# SO(4) multicriticality of two-dimensional Dirac fermions

Igor F. Herbut  Department of Physics, Simon Fraser University, Burnaby,
British Columbia, Canada V5A 1S6  Michael M. Scherer  Institute for
Theoretical Physics III, Ruhr-University Bochum, D-44801 Bochum, Germany

###### Abstract

We study quantum multicritical behavior in a (2+1)-dimensional Gross-Neveu-
Yukawa field theory with eight-component Dirac fermions coupled to two
triplets of order parameters that act as Dirac masses, and transform as
$(1,0)+(0,1)$ representation under the SO(4)$\simeq$ SO(3)$\times$SO(3)
symmetry group. This field theory is relevant to spin-1/2 fermions on
honeycomb or $\pi$-flux lattices, for example, near the transition points
between an $s$-wave superconductor and a charge-density wave, on one side, and
Néel order, on the other. Two triplets of such order parameters always allow
for a common pair of two other order parameters that would complete them to
the maximal set of compatible (anticommuting) orders of five. We first derive
a unitary transformation in the Nambu (particle-hole) space which maps any two
such triplets, possibly containing some superconducting orders, onto purely
insulating order parameters. This allows one to consider a universal SO(4)
Gross-Neveu-Yukawa description of the multicriticality without any Nambu
doubling. We then proceed to derive the renormalization-group flow of the
coupling constants at one-loop order in $4-\epsilon$ space-time dimensions,
allowing also a more general set of order parameters transforming under
SO($n_{a}$)$\times$SO($n_{b}$). While for $n_{a}=n_{b}>2$ in the bosonic
sector and with fermions decoupled there is a stable fixed point of the flow,
the Yukawa coupling to fermions quickly leads to its elimination by a generic
fixed-point collision in the relevant range of fermion flavor numbers $N_{f}$.
This suggests the replacement of the critical behavior by a runaway flow in
the physical case $n_{a}=n_{b}=3$. The structure of the RG flow at $n_{a}\neq
n_{b}$ is also discussed, and some non-perturbative arguments in favor of the
stability of the decoupled critical point when $n_{a}=3$ and $n_{b}=1$ in
$D=2+1$ are provided.

##  I Introduction

Two-dimensional quasi-relativistic Dirac electrons appear as low-energy
excitations in many electronic systems, graphene arguably being the best-known
example. Interactions between electrons can lead to various phase transitions
at strong coupling, with the concomitant quantum critical behavior defining
new generalized Gross-Neveu universality classes prl06 ; juricic09 ; vafek09 .
The latter have been a subject of numerous analytic fritz ; mesterhazy ;
janssen12 ; roy13 ; janssen14 ; knorr16 ; scherer16 ; mihaila17 ; li17 ;
classen17 ; zerf ; knorr18 ; torres18 ; wamer ; ihrig18 ; gracey18 ; yin18 ;
liu21 ; yerzhakov21 ; boyack21  and numerical studies sorella92 ; sorella12 ;
assaad ; toldin ; otsuka16 ; buividovich ; hesselmann ; tang ; lang19 ; li20 ;
ostmeyer20 ; otsuka20 ; huffman20 ; ostmeyer21 ; mondaini22 . Also, the
multicritical behavior near points of the phase diagram where three (or more)
phases meet has attracted considerable interest roy11 ; classen ; sato2017 ;
janssen ; roy18 ; torres2020 . Recent quantum Monte Carlo simulations liu ,
for example, have explored a situation where $N_{f}=2$ four-component Dirac
fermions are coupled to two order parameters, each appearing as three
relativistic mutually anticommuting masses in the Dirac Hamiltonian, and
transforming as a vector under an SO(3) symmetry group: the first internal
SO(3) rotates between the charge-density-wave (CDW) and two components of the
$s$-wave superconductor ($s$SC1, $s$SC2), and the second SO(3) rotates the
three spatial components of the Néel order (Néel1, Néel2, Néel3).

Here, we further investigate this multicriticality within a field-theoretic
approach. We extend the previous analysis to a generalized Dirac system with
coupled order parameters that transform under SO($n_{a}$)$\times$SO($n_{b}$).
For the specific value of $n_{a}=n_{b}=3$, we come to a different conclusion
about the nature of multicritical behavior than Ref. liu, .

It is useful to first recall some basic facts about the mass-like order
parameters for the two-dimensional Dirac fermions in question. Two-dimensional
spin-1/2 electrons on honeycomb lattice, to take a specific example, are
described by eight-component Dirac fermions, due to having two spin, two
valley, and two sublattice degrees of freedom. Its standard relativistic Dirac
Hamiltonian therefore features an SU(4) symmetry that unifies and extends the
spin and the valley (or chiral) symmetries, which by themselves would only
yield the SO(3)$\times$SO(3)$\simeq$ SO(4)$\subset\,$SU(4) group of symmetry.
One pair of triplets of Dirac masses (which throughout this paper will be used
interchangeably with order parameters) that would transform as $(1,0)+(0,1)$
representation of this spin-valley SO(4) subgroup is the three components of
the quantum spin Hall (QSH1, QSH2, QSH3) insulator as the first, and the CDW
and the two Kekulé bond-density-wave (BDW1, BDW2) insulators as the second
triplet.

The inclusion of superconducting mass order parameters requires, however, the
usual Nambu-doubling of Dirac Hamiltonian into a sixteen-dimensional particle-
hole-symmetric Dirac Bogoliubov-de Gennes (BdG) Hamiltonian. The crucial
observation then is that the maximal number of “compatible” order parameters,
insulating and/or superconducting, which are represented by mutually
anticommuting Dirac masses in the Nambu-doubled representation is five ryu .
This number is ultimately determined by the antiunitary particle-hole symmetry
inherent to Nambu’s construction herbut12 . All 56 such quintuplets of
compatible Dirac masses have been listed with the help of computer algebra in
Ref. ryu, . Each triplet of mutually anticommuting masses can be found in only
one quintuplet, but any pair of mutually anticommuting masses, remarkably,
belongs to two different quintuplets herbut12 . This last feature will be of
crucial importance for our further discussion.

To quickly preview our results, we note that there are two main aspects of the
problem. First, we observe that any two triplets of Dirac masses that mutually
anticommute within the same triplet, but commute when from different triplets,
111That is the two triplets transform as $(1,0)+(0,1)$ representation of some
SO(3)$\times$ SO(3) subgroup of symmetry of the Dirac BdG Hamiltonian, can
always be complemented by the same pair of some other two Dirac masses to
their respective quintuplets. In other words, the only realization of the
$(1,0)+(0,1)$ mass order parameters is provided by two quintuplets that share
a common pair of masses. The two quintuplets (CDW, $s$SC1, $s$SC2, BDW1, BDW2)
and (Néel1, Néel2, Néel3, BDW1, BDW2) are one such example, which would
correspond to the two triplets studied in Ref. liu, . Further, (CDW, $s$SC1,
$s$SC2, BDW1, BDW2) and ($s$SC1, $s$SC2, QSH1, QSH2, QSH3) form another pair
of quintuplets that would correspond to our example above ghaemi ; herbut10 .
Building on this result, we construct a unitary transformation in the Nambu
space that changes the common pair of the two quintuplets, allowing us to
exchange the first pair of above quintuplets for the second. This is of
practical use, since the second pair of quintuplets has the two
superconducting orders in common, with all the other orders appearing being
insulators. One can therefore in complete generality consider a fully
insulating realization of the $(1,0)+(0,1)$ order parameters, and remain in
the original (pre-Nambu-doubled) eight-dimensional representation.

Second, multicritical behavior has been an important part of general studies
of critical phenomena from the beginning. Given two order parameters that
transform as vectors under SO($n_{a}$)$\times$SO($n_{b}$) symmetry, based on
high-order perturbative calculations and some non-perturbative arguments
calabrese , it is believed today that the bosonic field theory in three
dimensions has a decoupled critical point for $n_{a}+n_{b}>3$, an
SO($n_{a}+n_{b}$)-symmetric critical point for $n_{a}=n_{b}=1$, and a mixed-
symmetry (“biconal”) critical point for $n_{a}=1$ , $n_{b}=2$. To the first
order in $4-\epsilon$ dimension, however, the decoupled fixed point is stable
only when $n_{a}+n_{b}>8$, and the SO($n_{a}+n_{b}$)-symmetric fixed point is
stable for $n_{a}+n_{b}<4$ book . The main point is that there always exists a
stable (critical) fixed point in the purely bosonic theory for the order
parameters. This is in spite of the fixed points moving around with the change
of parameters such as $\epsilon$, $n_{a}$, and $n_{b}$, and colliding with
each other. The usual annihilation of fixed points and their transfer to the
complex plane, relevant to many other field theories halperin ; nauenberg ;
kubota ; herbut06 ; gies ; book1 ; kaplan ; herbut14 ; herbut16 ; gorbenko ;
ihrig19 ; faedo , does not occur here.

The reason for this lies in the structure of the theory itself: (a) when one
of the two quartic couplings vanishes and the order parameters decouple, they
remain decoupled at all stages of the renormalization group (RG)
transformation. This condition defines therefore an RG-invariant line, along
which there always exists a Wilson-Fisher fixed point. (b) When the two
quartic couplings have a particular ratio so that the theory acquires an
enlarged SO($n_{a}+n_{b}$) symmetry, it retains it during the RG, so another
Wilson-Fisher fixed point always exists on this invariant line as well. Due to
the existence of these two invariant lines with concomitant Wilson-Fisher
fixed points on them, the only change with varying the parameters such as the
dimension, $n_{a}$, and $n_{b}$ in the theory is the exchange of stability of
the two Wilson-Fisher fixed points with the third (biconal) fixed point, the
location of which is not fixed by any symmetry. This is in contrast to the
generic scenario; in general, when two fixed points collide they develop an
imaginary part, and in the space of real (physical) couplings there remains
only a runaway flow.book1  This is avoided in the
SO($n_{a}$)$\times$SO($n_{b}$)-symmetric bosonic theory only because of the
existence of two RG-invariant lines from which the Wilson-Fisher fixed points
cannot escape.

This crucial feature of the field theory is lost upon coupling to fermions.
Since the Dirac fermions couple to both order parameters, the integration over
any number of fermions unavoidably couples them. Furthermore, the
SO($n_{a}+n_{b}$) symmetry for special values of the two quartic couplings is
immediately violated as well. Hence, there is nothing to protect the fixed
points from complexification after a collision, and that is precisely what we
find. All four fixed points (including the Gaussian one) move when the
parameter $n_{a}=n_{b}$ is changed, and in particular the stable (biconal)
fixed point collides with the SO($n_{a}+n_{b}$) fixed point at fairly low
number of fermions.

The rest of the paper is organized as follows. In the next section we discuss
the algebra of Dirac mass order parameters, and perform the unitary
transformation onto insulators. In Sec. III we introduce the Gross-Neveu-
Yukawa field theory which features the two triplets of order parameters
coupled to Dirac fermions. The RG flow in this theory computed to one-loop
near four space-time dimensions is given in Sec. IV. The analysis of the fixed
points and their collisions with the change of number of Dirac fermions is
provided in Sec. V. Discussion of our results and the summary are given in
Sec. VI. Mathematical details of the Clifford-algebraic structure of the mass
order parameters of two-dimensional Dirac fermions and the alternative
derivation of the RG beta-functions are given in the Appendices.

