justified

# Gross-Neveu-Heisenberg criticality from $\boldsymbol{2+\epsilon}$ expansion

Konstantinos Ladovrechis  Shouryya Ray  Tobias Meng  Lukas Janssen  Institut
für Theoretische Physik and Würzburg-Dresden Cluster of Excellence ct.qmat, TU
Dresden, 01062 Dresden, Germany

###### Abstract

The Gross-Neveu-Heisenberg universality class describes a continuous quantum
phase transition between a Dirac semimetal and an antiferromagnetic insulator.
Such quantum critical points have originally been discussed in the context of
Hubbard models on $\pi$-flux and honeycomb lattices, but more recently also in
Bernal-stacked bilayer models, of potential relevance for bilayer graphene.
Here, we demonstrate how the critical behavior of this fermionic universality
class can be computed within an $\epsilon$ expansion around the lower critical
space-time dimension of two. This approach is complementary to the previously
studied expansion around the upper critical dimension of four. The crucial
technical novelty near the lower critical dimension is the presence of
different four-fermion interaction channels at the critical point, which we
take into account in a Fierz-complete way. By interpolating between the lower
and upper critical dimensions, we obtain improved estimates for the critical
exponents in 2+1 space-time dimensions. For the situation relevant to single-
layer graphene, we find an unusually small leading-correction-to-scaling
exponent, arising from the competition between different interaction channels.
This suggests that corrections to scaling may need to be taken into account
when comparing analytical estimates with numerical data from finite-size
extrapolations.

##  I Introduction

Fermionic quantum critical points are continuous quantum phase transitions
that are driven by interactions between gapless fermionic degrees of freedom.
They can be viewed as the simplest examples of quantum critical points that do
not exhibit classical analogs. Such transitions were originally discussed in
the context of toy models, mimicking aspects relevant to high-energy physics,
such as chiral symmetry breaking and spontaneous mass generation [1], non-
perturbative renormalizability [2, 3, 4], and asymptotic safety [5]. The
poster child of fermionic quantum criticality is embodied by the
$(2+1)$-dimensional Gross-Neveu-Ising transition, across which massless Dirac
fermions in two spatial dimensions acquire an interaction-induced mass gap as
a consequence of a spontaneous $\mathbb{Z}_{2}$ symmetry breaking [6, 7, 8, 9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. From a field-theoretical
viewpoint, a crucial simplicity of the Gross-Neveu-Ising transition is the
absence at criticality of any further four-fermion interaction channel at all
orders in perturbation theory [22]. This allows one to compute loop
corrections to high orders not only in the vicinity of the upper critical
space-time dimension of four [13], but also near the lower critical space-time
dimension of two [19].

In many physically relevant lattice realizations of fermion quantum
criticality, however, the symmetry that spontaneously breaks across the
transition is continuous, and governed by a vector order parameter. A well-
known example is given by the Hubbard model on the honeycomb lattice [23, 24,
25], which realizes as function of on-site repulsion $U$ a direct and
continuous transition between a symmetric Dirac semimetal at small $U$ and an
antiferromagnetic insulator at large $U$. In the strong-coupling phase, the
fermionic spectrum is gapped out and the SU(2) spin symmetry is spontaneously
broken. The transition is expected to fall into the Gross-Neveu-Heisenberg
universality class [26], which has been heavily discussed in recent years [27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 16, 38, 13, 39]. The transition
between nematic and coexistent nematic-antiferromagnetic orders on the Bernal-
stacked honeycomb bilayer has recently been identified as another potential
realization of Gross-Neveu-Heisenberg criticality, albeit with the number of
fermion degrees of freedom doubled in comparison with the single-layer case
[40]. Another possible realization on the Bernal-stacked honeycomb bilayer has
been proposed for the transition between the trigonal-warping-induced Dirac
semimetal and the antiferromagnetic insulator [41, 42]. In this latter case,
each quadratic band touching point present in the noninteracting limit for
vanishing trigonal warping splits into four Dirac cones, leading to a
quadrupled number of fermion degrees of freedom in comparison with the the
single-layer case [42].

Similar to classical universality classes, each fermionic quantum universality
class is characterized by a unique set of universal critical exponents. For
the relativistic Gross-Neveu-type criticalities, the dynamical critical
exponent $z=1$, leaving three independent quantities $1/\nu$, $\eta_{\phi}$,
and $\eta_{\psi}$, corresponding to the correlation-length exponent and the
order-parameter and fermion anomalous dimensions. While for the Gross-Neveu-
Ising case provisional convergence of the predictions of the different methods
appears within reach [43, 11], the disagreement among the different literature
results in case of the Gross-Neveu-Heisenberg criticality remains significant,
see, e.g., Ref. [29] for a recent overview. On the field-theoretical side, a
major challenge is the fact that the perturbative series are at best
asymptotically convergent, requiring appropriate resummation schemes. In the
Gross-Neveu-Ising case, a significant step forward has been the possibility to
employ interpolational schemes that make use of the known expansions near the
lower and upper critical dimensions simultaneously [16, 43]. In the Gross-
Neveu-Heisenberg case, however, different four-fermion interaction channels
are generically present at criticality, necessitating a Fierz-complete study
that deals with these channels in an unbiased way [44].

In this work, we provide such an analysis. We study the renormalization group
(RG) fixed-point structure of the theory space defined by the symmetries of
the Gross-Neveu-Heisenberg field theory within an expansion around the lower
critical space-time dimension of two. After Fierz reduction, this space is
spanned by six independent four-fermion couplings. We identify the fixed point
corresponding to Gross-Neveu-Heisenberg criticality and determine the
corresponding quantum critical behavior in terms of the correlation-length
exponent $1/\nu$ and the order-parameter anomalous dimensions $\eta_{\phi}$ to
one-loop order. The fermion anomalous dimension $\eta_{\psi}$ is computed to
two-loop order. To arrive at these results, we derive general formulae for the
order-parameter and fermion anomalous dimensions that we expect to be of use
for the community also beyond the particular theory space studied in this
work. Our results near the lower critical dimension allow us to provide
improved estimates for the exponents in the physical situation in 2+1 space-
time dimensions by employing an interpolational resummation scheme that takes
also previous results near the upper critical dimension [13] into account.

The remainder of this paper is organized as follows: In the next section, we
determine the theory space of the Gross-Neveu-Heisenberg model, its
symmetries, and a corresponding Fierz-complete basis. The RG flow and the
fixed-point structure are discussed in Sec. III. In Sec. IV, we derive the
general formulae to compute the order-parameter and fermion anomalous
dimensions $\eta_{\phi}$ and $\eta_{\psi}$ to one-loop and two-loop order,
respectively, and use these to provide a complete set of critical exponents
near the lower critical dimension. Estimates for the exponents in 2+1 space-
time dimensions using an interpolational resummation scheme are presented in
Sec. V. Section VI contains our conclusions and outlook.

##  II Gross-Neveu-Heisenberg theory space

###  II.1 Microscopic model

Within a purely fermionic formulation, the Gross-Neveu-Heisenberg model can be
defined via the Euclidean microscopic action [12, 26, 16, 39, 40]

$\displaystyle S_{\text{GNH}}$
$\displaystyle=\int\\!\mathrm{d}^{D}{x}\Bigl{\\{}\bar{\psi}^{\alpha}(\gamma_{\mu}\otimes\mathbb{1}_{2})\partial_{\mu}\psi^{\alpha}$
$\displaystyle\quad-\frac{{g}_{1}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\mathbb{1}_{2}\otimes\vec{\sigma})\psi^{\alpha}\right]^{2}\Bigr{\\}}$
(1)

in $D=2+\epsilon$ space-time dimensions. In the above equation, the space-time
index $\mu=0,\dots,D-1$, the flavor index $\alpha=1,\dots,N_{\mathrm{f}}$,
with $N_{\mathrm{f}}$ counting the number of four-component Dirac fermions,
and $\vec{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ denotes the vector of
$2\times 2$ Pauli matrices. Here and in what follows, if not stated otherwise,
summation over repeated indices is implicitly assumed. In two dimensions, we
employ an irreducible two-dimensional representation of the Clifford algebra
$\\{\gamma_{\mu},\gamma_{\nu}\\}=2\delta_{\mu\nu}\mathbb{1}_{2}$, such as

$\gamma_{0}=\begin{pmatrix}0&-\mathrm{i}\\\
\mathrm{i}&0\end{pmatrix}\qquad\text{and}\qquad\gamma_{1}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix}.$ (2)

The Dirac conjugate field is defined as
$\bar{\psi}\equiv\psi^{\dagger}(\gamma_{0}\otimes\mathbb{1}_{2})$. The Gross-
Neveu-Heisenberg coupling $g_{1}$ has mass dimension $[g_{1}]=2-D$. The
coupling becomes dimensionless, and the theory perturbatively renormalizable,
at $D=2$. Thus, $D=2$ defines a critical space-time dimension that can further
be identified as the lower critical dimension, around which we expand. This
approach is complementary to the previously studied expansion around the upper
critical dimension of four [12, 26, 13].

Within the large-$N_{\mathrm{f}}$ expansion, an ultraviolet completion of the
above model exists for all dimensions $2<D<4$ [39], and we assume this
property to hold also for finite $N_{\mathrm{f}}$. A potential complication is
the possibility of emergent evanescent operators that may in principle be
generated within the perturbative expansion [45, 46, 19]. At the one-loop
level, these induce shifts in the $\beta$ functions, which cancel with
corresponding contributions from the two-loop diagrams [47]. As for the
critical exponents, we therefore expect the emergence of evanescent operators
to play a role only beyond the leading-order estimates computed in the present
work.

In the (2+1)-dimensional realization of the model on the single-layer
honeycomb lattice [24, 26], the number of four-component fermion flavors is
$N_{\mathrm{f}}=2$. The criticality between nematic and coexistent nematic-
antiferromagnetic orders in bilayer graphene corresponds to $N_{\mathrm{f}}=4$
[40]. At the transition between the warping-induced spin-$1/2$ Dirac semimetal
and the antiferromagnetic insulator on the honeycomb bilayer, the number of
gapless four-component flavors is $N_{\mathrm{f}}=8$ [42].

###  II.2 Symmetries

In contrast to the Gross-Neveu-Ising case [19], the microscopic action defined
in Eq. (II.1) is not closed under RG transformations. Already at the one-loop
order, fluctuations induce new interaction channels that need to be taken into
account in a consistent way. The possible newly generated terms are, however,
strongly constrained by the symmetries of the microscopic action. These
symmetries are:

#### Lorentz invariance:

In two Euclidean space-time dimensions, the Dirac spinors transform as

$\psi(x)\mapsto
e^{-\mathrm{i}\varepsilon\frac{\mathrm{i}}{4}[\gamma_{0},\gamma_{1}]\otimes\mathbb{1}_{2}}\psi(x^{\prime}),$
(3)

where we have suppressed the flavor index for simplicity. The space-time
coordinate $x=(x^{\mu})$ transforms as
$x^{\mu}\mapsto{x^{\prime}}^{\mu}=\Lambda^{\mu}{}_{\nu}x^{\nu}$, with rotation
matrix $(\Lambda^{\mu}{}_{\nu})\in\mathrm{O}(2)$.

#### Flavor symmetry:

$\psi^{\alpha}\mapsto U^{\alpha\beta}\psi^{\beta},$ (4)

with unitary matrix $U\in\operatorname{U}(N_{\mathrm{f}})$.

####  $\mathbbm{Z}_{2}$ chiral symmetry:

$\psi\mapsto(\gamma_{5}\otimes\mathbb{1}_{2})\psi,$ (5)

where $\gamma_{5}=\mathrm{i}\gamma_{0}\gamma_{1}$ denotes the chiral matrix,
which is diagonal in the representation of Eq. (2). Note that
$\bar{\psi}\mapsto-\bar{\psi}(\gamma_{5}\otimes\mathbb{1}_{2})$, such that the
mass bilinears $\bar{\psi}\psi$ and
$\mathrm{i}\bar{\psi}(\gamma_{5}\otimes\mathbb{1}_{2})\psi$ are odd under
chiral symmetry.

#### SU(2) spin symmetry:

$\psi^{\alpha}\mapsto
e^{\mathrm{i}\phi\vec{n}\cdot(\mathbb{1}_{2}\otimes\vec{\sigma})}\psi^{\alpha},$
(6)

