Torus Spectroscopy of the Gross-Neveu-Yukawa Quantum Field Theory: Free Dirac versus Chiral Ising Fixed Point

TL;DR

This paper characterizes the low-energy spectra of the Gross-Neveu-Yukawa quantum field theory at different fixed points using numerical and analytical methods, revealing insights into quantum critical behavior of Dirac systems.

Contribution

It provides the first comprehensive comparison of numerical and analytical torus spectra at the free Dirac and chiral Ising fixed points in 2+1 dimensions.

Findings

Numerical spectra match analytical predictions, confirming the chiral Ising fixed point.

Identified spectral features indicating emergent symmetry at criticality.

Analyzed crossover behavior between fixed points and corrected previous Fermi velocity estimates.

Abstract

We establish the universal torus low-energy spectra at the free Dirac fixed point and at the strongly coupled chiral Ising fixed point and their subtle crossover behaviour in the Gross-Neuveu-Yukawa field theory with ${n_\text{D}=4}$ component Dirac spinors in $D=(2+1)$ dimensions. These fixed points and the field theories are directly relevant for the long-wavelength physics of certain interacting Dirac systems, such as repulsive spinless fermions on the honeycomb lattice or $\pi$-flux square lattice. The torus energy spectrum has been shown previously to serve as a characteristic fingerprint of relativistic fixed points and is a powerful tool to discriminate quantum critical behaviour in numerical simulations. Here, we use a combination of exact diagonalization and quantum Monte Carlo simulations of strongly interacting fermionic lattice models, to compute the critical torus energy spectrum on finite-size clusters with periodic boundaries and extrapolate them to the thermodynamic limit. Additionally, we compute the torus energy spectrum analytically using the perturbative expansion in ${\epsilon = 4 - D}$, which is in good agreement with the numerical results, thereby validating the presence of the chiral Ising fixed point in the lattice models at hand. We show that the strong interaction between the spinor field and the scalar order-parameter field strongly influences the critical torus energy spectrum and we observe prominent multiplicity features related to an emergent symmetry predicted from the quantum field theory. Building on these results we are able to address the subtle crossover physics of the low-energy spectrum flowing from the chiral Ising fixed point to the Dirac fixed point, and analyze earlier flawed attempts to extract Fermi velocity renormalizations from the low-energy spectrum.

Tables

Table 1. Table 1: Quantum numbers of the most relevant energy levels in the torus spectrum of the chiral Ising CFT. The table denotes the particle-hole quantum number PH, the fermion sector relative to half filling n f subscript 𝑛 𝑓 n_{f} , the momentum 𝐤 𝐤 \mathbf{k} , and the reduced momentum κ 𝜅 \kappa . For simplicity, we here omit showing the irreducible representations under the lattice point-group symmetry.

Levels	PH	$I_x$	$n_f$	$𝐤$	$\kappa$
$1_T,{ϵ}_T$	$1$	$1$	$0$	$0$	$0$
${\sigma }_T,{\sigma }_T^{\prime }$	$-1$	$-1$	$0$	$0$	$0$
Kekule	$1$	$1$	$0$	$𝐊$	$0$
Chern	$-1$	$-1$	$0$	$0$	$0$
${\psi }_T$	–	–	$\pm 1$	$𝐊$	$0$
${\left(\psi +d\right)}_T,{\left(\psi +d\right)}_T^{\prime }$	–	–	$\pm 1$	${𝐤}_{\min }$	$1$
${\left(2\psi \right)}_T$	–	–	$\pm 2$	$0$	$0$

Table 2. Table 2: Low-lying spectrum of the N = 4 𝑁 4 N=4 GNY model on the torus, from numerically computing the strong-coupling expansion in the text. We measure energies with respect to the ground state. The states are labeled by the eigenvalue Q 𝑄 Q defined in the text, and their degeneracy (deg.) is given. The labelling of the levels is used to agree with that used in previous sections, where C and K refer to Chern and Kekule respectively.

Level	$Q$	deg.	$\sqrt{𝒜}E$
$1_T$, ${\sigma }_T$	$\pm 2$	2	0
${\psi }_T$	$\pm 1$	8	$2.6033{ϵ}^{1/3}+0.3287\sqrt{{\tau }_2}{f}_{1/2}^{\left(3\right)}\left(\tau \right){ϵ}^{2/3}$
C, K, $\left(2{\psi }_T\right)$	$0$	6	$3.7443{ϵ}^{1/3}+0.5698\sqrt{{\tau }_2}{f}_{1/2}^{\left(3\right)}\left(\tau \right){ϵ}^{2/3}$
${ϵ}_T$, ${\sigma }_T^{\prime }$	$\pm 2$	2	$4.4713{ϵ}^{1/3}+0.1183\sqrt{{\tau }_2}{f}_{1/2}^{\left(3\right)}\left(\tau \right){ϵ}^{2/3}$

Table 3. Table 3: Low-lying spectrum of the N = 4 𝑁 4 N=4 GNY model for particular torus shapes considered in numerics, written in terms of the ξ CI superscript 𝜉 CI \xi^{\text{CI}} defined in Section III.2 . We measure energies with respect to the ground state, extrapolated to ϵ = 1 italic-ϵ 1 \epsilon=1 . Here we give the energy spectrum for modular parameters τ = exp ⁡ ( i ​ π / 3 ) 𝜏 𝑖 𝜋 3 \tau=\exp(i\pi/3) and τ = i 𝜏 𝑖 \tau=i appropriate to the honeycomb and square lattices respectively.

Level	$Q$	deg.	$\frac{3}{4\pi }{\xi }_h^{\text{CI}}$	$\frac{3}{4\pi }{\xi }_s^{\text{CI}}$
$1_T$, ${\sigma }_T$	$\pm 2$	2	0	0
${\psi }_T$	$\pm 1$	8	0.246	0.247
C, K, $\left(2{\psi }_T\right)$	$0$	6	0.312	0.315
${ϵ}_T$, ${\sigma }_T^{\prime }$	$\pm 2$	2	0.612	0.613

Table 4. Table 4: Symmetry properties of the operators connecting the ground state to the relevant excited states for the honeycomb lattice. The tables denotes the action under particle-hole transformation PH, the fermion sector n f subscript 𝑛 𝑓 n_{f} and the momentum 𝐤 𝐤 \mathbf{k} , as well as the action under rotations about 60 ∘ superscript 60 60^{\circ} ( C 6 subscript 𝐶 6 C_{6} ) and 120 ∘ superscript 120 120^{\circ} ( C 3 subscript 𝐶 3 C_{3} ), lattice inversion ( C 2 subscript 𝐶 2 C_{2} ), and vertical ( I x subscript 𝐼 𝑥 I_{x} ) and horizontal ( I y subscript 𝐼 𝑦 I_{y} ) mirror reflection. 1 ​ ( − 1 ) 1 1 1\ (-1) indicates a vanishing commutator (anti-commutator), and – corresponds to an otherwise broken symmetry. In the case of degenerate excited states operators with different symmetries can be allowed, which connect to superimposed states of the degenerate subspace.

Level	PH	$n_f$	$𝐤$	$C_6$	$C_3$	$C_2$	$I_x$	$I_y$
$1_T,{ϵ}_T$	$1$	$0$	$0$	$1$	$1$	$1$	$1$	$1$
${\sigma }_T,{\sigma }_T^{\prime }$	$-1$	$0$	$0$	$1$	$-1$	$-1$	$-1$	$1$
Kekule	$1$	$0$	$𝐊$	$1$	–	–	$1$	–
Chern	$-1$	$0$	$0$	–	–	$-1$	$-1$	$1$
${\psi }_T$	–	$1$	$𝐊$	–	–	–	–	–
${\left(2\psi \right)}_T$	–	$2$	$0$	–	–	–	–	–

Equations

$${\mathcal{L}}_{\text{GNY}}=-{\overline{\Psi }}^j\left(\cancel{\partial }+g_{\mathrm{Y}}\varphi \right){\Psi }^j+\frac{1}{2}\varphi \left(s-{\partial }^2\right)\varphi +\frac{\lambda }{4!}{\varphi }^4,$$

$$\Psi \to {\gamma }^1\Psi ,\overline{\Psi }\to -\overline{\Psi }{\gamma }^1,\varphi \to -\varphi .$$

$${\mathcal{L}}_{\text{D}}=-{\overline{\Psi }}^j\cancel{\partial }{\Psi }^j.$$

$${\mathcal{L}}_{\text{GN}}=-{\overline{\Psi }}^j\cancel{\partial }{\Psi }^j-\frac{g}{2}{\left({\overline{\Psi }}^j{\Psi }^j\right)}^2,$$

$$H_{\mathrm{h}}=-t\sum _{⟨i,j⟩}\left(c_i^{†}{c}_j^{}+\text{h.c.}\right)+V\sum _{⟨i,j⟩}\left(n_i-\frac{1}{2}\right)\left(n_j-\frac{1}{2}\right).$$

$$H_{\mathrm{s}}=-t\sum _{⟨i,j⟩}\left({\mathrm{e}}^{\mathrm{i}{\theta }_{ij}}{c}_i^{†}{c}_j^{}+\text{h.c.}\right)+V\sum _{⟨i,j⟩}\left(n_i-\frac{1}{2}\right)\left(n_j-\frac{1}{2}\right).$$

$$T:c_{\alpha }(x_i,y_i)\to c_{\alpha }(x_i,y_i).$$

$$I_x:\left(\begin{array}{c}{c}_A(x_i,y_i)\\ c_B(x_i,y_i)\end{array}\right)\to \left(\begin{array}{c}{c}_B(-x_i,y_i)\\ c_A(-x_i,y_i)\end{array}\right),$$

$$I_y:\left(\begin{array}{c}{c}_A(x_i,y_i)\\ c_B(x_i,y_i)\end{array}\right)\to \left(\begin{array}{c}{c}_A(x_i,-y_i)\\ c_B(x_i,-y_i)\end{array}\right).$$

$$\mathcal{C}:\left(\begin{array}{c}{c}_A(x_i,y_i)\\ c_B(x_i,y_i)\end{array}\right)\to \left(\begin{array}{c}{c}_A^{†}(x_i,y_i)\\ -c_B^{†}(x_i,y_i)\end{array}\right),$$

$${\gamma }^0=\left(\begin{array}{cc}{\sigma }^z& 0\\ 0& {\sigma }^z\end{array}\right),{\gamma }^1=\left(\begin{array}{cc}{\sigma }^x& 0\\ 0& {\sigma }^x\end{array}\right),{\gamma }^2=\left(\begin{array}{cc}{\sigma }^y& 0\\ 0& -{\sigma }^y\end{array}\right).$$

$${\gamma }^3=\left(\begin{array}{cc}0& {\sigma }^y\\ {\sigma }^y& 0\end{array}\right),{\gamma }^4=\left(\begin{array}{cc}0& -\mathrm{i}{\sigma }^y\\ \mathrm{i}{\sigma }^y& 0\end{array}\right).$$

$$T:$$

$$I_x:$$

$$I_y:$$

$$\mathcal{C}:$$

$${\Delta }_i=\frac{v_{\text{F}}(V)}{L}{\xi }_i^{\text{D}},$$

$$\kappa =\frac{3}{4\pi }L|\mathbf{k}\mathrm{mod}{\mathbf{K}}^{(\prime )}|,$$

$${\Delta }_i=\frac{v_{\text{F}}^c}{L}{\xi }_i^{\text{CI}},$$

$${\Delta }_i$$

$$\delta \ll L\ll |V-V_c{|}^{-\nu }.$$

$$v_{\text{F}}(V)=\underset{L\to \infty }{\lim }\frac{{\Delta }_{i}L}{{\xi }_i^{\text{D}}},(V<V_c).$$

$$v_{\text{F}}^{\text{slope}}(V)=\underset{L\to \infty }{\lim }\frac{\Delta ({\mathbf{k}}_{\min })-\Delta (\mathbf{K})}{|{\mathbf{k}}_{\min }-\mathbf{K}|},$$

$$\frac{⟨{\Psi }_0|\hat{O}|{\Psi }_0⟩}{⟨{\Psi }_0|{\Psi }_0⟩}=\underset{\Theta \to \infty }{\lim }\frac{⟨{\Psi }_T|{\mathrm{e}}^{-\frac{\Theta }{2}H}\hat{O}{\mathrm{e}}^{-\frac{\Theta }{2}H}|{\Psi }_T⟩}{⟨{\Psi }_T|{\mathrm{e}}^{-\Theta H}|{\Psi }_T⟩},$$

$$⟨\hat{O}(\tau ){\hat{O}}^{†}⟩\sim {\left|⟨{\Psi }_0|\hat{O}|{\Psi }_{{\Delta }_{\hat{O}}}⟩\right|}^2\exp (-{\Delta }_{\hat{O}}\tau ),$$

$${\Delta }_i\to {\hat{\Delta }}_i=\frac{v_{\text{F}}^0\left|\mathbf{D}-{\mathbf{k}}_{\text{min}}\right|}{{\epsilon }^0\left({\mathbf{k}}_{\text{min}}\right)}{\Delta }_i,$$

$$H_{\text{GNY}}$$

$$\left[\varphi (x),\Pi (x^{\prime })\right]$$

$$\left\{{\Psi }_a^j(x),{\Psi }_b^{j†}(x^{\prime })\right\}$$

$$g_{\mathrm{Y},ren.}^{\ast 2}$$

$${\lambda }_{\text{ren}}^{\ast }$$

$$\varphi (x)$$

$$\Pi (x)$$

Full text

I Introduction