##  II Algebra of Dirac masses

We are interested in the relativistically invariant mass-terms that can be
added to the standard two-dimensional Dirac Hamiltonian

$H_{0}(\vec{p})=p_{1}\alpha_{1}+p_{2}\alpha_{2},$ (1)

where $\alpha_{i}$, $i=1,2$ are two $8\times 8$ Hermitian matrices satisfying
the anticommutation relation $\\{\alpha_{i},\alpha_{j}\\}=2\delta_{ij}$. This
Hamiltonian represents the low-energy degrees of freedom of spin-1/2 electrons
near two inequivalent valleys on graphene’s honeycomb lattice, for example.

To have a unified treatment of both insulating (particle-number preserving)
and superconducting (particle-number violating) mass-like order parameters we
consider the usual Nambu-doubled Bogoliubov-de Gennes (BdG) Hamiltonian

$H_{\mathrm{BdG}}=H_{0}(\vec{p})\oplus(-H_{0}^{T}(-\vec{p}))=p_{1}\Gamma_{1}+p_{2}\Gamma_{2},$
(2)

where $\Gamma_{i}$ are now $16\times 16$ Hermitian matrices, which still
satisfy $\\{\Gamma_{i},\Gamma_{j}\\}=2\delta_{ij}$.

The algebraic structure of the mass-like orders follows from the following
observation: add two compatible (anticommuting) mass-terms of the above BdG
Hamiltonian, with two $16\times 16$ Hermitian matrices $M_{1}$ and $M_{2}$,
such that $\\{M_{i},M_{j}\\}=2\delta_{ij}$, and $\\{M_{i},\Gamma_{j}\\}=0$.
The Hamiltonian is then

$H=H_{\mathrm{BdG}}+m_{1}M_{1}+m_{2}M_{2},$ (3)

with the parameters (masses) $m_{1,2}$ real. The masses can even be taken to
be arbitrary functions of space. The group of symmetry of the above
Hamiltonian is SO(4)$\equiv$SO(3)$\times$SO(3), generated by six $16\times 16$
Hermitian generators $A_{i}$ and $B_{i}$, $i=1,2,3$, such that
$\\{A_{i},A_{j}\\}=\\{B_{i},B_{j}\\}=2\delta_{ij}$, and $[A_{i},B_{j}]=0$, cf.
Ref. herbut12, . In other words,

$[A_{i},\Gamma_{j}]=[A_{i},M_{j}]=0,$ (4)

and

$[B_{i},\Gamma_{j}]=[B_{i},M_{j}]=0,$ (5)

$i=1,2,3$, $j=1,2$. This result follows from the theory of real
representations of Clifford algebras okubo  and was derived in Ref. herbut12,
, where it was used to reveal the internal structure of the vortex
configuration in the masses $m_{1,2}$. It immediately follows that the two
mass-matrices $M_{1,2}$ can be embedded into the maximal set of five mutually
anticommuting mass-matrices herbut12  in two different ways: (1)
$M_{A}=(M_{1},M_{2},A_{i}X)$, and (2) $M_{B}=(M_{1},M_{2},B_{i}X)$, $i=1,2,3$.
Here, the matrix $X$ is defined as $X=\Gamma_{1}\Gamma_{2}M_{1}M_{2}$, so that
it anticommutes with the matrices present in $H$:
$\\{X,M_{i}\\}=\\{X,\Gamma_{i}\\}=0$. This in turn implies that the Clifford
algebra generated by the seven matrices $C_{A}=(\Gamma_{i},M_{A})$
($C_{B}=(\Gamma_{i},M_{B})$) is quaternionic okubo , i.e. that besides the
unity matrix it has three additional Casimir operators $B_{i}$ ($A_{i}$). More
details of this construction were given in Ref. herbut12, , and are here
provided in App. A.

For our current purposes the crucial fact is that the two sets of mass-
matrices $A_{i}X$ and $B_{i}X$, which share a common pair $(M_{1},M_{2})$ in
their respective quintuplets $M_{A}$ and $M_{B}$, commute, i.e.

$[A_{i}X,B_{j}X]=0,$ (6)

and transform as $(1,0)+(0,1)$ representation of the above group SO(4). The
inverse is also true: given two triplets of mutually anticommuting mass
matrices such that any two matrices from different triplets commute, they
always are completed by the same and unique pair of mass matrices to their
respective quintuplets. The above algebraic structure of the mass terms is
thus a universal characteristic of two-dimensional Dirac equation for spin-1/2
particles on a bipartite lattice.

An important corollary is that the six order parameters in the $(1,0)+(0,1)$
representation can always be assumed to be insulating. To show this, let us be
specific, and choose

$C_{A}=(\Gamma_{1},\Gamma_{2},K_{1},K_{2},C,S_{1},S_{2}),$ (7)

with $K_{1,2}$ being two matrices for the Kekulé bond-density-wave orders, $C$
for the charge-density-wave, and $S_{1,2}$ for the two components of the
$s$-wave superconducting order. The main assumption is that two ($S_{1}$ and
$S_{2}$) out of five mass matrices in the algebra are two components of the
same superconducting mass order parameter. The second algebra that contains
$(\Gamma_{1},\Gamma_{2},K_{1},K_{2})$ as a sub-algebra is then

$C_{B}=(\Gamma_{1},\Gamma_{2},K_{1},K_{2},N_{1},N_{2},N_{3}).$ (8)

We can identify one of the Casimir operators of $C_{B}$ as
$A_{1}=iS_{1}S_{2}$, which rotates between two components of the
superconducting order parameter, and thus is nothing but the particle number.
Therefore, the mass-matrices $N_{i}$, $i=1,2,3$ satisfy $[N_{i},A_{1}]=0$, and
represent insulators 222In fact $N_{i}$ happen to be the three components of
the Néel order parameter..

Let us now perform the unitary transformation using

$\displaystyle U=e^{i\frac{\pi}{4}\left[iK_{1}S_{1}+iK_{2}S_{2}\right]},$ (9)

which simply reshuffles the elements of the first Clifford algebra

$\displaystyle C_{A}\rightarrow
UC_{A}U^{\dagger}=(\Gamma_{1},\Gamma_{2},S_{1},S_{2},C,K_{1},K_{2}).$

The second Clifford algebra then transforms as

$\displaystyle C_{B}\rightarrow
UC_{B}U^{\dagger}=(\Gamma_{1},\Gamma_{2},S_{1},S_{2},Y_{1},Y_{2},Y_{3}),$

with some new mass-matrices $Y_{i}$, given by $Y_{i}=UN_{i}U^{\dagger}$. Since
$\\{Y_{i},S_{1}\\}=\\{Y_{i},S_{2}\\}=0$, it readily follows that

$\displaystyle[Y_{i},iS_{1}S_{2}]=0\,.$ (10)

Since the matrix $iS_{1}S_{2}$ is the particle number operator, however, the
three masses $Y_{i}$ evidently represent insulating order parameters 333They
are actually standing for the three spin-components of the quantum spin-Hall
insulator.. The upshot is that in the $16\times 16$ BdG representation,

$\displaystyle
U(H_{\mathrm{BdG}}+a_{1}C+a_{2}S_{1}+a_{3}S_{2}+b_{i}N_{i})U^{\dagger}$
$\displaystyle\quad\quad\quad\quad=H_{\mathrm{BdG}}+a_{1}C+a_{2}K_{1}+a_{3}K_{2}+b_{i}Y_{i},$
(11)

with all six order parameters on the right-hand-side of the last equation
representing insulators. The right-hand-side is then block-diagonal, and the
Nambu doubling can be disposed of. The general two-dimensional Dirac
Hamiltonian with two sets of mass-order-parameters that transform as
$(1,0)+(0,1)$ representation of the SO(4) can therefore simply be taken to be
$8\times 8$ matrix.

##  III Field theory

The above reasoning motivates one to consider a Lagrangian of the form

$\displaystyle\mathcal{L}=\mathcal{L}_{\psi}+\mathcal{L}_{ab}+\mathcal{L}_{\psi
ab}\,,$ (12)

where $\mathcal{L}_{\psi}$ is the kinetic term of the Dirac fermions,
$\mathcal{L}_{ab}$ includes all purely bosonic (order-parameter) terms, and
$\mathcal{L}_{\psi ab}$ represents the Yukawa interactions. Explicitly, the
Dirac Lagrangian $\mathcal{L}_{\psi}$ in general dimension is defined as

$\displaystyle\mathcal{L}_{\psi}=\psi^{\dagger}\partial_{\tau}\psi+\psi^{\dagger}\Big{[}\mathbb{I}_{4\times
4}\otimes\sum_{i=1}^{d}\alpha_{i}\big{(}-i\frac{\partial}{\partial
x_{i}}\big{)}\Big{]}\psi\,,$ (13)

where $\tau$ is the imaginary time, and $x_{i}$ the spatial coordinates.
$\psi=\psi(x_{i},\tau)$ is a $4\cdot 2^{d-1}$-component Grassmann field and
$d$ is the spatial dimension. For graphene, $d=2$. The $\alpha_{i},\
i\in\\{1,\ldots,d\\}$ are $2^{d-1}\times 2^{d-1}$-dimensional matrices, which
satisfy the Clifford algebra
$\left\\{\alpha_{i},\alpha_{j}\right\\}=2\delta_{ij}$. This setup allows for
an introduction of an additional matrix $\beta$ that obeys

$\displaystyle\left\\{\beta,\alpha_{i}\right\\}=0\
\forall\,i\,,\quad\text{and}\quad\beta^{2}=1\,,$ (14)

and we define all matrices $\alpha_{i},\beta$ as Hermitian. For example, an
explicit choice in $d=2$ may be
$\alpha_{1}=\sigma_{1},\alpha_{2}=\sigma_{2},\beta=\sigma_{3}$ where the
$\sigma_{i}$ are the conventional Pauli matrices. The bosonic part
$\mathcal{L}_{ab}$ includes two coupled and real-valued order parameters
$\vec{a}$ and $\vec{b}$ with components $a_{i},i=1,\ldots,n_{a}$ and
$b_{j},j=1,\ldots,n_{b}$ which transform as vectors under SO($n_{a}$) and
SO($n_{b}$), respectively. The Lagrangian $\mathcal{L}_{ab}$ then reads

$\displaystyle\mathcal{L}_{ab}=$
$\displaystyle\left[(\partial_{\mu}a_{i})(\partial^{\mu}a_{i})+(\partial_{\mu}b_{j})(\partial^{\mu}b_{j})\right]+r_{a}a_{i}a_{i}+r_{b}b_{j}b_{j}$
$\displaystyle+\lambda_{a}\left(a_{i}a_{i}\right)^{2}+\lambda_{b}\left(b_{j}b_{j}\right)^{2}+2\lambda_{ab}a_{i}a_{i}b_{j}b_{j}\,,$
(15)