under which the Heisenberg bilinear
$\bar{\psi}(\mathbb{1}_{2}\otimes\vec{\sigma})\psi$ transforms as a vector.

#### Time-reversal symmetry:

$\psi(x)\mapsto\mathcal{T}\psi(x),$ (7)

with the time-reversal operator
$\mathcal{T}=(\gamma_{1}\otimes\sigma_{y})\mathcal{K}$ in Euclidean time,
where $\mathcal{K}$ denotes complex conjugation. We note that the Heisenberg
bilinear $\bar{\psi}(\mathbb{1}_{2}\otimes\vec{\sigma})\psi$ is odd under time
reversal, in agreement with its lattice realizations in 2+1 dimensions, where
it corresponds to antiferromagnetic orders [24, 26, 40]. The mass bilinears
$\bar{\psi}\psi$ and
$\mathrm{i}\bar{\psi}(\gamma_{5}\otimes\mathbb{1}_{2})\psi$ are even under
time reversal.

#### Inversion symmetry:

$\psi(x)\mapsto\mathcal{I}\psi(x^{\prime}),$ (8)

with the inversion operator $\mathcal{I}=\gamma_{0}\otimes\mathbb{1}_{2}$, and
the space-time coordinate $x=(x^{0},x^{1})$ transforming as $x\mapsto
x^{\prime}=(x^{0},-x^{1})$. While the standard mass term $\bar{\psi}\psi$ is
inversion symmetric, the bilinear $\mathrm{i}\bar{\psi}\gamma_{5}\psi$ is odd
under inversion.

###  II.3 Classification of four-fermion operators

The above symmetries forbid mass terms or other bilinears in the effective
action obtained by integrating out fermionic fluctuations. However, in
addition to the Gross-Neveu-Heisenberg four-fermion interaction already
present in the microscopic model, Eq. (II.1), there exist other four-fermion
terms that feature the same symmetries and thus will generically be generated
under the RG. In order to identify a basis of the theory space, we now
classify all possible four-fermion terms according to their symmetries.

Flavor symmetry allows for two different types of four-fermion terms: those
having singlet flavor structure
$(\bar{\psi}^{\alpha}\mathcal{O}\psi^{\alpha})(\bar{\psi}^{\beta}\mathcal{Q}\psi^{\beta})$
and those with non-singlet flavor structure
$(\bar{\psi}^{\alpha}\mathcal{O}\psi^{\beta})(\bar{\psi}^{\beta}\mathcal{Q}\psi^{\alpha})$.
Here, $\mathcal{O}$ and $\mathcal{Q}$ denote $4\times 4$ matrices that act on
the spinor indices of $\psi$ and $\bar{\psi}$. The singlet and non-singlet
terms are related to each other via Fierz identities [48, 22], and it is
therefore always possible to write the latter as a linear combination of the
former. Similarly, terms of the form
$(\psi^{\alpha\top}\mathcal{O}\psi^{\beta})(\bar{\psi}^{\beta}\mathcal{Q}\bar{\psi}^{\alpha\top})$
are related to the above terms via Fierz identities, and do not lead to any
new independent interaction channels. For our purposes, it thus suffices to
determine the invariant four-fermion terms with flavor singlet structure.
These can be constructed from the symmetry transformation properties of the
bilinears

$\bar{\psi}^{\alpha}\mathcal{O}\psi^{\alpha},$ (9)

with $\mathcal{O}$ being a $4\times 4$ matrix acting on the spinor indices of
$\psi$ and $\bar{\psi}$. A basis in the sixteen-dimensional space of $4\times
4$ matrices is given by the direct product of the basis matrices of the charge
and the spin sectors,

$\\{\mathbb{1}_{2},\gamma_{0},\gamma_{1},\gamma_{5}\\}\otimes\\{\mathbb{1}_{2},\vec{\sigma}\\},$
(10)

and any $\mathcal{O}$ can hence be written as a linear combination of these.

We now discuss the transformation properties of these basis matrices. They can
be divided into two groups of eight matrices each,
$A=\\{\mathbb{1}_{2},\gamma_{5}\\}\otimes\\{\mathbb{1}_{2},\vec{\sigma}\\}$
and $B=\\{\gamma_{0},\gamma_{1}\\}\otimes\\{\mathbb{1}_{2},\vec{\sigma}\\}$,
which commute and anticommute, respectively, with the chiral matrix
$\gamma_{5}$. Each group can be further split into sets of spin SU(2) scalars
and vectors, respectively, viz.,
$A_{\text{S}}=(\mathbb{1}_{2},\gamma_{5})\otimes\mathbb{1}_{2}$,
$A_{\text{V}}=(\mathbb{1}_{2},\gamma_{5})\otimes\vec{\sigma}$,
$B_{\text{S}}=(\gamma_{0},\gamma_{1})\otimes\mathbb{1}_{2}$, and
$B_{\text{V}}=(\gamma_{0},\gamma_{1})\otimes\vec{\sigma}$. Any flavor-singlet
four-fermion term invariant under both chiral and SU(2) spin symmetries can
therefore be written as
$\bar{\psi}^{\alpha}\mathcal{O}\psi^{\alpha}\bar{\psi}^{\beta}\mathcal{Q}\psi^{\beta}$
with $\mathcal{O}$ and $\mathcal{Q}$ from the same set $A_{\text{S}}$,
$A_{\text{V}}$, $B_{\text{S}}$, or $B_{\text{V}}$. Invariance under time
reversal and inversion then implies that $\mathcal{O}=\mathcal{Q}$. Finally,
Lorentz invariance implies that the two different four-fermion terms with
$\mathcal{O}=\mathcal{Q}\in B_{\text{S}}$ appear symmetrically in the
Lagrangian with the same coefficients, and equivalently for
$\mathcal{O}=\mathcal{Q}\in B_{\text{V}}$. Assuming $N_{\mathrm{f}}>1$, a
Fierz-complete basis of the Gross-Neveu-Heisenberg theory space therefore
contains six four-fermion terms,

$\displaystyle\mathcal{L}_{\text{int}}$
$\displaystyle=-\frac{g_{1}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\mathbb{1}_{2}\otimes\vec{\sigma})\psi^{\alpha}\right]^{2}-\frac{g_{2}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\gamma_{\mu}\otimes\vec{\sigma})\psi^{\alpha}\right]^{2}$
$\displaystyle\quad-\frac{g_{3}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\gamma_{5}\otimes\vec{\sigma})\psi^{\alpha}\right]^{2}-\frac{g_{4}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\mathbb{1}_{2}\otimes\mathbb{1}_{2})\psi^{\alpha}\right]^{2}$
$\displaystyle\quad-\frac{g_{5}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\gamma_{\mu}\otimes\mathbb{1}_{2})\psi^{\alpha}\right]^{2}-\frac{g_{6}}{2N_{\mathrm{f}}}\left[\bar{\psi}^{\alpha}(\gamma_{5}\otimes\mathbb{1}_{2})\psi^{\alpha}\right]^{2},$
(11)

parametrized by the six couplings $\boldsymbol{g}=(g_{1},\dots,g_{6})$.111For
$N_{\mathrm{f}}=1$, there exist further Fierz identities that may reduce the
number of indepdenent four-fermion terms. Here, $g_{1}$ corresponds to the
Gross-Neveu-Heisenberg coupling [12, 26, 16, 39, 40], $g_{4}$ corresponds to
the Gross-Neveu-Ising coupling [3, 4, 6, 5, 22, 19], and $g_{5}$ corresponds
to the Thirring coupling [47, 48, 49]. In following, we study the flow of the
full effective action

$S=\int\\!\mathrm{d}^{D}x\left[\bar{\psi}^{\alpha}(\gamma_{\mu}\otimes\mathbb{1}_{2})\partial_{\mu}\psi^{\alpha}+\mathcal{L}_{\text{int}}\right],$
(12)

out of which $S_{\text{GNH}}$ is a subspace that we explicitly show to be not
closed under the RG.

##  III Renormalization group flow

###  III.1 Flow equations

In order to compute the RG flow equations for the six four-fermion couplings
$\boldsymbol{g}=(g_{1},\dots,g_{6})$, we employ the general one-loop formula
derived in [22]. A straightforward evaluation of the matrix algebra occurring
in this formula, using standard computer algebra software, yields the
following flow equations, valid for arbitrary $N_{\mathrm{f}}$,

$\displaystyle\beta_{1}$
$\displaystyle=\left[\epsilon+\frac{2(2g_{2}-g_{3}+g_{4}+2g_{5}+g_{6})}{N_{\mathrm{f}}}\right]g_{1}$
$\displaystyle\quad-\frac{2(2N_{\mathrm{f}}+1)}{N_{\mathrm{f}}}g^{2}_{1}+\frac{4(g_{3}g_{5}+g_{2}g_{6})}{N_{\mathrm{f}}},$
(13) $\displaystyle\beta_{2}$ $\displaystyle=\epsilon
g_{2}+\frac{8g^{2}_{2}}{N_{\mathrm{f}}}+\frac{2[g_{3}(g_{3}+g_{4})+g_{1}(g_{1}+g_{6})]}{N_{\mathrm{f}}},$
(14) $\displaystyle\beta_{3}$
$\displaystyle=\left[\epsilon+\frac{2(g_{1}+2g_{2}-g_{4}+2g_{5}-g_{6})}{N_{\mathrm{f}}}\right]g_{3}$
$\displaystyle\quad+\frac{2(2N_{\mathrm{f}}+1)}{N_{\mathrm{f}}}g_{3}^{2}+\frac{4(g_{2}g_{4}+g_{1}g_{5})}{N_{\mathrm{f}}},$
(15) $\displaystyle\beta_{4}$
$\displaystyle=\left[\epsilon+\frac{2(3g_{1}+6g_{2}+3g_{3}+2g_{5}+g_{6})}{N_{\mathrm{f}}}\right]g_{4}$
$\displaystyle\quad-\frac{2(2N_{\mathrm{f}}-1)}{N_{\mathrm{f}}}g_{4}^{2}+\frac{4(3g_{2}g_{3}+g_{5}g_{6})}{N_{\mathrm{f}}},$
(16) $\displaystyle\beta_{5}$ $\displaystyle=\epsilon
g_{5}+\frac{2(3g_{1}g_{3}+g_{4}g_{6})}{N_{\mathrm{f}}},$ (17)
$\displaystyle\beta_{6}$
$\displaystyle=\left[\epsilon-\frac{2(3g_{1}-6g_{2}+3g_{3}+g_{4}-2g_{5})}{N_{\mathrm{f}}}\right]g_{6}$
$\displaystyle\quad+\frac{2(2N_{\mathrm{f}}-1)}{N_{\mathrm{f}}}g_{6}^{2}+\frac{4(3g_{1}g_{2}+g_{4}g_{5})}{N_{\mathrm{f}}},$
(18)

where $\epsilon=D-2$ and we have rescaled the couplings as
$\ell_{1}^{\text{(F)}}g_{i}/2\pi\mapsto g_{i}$ for $i=1,\dots,6$, with
$\ell_{1}^{\text{(F)}}$ a dimensionless regulator-dependent constant. The sign
of the $\beta$ functions is defined such that a coupling $g_{i}$ decreases
(increases) in the flow towards the infrared if $\beta_{i}>0$ ($\beta_{i}<0$).
For $g_{1}=g_{2}=g_{3}=0$, the above flow equations are consistent with those
of Ref. [47] upon identifying $\ell^{\text{(F)}}_{1}=1$ for the minimal
subtraction scheme employed therein. We note that the flow equations (13)–(18)
are invariant under the exchange of the couplings as
$(g_{1},g_{2},g_{3},g_{4},g_{5},g_{6})\leftrightarrow(-g_{3},g_{2},-g_{1},-g_{6},g_{5},-g_{4})$.
This property arises from the fact that the chiral matrix $\gamma_{5}$
anticommutes with the fermion propagator and squares to one.

###  III.2 Large-Nf fixed-point structure

