Fermion-induced quantum criticality in two-dimensional Dirac semimetals: Non-perturbative flow equations, fixed points and critical exponents

TL;DR

This paper investigates how quantum fluctuations of Dirac fermions influence the nature of phase transitions in 2D Dirac semimetals, revealing conditions under which transitions change from first to second order and calculating associated critical exponents.

Contribution

It provides a non-perturbative analysis of fermion-induced quantum criticality in 2+1 dimensions, identifying the critical fermion number where the transition's nature changes and computing critical exponents.

Findings

Quantum fluctuations can turn a first-order transition into a continuous one.

Critical fermion number for transition change is determined.

Large N_f results show the transition change occurs in 2+1 dimensions.

Abstract

We establish a scenario where fluctuations of new degrees of freedom at a quantum phase transition change the nature of a transition beyond the standard Landau-Ginzburg paradigm. To this end we study the quantum phase transition of gapless Dirac fermions coupled to a $\mathbb{Z}_3$ symmetric order parameter within a Gross-Neveu-Yukawa model in 2+1 dimensions, appropriate for the Kekul\'e transition in honeycomb lattice materials. For this model the standard Landau-Ginzburg approach suggests a first order transition due to the symmetry-allowed cubic terms in the action. At zero temperature, however, quantum fluctuations of the massless Dirac fermions have to be included. We show that they reduce the putative first-order character of the transition and can even render it continuous, depending on the number of Dirac fermions $N_f$. A non-perturbative functional renormalization group approach is employed to investigate the phase transition for a wide range of fermion numbers. For the first time we obtain the critical $N_f$, where the nature of the transition changes. Furthermore, it is shown that for large $N_f$ the change from the first to second order of the transition as a function of dimension occurs exactly in the physical 2+1 dimensions. We compute the critical exponents and predict sizable corrections to scaling for $N_f =2$.

Figures

Tables

Table 1. Table 1: Numerical values for the largest three critical exponents for different N f subscript 𝑁 𝑓 N_{f} in D = 2 + 1 𝐷 2 1 D=2+1 in LPA 12 ′ superscript 12 ′ 12^{\prime} . The inverse correlation length exponent is given by θ 1 = ν − 1 subscript 𝜃 1 superscript 𝜈 1 \theta_{1}=\nu^{-1} .The critical exponent deciding over stability θ 2 subscript 𝜃 2 \theta_{2} is printed in boldface and also determines the corrections to scaling in the quantum critical regime.

$N_f$	${\theta }_1$	${\theta }_2$	${\theta }_3$
1/2	1.0620	+0.5130	-0.9021
1	0.8248	+0.1574	-0.8044
2	0.8623	-0.00497	-0.8779
3	0.9017	-0.02398	-0.9203
4	0.9237	-0.02646	-0.9413
5	0.9377	-0.02560	-0.9536
10	0.9672	-0.01793	-0.9775
20	0.9831	-0.01054	-0.9889
50	0.9931	-0.00465	-0.9956
$\infty$	1	0	-1

Table 2. Table 2: Correlation length exponent ν = θ 1 − 1 𝜈 superscript subscript 𝜃 1 1 \nu=\theta_{1}^{-1} and anomalous dimensions for the emergent SUSY scenario, i.e. N f = 1 / 2 subscript 𝑁 𝑓 1 2 N_{f}=1/2 in the U ​ ( 1 ) 𝑈 1 U(1) symmetric model. The epsilon expansion results ( ϵ 3 superscript italic-ϵ 3 \epsilon^{3} ) to three-loop order are taken from Ref.  Zerf2016, with direct substitution ϵ = 1 italic-ϵ 1 \epsilon=1 . The conformal bootstrap (cBS) results have been calculated in Ref.  Bobev:2015vsa, . The FRG results are obtained within LPA12 ′ .

$N_f=1/2$	$\nu$	${\eta }_{\varphi }$	${\eta }_{\psi }$
${ϵ}^3$	0.985	1/3	1/3
cBS	0.917	1/3	1/3
FRG	0.941	0.354	0.323

Table 3. Table 3: Critical exponents of the Kekulé transition for different N f subscript 𝑁 𝑓 N_{f} in D = 2 + 1 𝐷 2 1 D=2+1 . The QMC values and 1-loop RG results are taken from Ref.  li2015, , the FRG results are calculated in LPA12 ′ . We added our calculation of η ψ subscript 𝜂 𝜓 \eta_{\psi} for completion.

$N_f$		$\nu$			${\eta }_{\varphi }$		${\eta }_{\psi }$
	QMC	${ϵ}^1$	FRG	QMC	${ϵ}^1$	FRG	FRG
$2$	1.04	1.25	1.16	0.71(3)	0.67	0.88	0.062
$3$	1.05	1.26	1.11	0.77(2)	0.75	0.92	0.038
$4$	1.12	1.25	1.08	0.80(4)	0.80	0.95	0.027
$5$	1.08	1.23	1.07	0.85(4)	0.83	0.96	0.021
$6$	1.07	1.22	1.06	0.89(4)	0.86	0.97	0.017

Equations

$${\mathcal{L}}_{\psi }=\bar{\psi }{\gamma }_{\mu }{\partial }_{\mu }\psi ,$$

$${\gamma }_0=\text{mathds}{1}_2\otimes {\sigma }_z,{\gamma }_1={\sigma }_z\otimes {\sigma }_y,{\gamma }_2=\text{mathds}{1}_2\otimes {\sigma }_x,$$

$${\gamma }_3={\sigma }_x\otimes {\sigma }_y,{\gamma }_5={\sigma }_y\otimes {\sigma }_y,$$

$${\mathcal{L}}_{\psi \varphi }=i\bar{h}\bar{\psi }\left({\varphi }_1{\gamma }_3+{\varphi }_2{\gamma }_5\right)\psi .$$

$$\psi \to e^{i\theta {\gamma }_{35}/2}\psi ,\bar{\psi }\to \bar{\psi }{e}^{-i\theta {\gamma }_{35}/2},\varphi \to e^{i\theta }\varphi ,$$

$${\mathcal{L}}_{\varphi }$$

$$+{\bar{g}}_5({\varphi }^3+{\varphi }^{\ast 3}){\left|\varphi \right|}^2+{\bar{g}}_6({\varphi }^3+{\varphi }^{\ast 3}{)}^2+{\bar{\lambda }}_6{\left|\varphi \right|}^6+\dots$$

$$[{\bar{m}}^2]$$

$$[{\bar{g}}_5]$$

$$[\bar{g}{]}_{\text{NGFP}}=1+\frac{ϵ}{2}+\text{fluctuations}\gtrless 0.$$

$${\beta }_i({\alpha }_1^{\ast },{\alpha }_2^{\ast },\dots )=0\forall i.$$

$${\partial }_t{\alpha }_i=$$

$${\theta }_1\ge {\theta }_2\ge ....$$

$${\partial }_t{\Gamma }_k=\frac{1}{2}\text{STr}\left[({\Gamma }_k^{(2)}+R_k{)}^{-1}{\partial }_t{R}_k\right],$$

$$\displaystyle\Big{(}\Gamma_{k}^{(2)}\Big{)}(p,q)=\frac{\overrightarrow{\delta}}{\delta\Phi(-p)^{T}}\Gamma_{k}\frac{\overleftarrow{\delta}}{\delta\Phi(q)}\,.$$

$${\Gamma }_k$$

$$\displaystyle\hskip 25.6073pt-\frac{1}{2}Z_{\phi,k}(\phi_{1}\partial_{\mu}^{2}\phi_{1}+\phi_{2}\partial_{\mu}^{2}{\phi}_{2})+U_{k}(\bar{\rho},\bar{\tau})\Big{]}\,.$$

$$\bar{\rho }$$

$$\bar{\tau }$$

$$u(\rho ,\tau )=k^{-D}U\left(\bar{\rho },\bar{\tau }\right),h^2=\frac{k^{D-4}}{Z_{\varphi ,k}{Z}_{\psi ,k}^2}{\bar{h}}_k^2,$$

$${\eta }_{\varphi }=-\frac{{\partial }_t{Z}_{\varphi ,k}}{Z_{\varphi ,k}},{\eta }_{\psi }=-\frac{{\partial }_t{Z}_{\psi ,k}}{Z_{\psi ,k}}.$$

$${\partial }_t{h}^2$$

$${\partial }_{t}u$$