where summation convention over repeated indices is adopted, and
$\mu\\!=\\!0,\ldots,d$. For general quadratic and quartic interaction
parameters $\\{r_{a},r_{b},\lambda_{a},\lambda_{b},\lambda_{ab}\\}$, the
Lagrangian $\mathcal{L}_{ab}$ has SO($n_{a}$)$\times$SO($n_{b}$) symmetry. For
$\lambda_{ab}=0$, the order parameters $\vec{a}$ and $\vec{b}$ are decoupled
in $\mathcal{L}_{ab}$, and for $r_{a}=r_{b}$ and
$\lambda_{a}=\lambda_{b}=\lambda_{ab}$, the Lagrangian $\mathcal{L}_{ab}$
acquires a larger SO($n_{a}+n_{b}$) symmetry.

For the remainder of this work, we restrict the number of bosonic field
components in each order parameter to maximally three, i.e. $n_{a},n_{b}\leq
3$, unless explicitly stated otherwise. We will also be particularly
interested in the case $n_{a}=n_{b}=3$.

Finally, we write down the Yukawa interaction

$\displaystyle\mathcal{L}_{\psi
ab}=\psi^{\dagger}\left[g_{a}(a_{i}\sigma_{i})\\!\otimes\\!\mathbbm{1}_{2}+g_{b}\mathbbm{1}_{2}\otimes(b_{j}\sigma_{j})\right]\otimes\beta\,\psi\,,$
(16)

introducing the Yukawa couplings $g_{a}$ and $g_{b}$. Importantly, as we will
show below, the full model, cf. Eq. (12), with maximal $n_{a}=n_{b}=3$, can be
merely SO(3)$\times$SO(3) invariant, i.e. there is no larger SO(6) at special
values of $\lambda_{a}$, $\lambda_{a}$, and $\lambda_{ab}$. Furthermore, the
order parameters $\vec{a},\vec{b}$ are always coupled, i.e. there is no
decoupling even for $\lambda_{ab}=0$.

##  IV Renormalization group flow

We perform a one-loop Wilsonian integration over the fermionic and bosonic
modes within the momentum shell from $\Lambda/b$ to $\Lambda$, and expand in
$\epsilon=4-D$ with $D=d+1$. This leads to the RG flow of the squared Yukawa
coupling constants

$\displaystyle\dot{g}_{a}^{2}$ $\displaystyle=\epsilon
g_{a}^{2}+\left(n_{a}-2N_{f}-4\right)g_{a}^{4}-3n_{b}g_{a}^{2}g_{b}^{2}\,,$
(17) $\displaystyle\dot{g}_{b}^{2}$ $\displaystyle=\epsilon
g_{b}^{2}+\left(n_{b}-2N_{f}-4)\right)g_{b}^{4}-3n_{a}g_{a}^{2}g_{b}^{2}\,,$
(18)

where $\dot{g}_{i}^{2}\equiv dg_{i}^{2}/d\ln b,i=a,b$, and we rescaled the
couplings as $(g_{i}^{2}/2)(\Lambda^{D-4}S_{D}/(2\pi)^{D})\rightarrow
g_{i}^{2}$. Here, $S_{D}$ is the area of the unit sphere in $D$ dimensions.
The flow of the quartic couplings is

$\displaystyle\dot{\lambda}_{a}$
$\displaystyle=\\!\epsilon\lambda_{a}\\!-\\!4N_{\\!f}g_{a}^{2}\lambda_{a}\\!-\\!4(n_{a}\\!+\\!8)\lambda_{a}^{2}\\!-\\!4n_{b}\lambda_{ab}^{2}\\!+\\!N_{\\!f}g_{a}^{4}\,,$
(19) $\displaystyle\dot{\lambda}_{b}$
$\displaystyle=\epsilon\lambda_{b}\\!-\\!4N_{\\!f}g_{b}^{2}\lambda_{b}\\!-\\!4(n_{b}\\!+\\!8)\lambda_{b}^{2}\\!-\\!4n_{a}\lambda_{ab}^{2}\\!+\\!N_{\\!f}g_{b}^{4}\,,$
(20) $\displaystyle\dot{\lambda}_{ab}$
$\displaystyle=\epsilon\lambda_{ab}\\!-\\!2N_{\\!f}\left(g_{a}^{2}\\!+\\!g_{b}^{2}\right)\lambda_{ab}\\!-\\!16\lambda_{ab}^{2}\\!-\\!4(n_{a}\\!+\\!2)\lambda_{a}\lambda_{ab}$
$\displaystyle\quad-4(n_{b}\\!+\\!2)\lambda_{b}\lambda_{ab}+3N_{\\!f}g_{a}^{2}g_{b}^{2}\,.$
(21)

with $(\lambda_{i}/4)(\Lambda^{D-4}S_{D}/(2\pi)^{D})\rightarrow\lambda_{i}$,
$i\in\\{a,b,ab\\}$. We have also generalized the model to include a number of
$N_{f}$ four-component Dirac fermions, with $N_{f}=2$ in graphene. The bosonic
and fermionic anomalous dimensions read as $\eta_{i}=2N_{f}g_{i}^{2}\,,\
i\in\\{a,b\\}$ and
$\eta_{\psi}=\frac{1}{2}n_{a}g_{a}^{2}+\frac{1}{2}n_{b}g_{b}^{2}$,
respectively.

We have checked that these flow equations are in agreement with Refs. classen,
and janssen,  in the respective limiting cases. A derivation of these flow
equations from the consideration of the limiting cases is presented in App. B.

##  V Fixed point analysis

We analyze the fixed points of the above RG flow in several steps of
generalization, starting with a purely bosonic system. To that end, we first
explore the constrained version of our full model, with $r_{a}=r_{b}$,
$\lambda_{a}=\lambda_{b}$, cf. Sec. V.1. Then we include a finite number of
fermions setting also $g_{a}=g_{b}$, and display the mechanism of fixed-point
appearance and collision, see Sec. V.2. Finally, we discuss the stable fixed-
point solutions of the full SO($n_{a}$)$\times$SO($n_{b}$) system without
constraints on the model parameters in Sec. V.3.

###  V.1 Purely bosonic case

We begin by recalling the situation for the bosonic SO(3)$\times$SO(3) sub-
sector without fermions, i.e. $n_{a}=n_{b}=3$ and $N_{f}=0$. As well known,
this limit features four different fixed points in the one-loop analysis, i.e.
the Gaussian fixed point (GFP), the decoupled fixed point (DFP), the SO(6)
isotropic fixed point (IFP), and the biconal fixed point (BFP)calabrese ;
kosterlitz1976 .

Figure 1:  RG flow at the critical surface (left panel) and the geometrical
solution for the fixed points (right panel) when $N_{f}=0$ and
$n_{a}=n_{b}=3$. The blue and the red lines are the points where the two beta-
functions vanish, and the dashed gray line denotes $\lambda^{\prime}=\lambda$,
where the theory becomes O(6)-symmetric.

Let us be on the critical surface $r_{a}=r_{b}=0$, and for simplicity of
presentation restrict the remaining quartic couplings
$\lambda=\lambda_{a}=\lambda_{b}$, and $\lambda^{\prime}=\lambda_{ab}$. The
fixed-point coordinates then read $\lambda=0,\lambda^{\prime}=0$ (GFP),
$\lambda=\frac{\epsilon}{44},\lambda^{\prime}=0$ (DFP),
$\lambda=\frac{\epsilon}{56},\lambda^{\prime}=\frac{\epsilon}{56}$ (IFP), and
$\lambda=\frac{3\epsilon}{136},\lambda^{\prime}=\frac{\epsilon}{136}$ (BFP),
see left panel of Fig. 1. To determine the stability of the fixed points, we
calculate the stability matrix and evaluate it at the fixed points. The
stability matrix of a fixed point $(\lambda_{\ast},\lambda^{\prime}_{\ast})$
is defined as

$\displaystyle\mathcal{S}=\begin{pmatrix}\frac{\partial\beta_{\lambda}}{\partial\lambda}&\frac{\partial\beta_{\lambda^{\prime}}}{\partial\lambda}\\\
\frac{\partial\beta_{\lambda}}{\partial\lambda^{\prime}}&\frac{\partial\beta_{\lambda^{\prime}}}{\partial\lambda^{\prime}}\end{pmatrix}\Big{|}_{(\lambda,\lambda^{\prime})\to(\lambda_{\ast},\lambda^{\prime}_{\ast})}\,,$
(22)

and its eigenvalues determine the stability of the respective fixed point,
i.e. a positive eigenvalue corresponds to an infrared repulsive (unstable,
relevant) direction, and a negative eigenvalue corresponds to an infrared
attractive (stable, irrelevant) direction of the RG flow. (The parameters
$r_{a}$ and $r_{b}$ that have been set to zero correspond of course to
relevant directions.) Therefore, a stable fixed point needs both eigenvalues
$(\theta_{1},\theta_{2})$ of its stability matrix $\mathcal{S}$ to be
negative.

We calculate the stability of the bosonic fixed points as listed above and
find that the BFP is the only stable fixed point at one-loop, with the
eigenvalues $\theta_{1}=-\epsilon,\theta_{2}=-\epsilon/17$. It is well known,
however, that the one-loop calculation significantly overestimates the regime
of stability of the BFP, as well as of the IFP for smaller number of bosonic
field components. In a higher-order calculation it eventually turns out that
in a O($n_{a}$)$\oplus$O($n_{b}$) theory of two coupled bosonic order
parameters, the DFP is stable for $n_{a}+n_{b}\gtrsim 4$, cf. Ref. calabrese,
.

We now take a closer look at the fixed points in the purely bosonic theory,
which will be helpful once we include the fermions. The fixed points are
defined by the equations

$\displaystyle\beta_{\lambda}=0\,,\quad\text{and}\quad\beta_{\lambda^{\prime}}=0\,.$
(23)

The first equation defines an ellipse, shown in red in the right panel of Fig.
1. The second equation can be readily factorized to read

$\displaystyle\beta_{\lambda^{\prime}}=\lambda^{\prime}(\epsilon-16\lambda^{\prime}-40\lambda)\,,$
(24)

which defines two straight lines, i.e.

$\displaystyle\lambda^{\prime}=0,\quad\lambda^{\prime}=(\epsilon-40\lambda)/16\,,$
(25)

as shown in blue in Fig. 1. The fixed points to this order of calculation are
therefore given by the four intersections of the two straight lines with the
ellipse. Furthermore, one intersection of the second straight line and the
ellipse also lies on the line $\lambda=\lambda^{\prime}$ (dashed gray line in
Fig. 1).

Higher-loop corrections to the flow would deform the ellipse and the second
straight line, but still keep one of their intersections at the
$\lambda=\lambda^{\prime}$ SO(6)-symmetric line. $\lambda^{\prime}=0$ also
remains a solution of $\beta_{\lambda^{\prime}}=0$ at all orders of
calculation. So the IFP and DFP always exist somewhere on these two straight
RG-invariant lines. The only fixed point not constrained by the symmetries is
the BFP, which can pass through IFP and DFP with a change of parameters, but
cannot annihilate them.