The topology of the RG flow is determined by the solutions of the fixed-point
equations $\partial_{t}g_{i}|_{\boldsymbol{g}^{\star}}=0$. Any nontrivial
solution $\boldsymbol{g}^{\star}\neq 0$ is located at finite distance of order
$\mathcal{O}(\epsilon)$ from the Gaussian fixed point
$\boldsymbol{g}^{\star}=0$. While the Gaussian fixed point is fully infrared
stable for $\epsilon>0$, all interacting fixed points have at least one
infrared relevant direction. For initial values of the couplings beyond a
certain finite threshold of order $\mathcal{O}(\epsilon)$, the flow is
unstable and diverges at finite RG time, indicating the onset of spontaneous
symmetry breaking. Critical fixed points, which govern the quantum critical
behavior of continuous phase transitions, are those with precisely one
infrared relevant direction. Among the different critical fixed points, we are
looking for the one corresponding to the onset of spin symmetry breaking via
condensation of the SU(2) vector

$\langle\bar{\psi}^{\alpha}(\mathbbm{1}_{2}\otimes\vec{\sigma})\psi^{\alpha}\rangle\neq
0.$ (19)

This fixed point is readily identified in the large-$N_{\mathrm{f}}$ limit. In
this limit, the individual flow equations for $g_{1},g_{2},\dots,g_{6}$
decouple and the fixed-point structure can be computed analytically. For the
set of six quadratic fixed-point equations, there may be at most $2^{6}=64$
possibly degenerate and/or complex solutions. However, in the
large-$N_{\mathrm{f}}$ limit, the quadratic term $\propto g_{2}^{2}$ vanishes
in the flow equation for $g_{2}$. The same is true, albeit also beyond the
large-$N_{\mathrm{f}}$ limit, for the quadratic term $\propto g_{5}^{2}$ in
the flow equation of the Thirring coupling $g_{5}$ [47]. The absence of such
quadratic term can be understood as a fixed point located at infinite coupling
$g_{2}^{\star}\to\infty$ and $g_{5}^{\star}\to\infty$, respectively. In fact,
new fixed points located at $g_{2}^{\star}\propto N_{\mathrm{f}}$ emerge upon
the inclusion of the $\mathcal{O}(1/N_{\mathrm{f}})$ corrections in the flow
equation for $g_{2}$. However, the Thirring coupling $g_{5}^{\star}$ vanishes
at any fixed point also beyond the large-$N_{\mathrm{f}}$ limit, as a
consequence of the absence of a $g_{5}^{2}$ term also for finite
$N_{\mathrm{f}}$. For large, but finite, $N_{\mathrm{f}}$, in addition to the
Gaussian fixed point at vanishing couplings, there are therefore $2^{5}-1=31$
possibly degenerate and/or complex fixed points at finite couplings. In the
large-$N_{\mathrm{f}}$ limit, the critical fixed points are located on the
coordinate axes of our flavor-singlet basis. The symmetry breaking pattern
corresponding to each of the critical fixed points can then be simply
identified on the basis of the mean-field decoupling of the corresponding
flavor-singlet four-fermion term, which is controlled in the
large-$N_{\mathrm{f}}$ limit. This implies that the fixed point located at

$\boldsymbol{g}^{\star}_{\text{GNH}}=[1/4+\mathcal{O}(1/N_{\mathrm{f}}),\mathcal{O}(1/N_{\mathrm{f}}),0,0,0,\mathcal{O}(1/N_{\mathrm{f}}^{2})]\epsilon$
(20)

corresponds to Gross-Neveu-Heisenberg quantum criticality, with SU(2) order
parameter
$\langle\bar{\psi}^{\alpha}(\mathbb{1}_{2}\otimes\vec{\sigma})\psi^{\alpha}\rangle$,
and the fixed point located at

$\boldsymbol{g}^{\star}_{\text{GNI}}=[0,0,0,1/4+\mathcal{O}(1/N_{\mathrm{f}}),0,0]\epsilon$
(21)

corresponds to Gross-Neveu-Ising quantum criticality, with $\mathbbm{Z}_{2}$
order parameter
$\langle\bar{\psi}^{\alpha}(\mathbb{1}_{2}\otimes\mathbb{1}_{2})\psi^{\alpha}\rangle$.
Another fixed point located at

$\boldsymbol{g}^{\star}_{\text{GNI}^{\prime}}=[0,0,0,0,0,1/4+\mathcal{O}(1/N_{\mathrm{f}})]\epsilon$
(22)

also corresponds to Gross-Neveu-Ising quantum criticality, with
$\mathbbm{Z}_{2}$ order parameter
$\langle\bar{\psi}^{\alpha}(\gamma_{5}\otimes\mathbb{1}_{2})\psi^{\alpha}\rangle$,
breaking not only chiral symmetry but also inversion symmetry. This fixed
point is completely equivalent to the one at
$\boldsymbol{g}_{\text{GNI}}^{\star}$ due to the invariance of the flow
equations under exchange of the couplings as
$(g_{1},g_{2},g_{3},g_{4},g_{5},g_{6})\leftrightarrow(-g_{3},g_{2},-g_{1},-g_{6},g_{5},-g_{4})$.
In particular, the set of eigenvalues of the stability matrix at
$\boldsymbol{g}_{\text{GNI}}^{\star}$ and
$\boldsymbol{g}_{\text{GNI}^{\prime}}^{\star}$ are equal. This statement holds
for arbitrary $N_{\mathrm{f}}$ within the one-loop approximation.

Note that the only flow equation containing a term $\propto g_{4}^{2}$ is the
one for $g_{4}$ itself, such that upon starting the flow on the axis
$\boldsymbol{g}\propto\boldsymbol{g}^{\star}_{\text{GNI}}$ in the ultraviolet,
no other couplings $g_{i\neq 4}$ are generated in the infrared, in agreement
with earlier work [47, 19]. Similarly, the flow equations for $g_{3}$,
$g_{4}$, and $g_{5}$ do not contain terms $\propto g_{i}g_{j}$ with
$i,j\in\\{1,2,6\\}$, such that $g_{3}$, $g_{4}$, and $g_{5}$ are not generated
under the RG if simultaneously absent initially. This is again a statement
that holds for arbitrary $N_{\mathrm{f}}$ within the one-loop approximation.
Put differently, the space spanned by the couplings $g_{1}$, $g_{2}$, and
$g_{6}$ is invariant under the one-loop RG. As this subspace of the full
theory space contains the Gross-Neveu-Heisenberg fixed point, we denote it as
“Gross-Neveu-Heisenberg subspace”. For the analysis of the Gross-Neveu-
Heisenberg criticality at the one-loop order, it is hence sufficient to
consider the flow only within this subspace. We emphasize, however, that the
RG invariance of the Gross-Neveu-Heisenberg subspace is _not_ symmetry
protected. At higher loop orders, there may be other terms that will require
one to consider the full six-dimensional theory space.

Table 1: Fixed points in Gross-Neveu-Heisenberg subspace, their locations, number of relevant directions, and collisions as function of $N_{\mathrm{f}}$. | $(g_{1}^{\star},g_{2}^{\star},g_{6}^{\star})/\epsilon$ | $\\#(\Theta_{i}>0)$ |
---|---|---|---
$\mathcal{O}$ | $(0,0,0)$ | 0 |  Gaussian
$\mathcal{H}$ | $(h_{1},h_{2},h_{6})$ | 1 |  Gross-Neveu-Heisenberg
$\mathcal{I}$ | $(0,0,-\frac{N_{\mathrm{f}}}{4N_{\mathrm{f}}-2})$ | $\begin{cases}1,&\text{for }N_{\mathrm{f}}>N_{\mathrm{f}}^{(1)}\\\ 2,&\text{for }N_{\mathrm{f}}<N_{\mathrm{f}}^{(1)}\\\ \end{cases}$ |  Gross-Neveu-Ising′, collides with $\mathcal{D}$ for $N_{\mathrm{f}}\to N_{\mathrm{f}}^{(1)}={3}/{2}$
$\mathcal{A}$ | $(0,-\frac{N_{\mathrm{f}}}{8},0)$ | 2 |  part of RG invariant plane spanned by $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$
$\mathcal{B}$ | $(\frac{N_{\mathrm{f}}}{4N_{\mathrm{f}}+4},0,-\frac{N_{\mathrm{f}}}{4N_{\mathrm{f}}+4})$ | 2 |  part of RG invariant plane spanned by $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$
$\mathcal{C}$ | $(\frac{N_{\mathrm{f}}}{4N_{\mathrm{f}}+4},-\frac{N_{\mathrm{f}}}{8},-\frac{N_{\mathrm{f}}}{4N_{\mathrm{f}}+4})$ | 3 |  part of RG invariant plane spanned by $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$
$\mathcal{D}$ | $(e_{1},e_{2},e_{6})$ | $\begin{cases}\text{complex},&\text{for }N_{\mathrm{f}}>N_{\mathrm{f}}^{(2)}\\\ 2,&\text{for }N_{\mathrm{f}}^{(1)}<N_{\mathrm{f}}<N_{\mathrm{f}}^{(2)}\\\ 1,&\text{for }N_{\mathrm{f}}<N_{\mathrm{f}}^{(1)}\\\ \end{cases}$ |  annihilates with $\mathcal{E}$ for $N_{\mathrm{f}}\nearrow N_{\mathrm{f}}^{(2)}=1.5146\dots$, collides with $\mathcal{I}$ for $N_{\mathrm{f}}\to N_{\mathrm{f}}^{(1)}={3}/{2}$
$\mathcal{E}$ | $(f_{1},f_{2},f_{6})$ | $\begin{cases}\text{complex},&\text{for }N_{\mathrm{f}}>N_{\mathrm{f}}^{(2)}\\\ 1,&\text{for }N_{\mathrm{f}}<N_{\mathrm{f}}^{(2)}\\\ \end{cases}$ |  annihilates with $\mathcal{D}$ for $N_{\mathrm{f}}\nearrow N_{\mathrm{f}}^{(2)}=1.5146\dots$

###  III.3 Gross-Neveu-Heisenberg subspace

The Gross-Neveu-Heisenberg subspace is defined by

$\boldsymbol{g}=(g_{1},g_{2},0,0,0,g_{6}),$ (23)

which is the smallest RG invariant subspace containing the Gross-Neveu-
Heisenberg fixed point. The latter is located at
$\boldsymbol{g}_{\text{GNH}}^{\star}=[h_{1}(N_{\mathrm{f}}),h_{2}(N_{\mathrm{f}}),0,0,0,h_{6}(N_{\mathrm{f}})]\epsilon$,
with real functions $h_{1}(N_{\mathrm{f}})>0$, $h_{2}(N_{\mathrm{f}})<0$, and
$h_{6}(N_{\mathrm{f}})>0$. The Gross-Neveu-Heisenberg subspace also contains
the Gross-Neveu-Ising′ fixed point at
$\boldsymbol{g}_{\text{GNI}^{\prime}}^{\star}=[0,0,0,0,0,-N_{\mathrm{f}}/(4N_{\mathrm{f}}-2)]\epsilon$,
as well as six additional real or complex fixed points. Table 1 shows the
locations of all fixed points in the Gross-Neveu-Heisenberg subspace, together
with their numbers of infrared relevant directions. The latter are obtained
from the eigenvalues $\Theta_{i}$ of the stability matrix
$(-\partial\beta_{i}/\partial g_{j})$ at the respective fixed point.
Importantly, for all values of $N_{\mathrm{f}}$, the Gross-Neveu-Heisenberg
fixed point features a single infrared relevant direction, corresponding to a
critical fixed point. For finite $N_{\mathrm{f}}<\infty$, it features finite
fixed-point couplings in all three channels
$-g_{1}/(2N_{\mathrm{f}})[\bar{\psi}(\mathbb{1}\otimes\vec{\sigma})\psi]^{2}$,
$-g_{2}/(2N_{\mathrm{f}})[\bar{\psi}(\gamma_{\mu}\otimes\vec{\sigma})\psi]^{2}$,
and
$-g_{6}/(2N_{\mathrm{f}})[\bar{\psi}(\gamma_{5}\otimes\mathbb{1}_{2})\psi]^{2}$
of the Gross-Neveu-Heisenberg subspace. The evolution of these fixed-point
couplings as function of $N_{\mathrm{f}}$ is depicted in Fig. 1.

Figure 1: Evolution of Gross-Neveu-Heisenberg fixed-point couplings
$\boldsymbol{g}^{\star}_{\text{GNH}}=[h_{1}(N_{\mathrm{f}}),h_{2}(N_{\mathrm{f}}),0,0,0,h_{6}(N_{\mathrm{f}})]\epsilon$
as function of flavor number $N_{\mathrm{f}}$.