Dirac fermions with a quasi-relativistic, i.e., gapless and linear, dispersion relation arise as low-energy quasi-particles in many condensed matter systems such as graphene and $d$-wave superconductors [1, 2, 3]. Interactions between the fermions can drive the system from the Dirac semimetallic (SM) phase through a quantum critical point (QCP) into various symmetry broken phases. Spinless fermions on the honeycomb or $\pi$-flux square lattice with repulsive nearest-neighbour interactions, for example, exhibit a SM to Mott insulator transition, where the ground state is charge ordered and spontaneously breaks discrete sublattice exchange symmetries [2, 4]. The quantum critical point of such phase transitions involves fermionic degrees of freedom and is believed to be described by the chiral Ising fixed point of the $D=(2+1)$ dimensional Gross-Neveu-Yukawa (GNY) field theory which features strong coupling between fermionic spinors and a scalar field [5, 6, 3, 7, 8]. Such universality classes do not have classical Landau-Ginzburg-Wilson analogues and are thus of particular interest. Recently, many attempts have been made to precisely measure the scaling dimensions of the operators of chiral QCPs in GNY field theories [4, 7, 9, 10, 11, 12, 13, 14, 15], which are a unique identifier of the universality class and directly related to the critical exponents of the phase transition. This task turned out to be particularly challenging and different methods, which were successful for charting more common critical points in the past, could not yet obtain completely consistent scaling dimensions or critical exponents. This summarizes the current situation for both the chiral Ising, and even more severely, the chiral Heisenberg universality class [10, 16].

Here, we strike a new path to tackle this problem. In fact, another way to identify and chart universality classes is to measure their critical torus energy spectrum as it was shown in Refs. [17, 18, 19, 20] for Wilson-Fisher and topological phase transitions. The low-energy gaps at a relativistic critical point on the torus are given, up to a non-universal factor $v$ describing the effective speed of light, by universal numbers $\xi_{i}$ times $1/L$, where $L$ is the linear extent of the cluster 111More generally, for a QCP with dynamic critical exponent $z$, the finite size spectrum takes the form $E_{i}=v\xi_{i}/L^{z}$, where $v$ is a non-universal pre-factor and the dimensionless numbers $\xi_{i}$ are universal constants. In this paper, we specialize to the relativistic case $z=1$, where $v$ is equal to the velocity of low-energy excitations. The order and degeneracy of the $\xi_{i}$ together with quantum numbers of the corresponding eigenstate (e.g. momentum, fermion number) provide a unique fingerprint of the universality class and can be obtained by many complementary numerical and analytical techniques.

Given the reported tension in the literature we want to shed light on the nature of the universality class from a different angle by confronting and comparing numerical torus energy spectra with analytical results. In this work, we use exact diagonalization (ED) and quantum Monte Carlo (QMC) simulations of fermionic lattice models, as well as $\epsilon$ expansion calculations of the effective low-energy field theory, to compute the critical torus energy spectrum for the chiral Ising universality class with an $(n_{\text{D}}=4)$-component spinor field. Although the $\epsilon$ expansion is only performed to low order, it provides important, exact statements about non-trivial multiplicities and quantum numbers of the low-energy spectrum in the scaling limit, while we provide high-quality data of the $\xi_{i}$ from numerical simulations.

Also, we want to emphasize that, while the $\xi_{i}$ are universal numbers, it is not only their precise values but typically the sequence of the low-energy levels with their degeneracies and quantum numbers which make the critical torus spectrum a universal fingerprint. In particular, changing the nature of the QCP by using another $n_{D}\neq 4$ would lead to a different multiplicity structure, i.e. a qualitative change. This is one of the main advantages of the critical torus spectrum compared to measuring critical exponents, where often very high precision data is necessary to distinguish different universality classes.

Furthermore, the study of the GNY field theory and the corresponding microscopic models allows us to clarify the crossover flow of the torus energy spectrum between two different infrared (IR) fixed points. This advance enables us to reliably measure the Fermi velocity, which is the condensed matter analogy of the speed of light of the GNY field theory, a topic of recent controversy [22, 23, 24].

Our work also completes an important intermediate step towards a quantitative understanding of massless Dirac fermions coupled to a $U(1)$ gauge field (QED3), which are of paramount importance for many quantum spin liquid candidates and exotic quantum phase transitions [25, 26, 27, 28, 29].

The paper is organized as follows. In Sec. II we introduce the GNY field theory, as well as the related Gross-Neveu (GN) field theory, which is formulated purely in terms of interacting fermions and introduce the two distinct renormalization group fixed points under consideration in this work. We present the fermion lattice models used to compute the chiral Ising critical torus spectrum, and establish the GNY field theory as a low-energy effective description of the lattice models. In Sec. III we provide a brief overview of our results: We discuss the different structures of the torus spectra of the free Dirac conformal field theory (CFT) and chiral Ising CFT, the crossover behaviour in finite volume and their impact on the renormalization of the Fermi velocity. In Sec. IV we give a more detailed analysis of the critical torus energy spectrum obtained from numerics. We show energy gaps from both ED and QMC simulations, and give details on the extrapolation to the thermodynamic limit. In Sec. V we present the $\epsilon$ expansion of the GNY field theory, and compare the results to the numerical spectra. Finally, in Sec. VI we conclude our results by comparing the different torus geometries among each other and discuss possible future perspectives.

II Field Theories and Model Hamiltonians

This section provides a concise introduction of the GNY field theory, the important infrared fixed points along with symmetry aspects that are relevant for the analysis of the torus spectrum. We also introduce the microscopic quantum many-body lattice models that exhibit fermionic quantum critical points, and which we examine by our numerical methods.

II.1 Quantum Field Theories

The fermionic quantum field theories (QFTs) that we explore in this work can be described by the GNY theory of fermionic fields coupled to a $\mathbb{Z}_{2}$ order parameter, i.e., a real, one-component scalar field [2, 30, 8] in $D=3$ (space-time) dimensions. Depending on the value of a tuning parameter $s$, the GNY theory describes a SM of non-interacting Dirac fermions, a symmetry broken phase with finite order parameter, and, in between those, a critical point belonging to the chiral Ising universality class [5, 6]. The most general form of the imaginary-time GNY Lagrangian is

[TABLE]

where $\Psi^{j}$ is an $n_{\text{D}}$-component Dirac spinor with ${j=1,...,{N_{f}}}$ flavors, so the total number of fermionic degrees of freedom is ${N=n_{\text{D}}N_{f}}$. The real scalar field is denoted by $\phi$, and $g_{\rm Y}$ is the Yukawa coupling strength between the spinor and scalar fields. We use the standard notation ${\cancel{\partial}=\gamma^{\mu}\partial_{\mu}}$, ${\mu\in\{0,\dots D-1\}}$ and ${\overline{\Psi}=\Psi\gamma^{0}}$, where the $\gamma^{\mu}$ are $n_{\text{D}}\times n_{\text{D}}$ matrices satisfying the Clifford algebra, $\{\gamma^{\mu},\gamma^{\nu}\}=2\delta^{\mu\nu}$. In these expressions, we have set the speed of light to unity. In $D=3$, a Dirac spinor has a minimum of $n_{\text{D}}=2$ components; however, in applications to condensed matter systems, the number of two-component Dirac fermions in a bulk lattice system is always doubled due to fermion doubling arguments [1, 31], so the total number of fermionic degrees of freedom $N$ is always a multiple of four.

In $D=3$, there is a critical value of the tuning parameter, $s=s_{c}$, such that for ${s<s_{c}}$ the scalar order parameter acquires a finite expectation value, ${\langle\phi\rangle\neq 0}$. Such a finite expectation value spontaneously breaks the ($\mathbb{Z}_{2}$) parity symmetry of the theory, which is given by taking $(x^{0},x^{1},x^{2})\rightarrow(x^{0},-x^{1},x^{2})$ together with

[TABLE]

The finite expectation value of $\phi$ acts as a Dirac mass in Eq. (1), resulting in a massive spectrum of fermions above a two-fold degenerate ground state. In contrast, for ${s>s_{c}}$, the parameter $s$ flows to positive infinity while $g_{\mathrm{Y}}$ and $\lambda$ flow to zero. In this limit, we may ignore the gapped bosonic fields, and at long distances the theory describes a SM of non-interacting, massless Dirac fermions with the Euclidean Lagrangian

[TABLE]

We call this fixed point the Dirac CFT, and its properties are easily obtained since $\mathcal{L}_{\text{D}}$ is exactly solvable.

Directly at the QCP, $s=s_{c}$, the interaction couplings $g_{\rm Y}$ and $\lambda$ flow to non-zero values of an interacting fixed point, determined by the chiral Ising universality class. We hence denote the critical theory of this emerging interacting fixed point the chiral Ising CFT. This QCP is non-perturbative directly in $D=3$, but there exists a perturbative expansion in $\epsilon=4-D$, where $\lambda\sim g_{\rm Y}^{2}\sim O(\epsilon)$, and the universal properties of the QCP may be obtained after extrapolating to $\epsilon=1$. This will be our primary analytic tool for studying the finite-size torus spectrum as detailed in Sec. V.

One may alternatively describe the above QCP using a purely fermionic field theory, the Gross-Neveu (GN) model [32], whose imaginary-time Lagrangian is

[TABLE]

with a self-interaction of strength $g>0$. For $D=3$, the coupling $g$ is renormalization group (RG) irrelevant in perturbation theory, so a weak-coupling analysis always results in a stable massless Dirac SM phase with $\langle\overline{\Psi}\Psi\rangle=0$. However, there is ample evidence for a non-perturbative UV fixed point at some value $g=g_{c}$, where for $g>g_{c}$ the system flows to strong coupling. At strong coupling, the system dynamically generates a mass by acquiring an expectation value $\langle\overline{\Psi}\Psi\rangle\neq 0$, spontaneously breaking the ($\mathbb{Z}_{2}$) parity transformation of Eq. (2), which will be examined in more detail in Sec. II.3. By the principle of universality, this fixed point should also be in the chiral Ising universality class, although it is only analytically accessible in an expansion in $\epsilon^{\prime}=D-2$ [5].

We note that both the GNY and GN models may be studied directly in $D=3$ (and even in fractional dimensions $2<D<4$) by a perturbative expansion in $1/N$, where it may be shown that the fixed points of the two models are exactly equivalent in the scaling limit within this expansion [5]. We discuss the $1/N$ and $\epsilon^{\prime}$ expansions of the torus spectrum in Appendix D, where we additionally give checks that the three expansions all give consistent torus spectra to leading order.

In both the GNY and GN field theories, there is a global U($N_{f}$) symmetry obtained by taking $\Psi^{j}\rightarrow U^{jj^{\prime}}\Psi^{j^{\prime}}$, with $U\in\mathrm{U}(N_{f})$. When $n_{\text{D}}=2$, it is sometimes conventional to decompose each Dirac fermion into two Majorana fields, after which the theory is invariant under the larger group O($2N_{f}$). This Majorana formulation is conventionally used to define the O($N$)-invariant chiral Ising CFTs [7, 11]. For $n_{\text{D}}>2$, these field theories no longer have an explicit O($N$) symmetry; however, due to the structure of the perturbative expansion for the beta functions, the scaling dimensions of all operators turn out to only depend on $N$ rather than $n_{\text{D}}$ or $N_{f}$ separately, which leads to identical critical properties for all theories with the same total number of components $N=n_{\text{D}}N_{f}$. We will show that the torus spectrum also only depends on $N$. Therefore, we conjecture that all of the chiral Ising CFTs with the same number of total degrees of freedom $N$ flow to the same O($N$)-invariant CFTs irrespective of the smaller global symmetries present in the Lagrangians of Eqns. (1) and (4), and will also obey an identical critical torus energy spectrum.

This emergent O($N$) symmetry also results in a particular degeneracy structure of the low-energy eigenstates of the critical torus energy spectrum, as shown in Secs. III.2 and V.1.

II.2 Model Hamiltonians

Our perturbative analysis of the GNY torus spectrum will be complemented by the non-perturbative analysis of microscopic fermionic quantum lattice models that exhibit a QCP which is widely believed to belong to the chiral Ising universality class as described by the GNY field theory Eq. (1). In particular, we consider two models of spinless fermions with Dirac cones in the non-interacting limit at half filling.

The first model is defined on a honeycomb lattice [see Fig. 1(a)] with the Hamiltonian

[TABLE]

Here, $\langle i,j\rangle$ stands for nearest neighbour bonds on the honeycomb lattice, $c_{i}^{(\dagger)}$ denotes fermionic annihilation (creation) operators at site $i$, and $n_{i}=c_{i}^{\dagger}c^{\phantom{\dagger}}_{i}$ is the fermionic number operator. The first term in Eq. (5) is the tight-binding model description of the fermionic hopping between nearest neighbour sites. The free dispersion relation for positive frequencies is shown in Fig. 2(a) and has Dirac cones at the two non-equivalent Dirac points in the Brillouin zone (BZ), denoted $\mathbf{K}$ and $\mathbf{K}^{\prime}$, characteristic for a SM state. At half-filling the spectrum is particle-hole symmetric. The second term describes a density-density interaction between fermions on neighboring sites, driving the system, for $V\gg t$, into a charge density wave (CDW) state, in which the particle-hole symmetry, together with a sublattice exchange parity symmetry, (a $\mathbb{Z}_{2}$ symmetry group) is spontaneously broken and the local densities within the two sublattices differ in the thermodynamic limit. The position of the QCP was determined previously to be at $V_{c}/t\approx 1.355$ [33, 4, 34, 35].

The second model that we consider describes interacting, spinless fermions on the $\pi$-flux square lattice at half filling [see Fig. 1(b)], with the Hamiltonian

[TABLE]

While the phases $\theta_{ij}$ depend on the choice of gauge, the total flux per plaquette, $\phi=\sum_{\langle i,j\rangle\in\Box}\theta_{ij}$ is fixed to be $\phi=\pi$. In the non-interacting case, $V=0$, Eq. (6) is a tight-binding model with, in general, complex hopping amplitudes. Its dispersion relation (for the specific choice of gauge $\theta_{ij}=\pi/4$) is plotted in Fig. 2(b) and again shows two distinct Dirac cones (now at the $\mathbf{X}$ and $\mathbf{X}^{\prime}$ points) in the BZ. For large values of $V\gg t$, this model also exhibits a CDW phase with different densities on the two sublattices that are coupled by the repulsion term. The position of the QCP was estimated to be at $V_{c}/t\approx 1.279$ [33, 4, 36]. In the remainder of this paper, we will set the energy scale of the lattice models by fixing the hopping amplitude $t=1$.

