Fermionic multicriticality near Kekul\'{e} valence-bond ordering in honeycomb lattice

TL;DR

This paper investigates the stability of emergent quantum multicritical points with enlarged symmetry in honeycomb lattice systems, showing that certain cubic terms become irrelevant near two dimensions, leading to stable MCPs with enhanced symmetry.

Contribution

It demonstrates that quantum fluctuations render cubic terms irrelevant at MCPs with $O(2+N_b)$ symmetry, establishing their stability in honeycomb lattice models.

Findings

Cubic terms are relevant near three dimensions but become irrelevant near two dimensions.

Enlarged $O(2+N_b)$ symmetric MCPs are stable due to Yukawa interactions.

Quantum fluctuations suppress discrete symmetry breaking at MCPs.

Abstract

We analyze the possibility of emergent quantum multicritical points (MCPs) with enlarged chiral symmetry, when strongly interacting gapless Dirac fermions acquire comparable propensity toward the nucleation of Kekul\'{e} valence-bond solid (KVBS) and charge-density-wave ($N_b=1$) or $s$-wave pairing ($N_b=2$) or anti-ferromagnet ($N_b=3$) in honeycomb lattice, where $N_b$ counts the number of bosonic order parameter components. Besides the cubic terms present in the order parameter description of KVBS due to the breaking of a discrete $Z_3$ symmetry, quantum fluctuations generate new cubic vertices near the high symmetry MCPs. All cubic terms are strongly relevant at the bare level near three spatial dimensions, about which we perform a leading order renormalization group analysis of coupled Gross-Neveu-Yukawa field theory. We show that due to non-trivial Yukawa interactions among gapless bosonic and fermionic degrees of freedom, all cubic terms ultimately become irrelevant at an $O(2+N_b)$ symmetric MCP, near two spatial dimensions, where $N_b=1,2,3$. Therefore, MCPs with an enlarged $O(2+N_b)$ symmetry near KVBS ordering are stable.

Equations

$$\displaystyle L_{f}=\Psi^{\dagger}(\tau,{\bf r})\bigg{[}\partial_{\tau}-i\sum^{d}_{j=1}\Gamma_{j}\partial_{j}\bigg{]}\Psi(\tau,{\bf r}),$$

$$L_{bf}=\sum _{j=1}^2{g}_{{}_j}{\Phi }_j{\Psi }^{†}{M}_j\Psi +g_{{}_3}\sum _{k=1}^{N_b}{\chi }_{{}_k}{\Psi }^{†}{M}_{2+k}\Psi ,$$

$$L_{\mathrm{Kek}}$$

$$+\sum _{j=1}^2\frac{{\lambda }_j}{4!}{\Phi }_j^4+\frac{2{\lambda }_{12}}{4!}{\Phi }_1^2{\Phi }_2^2,$$

$$L_{\mathrm{Sym}}^{N_b}$$

$$L_{\mathrm{Mix}}$$

$$\frac{d{g}_{{}_1}^2}{d\ell }$$

$$\frac{d{g}_{{}_2}^2}{d\ell }$$

$$\frac{d{u}_{{}_1}}{d\ell }$$

$$\frac{d{u}_{{}_2}}{d\ell }$$

$$\frac{d{u}_{{}_3}}{d\ell }$$

$$\frac{d{\lambda }_3}{d\ell }$$

$$\frac{d{\lambda }_1}{d\ell }$$

$$-12\left(u_{{}_1}^4+u_{{}_2}^4+N_b{u}_{{}_3}^4\right),$$

$$\frac{d{\lambda }_2}{d\ell }$$

$$\frac{d{\lambda }_{12}}{d\ell }$$

$$+\frac{1}{6}\left(11u_{{}_1}^2+33u_{{}_2}^2+24u_{{}_1}{u}_{{}_2}-N_b{u}_{{}_3}^2\right){\lambda }_{12}+6u_{{}_2}^2\left({\lambda }_1+{\lambda }_2\right)+2N_b{u}_{{}_3}^2{\lambda }_{23}-24\left(u_{{}_2}^4+u_{{}_1}^2{u}_{{}_2}^2+u_{{}_1}{u}_{{}_2}^3\right),$$