We have explicitly verified that perturbations out of the Gross-Neveu-
Heisenberg subspace are infrared irrelevant in the vicinity of the Gross-
Neveu-Heisenberg fixed point. As an illustration of the RG flow near the
Gross-Neveu-Heisenberg fixed point, Fig. 2 presents the flow diagram within
the plane spanned by the couplings $g_{1}$ and $g_{6}$ for fixed $g_{2}\equiv
h_{2}(N_{\mathrm{f}})$, for different values of $N_{\mathrm{f}}$. Therein, the
Gross-Neveu-Heisenberg fixed point labeled by $\mathcal{H}$ is marked as red
dot. The gray dots labeled by $\mathcal{O}^{\prime}$, $\mathcal{B}^{\prime}$,
and $\mathcal{I}^{\prime}$ indicate points in parameter space in which the
flow is perpendicular to the plane $g_{2}\equiv h_{2}(N_{\mathrm{f}})$. In the
large-$N_{\mathrm{f}}$ limit, they represent projections of the fixed points
$\mathcal{O}$, $\mathcal{B}$, and $\mathcal{I}$, respectively, and are
adiabatically connected to these upon lowering $N_{\mathrm{f}}$.

(a) $N_{\mathrm{f}}=2$ (b) $N_{\mathrm{f}}=4$ (c) $N_{\mathrm{f}}=8$

Figure 2: RG flow in the plane spanned by $g_{1}$ and $g_{6}$ for fixed
$g_{2}=h_{2}(N_{\mathrm{f}})$ through the Gross-Neveu-Heisenberg fixed point
for (a) $N_{\mathrm{f}}=2$, (b) $N_{\mathrm{f}}=4$, and (c)
$N_{\mathrm{f}}=8$. The red dot denotes the position of the Gross-Neveu-
Heisenberg fixed point $\mathcal{H}$. Gray dots labeled by
$\mathcal{O}^{\prime}$, $\mathcal{B}^{\prime}$, and $\mathcal{I}^{\prime}$
indicate points in parameter space in which the flow is perpendicular to the
plane $g_{2}=h_{2}(N_{\mathrm{f}})$, which become projections of the Gaussian
fixed point $\mathcal{O}$, the critical fixed point $\mathcal{B}$, and the
bicritical fixed point $\mathcal{I}$ in the large-$N_{\mathrm{f}}$ limit. The
gray line indicates the projection of the RG invariant subspace spanned by the
fixed points $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$.

The Gross-Neveu-Ising′ fixed point at
$\boldsymbol{g}_{\text{GNI}^{\prime}}^{\star}$ features a single relevant
direction only for $N_{\mathrm{f}}$ above a critical flavor number
$N_{\mathrm{f}}^{(1)}=3/2+\mathcal{O}(\epsilon)$. For $N_{\mathrm{f}}\to
N_{\mathrm{f}}^{(1)}$, it collides with the bicritical fixed point
$\mathcal{D}$, and exchanges, for $N_{\mathrm{f}}<N_{\mathrm{f}}^{(1)}$, its
role with respect to RG stability with the latter. Note that due to the
symmetry of the flow equations, a simultaneous fixed-point collision occurs
away from the Gross-Neveu-Heisenberg subspace at the same flavor number,
involving the Gross-Neveu-Ising fixed point
$\boldsymbol{g}_{\text{GNI}}^{\star}=[0,0,0,N_{\mathrm{f}}/(4N_{\mathrm{f}}-2),0,0]\epsilon$
and a bicritical fixed point at $[0,e_{2},-e_{1},-e_{6},0,0]\epsilon$. Such
fixed-point-collision scenario, involving an exchange of stability between
Gross-Neveu-Ising and bicritical fixed points, has been found previously in a
one-loop analysis in fixed $D=2+1$ space-time dimensions [22].

Remarkably, for a second critical flavor number
$N_{\mathrm{f}}^{(2)}=1.5146\dots$, i.e., only slightly above
$N_{\mathrm{f}}^{(1)}$, the fixed point $\mathcal{D}$ is involved in another
fixed-point collision. In this case, it merges with the fixed point
$\mathcal{E}$, with both of them disappearing into the complex coupling plane
for $N_{\mathrm{f}}>N_{\mathrm{f}}^{(2)}$. Such fixed-point annihilation has
been observed in a variety of gauge theories in 2+1 [50, 51, 52, 53, 54, 55,
56, 57, 58] and higher [59, 60] dimensions, but also in non-gauge theories
[22, 61, 62, 63, 64, 65, 66, 67, 68].

##  IV Critical exponents

The universal critical exponents we determine here are the correlation-length
exponent $1/\nu$, the anomalous dimensions of the order-parameter and fermion
fields, $\eta_{\phi}$ and $\eta_{\psi}$, respectively, and the corrections-to-
scaling exponent $\omega$. The dynamical critical exponent is $z=1$ exactly,
as a consequence of Lorentz invariance, which has been shown to emerge at low
energy for a Gross-Neveu-Heisenberg quantum critical point [40].

###  IV.1 Correlation-length exponent $\boldsymbol{1/\nu}$

The correlation-length exponent $1/\nu$ determines the divergence of the
correlation length near the quantum critical point. It is given by the unique
positive eigenvalue $\Theta_{1}>0$ of the stability matrix
$(-\partial\beta_{i}/\partial g_{j})$ at the corresponding critical fixed
point. To leading order in the perturbative expansion, we find both for the
Gross-Neveu-Heisenberg and the Gross-Neveu Ising fixed points

$1/\nu=\epsilon+\mathcal{O}(\epsilon),$ (24)

in agreement with the general result valid for all critical four-fermion
models near the lower critical dimension [22].

###  IV.2 Order-parameter anomalous dimension $\boldsymbol{\eta_{\phi}}$

As there appears no dangerously irrelevant coupling in the problem, we assume
hyperscaling to hold. The order-parameter anomalous dimension $\eta_{\phi}$ is
then linked to the correlation-length exponent $1/\nu$ and the susceptibility
exponent $\gamma$ via the hyperscaling relation

$\eta_{\phi}=2-\gamma/\nu.$ (25)

Within our fermionic formulation, we can determine the susceptibility exponent
$\gamma$ using the scheme described in Ref. [56]. To this end, we add the
corresponding infinitesimal mass term to the effective Lagrangian as

$\displaystyle\mathcal{L}\mapsto\mathcal{L}+\Delta\,\bar{\psi}^{\alpha}\mathcal{M}\psi^{\alpha},$
(26)

with $\mathcal{M}=\mathbb{1}_{2}\otimes\mathbb{1}_{2}$ for Gross-Neveu-Ising
criticality and $\mathcal{M}=\mathbb{1}_{2}\otimes\vec{\sigma}$ for Gross-
Neveu-Heisenberg criticality. In the presence of the infinitesimal mass term,
the scaling form of the free energy density near criticality reads [56]

$f(\delta\boldsymbol{g},\Delta)=|\delta\boldsymbol{g}|^{D\nu}\mathcal{F}^{\pm}\left(\frac{\Delta}{|\delta\boldsymbol{g}|^{x\nu}}\right),$
(27)

with scaling function $\mathcal{F}^{\pm}$. In the above equation,
$\delta\boldsymbol{g}$ is the eigenvector associated with the RG relevant
direction, and $x$ denotes the eigenvalue associated with the RG flow of
$\Delta$,

$\beta_{\Delta}=-x\Delta+\mathcal{O}(\Delta^{2}).$ (28)

Differentiating twice with respect to the mass parameter $\Delta$ yields the
scaling of the susceptibility,

$\chi=-\frac{\partial^{2}f}{\partial\Delta^{2}}\propto|\delta\boldsymbol{g}|^{-\gamma},\quad\text{with}\quad\gamma=(2x-D)\nu.$
(29)

With the help of the hyperscaling relation (25), the order-parameter anomalous
dimension is then given by

$\eta_{\phi}=D+2(1-x).$ (30)

At one-loop order, the flow of the mass parameter $\Delta$ has the form

$\beta_{\Delta}=-\left(1+\sum_{i}c_{i}g_{i}\right)\Delta+\mathcal{O}(\Delta^{2}),$
(31)

with coefficients

$\displaystyle c_{i}$
$\displaystyle=\frac{1}{4DN_{\mathrm{f}}}\sum_{\mu}\Bigl{\\{}N_{\mathrm{f}}\operatorname{Tr}(\mathcal{M}\gamma_{\mu}\mathcal{O}_{i}\gamma_{\mu})\operatorname{Tr}(\mathcal{M}\mathcal{O}_{i})$
$\displaystyle\quad-\operatorname{Tr}\left[\mathcal{O}_{i}\gamma_{\mu}\mathcal{M}\gamma_{\mu}\mathcal{O}_{i}\mathcal{M}\right]\Bigr{\\}},$
(32)

where $\mathcal{O}_{i}$ denotes the $4\times 4$ matrix in the four-fermion
term parametrized by $g_{i}$, and for brevity we have omitted factors of
$\mathbb{1}_{2}$ in direct products with $\gamma$ matrices, i.e.,
$\gamma_{\mu}\equiv\gamma_{\mu}\otimes\mathbb{1}_{2}$. Note that in the above
equation, no summation over repeated indices $i$ is assumed, and we have
rescaled the couplings in the same way as described below Eqs. (13)–(18). The
order-parameter anomalous dimension at a fixed point
$\boldsymbol{g}^{\star}=(g_{i}^{\star})$ can then be obtained from

$\eta_{\phi}=2+\epsilon-2\sum_{i}c_{i}g_{i}^{\star}$ (33)

in $D=2+\epsilon$ dimensions.

Evaluating the matrix algebra for the Gross-Neveu-Ising mass
$\mathcal{M}=\mathbb{1}_{2}\otimes\mathbb{1}_{2}$ and
$\mathcal{O}_{4}=\mathbb{1}_{2}\otimes\mathbb{1}_{2}$, using the Gross-Neveu-
Ising fixed-point value
$\boldsymbol{g}_{\text{GNI}}^{\star}=[0,0,0,N_{\mathrm{f}}/(4N_{\mathrm{f}}-2),0,0]\epsilon$,
yields the order-parameter anomalous dimension for Gross-Neveu-Ising
criticality

$\eta_{\phi}^{\text{GNI}}=2-\frac{2N_{\mathrm{f}}}{2N_{\mathrm{f}}-1}\epsilon+\mathcal{O}(\epsilon^{2}).$
(34)

This agrees with the known results near the lower critical dimension [69, 70,
71, 19], thereby providing a first cross-check of our calculations.222Note
that the definition for $N$ used in Ref. [19] deviates from our definition for
$N_{\mathrm{f}}$ as
$N^{\text{(Ref.~{}\cite[cite]{[\@@bibref{Number}{gracey16}{}{}]})}}=2N_{\mathrm{f}}^{\text{(this
work)}}$.

For Gross-Neveu-Heisenberg criticality, we assume
$\mathcal{M}=\mathbb{1}_{2}\otimes\sigma_{z}$ without loss of generality, and
use the couplings
$\boldsymbol{g}_{\text{GNH}}^{\star}=[h_{1}(N_{\mathrm{f}}),h_{2}(N_{\mathrm{f}}),0,0,0,h_{6}(N_{\mathrm{f}})]\epsilon$
at the Gross-Neveu-Heisenberg fixed point. Evaluating the matrix algebra
yields

$\displaystyle\eta_{\phi}^{\text{GNH}}$
$\displaystyle=2+\left[1-2\frac{\left(1+4N_{\mathrm{f}}\right)h_{1}(N_{\mathrm{f}})+2h_{2}(N_{\mathrm{f}})-h_{6}(N_{\mathrm{f}})}{N_{\mathrm{f}}}\right]\epsilon$
$\displaystyle\quad+\mathcal{O}(\epsilon^{2})$ (35)

for general $N_{\mathrm{f}}$. It is instructive to further expand our
small-$\epsilon$ results for large $N_{\mathrm{f}}$,