The QPTs in both lattice models introduced above may be described by either the GNY or GN quantum field theories with a single ($N_{f}=1$) flavor of $n_{\text{D}}=4$-component spinors [2, 3], corresponding to a value of $N=4$, as reviewed in Appendix C. If not specified otherwise, we mean for $N$ to take on this value in the following. An important aspect of the correspondence between the microscopic lattice models and the effective QFT description is the manifestation of several global symmetries, which we will examine in the following section.

II.3 Symmetries

The model Hamiltonians $H_{\rm h}$ and $H_{\rm s}$ possess several global symmetries, some of which are spontaneously broken in the ordered phase. Here we review these relevant symmetries and discuss how they manifest in the GNY and GN field theories.

Since both model Hamiltonians are bipartite, we may label the fermion annihilation (creation) operators as ${c^{(\dagger)}_{\alpha}(x_{i},y_{i})}$, where ${\mathbf{r}_{i}=(x_{i},y_{i})}$ is the coordinate of the Bravais lattice, and $\alpha=A,B$ is the sublattice index. In this section, we take the convention from Fig. 1 that the $A$ and $B$ sites are connected along the $x$-axis, and that the unit cell is centered on the point directly equidistant between these two sites. Furthermore, with regards to $H_{\rm s}$, we consider here, for simplicity, a gauge, in which the phases $\mathrm{e}^{\mathrm{i}\theta_{ij}}$ in $H_{\rm s}$ are all real (for example, by choosing $\theta_{ij}=\pi$ on one link of each square plaquette and zero on the others); then these symmetries take on an identical form in both models. With these conventions, we can now describe the symmetries of both models.

First, we have the usual global U(1) symmetry from particle-number conservation, $c_{\alpha}(x_{i},y_{i})\rightarrow\mathrm{e}^{\mathrm{i}\varphi}c_{\alpha}(x_{i},y_{i})$. In addition, both models have an anti-unitary time-reversal symmetry that is given by complex conjugation in real space, leaving the fermionic operators unchanged,

[TABLE]

We may also define a parity flip across either the vertical or horizontal axes, $I_{x}:(x_{i},y_{i})\rightarrow(-x_{i},y_{i})$ and $I_{y}:(x_{i},y_{i})\rightarrow(x_{i},-y_{i})$. These parity symmetries are actually part of the larger point group symmetries of these models, e.g., the dihedral symmetry group $D_{6}$ on the honeycomb lattice. Besides a change in coordinates, our convention of parity implies that $I_{x}$ also exchanges the two sublattices, so we have

[TABLE]

[TABLE]

Finally, we also define a particle-hole transformation,

[TABLE]

under which the density on each site transforms as ${n_{i}\rightarrow(1-n_{i})}$. In the CDW phase, the densities of fermions on the sublattices $A$ and $B$ differ, such that both $I_{x}$ and $\mathcal{C}$ are spontaneously broken.

We now characterize the symmetries of the field theory. As shown in Appendix C, the GNY and GN field theories may be derived from our model Hamiltonians in the continuum and scaling limits, and we find that the models possess a single flavor of four-component Dirac fermions, i.e., $n_{\text{D}}=4$, $N_{f}=1$. In discussing this realization of the field theory, it is useful to introduce the following explicit representation of the gamma matrices, which arises naturally from the derivation:

[TABLE]

In this representation, the first two indices of the four-spinors represent fermions on the honeycomb (square) lattice with momentum near $\mathbf{K}$ ($\mathbf{X}$) at sublattice A and B respectively, while the bottom two components represent fermions near momentum $\mathbf{K}^{\prime}$ ($\mathbf{X}^{\prime}$) at sublattice A and B. It is useful to define two additional gamma matrices, which anticommute with the above,

[TABLE]

We also define ten Hermitian matrices $\gamma^{ab}\equiv\mathrm{i}\gamma^{a}\gamma^{b}$, where $0\leq a<b\leq 4$. This parametrization is useful because an arbitrary $4\times 4$ Hermitian matrix may be written as a linear combination of the sixteen matrices $\{\mathbb{I},\gamma^{a},\gamma^{ab}\}$ with real coefficients.

Given this explicit representation, we may follow the derivation of Appendix C to obtain how the microscopic symmetry transformations act on the fields of the QFT. In particular, we use the fact that the transformations $T$, $I_{y}$, and $\mathcal{C}$ exchange the Dirac points, while $I_{x}$ does not (see also Fig. 2). The U(1) symmetry simply takes the form $\Psi\rightarrow\mathrm{e}^{\mathrm{i}\varphi}\Psi$, while the discrete symmetries are given by

[TABLE]

where $T$ is anti-unitary. It is straightforward to show that Eqns. (1) and (4) are invariant under these transformations (where we simply ignore the transformation rules on $\phi$ in the GN case). We furthermore see that the order parameters for the $\mathbb{Z}_{2}$ symmetry breaking may be given by either $\langle\phi\rangle$ or $\langle\overline{\Psi}\Psi\rangle$, which are both odd under $I_{x}$ and $\mathcal{C}$ but even under the rest of the above symmetry transformations.

In addition to the symmetries inherited from the microscopic model, the field theory has additional symmetries which are not present in the lattice model. Importantly, it is invariant under the three-dimensional Lorentz group, which is generated by $\{\gamma^{01},\gamma^{02},\gamma^{12}\}$ 222In application to the imaginary-time Euclidean theories, these are SO(3) rotations, but they will correspond to SO(2,1) Lorentz transformations after a Wick rotation to real-time., and includes the discrete rotational symmetries of the lattice as a subgroup. Finally, we have an U(1) (“chiral”) symmetry given by $\Psi\rightarrow\exp(\mathrm{i}\gamma^{34})\Psi$, which corresponds to performing independent U(1) rotations at the two inequivalent Dirac points. We stress that these emergent symmetries of the field theory only apply to the lattice models in the strict scaling limit, and any degeneracies found in the torus spectrum of the field theory due to these particular symmetries are expected to be approximate in the lattice models due to the presence of additional irrelevant operators.

II.4 Torus Compactifications

In numerical studies of two-dimensional quantum lattice models, one typically considers finite-size clusters constructed from the underlying lattice, with periodic boundary conditions taken in both lattice directions. This way, one effectively studies a torus compactification of the original infinite lattice model. Here, we choose finite-size clusters which preserve the maximal six- (four-)fold rotational symmetry $C_{6}$ ($C_{4}$) of the honeycomb (square) lattice model. As mentioned in the introduction, our analysis serves the dual purpose of (i) examining the energy level structure of fermionic model systems on such torus geometries, which serves as a universal fingerprint for the corresponding QFT, as well as (ii) deriving appropriate estimators for other physical quantities, such as the effective Fermi velocity (the effective speed of light) of the interacting fermion models, which depend on such spectroscopic data.

The torus clusters for the two microscopic models that we consider here exhibit different overall rhombic shapes [see Fig. 1(c), (d)], which act as an infrared (IR) cutoff (irrespectively of the lattice discretization, i.e., the ultraviolet cutoff). The IR characteristics remain influential in the thermodynamic limit [17], and we account for them also in the analysis of the GNY field theory on finite tori. This is done as follows: In the continuum limit, one can use complex coordinates on the two-dimensional torus, $x=x_{1}+\mathrm{i}\,x_{2}$. The torus is then defined by two complex periods, $\omega_{1}$ and $\omega_{2}$, such that the points $x+n\omega_{1}+m\omega_{2}$ are equivalent for all $n,m\in\mathbb{Z}$. The torus shape is then characterized by its modular parameter, $\tau\equiv\omega_{2}/\omega_{1}=\tau_{1}+\mathrm{i}\tau_{2}$. In particular, the considered triangular-lattice based tori (such as for the honeycomb lattice model) have a value of $\tau=\exp(\mathrm{i}\pi/3)$, while the square-lattice based tori correspond to $\tau=\mathrm{i}$, respectively, c.f. Fig. 1 (c), (d). In the framework of the GNY field theory the torus corresponds to periodic boundary conditions for both the $\Psi$ and the $\phi$ fields.

The lattice models we consider all use periodic boundary conditions for their simulation clusters. Since the Dirac points in the considered models are not located at the $\Gamma$ point ($\mathbf{k}=(0,0)$), several families of clusters arise, which differ in their momentum discretization grid around the Dirac points, as illustrated in Fig. 3. For example, when considering finite torus clusters with $N_{s}=2L^{2}$ lattice sites of the honeycomb lattice (c.f. Fig. 1), those with linear size $L\!\!\mod 3=0$ feature the Dirac points ($\mathbf{K}$ and $\mathbf{K}^{\prime}$) in their momentum space (family H), in contrast to clusters with $L\!\!\mod 3\neq 0$, which do not feature the Dirac points (family H’), so that the spectrum is gapped already in the tight-binding limit, $V=0$. Hence, the finite-size torus spectrum is qualitatively different for those two families already in the non-interacting limit, and we observe characteristic differences also for the interacting case. The case $L\!\!\mod 3=0$ in the lattice models corresponds to the standard periodic boundary condition case in the GNY field theory discussed above. The second case $L\!\!\mod 3\neq 0$ corresponds to the GNY field theory with twisted boundary conditions, $\Psi(x+n\omega_{1}+m\omega_{2})=e^{i\theta}\Psi(x)$. We will present numerical results for those spectra later on, although we will only give a few comments about the structure of the $\epsilon$-expansion for twisted boundary conditions in Section V. Similar considerations apply to the fermionic model $H_{\rm s}$ on the $\pi$-flux square lattice, as detailed in Sec. IV.5.

III Overview of the Central Results

This section provides an overview of our analytical and numerical findings. Further details are provided in the subsequent sections of this paper. For all the models and the field theories that we consider, the phase diagram is divided by a strongly coupled QCP (chiral Ising CFT) into an extended SM regime of massless Dirac fermions (Dirac CFT) and a regime with spontaneous symmetry breaking and gapped fermionic excitations. While the spectroscopic properties of the strongly coupled chiral Ising fixed point are of particular interest, the excitation spectrum of the free Dirac CFT characterizes the SM regime, which thus exhibits distinctly different spectral characteristics. Additionally, there are important crossover effects between these two regimes which need to be treated with care. We therefore start by exploring the spectroscopic properties of the Dirac CFT, before discussing our main findings for the chiral Ising CFT. We then discuss the crossover between these two CFTs as well as the subtleties in obtaining the correct Fermi velocity renormalization. In this section, we concentrate on the case of the honeycomb lattice model $H_{\rm h}$ and torus clusters that contain the Dirac points (family H).

III.1 Dirac CFT Torus Spectrum

In the SM phase with $V<V_{c}$ the torus spectrum is characterized by the free, massless Dirac CFT, defined by Eq. (3). Its excitation energies can be readily calculated analytically, and they are directly related to the Fermi velocity of the Dirac fermions. For large finite clusters with $N_{s}=2L^{2}$ lattice sites the energy levels scale as $1/L$, and the torus spectrum in the SM regime is then given by

[TABLE]

where $\Delta_{i}$, $i=1,2,...$ denotes the set of finite-size energy gaps (relative to the ground state energy) which make up the low-energy torus spectrum. Here, $v_{\text{F}}(V)$ is the renormalized Fermi velocity at interaction strength $V$, and the $\xi^{\text{D}}_{i}$ is a set of universal numbers that characterize the Dirac CFT [17, 18, 19]. The values $\xi_{i}^{\text{D}}$ in the low-energy regime are shown in Fig. 4(a). Additional quantum numbers are attached to the corresponding eigenstate, in particular the fermion number $n_{f}$ relative to half filling, and its momentum $\mathbf{k}$. At half-filling, many-body states of vanishing finite-size gaps reside both at the Dirac points $\mathbf{K}$ and $\mathbf{K}^{\prime}$, as well as at total momentum $\mathbf{k}=0$ 333Note that in the field theory the soft energy levels are around $\mathbf{k}=0$, as the information about $\mathbf{K}$ and $\mathbf{K}^{\prime}$ has been absorbed into the spinor components.. For example, at $n_{f}=0$ we have to put two fermions into the upper or lower band of the Dirac cones at $\mathbf{K}$ and $\mathbf{K}^{\prime}$ to find in total six zero energy states. This can be done by either putting both fermions into the two-bands of a single Dirac cone which results in a total momentum of $\mathbf{k}=2\mathbf{K}^{(\prime)}\equiv-\mathbf{K}^{(\prime)}$ of the many-body state, or by putting each one in a different Dirac cone with a total momentum of $\mathbf{k}=\mathbf{K}+\mathbf{K}^{\prime}\equiv 0$. We thus introduce a reduced momentum variable

[TABLE]

taken modulo the momentum space lattice vectors spanned by $\mathbf{K}$ and $\mathbf{K}^{\prime}$. States with small value of $\kappa$ then map to the low-energy sector described by the QFT. The normalization factor $3/(4\pi)L$ is chosen such that momenta closest to $\mathbf{k}=0$ and to the Dirac points correspond to a value of $\kappa=1$ (for lattice constant $a=1$). The levels $\Delta_{i}$ together with their quantum numbers and multiplicities provide a characteristic fingerprint of the Dirac CFT. We show the low-energy part of the torus spectrum of the Dirac CFT as a function of the reduced momentum in Fig. 4(a), which also displays the corresponding multiplicities of each level.

In the following, we will refer to gaps in the torus spectrum as denoting finite differences between the rescaled energy gaps ${L\Delta_{i}\propto\xi_{i}}$. For example, in Fig. 4(a), the lowest level at $\kappa=1$ has a finite gap with respect to the lowest level at $\kappa=0$. The raw many body spectrum in the thermodynamic limit, however, shows no gaps for ${V\leq V_{c}}$. This notion of gaps in the torus spectrum will turn out particularly useful to quantify the differences between the torus spectrum of the Dirac CFT and the one at the chiral Ising critical point, which we discuss next.

III.2 Chiral Ising CFT Torus Spectrum

At the critical point $s=s_{c}$, or $V=V_{c}$, our models are described by the strongly interacting chiral Ising CFT, for which no exact analytical solutions are known. An analytical approach to the critical torus spectrum of the chiral Ising fixed point of the $N=4$ component GNY field theory is provided by the $\epsilon$ expansion as detailed in Sec. V. From these calculations and numerical simulations of the microscopic models, we find that the critical torus spectrum of the interacting chiral Ising CFT [i.e., at finite $g_{\rm Y}>0$ in Eq. (1)] is characterized by a (different) set of finite-size energy gaps $\Delta_{i}$ that scale as