###  V.2 Inclusion of fermions for $n_{a}=n_{b}$

The above geometrical picture behind the fixed points of the RG flow changes
when fermions are included, as we discuss next, still focusing on the the case
of SO(3)$\times$SO(3) symmetry, i.e. $n_{a}=n_{b}=3$. We also assume
$g_{a}^{2}=g_{b}^{2}=g^{2}$. In this case, the ellipse is slightly deformed,
but the main change is the avoidance of the crossing of the two straight lines
of the solutions of $\beta_{\lambda^{\prime}}(N_{f}=0)=0$ due to the term
$\propto N_{f}g^{4}_{i}$, cf. Eq. (21). The resulting picture is presented in
Fig. 2 for four different choices of $N_{f}$.

Increasing the number of fermions $N_{f}$ starting from $N_{f}=0$, we first
observe that the positions of the intersections resulting from Eq. (23), i.e.
the fixed points, approach each other until they collide at some $N_{f}$. This
holds for both pairs of intersections originating from the purely bosonic
fixed points, i.e. the pair $\\{$IFP, BFP$\\}$ and the pair $\\{$GFP,
DFP$\\}$. For the collision of the first pair, we numerically find
$N_{f,c}\approx 0.0164$, so that for $N_{f,0}>N_{f}>N_{f,c}$ there is no
stable fixed point. The second pair collides at $N_{f,0}\approx 1.485$.

Figure 2: Evolution of the fixed points with the increase of the number of
Dirac fermions. The faint dashed lines represent the case without fermions,
$N_{f}=0$, and are shown for comparison.

To further explore the landscape of fixed-point solutions of the model, we
treat the number of bosonic components $n_{a}=n_{b}=n$ as a free variable and
keep the parameter constraints, i.e.
$r_{a}=r_{b},\lambda_{a}=\lambda_{b},g_{a}=g_{b}$. To be explicit, we show the
RG $\beta$ functions that can be extracted directly from Eqs. (17) to (21) as

$\displaystyle\beta_{g^{2}}$ $\displaystyle=\epsilon
g^{2}-2(N_{f}+n+2)g^{4}\,,$ (26) $\displaystyle\beta_{\lambda}$
$\displaystyle=\epsilon\lambda-4(n+8)\lambda^{2}-4n\lambda^{\prime
2}-4N_{f}\lambda g^{2}+N_{f}g^{4}\,,$ (27)
$\displaystyle\beta_{\lambda^{\prime}}$
$\displaystyle=\epsilon\lambda^{\prime}-16\lambda^{\prime
2}\\!-\\!8(2\\!+\\!n)\lambda\lambda^{\prime}\\!-\\!4N_{f}\lambda^{\prime}g^{2}\\!+\\!3N_{f}g^{4}\,.$
(28)

The resulting function $N_{f,c}(n)$ is shown in Fig. 3. There is a small but
finite number of fermions that causes the first fixed-point collision after
which the stable fixed point’s coordinates acquire an imaginary part.
Therefore, no stable real fixed-point solution exists for the theory at
$N_{f}=2$, at least at leading order in the epsilon expansion. Consequently,
the system is not expected to feature true quantum critical behavior.

Our results are in conflict with the numerical findings and the
renormalization group analysis in Ref. liu, . The difference in the RG results
can be traced back to the fermionic contribution in Eq. (21), i.e. the term
$3N_{f}g_{a}^{2}g_{b}^{2}$ which is reported with a different prefactor in
Ref. liu, . We have carefully checked our present calculation, and also
confirmed that our results are in agreement with earlier works, cf. Refs.
classen,  and janssen, . An independent argument corroborating our finding is
also presented in App. B.

Figure 3: $N_{f,c}$ as a function of the number of boson components $n$. A
finite number of fermions tentatively causes a fixed-point collision at which
the stable fixed-point solution is annihilated.

We also note that for the special case of $n=1$ there is an admissible fixed
point for all $N_{f}$. Here, the Yukawa coupling has the fixed-point value
$g^{2}_{\ast}=\epsilon/(6+2N_{f})$ and the bosonic couplings have the values

$\displaystyle\lambda_{\ast}$
$\displaystyle=\frac{3-N_{f}+\sqrt{9+N_{f}(66+N_{f})}}{144(3+N_{f})}\epsilon\,,$
(29) $\displaystyle\lambda^{\prime}_{\ast}$
$\displaystyle=\frac{3-N_{f}+\sqrt{9+N_{f}(66+N_{f})}}{48(3+N_{f})}\epsilon\,.$
(30)

This fixed point is stable. For example, for $N_{f}=1$ we find
$\theta_{1}\approx-0.39\epsilon,\theta_{2}=-\epsilon,\theta_{3}\approx-2.18\epsilon$.
For $N_{f}=2$ we find
$\theta_{1}\approx-0.67\epsilon,\theta_{2}=-\epsilon,\theta_{3}\approx-2.41\epsilon$.

###  V.3 General SO($n_{a}$)$\times$SO($n_{b}$) symmetry

The full system with two coupled and competing order parameters has a larger
coupling constant space than the one we studied so far. In fact, it allows the
two Yukawa couplings and also the two quartic couplings to be different. In
that case the relevant RG equations are Eqs. (17) – (21), which feature a
stable fixed-point solution, in contrast to the constrained system.

To analyze the fixed points we start with the flow of the two Yukawa
couplings, since this subset of flow equations is independent of the quartic
couplings. The subset of Yukawa beta-functions has four different fixed-point
solutions $Y_{i},\ i\in\\{1,2,3,4\\}$, reading

$\displaystyle Y_{1}:\quad$ $\displaystyle g_{a}^{2}=0,\ g_{b}^{2}=0\,,$ (31)
$\displaystyle Y_{2}:\quad$ $\displaystyle g_{a}^{2}=0,\
g_{b}^{2}=\frac{1}{4+2N_{f}-n_{b}}\epsilon\,,$ (32) $\displaystyle
Y_{3}:\quad$ $\displaystyle g_{a}^{2}=\frac{1}{4+2N_{f}-n_{a}}\epsilon,\
g_{b}^{2}=0\,,$ (33)

and $Y_{4}$ with the Yukawa couplings

$\displaystyle g_{a}^{2}$
$\displaystyle=\frac{2n_{b}-N_{f}-2}{n_{a}(2\\!+\\!4n_{b}\\!+\\!N_{f})\\!+\\!(2\\!+\\!N_{f})(n_{b}\\!-\\!4\\!-\\!2N_{f})}\epsilon,$
(34) $\displaystyle g_{b}^{2}$
$\displaystyle=\frac{2n_{a}-N_{f}-2}{n_{a}(2\\!+\\!4n_{b}\\!+\\!N_{f})\\!+\\!(2\\!+\\!N_{f})(n_{b}\\!-\\!4\\!-\\!2N_{f})}\epsilon\,.$
(35)

Fixed point $Y_{4}$ is the one that has been studied before in the
SO($n$)$\times$SO($n$) system, i.e. when $n_{a}=n_{b}$ and we enforced
$g_{a}^{2}=g_{b}^{2}$. We can determine the stability matrix in the Yukawa
subsector and diagonalize it. There we find that for $Y_{4}$ the two
eigenvalues are

$\displaystyle\theta_{g,1}$ $\displaystyle=-\epsilon,$ (36)
$\displaystyle\theta_{g,2}$
$\displaystyle=\frac{2(2-2n_{a}+N_{f})(2-2n_{b}+N_{f})}{n_{a}(2\\!+\\!4n_{b}\\!+\\!N_{f})\\!+\\!(2\\!+\\!N_{f})(n_{b}\\!-\\!4\\!-\\!2N_{f})}\epsilon$
(37)

For $N_{f}=2$ and $n_{a}=n_{b}$, we find
$\theta_{g,2}=2\frac{n_{a}-2}{n_{a}+4}$, so for $n_{a}>2$ this fixed point is
necessarily unstable, as $\theta_{g,2}>0$.

We next leave $N_{f}$ and $n_{a}$ as free parameters, fix $n_{b}=3$, and
consider $Y_{2}$. The eigenvalues of the stability matrix are

$\displaystyle Y_{2}:\quad\theta_{g,1}=-\epsilon,\
\theta_{g,2}=2\frac{N_{f}-4}{1+2N_{f}}\epsilon\,.$ (38)

For $N_{f}<4$, we find $\theta_{g,2}<0$ and therefore this fixed point is
potentially stable, provided that the bosonic sector also admits a stable
solution. It turns out that this actually is the case and the fixed point
coordinates of this stable fixed point in the bosonic sector read

$\displaystyle\lambda_{a}$
$\displaystyle=\frac{\epsilon}{4(n_{a}+8)}\,,\quad\lambda_{ab}=0\,,$ (39)
$\displaystyle\lambda_{b}$
$\displaystyle=\frac{1-2N_{f}+\sqrt{1+4N_{f}(43+N_{f})}}{88(1+2N_{f})}\epsilon\,.$
(40)

We note that this is in complete agreement with the findings in Ref. classen,
upon setting $n_{a}=1$. Diagonalization of the complete stability matrix for
all five couplings provides three additional eigenvalues which are

$\displaystyle\theta_{\lambda,1}$
$\displaystyle=-\epsilon\,,\quad\theta_{\lambda,2}=-\frac{1+4N_{f}(43+N_{f})}{1+2N_{f}}\epsilon\,,$
(41) $\displaystyle\theta_{\lambda,3}$
$\displaystyle=\frac{6\epsilon}{n_{a}\\!+\\!8}-\frac{5\\!+\\!34N_{f}+5\sqrt{1\\!+\\!4N_{f}(43\\!+\\!N_{f})}}{22+44N_{f}}\epsilon\,.$
(42)

The first two are clearly negative. We evaluate the last one for $N_{f}=2$ to
find

$\displaystyle\theta_{\lambda,3}(n_{a})=(-\frac{84}{55}+\frac{6}{n_{a}+8})\epsilon\,,$
(43)

which is negative for $n_{a}=1,2,3$, i.e.
$\theta_{\lambda,3}(1)=-142/165,\theta_{\lambda,3}(2)=-51/55,\theta_{\lambda,3}(3)=-54/55$.
Therefore, this represents the stable fixed point of the RG flow, including
the case when $n_{a}=n_{b}=3$, i.e. the SO(3)$\times$SO(3) symmetry we are
most interested in.

Figure 4: Flow in the plane of the two Yukawa couplings in the unconstrained
SO($n$)$\times$SO($n$) model, for $n=1,2,3$. Left: $n=1$. Middle $n=2$. Right
$n=3$. For $n=2$ a line of fixed points emerges, indicated by the dashed line.