$\eta_{\phi}^{\text{GNH}}=2-\left[1-\frac{1}{2N_{\mathrm{f}}}-\frac{5}{4N_{\mathrm{f}}^{2}}+\frac{19}{8N_{\mathrm{f}}^{3}}+\mathcal{O}(1/N_{\mathrm{f}}^{4})\right]\epsilon+\mathcal{O}(\epsilon^{2}),$
(36)

which agrees, up to the order calculated, with the large-$N_{\mathrm{f}}$
exponents computed for arbitrary $2<D<4$ [39], upon expanding the latter for
small $\epsilon=D-2$.333Note that the definitions for $\epsilon$ and $N$ used
in Ref. [39] deviate from our definitions for $\epsilon$ and $N_{\mathrm{f}}$
as
$\epsilon^{\text{(Ref.~{}\cite[cite]{[\@@bibref{Number}{gracey18a}{}{}]})}}=-\epsilon^{\text{(this
work)}}/2$ and
$N^{\text{(Ref.~{}\cite[cite]{[\@@bibref{Number}{gracey18a}{}{}]})}}=2N_{\mathrm{f}}^{\text{(this
work)}}$. This furnishes another nontrivial cross-check of our calculations.
For the cases relevant for interacting electrons on the single-layer [24, 26,
25, 30] and bilayer [40, 41, 42] honeycomb lattices, we explicitly find from
Eq. (IV.2), i.e., without expanding in $1/N_{\mathrm{f}}$,

$\displaystyle\eta_{\phi}^{\text{GNH}}=\begin{cases}2-\frac{\epsilon}{6}+\mathcal{O}(\epsilon^{2}),&\text{for
$N_{\mathrm{f}}=2$},\\\
2-0.812333\epsilon+\mathcal{O}(\epsilon^{2}),&\text{for
$N_{\mathrm{f}}=4$},\\\
2-0.921305\epsilon+\mathcal{O}(\epsilon^{2}),&\text{for
$N_{\mathrm{f}}=8$}.\end{cases}$ (37)

Equation (37) represents one of the main results of this work.

###  IV.3 Fermion anomalous dimension $\boldsymbol{\eta_{\psi}}$

(a) Gross-Neveu-Ising (b) Gross-Neveu-Heisenberg

Figure 3: (a) Corrections-to-scaling exponent $\omega$ for the Gross-Neveu-
Ising fixed point as a function of $N_{\mathrm{f}}$. (b) Same as (a), but for
the Gross-Neveu-Heisenberg fixed point, featuring a distinct minimum near
$N_{\mathrm{f}}=2$, corresponding to a slow flow towards the fixed point. The
inset shows the angle $\varphi$ between $\boldsymbol{g}^{\star}_{\text{GNH}}$
and the surface normal $\boldsymbol{n}=(1,0,1)/\sqrt{2}$ of the invariant
subspace spanned by the fixed points $\mathcal{A}$, $\mathcal{B}$, and
$\mathcal{C}$, as function of $N_{\mathrm{f}}$.

As in any critical four-fermion model near the lower critical dimension, the
fermion anomalous dimension vanishes at one-loop order,
$\eta_{\psi}=0+\mathcal{O}(\epsilon^{2})$. This implies that knowledge of the
$\mathcal{O}(\epsilon)$ fixed-point values, together with the result of the
corresponding two-loop selfenergy diagram, is sufficient to compute
$\eta_{\psi}$ to order $\mathcal{O}(\epsilon^{2})$. While at the one-loop
order all regulator dependences can be factored out by appropriate rescalings
of the couplings, this may no longer be true at higher orders. We employ a
minimal subtraction scheme analogous to Ref. [47], with an infrared cutoff in
the form of a mass term
$m\bar{\psi}(\mathbb{1}_{2}\otimes\mathbb{1}_{2})\psi$, and an effective
fermion propagator

$G(p)=\frac{-\mathrm{i}p_{\mu}(\gamma_{\mu}\otimes\mathbb{1}_{2})}{p^{2}+m^{2}}.$
(38)

Note that we have omitted the mass term in the numerator of the effective
propagator, which gives no contribution to the pole in $1/\epsilon$, as a
consequence of the infrared finiteness of the theory [47]. Evaluating the
sunset diagram for the fermion selfenergy at a fixed point
$\boldsymbol{g}^{\star}=(g_{i}^{\star})$ yields

$\displaystyle\eta_{\psi}=g_{i}^{\star}H_{ij}g_{j}^{\star},$ (39)

with the matrix elements

$\displaystyle H_{ij}$
$\displaystyle=\dfrac{1}{32N_{\mathrm{f}}^{2}}\sum_{\mu,\nu,\lambda}(\delta_{\mu,0}\delta_{\nu\lambda}+\delta_{\nu,0}\delta_{\mu\lambda}+\delta_{\lambda,0}\delta_{\mu\nu})\bigl{\\{}$
$\displaystyle\quad
N_{\mathrm{f}}\operatorname{Tr}(\gamma_{0}\mathcal{O}_{i}\gamma_{\mu}\mathcal{O}_{j})\operatorname{Tr}(\gamma_{\nu}\mathcal{O}_{j}\gamma_{\lambda}\mathcal{O}_{i})$
$\displaystyle\quad-\operatorname{Tr}(\gamma_{0}\mathcal{O}_{i}\gamma_{\mu}\mathcal{O}_{j}\gamma_{\nu}\mathcal{O}_{i}\gamma_{\lambda}\mathcal{O}_{j})\bigr{\\}},$
(40)

where $\mathcal{O}_{i}$ again denotes the $4\times 4$ matrix in the four-
fermion term parametrized by $g_{i}$, and for brevity we have omitted factors
of $\mathbb{1}_{2}$ in direct products with $\gamma$ matrices, i.e.,
$\gamma_{\mu}\equiv\gamma_{\mu}\otimes\mathbb{1}_{2}$. Note that in the above
equation, no summation over repeated indices $i$ and $j$ is assumed, and we
have rescaled the couplings as $g_{i}/2\pi\mapsto g_{i}$, which agrees with
the rescaling below Eqs. (13)–(18) for the present regularization scheme.

Evaluating the matrix algebra for the Gross-Neveu-Ising fixed point
$\boldsymbol{g}_{\text{GNI}}^{\star}$ yields

$\eta_{\psi}^{\text{GNI}}=\frac{4N_{\mathrm{f}}-1}{8(2N_{\mathrm{f}}-1)^{2}}\epsilon^{2}+\mathcal{O}(\epsilon^{3}),$
(41)

in agreement with the literature results [69, 70, 71, 19], providing another
cross-check of our approach.

For the Gross-Neveu-Heisenberg fixed point
$\boldsymbol{g}_{\text{GNH}}^{\star}=[h_{1}(N_{\mathrm{f}}),h_{2}(N_{\mathrm{f}}),0,0,0,h_{6}(N_{\mathrm{f}})]\epsilon$,
we find

$\displaystyle\eta_{\psi}^{\text{GNH}}$
$\displaystyle=\Bigl{[}3(4N_{\mathrm{f}}+1)h_{1}^{2}+24N_{\mathrm{f}}h_{2}^{2}+(4N_{\mathrm{f}}-1)h_{6}^{2}$
$\displaystyle\quad+12h_{1}h_{2}-6h_{1}h_{6}+12h_{2}h_{6}\Bigr{]}\frac{\epsilon^{2}}{2N_{\mathrm{f}}^{2}}+\mathcal{O}(\epsilon^{3})$
(42)

for general $N_{\mathrm{f}}$, leading to

$\eta_{\psi}^{\text{GNH}}=\left[\frac{3}{2N_{\mathrm{f}}}-\frac{9}{8N_{\mathrm{f}}^{2}}-\frac{3}{4N_{\mathrm{f}}^{3}}+\mathcal{O}(1/N_{\mathrm{f}}^{4})\right]\epsilon^{2}+\mathcal{O}(\epsilon^{3})$
(43)

in the large-$N_{\mathrm{f}}$ limit. The first two terms agree with the
previous large-$N_{\mathrm{f}}$ calculation in fixed space-time dimension
$2<D<4$ [39], when expanding the latter for small $\epsilon=D-2$. The third
term $\propto 1/N_{\mathrm{f}}^{3}$ does not agree: We have found
$-{3\epsilon^{2}}/(4N_{\mathrm{f}}^{3})$, whereas Ref. [39] suggests
$-9\epsilon^{2}/(8N_{\mathrm{f}}^{3})$. However, this discrepancy can be
traced back to a term $-2/[3(\mu-1)]$ on the right-hand side of Eq. (6.6) of
Ref. [39], which should not be there. Without that term, the
large-$N_{\mathrm{f}}$ result, when expanded near two dimensions, fully agrees
with our Eq. (43).444We are grateful to John Gracey for pointing this out to
us. For the physically relevant cases [24, 26, 25, 30, 40, 41, 42], we find

$\eta_{\psi}^{\text{GNH}}=\begin{cases}\frac{7}{72}\epsilon^{2}+\mathcal{O}(\epsilon^{3}),&\text{for
}N_{\mathrm{f}}=2,\\\
0.0748866\epsilon^{2}+\mathcal{O}(\epsilon^{3}),&\text{for
}N_{\mathrm{f}}=4,\\\
0.0422519\epsilon^{2}+\mathcal{O}(\epsilon^{3}),&\text{for
}N_{\mathrm{f}}=8,\end{cases}$ (44)

which represents another important result of our work.

###  IV.4 Corrections-to-scaling exponent $\boldsymbol{\omega}$

The exponent $\omega$ determines the leading corrections to scaling near the
quantum critical point. It is given by the negative of the second-largest
eigenvalue $\Theta_{2}<0$ of the stability matrix
$(-\partial\beta_{i}/\partial g_{j})$ at the corresponding critical fixed
point. The leading-order results for the Gross-Neveu-Ising and Gross-Neveu-
Heisenberg fixed points are shown in Figs. 3(a) and (b), respectively.

In the Gross-Neveu-Ising case, the corrections-to-scaling exponent $\omega$
vanishes for $N_{\mathrm{f}}\to N_{\mathrm{f}}^{(1)}=3/2$. This is a direct
consequence of the fixed-point collision occurring at this value of
$N_{\mathrm{f}}$. For $N_{\mathrm{f}}<N_{\mathrm{f}}^{(1)}$, the Gross-Neveu-
Ising fixed point develops a second relevant direction, as discussed in Sec.
III.3.

Remarkably, in the Gross-Neveu-Heisenberg case, the exponent features a
distinct minimum of $\omega\approx 0.3\epsilon+\mathcal{O}(\epsilon^{2})$ near
$N_{\mathrm{f}}=2$. This can be understood to arise from the competition
between the different interaction channels in the vicinity of the critical
fixed point. In particular, the infrared relevant direction of the Gross-
Neveu-Heisenberg fixed point, which is parallel to the fixed-point vector
$\boldsymbol{g}_{\text{GNH}}^{\star}$ itself, has a large component
perpendicular to the RG invariant plane spanned by the fixed points
$\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$. This RG invariant plane is
characterized by an O(4) symmetry generated by
$(\mathbbm{1}_{2},\gamma_{5})\otimes\vec{\sigma}$, under which the bilinears
$\bar{\psi}(\mathbbm{1}_{2}\otimes\vec{\sigma})\psi$ and
$\mathrm{i}\bar{\psi}(\gamma_{5}\otimes\mathbbm{1}_{2})\psi$ transform as
components of an O(4) vector, and which is an enhancement of the spin SU(2)
symmetry defined in Eq. (6). This is illustrated in the inset of Fig. 3(b),
which shows the angle $\varphi$ between $\boldsymbol{g}_{\text{GNH}}^{\star}$
and the surface normal $\boldsymbol{n}=(1,0,1)/\sqrt{2}$ of the RG invariant
plane spanned by the fixed points $\mathcal{A}$, $\mathcal{B}$, and
$\mathcal{C}$, featuring a distinct minimum near $N_{\mathrm{f}}=2$. The
presence of an RG invariant plane in the vicinity of the Gross-Neveu-
Heisenberg fixed point perpendicular to the fixed point’s relevant direction
arguably leads to a slow flow on the critical surface, i.e., towards the
critical point. Near $N_{\mathrm{f}}=2$, the corrections-to-scaling exponent
$\omega$ is therefore relatively small, implying that fluctuations over a
comparatively large number of length scales need to be integrated out to
approach the ultimate infrared behavior. We emphasize that this result arises
from the competition between the different interaction channels within our
Fierz-complete basis, and could not have been obtained within standard
$4-\varepsilon$ or large-$N_{\mathrm{f}}$ approaches, which typically involve
the fluctuations within the respective condensation channel only.