[TABLE]

with $v^{c}_{\text{F}}=v_{\text{F}}(V_{c})$ the renormalized Fermi velocity at the critical interaction strength. Such a scaling form of the critical torus spectrum of an interacting fixed point has been obtained also in studies of purely bosonic quantum critical points [17, 18, 19]. It can be considered a mass spectrum of the quantum critical theory with a mass scale set by the IR cutoff, which is proportional to $v^{c}_{\text{F}}/L$. Here, the $\xi^{\text{CI}}_{i}$ are again a set of universal numbers which are, however, distinct from those of the Dirac CFT, and, together with the levels multiplicities and quantum numbers, identify the chiral Ising CFT. The $\epsilon$-expansion predicts a rich level multiplicity structure for the low-energy levels because of the before mentioned emergent O($N$) symmetry of the chiral Ising field theory.

In Fig. 4(b), we show the low-energy part of the chiral Ising torus spectrum for the $N=4$-component GNY theory, normalized by the Fermi velocity $v_{\text{F}}^{c}$ as a function of the reduced momentum $\kappa$. In this figure, we compare our estimates for the $\xi^{\text{CI}}_{i}$ as obtained from the $\epsilon$ expansion at $\kappa=0$ (left panel) with finite-size extrapolated gap data from the microscopic lattice model $H_{\rm h}$ (right panel). The structure of level multiplicities and the quantum numbers of the corresponding eigenstates as obtained from the perturbative $\epsilon$-expansion compare remarkably well to the numerical analysis of the microscopic lattice models (when comparing actual values, one needs to keep in mind, that in the $\epsilon$ expansion they result from a simple extrapolation of only the leading term). It is important to note, that the degenerate levels obtained from the perturbative analysis appear as quasi-degenerate levels in the numerical data, because of corrections to scaling present in any particular lattice implementation of the GNY field theory. A very prominent feature of the chiral Ising torus spectrum is the opening of large gaps in the torus spectrum between the two (quasi-) degenerate ground state levels with $n_{f}=0$ and the other 14 levels at $\kappa=0$ (with $n_{f}=0,\pm 1,\pm 2)$, which all contribute to the ground state manifold in the Dirac CFT. For the chiral Ising CFT, the $n_{f}=\pm 1$ states form an 8-fold degenerate level, while we expect the 2-fold degenerate $n_{f}=\pm 2$ and the two 2-fold degenerate $n_{f}=0$ levels to form a 6-fold (quasi-) degenerate level.

The other large degeneracies of the higher levels in the Dirac CFT similarly split up in the chiral Ising CFT. Furthermore, very characteristically, we observe a very low-lying two-fold, nearly degenerate set of levels with the same quantum numbers as the ground state levels, similar to what was also observed in Wilson-Fisher CFTs [17, 19]. The quasi-degenerate nature of the ground state levels also transfers to the lowest $n_{f}=\pm 1$ levels at $\kappa=1$. Their energies are pushed to values slightly above the linear dispersion relation with the velocity $v_{\text{F}}^{c}$, where our best estimate for $v_{\text{F}}^{c}$ was obtained as described in Sec. III.4. We expect that such a two-fold quasi-degeneracy appears for all the lowest $n_{f}=\pm 1$ levels at all $\kappa>0$, however, these were inaccessible because of finite-size restrictions in the numerical calculations. A comparison of the critical torus spectrum for the chiral Ising universality class in Fig. 4(b) with the critical torus spectrum for the (Wilson-Fisher) Ising universality class [17, 19] where, in particular, only non-degenerate levels appear, shows the strong influence of the fermionic degrees of freedom on the critical spectrum. This, once more, demonstrates the potential of the critical torus energy spectrum as a useful identifier for universality classes.

III.3 Crossover Effects near the QCP

For all values $V<V_{c}$ the microscopic models, in the thermodynamic limit, flow towards the Dirac CFT fixed point with massless fermionic excitations featuring a linear light-cone, and the spectrum is defined by the universal numbers $\xi_{i}^{\text{D}}$. The energy spectrum on finite clusters, however, is affected by a pronounced crossover effect: This derives from the fact that in the vicinity of the QCP at $V_{c}$, the RG flow is first attracted towards the chiral Ising CFT fixed point on intermediate length scales, and later crosses over to the asymptotic Dirac CFT fixed point only beyond an increasingly larger length scale $L_{c}$, which diverges upon approaching the QCP. For $V$ near $V_{c}$, sufficiently large system sizes $L\gg L_{c}$ are thus required in order to probe the asymptotic Dirac CFT fixed point. As a result, the values of the scaled excitation gaps $\Delta_{i}\times L$ exhibit a continuous crossover between the asymptotic values and those at the QCP, in particular for those levels, for which $\xi^{\text{D}}_{i}$ and $\xi^{\text{CI}}_{i}$ differ notably (in particular levels at $\kappa=0$). This crossover behavior is illustrated in Fig. 5.

This may be made precise by using the theory of finite-size scaling. Assume we have a $D$-dimensional CFT perturbed by a single relevant operator, with an associated correlation length exponent $\nu>0$. We also perturb our CFT with any number of irrelevant operators, and call the usual critical exponent governing the leading corrections to scaling $\omega>0$. The above conditions should describe a typical critical point with a single relevant direction being probed by an experiment or numerical simulation. With these definitions, the finite-size spectrum of the perturbed CFT on the torus with “speed of light” $v$ and linear extent $L$ is given by [39]

[TABLE]

In this expression, the dimensionless scaling functions $\xi_{i}[x]$ and $\zeta_{i}[x]$ are universal up to overall multiplicative factors and a normalization of their arguments. We have included the non-universal length scale $\delta$ associated with the addition of the leading irrelevant operator to the CFT, and we only show the leading non-analytic dependence on $L$.

In any lattice model, the spectrum calculated numerically will contain all of these terms, but in this paper we focus on extracting the constants $\xi_{i}=\xi_{i}[0]$, which completely characterize the torus spectrum of the unperturbed CFT. Therefore, we are really interested in taking the limit

[TABLE]

By comparing this limit to Eq. (17), we see that the first inequality ensures that only $\xi_{i}[x]$ contributes, while the second ensures that the limit $\xi_{i}[x\rightarrow 0]$ is taken. In analyzing numerical or experimental results, where the non-universal scale $\delta$ may be anomalously large or $V-V_{c}$ cannot be tuned with arbitrary precision, one should always check that these conditions hold 444We additionally note that if $\omega$ happens to be very small numerically, then corrections to scaling will remain important for large system sizes. This does not seem to be a problem for our cases of interest; for the Dirac CFT, $\omega_{\text{D}}=1$, while for the chiral Ising CFT, $\omega_{\text{CI}}\sim 0.8$ at four-loops in the $\epsilon$-expansion [10]..

Applying this reasoning to the chiral Ising CFT, we find that, even if we tune into the semimetal phase, $V<V_{c}$, our finite-size spectrum will continue to be that of the chiral Ising CFT provided the linear extent of the torus satisfies $L\ll L_{c}=|V-V_{c}|^{-\nu_{\text{CI}}}$. Alternatively, we may apply this analysis to the Dirac CFT, which does not have any relevant operators. Instead, the coupling $V$ is actually an irrelevant perturbation to the Dirac CFT, so the length scale $L_{c}$ should actually be associated with $\delta$ in Eq. (18), and the Dirac CFT spectrum is obtained when $L\gg L_{c}$. At intermediate length scales, $L\sim L_{c}$, the torus spectrum is described by the full function $\xi_{i}[x]$.

III.4 Quantifying the Fermi Velocity Renormalization

In the SM phase, the considered model systems feature massless fermionic excitations in the thermodynamic limit, with a linear single particle dispersion and a renormalized Fermi velocity $v_{\text{F}}(V)$ which depends on the interaction strength in a non-universal manner. This is due to a RG flow towards the Dirac CFT fixed point for all values $V<V_{c}$. Only exactly at the QCP will the system finally flow to the chiral Ising CFT fixed point. While the shape of the fermion dispersion relation does not change within the SM phase, the Fermi velocity $v_{\text{F}}(V)$ can differ greatly from the non-interacting value $v_{\text{F}}^{0}$ and is model dependent [see Fig. 5(b)].

Since the torus energy spectrum is a universal property for the underlying Dirac CFT, the energy levels are given by $\Delta_{i}L=v_{\text{F}}(V)\xi^{\text{D}}_{i}$ within the SM phase, where the $\xi_{i}^{\text{D}}$ denote universal numbers describing the SM phase and do not depend on the interaction strength $V<V_{c}$. Hence, it is possible to determine the renormalization of the Fermi velocity from single energy levels (with $\xi_{i}^{\text{D}}>0$) of the spectrum,

[TABLE]

Here, we obtain $v_{\text{F}}(V)$ by measuring the energy gap $\Delta_{n_{f}=\pm 1}(\mathbf{k}_{\text{min}})$ of the single-fermion excitation $n_{f}=1$ at the momentum $\mathbf{k}_{\text{min}}$ closest to the Dirac point (family H), which shows particularly small finite-size and crossover effects. This momentum corresponds to $\kappa=1$, and the corresponding value of $\xi^{\text{D}}_{i}$ is exactly known [see Fig. 4(a)]. Note, that Eq. (19) cannot be readily applied at $V=V_{c}$ to extract the critical Fermi velocity $v_{\text{F}}^{c}$, since the values of $\xi^{\text{CI}}_{i}$ are not a priori known, and cannot be measured independently of $v_{\text{F}}^{c}$ in numerical simulations.

The results obtained by the above analysis for the renormalized Fermi velocity $v_{\text{F}}(V)$ in the SM regime of the honeycomb lattice model $H_{\rm h}$ is shown in Fig. 6, along with a linear regression, which describes the numerical results remarkably well. Assuming a non-singular behaviour of the Fermi velocity across the critical point [3], we extrapolate the linear functions to $V_{c}$ to obtain an estimate for the critical Fermi velocity $v_{\text{F}}^{c}$ which is approximately $35\%$ larger than $v_{\text{F}}^{0}$. Note that in the GNY field theory the speed of light is the analogue of the Fermi velocity, and stays constant due to strict Lorentz invariance throughout all the phases.

Furthermore, we perform the same analysis for a gapped level with $n_{f}=\pm 1$ on clusters in family H’, which yields a velocity renormalization in very good agreement with the previous estimate [see Fig. 6].

Alternatively, one may be tempted to consider the slope of the dispersion in the vicinity of the Dirac point,

[TABLE]

to extract the Fermi velocity renormalization: In fact, based on Eq. (14), we obtain $v^{\text{slope}}_{\text{F}}$ in terms of the universal numbers $\xi^{\text{D}}_{i}$ for the $n_{f}=1$ excitations at $\kappa=0$ and $\kappa=1$, respectively. These can be readily extracted from Fig. 4(a), and we obtain $v^{\text{slope}}_{\text{F}}(V)=v_{\text{F}}(V)$ for $V<V_{c}$. This estimator is, however, strongly influenced by the crossover effects, mentioned in Sec. III.3, close to the QCP: As seen from Fig. 4, the single-particle excitation at the Dirac point (with $n_{f}=1$, $\kappa=0$), which enters the estimator for the Fermi velocity renormalization in Eq. (20), shows particularly different values of $\xi^{\text{D}}_{i}$ and $\xi^{\text{CI}}_{i}$. Hence, performing numerical simulations on insufficiently large tori (due to limitations in the accessible system sizes), will lead to a false estimate for the Fermi velocity renormalization by Eq. (20).

An extrapolation of the Fermi velocity based on slopes between the Dirac point and the closest momentum nearby, as in Eq. (20), is thus particularly dangerous in the vicinity of interacting quantum critical points and a very careful analysis including proper finite-size scaling is necessary [41, 17]. As pointed out recently [24], such a crossover due to enhanced finite-size shifts in the excitation gap at the Dirac point led the authors in Ref. [23] to drastically underestimate the Fermi velocity in the Hubbard model on the honeycomb lattice near its QCP. While the nature of the quantum critical point is different in Ref. [23] (chiral Heisenberg fixed point), the problem of the crossover length scale also applies there.

IV Numerical Results for the Lattice Models

In this section, we present our numerical results for the microscopic lattice models in more detail. We begin by providing an overview of the two methods that we used for our numerical calculations, exact diagonalization (ED) and quantum Monte Carlo (QMC), and explain how we corrected for warping effects from the lattice discretization, before providing details of the extrapolation to the thermodynamic limit.

IV.1 Exact Diagonalization

Exact diagonalization (ED) [42, 43] can be used to calculate all low-energy gaps directly and exactly for all parameters $V$ on finite size clusters. In addition, their quantum numbers according to fermion-number conservation and lattice symmetries can be directly identified using a symmetry-adapted basis. In particular, the ED spectrum is divided into $n_{f}$ sectors combined with a $\mathbb{Z}_{2}$ charge according to the particle-hole symmetry $C$ at half filling $n_{f}=0$, as well as the momentum quantum number $\mathbf{k}$, and an irreducible representation of the lattice’ point-group. This also allows for the identification of appropriate quantum many-body operators with non-vanishing matrix elements between the ground state and the various low-lying excitations. These operators can then be used to extract the corresponding energy gaps within the QMC simulations from the decay of the imaginary-time correlation function, as described in the next section.

IV.2 Quantum Monte Carlo

We employ the projector lattice continuous-time quantum Monte Carlo algorithm (LCT-INT) detailed in [44]. The ground state expectation value of an observable $\hat{O}$ is accessed upon projecting a trial wave function $|\Psi_{T}\rangle$,

[TABLE]

where $|\Psi_{0}\rangle$ denotes the ground state of the Hamiltonian $H$. For this work, the simulations were performed with a projection length of up to $\Theta=160/t$ to ensure convergence within the statistical uncertainty. Importantly, the LCT-INT formulation does not rely on a Trotter decomposition, but instead decomposes the projection operator using an interaction expansion directly in continuous time, thus eliminating the Trotter error completely.