For the further discussion, we proceed with the concrete case of
$n_{a}=n_{b}=3$ and recall that – at least at one-loop order – the stable
fixed point is the one with $Y_{2}$ in the Yukawa subsector and which is
decoupled in the bosonic sector. That means that at the fixed point, the
$\vec{a}$ field completely decouples from both, the fermions and the $\vec{b}$
field (this is true for any $n_{a}=1,2,3$). Therefore, this sector of the
model belongs to the purely bosonic Heisenberg universality class whereas the
sector with the fermions and the $\vec{b}$ field belongs to the chiral
Heisenberg universality class. In fact, the critical exponents for these two
universality classes are known beyond the leading-order $\epsilon$
expansionhasenbusch2011 ; kos2016 . The correlation length exponent for the
three-dimensional Heisenberg universality class is very precisely known to be
$\nu_{a}=0.712$, cf. Refs. book, . The value of the correlation length
exponent for the chiral Heisenberg universality class is less certain, with
estimates laying in the range $0.84<\nu_{b}<1.31$, cf. Ref. zerf, .

Employing Aharony’s scaling relation book ; calabrese  we can deduce the RG
eigenvalue for the mixed quartic coupling that would couple the $\vec{a}$ and
$\vec{b}$ fields at the decoupled fixed point. The scaling dimension of the
mixed quartic coupling reads

$\displaystyle\theta_{ab}=\frac{1}{\nu_{a}}+\frac{1}{\nu_{b}}-D\,,$ (44)

where $D$ is the space-time dimension. For the decoupled fixed point to be
stable $\theta_{ab}$ needs to be negative. In leading-order $\epsilon$
expansion, we find that
$\frac{1}{\nu_{a}}=2-\frac{n_{a}+2}{n_{a}+8}\epsilon=2-\frac{5}{11}\epsilon$
and $\frac{1}{\nu_{b}}=2-\frac{84}{55}\epsilon$. Insertion into the expression
above yields

$\displaystyle\theta_{ab}=2-\frac{5}{11}\epsilon+2-\frac{84}{55}\epsilon-(4-\epsilon)=-\frac{54}{55}\epsilon,$
(45)

at the leading order in $\epsilon$-expansion. As a check, we observe that this
value indeed coincides with $\theta_{\lambda,3}(n_{a}=3)$ from Eq. (43).

Using the best current estimates for $1/\nu_{a,b}$ from other methods, we find
in $D=3$

$\displaystyle\theta_{ab}\approx\frac{1}{0.712}+\frac{1}{\nu_{b}}-3\in[-0.83,-0.41]\,.$
(46)

The scaling dimension $\theta_{ab}$ seems therefore likely to be negative and
the coupling $\lambda_{ab}$ irrelevant even at higher orders, supporting the
scenario where the Heisenberg+chiral Heisenberg fixed point would remain
stable beyond one loop. What is still missing to corroborate the stability of
the Heisenberg+chiral Heisenberg fixed point is an argument for the RG scaling
of the (vanishing) Yukawa coupling $g_{a}^{2}$.

We note that for $n_{a}=n_{b}$, the situation is somewhat special: In that
case the fixed points $Y_{2}$ and $Y_{3}$ are identical upon exchange of the
labels $a\leftrightarrow b$. Therefore, $Y_{3}$ also has to be a stable fixed
point if $Y_{2}$ is stable. We explore this stability property in the Yukawa
couplings for various configurations of $n_{a}=n_{b}=n\in\\{1,2,3\\}$, see
Fig. 4. Interestingly, for $n=1$ the fixed point $Y_{4}$ is stable and for
$n=3$ the fixed points $Y_{2},Y_{3}$ are stable. For $n=2$, the fixed points
$Y_{2},Y_{3},Y_{4}$ all have a marginal direction in the RG flow.

For the case of $n_{a}=n_{b}=1,N_{f}=2$ it is therefore interesting to also
look at the possible fixed points and their stability in the bosonic sector.
There we find a stable fixed point with coordinates

$\displaystyle\lambda_{a}=\frac{1+\sqrt{145}}{720}\epsilon=\lambda_{b},\quad\lambda_{ab}=\frac{1+\sqrt{145}}{240}\epsilon\,.$
(47)

We note that this is also not a fixed point with symmetry enlargement. Also
this is in agreement with the findings discussed around Eq. (29).

##  VI Discussion

We have employed a one-loop renormalization-group approach to study the
quantum multicritical behavior in a generalized Gross-Neveu-Yukawa field
theory for Dirac fermions coupled to two order-parameter fields. This work was
partly motivated by the quantum Monte Carlo (QMC) study of the lattice
Hamiltonian presented in Ref. liu, , where the quantum phase transition
between a $N_{f}=2$ Dirac fermion phase and a massive phase including a Néel
state and a superconductor-CDW state was investigated. Therein, a scaling
collapse of the QMC data is found at critical coupling strength, providing
evidence for a continuous quantum phase transition with a concomitant quantum
critical point.

Universality suggests that this putative quantum critical point is described
by the constrained continuum field theory studied in our Sec. V.2. As we have
shown, however, there is no stable fixed point for $n=3$ and $N_{f}=2$, i.e.
the numerical findings appear to be incompatible with the universal continuum
field theory. We therefore conclude that either (1) the lattice Hamiltonian of
Ref. liu,  should be described by a different continuum field theory, (2) the
leading-order perturbative RG approach presented here is insufficient to
capture the underlying fixed-point structure, or (3) the QMC data actually
exhibits a (very) weak first-order transition, which only appears to be
continuous.

About point (1), we can say that the continuum theory is constructed such that
it shares the full O(4) symmetry of the lattice Hamiltonian, cf. Ref. liu, . A
different continuum theory should also fulfill that minimal requirement, which
imposes a strong constraint. One possibility could be that topological
contributions appear in the lattice realization, while they are not reflected
in the present continuum setup. Option (2) would also be interesting since
examples where the one-loop order fails to capture the general fixed-point
structure of a model are rare. Possibilities for further studies include
higher-order epsilon expansions, but also a large-$N_{f}$ study at leading
order. The latter seems promising since at large $N_{f}\sim O(10)$ fermionic
contributions in beta-function dominate and the critical point reappears
LHKC2022 . An unusually large order parameter’s anomalous dimension detected
in Ref. liu,  may be an indication of the appropriateness of the large-$N_{f}$
approach, in which $\eta_{\psi}=1+O(1/N_{f})$. Finally, it may be
computationally expensive to decide on the option (3), as weakly first-order
and continuous transitions are naturally difficult to distinguish. Weakly
first-order transitions are expected to appear in a number of spin systems. In
the present context it would seem to require that the critical number of
fermions for the disappearance of the stable fixed point $N_{f,c}$ is below
but closer to the physical number $N_{f}=2$ than what we found. Higher order
calculations may be able to shed more light on the matter.

###### Acknowledgements.

We thank Shailesh Chandrasekharan, Emilie Huffman, Ribhu Kaul, and Hanqing Liu
for helpful discussions. MMS acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) through the DFG
Heisenberg programme (project id 452976698) and SFB 1238 (project C02, project
id 277146847). IFH is supported by the NSERC of Canada.

##  Appendix A Unitary transformation between different SO(3)$\times$SO(3)
vector order parameters

The Hamiltonian $H$ in Eq. (3) is particle-hole symmetric, i.e. it
anticommutes with the antiunitary particle-hole transformation operator $P$.
Since $P$ has the property $P^{2}=+1$, there exists a representation where
$P=K$, with $K$ being complex conjugation, and $P$ with a trivial unitary
part. In this special basis, the $16\\!\times\\!16$ BdG Hamiltonian reads

$\displaystyle H=p_{1}R_{1}+p_{2}R_{2}\,,$ (48)

where $\vec{p}=(p_{1},p_{2})=-i\nabla$ is the usual momentum operator, odd
under $P$, and $R_{1},R_{2}$ are two real Hermitian (symmetric)
$16\\!\times\\!16$ matrices with the anticommuting property

$\displaystyle\\{R_{i},R_{j}\\}=2\delta_{ij}\,,\quad i,j\in\\{1,2\\}\,.$ (49)

The largest Clifford algebra that has a $16\\!\times\\!16$ representation is

$\displaystyle
C(5,4)=(R_{1},R_{2},R_{3},R_{4},R_{5},I_{1},I_{2},I_{3},I_{4})\,,$ (50)

where the $R_{i},i\in\\{1,2,3,4,5\\}$ are real and $I_{j},j\in\\{1,2,3,4\\}$
are fully imaginary (antisymmetric). All matrices $R_{i},I_{j}$ are hereafter
chosen to be Hermitian, and to square to $+1$. They all also mutually
anticommute.

The Clifford algebra $C(n,m)$ is here defined as a set of $n+m$ mutually
anticommuting Hermitian generators which all square to $+1$, and are such that
the first $n$ are real, and the remaining $m$ are imaginary. Note that this is
a slight modification of the standard definitions as, e.g., given in Ref.
okubo , where the fully real representation of $C(n,m)$ are studied, at the
cost of the last $m$ elements being multiplied with an imaginary unit, and
therefore taken as anti-Hermitian. In consequence the last $m$ matrices in our
definition square to $+1$ instead of $-1$, while preserving the nature of the
representation. The difference between the two conventions reflects the
difference between Minkowski and Euclidian spaces typically employed in
particle and condensed matter physics.

In the present context of the 16$\times$16 massive Dirac BdG Hamiltonian in
two spatial dimensions, we have $n=2$, and then $m=5$ sets the maximal number
of compatible (anticommuting) masses.herbut12  This is because one can
construct the Clifford algebra

$\displaystyle C(2,5)=(R_{1},R_{2},I_{1},I_{2},I_{3},I_{4},I_{5})\,,$ (51)

with $I_{5}=iR_{5}R_{4}R_{3}$, with $I_{5}^{\dagger}=I_{5}$ and
$I_{5}^{2}=+1$. This algebra can be used to write the massive BdG-Dirac
Hamiltonian with five mutually anticommuting mass terms, i.e.

$\displaystyle
H_{\mathrm{BdG}}=p_{1}R_{1}+p_{2}R_{2}+\sum_{i=1}^{5}m_{i}I_{i}\,,$ (52)

where the $m_{i}\in\mathbb{R}$, with the spectrum

$\displaystyle\epsilon_{\vec{p},\pm}=\pm\sqrt{p_{1}^{2}+p_{2}^{2}+m_{1}^{2}+m_{2}^{2}+m_{3}^{2}+m_{4}^{2}+m_{5}^{2}}\,.$
(53)

Note that this Hamiltonian anticommutes with the particle-hole operator $P$,
as it should, because the parameters $\\{m_{i}\\}$ are real and therefore even
under $P$, whereas the momentum $\vec{p}$ is odd. The full summary of the
dimensions and the nature of real irreducible representations of Clifford
algebras provided in Ref. herbut12,  implies then that there cannot be more
than five mutually anticommuting mass-terms that can be added to
$H_{\mathrm{BdG}}$. The $C(2,5)$ is thus and the maximal Clifford algebra
relevant to spin-1/2 Dirac fermions in two dimensions.

An interesting fact about this algebra is that it is quaternionic okubo . In
our convention this means that, besides unit matrix, it allows for three
additional Casimir operators

$\displaystyle(A_{1},A_{2},A_{3})=(iR_{3}R_{4},iR_{4}R_{5},iR_{5}R_{3})\,,$
(54)

which evidently commute with all seven matrices in $C(2,5)$. The three $A_{i}$
close an SO(3) Lie algebra, mutually anticommute, and are imaginary and
Hermitian.