##  V Estimates for Gross-Neveu-Heisenberg criticality in 2+1 dimensions

The knowledge of the critical exponents in $D=2+\epsilon$ space-time
dimensions, together with literature results for these exponents in
$D=4-\varepsilon$ dimensions [13], allows us to employ an interpolational
resummation scheme in order to obtain estimates for the critical exponents in
the physical cases for $D=2+1$. In the Gross-Neveu-Ising case, such an
approach has previously been shown to lead to a significant improvement [16,
43]. Here, we focus on the leading exponents $1/\nu$, $\eta_{\phi}$, and
$\eta_{\psi}$. We do not attempt an interpolation of the corrections-to-
scaling exponent $\omega$, since the leading corrections to scaling, as
obtained in the previous subsection, arise from the competition between
different interaction channels, which has not been included in $4-\varepsilon$
expansion approaches to date.

We employ a scheme based on two-sided Padé approximants defined as

$[m/n](D):=\frac{a_{0}+a_{1}D+\dots+a_{m}D^{m}}{1+b_{1}D+\dots+b_{n}D^{n}},$
(45)

with non-negative integers $m$ and $n$, and real coefficients
$a_{0},\dots,a_{m}$ and $b_{1},\dots,b_{n}$, chosen such that the Padé
approximant matches both the $2+\epsilon$ and $4-\epsilon$ results, when
expanding the approximant near the lower and upper critical dimensions,
respectively. The order $m+n$ of the Padé approximant is determined by the
number of constraints given by the $2+\epsilon$ and $4-\varepsilon$ results.
Near the upper critical dimension, all leading exponents are known up to
quartic order in $\varepsilon=4-D$. Near the lower critical dimension, we have
computed the correlation-length exponent $1/\nu$ and the boson anomalous
dimension $\eta_{\phi}$ to linear order, and the fermion anomalous dimension
to quadratic order in $\epsilon=D-2$. This implies that the orders of the
corresponding Padé approximants are $m+n=6$ for $1/\nu$ and $\eta_{\phi}$, and
$m+n=7$ for $\eta_{\psi}$, respectively. While in principle several choices
for $m$ and $n$ are possible, some of these cannot satisfy all constraints
near $D=2$ and $D=4$ for real coefficients. This applies to $m=0$ for
$\eta_{\phi}$ and $1/\nu$, as well as to $m=0,1,2$ for $\eta_{\psi}$.
Furthermore, some choices lead to singularities of the corresponding Padé
approximants between $2<D<4$.

Table 2: Critical exponents of the Gross-Neveu-Heisenberg universality class for $N_{\mathrm{f}}=2$ four-component Dirac fermions in $D=3$ space-time dimensions, relevant for the transition between the Dirac semimetal and the antiferromagnetic insulator in the Hubbard model on the honeycomb lattice [26, 25, 30]. Here, we have used different two-sided Padé approximants $[m/n]$, interpolating between the expansions near the lower and upper critical dimensions. In the upper (lower) part of the table, marked as $\mathcal{O}(\epsilon,\varepsilon^{4})$ [$\mathcal{O}(\epsilon^{2},\varepsilon^{4})$], we have employed the results to linear (quadratic) order in $\epsilon=D-2$ and to quartic order in $\varepsilon=4-D$. The latter are obtained from Ref. [13]. Approximants that cannot satisfy all constraints are marked as “n.e.”, those exhibiting singularities in $2<D<4$ dimensions are marked as “sing.” The dashes “$-$” signify approximants for which the required $\epsilon^{2}$ corrections are not yet available. $N_{\mathrm{f}}=2$ | $[m/n]$ | $1/\nu$ | $\eta_{\phi}$ | $\eta_{\psi}$
---|---|---|---|---
$\mathcal{O}(\epsilon,\varepsilon^{4})$ | $[1/5]$ | 0.83569 | 1.07386 | n.e.
| $[2/4]$ | 0.76888 | 1.02731 | n.e.
| $[3/3]$ | sing. | 0.97641 | sing.
| $[4/2]$ | sing. | 0.95000 | 0.15943
| $[5/1]$ | 0.75902 | 1.02035 | sing.
| $[6/0]$ | 0.94485 | 1.03755 | 0.15592
$\mathcal{O}(\epsilon^{2},\varepsilon^{4})$ | $[3/4]$ | $-$ | $-$ | 0.14750
| $[4/3]$ | $-$ | $-$ | 0.10199
| $[5/2]$ | $-$ | $-$ | sing.
| $[6/1]$ | $-$ | $-$ | sing.
| $[7/0]$ | $-$ | $-$ | 0.13960

Figure 4: Correlation-length exponent $1/\nu$ (left column), order-parameter
anomalous dimension $\eta_{\phi}$ (center column), and fermion anomalous
dimension $\eta_{\psi}$ (right column) of the Gross-Neveu-Heisenberg
universality class as function of space-time dimension $2<D<4$ for
$N_{\mathrm{f}}=2$ (first row), $N_{\mathrm{f}}=4$ (second row), and
$N_{\mathrm{f}}=8$ (last row). The different curves in each panel correspond
to different Padé approximants $[m/n]$, which interpolate between the series
expansions around the lower and upper critical dimensions $D=2$ and $D=4$,
respectively. Data points at $D=3$ refer to literature results from
$4-\varepsilon$ expansion ($\triangledown$) [13], $1/N_{\mathrm{f}}$ expansion
($\triangle$) [39, 40, 42], functional RG ($\circ$) [16, 38, 40], as well as
determinantal quantum Monte Carlo ($\diamond$) [27, 28, 29, 30, 31] and hybrid
Monte Carlo ($\square$) [33, 34, 35, 36] simulations.

Figure 4 shows the non-singular Padé approximants for the critical exponents
at the Gross-Neveu-Heisenberg fixed point as a function of space-time
dimension $2<D<4$ for the physically relevant cases on the single-layer [24,
26, 25, 30] and bilayer [40, 41, 42] honeycomb lattices, i.e., for
$N_{\mathrm{f}}=2$, $N_{\mathrm{f}}=4$, and $N_{\mathrm{f}}=8$. The numerical
values in $D=2+1$ space-time dimensions of the different Padé approximants are
given in Tables 2, 3, and 4, respectively.

Table 3: Same as Table 2, but for $N_{\mathrm{f}}=4$ four-component Dirac fermions, relevant for the transition between nematic and coexistent nematic-antiferromagnetic orders on the Bernal-stacked honeycomb bilayer [40]. $N_{\mathrm{f}}=4$ | $[m/n]$ | $1/\nu$ | $\eta_{\phi}$ | $\eta_{\psi}$
---|---|---|---|---
$\mathcal{O}(\epsilon,\varepsilon^{4})$ | $[1/5]$ | 0.86441 | 1.03391 | sing.
| $[2/4]$ | 0.84006 | 1.00147 | sing.
| $[3/3]$ | sing. | sing. | 0.06418
| $[4/2]$ | sing. | sing. | 0.05950
| $[5/1]$ | 0.83956 | 0.99998 | sing.
| $[6/0]$ | 0.89489 | 1.02706 | 0.05906
$\mathcal{O}(\epsilon^{2},\varepsilon^{4})$ | $[3/4]$ | $-$ | $-$ | 0.05848
| $[4/3]$ | $-$ | $-$ | 0.05570
| $[5/2]$ | $-$ | $-$ | 0.05776
| $[6/1]$ | $-$ | $-$ | 0.05887
| $[7/0]$ | $-$ | $-$ | 0.05886
Table 4: Same as Table 2, but for $N_{\mathrm{f}}=8$ four-component Dirac fermions, relevant for the transition between the trigonal-warping-induced Dirac semimetal and the antiferromagnetic insulator on the Bernal-stacked honeycomb bilayer [41, 42]. $N_{\mathrm{f}}=8$ | $[m/n]$ | $1/\nu$ | $\eta_{\phi}$ | $\eta_{\psi}$
---|---|---|---|---
$\mathcal{O}(\epsilon,\varepsilon^{4})$ | $[1/5]$ | 0.92445 | 1.01040 | n.e.
| $[2/4]$ | 0.90752 | sing. | n.e.
| $[3/3]$ | 0.89871 | sing. | 0.02769
| $[4/2]$ | 0.90389 | sing. | 0.02679
| $[5/1]$ | 0.90865 | sing. | sing.
| $[6/0]$ | 0.93249 | 1.00979 | 0.02591
$\mathcal{O}(\epsilon^{2},\varepsilon^{4})$ | $[3/4]$ | $-$ | $-$ | 0.02682
| $[4/3]$ | $-$ | $-$ | 0.02657
| $[5/2]$ | $-$ | $-$ | 0.02660
| $[6/1]$ | $-$ | $-$ | 0.02718
| $[7/0]$ | $-$ | $-$ | 0.02676

The final best-guess estimates are determined via averaging over the results
from the different non-singular Padé approximants for the highest-order
expansion results available for each exponent. We thus arrive at

$1/\nu=\begin{cases}0.83(12),&\text{for }N_{\mathrm{f}}=2,\\\
0.86(4),&\text{for }N_{\mathrm{f}}=4,\\\ 0.913(20),&\text{for
}N_{\mathrm{f}}=8,\end{cases}$ (46)

for the correlation-length exponent, as well as

$\eta_{\phi}=\begin{cases}1.014(64),&\text{for }N_{\mathrm{f}}=2,\\\
1.016(18),&\text{for }N_{\mathrm{f}}=4,\\\ 1.0101(3),&\text{for
}N_{\mathrm{f}}=8,\end{cases}$ (47)

and

$\eta_{\psi}=\begin{cases}0.130(28),&\text{for }N_{\mathrm{f}}=2,\\\
0.0579(22),&\text{for }N_{\mathrm{f}}=4,\\\ 0.0268(4),&\text{for
}N_{\mathrm{f}}=8,\end{cases}$ (48)

for the boson and fermion anomalous dimensions, respectively. In the above
equations, the numbers in parentheses correspond to the maximal deviations
from the mean values among the different approximants, which can be understood
as a lower bound for the uncertainty of our best-guess estimates. We emphasize
that the true systematic error is hard to quantify and may be significantly
larger than this lower bound. This is particularly true for cases in which
only few non-singular Padé approximants exist. Nevertheless, we find it
reassuring that for the case of $N_{\mathrm{f}}=8$, for which the
large-$N_{\mathrm{f}}$ expansion is expected to yield reliable results, our
estimates for $1/\nu$ and $\eta_{\phi}$ are within error bars fully consistent
with the large-$N_{\mathrm{f}}$ results quoted in Ref. [42], and our estimate
for $\eta_{\psi}$ is within error bars almost consistent with those of Ref.
[42].555Note that the definition for $N$ used in Ref. [42] deviates from our
definition for $N_{\mathrm{f}}$ as
$2N^{\text{(Ref.~{}\cite[cite]{[\@@bibref{Number}{ray18}{}{}]})}}=N_{\mathrm{f}}^{\text{(this
work)}}$.

Our estimates are compared with a variety of literature results available for
$N_{\mathrm{f}}=2$ from $4-\varepsilon$ expansion [13], $1/N_{\mathrm{f}}$
expansion [39], functional RG [16, 38], as well as determinantal quantum Monte
Carlo [27, 28, 29, 30, 31, 32] and hybrid Monte Carlo [33, 34, 35, 36]
simulations in Table 5. Available literature results for different
$N_{\mathrm{f}}$ are included as black dots in Fig. 4. For $N_{\mathrm{f}}=2$,
the deviations between the results of the different methods are considerable:
For $1/\nu$, the analytical estimates are typically significantly smaller than
those of determinantal quantum Monte Carlo calculations; the hybrid Monte
Carlo estimates lie roughly between these two. For $\eta_{\phi}$, on the other
hand, the analytical estimates are significantly larger than those of most of
the quantum Monte Carlo simulations. In the case of $\eta_{\psi}$, analytical
estimates are again significantly smaller than those of determinantal quantum
Monte Carlo calculations; hybrid Monte Carlo estimates for $\eta_{\psi}$ are
not available at present. For $N_{\mathrm{f}}=4$, literature results from
$4-\varepsilon$ expansion [13], $1/N_{\mathrm{f}}$ expansion [39], and
functional RG [40], all of which as compiled in Ref. [40], agree very well
with our results in the case of $1/\nu$ and $\eta_{\phi}$; some deviations, in
particular from the functional RG estimate, are present in the case of
$\eta_{\psi}$.