$$u^{(n,m)}=\frac{{\partial }^{n+m}}{\partial {\rho }^n\partial {\tau }^m}u.$$

$$u(\rho ,\tau )=$$

$$m^2={\lambda }_{1,0},g={\lambda }_{0,1},8\lambda ={\lambda }_{2,0}.$$

$$\displaystyle\partial_{t}\lambda_{ij}=\left(\frac{\partial^{i+j}}{\partial\rho^{i}\tau^{j}}\partial_{t}u(\rho,\tau)\right)\Big{|}_{\rho=\tau=0}\,.$$

$${\partial }_t{h}^2=$$

$${\partial }_t{m}^2$$

$$+\frac{32v_D}{D}\left(1-\frac{{\eta }_{\varphi }}{2+D}\right)\left(\frac{9g^2}{(1+m^2{)}^3}-\frac{4\lambda }{(1+m^2{)}^2}\right),$$

$${\partial }_{t}g$$

$$+\frac{32v_D}{D}\left(1-\frac{{\eta }_{\varphi }}{2+D}\right)\left(\frac{12\lambda g}{(1+m^2{)}^3}-\frac{{\lambda }_{1,1}}{(1+m^2{)}^2}\right),$$

$${\partial }_t\lambda$$

$$\displaystyle+\frac{16v_{D}}{D}\left(1-\frac{\eta_{\phi}}{2+D}\right)\Bigg{(}\frac{162g^{4}}{(1+m^{2})^{5}}-\frac{216g^{2}\lambda}{(1+m^{2})^{4}}$$

$$+\frac{40{\lambda }^2}{(1+m^2{)}^3}+\frac{18g{\lambda }_{11}}{(1+m^2{)}^3}-\frac{9}{8}\frac{{\lambda }_{02}}{(1+m^2{)}^2}$$

$$\displaystyle-\frac{3}{16}\frac{\lambda_{30}}{(1+m^{2})^{2}}\Bigg{)}\,,$$

Full text

I Introduction

I.1 Quantum Phase Transitions

Correlated many-body systems display a variety of complex collective phenomena that are related to the appearance of phase transitions. In particular, quantum phase transitions (QPT)sachdev2011 , i.e. zero-temperature phase transitions driven by quantum fluctuations, are believed to be key to the understanding of unconventional properties of correlated many-body systems, even at finite temperature. For example QPTs are considered to play an important role in the phase diagram of high-temperature superconductors sachdev2003 . In general, the established theoretical framework for the description of phase transitions is based on the identification of an appropriate order parameter $\phi$ distinguishing a disordered phase with expectation value $\langle\phi\rangle=0$, from a phase with a finite symmetry-breaking order $\langle\phi\rangle>0$, as a function of some tuning parameterlandau2013statistical . In the case of a QPT, the tuning parameter is a quantity which controls the strength of the quantum fluctuations, e.g., an interaction strength, a doping level or – more formally – the value of Planck’s constant.

In a phase transition of second order, the expectation value of the order parameter $\langle\phi\rangle$ is a continuous function of the tuning parameter, and the transition point is referred to as a quantum critical point (QCP). In the vicinity of such a QCP an universal critical behavior emerges. In many cases, the critical behavior can successfully be accessed by Landau-Ginzburg-Wilson (LGW) theoryWegner:1972ih ; Wilson:1973jj . It describes the phase transition by a continuum field theory formulated exclusively in terms of the order parameter in the appropriate space-time and restricted by the symmetries of the underlying system. The overwhelming success of this theoretical paradigm notwithstanding, it appears that there are QPTs which demand concepts beyond the LGW paradigmsenthil2004 , deconfined quantum critical pointsSenthil1490 being a prime example. These should describe a second-order QPT between two distinct symmetry-broken phases, in contrast to the LGW theory that would suggest either a first order transition or a coexistence region, with both order parameters having a finite expectation value. The reason why the LGW theory could fail is that the deconfined QCP cannot be described purely in terms of order parameter fluctuations. Instead, new degrees of freedom such as spinons emerge right at the critical point, and a more natural formulation of the transition is provided directly in terms of these.

I.2 Fluctuation-induced critical points

Similarly to the emergent spinons at a deconfined QCP, gapless fermions can provide additional degrees of freedom specific to a given quantum phase transition, so that it becomes necessary to go beyond a description purely in terms of the order parameter fields. For example, the presence of gapless fermion fluctuations in continuous quantum phase transitions towards an order parameter with $O(N)$ symmetry severely modifies the quantitative estimates for its critical exponents and defines the fermionic or chiral universality classesRosenstein:1992 ; Rosenstein:1993zf . Moreover, the impact of the fermion fluctuations could be so strong that a transition becomes second order, although the pure LGW theory would predict a first order transition. In this case we speak of a fermion-induced quantum critical pointli2015 . More generally, a fluctuation-induced critical point occurs, when quantum or thermal fluctuations render a putatively first-order transition continuous. Evidence for a fermion-induced second order transition has been put forward in Refs. li2015, ; Scherer:2016zwz, ; 2016arXiv161007603J, in the context of the 2+1 dimensional honeycomb electrons near the transition to a Kekulé valence bond solid (VBS). A closely related scenario was also proposed in 3D double-Weyl semimetals for the transition to a nodal-nematic order2016arXiv160906313J . The formation of a Kekulé VBS has been observed in artificial graphene gomes2012 and graphene on a Copper substrate Gutierrez2016 . Theoretically, the Kekulé VBS and the nodal-nematic order are described in terms of a $\mathbb{Z}_{3}$-symmetric complex-valued order parameter field $\phi$. As a consequence of the discrete $\mathbb{Z}_{3}$ symmetry, a cubic term $\propto(\phi^{3}+\phi^{*3})$ is allowed in the LG free energy $F_{\text{LG}}$. If such a term is present, the minimum of $F_{\text{LG}}$ jumps discontinuously from $\phi=0$ to $\phi\neq 0$ at the phase transition, see Fig. 1. That is this Landau-type mean-field criterion suggests to expect a first-order phase transition in the considered systemPhysRevB.8.3419 . Indeed, for the three-state Potts model, which exhibits such a $\mathbb{Z}_{3}$ symmetry, the first-order behavior of the thermal phase transition has been corroborated for three and higher spatial dimensionsRevModPhys.54.235 . On the other hand, an exact result for the three-state Potts model in two dimensions suggests a second-order phase transition in contrast to the argument abovebaxter1973 . In this case the strength of fluctuations is increased by decreasing the dimensionality of the system.

In the case of the Kekulé VBS and the nodal-nematic order in Dirac materials, the formulation in terms of a $\mathbb{Z}_{3}$ order parameter field is insufficient. The strength of fluctuations is increased by the inclusion of additional critical degrees of freedom, such as gapless Dirac fermions. In fact, it is argued in Refs. li2015, ; Scherer:2016zwz, (Kekulé) and Ref. 2016arXiv160906313J, (nodal-nematic) that the first-order behavior breaks down and a critical point is induced if the number of fermions $N_{f}$, corresponding to the number of pairs of Dirac points, is sufficiently large (cf. also Ref. Rosenstein:1992, ). That is, there should be a critical value $N_{f,c}$ so that the transition is of first order for $N_{f}<N_{f,c}$ and of second order for $N_{f}>N_{f,c}$. However, the precise value of the critical fermion number is presently not known. Moreover, so far there is no non-perturbative calculation which would access both the first and the second order regimes. In Ref. 2016arXiv161007603J, Majorana quantum Monte Carlo simulations have been performed for $N_{f}\in\{2,\ldots,6\}$, which all yield a second order transition. On the other hand, Refs. Scherer:2016zwz, and 2016arXiv160906313J, used perturbative expansions around the upper critical dimension and large $N_{f}$ to study the fermion-induced critical point.

The question whether fluctuations can indeed alter the nature of a phase transition relative to the LGW prediction is inherently non-perturbative, as it requires sizable fluctuation effects. Thus a non-perturbative study which would access the full range of $N_{f}$ is required in order to establish the possible occurrence of a fermion-induced QCP. Here we therefore perform such a study through an investigation of the zero-temperature Kekulé transition in 2+1 dimensional Dirac materials with the help of the non-perturbative functional renormalization group (FRG)wetterich1993 . The FRG has been proven to be a suitable tool in the study of non-perturbative aspects of quantum and statistical field theories, see, e.g., Refs. Berges:2000ew, ; Polonyi:2001se, ; Pawlowski:2005xe, ; Gies:2006wv, ; Delamotte:2007pf, ; Braun:2011pp, ; metzner2011, for reviews and Sec. III for important aspects related to this work. Once the appearance of a fermion-induced QCP is established, the question arises if the vicinity to a first order transition affects the critical behavior in the second order transition regime close to $N_{f,c}$. We therefore compute the critical exponents characterizing the second order transition, and determine and discuss the corrections to scaling.