$$\frac{d{\lambda }_{13}}{d\ell }$$

$$+2\left[\left(N_b+2\right)u_{{}_3}^2{\lambda }_3+u_{{}_1}^2{\lambda }_{13}+u_{{}_2}^2{\lambda }_{23}+3u_{{}_3}^2\left({\lambda }_1+{\lambda }_{13}\right)+2u_{{}_1}{u}_{{}_3}{\lambda }_{13}\right]-24\left(u_{{}_1}^2{u}_{{}_3}^2+u_{{}_1}{u}_{{}_3}^3+u_{{}_3}^4\right),$$

$$\frac{d{\lambda }_{23}}{d\ell }$$

$$+\left[2\left(u_{{}_2}^2{\lambda }_{13}+u_{{}_3}^2{\lambda }_{12}\right)+3{\left(u_{{}_2}+u_{{}_3}\right)}^2{\lambda }_{23}\right]-24u_{{}_2}^2{u}_{{}_3}^2,$$

$$\left(g_{\ast }^2,{\lambda }_{\ast },u_{\ast }\right)=\left(ϵ,3\frac{H(N_f,N_b)}{10+N_b}ϵ,0\right)\frac{1}{2(N_f+1)-N_b},$$

$$M_{ij}\left(\{C\}\right)=\frac{d}{d{C}_j}\left(\frac{d{C}_i}{d\ell }\right),$$

$$\nu =\frac{1}{2}+\frac{N_f}{2}{g}_{\ast }^2+\frac{4+N_b}{24}{\lambda }_{\ast }.$$

Full text

Introduction: Gapless Dirac fermions constitute an ideal arena to explore the effects of electronic interactions and emergent quantum critical phenomena of itinerant systems. Typically at low enough temperatures and for sufficiently strong interactions, nodal Dirac fermions become susceptible toward a gap opening (mass generation), leading to a maximal gain of the condensation energy. The effective field theory describing the associated quantum phase transition assumes the form of a Gross-Neveu-Yukawa (GNY) model, which besides capturing the dynamics of fermionic and bosonic (order-parameter) fields, also accounts for the Yukawa coupling between them zinn-justin . Traditionally, the GNY theory is analyzed using a perturbative $\epsilon$ expansion, controlled by a parameter $\epsilon=3-d$, measuring the deviation from the upper critical three spatial dimensions. Close to the quantum phase transitions, the notion of sharp fermionic or bosonic excitations becomes moot, and the system accommodates a strongly coupled relativistic (due to an emergent Lorentz symmetry anber-lorentz ; roy-lorentz ; roy-kennett-yang-juricic ) ‘soup’ of these degrees of freedom, constituting a non-Fermi liquid.

Such rich field theory predictions recently became relevant in the context of condensed matter physics due to the possible realization of symmetry protected emergent Dirac excitations from concrete lattice models as, for example, in honeycomb graphene-review and $\pi$-flux square marston-affleck lattices. The associated quantum phase transitions are succinctly captured by minimal Hubbardlike models, containing only the finite-range components of the Coulomb interaction herbut-solo ; raghu-honerkamp ; honerpkapm ; herbut-juricic-roy ; gonzalez ; dagofer-hohendler ; juricic-roy-TBLG . Simplicity of these lattice models also permits numerical demonstration of quantum criticality using, for example, quantum Monte Carlo simulations sorella-1 ; herbut-assaad-1 ; chandrasekharan ; troyer-honeycomb ; herbut-assaad-2 ; sorella-2 ; hong-yao-NN-honeycomb ; kaul-itinerant , besides field theoretic analyses rosenstein ; herbut-juricic-vafek ; herbut-juricic-roy-SC ; roy-yang ; sslee ; hong-yao-1 ; machiejko-zarf ; klebanov ; knorr ; roy-juricic-PRL . An intriguing outcome in this context is the following. When Dirac quasiparticles acquire comparable propensity toward the formation of more than one (typically two) competing phases, respectively breaking O$(N_{1})$ and O$(N_{2})$ symmetries, such that the corresponding order-parameters can be rotated into each other by the generators of the emergent chiral symmetry, then it is conceivable to find a stable itinerant multicritical point (MCP), possessing an enlarged O$(N)$ symmetry, where $N=N_{1}+N_{2}$, but $N\leq 5$ in graphenelike systems, described by an eight-component Dirac spinor roy-MCP ; royjuricic-MCP ; herbut-MCP ; roygoswmaijericic-MCP . In particular, emergent MCPs with O$(N)$ symmetry with $3\leq N\leq 5$ lack any analogue in pure bosonic systems calabrese and can only be found in strongly interacting Dirac materials.