Table 5: Gross-Neveu-Heisenberg critical exponents for $N_{\mathrm{f}}=2$ four-component Dirac fermions from interpolation between series expansions near lower and upper critical dimensions (this work) in comparison with previous results from fourth-order $4-\varepsilon$ expansion [13], second-order (third-order for $\eta_{\psi}$) $1/N_{\mathrm{f}}$ expansion [39], functional RG in local potential approximation (LPA’) [16] and next-to-leading order derivative expansion (NLO) [38], as well as determinantal quantum Monte Carlo (DQMC) [27, 28, 29, 30, 31, 32] and hybrid Monte Carlo (HMC) [33, 34, 35, 36] simulations for different linear lattice sizes $L$ and inverse temperatures $\beta$ or projection times $\tau$. In cases where $\eta_{\phi}$ and/or $1/\nu$ were not computed directly, we have employed appropriate hyperscaling relations to obtain these. In cases where results from different Padé approximants ($4-\varepsilon$ expansion), different regulators (functional RG), or different lattices (Monte Carlo simulations) are available within the same work, we show the corresponding mean values. $N_{\mathrm{f}}=2$ | Year | $1/\nu$ | $\eta_{\phi}$ | $\eta_{\psi}$
---|---|---|---|---
Interpolation (this work) | 2022 | 0.83(12) | 1.01(6) | 0.13(3)
$4-\varepsilon$ expansion, $\mathcal{O}(\varepsilon^{4})$ [13] | 2017 | 0.64 | 0.98 | 0.17
$1/N_{\mathrm{f}}$ expansion, $\mathcal{O}(1/N_{\mathrm{f}}^{2,3})$ [39] | 2018 | 0.85 | 1.18 | 0.11
functional RG, NLO [38] | 2018 | 0.80 | 1.03 | 0.07
functional RG, LPA’ [16] | 2014 | 0.77 | 1.01 | 0.08
DQMC, $\tau\sim L\leq 40$ [31] | 2020 | 0.95(5) | 0.75(4) | 0.23(4)
DQMC, $\tau\sim L\leq 40$ [30] | 2016 | 0.98(1) | 0.47(7) | 0.22(2)
DQMC, $\beta=L\leq 24$ [29] | 2021 | 1.11(4) | 0.80(9) | 0.29(2)
DQMC, $\beta=L\leq 21$ [28] | 2019 | 1.14(9) | 0.79(5) | –
DQMC, $\tau=60$, $L\leq 18$ [27] | 2015 | 1.19(6) | 0.70(15) | –
DQMC, $\beta=L\leq 20$ [32] | 2021 | 1.01(8) | 0.55(2) | –
HMC, $\beta\leq 12$, $L\leq 102$ [36] | 2021 | 0.84(4) | 0.52(1) | –
HMC, $\beta\leq 12$, $L\leq 102$ [35] | 2020 | 0.84(4) | 0.85(13) | –
HMC, $\beta=21,L\leq 24$ [34] | 2019 | 1.08 | 0.62 | –
HMC, $\beta=21,L\leq 18$ [33] | 2018 | 0.86 | 0.87(2) | –

##  VI Conclusions

To conclude, we have determined the critical behavior of the Gross-Neveu-
Heisenberg universality class within an $\epsilon$ expansion around the lower
critical space-time dimension of two. In contrast to the Gross-Neveu-Ising
case [19], the critical fixed point associated with the Gross-Neveu-Heisenberg
universality class is characterized by a combination of different four-fermion
interaction channels, requiring an approach that takes these channels into
account in an unbiased way. For the Gross-Neveu-Heisenberg case, a Fierz-
complete basis of the theory space compatible with the symmetries of the model
comprises six four-fermion interaction terms. Applying the general formula
derived in Ref. [22] to this system has allowed us to derive the flow
equations of this six-dimensional theory space. By making use of hyperscaling
relations and the flow of an infinitesimal symmetry breaking fermion bilinear,
we have demonstrated how to compute the full set of critical exponents within
the fermionic language. Applying this scheme to the Gross-Neveu-Ising fixed
point, for which various literature results are available, facilitates a
nontrivial cross-check of our approach. Our results for the leading-order
order-parameter anomalous dimension $\eta_{\phi}$ and the next-to-leading
order fermion anomalous dimension $\eta_{\psi}$ at the Gross-Neveu-Heisenberg
fixed point are original.

These results have allowed us to obtain improved estimates for the critical
exponents in $D=2+1$ space-time dimensions, as relevant for interacting
fermion models on the honeycomb and bilayer honeycomb lattices. Here, we have
employed a resummation scheme that takes the expansions near the lower and
upper critical dimensions simultaneously into account. For the Gross-Neveu-
Ising case, such an interpolational approach has previously been shown to
provide significantly more reliable estimates in comparison with standard
extrapolation schemes [16, 43]. In the Gross-Neveu-Heisenberg case, our
results for $N_{\mathrm{f}}=8$, relevant for the transition between the
trigonal-warping-induced semimetal and the antiferromagnetic insulator on the
Bernal-stacked honeycomb bilayer [41, 42], agree with previous
large-$N_{\mathrm{f}}$ estimates [39] within an uncertainty on the level of
$\lesssim 3\%$. For $N_{\mathrm{f}}=4$, relevant for the nematic-to-
coexistence transition on the honeycomb bilayer [40], the deviations between
our estimates and the large-$N_{\mathrm{f}}$ results [39], upon appropriate
resummation of the latter [40], are only slightly larger as compared with the
$N_{\mathrm{f}}=8$ case, with the largest relative difference of $\simeq 8\%$
occurring for the fermion anomalous dimension.

Interestingly, for $N_{\mathrm{f}}=2$, relevant for the semimetal-to-
antiferromagnet transition in the honeycomb-lattice Hubbard model [26, 25,
30], we have found that the Gross-Neveu-Heisenberg fixed point is
characterized by a slow flow towards criticality, corresponding to a small
corresponding exponent $\omega$ and generically sizable corrections to
scaling. This result can be understood to arise from the competition between
different interaction channels present at the Gross-Neveu-Heisenberg fixed
point. This is in contrast to the Gross-Neveu-Ising fixed point, which is
characterized by a single and uniquely identifiable interaction channel. The
critical point of a lattice model of, e.g., spinless fermions interacting via
a repulsive nearest-neighbor density-density interaction, can therefore be
close in theory space to the Gross-Neveu-Ising fixed point, leading to small
scaling corrections. The generically large scaling corrections in the Gross-
Neveu-Heisenberg case for $N_{\mathrm{f}}=2$ might explain the significant
spread between the estimates from the various numerical and analytical
approaches, cf. Table 5. To track down the origin of these discrepancies, it
would be interesting to test whether the data obtained in the simulations are
in principle compatible with a small corrections-to-scaling exponent $\omega$.
Within our one-loop analysis, we estimate $\omega\approx 0.3$ for
$N_{\mathrm{f}}=2$; however, a more accurate estimate, obtained from, e.g., a
full two-loop analysis around the lower critical dimension, or an
interpolation between the lower and upper critical dimensions, would certainly
be highly desirable. An interpolational approach to estimate $\omega$ would
require to compute the scaling dimensions of the different four-fermion terms
within the $4-\varepsilon$ expansion, which might be an interesting direction
for future work.

On more general grounds, our work demonstrates how to determine the critical
behavior of fermion models in cases where the corresponding critical fixed
point is characterized by different four-fermion interaction channels. This
should be of relevance for other fermionic universality classes as well. In
particular, our general formulas for the order-parameter anomalous dimension
$\eta_{\phi}$ to linear order in $\epsilon=D-2$, see Eqs. (32) and (33), and
the fermion anomalous dimension $\eta_{\psi}$ to quadratic order in
$\epsilon$, see Eqs. (39) and (40), together with the general formula for the
flow equations of relativistic four-fermion models [22], could be immediately
applied to other relativistic universality classes, such as Gross-Neveu-XY
[72, 73, 74], Gross-Neveu-SO($N$) [75, 76, 77], or nematic [78, 79]
transitions. These may host even more interesting phenomena, such as emergent
supersymmetry [80, 81, 82], fixed-point annihilation and complexification [22,
77], or quasiuniversal behavior [79].

###### Acknowledgements.

We thank John Gracey for very valuable discussions and comments on the
manuscript, and Michael Scherer for collaborations on related topics. This
work has been supported by the Deutsche Forschungsgemeinschaft (DFG) through
SFB 1143 (A04 and A07, Project No. 247310070), the Würzburg-Dresden Cluster of
Excellence ct.qmat (EXC 2147, Project No. 390858490), and the Emmy Noether
program (ME4844/1-1, Project No. 327807255, and JA2306/4-1, Project No.
411750675).