I.3 Main results

In this work we study whether strong fluctuations can induce a second order phase transition and consequently critical behavior. To that end, we consider the Kekulé ordering transition at zero temperature in two-dimensional Dirac semimetals, where the conventional Landau-Ginzburg paradigm of phase transitions would predict a discontinuous transition. By means of the non-perturbative functional renormalization group, we investigate a full range of the number of Dirac fermions $N_{f}$ for this model. We indeed find that fluctuations render the phase transition continuous in $D=2+1$ dimensions if $N_{f}$ is sufficiently large. This is sketched in Fig. 1. We also obtain the first direct estimate for the critical fermion flavor number, where the nature of the transition changes. We find $N_{f,c}\approx 1.9$, which lies in the limits provided by quantum Monte Carlo ($N_{f}<2$) and SUSY ($N_{f,c}>1/2$) calculations. In particular, spin-1/2 fermions on the honeycomb lattice, $N_{f}=2$, show a second order transition with quantum critical behavior. In the critical regime $N_{f}>1.9$, we determine the critical exponents characterizing the scaling of correlation functions. We find the second largest critical exponent to be negative but small in magnitude, which implies a nearly marginal RG direction in the theory. In the limit of a large number of Dirac fermions, the corresponding coupling becomes exactly marginal. This leads to sizable corrections to scaling close to the quantum critical point for the whole $N_{f}$-regime. We expect these corrections to be detectable by QMC simulations.

I.4 Outline

We use a Gross-Neveu-Yukawa theory to model the quantum phase transition of massless Dirac fermions towards the Kekulé valence bond solid (VBS), which we motivate in Sec. II. The Kekulé VBS is described in terms of a complex order parameter field with $\mathbb{Z}_{3}$ symmetry, Yukawa-coupled to Dirac fermions. In Sec. III, we explain the renormalization group picture of phase transitions and introduce renormalization group fixed points and critical exponents. This allows a more specific notion about what sizable fluctuations mean in the considered context. We also present our functional RG approach and discuss its key features. In Sec. IV, we discuss our results for the fermion-induced quantum critical point. We carefully establish convergence within our truncation scheme and compare to previous studies if possible. In particular, we also highlight the advantages of the approach presented here. We give our best estimate for the fermion number, where the nature of the transition changes, and determine the critical behavior in the second order regime. Finally, we draw conclusions in Sec. V. Technical details are relegated to the appendices.

II Semimetal to Kekulé VBS transition

II.1 Honeycomb fermions and Dirac semimetals

Since the discovery of graphene, considerable interest has been generated by Dirac materials, that impact a broad range of research from fundamental physics to concrete technological applicationsvafek2013 ; wehling2014 . More specifically, charge-neutral grapheneNovoselov2005 ; Geim2007 ; castroneto2009 serves as a the prototypical example for a class of materials where low-energy Dirac excitations emerge as a result from the underlying honeycomb lattice structure. A wide variety of possible ordering patterns have been proposed to become relevant to the phase diagram of interacting Dirac electrons on the honeycomb latticesorella1992 ; herbut2006 ; hou2007 ; honerkamp2008 ; raghu2008 ; herbut2009 ; ryu2009 ; roy2010c ; weeks2010 ; classen2014 ; volpez2016 ; sanchez2016 . Among these ordered states, the Kekulé valence bond solid (VBS) hou2007 ; ryu2009 ; roy2010c ; classen2014 – a particular dimerization pattern of the fermions on the honeycomb lattice – has recently attracted special attention due to its experimental realization gomes2012 ; Gutierrez2016 .

Explicitly, we consider two-dimensional materials whose low-energy effective theory can be expressed in terms of a free Dirac LagrangianSemenoff1984

[TABLE]

where the conjugate field is defined as $\bar{\psi}=\psi^{\dagger}\gamma_{0}$ and $\partial_{\mu}=(\partial_{\tau},\vec{\nabla})$ denotes the imaginary-time and space derivative. We have set the Fermi velocity to unity, $v_{F}=1$. The $\gamma$ matrices satisfy the Clifford algebra $\{\gamma_{\mu},\gamma_{\nu}\}=2\delta_{\mu\nu}$ and their four-dimensional representation in $D=2+1$ dimensions reads

[TABLE]

with the Pauli matrices $\sigma_{i}$. In momentum space, we have $\psi(x)=\int d^{D}qe^{iqx}\psi(q)$ with momentum vectors $q=(\omega,\vec{q})$ gathering Matsubara frequency $\omega$ and wavevector $\vec{q}$. The four spinor components of the Dirac field $\psi$ can be related to the two sublattice and the $K,-K$ valley degrees of freedom of spinless fermions on the honeycomb lattice. Furthermore we define

[TABLE]

which anticommute with all $\gamma_{\mu}$. Their combination given by $\gamma_{35}=-i\gamma_{3}\gamma_{5}$ commutes with all $\gamma_{\mu}$ and anticommutes with $\gamma_{3}$ and $\gamma_{5}$. As a convenient generalization of this model, we introduce the number of fermion flavors of four-component spinors by $\psi^{\dagger}\rightarrow\psi_{i}^{\dagger}$ with $i\in\{1,\ldots,N_{f}\}$. This corresponds to an arbitrary number of flavors or of pairs of Dirac points $\pm\vec{K}_{i}$ in the fermion spectrum. $N_{f}=2$ then corresponds to an eight-component spinor, which describes the “graphene case” of spin-1/2 fermions on the honeycomb lattice with spin-rotation invariance.

II.2 Kekulé valence bond solid

We are interested in the phase transition from the above Dirac-semimetal phase into the Kekulé valence bond solidhou2007 . A sketch of this state on the honeycomb lattice is shown in Fig. 1, exhibiting varying bond strengths between nearest-neighbors. In other words, the uniform nearest-neighbor hopping amplitude of the symmetric phase receives two different, but real values in the symmetry-broken phase, arranged in the specific pattern shown in the figure. This pattern breaks the translational symmetry of the original lattice and reduces the sixfold rotation symmetry to a threefold one. In the low-energy effective theory, the modulation of the nearest-neighbor hopping is described by a complex order parameter $\phi=(\phi_{1}+i\phi_{2})/\sqrt{2}$, whose phase controls the angles of the dimerization patternryu2009 and in our case has to be $\mathbb{Z}_{3}$ symmetric. This is related to the presence of three degenerate Kekulé ground states as obtained by translation of the pattern shown in Fig. 1 along the lattice vectors. The $\mathbb{Z}_{3}$ symmetric order parameter couples to the fermions via a Yukawa term of the form roy2010c

[TABLE]

We note that this part of the action is still invariant under the general “chiral” continuous $U(1)$ transformation

[TABLE]

corresponding to translations on the honeycomb lattice. This property will be violated by the bosonic action

[TABLE]

in which the terms $\propto(\phi^{3}+\phi^{*3})$ are responsible for the reduction of the continuous $U(1)$ symmetry down to $\mathbb{Z}_{3}$. In general higher order couplings can also occur as denoted by the ellipsis. Here, we have added couplings up to the order $\phi^{6}$. Dimensional analysis yields the following canonical dimensions for the bosonic couplings

[TABLE]

as well as $[\bar{h}]=(4-D)/2$ for the Yukawa coupling. Accordingly, couplings of higher order than $\phi^{4}$ are omitted in perturbative calculations close to $D=3+1$ dimensions as they become irrelevant (see next section). In three spacetime dimensions, however, couplings up to $\phi^{6}$ are relevant, and we will indeed find that they do play a role in deciding the nature of the phase transition.