For the honeycomb lattice model, this observation supports the scenarios with the MCPs of following symmetries: (a) O(3), arising from the competition between Kekulé valence bond solid (KVBS) or $s$-wave pairing and charge-density-wave, (b) O(4), where KVBS and $s$-wave pairing form a chiral supervector, and (c) O(5), involving anti-ferromagnet and KVBS or quantum spin-Hall insulator and $s$-wave pairing. But, the emergent symmetry at an MCP becomes a subtle issue when one of the competing phases is the discrete three-fold rotational ($Z_{3}$) symmetry breaking KVBS Chamon-Kekule ; Chamon_SO5 ; roy-herbut-Kekule , since the corresponding cubic term in the bosonic order-parameter theory is strongly relevant near $d=3$ KekuleZ3_1 ; KekuleZ3_2 and responsible for first-order transitions Golner ; RMP-Potts . Therefore, question arises regarding the ultimate stability of the MCPs possessing a seemingly putative enlarged O$(2+N_{b})$ symmetry, with $N_{b}=1,2,3$, respectively, for charge-density-wave, $s$-wave pairing and anti-ferromagnet, in the vicinity of a KVBS ordering. This is the question we address in this work, and arrive at the following conclusions.

We show that when KVBS order, described by a two-component order parameter ${\bm{\Phi}}=(\Phi_{1},\Phi_{2})$, is coupled to an O$(N_{b})$ symmetry breaking field (${\bm{\chi}}$), a new cubic vertex $\Phi_{1}|{\bm{\chi}}|^{2}$ gets generated near a MCP through quantum corrections, besides $\Phi^{3}_{1}$ and $\Phi_{1}\Phi^{2}_{2}$ vertices, see Fig. 1. However, due to the Yukawa coupling between gapless fermions and bosons, all the cubic terms ultimately become irrelevant close to the O$(2+N_{b})$ symmetric MCPs at least near $d=2$, see Fig. 2. Therefore, MCPs with an enlarged O$(2+N_{b})$ symmetry near KVBS ordering are stable, and control continuous quantum phase transitions from an interacting Dirac liquid to (1) O$(N_{b})$ and (2) O$(N_{b}+2)$ symmetry breaking insulators, (3) KVBS, as well as (4) the direct transition between KVBS and O$(N_{b})$ symmetry breaking insulator.

Effective field theory: We begin the discussion by introducing the GNY model in the presence of all symmetry allowed cubic terms close to an $O(2+N_{b})$ symmetric quantum MCP. The corresponding imaginary time ($\tau$) action reads as ${\mathcal{S}}=\int d\tau d^{d}{\bm{r}}L$, with $L=L_{f}+L_{bf}+L_{b}$ and ${\bm{r}}=\left(x_{1},\cdots,x_{d}\right)$ is the spatial coordinate. The dynamics of massless Dirac fermions is captured by

[TABLE]

where $\Gamma_{j}$s are mutually anticommuting eight-dimensional Hermitian matrices, and $\Psi^{\dagger}$ and $\Psi$ are independent eight-component Grassmann variables. The following discussion is, however, impervious to specific matrix and spinor representations. The coupling between massless Dirac fermions and bosonic order parameter fields reads

[TABLE]

where $X\equiv X(\tau,{\bf r})$ for $X=\Psi^{\dagger},\Psi,{\bm{\Phi}},{\bm{\chi}}$. Mutually anticommuting eight-dimensional Hermitian mass matrices $M_{j}$s, satisfy $\{\Gamma_{i},M_{j}\}=0$, with $N_{b}\leq 5-d$. The $O(2+N_{b})$ chiral rotations among the mass matrices are generated by $G_{ij}=[M_{i},M_{j}]/(2i)$, where $i,j=1,\cdots,N_{b}+2$. Note $G_{ij}$s are the generators of the chiral symmetry for massless Dirac fermions as $[\Gamma_{i},G_{jk}]=0$.