This observation can be used to derive the central result quoted in the main
text as follows. Take a subalgebra of the above $C(2,5)$ to be

$\displaystyle C(2,2)=(R_{1},R_{2},I_{1},I_{2})\,.$ (55)

Then there exists precisely one other triplet besides $(I_{3},I_{4},I_{5})$
which may be used to complete this subalgebra to the maximal, but different,
$C(2,5)$. This new triplet will transform as a vector under the above SO(3)
algebra of Casimir operators, and it reads

$\displaystyle(A_{1}X,A_{2}X,A_{3}X)\,\ \text{with}\
X=R_{1}R_{2}I_{1}I_{2}\,.$

These three matrices are also all imaginary, Hermitian, and square to $+1$. We
thus found another maximal $C(2,5)$ algebra which contains the same subalgebra
$\\{R_{1},R_{2},I_{1},I_{2}\\}$ as the first one. Explicitly, it reads

$\displaystyle C_{A}(2,5)=(R_{1},R_{2},I_{1},I_{2},A_{1}X,A_{2}X,A_{3}X)\,.$

To distinguish it from the first one we added an index $A$ in $C_{A}(2,5)$,
and we refer to the first one as $C_{B}(2,5)$ hereafter. $C_{A}(2,5)$ and
$C_{B}(2,5)$ share the same subalgebra $C(2,2)=(R_{1},R_{2},I_{1},I_{2})$.
$C_{A}(2,5)$ is of course also quaternionic, and the three non-trivial
Casimirs are simply

$\displaystyle(B_{1},B_{2},B_{3})=(iI_{3}I_{4},iI_{4}I_{5},iI_{5}I_{3})\,$
(56)

$B_{i}$ are also Hermitian and mutually anticommuting, and provide another
imaginary representation of SO(3) Lie algebra.

In sum, given the 16$\times$16 representation of $C(2,2)$ as
$(R_{1},R_{2},I_{1},I_{2})$, there are two different ways to complete it to
the maximal Clifford algebra $C(2,5)$:

  1. 1.

$C_{B}(2,5)=(R_{1},R_{2},I_{1},I_{2},I_{3},I_{4},I_{5})$, with Casimir
operators $A_{i}$, $[A_{i},C_{B}(2,5)]=0,i\in\\{1,2,3\\}$,

  2. 2.

$C_{A}(2,5)\\!=\\!(R_{1},R_{2},I_{1},I_{2},A_{1}X,A_{2}X,A_{3}X)$ with Casimir
operators $B_{i}$, $[B_{i},C_{A}(2,5)]\\!=\\!0,\\!i\\!\in\\{1,2,3\\}$.

The two vectors $(A_{1}X,A_{2}X,A_{3}X)$ and $(I_{3},I_{4},I_{5})$ are a
$(1,0)+(0,1)$ representation of the SO(3)$\times$SO(3) algebra of Casimirs,
which are

$\displaystyle(iR_{3}R_{4},iR_{4}R_{5},$ $\displaystyle
iR_{5}R_{3},iI_{3}I_{4},iI_{4}I_{5},iI_{5}I_{3})$
$\displaystyle\quad=(A_{1},A_{2},A_{3},B_{1},B_{2},B_{3})\,,$ (57)

with $[A_{i},B_{j}]=0$. In other words, the algebra
$C(2,2)=(R_{1},R_{2},I_{1},I_{2})$ has the SO(4)$\simeq$SO(3)$\times$SO(3)
algebra of Casimirs.

One can also summarize the above discussion by noting that

$\displaystyle R_{5}=$ $\displaystyle\
R_{1}R_{2}R_{3}R_{4}I_{1}I_{2}I_{3}I_{4}$ $\displaystyle\Rightarrow
B_{1}X=I_{5},B_{2}X=I_{3},B_{3}X=I_{4}\,.$ (58)

So the Clifford algebra $C(2,2)\\!=\\!(R_{1},R_{2},I_{1},I_{2})$, has the
SO(4) Lie algebra of Casimir operators

$\displaystyle(\vec{A},\vec{B})$
$\displaystyle=(iR_{3}R_{4},iR_{4}R_{5},iR_{5}R_{3},iI_{3}I_{4},iI_{4}I_{5},iI_{5}I_{3})$
$\displaystyle=\mathrm{SO}(3)\times\mathrm{SO}(3)\,.$ (59)

$(\vec{A}X,\vec{B}X)$ with $X=R_{1}R_{2}I_{1}I_{2}$ then transforms as
$(1,0)+(0,1)$ representation of this SO(4) Lie algebra. Finally,

  * $C_{A}(2,5)=C(2,2)+(\vec{A}X)$,

  * $C_{B}(2,5)=C(2,2)+(\vec{B}X)$,

are the two maximal Clifford algebras that contain the original $C(2,2)$. This
was quoted in the text.

### Transformation

Say we are given the two maximal Clifford algebras $C_{A}(2,5)$ and
$C_{B}(2,5)$ with the common subalgebra $C(2,2)=(R_{1},R_{2},I_{1},I_{2})$. We
can construct two Clifford algebras with any other pair of imaginary matrices
from $C_{A}(2,5)$ (or $C_{B}(2,5)$) being common as follows:

Define

$\displaystyle
U=e^{i\frac{\pi}{4}\left[iI_{1}(A_{1}X)+iI_{2}(A_{2}X)\right]}\,.$ (60)

Its action on $C_{A}(2,5)$ is as follows:

  1. 1.

$UR_{i}U^{\dagger}=R_{i},\ i\in\\{1,2\\}$

  2. 2.

$UI_{1}U^{\dagger}=e^{i\frac{\pi}{4}iI_{1}A_{1}X}I_{1}e^{-i\frac{\pi}{4}iI_{1}A_{1}X}\\\
{}\quad\quad\ \ \
=e^{i\frac{\pi}{2}iI_{1}A_{1}X}I_{1}=\\!-\\!\left(\sin\frac{\pi}{2}\right)I_{1}A_{1}XI_{1}=A_{1}X$

  3. 3.

$UI_{2}U^{\dagger}=A_{2}X$

  4. 4.

$UA_{1}XU^{\dagger}=-I_{1}$

  5. 5.

$UA_{2}XU^{\dagger}=-I_{2}$

  6. 6.

$UA_{3}XU^{\dagger}=A_{3}X$

The net result is the same set of matrices but with exchanges of
$I_{1}\leftrightarrow A_{1}X,I_{2}\leftrightarrow A_{2}X$. The action on the
matrices in $C_{B}(2,5)$, on the other hand, is

  1. 1.

$UR_{i}U^{\dagger}=R_{i},\ i\in\\{1,2\\}$

  2. 2.

$UI_{1}U^{\dagger}=A_{1}X$

  3. 3.

$UI_{2}U^{\dagger}=A_{2}X$

  4. 4.

$UB_{i}XU^{\dagger}\\\
{}=\\!e^{i\\!\frac{\pi}{4}(iI_{1}\\!A_{1}\\!X)}e^{i\frac{\pi}{4}(iI_{2}\\!A_{2}\\!X)}\\!B_{i}Xe^{-i\\!\frac{\pi}{4}(iI_{2}\\!A_{2}\\!X)}e^{-i\frac{\pi}{4}(iI_{1}\\!A_{1}\\!X)}\\\
{}=\\!e^{i\frac{\pi}{4}(iI_{1}A_{1}X)}\left[-I_{2}A_{2}XB_{i}X\right]e^{-i\frac{\pi}{4}(iI_{1}A_{1}X)}\\\
{}=-I_{1}(A_{1}X)(-)I_{2}A_{2}XB_{i}X\\\
{}=I_{1}A_{1}XI_{2}A_{2}B_{i}(X^{2})\\\ {}=-I_{1}A_{1}I_{2}A_{2}B_{i}X\\\
{}=-I_{1}I_{2}A_{1}A_{2}B_{i}X\\\ {}=-iI_{1}I_{2}A_{3}B_{i}X$.

So, we obtain

$\displaystyle C_{B}(2,5)\rightarrow\ $ $\displaystyle UC_{B}(2,5)U^{\dagger}$
$\displaystyle=(R_{1},R_{2},A_{1}X,A_{2}X,iI_{1}I_{2}A_{3}B_{i}X)\,,$ (61)

with $i\in\\{1,2,3\\}$.

Let us make couple of remarks.

  1. 1.

The matrices $iI_{1}I_{2}A_{3}B_{i}X$ are imaginary, since $X$ is real and
$I_{1},I_{2},A_{3},B_{i}$ are imaginary.

  2. 2.

$\\{R_{i},iI_{1}I_{2}A_{3}B_{i}X\\}=0$, because $R_{i},\ i\in\\{1,2\\}$
commutes with $iI_{1}I_{2}$, $A_{3}$, and $B_{i}$, but anticommutes with
$X=R_{1}R_{2}I_{1}I_{2}$.

  3. 3.

$\\{A_{i}X,iI_{1}I_{2}A_{3}B_{i}X\\}=0$, since $A_{i}X,\ i\in\\{1,2\\}$
commutes with $iI_{1}I_{2}$ and $B_{i}X$, but anticommutes with $A_{3}$.

  4. 4.

$iI_{1}I_{2}A_{3}B_{i}X$ mutually anticommute, $i\in\\{1,2,3\\}$.

The final result of the transformation is that

$\displaystyle C_{A}(2,5)\rightarrow$ $\displaystyle\
UC_{A}(2,5)U^{\dagger}=C_{A}(2,5)$
$\displaystyle=(R_{1},R_{2},I_{1},I_{2},A_{i}X)\,,$ (62) $\displaystyle
C_{B}(2,5)\rightarrow$ $\displaystyle\ UC_{B}(2,5)U^{\dagger}$
$\displaystyle=(R_{1},R_{2},A_{1}X,A_{2}X,iI_{1}I_{2}A_{3}B_{i}X)\,,$ (63)

so that the transformed Cliffod algebras have the imaginary (mass) matrices
$(A_{1}X,A_{2}X)$ in common, besides the real matrices $(R_{1},R_{2})$, of
course. Again, the matrices from the two different remaining triplets
$(I_{1},I_{2},A_{3})$ and $(iI_{1}I_{2}A_{3}B_{i}X)$ commute between
themselves.

In the original notation, we obtain the relation: if
$C_{B}(2,5)=(R_{1},R_{2},I_{1},I_{2},I_{3},I_{4},I_{5})$ and
$UC_{B}(2,5)U^{\dagger}=(R_{1},R_{2},I_{3},I_{4},I_{1},I_{2},I_{5})$, then for
$UC_{B}(2,5)U^{\dagger}$ the Casimir $SO(3)$ algebra is
$\vec{A}=(iR_{3}R_{4},iR_{4}R_{5},iR_{5}R_{3})$. Since $X\rightarrow
UXU^{\dagger}=R_{1}R_{2}I_{3}I_{4}$, so that the matrices
$\vec{A}X=(iR_{k}I_{1}I_{2}),\ k\in\\{5,3,4\\}$. So, we find that for
$UC_{A}(2,5)U^{\dagger}=\\{R_{1},R_{2},I_{3},I_{4},iR_{k}I_{1}I_{2}\\},\
k\in\\{5,3,4\\}$ and the new $(1,0)+(0,1)$ representation is