In equation (6), we employed the low-energy effective particle-hole symmetry of the model to exclude a linear term in the time-derivative of the order parameter Lagrangian. Furthermore, we have written down a Lorentz-symmetric kinetic term, i.e. we have set the boson velocity to the same value as the Fermi velocity $v_{B}=v_{F}=1$. This Lorentz-symmetric form of the total action $S=\int d\tau d^{D-1}x\,(\mathcal{L}_{\psi}+\mathcal{L}_{\psi\phi}+\mathcal{L}_{\phi})$ however, is not dictated a priori; rather, the Lorentz symmetry near fermionic quantum critical points with equal velocities for fermions and bosons has been argued to emerge naturally in the deep infrared regime in a large class of Yukawa theories of the same kindRoy:2015zna , even if $v_{F}\neq v_{B}$ on intermediate scales. Lorentz symmetry also emerges in our case, as the quantum critical point which we will find is located at vanishing $U(1)$-breaking couplings, i.e. $\bar{g}^{\ast}=\bar{g}_{5}^{\ast}=\bar{g}_{6}^{\ast}=...=0$ and therefore the analysis given in Ref. Roy:2015zna, is directly applicable. Away from the QCP, and in particular in the symmetry broken phase, however, it could also be interesting to consider the case with two different velocities. We will leave this issue for future work.

III Renormalization Group approach

III.1 Renormalization group scaling argument

Before we explain the details of our RG approach, we would like to translate the picture from the saddle-point approximation of the free energy into a renormalization group language. In the RG picture, a system is considered scale-dependent, i.e. the couplings describing it evolve with the scale when, for example, the energy scale is lowered. Critical behavior near a second-order phase transition is related to the presence of an infrared attractive renormalization group fixed point of this scale evolution. At a RG fixed point, the system becomes scale invariant by definition. Consequently, the free energy adopts a scaling form accompanied by universal critical exponents. In the standard scenario, the critical point is approached by fine-tuning of a parameter, e.g., the temperature $T\to T_{c}$ for a thermal transition and, e.g., a doping level or a coupling strength for a quantum transition. In the renormalization group a tuning parameter translates to the presence of an RG relevant direction, which is infrared repulsive, i.e. the corresponding direction also has to be fine-tuned to approach an otherwise infrared-attractive RG fixed point. Therefore, in order to find a scaling solution and observe critical behavior in a system, an identical number of physical tuning parameters and RG relevant directions is required.

Dimensional analysis, or power counting, gives a first hint about the RG relevance of the parameters of a model and is valid near the Gaußian fixed point. At a non-Gaussian fixed point, however, the predictions from dimensional analysis about the RG scaling of parameters receive fluctuation corrections. They can affect the relevance of the RG directions. A positive (negative) power counting dimension suggests a relevant (irrelevant) parameter. Parameters with vanishing power counting dimension are called marginal. In our case, the canonical dimensions for the lowest order couplings appearing in the Lagrangian in Eq. (6) are presented in Eq. (7). Here, $D$ specifies the spacetime dimension of the system. The power counting dimension of the parameter $m^{2}$ is strongly relevant and certainly needs to be fine-tuned to approach the critical point. In contrast, the other power counting dimensions depend on $D$, so the consideration is more subtle: For example, in $D=4-\epsilon$ dimensions, we observe $[\bar{\lambda}]=\epsilon$ and the coupling $\bar{\lambda}$ therefore corresponds to a slightly relevant direction controlled by the size of $\epsilon$. On the other hand, the critical behavior of this type of continuous field theories below four dimensions is typically governed by non-Gaußian fixed points (NGFP)wilson1972 and the dimensional counting is modified. The sign and the magnitude of the fluctuation contribution generally depends on the dimensionality and possibly on the coupling to other degrees of freedom. When these modifications are large enough, they can change the sign of the power counting dimension of a canonically relevant coupling. In fact, at the non-Gaussian $O(2)$ Wilson-Fisher fixed point, the $\lambda$ direction is turned irrelevant by a negative contribution from the fluctuations, i.e. $[\lambda]_{\text{NGFP}}<0$, and therefore $\lambda$ automatically approaches its fixed point value in the infrared. It turns out that this reasoning continues to apply even to the case of three dimensionszinnjustin1998 , when $\epsilon=1$.

In the $\mathbb{Z}_{3}$ symmetric scenario, however, we have another term in the free energy, i.e. the cubic order parameter term $\propto\bar{g}$. It generically introduces an additional RG relevant direction with a larger power-counting dimension, e.g., in $4-\epsilon$ dimensions we have $[\bar{g}]=1+\epsilon/2$. Rendering this direction irrelevant to obtain a fluctuation-induced critical point requires a more dominant role of fluctuations. This cannot be controlled by a small $\epsilon$ as

[TABLE]

Therefore, this scenario needs the fluctuations to provide a large negative contribution with $|\text{fluctuations}|>3/2$ in the case of $D=2+1$ independent of their origin. This constitutes an intrinsically non-perturbative scenario and motivates us to use the functional renormalization group in the following.

III.2 Fixed points and critical exponents

More explicitly, the renormalization group theory describes the scale dependence of a physical system by providing $\beta$ functions for the different couplings of a modelherbutbook . The $\beta$ functions are differential equations that encode the evolution of the system with respect to the energy (or momentum) scale $k$. Starting from a “microscopic” model for a system at some ultraviolet (UV) cutoff scale $k=\Lambda$, one can then infer the low-energy, or infrared (IR), characteristics in terms of the solution of the $\beta$ functions. In our case, the UV scale $\Lambda$ corresponds to the scale at which our effective model is valid, e.g., the largest energy scale where a description of the material in terms of Dirac excitations is justified. This is much smaller than the system’s bandwidth and going beyond that energy scale would require an appropriate lattice description.

To formulate the renormalization group approach, we introduce the generalized set of dimensionless couplings for the theory by $\alpha_{i}$, $i\in\{1,2,\ldots\}$ which will be specified below. The $\beta$ functions can be written in the form $\partial_{t}\alpha_{i}=\beta_{i}(\alpha_{1},\alpha_{2},\ldots)$, where the change in scale is written in terms of the renormalization group time $t=\ln(k/\Lambda)\leq 0$. A renormalization group fixed point $\alpha^{*}$ is defined by a simultaneous vanishing of all beta functions of the model,

[TABLE]

The critical properties and scaling behavior near such a transition are encoded in the linearized RG flow close to the fixed point $\alpha^{*}$,

[TABLE]

where $\mathcal{B}_{i,j}=(\partial\beta_{i}/\partial\alpha_{j})|_{\alpha=\alpha^{*}}$ is the stability matrix. The eigenvalues $\theta_{i}$ of $(-\mathcal{B}_{i,j})$ are called the critical exponents and are universal quantities that characterize the scaling laws at the putative continuous phase transition. For the following discussion, we introduce an ordering according to the size of the critical exponents, i.e.

[TABLE]

All positive critical exponents $\theta_{i}$ correspond to RG-relevant directions, i.e., the fixed point repels the flow in that direction. In turn, negative $\theta_{i}$ are RG irrelevant and correspond to attractive directions. Fixed points with no more than one relevant direction (corresponding to no more than one positive critical exponent) can be accessed by tuning a single parameter. In the present context, the system features a single tuning parameter. Therefore, a stable infrared-attractive fixed point which gives rise to physically accessible quantum critical behavior is specified by a set of critical exponents where only one eigenvalue is positive, $\theta_{1}>0$, and all others are negative, $\theta_{i}<0$ for $i>1$. Thus, the sign of the second largest critical exponent $\theta_{2}$ as calculated from the stability matrix decides over the stability of a fixed point.

III.3 Functional RG approach in a nutshell

For our study on the quantum transition from the semi-metallic state of two-dimensional Dirac fermions to a Kekulé VBS state, we employ the functional renormalization group (FRG) approachwetterich1993 . It generalizes Wilson’s momentum-shell RG approach by modifying the functional integral representation of the partition function $Z=\int_{\Lambda}\mathcal{D}\Phi\,\exp(-S[\Phi])$ by a regulator insertion in the microscopic action, i.e. $S\rightarrow S+\int_{p}\frac{1}{2}\Phi(-p)R_{k}(p)\Phi(p)$. Here, $\Phi$ represents a collective field variable for all field degrees of freedom of a specific model. In our case it includes both, Dirac fermions as well as complex bosons. The regulator insertion controls how high energy degrees of freedom are integrated out in the RG procedure and the common Wilsonian momentum-shell RG integration would correspond to a specific choice of $R_{k}(p)$. In general the regulator $R_{k}(p)$ depends on an IR cutoff scale $k$ and the momentum $p$ of the field configurations that are integrated over in the partition function. In particular, $R_{k}(p)$ is designed to suppress the low-momentum fluctuations by choosing $R_{k}(p)>0$ for $p^{2}<k^{2}$. Further, the regulator function satisfies $R_{k}(p)\rightarrow\infty$ for $k\rightarrow\Lambda\rightarrow\infty$ and $R_{k}(p)\rightarrow 0$ for $k/|p|\rightarrow 0$. We then define the generating functional for the one-particle irreducible correlation functions, i.e., the scale-dependent effective action $\Gamma_{k}$ as the (modified) Legendre transform of the Schwinger functional $\ln Z_{k}$. By choice of the regulator properties, $\Gamma_{k}$ then interpolates between the microscopic action $S$ and the quantum effective action $\Gamma_{k\rightarrow 0}$ in the IR.