## References

  * Gross and Neveu [1974] D. J. Gross and A. Neveu, Dynamical symmetry breaking in asymptotically free field theories, Phys. Rev. D 10, 3235 (1974).
  * Gawȩdzki and Kupiainen [1985] K. Gawȩdzki and A. Kupiainen, Renormalizing the nonrenormalizable, Phys. Rev. Lett. 55, 363 (1985).
  * Rosenstein  _et al._ [1989] B. Rosenstein, B. J. Warr, and S. H. Park, Four-fermion theory is renormalizable in 2+1 dimensions, Phys. Rev. Lett. 62, 1433 (1989).
  * Zinn-Justin [1991] J. Zinn-Justin, Four-fermion interaction near four dimensions, Nucl. Phys. B 367, 105 (1991).
  * Braun  _et al._ [2011] J. Braun, H. Gies, and D. D. Scherer, Asymptotic safety: A simple example, Phys. Rev. D 83, 085012 (2011).
  * Hands  _et al._ [1993] S. Hands, A. Kocic, and J. Kogut, Four-Fermi Theories in Fewer Than Four Dimensions, Ann. Phys. (N. Y.) 224, 29 (1993).
  * Wang  _et al._ [2014] L. Wang, P. Corboz, and M. Troyer, Fermionic quantum critical point of spinless fermions on a honeycomb lattice, New J. Phys. 16, 103008 (2014).
  * Li  _et al._ [2015] Z.-X. Li, Y.-F. Jiang, and H. Yao, Fermion-sign-free Majarana-quantum-Monte-Carlo studies of quantum critical phenomena of Dirac fermions in two dimensions, New J. Phys. 17, 085003 (2015).
  * Hesselmann and Wessel [2016] S. Hesselmann and S. Wessel, Thermal Ising transitions in the vicinity of two-dimensional quantum critical points, Phys. Rev. B 93, 155157 (2016).
  * Huffman and Chandrasekharan [2017] E. Huffman and S. Chandrasekharan, Fermion bag approach to Hamiltonian lattice field theories in continuous time, Phys. Rev. D 96, 114502 (2017).
  * Huffman and Chandrasekharan [2020] E. Huffman and S. Chandrasekharan, Fermion-bag inspired Hamiltonian lattice field theory for fermionic quantum criticality, Phys. Rev. D 101, 074501 (2020).
  * Rosenstein  _et al._ [1993] B. Rosenstein, H.-L. Yu, and A. Kovner, Critical exponents of new universality classes, Phys. Lett. B 314, 381 (1993).
  * Zerf  _et al._ [2017] 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. D 96, 096010 (2017).
  * Rosa  _et al._ [2001] L. Rosa, P. Vitale, and C. Wetterich, Critical Exponents of the Gross-Neveu Model from the Effective Average Action, Phys. Rev. Lett. 86, 958 (2001).
  * Höfling  _et al._ [2002] F. Höfling, C. Nowak, and C. Wetterich, Phase transition and critical behavior of the $d=3$ Gross-Neveu model, Phys. Rev. B 66, 205111 (2002).
  * Janssen and Herbut [2014] 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).
  * Vacca and Zambelli [2015] G. P. Vacca and L. Zambelli, Multimeson Yukawa interactions at criticality, Phys. Rev. D 91, 125003 (2015).
  * Knorr [2016] B. Knorr, Ising and Gross-Neveu model in next-to-leading order, Phys. Rev. B 94, 245102 (2016).
  * Gracey  _et al._ [2016] J. A. Gracey, T. Luthe, and Y. Schröder, Four loop renormalization of the Gross-Neveu model, Phys. Rev. D 94, 125028 (2016).
  * Gracey [2018a] J. A. Gracey, Large $N_{f}$ quantum field theory, Int. J. Mod. Phys. A 33, 1830032 (2018a).
  * Iliesiu  _et al._ [2018] L. Iliesiu, F. Kos, D. Poland, S. S. Pufu, and D. Simmons-Duffin, Bootstrapping 3D fermions with global symmetries, J. High Energy Phys. 1 (2018) 36.
  * Gehring  _et al._ [2015] F. Gehring, H. Gies, and L. Janssen, Fixed-point structure of low-dimensional relativistic fermion field theories: Universality classes and emergent symmetry, Phys. Rev. D 92, 085046 (2015).
  * Sorella and Tosatti [1992] S. Sorella and E. Tosatti, Semi-Metal-Insulator Transition of the Hubbard Model in the Honeycomb Lattice, Europhys. Lett. 19, 699 (1992).
  * Herbut [2006] I. F. Herbut, Interactions and Phase Transitions on Graphene’s Honeycomb Lattice, Phys. Rev. Lett. 97, 146401 (2006).
  * Assaad and Herbut [2013] 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).
  * Herbut  _et al._ [2009] I. F. Herbut, V. Juričić, and B. Roy, Theory of interacting electrons on the honeycomb lattice, Phys. Rev. B 79, 085116 (2009).
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb and $\pi$-flux Hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum Monte Carlo, Phys. Rev. B 91, 165108 (2015).
  * Liu  _et al._ [2019] Y. Liu, Z. Wang, T. Sato, M. Hohenadler, C. Wang, W. Guo, and F. F. Assaad, Superconductivity from the condensation of topological defects in a quantum spin-Hall insulator, Nat. Commun. 10, 2658 (2019).
  * Liu  _et al._ [2021] Y. Liu, Z. Wang, T. Sato, W. Guo, and F. F. Assaad, Gross-Neveu Heisenberg criticality: Dynamical generation of quantum spin Hall masses, Phys. Rev. B 104, 035107 (2021).
  * Otsuka  _et al._ [2016] 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).
  * Otsuka  _et al._ [2020] 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).
  * Xu and Grover [2021] X. Y. Xu and T. Grover, Competing Nodal $d$-Wave Superconductivity and Antiferromagnetism, Phys. Rev. Lett. 126, 217002 (2021).
  * Buividovich  _et al._ [2018] 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).
  * Buividovich  _et al._ [2019] P. Buividovich, D. Smith, M. Ulybyshev, and L. von Smekal, Numerical evidence of conformal phase transition in graphene with long-range interactions, Phys. Rev. B 99, 205434 (2019).
  * Ostmeyer  _et al._ [2020] 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).
  * Ostmeyer  _et al._ [2021] 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).
  * Lang and Läuchli [2019] 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).
  * Knorr [2018] B. Knorr, Critical chiral Heisenberg model with the functional renormalization group, Phys. Rev. B 97, 075129 (2018).
  * Gracey [2018b] J. A. Gracey, Large $N$ critical exponents for the chiral Heisenberg Gross-Neveu universality class, Phys. Rev. D 97, 105009 (2018b).
  * Ray and Janssen [2021] S. Ray and L. Janssen, Gross-Neveu-Heisenberg criticality from competing nematic and antiferromagnetic orders in bilayer graphene, Phys. Rev. B 104, 045101 (2021).
  * Pujari  _et al._ [2016] S. Pujari, T. C. Lang, G. Murthy, and R. K. Kaul, Interaction-Induced Dirac Fermions from Quadratic Band Touching in Bilayer Graphene, Phys. Rev. Lett. 117, 086404 (2016).
  * Ray  _et al._ [2018] S. Ray, M. Vojta, and L. Janssen, Quantum critical behavior of two-dimensional Fermi systems with quadratic band touching, Phys. Rev. B 98, 245128 (2018).
  * Ihrig  _et al._ [2018] B. Ihrig, L. N. Mihaila, and M. M. Scherer, Critical behavior of Dirac fermions from perturbative renormalization, Phys. Rev. B 98, 125109 (2018).
  * Jaeckel and Wetterich [2003] J. Jaeckel and C. Wetterich, Flow equations without mean field ambiguity, Phys. Rev. D 68, 025020 (2003).
  * Bennett and Gracey [1999] J. F. Bennett and J. A. Gracey, Three-loop renormalization of the $SU(N_{c})$ non-abelian Thirring model, Nucl. Phys. B 563, 390 (1999).
  * Ali and Gracey [2001] D. B. Ali and J. A. Gracey, Four loop wave function renormalization in the non-abelian Thirring model, Nucl. Phys. B 605, 337 (2001).
  * Bondi  _et al._ [1990] A. Bondi, G. Curci, G. Paffuti, and P. Rossi, Metric and central charge in the perturbative approach to two dimensional fermionic models, Ann. Phys. (N. Y.) 199, 268 (1990).
  * Gies and Janssen [2010] H. Gies and L. Janssen, UV fixed-point structure of the three-dimensional Thirring model, Phys. Rev. D 82, 085018 (2010).
  * Janssen and Gies [2012] L. Janssen and H. Gies, Critical behavior of the ($2+1$)-dimensional Thirring model, Phys. Rev. D 86, 105007 (2012).
  * Halperin  _et al._ [1974] 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).
  * Nahum  _et al._ [2015] A. Nahum, J. T. Chalker, P. Serna, M. Ortuño, and A. M. Somoza, Deconfined Quantum Criticality, Scaling Violations, and Classical Loop Models, Phys. Rev. X 5, 041048 (2015).
  * Ihrig  _et al._ [2019] 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).
  * Gies and Jaeckel [2006] H. Gies and J. Jaeckel, Chiral phase structure of QCD with many flavors, Eur. Phys. J. C 46, 433 (2006).
  * Kaplan  _et al._ [2009] D. B. Kaplan, J.-W. Lee, D. T. Son, and M. A. Stephanov, Conformality lost, Phys. Rev. D 80, 125005 (2009).
  * Braun  _et al._ [2014] J. Braun, H. Gies, L. Janssen, and D. Roscher, Phase structure of many-flavor ${\mathrm{QED}}_{3}$, Phys. Rev. D 90, 036002 (2014).
  * Janssen [2016] L. Janssen, Spontaneous breaking of Lorentz symmetry in ($2+\epsilon$)-dimensional QED, Phys. Rev. D 94, 094013 (2016).
  * Herbut [2016] I. F. Herbut, Chiral symmetry breaking in three-dimensional quantum electrodynamics as fixed point annihilation, Phys. Rev. D 94, 025036 (2016).
  * Gukov [2017] S. Gukov, RG flows and bifurcations, Nucl. Phys. B 919, 583 (2017).
  * Herbut and Janssen [2014] I. F. Herbut and L. Janssen, Topological Mott Insulator in Three-Dimensional Systems with Quadratic Band Touching, Phys. Rev. Lett. 113, 106401 (2014).
  * Janssen and Herbut [2017] L. Janssen and I. F. Herbut, Phase diagram of electronic systems with quadratic Fermi nodes in $2<d<4$: $2+\epsilon$ expansion, $4-\epsilon$ expansion, and functional renormalization group, Phys. Rev. B 95, 075101 (2017).
  * Herbut and Janssen [2016] I. F. Herbut and L. Janssen, Critical $O(2)$ and $O(3)$ ${\varphi}^{4}$ theories near six dimensions, Phys. Rev. D 93, 085005 (2016).
  * Gracey  _et al._ [2018] J. A. Gracey, I. F. Herbut, and D. Roscher, Tensor $O(N)$ model near six dimensions: Fixed points and conformal windows from four loops, Phys. Rev. D 98, 096014 (2018).
  * Gorbenko  _et al._ [2018] V. Gorbenko, S. Rychkov, and B. Zan, Walking, Weak first-order transitions, and Complex CFTs II. Two-dimensional Potts model at $Q>4$, SciPost Phys. 5, 50 (2018).
  * Ma and He [2019] H. Ma and Y.-C. He, Shadow of complex fixed point: Approximate conformality of $Q>4$ Potts model, Phys. Rev. B 99, 195130 (2019).
  * Ma and Wang [2020] R. Ma and C. Wang, Theory of deconfined pseudocriticality, Phys. Rev. B 102, 020407 (2020).
  * Nahum [2020] A. Nahum, Note on Wess-Zumino-Witten models and quasiuniversality in $2+1$ dimensions, Phys. Rev. B 102, 201116 (2020).
  * [67] M. Weber and M. Vojta, SU(2)-symmetric spin-boson model: Quantum criticality, fixed-point annihilation, and duality, arXiv:2203.02518 .
  * [68] H. Hu and Q. Si, Kondo destruction and fixed-point annihilation in a Bose-Fermi Kondo model, arXiv:2207.08744 .
  * Gracey [1990] J. Gracey, Three-loop calculations in the O($N$) Gross-Neveu model, Nucl. Phys. B 341, 403 (1990).
  * Gracey [1991] J. Gracey, Computation of the three-loop $\beta$-function of the O($N$) Gross-Neveu model in minimal subtraction, Nucl. Phys. B 367, 657 (1991).
  * Gracey [2008] J. Gracey, Four loop $\overline{\text{MS}}$ mass anomalous dimension in the Gross-Neveu model, Nucl. Phys. B 802, 330 (2008).
  * Bobev  _et al._ [2015] N. Bobev, S. El-Showk, D. Mazáč, and M. F. Paulos, Bootstrapping the Three Dimensional Supersymmetric Ising Model, Phys. Rev. Lett. 115, 051601 (2015).
  * Li  _et al._ [2017] Z.-X. Li, Y.-F. Jiang, S.-K. Jian, and H. Yao, Fermion-induced quantum critical points, Nat. Commun. 8, 314 (2017).
  * Classen  _et al._ [2017] L. Classen, I. F. Herbut, and M. M. Scherer, Fluctuation-induced continuous transition and quantum criticality in Dirac semimetals, Phys. Rev. B 96, 115132 (2017).
  * Seifert  _et al._ [2020] U. F. P. Seifert, X.-Y. Dong, S. Chulliparambil, M. Vojta, H.-H. Tu, and L. Janssen, Fractionalized Fermionic Quantum Criticality in Spin-Orbital Mott Insulators, Phys. Rev. Lett. 125, 257202 (2020).
  * Ray  _et al._ [2021] S. Ray, B. Ihrig, D. Kruti, J. A. Gracey, M. M. Scherer, and L. Janssen, Fractionalized quantum criticality in spin-orbital liquids from field theory beyond the leading order, Phys. Rev. B 103, 155160 (2021).
  * Janssen and Seifert [2022] L. Janssen and U. F. P. Seifert, Phase diagrams of SO($N$) Majorana-Hubbard models: Dimerization, internal symmetry breaking, and fluctuation-induced first-order transitions, Phys. Rev. B 105, 045120 (2022).
  * Vojta  _et al._ [2000] M. Vojta, Y. Zhang, and S. Sachdev, Quantum Phase Transitions in $\mathit{d}$-Wave Superconductors, Phys. Rev. Lett. 85, 4940 (2000).
  * Schwab  _et al._ [2022] J. Schwab, L. Janssen, K. Sun, Z. Y. Meng, I. F. Herbut, M. Vojta, and F. F. Assaad, Nematic Quantum Criticality in Dirac Systems, Phys. Rev. Lett. 128, 157203 (2022).
  * Lee [2007] S.-S. Lee, Emergence of supersymmetry at a critical point of a lattice model, Phys. Rev. B 76, 075103 (2007).
  * Jian  _et al._ [2015] S.-K. Jian, Y.-F. Jiang, and H. Yao, Emergent Spacetime Supersymmetry in 3D Weyl Semimetals and 2D Dirac Semimetals, Phys. Rev. Lett. 114, 237001 (2015).
  * Gies  _et al._ [2017] H. Gies, T. Hellwig, A. Wipf, and O. Zanusso, A functional perspective on emergent supersymmetry, J. High Energy Phys. 12 (2017) 132.