The pure bosonic part of theory can be decomposed according to $L_{b}=L_{\rm Kek}+L^{N_{b}}_{\rm Sym}+L_{\rm Mix}$, with

[TABLE]

Therefore, the effective field theory contains 12 coupling constants. Even though three coupling constants appearing in $L_{\rm Mix}$ are absent at the bare level, $u_{{}_{3}}$, $\lambda_{13}$ and $\lambda_{23}$ get generated through quantum corrections, see Fig. 1, and thus have to be included from the outset.

Near a pure KVBS ordering $L^{N_{b}}_{\rm Sym}=L_{\rm Mix}=0$ and $g_{{}_{3}}=0$. In addition, $u_{1}=-u_{2}$, $g_{{}_{1}}=g_{{}_{2}}$ and $\lambda_{1}=\lambda_{2}=\lambda_{12}$. However, such a symmetry is broken, when the system acquires comparable propensity toward the formation of an O$(N_{b})$ symmetry breaking phase. Nevertheless, the above effective field theory enjoys a hidden O$(N_{b}+1)$ symmetry, which we identify by constructing a composite bosonic field according to $(\Phi_{2},{\bm{\chi}})$, and setting $g_{{}_{2}}=g_{{}_{3}}$, $u_{{}_{2}}=u_{{}_{3}}$, $\lambda_{2}=\lambda_{3}=\lambda_{23}$ and $\lambda_{12}=\lambda_{13}$. The effective field theory then describes a $Z_{2}\otimes{\rm O}(N_{b}+1)$ symmetric GNY model for massless Dirac fermions in the presence of cubic terms and contains 7 coupling constants. The renormalization group (RG) flow equations also reflect this symmetry, about which more in a moment. Anticipating the outcome, we set the Fermi ($v_{{}_{F}}$) and bosonic ($v_{{}_{B}}$) velocities to be equal (due to an emergent Lorentz symmetry at MCP) roy-lorentz , and $v_{{}_{F}}=v_{{}_{B}}=1$ for simplicity.

$\epsilon$ expansion: The bare (engineering) scaling dimension of the fermionic and bosonic fields, respectively, $D[\Psi]=d/2$ and $D[{\bm{\Phi}}]=D[{\bm{\chi}}]=(d-1)/2$, follow from the corresponding noninteracting parts of the action, the scaling dimension of momentum $D[k]=1$ and imaginary time $D[\tau]=-1$. Consequently, the scaling dimension of all three Yukawa couplings and six four-boson couplings is $D[\{g^{2}\}]=D[\{\lambda\}]=3-d$, and that for the bosonic mass parameters (controlling the transition from a Dirac semimetal to ordered phases) $D[\{m^{2}\}]=2$ footnote . Our focus here is restricted on the critical hyperplane, defined by $\{m^{2}\}=0$. Therefore, in the absence of any cubic terms, the GNY model near fermionic MCPs can be addressed in terms of a perturbative $\epsilon$ expansion about three spatial dimensions, with $\epsilon=3-d$ roy-MCP ; royjuricic-MCP ; herbut-MCP ; roygoswmaijericic-MCP .

On the other hand, the scaling dimension of all three cubic terms is $D[\{u\}]=(2+\epsilon)/2$, implying that they are strongly relevant perturbations at the bare level in $d=3$ KekuleZ3_1 ; KekuleZ3_2 . Therefore, the entire theory cannot be controlled in terms of a single expansion parameter $\epsilon$. Alternatively, one may pursue a double $\epsilon$ expansion, with $D[\{g^{2}\}]=D[\{\lambda\}]=\epsilon_{1}$ and $D[\{u\}]=\epsilon_{2}/2$, where $\epsilon_{1}=3-d$ and $\epsilon_{2}=5-d$, somewhat similar to the ones employed in interacting and disordered systems boyanovsky-cardy ; roy-dassarma ; yerzhakov-maciejko . However, such double $\epsilon$ expansion does not yield any additional control over the perturbative analysis. So, we abandon it for rest of the discussion. Nonetheless, the cubic vertices receive perturbative corrections, and the question we seek to answer in rest of the paper is whether these corrections (in particular, the ones arising due to the Yukawa couplings) can provide sufficiently negative scaling dimensions to the cubic vertices to turn them irrelevant near the fermionic MCPs, located at $\{g^{2}\}\sim\epsilon$, $\{\lambda\}\sim\epsilon$ and $\{u\}=0$, at least when $\epsilon\to 1$ or as $d\to 2$.