The central equation of the FRG is the Wetterich equation: an exact functional differential equation which describes the RG evolution of $\Gamma_{k}$. It provides the explicit interpolation procedure between $S$ and $\Gamma$. The Wetterich equation readswetterich1993

[TABLE]

where $\partial_{t}=k\partial_{k}$ and we have introduced the Hessian $\Gamma_{k}^{(2)}$

[TABLE]

In terms of the exact Wetterich equation, the functional RG approach allows to devise suitable truncations for a practical implementation of Wilson’s idea of successively integrating out degrees of freedom in the functional integral formalism.

III.4 Key features of the FRG

In general, the FRG provides a unified framework to describe universal and non-universal properties of quantum and statistical field theories within and beyond the scope of perturbation theory. That means, given an appropriate model, it can not only access the universal behavior in the vicinity of a continuous phase transition, but also the system’s physical properties away from the transition. In particular, first-order phase transitions can also be studiedtissier2000 ; jakubczyk2009 ; braun2010 ; schaefer2006 ; rennecke2014 .

The Wetterich equation has a one-loop structure, resembling the one-loop functional determinant. Therefore, a proliferation of diagrams as they arise beyond leading order in higher-order loop expansions is avoided with the FRG. On the other hand, the presence of the Hessian in the denominator of the one-loop expression $\Gamma_{k}^{(2)}$ encodes higher-order effects within the one-loop structure. For example, the threshold effects introduced by this non-perturbative propagator improve conventional one-loop results for critical exponents while maintaining the simple one-loop form. Furthermore, the FRG can be applied in arbitrary dimension and does not require the identification of a small expansion parameter. Instead, approximations in the FRG approach result from of a particular truncation of the space of operators that generate the RG flow, see Sec. IV.2. In contrast to perturbative RG approaches where typically highly non-trivial resummation techniques are required in order to obtain quantitative results for critical exponents, the FRG method already incorporates a specific resummation through the threshold effectsdelamotte2005 ; delamotte2010b ; delamotte2010c . Moreover, the inclusion of chiral fermions is straightforward within the FRG and has been successfully employed in the context of Yukawa models for Dirac materialsRosa:2000ju ; Hofling:2002hj ; Gies:2009da ; Janssen:2012pq ; mesterhazy2012 ; janssen2014 ; Classen:2015mar . In contrast to lattice field theory, we employ a continuum formulation here, which means that our results are unaffected by finite-volume or discretization artifacts. Thus, the FRG constitutes a well-suited tool to investigate interacting RG fixed points, in particular in regimes where the validity of the perturbative RG approach is called into question. From a more general point of view, results from the perturbative RG, functional RG and numerical methods, as, e.g., (quantum) Monte Carlo calculations, complement each other to provide a more comprehensive understanding of physical systems Mihaila:2017ble ; otsuka2016 .

IV Fermion-induced QCP from FRG

IV.1 FRG flow equations

The scale-dependent effective action $\Gamma_{k}$ generally contains all symmetry-compatible operators of the considered theory. In practice, the application of Wetterich’s equation requires truncating the infinite-dimensional space of RG directions or couplings to a closed and tractable subspace. Here, we employ a scheme following a derivative expansion which we truncate after the leading order. It is based on the construction of our model in Sec. II, and the ansatz for $\Gamma_{k}$ explicitly reads

[TABLE]

with the spacetime dimension $D$. The Dirac spinors $\bar{\psi},\psi$ come with a number of $N_{f}d_{\gamma}$ spinor components, where $d_{\gamma}$ denotes the dimension of the Clifford algebra. $\phi_{1}$ ($\phi_{2}$) is the real (imaginary) part of a complex scalar field $\phi=(\phi_{1}+i\phi_{2})/\sqrt{2}$. In the first and second line, we gather the kinetic terms of the boson and fermion fields which are equipped with a scale-dependent, uniform wave function renormalization $Z_{\psi/\phi,k}$. Additionally, in the first line, we list the Yukawa coupling $\bar{h}_{k}$ which also carries a scale-dependence. Finally, we have introduced the scale-dependent boson effective potential $U_{k}$ which includes all symmetry-allowed boson self-couplings. We write it as a function of two $\mathbb{Z}_{3}$ invariants

[TABLE]

For the determination of the renormalization group fixed point and its properties it is particularly convenient to introduce dimensionless quantities. To that end, we define the dimensionless effective potential and the dimensionless Yukawa coupling as

[TABLE]

with $\bar{\rho}=k^{D-2}Z_{\phi,k}^{-1}\rho$ and $\bar{\tau}=k^{3(D-2)/2}Z_{\phi,k}^{-3/2}\tau$. Dimensionfull quantities are denoted with a bar, dimensionless without. Further, we introduce the boson and fermion anomalous dimensions

[TABLE]

We proceed with the evaluation of the Wetterich equation to calculate explicit expressions for the flow of the effective potential $u$, the squared Yukawa coupling $h^{2}$ and the anomalous dimensions $\eta_{\phi},\eta_{\psi}$. To do so, we substitute our ansatz for $\Gamma_{k}$, Eq. (15), into the Wetterich equation, Eq. (13), and project it onto the different couplings. This leads to a set of coupled differential equations of the form

[TABLE]

together with two algebraic equations for $\eta_{\phi}$ and $\eta_{\psi}$. Here, $u^{(n,m)}$ denotes the derivative of the effective potential with respect to the invariants, i.e.

[TABLE]

For the definition of the scale-dependent boson couplings from the effective potential, we expand $u$ in powers of $\rho$ and $\tau$. Generally, to properly access the fixed-point properties, we expand the effective fixed-point potential about its scale-dependent minimum. In the considered model, it turns out that the minimum of the fixed-point potential typically lies at the origin $(\rho,\tau)=(0,0)$ except for very small $N_{f}\lesssim 1/2$. We therefore choose an expansion corresponding to a vanishing vacuum expectation value of the bosonic field, which we refer to as the expansion in the symmetric (SYM) regime,

[TABLE]

This expansion includes the (dimensionless versions of the) previous couplings

[TABLE]

The $\beta$-functions for the expansion parameters $\lambda_{i,j}$ in the SYM regime can then be obtained by appropriate projection prescriptions for the flow of the dimensionless effective potential $\partial_{t}u$. That is we obtain the flow of the bosonic couplings in the SYM regime from

[TABLE]

Accordingly, we set $\rho=\tau=0$ in the expressions for $\partial_{t}h^{2}$, $\eta_{\phi}$ and $\eta_{\psi}$. For the explicit fixed-point solutions, we expand the effective potential up to a finite order in the boson fields. That is, it includes all symmetry-allowed couplings up to order $\phi^{N}$ where, in practice, $N\in\{4,6,8,12\}$. We can then also study the dependence of our results on $N$ to establish a convergence within the potential expansion. We refer to the approximation including the effective potential up to order $N$, together with the running wave function renormalizations and Yukawa coupling as LPA$N^{\prime}$ in the following. We explain further details of the derivation for the flow equations including projection prescriptions and regulator choices in App. A. We also present the explicit and fully analytical expressions for the $\beta$ functions for an expansion with general vacuum expectation value of the field in App. A.

IV.2 FRG flow equations in the symmetric regime

For vanishing expectation value the flow of the Yukawa coupling is completely given by its dimensional running

[TABLE]

because the loops with $\phi_{1}$ and $\phi_{2}$ propagator exactly cancel each other.

As example of the flow of the boson couplings, we list the flow equations of the lowest order couplings of Eq. (23)

[TABLE]

Their graphical representation is given in Fig. 2. At these expressions, one can recognize the FRG enhancement by the dressing of the propagators with the anomalous dimensions. The anomalous dimensions are given by the solution of

[TABLE]

IV.3 Emergent $U(1)$ symmetry