$\displaystyle(I_{1},I_{2},I_{5})+(iR_{k}I_{1}I_{2};k=3,4,5)\,.$ (64)

Again, all matrices in the first set commute with matrices in the second set,
as it should be.

##  Appendix B Derivation of RG flow equations from known limiting cases

The RG equations for a Dirac system with fermion field $\psi$ coupled to two
order parameters with real field components $\phi_{i},\
i\in\\{1,\ldots,n_{a}\\}$ and $\chi_{i},\ i\in\\{1,\ldots,n_{b}\\}$ have a
general structure following from the consideration of rescalings of the
couplings and loop diagrams. The general structure for the $\beta$ functions
of $g_{i}^{2},\lambda_{i},\lambda_{ab}$ where $i\in\\{a,b\\}$ reads

$\displaystyle\beta_{g_{i}^{2}}$
$\displaystyle=\left(\epsilon-\eta_{i}-2\eta_{\psi}\right)g_{i}^{2}+c_{g_{i},1}g_{i}^{4}+c_{g_{i},2}g_{a}^{2}g_{b}^{2}\,,$
(65) $\displaystyle\beta_{\lambda_{i}}$
$\displaystyle=\left(\epsilon-2\eta_{i}\right)\lambda_{i}\\!+\\!c_{\lambda_{i},1}\lambda_{i}^{2}\\!+\\!c_{\lambda_{i},2}\lambda_{ab}^{2}\\!+\\!c_{\lambda_{i},3}N_{f}g_{i}^{4}\,,$
(66) $\displaystyle\beta_{\lambda_{ab}}$
$\displaystyle=\left(\epsilon-\eta_{a}-\eta_{b}\right)\lambda_{ab}+c_{\lambda_{ab},1}\lambda_{ab}^{2}$
$\displaystyle\quad\\!+\\!c_{\lambda_{ab},2}\lambda_{a}\lambda_{ab}\\!+\\!c_{\lambda_{ab},3}\lambda_{b}\lambda_{ab}\\!+\\!c_{\lambda_{ab},4}N_{f}g_{a}^{2}g_{b}^{2}\,,$
(67)

where, from the general point of view, the coefficients $c_{i,j}$ and the
rescaling coefficients $\eta_{\psi},\eta_{i},\ i\in\\{a,b\\}$ are undetermined
at first.

For a model where the bosonic sector features a O($n_{a}$)$\oplus$O($n_{b}$)
symmetry, the coefficients $c_{\lambda_{i},j},i\in\\{a,b\\},j\in\\{1,2\\}$ and
$c_{\lambda_{ab},i},i\in\\{1,2,3\\}$ are well known from statistical models,
see, e.g., Ref. calabrese, . They read

$\displaystyle c_{\lambda_{a},1}$ $\displaystyle=-4(n_{a}+8),\quad
c_{\lambda_{b},1}=-4(n_{b}+8)\,,$ (68) $\displaystyle c_{\lambda_{a},2}$
$\displaystyle=-4n_{b},\quad c_{\lambda_{b},2}=-4n_{a}\,,$ (69) $\displaystyle
c_{\lambda_{ab},1}$ $\displaystyle=-16,\quad
c_{\lambda_{ab},2}=-4(n_{a}+2)\,,$ (70) $\displaystyle c_{\lambda_{ab},3}$
$\displaystyle=-4(n_{b}+2)\,.$ (71)

We now restrict ourselves to the case where the Dirac mass terms related to
the bosonic fields $\phi_{i}$ have anticommuting components and the same for
$\chi_{j}$, but $\phi_{i}$ and $\chi_{j}$ do not have to be mutually
anticommuting. The rescaling coefficients $\eta_{i},\ i\in\\{a,b\\}$ are
related to the wave function renormalizations of the bosons and at a fixed
point they correspond to the boson anomalous dimensions. Diagrammatically,
they are given by the one-loop diagram of the boson propagator with the closed
fermion loop providing a factor of $N_{f}$. The expression for $\eta_{i}$ is
well known, cf. Ref. janssen, , again

$\displaystyle\eta_{i}=2N_{f}g_{i}^{2}\,,\quad i\in\\{a,b\\}$ (72)

Also $\eta_{\psi}$ is known from the one-loop diagram of the fermion
propagator janssen . Here, all boson components contribute additively

$\displaystyle\eta_{\psi}=\frac{1}{2}n_{a}g_{a}^{2}+\frac{1}{2}n_{b}g_{b}^{2}\,.$
(73)

Next, we consider the Yukawa couplings and the corresponding coefficients
$c_{g_{i},1}$ and $c_{g_{i},2}$. The first coefficient, $c_{g,1}$, is the one
obtained from consideration of a Yukawa system with a coupling to a single
order parameter, e.g., a Néel state or a CDW. It can also be inferred from
Ref. janssen,  and it reads

$\displaystyle c_{g_{i},1}=2(n_{i}-2)\,\quad i\in\\{a,b\\}\,.$ (74)

Furthermore, we can use the RG flow equations from the single order parameter
case to extract $c_{\lambda_{i},3}$,

$\displaystyle c_{\lambda_{i},3}=1\,\quad i\in\\{a,b\\}\,.$ (75)

Now, the second Yukawa coefficients $c_{g_{i},2}$ are not that obvious because
they depend on the mutual anticommuting properties of the two different order
parameters $\phi$ and $\chi$. Here, to connect to the work in Ref. liu, , we
can exploit our own RG flow equations from Ref. classen, , where we considered
the competition between the Néel and the CDW state which corresponds to the
case $n_{a}=3,n_{b}=1$. The diagram corresponding to the coefficients
$c_{g_{i},2}$ is the one where the internal bosonic line is contributed from
the order parameter that is not the external bosonic line of the considered
Yukawa vertex. Therefore, the diagram for vertex $g_{a}$ comes with
multiplicity $n_{b}$ and vice versa. The remaining factor is fixed by
consideration of the contribution to the Néel and CDW Yukawa coupling in Ref.
classen, . We find

$\displaystyle c_{g_{a},2}=-2n_{b},\quad c_{g_{b},2}=-2n_{a}\,.$ (76)

The only coefficient left to fix is now $c_{\lambda_{ab},4}$. It follows from
the corresponding equation (Eq. (16)) in Ref. classen,  because of the same
anticommuting properties of the Yukawa vertices for the considered bosonic
components. We find

$\displaystyle c_{\lambda_{ab},4}=3\,.$ (77)

With this, all the coefficients are determined and the set of flow equations
is complete. Putting all that together and slightly rearranging for better
comparison, we obtain the following set of equations given as Eqs. (17) to
(21) in the main text.

Finally, we note that the beta functions for the quartic couplings differ in
the fermionic loop contribution $\propto g^{4}$, which is particularly obvious
when considering the constrained system where $r_{a}=r_{b}$,
$\lambda=\lambda_{a}=\lambda_{b}$, $\lambda^{\prime}=\lambda_{ab}$, and
$g_{a}^{2}=g_{b}^{2}=g^{2}$, cf. Eqs. (27) – (28). As a consequence, in the
presence of fermions, the theory does not have an SO(6) symmetry even at one-
loop level. To see this more directly, we can diagonalize the matrix

$\displaystyle\mathcal{Y}=\vec{a}\cdot\vec{\sigma}\otimes\mathbb{I}+\mathbb{I}\otimes\vec{b}\cdot\vec{\sigma}\,,$
(78)

which appears in the Yukawa coupling and which determines the loop
contributions to the quartic couplings. The eigenvalues of $\mathcal{Y}$ are
$\pm(a\pm b)$. Expanding the fermion determinant for each component separately
leads to a quartic term in the form

$\displaystyle\sim(a+b)^{4}+(a-b)^{4}=\ $
$\displaystyle(a^{4}+4a^{3}b+6a^{2}b^{2}+4ab^{3}+b^{4})$
$\displaystyle+\\!(a^{4}\\!-\\!4a^{3}b+6a^{2}b^{2}\\!-\\!4ab^{3}\\!+\\!b^{4})$
$\displaystyle=\ $ $\displaystyle 2(a^{4}+b^{4}+6a^{2}b^{2})\,.$

Here, we see how a relative factor of 6 arises in the beta-functions for the
couplings $\lambda$ and $\lambda^{\prime}$.