The trial wave function is chosen as a zero momentum, particle-hole (anti-)symmetric ground state of the free Hamiltonian, and is represented by a Slater determinant $|\Psi_{T}\rangle$. Furthermore, the invariance of the Hamiltonian under the reflection symmetry $I_{x}$ can be used to separate the two quasi-degenerate ground states of the interacting system. We therefore consider trial wave functions $|\Psi_{T}\rangle_{\pm}$ with $I_{x}$-eigenvalue $\pm 1$ in order to project onto the ground state of the corresponding symmetry sectors.

The lowest energy gaps $\Delta_{\hat{O}}$ are extracted from the asymptotic decay of imaginary-time correlation functions, dominated by

[TABLE]

where for sufficiently large $\tau$, with $1\ll\tau\ll\Theta/2$, this leading exponential decay is dominated by the smallest gap $\Delta_{\hat{O}}$ accessible by the operator $\hat{O}$. The relevant energy gaps correspond to excited states in different symmetry sectors, which are determined by their fermion-number, momentum, particle-hole symmetry, and an irreducible representation of the lattice point-group. The operators $\hat{O}$ have to connect the ground state $|\Psi_{0}\rangle$ to the desired excited state $|\Psi_{\Delta_{\hat{O}}}\rangle$ such that the overlap $|\langle\Psi_{0}|\hat{O}|\Psi_{\Delta_{\hat{O}}}\rangle|$ is finite. Feasible operators $\hat{O}$ can be categorized by their action under the various symmetry operations. For example, states with opposite parity under $I_{x}$ are connected by operators for which the anti-commutator with the reflection operator $I_{x}$ vanishes, $\{\hat{O},I_{x}\}=0$. A detailed list of the symmetry properties of the operators connecting the ground state to the various relevant excited states can be found in Table 4. The explicit expressions of possible operator implementations can be found in Appendix A.

IV.3 Warping Corrections

It is well known that the Dirac cones for the non-interacting lattice models, $V=0$, are strongly modified by warping effects from the asymptotic dispersion relations, see Fig. 2. As a result, rather large systems are required to directly probe the linear dispersion regime, especially for those clusters that do not contain the Dirac points. In order to reduce the finite-size effects for the energy levels measured on clusters without Dirac points, we appropriately rescale the finite-energy spectra at any momentum $\mathbf{k}$ in the vicinity of the Dirac points for all interaction strengths $V$ as

[TABLE]

where $v_{\text{F}}^{0}$ denotes the Fermi velocity for $V=0$, $\varepsilon^{0}(\mathbf{k})$ the dispersion relation of the non-interacting system, and $\mathbf{D}$ is equal to $\mathbf{K}$ or $\mathbf{X}$, for $H_{\rm h}$ and $H_{\rm s}$, respectively. This unwarping of the energy spectrum is illustrated in Fig. 7. In the following, $\hat{\Delta}_{i}$ will always indicate that an unwarping according to the above equation has been performed.

IV.4 Results - Evolution of the Torus Spectrum of $H_{\rm h}$ with $V$

In this section, we provide a detailed analysis of the spinless fermion $t-V$ model on the honeycomb lattice, described by Eq. (5) at half-filling. The underlying structure of the energy level spectrum is uncovered upon appropriately rescaling the finite-size energy gaps by a factor of the linear system scale, as quantified by $L=\sqrt{N_{s}/2}$. We, therefore, consider in the following the low-energy gaps rescaled as $\Delta_{i}\times L$, which we call the spectrum. We first consider the evolution of the full low-energy spectrum with the interaction strength $V$, including the free system $V=0$ and the QCP at $V_{c}\approx 1.355$ [33, 4], based on ED calculations on clusters of a few ten sites. This is very instructive in order to identify the qualitative features of a QCP that remain valid in the thermodynamic limit, even though the quantitative values of the energy gaps may be subject to substantial finite-size effects.

Fig. 8(a) shows the evolution of the spectrum on the ${N_{s}=24}$ sites cluster, which is characteristic of tori that contain the Dirac points (family H). In the non-interacting case, ${V=0}$, we identify a six-fold degenerate ground state in the half-filled sector, ${n_{f}=0}$. The lowest single-fermion level ${n_{f}=\pm 1}$ 555Due to particle-hole symmetry the spectrum is identical for ${n_{f}=\pm\nu,\,\nu\in\mathbb{Z}_{0}}$, i.e., if we add $\nu$ fermions or holes. Therefore, we only show the $n_{f}\geq 0$ sectors. (8-fold degenerate) and the two-fermion sector ${n_{f}=\pm 2}$ (two-fold degenerate; not shown) are also gapless, as fermions can be created at the gapless Dirac points. The finite-size spectrum immediately gaps out for finite ${V>0}$, and the system undergoes a transition into the CDW phase for ${V>V_{c}}$, where a two-fold (quasi-)degenerate ground state of a $\mathbb{Z}_{2}$ even and a $\mathbb{Z}_{2}$ odd level is observed, while all fermionic excitations ${n_{f}\neq 0}$ are gapped.

According to their quantum numbers we can label the energy levels and relate them to the known instabilities of the Dirac phase: Mass gaps for spinless Dirac fermions can be generated by breaking the sublattice symmetry ($\sigma_{T}$), by breaking the time-reversal symmetry with zero net magnetic flux through the honeycomb unit cell (Chern), or by a Kekulé dimerization which creates two distinct real masses [46, 47, 48, 49].

Another prominent level near the quantum critical point is the state corresponding to a detuning from the quantum critical point ($\epsilon_{T}$). This level lies in the same symmetry-sector as the ground state ($1_{T}$) and typically shows a characteristic shape with a minimum around the QCP. This level is also the leading contribution to the fidelity susceptibility at the quantum critical point [50, 51].

The critical torus spectrum in the past showed qualitatively similar structures as the operator content of the corresponding field theory, i.e., the scaling dimensions of the fields of the GNY CFT [17, 19]. We have, therefore, chosen the labels $\sigma_{T}$ ($\epsilon_{T}$) as the torus analogues of the lowest particle-hole odd (even) scalar fields and $\psi_{T}$ as the torus analogue of the lowest vector field (the lowest single-fermion excitation ${n_{f}=\pm 1}$). Furthermore, we label the lowest two-fermion excitation ${n_{f}=\pm 2}$ as $(2\psi)_{T}$ and the torus analogue of the fermionic descendent field, i.e., the lowest fermionic excitation at $\mathbf{k}_{\text{min}}$, as $(\psi+d)_{T}$. A prime on a level symbol is used to indicate the second-lowest level in the same symmetry sectors as the corresponding unprimed level. In Tab. 1 we list the most important quantum numbers for these levels.

The finite-size spectrum for the family of clusters that do not include the Dirac points is structurally different in the SM phase [see Fig. 8(b)]. The ground state is unique even for ${V=0}$ and the fermion excitations ${n_{f}=\pm 1}$ and ${n_{f}=\pm 2}$ are gapped. The $\sigma_{T}$ field strongly decreases in energy with a finite value at the critical point and constitutes the second state in the two-fold degenerate ground state manifold for ${V\gg V_{c}}$, but, in contrast to the clusters with Dirac points has a large gap at $V_{c}$. Again, the $\epsilon_{T}$ field shows a very characteristic shape with a strongly reduced gap only around the critical point.

IV.5 Results - Torus Spectrum at Criticality

We next examine in more detail the spectrum at criticality for the two microscopic lattice models, $H_{\rm h}$ and $H_{\rm s}$. In Fig. 9 we show the critical torus energy spectra, as obtained from the different cluster geometries. In this figure, the critical torus energy spectra are rescaled by the critical Fermi velocity, which, up to a global factor, identifies the universal numbers $\xi_{i}^{\text{GNY}}$ for the $N=4$ chiral Ising GNY universality class in $D=(2+1)$ dimensions. As mentioned in Sec. II.4, one furthermore has to distinguish for each model different families of finite size clusters. For the honeycomb lattice, the families H (H’) are distinguish by the presence (absence) of the Dirac points in their momentum space [c.f. Fig. 3(a)]. For the square lattice, we distinguish three families: Clusters in family S contain the two Dirac points among the lattice momenta, while for those belonging to family S’ (S”), the Dirac points are (not) located at the center between four lattice momenta, respectively [c.f. Fig. 3(b)].

To obtain the critical torus spectrum in the thermodynamic limit we extrapolated the finite-size results $\Delta_{i}\times L$ for the different levels [c.f. also Fig. 8] obtained from ED and QMC to $N_{s}\rightarrow\infty$, as shown in Fig. 10 for the model $H_{\rm h}$, on which we focus in the remainder of this section. The details of the corresponding finite-size analysis of the torus spectrum at the QCP of the square lattice model $H_{\rm s}$ are provided in Appendix B. While the quality of the extrapolations is not equally good for all levels since not all energy gaps could be measured with QMC, it is important to note that the qualitative structure of the low-energy levels, their quantum numbers and (approximate) multiplicities are already present on the smaller clusters with a few tens of sites, and it is this qualitative structure that serves as a fingerprint for the chiral Ising universality class. For the same reason we also omit to give error bars on the extrapolated levels.

The critical torus energy spectrum for tori that contain the Dirac points (family H) show very characteristic features [c.f. Figs. 9, 10]: The gap for the $\sigma_{T}$ field is remarkably small and appears to form a two-fold degenerate state together with the ground state, i.e., the vacuum level $1_{T}$ in the thermodynamic limit. The $\epsilon_{T}$ level and the $\sigma_{T}^{\prime}$ level are also close to each other and build a second copy of such a two-fold (nearly) degenerate level in the thermodynamic limit. The fermion mode $\psi_{T}$ is the next lowest level above $\sigma_{T}$. The subsequent Kekule, Chern and $(2\psi)_{T}$ levels are very close to each other and build a six-fold (nearly) degenerate level. Two single-fermion levels $(\psi+d)_{T}^{(\prime)}$ with momentum $\mathbf{k}_{\text{min}}$ are also found to build a nearly degenerate set of levels with energy comparable to the $\epsilon_{T}$ and $\sigma_{T}^{\prime}$ states. The characteristic two-fold nearly degenerate levels are a result of the two-fold nearly degenerate ground state.

For tori that do not contain the Dirac points (family H’), the critical torus energy spectrum is strongly altered, with a much larger $\sigma_{T}$ gap, followed by the $\psi_{T}$ field. The opening of the $\sigma_{T}$ gap also leads to a strong splitting of the other two-fold nearly degenerate levels observed on the clusters with Dirac points. The $(2\psi)_{T}$ and Chern fields are, again, very close to each other, while the Kekule field is not clearly defined. The $\epsilon_{T}$ level seems to be only slightly influenced by the choice of the torus shape.

IV.6 Results - Estimation of the Fermi Velocity Renormalization

We already discussed our approach to estimate the Fermi velocity renormalization and the subtle crossover effects near the QCP in Secs. III.4 and III.3. Here, we provide further details on this procedure for $H_{\rm h}$, while the case of $H_{\rm s}$ is treated in Appendix B. In particular, we consider using different levels to access $v_{\text{F}}(V)$ based on the general formula Eq. (19) within the SM phase. Here, we focus on the Hamiltonian $H_{\rm h}$, for which larger-system QMC data is available.

In Fig. 11(a) we show various renormalized energy levels $\Delta_{i}\times L$ within the SM phase of the honeycomb lattice model $H_{\rm h}$ from finite clusters that do not contain the Dirac points, and also extrapolated these values to the thermodynamic limit $N_{s}\rightarrow\infty$. Because of the very small finite-size effects of the $n_{f}=\pm 1$ level ($\psi_{T}$), the Fermi velocity renormalization is best read off from this level and its $V$-dependence can be well approximated by a linear function [see center panel in Fig. 11(a)]. Furthermore, the renormalization of other levels with stronger finite-size effects [see left and right panels in Fig. 11(a)] is well approximated by the same function in the SM phase. It is important to note, that the strong drop of the $\sigma_{T}$ level close to the critical point is not due to a sudden, strong decrease of the Fermi velocity, but mainly because of the strong difference of the universal numbers $\xi^{\text{D}}_{\sigma_{T}}$ in the SM phase and $\xi^{\text{GNY}}_{\sigma_{T}}$ at the chiral Ising critical point. This level, thus, provides us with another dramatic example of the crossover effects that have been discussed in Sec. III.3.

For clusters that contain the Dirac points, many gaps, such as the single particle gap at the Dirac point ($\psi_{T}$), vanish faster than $1/L$ in the SM phase, i.e., these levels have a vanishing value of $\xi^{\text{D}}_{i}=0$, and they, thus, cannot be used to extract $v_{\text{F}}(V)$ based on Eq. (19). On these clusters, we, therefore, measure the gap of a single fermion excitation $n_{f}=1$ with $\kappa=1$, labeled $(\psi+d)_{T}$, which is gapped, and obtain the $v_{\text{F}}(V)$ renormalization from this value [see Fig. 11(b)]. The extrapolated values and the linear regression function is also shown in Fig. 6 and was further discussed in Sec. III.4. Again, we obtain an approximately linear behaviour and, importantly, the regression functions agree very well among the two families of finite-size clusters. The agreement of $v_{\text{F}}(V)$ as extracted using different levels and cluster families demonstrates that this approach of computing the Fermi velocity renormalization within the SM phase is quite reliable.

V Torus Spectrum in the $\epsilon$ Expansion

In this section, we provide the details of the analytical $\epsilon$ expansion to extract the torus spectrum for the GNY field theory.

V.1 General Structure of the Expansion

We want to examine the finite-size spectrum of the GNY field theory given in Eq. (1). For this purpose, we will use a real-time Hamiltonian formulation, with

[TABLE]

These operators satisfy the equal-time commutation relations

[TABLE]

where the Dirac field $\Psi_{a}^{j}$ has spinor index $a=1,...,n_{\text{D}}$ and flavor index $j=1,...,N_{f}$, and we define $N=n_{\text{D}}N_{f}$ to be the total number of degrees of freedom. In Eq. (24) we have assumed that the tuning parameter $s$ has already been set to its critical value $s_{c}$, and used the fact that $s_{c}=0$ in dimensional regularization. At leading order in $\epsilon=3-d$, the interaction couplings flow to the fixed point values [5]