We can now directly evaluate the renormalization group equations for arbitrary dimension $D$ and fermion flavor number $N_{f}$. This allows us to determine the fixed point structure and the critical behavior. The fixed point corresponding to the potentially continuous Kekulé transition exhibits an enhanced $U(1)$ symmetry. This has also been observed in the perturbative RG calculationsli2015 ; Scherer:2016zwz . Therefore, all couplings that reduce the symmetry to $\mathbb{Z}_{3}$ vanish at the fixed point. More specifically, this allows us to reduce the set of couplings to be solved for in the fixed point equations. We can directly use

[TABLE]

and solve for the NGFP solution in the $\lambda_{i,0}$ couplings. For the stability matrix, however, all the couplings have to be taken into account as the RG scaling generally is different from zero, also for the vanishing couplings. The diagonalization of the stability matrix is simplified by an effective decoupling of the $U(1)$ symmetric and the $U(1)$-symmetry-broken sector at the fixed point leading to a block diagonal structure.

IV.4 Connection to perturbative RG

In order to establish a connection to the perturbative RG approach, we compare our results for the universal critical exponents to a perturbative expansion in $4-\epsilon$ dimensions to first order in $\epsilon$, cf. Ref. Scherer:2016zwz, . Since this expansion becomes exact close to the upper critical dimension, our FRG results should reduce to this limitjanssen2014 ; Classen:2015mar ; eichhorn2016 . To confirm this, we have calculated the universal critical exponents for dimensions $D\in[3,4]$. In Fig. 3, we present the cases $N_{f}=1$ and $N_{f}=10$, showing that the FRG results indeed approach the values of the $\epsilon$-expansion in the perturbative regime. However, for dimensions close to three, which is relevant for two-dimensional Dirac materials, deviations increase. In fact, the second critical exponent remains positive for $N_{f}=1$ in contrast to the $\epsilon$-expansion result. For $N_{f}=10$, $\theta_{2}$ changes sign, but not before $D\approx 3.01$. In this case, the $U(1)$ fixed point becomes stable in $D=3$ and renders the putative discontinuous Kekulé transition second order. However in comparison to first order in $\epsilon$, the tendency to a fluctuation induced second order transition is still reduced. This general trend can be attributed to the threshold effects from the FRG approach, i.e. the mass contributions in the propagators. Thus the impact of fluctuations is reduced by the diminution of loop contributions. In contrast, these effects are not included in a perturbative expansion to one-loop order, so that fluctuations there tend to be overestimated.

IV.5 Fermion-induced criticality

IV.5.1 Fixed point stability for varying numbers of fermions

To determine if the nature of the Kekulé phase transition can be changed from first to second order in the case of two-dimensional Dirac materials as graphene, we calculate the critical exponents characterizing the $U(1)$-symmetric fixed point in $D=2+1$ for different numbers of pairs of Dirac points $N_{f}$. As we explained in Sec. III, we determine the second order regime by analyzing the sign of the second largest critical exponent, which decides over the stability of a fixed point. If it becomes negative, the phase transition becomes second order and universal critical behavior emerges. The largest critical exponent is related to the correlation length exponent $\theta_{1}=\nu^{-1}$. It is always positive. With respect to the second largest exponent, we find that it changes sign at a critical $N_{f,c}$. Below $N_{f,c}$, it is positive in our calculations so that a first order transition takes place. Above the critical number of fermions, $\theta_{2}$ is always negative. In this regime, fluctuations are strong enough to change the sign of $\theta_{2}$ and make the $U(1)$ fixed point stable. All further critical exponents of the $U(1)$-symmetric fixed point are negative in our calculation. Our best estimate for the critical fermion flavor number is $N_{f,c}\approx 1.9$. This lies in the range of upper and lower bounds obtained by QMC calculationsli2015 , $N_{f,c}<2$, and an emergent SUSY theory2016arXiv161007603J ; Zerf2016 , $N_{f,c}>1/2$. Further, we show the largest critical exponents $\theta_{1}$ and $\theta_{2}$ as function of the number of fermions in Fig. 4.

To investigate the effect of higher order couplings, we have determined the critical behavior for different expansion orders of the order parameter potential. More specifically we have taken into account all symmetry-allowed couplings up to $\phi^{N}$ for $N\in\{4,6,8,12\}$. Interestingly, we find that for small $N_{f}$ a $\phi^{4}$ expansion as usually employed in perturbative calculations is not sufficient. It differs substantially from higher-order expansions as can be cleary seen in Fig. 4. In contrast, we observe a convergence of our results for orders higher than $\phi^{6}$. These findings confirm our expectations that in $D=2+1$ dimensions and for not too large fermion flavor numbers, the contribution of higher order couplings is not negligible, because couplings up to $\phi^{6}$ are relevant at the microscopic level. Although at the interacting fixed point higher order eigenperturbations become irrelevant, their scaling still contributes to the second critical exponent through non-diagonal entries in the stability matrix. Finally, for increasing $N_{f}$, boson loops become sub-leading compared to fermion contributions so that the difference between the expansion orders of the boson potential stops playing a role.

IV.5.2 Large $N_{f}$

For large $N_{f}$, we can solve the fixed point equations analytically. This allows us to obtain exact results for the critical exponents. In the symmetric regime for large $N_{f}$, the FRG equations simplify to

[TABLE]

where we have rescaled $u\rightarrow u/N_{f}$, $Z_{\phi}\rightarrow Z_{\phi}/N_{f}$. The fixed point that we are interested in has nonzero $h^{2}$. Note, that the solution with $h^{2}=0$ is unstable, because it has an additional relevant eigenperturbation in the direction of $h^{2}$. Consequently, from $\partial_{t}h^{2}=0$ follows that $\eta_{\phi}=4-D$. Furthermore, $\partial_{t}h^{2}$ only depends on $h^{2}$ and

[TABLE]

only depends on $u^{(n,m)}$ and $h^{2}$. Therefore, the stability matrix becomes triangular and we can read off its eigenvalues from the diagonal entries. That means the critical exponents are given by the scaling dimensions of the couplings at the interacting fixed point. It follows for the boson couplings $\theta_{i}=[u^{(n,m)}]=D-(2n+3m)[\phi]=D-(2n+3m)\frac{1}{2}(D-2+\eta_{\phi})$ and we obtain

[TABLE]

with $\eta_{\phi}=4-D$. We show how the largest critical exponents approach these values in Fig. 4. In particular, we find for the cubic coupling

[TABLE]

Thus for large $N_{f}$, the cubic coupling is exactly marginal at the interacting fixed point leading to a line of fixed points for arbitrary $\lambda_{0,1}$. We can also see this at the exact solution of Eq. (33)

[TABLE]

with the hypergeometric function ${{}_{2}F_{1}}$ and $h^{2}_{\ast}=D(8-6D+D^{2})/(2v_{D}d_{\gamma}(4-3D))$. The arbitrary function $f(\tau/\rho^{3/2})$ must be determined by boundary conditions. But we see that the cubic term $\lambda_{0,1}\tau=\lambda_{0,1}\rho^{3/2}\frac{\tau}{\rho^{3/2}}$ is of the form $\rho^{D/2}f(\tau/\rho^{3/2})$ in $D=3$. That is it satisfies the fixed point equation for every $\lambda_{0,1}$. As a consequence of Eq. (39), a second order phase transition can only occur for $D\leq 3$ for large $N_{f}$.

IV.5.3 Nature of the phase transition for general $D$ and $N_{f}$

We have also performed the stability analysis of the fixed point which describes the phase transition to the $\mathbb{Z}_{3}$-symmetric Kekulé state for varying spacetime dimensions. This allows to estimate how close the theory in 2+1 spacetime dimensions is to the other regime, i.e. if the transition is weakly first or second order, respectively. Our result is shown in Fig. 5. The second order regime appears close to three space-time dimensions for $N_{f}\gtrsim 1.9$. Compared to the first order $\epsilon$-expansion Scherer:2016zwz , it is severely reduced. As we explained above, the reason is that higher order couplings play a role and that threshold effects coming from the full propagators in the fRG equations reduce fluctuation effects. Nevertheless, systems with $N_{f}=2$ and $D=2+1$ still lie in the range with a stable fixed point. That means that, e.g., spin-1/2 fermions on the two-dimensional honeycomb lattice show a second order transition with critical scaling. The scaling is determined by the correlation length exponent $\theta_{1}$, the exponent of the cubic anisotropy $\theta_{2}$ and the anomalous dimensions $\eta_{i}$ see tables 1 and 3. However, as a result of being near to the first order transition, we can expect substantial corrections to this scaling in the whole second order regime as explained in the next section.