## References

  * (1) I. F. Herbut, Interactions and Phase Transitions on Graphene’s Honeycomb Lattice, Phys. Rev. Lett. 97, 146401, (2006).
  * (2) I. F. Herbut, V. Juričić, and B. Roy, Theory of interacting electrons on honeycomb lattice, Phys. Rev. B 79 , 085116 (2009).
  * (3) I. F. Herbut, V. Juričić, and O. Vafek, Relativsitc Mott criticality in graphene, Phys. Rev. B 80, 075432 (2009).
  * (4) L. Fritz, Quantum critical transport at a semimetal-to-insulator transition on the honeycomb lattice, Phys. Rev. B 83, 035125 (2011).
  * (5) D. Mesterházy, J. Berges, and L. von Smekal, Effect of short-range interactions on the quantum critical behavior of spinless fermions on the honeycomb lattice, Phys. Rev. B 86 245431 (2012).
  * (6) L. Janssen and H. Gies, Critical behavior of the (2+1)-dimensional Thirring model, Phys. Rev. D 86, 105007 (2012).
  * (7) B. Roy, V. Juričić, and I. F. Herbut, Quantum superconducting criticality in graphene and topological insulators, Phys. Rev. B 87, 041401 (2013).
  * (8) L. Janssen and I. F. Herbut, Antiferromagnetic critical point on graphene’s honeycomb lattice: A functional renormalization group approach, Phys. Rev. B 89, 205403 (2014).
  * (9) B. Knorr, Ising and Gross-Neveu model in next-to-leading order, Phys. Rev. B 94, 245102 (2016).
  * (10) M. M. Scherer and I. F. Herbut, Gauge-field-assisted Kekulé quantum criticality, Phys. Rev. B 94 , 205136 (2016).
  * (11) L. N. Mihaila, N. Zerf, B. Ihrig, I. F. Herbut, and M. M. Scherer, Gross-Neveu-Yukawa model at three loops and Ising critical behavior of Dirac systems, Phys. Rev. B 96, 165133 (2017).
  * (12) Z.-X. Li, Y.-F. Jiang, S.-K. Jian, H. Yao, Fermion-induced quantum critical points, Nat. Comm. 8, 314 (2017).
  * (13) L. Classen, I. F. Herbut, M. M. Scherer, Fluctuation-induced continuous transition and quantum criticality in Dirac semimetals, Phys. Rev. B 96, 115132 (2017).
  * (14) N. Zerf, L. N. Mihaila, P. Marquard, I. F. Herbut, and M. M. Scherer, Four-loop critical exponents for the Gross-Neveu-Yukawa models, Phys. Rev. B 96, 096010 (2017).
  * (15) B. Knorr, Critical chiral Heisenberg model with the functional renormalization group, Phys. Rev. B 97, 075129 (2018).
  * (16) E. Torres, L. Classen, I. F. Herbut, and M. M. Scherer, Fermion-induced quantum criticality with two length scales in Dirac systems, Phys. Rev. B 97, 125137 (2018).
  * (17) K. Wamer and I. Affleck, Renormalization group analysis of phase transitions in the two-dimensional Majorana-Hubbard model, Phys. Rev. B 98, 245120 (2018).
  * (18) B. Ihrig, L. N. Mihaila, and M. M. Scherer, Critical behavior of Dirac fermions from perturbative renormalization, Phys. Rev. B 98 125109 (2018).
  * (19) J. A. Gracey, Large N critical exponents for the chiral Heisenberg Gross-Neveu universality class, Phys. Rev. D 97, 105009 (2018).
  * (20) S. Yin, S.-K. Jian, and H. Yao, Chiral tricritical point: A new universality class in Dirac systems, Phys. Rev. Lett. 120, 215702 (2018).
  * (21) Y. Liu, Z.- Y. Meng, and S. Yin, Fermion-enhanced first-order phase transition and chiral Gross-Neveu tricritical point, Phys. Rev. B 103 075147 (2021).
  * (22) H. Yerzhakov and J. Maciejko, Random-mass disorder in the critical Gross-Neveu-Yukawa models,, Nucl. Phys. B 962, 115241 (2021).
  * (23) R. Boyack, H. Yerzhakov, and J. Maciejko, Quantum phase transitions in Dirac fermion systems,, Eur. Phys. J. Spec. Top. 230, 979 (2021).
  * (24) S. Sorella and E. Tosatti, Semi-Metal-Insulator transition of the Hubbard model in the honeycomb lattice, Europhys. Lett. 19, 699 (1992).
  * (25) S. Sorella, Y. Otsuka, and S. Yunoki, Absence of a spin liquid phase in the Hubbard Model on the honeycomb lattice, Sci. Rep. 2, 992 (2012).
  * (26) F. F. Assaad and I. F. Herbut, Pinning the Order: The Nature of Quantum Criticality in the Hubbard Model on Honeycomb Lattice, Phys. Rev. X 3, 031010 (2013).
  * (27) F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb andp-flux Hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum Monte Carlo, Phys. Rev. B 91, 165108 (2015).
  * (28) Y. Otsuka, S. Yunoki, and S. Sorella, Universal quantum criticality in the metal-insulator transition of two-dimensional interacting Dirac electrons, Phys. Rev. X 6, 011029 (2016).
  * (29) P. Buividovich, D. Smith, M. Ulybyshev, and L. von Smekal, Hybrid Monte Carlo study of competing order in the extended fermionic Hubbard model on the hexagonal lattice, Phys. Rev. B 98, 235129 (2018).
  * (30) S. Hesselmann, D. D. Scherer, M. M. Scherer, and S. Wessel, Bond-ordered states and f -wave pairing of spinless fermions on the honeycomb lattice, Phys. Rev. B 98, 045142 (2018).
  * (31) H.-K. Tang, J. N. Leaw, J. N. B. Rodrigues, I. F. Herbut, P. Sengupta, F. F. Assaad, and S. Adam, The role of electron-electron interactions in two-dimensional Dirac fermions, Science 361, 570-574 (2018).
  * (32) T. C. Lang and A. M. Läuchli, Quantum Monte Carlo simulation of the chiral Heisenberg Gross-Neveu-Yukawa phase transition with a single Dirac cone, Phys. Rev. Lett. 123, 137602 (2019).
  * (33) B.-H. Li, Z.-X. Li, and H. Yao, Fermion-induced quantum critical point in Dirac semimetals: A sign-problem-free quantum Monte Carlo study, Phys. Rev. B 101, 085105 (2020).
  * (34) J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Semimetal Mott insulator quantum phase transition of the Hubbard model on the honeycomb lattice, Phys. Rev. B 102, 245105 (2020).
  * (35) Y. Otsuka, K. Seki, S. Sorella, and S. Yunoki, Dirac electrons in the square-lattice Hubbard model with a d-wave pairing field: The chiral Heisenberg universality class revisited, Phys. Rev. B 102, 235105 (2020).
  * (36) E. Huffman and S. Chandrasekharan, Fermion-bag inspired Hamiltonian lattice field theory for fermionic quantum criticality, Phys. Rev. D 101, 074501 (2020).
  * (37) J. Ostmeyer, E. Berkowitz, S. Krieg, T. A. Lähde, T. Luu, and C. Urbach, Antiferromagnetic character of the quantum phase transition in the Hubbard model on the honeycomb lattice, Phys. Rev. B 104, 155142 (2021).
  * (38) S. Tarat, Bo Xiao, R. Mondaini, and R. T. Scalettar, Deconvolving the components of the sign problem, Phys. Rev. B 105, 04510 (2022).
  * (39) B. Roy, Multicritical behavior of $Z_{2}\times O(2)$ Gross-Neveu-Yukawa theory in graphene, Phys. Rev. B 84, 113404 (2011).
  * (40) L. Classen, I. F. Herbut, L. Janssen, and M.M. Scherer, Mott multicriticality of Dirac electrons in graphene, Phys. Rev. B 92, 035429 (2015); Competition of density waves and quantum multicritical behavior in Dirac materials from functional renormalization, Phys. Rev. B 93, 125119 (2016).
  * (41) T. Sato, M. Hohenadler, F. F. Assaad, Dirac Fermions with Competing Orders: Non-Landau Transition with Emergent Symmetry, Phys. Rev. Lett. 119, 197203 (2017).
  * (42) L. Janssen, I. F. Herbut, and M. M. Scherer, Compatible orders and fermion-induced emergent symmetry in Dirac systems, Phys. Rev. B 97, 041117(R) (2018).
  * (43) B. Roy, P. Goswami, and V. Juričić, Itinerant quantum multicriticality of two-dimensional Dirac fermions, Phys. Rev. B 97, 205117(R) (2018).
  * (44) E. Torres, L. Weber, L. Janssen, S. Wessel, M. M. Scherer, Emergent symmetries and coexisting orders in Dirac fermion systems, Phys. Rev. Research 2, 022005 (2020).
  * (45) H. Liu, E. Huffman, S. Chandrasekharan, R. K. Kaul, Quantum Criticality of Antiferromagnetism and Superconductivity with Relativity, Phys. Rev. Lett 128, 117202 (2022).
  * (46) S. Ryu, C. Mudry, C-Y. Hou, and C. Chamon, Masses in graphene-like two-dimensional electronic systems: topological defects in order parameters and their fractional exchange statistics, Phys. Rev. B 80, 205319 (2009).
  * (47) I. F. Herbut, Isospin of topological defects in Dirac systems, Phys. Rev. B 85, 085304 (2012).
  * (48) P. Ghaemi, S. Ryu, D.-H. Lee, Quantum valley Hall effect in proximity-induced superconducting graphene: An experimental window for deconfined quantum criticality, Phys. Rev. B 81, 081403(R) (2010).
  * (49) I. F. Herbut, Topological Insulator in the Core of the Superconducting Vortex in Graphene, Phys. Rev. Lett. 104, 066404 (2010).
  * (50) P. Calabrese, A. Pelissetto, E. Vicari, Multicritical phenomena in $O(n_{1})\oplus O(n_{2})$-symmetric theories, Phys. Rev. B 67, 054505, (2003).
  * (51) I. Herbut, Modern Approach to Critical Phenomena (Cambridge University Press, Cambridge, England, 2007); Problems 3.8 and 3.10.
  * (52) B. I. Halperin, T. C. Lubensky, and S.-k. Ma, First-order phase transitions in superconductors and smectic-A liquid crystals, Phys. Rev. Lett. 32, 292 (1974).
  * (53) D. Nauenberg, D. J. Scalapino, Singularities and scaling functions at the Potts-model multicritical point, Phys. Rev. Lett. 44, 837 (1980).
  * (54) K. Kubota, H. Terao, Dynamical symmetry breaking in $QED_{3}$ from the Wilson RG point of view, Prog. Th. Phys. 105 809, (2001).
  * (55) I. F. Herbut, K. Kaveh, Chiral symmetry breaking in $QED_{3}$ in presence of irrelevant interactions: a renormalization group study, Phys. Rev. B 71, 184519 (2005).
  * (56) H. Gies and J. Jaeckel, Chiral phase structure of QCD with many flavors, Eur. Phys. J. C 46, 433 (2006).
  * (57) I. Herbut, Modern Approach to Critical Phenomena (Cambridge University Press, Cambridge, England, 2007); Chapter 4.3.
  * (58) D. B. Kaplan, J.-W. Lee, D. T. Son, and M. A. Stephanov, Conformality lost, Phys. Rev. D 80, 125005 (2009).
  * (59) I. F. Herbut, L. Janssen, Topological Mott insulator in three-dimensional systems with quadratic band touching, Phys. Rev. Lett. 113, 106401 (2014).
  * (60) I. F. Herbut, Chiral symmetry breaking in three-dimensional quantum electrodynamics as fixed point annihilation, Phys. Rev. D 94 025036 (2016).
  * (61) V. Gorbenko, S. Rychkov, B. Zan, Walking, weak first-order transitions, and complex CFTs, J. High. En. Phys. 108, (2018); Walking, weak first-order transitions, and Complex CFTs II. Two-dimensional Potts model at $Q>4$, SciPost Phys. 5, 050 (2018).
  * (62) B. Ihrig, N. Zerf, P. Marquard, I. F. Herbut, and M. M. Scherer, Abelian Higgs model at four loops, fixed-point collision, and deconfined criticality, Phys. Rev. B 100, 134507 (2019).
  * (63) A. F. Faedo, C. Hoyos, D. Mateos, and J. G. Subils, Holographic complex conformal field theories, Phys. Rev. Lett. 124, 161601 (2020).
  * (64) S. Okubo, Real representations of finite Clifford algebras 1: classification, J. Math. Phys. 32, 1657 (1991); Real representations of finite Clifford algebras 2: explicit construction and pseudo-octonion, J. Math Phys. 32, 1669 (1991).
  * (65) J. M. Kosterlitz, D. Nelson and M. E. Fisher, Bicritical and tetracritical points in anisotropic antiferromagnetic systems, Phys. Rev. B 13, 412 (1976).
  * (66) M. Hasenbusch and E. Vicari, Anisotropic perturbations in three-dimensional O(N)-symmetric vector models, Phys. Rev. B 84, 125136 (2011).
  * (67) F. Kos, D. Poland, D. Simmons-Duffin, and A. Vichi, Precision islands in the Ising and O(N) models, J. of High En. Phys. 2016, 36 (2016).
  * (68) H. Liu, E. Huffman, S. Chandrasekharan, R. K. Kaul, to appear.
  * Note (1) That is the two triplets transform as $(1,0)+(0,1)$ representation of some SO(3)$\times$ SO(3) subgroup of symmetry of the Dirac BdG Hamiltonian.
  * Note (2) In fact $N_{i}$ happen to be the three components of the Néel order parameter.
  * Note (3) They are actually standing for the three spin-components of the quantum spin-Hall insulator.