[TABLE]

As a reminder of our notation, we use complex coordinates for the two-dimensional torus, $x=x_{1}+ix_{2}$, and then define the torus by two complex periods, $\omega_{1}$ and $\omega_{2}$, such that the points $x+n\omega_{1}+m\omega_{2}$ are equivalent for all $n,m\in\mathbb{Z}$. The torus is characterized by its modular parameter, $\tau\equiv\omega_{2}/\omega_{1}=\tau_{1}+i\tau_{2}$, and its area is given by $\mathcal{A}=\mathrm{Im}\left(\omega_{2}^{\ast}\omega_{1}\right)$. In Eq. (24), we take the spatial integral to be over $d/2$ copies of the two-dimensional torus with modular parameter $\tau$, which preserves point group symmetries at all steps of the calculation, and does not introduce any extra unphysical parameters. In addition to the other symmetries mentioned in this paper, we note that the full torus spectrum is also invariant under the modular transformations, $\tau\rightarrow\tau+1$ and $\tau\rightarrow-1/\tau$, under which the torus area $\mathcal{A}$ is also invariant. The length scale $L$ introduced in earlier sections to define the universal numbers $\xi_{i}^{\text{CI}}$ are related to the area of the torus by $\mathcal{A}=|\mathbf{a}_{1}|\tau_{2}L^{2}$ (see Figure 1).

As mentioned in Section II.4, the structure of the $\epsilon$ expansion on the torus turns out to be very different if we allow twisted boundary conditions. In particular, if we consider torus clusters without Dirac points, this corresponds to a boundary condition $\Psi(x+n\omega_{1}+m\omega_{2})=e^{i\theta}\Psi(x)$ on the fermions, while the bosonic field remains fully periodic. The twisted boundary condition results in a finite-size mass gap for the fermions proportional to $\theta$, so one does not need to separate out the fermionic zero modes. Then following the arguments in Section III.A. of Ref. 19, whenever $\theta^{2}\gtrsim\epsilon^{1/3}$, the torus spectrum is given by an effective Hamiltonian which only involves the bosonic zero modes, implying that the torus clusters which do not have Dirac points will have a dramatically different spectrum from the $L\!\!\mod 3=0$ case corresponding to the field theory with periodic boundary conditions. In the remainder of this section, we will only focus on the fully periodic setup.

As emphasized in previous work [17, 19], the presence of massless bosonic fields invalidates naive perturbation theory in a finite volume. To obtain the finite-size spectrum, one must separate out the zero-momentum part of the fields, and subsequently treat the interactions of these zero-modes exactly. The non-perturbative treatment of the zero modes results in a torus spectrum which is dramatically different from the particle-like Fock spectrum of the Dirac case. To this end, we write the mode expansions for our fields as

[TABLE]

Here, $a^{\dagger}$, $b_{s}^{\,j\dagger}$, and $c_{s}^{\,j\dagger}$ create Fock states for bosons, fermions, and anti-fermions respectively (and the fermions have pseudospin and flavor indices $s$ and $j$). Dot products for complex coordinates are defined as $k\cdot x\equiv\mathrm{Re}\left(kx^{\ast}\right)$. The momentum sums are performed over the reciprocal lattice, which is given by

[TABLE]

where $k_{1}=-i\omega_{2}/\mathcal{A}$, $k_{1}=i\omega_{1}/\mathcal{A}$. The commutation relations for the zero-mode parts are

[TABLE]

We now place the mode expansion into Eq. (24) and separate the Hamiltonian into an unperturbed and interacting part, $H_{\text{GNY}}=H_{0}+V$, where we insist that all zero-mode operators are included in $V$. Then $H_{0}$ is a free Fock Hamiltonian,

[TABLE]

where $E_{0}=-(3/2)\sum_{k\neq 0}|k|$ is the leading contribution to the ground state energy, and repeated flavor indices are always summed from $j=1,...,N_{f}$. It is possible to compute the universal part of the ground state energy (the calculation for the Wilson-Fisher CFT is given in Ref. 19), but in this paper we will only compute the energy splittings from the ground state, and hereafter we subtract the ground state energy from $H_{\text{GNY}}$. The rest of the Hamiltonian is

[TABLE]

where we only show the terms needed to obtain the leading loop corrections to the spectrum.

We now treat $V$ as a perturbation to $H_{0}$. The spectrum of $H_{0}$ is just the Fock spectrum, but crucially, every state in the unperturbed spectrum is infinitely degenerate. This is because the zero mode $\varphi$ does not appear, so we may multiply each eigenstate by an arbitrary normalizable function of $\varphi$ without changing the energy. Each state additionally has a $2^{N}$-fold degeneracy due to the fermionic zero modes, since we may arbitrarily choose $\psi^{j\dagger}_{a}\psi^{j}_{a}=0,1$ for each value of $a$ and $j$.

We use an effective Hamiltonian method to treat $V$, which will describe the splitting of each Fock state due to interactions between the zero modes. We consider a degenerate subspace of $H_{0}$ with energy $\epsilon_{0}$, i.e., the set of states satisfying $H_{0}|\alpha_{0}\rangle=\epsilon_{0}|\alpha_{0}\rangle$. Then we construct an effective Hamiltonian which acts on this subspace, but whose eigenvalues are the exact eigenvalues, $H_{\text{eff}}|\alpha\rangle=E_{\alpha}|\alpha\rangle$, where $E_{\alpha}=\epsilon_{0}+O(V)$ are the exact eigenvalues of $H$. This effective Hamiltonian may be obtained perturbatively in $V$, and at leading order it is given by [52]

[TABLE]

where $P_{0}$ is the projection operator onto the degenerate subspace of interest. In this paper, we will only compute the effective Hamiltonian for the Fock vacuum, $|0\rangle$. In principle it is possible to obtain the effective Hamiltonian for any Fock state, but their structure becomes increasingly intricate at higher energies [19].

Taking $P_{0}=|0\rangle\langle 0|$, and combining Eqns (30), (31), and (32), we obtain $H_{\text{eff}}=|0\rangle\langle 0|h_{k=0}$, with

[TABLE]

The sum in the second line of this expression is ultraviolet divergent. The evaluation of sums of this form using dimensional regularization is treated at length in Appendix C of Reference 19, so we simply quote the result:

[TABLE]

Here, we define the function

[TABLE]

where the special function $\Theta(\lambda,\mathbf{\Omega})$, known as the two-dimensional Riemann Theta function, is defined as

[TABLE]

for a $2\times 2$ matrix $\mathbf{\Omega}$. The matrices appearing in Eq. (35) are

[TABLE]

We now discuss the spectrum of $h_{k=0}$, which acts on the space of zero modes. The eigenfunctions are a product of a bosonic and fermionic part, $F[\varphi]\otimes|n_{a}^{j}\rangle$, where the zero-mode operators act as (temporarily using hats to distinguish operators from their eigenvalues)

[TABLE]

Focusing on the fermionic part of the Hilbert space, we can show that the effective Hamiltonian is symmetric under a full U($N$) symmetry group. This can be seen by choosing a basis such that $\gamma^{0}=\mathrm{diag}(\mathbb{I},-\mathbb{I})$, after which the fermionic part of the Hamiltonian may be written

[TABLE]

Then by performing the transformation $\tilde{\psi}^{j}_{a+n_{\text{D}}/2}=\psi^{j\,\dagger}_{a+n_{\text{D}}/2}$, $\tilde{\psi}^{j\,\dagger}_{a+n_{\text{D}}/2}=\psi^{j}_{a+n_{\text{D}}/2}$, on the second term in the parenthesis, we have

[TABLE]

which is manifestly invariant under transformations of the ($\tilde{\psi}$) $\tilde{\psi}^{\dagger}$ fields as (anti-)fundamental vectors of U($N$). The emergent O($N$) symmetry of the chiral Ising CFTs noted at the end of Sec II.1 is the subgroup of this U($N$) obtained by taking purely real elements of the Lie group. The enlarged symmetry of the effective Hamiltonian compared to that in Eq. (24) occurs because the zero mode does not appear in the kinetic term, $\overline{\Psi}^{j}\cancel{\nabla}\Psi$, so the finite-momentum parts of $\Psi$ have less symmetry than the zero momentum part. Thus, we expect this extra symmetry of the zero-mode Hamiltonian to hold at all orders in perturbation theory. The results of Appendix D show that this emergent symmetry occurs in the $1/N$ and $\epsilon^{\prime}=D-2$ expansions as well.

Proceeding, we denote the above operator by $\hat{Q}\equiv\psi^{j\dagger}\gamma^{0}\psi^{j}$, and its eigenvalues $Q$ take integer values in the range $Q\in[-N/2,N/2]$. The degeneracy of the eigenvalue $Q$ is

[TABLE]

We now use Eq. (26) to write the critical couplings as $g_{\mathrm{Y},\mathrm{ren}}^{\ast 2}=Y\epsilon$ and $\lambda_{\text{ren}}^{\ast}=U\epsilon$, where $Y$ and $U$ only depend on $N$. After the canonical transformation $\varphi\rightarrow\epsilon^{-1/6}$ and $\pi\rightarrow\epsilon^{1/6}$, our final form for the effective Hamiltonian is

[TABLE]

This is the final version of the effective Hamiltonian that we will work with. The purpose of the canonical transformation was to make the $\epsilon$ dependence of the spectrum clear: the first line of Eq. (42) gives the leading $O(\epsilon^{1/3})$ contribution to the energy spectrum, and the second line (which required the computation of a one-loop diagram), gives the $O(\epsilon^{2/3})$ correction. The omitted terms in Eqns. (31) and (32) can be shown to contribute only at higher orders in $\epsilon$ (we direct the interested reader to Ref. 19 for details on this point).

The lowest energies of the critical GNY torus spectrum are given by numerically solving the Hamiltonians in Eq. (42) for each $Q$, and for a given $Q$ the set of states obtained have a degeneracy given by Eq. (41). Additionally, the spectrum of $h^{(Q)}_{k=0}$ is identical to the spectrum of $h^{(-Q)}_{k=0}$, where the bosonic part of the eigenfunctions are related by $F^{(Q)}[\varphi]=F^{(-Q)}[-\varphi]$. Thus, for a given $N$, we need to numerically solve the effective Hamiltonians for $Q=0,1,...,N/2$. Since larger values of $|Q|$ lower the minimum of the potential, we expect that the ground states of the system to be given by the ground states of the $Q=\pm N/2$ sectors. From Eq. (41), these two sectors are individually non-degenerate, and the resulting ground state is always exactly two-fold degenerate. We write the ground states as

[TABLE]

where the $\pm$ index indicates the eigenvalue of this state under the $\mathbb{Z}_{2}$ symmetry $\varphi\rightarrow-\varphi$. The states $|\mathrm{GS},+\rangle$ and $|\mathrm{GS},-\rangle$ correspond respectively to the levels denoted $1_{T}$ and $\sigma_{T}$ in earlier sections. The two ground states are exactly degenerate in the scaling limit, although we expect the state $|\mathrm{GS},-\rangle$ to acquire a gap in any lattice realization of the transition due to non-universal corrections to scaling.

The first excited state then corresponds to the effective Hamiltonians with $Q=N/2-1$ and $Q=-N/2+1$. From Eq. (41), this state has a total degeneracy of $2N$. The fermionic part of this state is obtained by acting on the ground state either with $\psi^{\dagger}_{a}$ for $1\leq a\leq n_{\text{D}}/2$, or with $\psi_{a}$ for $n_{\text{D}}/2+1\leq a\leq n_{\text{D}}$. These can be considered the particles ($n_{f}=1$) or holes ($n_{f}=-1$) in the language of previous sections, and these states clearly correspond to those labelled $\psi_{T}$ earlier.

We note that the lowest-lying finite momentum states are obtained by constructing an effective Hamiltonian around the zeroth order finite momentum Fock states [17, 19]. In the present model, the energy of these states are not shifted from the zeroth-order value $|\mathbf{k}|$ ($v_{F}=1$ here) until order $\epsilon$, so at the order we are working they are unchanged compared to the Dirac CFT. This is in agreement with the small shift of this level seen in numerics, see Figure 4.

We now detail the low-energy zero momentum spectrum for the $N=4$ case relevant to the model Hamiltonians $H_{\rm h}$ and $H_{\rm s}$, where we may relate each individual state with those obtained in numerics. A similar analysis of the spectrum may be done for any value of $N$.

V.2 The case of $N=4$

Specializing to the case with four degrees of freedom, we need to solve for the lowest eigenvalues of Eq. (42) for $Q=0,1,2$. We obtain the spectrum by first computing the eigenfunctions and eigenvalues of the $O(\epsilon^{1/3})$ part of the spectrum numerically, giving us the leading-order contribution to the energy and eigenfunctions. We then compute the $O(\epsilon^{2/3})$ contribution from these eigenfunctions using ordinary first-order perturbation theory. The results of this computation are shown in Tab. 2. In Tab. 3 we list the explicit numbers for the torus shapes considered in this manuscript.

By looking at the transformation of the eigenfunctions under the symmetries in Section II.3, we can explicitly relate these states to the $\kappa=0$ states enumerated in the numerical simulations of previous sections. In Fig. 4(b), we compare the resulting spectrum from $\epsilon$-expansion to the one obtained from numerics. We observe a favorable agreement between the two methods. In particular, the sequence of the eigenstates’ quantum numbers and degeneracies (quasi-degeneracies in numerics, see below) are identical. Also, the relative energy gaps in the $\kappa=0$ sector are similar, i.e. we observe large gaps between the $Q=\pm 2$ levels to the other states, while the $Q=0$ and $Q=\pm 1$ levels are very close in energy. Quantitatively, the $\epsilon$-expansion underestimates the gaps observed in numerics, which we relate to be mainly an artefact of the low-order expansion.

The degeneracies of the torus spectrum levels because of the emergent O($N$) symmetry of the chiral Ising CFT [see Secs. V.1, II.1] is a highly nontrivial prediction of the field theory, suggesting that the further splitting seen between these levels in numerics is non-universal and an artefact of corrections to scaling in explicit lattice realizations.

We also note that the energy depends on the shape of the torus rather weakly [Tab. 3], but that the levels for the square torus are slightly higher than those for the triangular torus.