RG analysis: To answer this question, we here restrict ourselves to the leading order $\epsilon$ expansion for $\{g\}$ and $\{\lambda\}$, and account for one-loop corrections to $\{u\}$. To this end we integrate over the Matsubara frequency $-\infty\leq\omega\leq\infty$ and the fast Fourier modes, living within a thin Wilsonian momentum shell $\Lambda e^{-\ell}\leq|{\bf k}|\leq\Lambda$. We perform the matrix algebra in $d=2$ and subsequently carry the momentum integral in dimension $d=3-\epsilon$. The resulting coupled RG flow equations are

[TABLE]

in terms of dimensionless coupling constants, defined as $X\Lambda^{-\epsilon}/(8\pi^{2})\to X$, where $X=\{g^{2}\},\{\lambda\}$ and $\{u\}\Lambda^{-\frac{2+\epsilon}{2}}/(8\pi^{2})\to\{u\}$. Here, $N_{f}$ is the number of four-component fermion flavors, and hence for graphenelike systems $N_{f}=2$. The details of the RG calculation are presented in the Supplemental Materials supple . The underlying $Z_{2}\otimes{\rm O}(N_{b}+1)$ symmetry of the GNY theory can be appreciated in the following way. If we set $g_{{}_{2}}=g_{{}_{3}}$, $u_{{}_{2}}=u_{{}_{3}}$, $\lambda_{2}=\lambda_{3}=\lambda_{23}$ and $\lambda_{12}=\lambda_{13}$, then RG flow equations of the following couplings are identical: (1) $g_{{}_{2}}$ and $g_{{}_{3}}$, (2) $u_{{}_{2}}$ and $u_{{}_{3}}$, (3) $\lambda_{2}$, $\lambda_{3}$ and $\lambda_{23}$, and (4) $\lambda_{12}$ and $\lambda_{13}$. It is quite challenging to find all possible solutions of above 12 coupled RG flow equations, so we instead focus on the specific and relevant case, the O$(2+N_{b})$ symmetric MCP, located at $g^{2}_{{}_{1}}=g^{2}_{{}_{2}}=g^{2}_{{}_{3}}=g^{2}_{\ast}$, $\lambda_{1}=\lambda_{2}=\lambda_{12}=\lambda_{3}=\lambda_{13}=\lambda_{23}=\lambda_{\ast}$, and $u_{{}_{1}}=u_{{}_{2}}=u_{{}_{3}}=u_{\ast}$, where

[TABLE]

and $H(x,y)=2-2x-y+[4x^{2}+(y-2)^{2}+4x(38+5x)]^{1/2}$.

To analyze the stability of such fixed points, we compute the stability matrix (${\bm{M}}$), defined as

[TABLE]

and its eigenvalues in its vicinity. Here $\{C\}$ is the set of 12 coupling constants, and thus $i,j=1,\cdots,12$. The results are displayed in Fig. 2. Sufficiently close to two spatial dimensions (as $\epsilon\to 1$), all 12 eigenvalues of the stability matrix are negative for any value of $N_{f}$ and $N_{b}$. Hence, all cubic terms (namely, $u_{{}_{1}},u_{{}_{2}}$ and $u_{{}_{3}}$) become irrelevant in the close vicinity of the O$(2+N_{b})$ symmetric MCPs, indicating their stability in two dimensions. Note that irrelevance of the cubic terms is solely introduced by nontrivial Yukawa coupling between gapless bosonic and fermionic degrees of freedom. Hence, such a quantum MCP can only be realized in strongly interacting Dirac systems. We also note that with increasing number of the order-parameter components $N_{b}$, the regime of stability (the shaded region in Fig. 2) of the MCPs increases footnote_2d . The fact that the MCPs are stable over a range of $\epsilon$ (see Fig. 2), suggests that their stability in the presence of cubic couplings is possibly nonperturbative in nature. Only the range of $\epsilon$, over which the MCPs are stable, can be renormalized at each order in a perturbative expansion.