IV.6 Critical behavior in 2D Dirac materials

IV.6.1 Corrections to scaling

We list our best numerical values for the critical exponents, i.e. the eigenvalues of the stability matrix, as obtained within LPA$12^{\prime}$ in Tab. 1. In the physical case of 2+1 dimensions, we see from Tab. 1 and Fig. 4 that the magnitude of $\theta_{2}$ is generally small, i.e. $|\theta_{2}|\ll 1$. This implies that there is an RG direction which flows very slowly, i.e. almost logarithmically, towards its quantum critical point. Therefore, we expect that the corrections to scaling should be easily visible in the vicinity of the quantum critical point.

Generally, the leading correction to scaling at a continuous phase transition is induced by the least irrelevant RG direction at the critical point. Here, this RG direction is related to the cubic coupling $g$ and indeed represented by the critical exponent $\theta_{2}$ which gives rise to a power-law correction to the scaling. For example, the scaling of the correlation length should follow a behavior of the formwegner1972b ; herbutbook ; assaad2013

[TABLE]

where $A,B$ are non-universal amplitudes, and in particular $B\propto g$. The ellipsis indicate further corrections due to the stronger irrelevant directions, which vanish very quickly; see the values of $\theta_{3}$ in Tab. 1. $\Delta$ marks the distance to the critical point. In general, it represents the difference of the tuning parameter and its critical value, e.g. a coupling strength in a lattice description or the reduced temperature at a thermal phase transition.

Eq. (41) shows that for small negative $\theta_{2}$ the term proportional to the cubic coupling $\propto g$ will vanish only very slowly. This behavior should be clearly detectable in a Quantum Monte Carlo simulation. The corrections to scaling will become particularly important in the limits $N_{f}\gtrsim N_{f,c}$ and $N_{f}\to\infty$ as in both cases $\theta_{2}\to 0^{-}$.

IV.6.2 Emergent SUSY fixed point

To obtain an estimate of the accuracy of our results, we compare them to several other methods. In particular, there exist exact results in the limit where we set $N_{f}=1/2$ and restrict the symmetry to $U(1)$, which is the emergent symmetry at the considered fixed point. In that case, the model under investigation exhibits the same universal critical behavior as the Lagrangian describing the semimetal-superconductor quantum phase transition of a two-component Dirac fermion on the 2D surface of a 3D topological insulatorhasan2010 ; qi2011 reading Roy:2012wz ; Zerf2016

[TABLE]

where the gamma matrices are defined using the Pauli matrices, i.e. $\gamma_{0}=\sigma_{3},\gamma_{1}=\sigma_{1},\gamma_{2}=\sigma_{2}$, and $\bar{\psi}=-i\psi^{\dagger}\gamma_{0}$. For the model in Eq. (42), an emergent SUSY scenario described by Wess-Zumino theorywess1974 at the critical point was discussedlee2007 . Quantitative estimates for the critical exponents have been suggested by means of the conformal bootstrapBobev:2015vsa and within the $4-\epsilon$ expansion up to three loop orderZerf2016 . We present the corresponding comparison of the correlation length exponent and the anomalous dimensions in Tab. 2. Our results for the correlation length exponent lies within 5% of the conformal bootstrap and the three-loop calculations. The anomalous dimensions agree within 6% to the exact result of $\eta_{\psi}=\eta_{\phi}=1/3$. Let us also note a technical aspect here. Small values of $N_{f}$ reduce the impact of fermion fluctuations, which will shift the fixed point value of the $\left|\phi\right|^{2}$ coupling to negative values. This occurs around $N_{f}\approx 0.6$ in our calculations. In this case one can adapt the expansion of the effective potential to the new minimum appearing away from $\rho=\tau=0$ to obtain better estimates of critical exponents. We refrain from doing this here, because it only affects the $N_{f}=1/2$ case, where we nevertheless find reasonable critical exponents as demonstrated above.

IV.6.3 Comparison with MQMC

For larger fermion flavor numbers, we can compare our results with quantum Monte Carlo calculationsli2015 in $D=2+1$. We list the correlation length exponent and the boson anomalous dimension in Tab. 3. In addition, we quote the results of a first order expansion around the upper critical dimension. Our anomalous dimensions lie within 19% for $N_{f}=2$ to 8% for $N_{f}=6$ of the MC results. In this respect, it is known, however, that the FRG at this truncation order does not lead to very accurate boson anomalous dimensions. In contrast, we find that the correlation length exponents of our FRG calculation and the QMC result agree within 11% for $N_{f}=2$, which is improved to 2% agreement for $N_{f}=6$. Let us note here, that our prediction of substantial corrections to scaling have not yet been considered in the QMC simulation so that the agreement between the critical exponents may even improve.

In conclusion of the comparison, we believe that we can give good estimates for the critical behavior for a broad range of $D$ and $N_{f}$. We reproduce the perturbative RG results close to $D=4$. Furthermore, our critical exponents are comparable to available quantitative approaches in $D=3$. Thus we are confident that our analysis provides valid results for the whole range of $D$ and $N_{f}$, which has not been available to this extent before.

V Conclusions

In this work, we have studied the quantum phase transition to the Kekulé phase in two-dimensional Dirac materials. The condensation of the Kekulé order parameter reduces the “chiral” $U(1)$ symmetry of Dirac materials Eq. (5) to $\mathbb{Z}_{3}$. For electrons on the honeycomb lattice, it corresponds to a structural transition, where the translational invariance is broken because a valence bond solid is formed. The Kekulé VBS has already been observed in molecular graphenegomes2012 and graphene on a Copper substrateGutierrez2016 . With respect to graphene on a substrate, it is interesting to consider the effect of doping. A chemical potential introduces a new scale to the system and leads to a finite density of states, in favor of chiral symmetry breaking. It also preserves the chiral symmetry and the nesting close to the two inequivalent Dirac cones. But the nesting vector is not commensurate to the lattice anymore, which is disadvantageous for the Kekulé order. Furthermore, a finite chemical potential destroys Lorentz symmetry. However, universal critical behavior is determined by the scaling of the Fermi surface, not by its absolute magnitude, and the behavior at zero doping is continuously connected to the one at finite doping. Thus we expect that the infrared behavior of the weakly doped case is still governed by the same fixed point as the neutral system as long as long-range order can still develop.

Within a conventional Landau-Ginzburg description in terms of an order parameter field with small fluctuations, the corresponding phase transition would be of first order. The reason is that cubic terms of the order parameter field are allowed by symmetry and consequently appear in the free energy. Interestingly, the order of the phase transition can be changed when fluctuations become strong. This happens for example in the closely related three-states Potts model in 1+1 dimensionsbaxter1973 ; RevModPhys.54.235 , where the effect of fluctuations is increased due to the reduced dimensionality. Here, we have considered another intriguing mechanism to induce a second order transition, which is specific to the quantum character of the phase transition. At zero temperature, the Dirac fermions, which are massless in the symmetric phase, provide an additional critical mode at the phase transition. That means that fluctuations are increased due to the presence of further degrees of freedom, besides the order parameter field. We have investigated if this can indeed lead to a change of the nature of the Kekulé transition in Dirac materials.

We have modeled the system in terms of a Gross-Neveu-Yukawa theory with a generalized number of Dirac fermions $N_{f}$. Thereby, for example, $N_{f}=2$ describes graphene, and $N_{f}=1$ corresponds to spinless fermions on the honeycomb lattice. Further, we have discussed the emergent SUSY scenario for the case of $N_{f}=1/2$. With the help of the non-perturbative functional renormalization group, we have derived flow equations for the parameters of the model and investigated the fixed point structure. The fixed point corresponding to the Kekulé transition shows an emergent $U(1)$ symmetry and its stability depends on the number of Dirac fermions. Stability means the ability to tune to the critical point, which in turn decides if the transition is of first or second order. We found that the nature of the transition changes when $N_{f}\approx 1.9$. That means that the physically interesting case of graphene, or graphene-like systems, displays a second order phase transition. Hence, if the Kekulé transition becomes experimentally accessible, critical behavior should be observable.

The critical behavior at a second order phase transition is revealed by the scaling of correlation functions, characterized by critical exponents. We have calculated these critical exponents for the second order regime. Due to the proximity to the first order transition, the second largest critical exponent is, albeit irrelevant, close to zero. For large $N_{f}$, it further tends to exactly zero, i.e. the associated coupling, which corresponds to the cubic term of the order parameter field, is exactly marginal. We expect that this affects the scaling laws close to the phase transition, because they are modified by the scaling of the (almost) marginal coupling. The effect should be clearly visible in QMC simulations.