As discussed earlier, the chiral Ising CFT with $N$ degrees of freedom appears to always flow to a fixed point with full O($N$) symmetry in perturbation theory, even when the original field theory does not possess this symmetry. We have already noted how this symmetry appears in the torus spectrum below Eq. (39), where it is a subgroup of a larger SU($N$) symmetry. Therefore, we may classify the states in Table 2 by their representations under these symmetry groups. From this perspective, the large degeneracies of the torus spectrum may be related to the large emergent symmetry of the CFT. Obtaining the relevant irreducible representations for a given state is easiest when $\hat{Q}$ is written in the form of Eq. (40), where the states are given by acting on the lowest-$Q$ state by antisymmetrized products of the SU($N$) vectors $\tilde{\psi}^{j\,\dagger}_{a}$. In this way, we see that the $Q=\pm 2$ states are SU(4) singlets, the eight $Q=\pm 1$ states are two inequivalent SU(4) vectors, and the six-fold degenerate $Q=0$ states transform into each other as an antisymmetric SU(4) tensor. The enumerations of multiplets and their degeneracies is not altered if we instead consider the O(4) subgroup of SU(4).

VI Conclusions & Outlook

In this manuscript we have shown how to calculate the critical torus energy spectrum for the chiral Ising universality class of the GNY theory with $n_{\mathrm{D}}=4$ spinor components in $D=(2+1)$ dimensions from the investigation of strongly interacting fermionic tight-binding models. We have computed the low-energy spectrum on finite-size clusters on different spatial torus geometries (honeycomb vs. square models) using exact diagonalization and quantum Monte Carlo approaches, which complement each other particularly well for this task. We have extrapolated the finite-size results to the thermodynamic limit to obtain the critical torus energy spectrum which serves as a unique fingerprint of the QCP’s universality class.

Furthermore, we have calculated the critical torus energy spectrum for the chiral Ising universality class using the perturbative expansion in $\epsilon=4-D$. This analytical approach shows a good qualitative agreement with our numerical results. In particular it predicts non-trivial degeneracies of levels which we also observe in numerics after extrapolation to the thermodynamic limit. This validates the description of the QCPs in the lattice models as chiral Ising critical points of the GNY theory. The $\epsilon$-expansion results also suggest, that the critical torus spectrum of GNY theories only depends on the total number of fermionic degrees of freedom $N=n_{\text{D}}N_{f}$, instead of depending on $n_{\text{D}}$ and $N_{f}$ individually.

We also observe that the finite-size clusters used to approach the thermodynamic limit split into families with distinct critical torus spectra. These families can be distinguished by the properties of the clusters momentum space; one family of clusters has the Dirac points in their momentum space and, in the field theory, correspond to periodic boundary conditions of the fermionic and bosonic fields. The other family does not have the Dirac points and describes twisted boundary conditions for the fermionic fields, while the bosonic field remains unaltered. The critical torus spectrum is very different for the different families.

Furthermore, we have computed the renormalization of the Fermi velocity due to the interactions between the fermions in the SM phase, which we derive from the renormalization of the torus energy level of a single fermion mode. We have shown that the so-obtained approximately linear velocity renormalization also describes the behaviour of other energy levels in the SM phase. Assuming that the Fermi velocity behaves continuously also at the critical point, we extrapolate the linear renormalization to the critical point, and obtain an approximately 35% increase at criticality compared to the non-interacting case. Additionally, we have investigated the crossover behaviour between the chiral Ising and the Dirac critical points for finite-size systems. We point out, that this crossover behaviour can lead to bad estimators for the Fermi velocity renormalization.

We hope that this work further strengthens the interpretation of the critical torus energy spectrum as a universal fingerprint of quantum critical points, that is, as we have seen here, capable of detecting the coupling of bosonic fields to fermionic spinors. We anticipate that this work inspires future research on the critical torus spectrum for chiral Ising models using different methods, a different number of spinor components $N$, or on chiral XY and chiral Heisenberg models where the spinors are coupled to a continuous O(2)/O(3)-order parameter. Recently, it was shown how to create fermionic tight-binding models with a single Dirac cone that is exactly linear in the entire Brillouin zone of the finite-size system [22]. Such systems could also be very beneficial to study chiral universality classes, because the portion of the Brillouin zone showing Dirac physics is much larger than in the models considered within this chapter, and finite-size extrapolation might become easier. Finally we believe our results complete an important step towards a quantitative understanding of the torus energy spectrum of QED3-like theories, believed to describe quantum spin liquids on the Kagome and the triangular lattice [25, 26, 27, 28].

Acknowledgements.

MS, TCL, and AML acknowledge support by the Austrian Science Fund for project SFB FoQus (F-4018) and DFG-FOR1807 (I-2868). MS acknowledges support by the Austrian Science Fund (FWF) through Grant No. P 31701-N27. SH and StW acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) under project number RGT 1995. SeW acknowledges support from the NIST NRC Postdoctoral Associateship award. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC). This work was supported by the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Focal Point Scientific Computing at the University of Innsbruck. We thank the IT Center at RWTH Aachen University and the JSC Jülich for access to computing time through JARA-HPC.

Appendix

Appendix A Operators for Gap Estimation with QMC simulations

The LCT-INT algorithm [44] used in this work makes explicit use of a weak coupling expansion of the partition function in terms of the interacting part of the Hamiltonian. As a result imaginary time displaced expectations values of the form in Eq. (22) can be expressed as expectation values of the free Hamiltonian. One can therefore employ the Wick theorem to calculate a Monte Carlo estimator of Eq. (22) using Green functions of the form

[TABLE]

which are calculated during the LCT-INT sampling process.

A.1 The $\sigma_{T}$ Level

The staggered density operator

[TABLE]

corresponds to the order parameter of the commensurate charge-density-wave of the bipartite lattice [see Fig. 12(b)]. The operator is anti-symmetric under particle-hole transformation and connects the $\mathbb{Z}_{2}$-even and -odd quasi-degenerate lowest energy states with each other. Furthermore, this operator provides an overlap of the ground state with the energetically higher $\sigma_{T}^{\prime}$ level.

A.2 The $\psi_{T}$ and $(\Psi+d)_{T}$ Levels

The $n_{f}=0$ sector can be connected to the $n_{f}=1$ ($n_{f}=-1$) sector by operators that create (annihilate) a fermion with a certain momentum. For lattice clusters that contain the Dirac point, the lowest excited state of the $n_{f}=1$ sector is connected to the $n_{f}=0$ ground state by the operator

[TABLE]

which creates a fermion at the Dirac point $\mathbf{K}$. Note that the eight-fold degeneracy of the $\sigma_{T}$ level follows from the valley, orbital and particle-hole degeneracy, $2_{\text{valley}}\times 2_{\text{orbital}}\times 2_{\text{PH}}=8$.

Overlap to higher excited states can then be achieved by creating a fermion with the $n$-th closest momentum to the Dirac point,

[TABLE]

where the momentum $\mathbf{q}$ is chosen accordingly on the lattice cluster.

A.3 The $2\psi_{T}$ Level

The $(2\psi)_{T}$ level corresponds to the lowest excited state in the $n_{f}=\pm 2$ sectors. The ground state can be connected to them by operators that create two fermions, one with momentum $\mathbf{K}$ and one with $-\mathbf{K}$,

[TABLE]

Note that the total momentum is zero. In this case, the valley and orbital degeneracies become redundant, and the two-fold degeneracy of the $2\psi_{T}$ level follows directly from particle-hole symmetry.

A.4 The $\epsilon_{T}$ Level

The $\epsilon_{T}$ level corresponds to the first excited state with identical symmetries as the ground state. In order to connect states within the same symmetry sector, possible operators must commute with all symmetry operations. A suitable choice is therefore given by either part of the Hamiltonian,

[TABLE]

Since both operators have a finite ground state expectation value, one has to extract the gap $\Delta_{\epsilon_{T}}$ using the formula

[TABLE]

Note that the operator $\hat{O}_{\epsilon_{T}}^{(V)}$ is related to the weak coupling expansion used in LCT-INT. One can therefore calculate its correlation function, as well as the fidelity susceptibility, from the distribution of interaction vertices during the Monte Carlo sampling [44].

A.5 The Kekule Level

The Kekule level corresponds to the lowest excited states in the $n_{f}=0$ sector with momentum $\mathbf{K}$ and identical particle-hole parity as the ground state. This level is two-fold degenerate due to the valley degeneracy. Possible operators can be constructed from the Kekule bond pattern, which itself is 3-fold degenerate on the honeycomb lattice ($K_{1}$, $K_{2}$ and $K_{3}$) [see Fig. 12(c)]. The Kekule pattern features an enlarged unit cell, which in reciprocal space corresponds to the momentum at the Dirac point. Because of the reduced lattice symmetry at finite momenta, states with momentum $\mathbf{K}$ do not have a well defined inversion parity. Nevertheless, one can choose to construct the Kekule operators such that they are (anti-)symmetric under lattice inversion. In this case the operators do not have well defined momenta, and provide overlap of the ground state with states of momentum $\mathbf{K}$ as well as $-\mathbf{K}$,

[TABLE]

Note that both $\hat{O}_{\text{Kekule}}^{(+)}$ and $\hat{O}_{\text{Kekule}}^{(-)}$ transform symmetric under the $\mathbb{Z}_{2}$ particle-hole transformation.

A.6 The Chern Level

The Chern level denotes the lowest excited states in the $n_{f}=0$ sector with zero momentum and opposite particle-hole parity as the ground state. These states are two-fold degenerate and transform according to a two-dimensional irreducible representation of the point group. They can be connected to the ground state by current operators that are anti-symmetric under particle-hole transformation and break rotational symmetry [see Fig. 12(d)], such as

[TABLE]

Appendix B Critical Torus Spectrum for $H_{\rm s}$

In this appendix we analyze the chiral Ising QCP of the model Eq. (6) on the $\pi$-flux square lattice. For this case, we accessed finite-size data only from ED and we show the finite-size extrapolations of the most important low-energy levels in Fig. 13. For the ED calculations for $H_{\rm s}$, we used the uniform gauge choice $\theta_{ij}=\pi/4$, $\theta_{ji}=-\theta_{ij}$, in order to assure a four-fold rotational symmetry of the Hamiltonian. The extrapolations of model Eq. (5) on the honeycomb lattice [see Fig. 10] show that the extrapolations of ED data alone typically give rather good agreement with the extrapolations of much larger clusters from QMC data. We, thus, assume that the extrapolations of the square-lattice model also give satisfactory qualitative estimates for the critical torus spectrum on the $\pi$-flux square lattice.

The critical torus spectrum for the clusters that contain the Dirac points [see left panel in Fig. 13] shows a very similar structure to the one of the honeycomb lattice, as it is also illustrated in Fig. 9 in the main text. Again, the $\sigma_{T}$ level shows a very low energy gap, and the Kekule, Chern and $(2\psi)_{T}$ levels are very close to each other, forming a $6$-fold (nearly) degenerate level in the thermodynamic limit. Furthermore, both the $\epsilon_{T}$ and $\sigma_{T}^{\prime}$ levels, as well as the $(\psi+d)_{T}$ and $(\psi+d)_{T}^{\prime}$ levels, become nearly degenerate, as was observed on the honeycomb lattice, and suggested from the $\epsilon$-expansion.

In contrast to the honeycomb case, the sequence of square lattices without Dirac points (and full four-fold rotational $C_{4}$ symmetry) separates into two families, which we denote as S’ and S”, respectively [see right panel in Fig. 13]. These two families are distinguished by the distance $||\mathbf{X}-\mathbf{k}_{\text{min}}||L$ of the closest momentum point $\mathbf{k}_{\text{min}}$ to the Dirac point $\mathbf{X}$ [see Fig. 3]. The scaled energy gap $\Delta_{i}\times L$ of, for example, the single particle mode $n_{f}=\pm 1$ is, thus, different for both families in the non-interacting case, $V=0$, which develops into a distinct critical spectrum for clusters of family S’ and S”. While the finite-size extrapolated $\epsilon_{T}$ levels are very similar for both of those families, the other levels are rescaled, but the sequence remains almost unchanged and is comparable with the results from the honeycomb lattice.

We can also attempt to compute the renormalization of the Fermi velocity with the interaction strength $V$ in the SM phase $V<V_{c}$ for the $\pi$-flux square lattice model in an analogous way to the honeycomb case [see Section IV.6]. We again choose the single-fermion level at $\mathbf{k}_{\text{min}}$ to derive the velocity renormalization, because it shows the smallest finite-size effects. We here use the largest finite size cluster to estimate the Fermi-velocity renormalization, since an extrapolation of the energy gaps to the thermodynamic limit for each $V<V_{c}$ could not be reliably obtained for this model. Like in the case of the honeycomb lattice, the renormalization can be well approximated by a linear function, and the behaviour of other characteristic levels, like the two-fermion level $n_{f}=\pm 2$, is consistent [see right panel in Fig. 14(a)]. The lowest $n_{f}=0$ excitation ($\sigma_{T}$) shows strong finite-size effects since it is strongly influenced by the crossover effects described in Sec. III.3, and, without QMC data, the available system-sizes are too small to read off the Fermi velocity renormalization from this level [see left panel in Fig. 14(a)]. The family of clusters that contain the Dirac points also yields a compatible renormalization of the Fermi velocity [see Fig. 14(b)]. At criticality, the Fermi velocity is increased by approximately $30\%$ compared to the non-interacting value $v_{\text{F}}^{0}$.

Appendix C Derivation of Quantum Field Theories from Model Hamiltonians

In this Appendix, we explicitly derive the GNY and GN field theories directly from the microscopic model Hamiltonians of Section II.2. We will mostly focus on the honeycomb lattice Hamiltonian $H_{\rm h}$; the corresponding derivation for $H_{\rm s}$ follows identical steps, and the only differences are unimportant numerical factors. This explicit computation is done to aid mapping the microscopic symmetries of our model Hamiltonians to those in the GNY and GN models, outlined in Section II.3.

We first recall the exact solution of the $V=0$ limit of the honeycomb model $H_{\rm h}$ [46] in order to fix notation before adding interactions. The honeycomb model has two sites per unit cell, and one may write the quadratic part of the Hamiltonian Eq. (5) as