We should also mention that a leading order $\epsilon$ expansion in a purely bosonic system suggests a putative O(3) symmetric MCP, which, however, looses stability once the higher order corrections are taken into account calabrese . On the other hand, there exists neither O(4) nor O(5) symmetric MCP in a purely bosonic theory. Therefore, appearance of O$(2+N_{b})$ symmetric quantum MCPs, with $N_{b}=1,2,3$, and their stability against the cubic perturbations, when the system resides at the brink of Kekulé ordering in honeycomb lattice, are purely fermion driven phenomena in Dirac materials. Also a pure Kekulé O$(2)$ quantum critical point (in the absence of the ${\bm{\chi}}$ field) is stabilized due to gapless Dirac fermions KekuleZ3_1 ; KekuleZ3_2 ; comment-Kekule ; footnote_rosenstein ; footnote_2d .

Besides the O$(2+N_{b})$ symmetric MCP, there exists two more interacting fixed points, possessing $O(2)$ and $O(N_{b})$ symmetries. Respectively, they control transition from a Dirac semimetal to KVBS and O$(N_{b})$ symmetry breaking phase. The cubic terms are irrelevant at these two fixed points. However, they ultimately become unstable toward the MCP. We could not find (numerically) any fixed point at finite $\{u_{i}\}$, when $\epsilon$ is close to 1.

Discussion: To summarize, we address the stability of quantum MCPs with enlarged $O(2+N_{b})$ symmetry, when a correlated Dirac liquid, realized on a honeycomb lattice, acquires comparable propensity toward the nucleation of KVBS and charge-density-wave ($N_{b}=1$) or $s$-wave pairing ($N_{b}=2$) or antiferromagnet ($N_{b}=3$). We show that quantum corrections generate new cubic vertices (see Fig. 1) near such MCPs, besides the ones for KVBS due to the breaking of discrete $Z_{3}$ rotational symmetry in honeycomb lattice. All cubic terms are strongly relevant at the bare level and responsible for a generic first-order transition in pure bosonic systems Golner ; RMP-Potts . However, due to nontrivial boson-fermion Yukawa couplings, they all become irrelevant near high symmetric MCPs, close to two spatial dimensions, see Fig. 2. Therefore, $O(2+N_{b})$ symmetric quantum MCPs are expected to be stable in honeycomb Dirac systems. At this MCP the bosonic and fermionic anomalous dimensions are respectively given by $\eta_{b}=2N_{f}g^{2}_{\ast}$ and $\eta_{\Psi}=(2+N_{b})g^{2}_{\ast}/2$, and the correlation length exponent is

[TABLE]

Together they determine the universality class of continuous transitions from a Dirac semimetal to (a) KVBS, (b) $O(N_{b})$ and (c) $O(2+N_{b})$ symmetry breaking orders, and (d) the direct transition from KVBS to an $O(N_{b})$ symmetry breaking order through itinerant MCP.

Our results should be germane in the context of recent quantum Monte Carlo simulations on correlated Dirac liquid, in the presence of competing orderings sato-hohendler-assaad ; HongYao-2019 . So far, emergence of high symmetry has only been reported slightly away from the itinerant MCP. We hope that our results will motivate future works to explore symmetry enlargement in the proximity to the MCP and KVBS ordering. In addition, our results can also be relevant for slow (due to sufficiently small Fermi velocity) and strongly interacting (due to substantial bandwidth suppression) Dirac fermions in twisted bilayer graphene near magic angle (MA-TBLG) TBLG-1 ; TBLG-2 , and in correlated organic Dirac material pressured $\alpha$-(BEDT-TTF)2I3 organicDirac . A recent experiment reported the existence of an insulating phase near the charge-neutrality point in the former system TBLG-3 . With the application of suitable non-thermal tuning parameters (for example, pressure, strain, twist angle etc.) it is, at least in principle, conceivable to drive such a Dirac insulator through a quantum MCP, where the present discussion can be pertinent Kekule-MA-TBLG .