We also showed that for an accurate determination of the critical $N_{f}$ and the critical exponents, it is crucial to account for higher order couplings in the continuum description. In perturbative calculations, couplings of higher order than $\phi^{4}$ are usually omitted. But in 2+1 dimensions, couplings up to $\phi^{6}$ are relevant at the bare level. Although they turn irrelevant at the interacting fixed point, which describes the Kekulé transition, their scaling affects the critical exponents.

A further aspect of the critical behavior of the Kekulé transition is that through the breaking of the continuous $U(1)$ symmetry to $\mathbb{Z}_{3}$ a second scale is introduced. The reason is that the Goldstone mode of the $U(1)$ symmetry acquires a dynamical mass in the symmetry-broken phase, which is different from the longitudinal mass. Hence, the scaling of correlations is different on the two sides of the phase transition Leonard:2015wyg ; 2016arXiv161007603J . Interestingly, the FRG allows to study the dynamical generation of a mass gap in the spontaneously broken regime. We leave it for future work to demonstrate the different scaling in the Kekulé and the Dirac semimetal phase.

Acknowledgements.

The authors are grateful to E. Torres and Andrés Goens for discussions. IFH is supported by the NSERC of Canada. This work has been supported by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center SFB 1238, TP C04. LC is supported by the U.S. Department of Energy (DOE), Division of Condensed Matter Physics and Materials Science, under Contract No. DE-AC02-98CH10886.

Appendix A FRG flow equations

In this appendix, we present the derivation and the full expressions for the FRG flow equations. That means in our truncation for the effective action Eq. (15), we obtain flow equations for $Z_{\psi,k},Z_{\phi,k},\bar{h}_{k}$ and $U_{k}$. We express them in terms of the corresponding dimensionless quantities defined in Eqs. (18)-(19). For convenience, here and in the following, we suppress the index $k$.

A.1 Effective potential

The flow of the effective potential is obtained by evaluating Eq. (13) for constant bosonic field $\phi$ and vanishing fermion field $\psi$. This allows us to present a closed form for the flow of the full effective potential within the given truncation. It reads

[TABLE]

with the volume element $v_{D}^{-1}=2^{D+1}\pi^{D/2}\Gamma[D/2]$. We have denoted the derivatives with respect to the invariants by $u^{(i,j)}=\frac{\partial^{i+j}}{\partial\rho^{i}\partial\tau^{j}}u$. We have furhter introduced the longitudinal and transverse masses $m^{2}_{L/T}$ and the threshold functions $l_{0}^{(B)}$ and $I_{0}^{(F)}$. The threshold functions involve the loop integrations and the corresponding regulator dependence. We give the exact expressions in the App. B. The longitudinal and transversal mass $m_{L}^{2}$ and $m_{T}^{2}$ are obtained by diagonalizing $(\Gamma_{k}^{(2),B}+R_{k}^{B})^{-1}$ and read $m_{L/T}^{2}=(u_{11}+u_{22})/2\pm\sqrt{4u_{11}^{2}+(u_{11}-u_{22})^{2}}/2$ as function of $\phi_{1},\phi_{2}$ with $u^{(i,j)}=\frac{\partial^{i}}{\partial\rho^{i}}\frac{\partial^{j}}{\partial\tau^{j}}u$. When expressed in terms of $\rho,\tau$ the expressions become

[TABLE]

As we have explained in the main text, we expand $u$ in powers of $\rho$ and $\tau$ about its scale-dependent minimum. In our case the minimum of the fixed-point potential lies at the origin $(\rho,\tau)=(0,0)$ (except for very small $N_{f}\lesssim 1/2$).

A.2 Yukawa coupling

To extract the flow equation of the Yukawa coupling we decompose the two-point function into its fluctuation dependent and independent parts $\Gamma_{k,0}^{(2)}=\Gamma_{k}^{(2)}|_{\Delta\phi=\psi=0},$ and $\Delta\Gamma_{k}^{(2)}=\Gamma_{k}^{(2)}-\Gamma_{k,0}^{(2)}\,$ with $\phi=\phi_{0}+\Delta\phi$. Then, we expand the Wetterich equation as

[TABLE]

where we have defined the scale derivative $\tilde{\partial}_{t}$ which acts only on the $t$-dependence of the regulator (see also Eq. (62)). The real and imaginary components of the bosonic field are also divided into their vacuum expectation value and a fluctuating part, $\phi_{1}=\bar{\phi}_{1}+\Delta\phi_{1}$ and $\phi_{2}=\bar{\phi}_{2}+\Delta\phi_{2}$. To project the flow equation onto the Yukawa coupling, we then use

[TABLE]

evaluated at $\bar{\psi}=\psi=0$, $\Delta\phi_{1}=\Delta\phi_{2}=0$ and $p^{\prime}=p=q=0$. The flow of the squared Yukawa coupling can then be expressed in terms of field expectations values $\bar{\phi}_{1},\bar{\phi}_{2}$. We rewrite them in terms of the invariants at the minimum of the effective potential $\rho_{0}=(\bar{\phi}_{1}^{2}+\bar{\phi}_{2}^{2})/2$, $\tau_{0}=(\bar{\phi}_{1}^{3}-3\bar{\phi}_{1}\bar{\phi}_{2}^{2})/\sqrt{2}$. Therefore we choose our coordinate system so that the minimum lies on the $\phi_{1}$ axis, i.e. $\bar{\phi}_{2}=0$. Then, if we denote the value of $\rho$ at the minimum as $\rho_{0}=\kappa$, we have the relation $\kappa=\bar{\phi}_{1}^{2}/2$ and $\tau_{0}=\sqrt{2\kappa^{3}}$. With this the $\beta$-function of the squared Yukawa coupling becomes

[TABLE]

with \omega_{ijk}=\frac{\partial^{3}}{\partial\phi_{i}\partial\phi_{j}\partial\phi_{k}}u\big{|}_{min} and m_{L/T,0}^{2}=m_{L/T}^{2}\big{|}_{\text{min}} evaluated at the minimum of the potential. They read

[TABLE]

In our case these expressions simplify even further, because $\kappa=0$ in the SYM regime. The threshold functions $l_{ij}^{(FB)}$ required for the Yukawa coupling are also displayed in App. B.

A.3 Anomalous dimensions

To determine the expressions for the anomalous dimensions of the boson and fermion fields, we first need appropriate projection prescriptions for the wave function renormalizations. To that end, we evaluate the Wetterich equation for momentum-dependent fields and choose the following prescriptions

[TABLE]

with short-hand notation $\int_{q}=\int\frac{d^{D}q}{(2\pi)^{D}}$. The evaluation of these prescriptions yields a fermion anomalous dimension reading

[TABLE]

Similarly, the boson anomalous dimension is given by:

[TABLE]

and the threshold functions $m_{12}^{(FB)},m_{22}^{(B)},m_{4}^{(F)}$ and $m_{2}^{(F)}$ are again listed in App. B. This completes our list of FRG equations.

A.4 Flow equations in the symmetric regime

In the main text, we have stated the above expressions in the symmetric regime, where they simplify significantly. To obtain them, we have set $\kappa=0$ and $m_{L,0}^{2}=m_{T,0}^{2}=u^{(1,0)}=m^{2}$ and used the expressions for the threshold functions with the linear cutoff as given in the next section.

Appendix B Threshold functions

In the above flow equations, we have abbreviated the result of the loop integrations in terms of threshold functions. These functions also contain the regulator dependence of the FRG equations. They are defined as

[TABLE]

where we used regulators of the form $R_{\phi}=Z_{\phi}p^{2}r_{\phi}$ and $R_{\psi}=Z_{\psi}i\gamma_{\mu}p_{\mu}r_{\psi}$ and we defined

[TABLE]

To calculate the functional renormalization group flow equations for the running couplings/functionals, we further have to specify the regulator functions $r_{\psi}$ and $r_{\phi}$. Here, we choose linear regulators, which allow for a fully analytical calculation of the beta functions, i.e.

[TABLE]

with the Heavyside step function $\Theta(x)$. Furthermore they optimize the convergence to the physical solution litim1 ; litim2 ; litim3 ; litim4 . Then the threshold functions read

[TABLE]