[TABLE]

Here, the sum $\sum_{\mathbf{r}_{i}}$ is over all unit cells in the lattice, and ${\mathbf{s}_{1}=(1,0)}$, ${\mathbf{s}_{2}=(-1/2,\sqrt{3}/2)}$, and ${\mathbf{s}_{3}=(-1/2,-\sqrt{3}/2)}$ are the unit vectors connecting nearest-neighbors on the honeycomb lattice (we take the lattice spacing and the coupling $t$ to unity). We define the Fourier transform on the infinite honeycomb lattice by

[TABLE]

where the integral is over the Brillouin zone, and $\mathcal{A}_{\mathrm{BZ}}=8\pi^{2}/3\sqrt{3}$ is its area. This transforms the Hamiltonian to

[TABLE]

with $h(\mathbf{k})=-\sum_{i=1}^{3}\exp(\mathrm{i}\mathbf{k}\cdot(\mathbf{s}_{i}+\mathbf{s}_{1}))$. The resulting tight-binding spectrum is $E(\mathbf{k})=\pm|h(\mathbf{k})|$, and it is pictured in Fig. 2. It vanishes at two inequivalent points in the Brillouin zone, $\mathbf{K}$ and $\mathbf{K}^{\prime}$. Defining an appropriate four-spinor,

[TABLE]

with $\kappa=(2\pi)\mathrm{e}^{-i\pi/6}/\sqrt{\mathcal{A}_{\mathrm{BZ}}}$, the low energy effective Hamiltonian for $V=0$ becomes

[TABLE]

where $v_{\text{F}}=3/2$, and the explicit representation of the gamma matrices being used is given in Eq. (11).

We now consider interactions, $V>0$. We first define the scalar order parameter:

[TABLE]

This measures the difference in density on the two sublattices, so its expectation value should vanish in the SM phase, while it saturates to $\pm 1$ deep in the ordered phase. We note the following identity,

[TABLE]

We compare this to the interaction term in $H_{\rm h}$:

[TABLE]

where $h(\mathbf{q})$ is the same function appearing in Eq. (57), and we have dropped a constant. In the low-energy limit, we are only interested in the regions of these integrals where the fermions may be expanded in terms of the $\Psi$ fields of Eq. (58). Analyzing the first line of Eq. (62), we find that the only nonzero values of $\mathbf{q}$ which allow this expansion are $\mathbf{q}=\mathbf{K}$, $\mathbf{K}^{\prime}$, but these are precisely where $h(\mathbf{q})=0$, so these regions do not contribute. Thus, it suffices to expand the integrand near $\mathbf{q}=0$, where $h(\mathbf{q})\approx 3$. Then using Eq. (61), we find that

[TABLE]

We proceed by decoupling this term with a scalar field by a Hubbard-Stratonovich transformation, using the identity

[TABLE]

for every $\mathbf{r}_{i}$ after expressing our theory as a path integral. This implies that we may replace the interaction Hamiltonian by

[TABLE]

We note that $\phi$ transforms in the same way as the order parameter $\Phi$, and it will act as the order parameter in what follows. We also expect that $\phi$ has a uniform condensate in the ordered phase, so in the scaling limit we only need to consider its Fourier components close to zero momentum. Concentrating on the interaction term,

[TABLE]

the only regions of interest at low energy are $\mathbf{k}\sim 0$ and $\mathbf{k}^{\prime}\sim\mathbf{K},\mathbf{K}^{\prime}$. Using Eqs. (58)-(11) and going back to real space, we find

[TABLE]

where positive numerical constants have been absorbed into the definition of $V^{\prime}\propto V$.

To summarize, our full Hamiltonian is

[TABLE]

A renormalization group transformation will generate dynamical terms for the scalar field as well as all interactions allowed by symmetry, resulting in the GNY QFT of Eq. (1) modulo irrelevant terms. Alternatively, one may integrate out the $\phi$ field entirely at this stage and obtain the purely fermionic GN model of Eq. (4). This establishes that the phase transition of $H_{\rm h}$ lies in the chiral Ising universality class, and gives an explicit method for relating its symmetries to those in field theory.

Appendix D Spectrum in $1/N$ and $D-2$ Expansions

In this Appendix, we detail the spectrum of the chiral Ising universality class using the $1/N$ and $\epsilon^{\prime}=d-1$ expansions. A major benefit of the $1/N$ expansion is that one may easily obtain every state in the spectrum for both $g=g_{c}$ and for relevant deviations $g\neq g_{c}$ exactly in the $N=\infty$ limit. This allows us to compute crossovers from the critical point to the proximate phases. We will also show that the $N\rightarrow\infty$ limit of the torus spectrum commutes with the $\epsilon\rightarrow 0$ and $\epsilon^{\prime}\rightarrow 0$ limits to leading order.

D.1 $1/N$-expansion

We will find it convenient to work with the GN model, Eq. (4). After decoupling the interaction by a Hubbard-Stratonovich transformation, the imaginary-time Lagrangian is given by

[TABLE]

where we have defined $g=\tilde{g}/N$, expecting that $\tilde{g}$ remains finite for $N\rightarrow\infty$. As in the main text, we have $N_{f}$ flavors of $n_{\text{D}}$-component Dirac fermions, and define $N=n_{\text{D}}N_{f}$ to be the total number of degrees of freedom. We will take $n_{\text{D}}$ to be finite, taking $N_{f}$ (and therefore $N$) to infinity.

We may now integrate out the fermions exactly, obtaining the Euclidean action

[TABLE]

For $N\rightarrow\infty$, the path integral may be evaluated at its saddle-point. Defining $\Delta=\langle\phi\rangle$ in this limit, the saddle-point configuration is given by the gap equation

[TABLE]

which is to be solved for $\Delta$. Both sides of this equation are divergent and regularization-dependent, but if we consider deviations from the critical coupling, we obtain

[TABLE]

The left-hand side of Eq. (72) is finite for $1<d<3$. From this equation, we may solve for the energy gap $\Delta$ by writing

[TABLE]

where we define the special function

[TABLE]

and the Riemann Theta function $\Theta\left(\lambda,\mathbf{\Omega}(\tau)\right)$ was defined in Eqns. (36)-(37). We note that $\Delta$ is a monotonically increasing function of $g$, and that at $g=g_{c}$, it only depends on the shape of the torus. This gap equation is identical to that obtained in the large-$N$ limit of the Wilson-Fisher CFT [18].

We may now write the saddle-point Lagrangian as

[TABLE]

Because the gap equation only depends on $\Delta^{2}$, we have a sign ambiguity in choosing $\Delta=\langle\phi\rangle$. Unlike the Wilson-Fisher case [18], this sign difference is physical, and it has an effect on the torus spectrum. In particular, the full set of states are those obtained with positive $\Delta$ and those obtained for negative $\Delta$. But this sign has no effect on the spectrum of the theory, so we simply have two copies of a free massive Dirac spectrum with gap $|\Delta|$. Both sectors of the theory have the same ground-state energy, which may be calculated by temporarily taking a finite length $0<\tau<T$ in the Euclidean-time partition function $\mathcal{Z}_{\pm}(T)=\int\mathcal{D}\overline{\psi}\mathcal{D}\psi\exp(-\int_{0}^{T}d\tau\int d^{d}x\mathcal{L}_{\pm})$, and then calculating $E_{0}=-\lim_{T\rightarrow\infty}T^{-1}\ln\mathcal{Z}_{\pm}(T)$. This results in

[TABLE]

where in the second line we have used $g_{c}^{-1}=0$ in dimensional regularization, and the special function is defined in Eq. (74). We note that the ground state energy is actually ambiguous up to an overall constant, but our choice is such that the energy density vanishes at $g=g_{c}$ and $\mathcal{A}=\infty$.

To summarize, most of the states in the spectrum are described by the Hamiltonian

[TABLE]

where $s=1,...,n_{\text{D}}/2$ and $j=1,...,N_{f}$. In addition to the degeneracies for different values of $k$, $s$, and $j$, every state is doubled due to the different values of $\pm\langle\phi\rangle$. The ground state has energy $E_{0}$ and is two-fold degenerate, while the first excited state has energy $E_{0}+\Delta$ and degeneracy $2N$. We note that this Hamiltonian has an emergent SU($N$) symmetry, where the $N$ operators $b^{j}_{s}(k)$ and $c^{j}_{s}(k)$ together transform in the fundamental representation.

However, the Hamiltonian $H_{N}$ does not describe every state in the spectrum. The excited states which are singlets under the SU($N$) symmetry are instead described by fluctuations of the scalar mode $\phi$. We may obtain their energy by writing $\phi=\Delta+\tilde{\phi}/\sqrt{N}$ and expanding Eq. (70). To leading order in $1/N$, we obtain (ignoring a constant)

[TABLE]

where

[TABLE]

The function $\Pi(\omega,k)$ is the inverse Euclidean propagator of the $\tilde{\phi}$ field, so solving $\Pi(i\omega=E_{\phi}(k),k)=0$ for $E_{\phi}(k)$ gives the energy splitting of the states created by the scalar field from the ground state. The integrals and sums are convergent, and one can solve this equation numerically [19, 20]. Each solution obtained this way has a two-fold degeneracy corresponding to the sign choice $\pm\langle\phi\rangle$. Together with Eq. (77), this describes the complete spectrum at $N=\infty$.

D.2 $\epsilon^{\prime}$-expansion

We now consider the expansion in $d=1+\epsilon^{\prime}$ spatial dimensions, which we only study at leading order. The critical Hamiltonian associated with Eq. (4) is

[TABLE]

where $g^{\ast}=2\pi\epsilon^{\prime}/(N-2)$ to leading order in $\epsilon^{\prime}$ [5]. The computation proceeds as in Section V, but is simpler due to the lack of bosonic modes. We find the zero-mode effective Hamiltonian

[TABLE]

where $Q$ has the same definition as in Section V, leading to the same integer eigenvalues $Q\in[-N/2,N/2]$ and degeneracy Eq. (41). The Hamiltonian also has a $Q\rightarrow-Q$ symmetry, so the degeneracy analysis of the lowest two states is identical to that of the $\epsilon$ expansion. Extrapolating to $\epsilon^{\prime}=1$ for the $N=4$ case results in the predictions $\sqrt{\mathcal{A}}E_{1}=3\pi/2\approx 4.7$ and $\sqrt{\mathcal{A}}E_{2}=2\pi\approx 6.3$ for the energies of the first and second excited states respectively. States with higher energy have a spectrum given by effective Hamiltonians around higher Fock states.

The above expressions explicitly fail for the $N=2$ case. For $D=N=2$, our model is commonly referred to as the massless Thirring (Luttinger) model in the high energy (condensed matter) literature. This case is special as it does not have a phase transition; there is instead a single gapless Luttinger liquid phase [53]. It is also known that the $N=4$, $D=2$ model exhibits a Kosterlitz-Thouless transition at $g=0$ [54], which is not captured by the $\epsilon^{\prime}$ expansion. Therefore, as in the case for the O($N$) model, the expansion close to $D=2$ and small values of $N$ likely requires the consideration of non-perturbative effects to obtain the correct critical behavior, which would affect the accuracy of the $\epsilon^{\prime}$ expansion for smaller values of $N$ [55].

D.3 Correspondence between expansions

The large-$N$ expansion for the torus spectrum is much simpler than the $\epsilon$ or $\epsilon^{\prime}$ expansions in that the zero-momentum modes do not play a special role, so perturbation theory for the spectrum takes a similar form as perturbation theory for the critical exponents. Here, we show that all three expansions are compatible with each other, at least at leading order for low-lying states. We note that the large-$N$ gap equation, Eq. (73), may be solved in any dimension $d$, and in particular it may be solved analytically at leading order in either $\epsilon$ or $\epsilon^{\prime}$. An explicit computation finds the energy gap at criticality to be

[TABLE]

The rest of the spectrum is then given by Eq. (77) (with the exception of the singlet states, which we do not obtain analytically).

We now consider the large $N$ limit of the $\epsilon$ expansion derived in Section V. As argued there, the energy spectrum on the torus is given by solving the Hamiltonians $h_{k=0}^{(Q)}$ in Eq. (42) for $Q=0,1,...,N/2$, with the lowest energies given by values of $Q$ close to $N/2$. Therefore, to take the large-$N$ limit, it makes sense to take

[TABLE]

to concentrate on the low-energy spectrum.

At order $\epsilon^{1/3}$, the Hamiltonians $h_{k=0}^{(Q)}$ are shifted harmonic oscillators with a minimum at $\varphi=-(6\sqrt{Y}Q/U)^{1/3}$, which grows as $\sqrt{N}$ at large $N$. If we expand the Hamiltonians around this minimum with the assumption in Eq. (83), we find that they reduce to simple harmonic oscillators:

[TABLE]

Here, we have used the large-$N$ expansions of the couplings $U$ and $Y$. The anharmonic terms in the effective Hamiltonians vanish at $N=\infty$, so we may simply read off the low-energy spectrum from the well-known oscillator spectrum. The $q$-independent terms in the first line contribute to the universal part of the ground state energy, which is negative as expected for bosonic fields. With respect to the ground state, we see that we have two towers of free “particle-like” excitations, with energies given by

[TABLE]

The first excitation agrees with the large-$N$ gap given by solving the gap equation, Eq. (82). For the spectrum to match the large-$N$ results, this value of $\Delta_{2}$ must therefore be equal to the smallest solution of the equation $\Pi(i\omega=\Delta_{2},k=0)=0$, with $\Pi$ defined in Eq. (79).

The large-$N$ limit of the $\epsilon^{\prime}$ expansion is especially simple: we just take the limit Eq. (83) in Eq. (81). For $N\rightarrow\infty$, this gives a set of effective Hamiltonians

[TABLE]

This results in a tower of Fock states with mass $\sqrt{\mathcal{A}}\Delta=\pi\epsilon^{\prime}$, in agreement with Eq. (82). Since the $\epsilon^{\prime}$ expansion uses a purely fermionic theory, the bosonic singlet states will necessarily come from an effective Hamiltonian about a higher Fock state with an even number of fermions.