2. Nematic quantum criticality in Dirac systems

The results presented in this chapter are the outcome of a collaboration with Lukas Janssen, Kai Sun, Zi Yang Meng, Igor F. Herbut, Matthias Vojta, and Fakher F. Assaad. These findings have been published in [14], with significant portions reproduced here verbatim. My contributions to the project comprise the quantum Monte Carlo (QMC) and mean field calculations, including the implementation of the models in code and the creation of corresponding figures. The \(\epsilon\)-expansion has been performed by L. J., the models are designed by F. F. A. and K. S., while the interpretation of data and written text is a combined work of all authors.

2.1. Introduction

In a strongly correlated electron system, global symmetries, such as spin rotation, point group, or translational symmetries, can be spontaneously broken as a function of some external tuning parameter. This challenging problem has been studied extensively numerically and experimentally over the last years and impacts our understanding of quantum criticality [37] in cuprates [38] and heavy fermions [39]. The problem greatly simplifies when the Fermi surface reduces to isolated Fermi points in \(2+1\) dimensions and the critical point features emergent Lorentz symmetry. In this context, spin, time reversal, and translational symmetry breaking generically correspond to the dynamical generation of mass terms [40], and the semimetal-to-insulator transition belongs to one of the various Gross-Neveu universality classes [41, 42, 43, 44, 45, 46, 47].

Across nematic transitions, rotational symmetry is spontaneously broken [48, 49]. For continuum Dirac fermions with Hamiltonian \(H(\boldsymbol{k}) = v \left( k_x \sigma_x + k_y \sigma_y \right)\) in momentum space, where \(\boldsymbol{\sigma}\) are Pauli spin matrices and \(v\) is the Fermi velocity, nematic transitions correspond to the dynamical generation of nonmass terms, such as \(m \sigma_x\). They shift the position of the Dirac cone and as such break rotational, and therewith also Lorentz, symmetries. Such nematic transitions have been studied theoretically in the past in the context of \(d\)-wave superconductors [50, 51, 52, 53, 54] and bilayer graphene [55]. Fundamental questions pertaining to the very nature of the transition remain open: While initial renormalization group (RG) calculations based on the \(\epsilon\) expansion suggested a first-order transition [50, 51], a continuous transition has been found in large-N analyses [52, 53]. In this work, we use quantum Monte Carlo (QMC) simulations and a revised \(\epsilon\)-expansion analysis to study these transitions. We introduce two different models of Dirac fermions with twofold and fourfold lattice rotational symmetries, respectively. They are designed to feature a nematic transition, tuned by a parameter \(h\). The fermion dispersion of these models and the meandering of their Dirac cones are shown in Fig. 2.1(a). We demonstrate numerically and analytically that for both models this transition is continuous, realizing a new family of quantum universality classes in Dirac systems without emergent Lorentz invariance.

This chapter is organized as follows. In Sec. 2.2 we define the models and derive their symmetries. In Sec. 2.3 we perform a mean-field analysis of the models. In Sec. 2.4 we derive the continuum-field theories of the models. In Sec. 2.5 we perform an \(\epsilon\)-expansion analysis of these field theories, finding a continuous transition for both models, where App. A.1 contains more details on the analysis. In Sec. 2.6 we introduce the numerical method used for the QMC simulations and prove the absence of the negative sign problem. In Sec. 2.8 we define the QMC observables and show numerical results confirming a continuous transition for both models. In Sec. 2.9 we summarize our results. Furthermore, App. A.2 demonstrates how to use pyALF (cf. Section 4) to reproduce the QMC results of this chapter. App. A.3 and App. A.4 display the source code used for data collapse and retrieving the fermion dispersions, respectively. Finally, App. A.5 shows that the results from Sec. 2.8 do not change qualitatively when varying the number of spin degrees of freedom \(N_\sigma\) or the coupling parameter \(\xi\).

Fig. 2.1 (a) Contour plot of the fermion dispersion in the disordered phase from mean-field theory. The plotted area indicates the Brillouin zone. The green lines and arrows indicate the point group symmetries. Black (gray) dots sketch the meandering of the Dirac cones in the nematic phase for \(\langle \hat{s}_{\boldsymbol{R}}^z \rangle > 0\) (\(<0\)). (b)-(d) Fermion dispersion from QMC at \(L = 20\) for (b) \(h = 5.0 > h_\mathrm{c}\) featuring isotropic Dirac cones (c) \(h \simeq h_\mathrm{c}\), \(h \approx 3.27\) (left) and \(h \approx 3.65\) (right) featuring anisotropic Dirac cones, and (d) at \(h= 1.0 < h_\mathrm{c}\) featuring broken point-group symmetries. Color scale applies to all plots.

2.2. Models

Inspired from Refs. [49, 56, 57], we design two models of (2+1)-dimensional Dirac fermions, \(\mathcal{H}_0\), coupled to a transverse-field Ising model (TFIM),

(2.1)\[\begin{align} \mathcal{H}_\text{Ising} = -J \sum_{\langle \boldsymbol{R},\boldsymbol{R}' \rangle} \hat{s}_{\boldsymbol{R}}^z \hat{s}_{\boldsymbol{R}'}^z -h \sum_{\boldsymbol{R}} \hat{s}_{\boldsymbol{R}}^x, \end{align}\]

where \(\boldsymbol{R}\) denotes a unit cell and \(\langle \boldsymbol{R},\boldsymbol{R}' \rangle\) runs over adjacent unit cells. A Yukawa coupling, \(\mathcal{H}_\text{Yuk}\), between the Ising field and a nematic fermion bilinear yields the desired models, \(\mathcal{H}= \mathcal{H}_0 + \mathcal{H}_{\text{Ising}} +\mathcal{H}_{\text{Yuk}}\), that correspond to one of many possible lattice regularizations of the continuum field theories of Eqs. (2.18) and (2.19).

In the \(C_{2v}\) model, depicted in Fig. 2.2(a), we employ a \(\pi\)-flux Hamiltonian on the square lattice as

(2.2)\[\begin{split} \mathcal{H}_0^{C_{2v}} & = -t\sum_{\boldsymbol{R}} \sum_{\sigma=1}^{N_\sigma} \hat{a}_{\boldsymbol{R},\sigma}^\dag \Big( \hat{b}_{\boldsymbol{R},\sigma} e^{-i\tfrac{\pi}{4}} + \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_-,\sigma} e^{ i\tfrac{\pi}{4}} \\ &\quad + \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_--\boldsymbol{e}_+,\sigma} e^{-i\tfrac{\pi}{4}} + \hat{b}_{\boldsymbol{R}-\boldsymbol{e}_+,\sigma} e^{ i\tfrac{\pi}{4}} \Big) + \mathrm{H.c.},\end{split}\]

where \(\hat a\) and \(\hat b\) with spin index \(\sigma\) are fermion annihilation operators on the two sublattices, \(t\) is the hopping parameter, and \(N_\sigma=2\) is the number of spin degrees of freedom. \(\mathcal{H}_0\) features two inequivalent Dirac points per spin component in the Brillouin zone (BZ). The Ising spins \(\hat{s}_{\boldsymbol{R}}\) couple, with the sign structure indicated in Fig. 2.2(a), to the nearest-neighbor fermion hopping terms,

(2.3)\[\begin{split} \mathcal{H}_\text{Yuk}^{C_{2v}} & = -\xi \sum_{\boldsymbol{R}} \sum_{\sigma=1}^{N_\sigma} \hat{s}_{\boldsymbol{R}}^z \hat{a}^\dag_{\boldsymbol{R},\sigma} \Big( \hat{b}_{\boldsymbol{R},\sigma} e^{-i\tfrac{\pi}{4}} - \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_-,\sigma} e^{ i\tfrac{\pi}{4}} \\ &\quad - \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_--\boldsymbol{e}_+,\sigma} e^{-i\tfrac{\pi}{4}} + \hat{b}_{\boldsymbol{R}-\boldsymbol{e}_+,\sigma} e^{ i\tfrac{\pi}{4}} \Big) + \mathrm{H.c.},\end{split}\]

where \(\xi\) denotes the coupling strength.

Fig. 2.2 Sketch of (a) \(C_{2v}\) and (b) \(C_{4v}\) models, defined on \(\pi\)-flux single-layer and bilayer square lattices, with lattice vectors \(\boldsymbol{e}_{\nicefrac{+}{-}}\) and \(\boldsymbol{e}_{\nicefrac{x}{y}}\), respectively. Dark pink regions indicate unit cells, containing two orbitals (\(a\) and \(b\)) and one Ising spin (green arrow) in both cases. Fermions hop along the red lines and acquire a phase factor \(e^{i \pi/4}\) when following the direction of the arrow. Red and blue squares in (a) indicate the sign structure in the Yukawa coupling of the \(C_{2v}\) model.

The model has a \(C_{2v}\) point group symmetry, composed of reflections, \(\hat{T}_{\pm}\) on the \(\boldsymbol{e}_\pm = \boldsymbol{e}_x \pm \boldsymbol{e}_y\) axis. \(\hat{T}_{\pm}\) pins the Dirac cones to the \(\boldsymbol{K}_{\pm} = (\pi/2, \pm \pi/2)\) points in the BZ. Aside from the above reflections, \(\pi\) rotations about the \(z\) axis are obtained as \(\hat{T}_+ \hat{T}_- \). Further, the model exhibits an explicit \(\text{SU}(N_\sigma)\) spin symmetry that is enlarged to \(\text{O}(2N_\sigma)\) (cf. Section 2.6.1.1).

The \(C_{4v}\) model corresponds to a bilayer \(\pi\)-flux model, in which the Ising spins are located on the rungs, Fig. 2.2(b). The fermion hopping Hamiltonian is

(2.4)\[\begin{split} \mathcal{H}_0^{C_{4v}} & = -t\sum_{\boldsymbol{R}} \sum_{\sigma=1}^{N_\sigma} \hat{a}_{\boldsymbol{R},\sigma}^\dag \Big( \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_x,\sigma} e^{i\tfrac{\pi}{4}} + \hat{b}_{\boldsymbol{R}-\boldsymbol{e}_x,\sigma} e^{ i\tfrac{\pi}{4}} \\ &\quad + \hat{b}_{\boldsymbol{R}+\boldsymbol{e}_y,\sigma} e^{-i\tfrac{\pi}{4}} + \hat{b}_{\boldsymbol{R}-\boldsymbol{e}_y,\sigma} e^{ -i\tfrac{\pi}{4}} \Big) + \mathrm{H.c.},\end{split}\]

featuring four Dirac cones per spin component. The Yukawa coupling reads

(2.5)\[ \mathcal{H}_\text{Yuk}^{C_{4v}} = -\xi \sum_{\boldsymbol{R}} \sum_{\sigma=1}^{N_\sigma} i \hat{s}_{\boldsymbol{R}}^z \hat{a}^\dag_{\boldsymbol{R},\sigma} \hat{b}_{\boldsymbol{R},\sigma} + \mathrm{H.c.},\]

amounting to a coupling of the Ising spins to the interlayer fermion current. The \(C_{4v}\) Hamiltonian commutes with \(\hat{T}_{\pi/2}\), corresponding to \(\pi/2\) rotation about the \(z\) axis. The model is invariant under reflections \(\hat T_x\) and \(\hat T_y\) along the \(x\) and \(y\) axes, respectively. Reflections along \(\boldsymbol{e}_{\pm} = \boldsymbol{e}_x \pm \boldsymbol{e}_y\), denoted by \(\hat{T}_{\pm}\), can be derived from \(\hat{T}_{\pi/2}\), \(\hat{T}_x\), and \(\hat{T}_{y}\), and therefore also leave the model invariant. The model hence has a \(C_{4v}\) symmetry. Particle-hole symmetry, imposes \(A(\boldsymbol{k},\omega) = A(-\boldsymbol{k} + \boldsymbol{Q},-\omega)\), where \(\boldsymbol{Q}=(\pi,\pi)\) such that alongside with the \(C_{4v}\) symmetry the Dirac cones are pinned to the \(\pm\boldsymbol{K}_{\pm}\) points in the BZ (cf. Section 2.2.2).

2.2.1. Fourier transformed models

We define the Fourier transformation as:

(2.6)\[\begin{split}\begin{pmatrix} \hat{a}_{\ve{R},\sigma}^\dag \\ \hat{b}_{\ve{R},\sigma}^\dag \end{pmatrix} &= \frac{1}{\sqrt{N}} \sum_{\ve{k}} e^{-i\ve{k}\ve{R}} \begin{pmatrix} \hat{a}_{\ve{k},\sigma}^\dag \\ \hat{b}_{\ve{k},\sigma}^\dag \end{pmatrix} \\ s^z_{\ve{R}} &= \frac{1}{N}\sum_{\ve{q}} s^z_{\ve{q}} e^{i\ve{q}\ve{R}}\end{split}\]

with this definition, both models take the form

(2.7)\[ \mc{H} = -\sum_{\sigma=1}^{N_\sigma} \sum_{\ve{k}} \hat{a}_{\ve{k},\sigma}^\dag \Bigg( \hat{b}_{\ve{k},\sigma} Z_0(\ve{k}) + \frac{\xi}{N}\sum_{\ve{q}} \hat{b}_{\ve{k}-\ve{q},\sigma} s^z_{\ve{q}} Z_\text{Yuk}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising}\]

with

\[\begin{split} Z_0(\ve{k}) &= \left\{ \begin{array}{ll} 2 t \left( e^{i\frac{\pi}{4}} \cos k_x + e^{-i\frac{\pi}{4}} \cos k_y \right) e^{-i\ve{k}_y} & C_{2v} \text{ model} \\ 2 t \left( e^{i\frac{\pi}{4}} \cos k_x + e^{-i\frac{\pi}{4}} \cos k_y \right) & C_{4v} \text{ model} \end{array} \right. \\ Z_\text{Yuk}(\ve{k}) &= \left\{ \begin{array}{ll} i 2 \left( -e^{i\frac{\pi}{4}} \sin k_x + e^{-i\frac{\pi}{4}} \sin k_y \right) e^{-i\ve{k}_y} & C_{2v} \text{ model} \\ i & C_{4v} \text{ model.} \end{array} \right.\end{split}\]

2.2.2. Symmetries

2.2.2.1. The \(C_{2v}\) model

The First model has a \(C_{2v}\) symmetry, consisting of two reflections: \(T_+\) and \(T_-\) on \(\ve{e}_\pm = \ve{e}_x \pm \ve{e}_y\), the \(\pi\) rotation needed by the point group can be obtained as \(T_\pi = T_+ \cdot T_-\). The \(T_-\) invariance hinges on the Z\(_2\) Ising symmetry, \(s^z \rightarrow -s^z\), and is therefore broken in the ordered phase.

The \(C_{2v}\) symmetry pins the Dirac points (up to a gauge choice) to \((\pi/2, \pm \pi/2)\), while in the ordered phase, meandering parallel to \(\ve{e}_+\) is possible.

To show this symmetry, we expand the momentum- and real-space vectors as: \(\ve{k} = k_+\ve{e}_+ + k_-\ve{e}_-\) and \(\ve{R} = R_+\ve{e}_+ + R_-\ve{e}_-\).

The first reflection \(T_+\) reads:

(2.8)\[\begin{split}\hat{T}_+^{-1} \begin{pmatrix} \hat{a}^\dag_{\ve{k},\sigma} \\ \hat{b}^\dag_{\ve{k},\sigma} \end{pmatrix} \hat{T}_+ &= \begin{pmatrix} \hat{b}^\dag_{(k_+,-k_-),\sigma} \phantom{\^{-i k_+}} \\ \hat{a}^\dag_{(k_+,-k_-),\sigma} e^{-i k_+} \end{pmatrix} \\ \hat{T}_+^{-1} s^z_{\ve{q}} \hat{T}_+ &= s^z_{(q_+,-q_-)}\end{split}\]

Inserting the above in Eq. (2.7), we obtain:

\[\begin{align*} &\hat{T}_+^{-1} \mc{H}^{C_{2v}} \hat{T}_+ = \\ &= -\sum_{\ve{k}} e^{i k_+} \hat{b}_{(k_+,-k_-),\sigma}^\dag \Bigg( \hat{a}_{(k_+,-k_-),\sigma} Z_0^{C_{2v}}(\ve{k}) + \xi\sum_{\ve{q}} \hat{a}_{(k_+-q_+,-k_-+q_-),\sigma} \hat{s}^z_{(q_+,-q_-)} Z_\text{Yuk}^{C_{2v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &= -\sum_{\ve{k}} e^{-i k_+} \hat{a}_{(k_+,k_-),\sigma}^\dag \Bigg( \hat{b}_{(k_+,k_-),\sigma} \bar{Z}_0^{C_{2v}}(k_+,-k_-) + \xi\sum_{\ve{q}} \hat{b}_{(k_+-q_+,k_--q_-),\sigma} \hat{s}^z_{(q_+,q_-)} \bar{Z}_\text{Yuk}^{C_{2v}}(k_+,-k_-) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ & \qquad \left[ \begin{array}{rl} \bar{Z}_0^{C_{2v}}(k_+,-k_-) &= Z_0^{C_{2v}}(\ve{k}) e^{i k_+} \\ \bar{Z}_\text{Yuk}^{C_{2v}}(k_+,-k_-) &= Z_\text{Yuk}^{C_{2v}}(\ve{k}) e^{i k_+} \end{array} \right] \\ &= -\sum_{\ve{k}} \hat{a}_{(k_+,k_-),\sigma}^\dag \Bigg( \hat{b}_{(k_+,k_-),\sigma} Z_0^{C_{2v}}(\ve{k}) + \xi\sum_{\ve{q}} \hat{b}_{\ve{k}-\ve{q},\sigma} \hat{s}^z_{\ve{q}} Z_\text{Yuk}^{C_{2v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &=\mc{H}^{C_{2v}} \end{align*}\]

In real space, \(\hat{T}_+\) translates to:

\[\begin{split}\hat{T}_+^{-1} \begin{pmatrix} \hat{a}^\dag_{\ve{R},\sigma} \\ \hat{b}^\dag_{\ve{R},\sigma} \end{pmatrix} \hat{T}_+ &= \begin{pmatrix} \hat{b}^\dag_{(R_+ ,-R_-),\sigma} \\ \hat{a}^\dag_{(R_++1,-R_-),\sigma} \end{pmatrix} \\ \hat{T}_+^{-1} s^z_{\ve{R}} \hat{T}_+ &= s^z_{(R_+,-R_-)}\end{split}\]

The second reflection \(T_-\) can be expressed as:

(2.9)\[\begin{split}\hat{T}_-^{-1} \begin{pmatrix} \hat{a}^\dag_{\ve{k},\sigma} \\ \hat{b}^\dag_{\ve{k},\sigma} \end{pmatrix} \hat{T}_- &= \begin{pmatrix} \hat{b}^\dag_{(-k_+,k_-),\sigma} \phantom{e^{i k_-}}\\ \hat{a}^\dag_{(-k_+,k_-),\sigma} e^{i k_-} \end{pmatrix} \\ \hat{T}_-^{-1} s^z_{\ve{q}} \hat{T}_- &= -s^z_{(-q_+,q_-)}\end{split}\]

Inserting this into Eq. (2.7), we obtain:

\[\begin{align*} &\hat{T}_-^{-1} \mc{H}^{C_{2v}} \hat{T}_- = \\ &= -\sum_{\ve{k}} e^{-i k_-} \hat{b}_{(-k_+,k_-),\sigma}^\dag \Bigg( \hat{a}_{(-k_+,k_-),\sigma} Z_0(\ve{k}) + \xi\sum_{\ve{q}} \hat{a}_{(-k_++q_+,k_--q_-),\sigma} (-\hat{s}^z_{(-q_+,q_-)}) Z_\text{I}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &= -\sum_{\ve{k}} e^{i k_-} \hat{a}_{(k_+,k_-),\sigma}^\dag \Bigg( \hat{b}_{(k_+,k_-),\sigma} \bar{Z}_0(-k_+,k_-) + \xi\sum_{\ve{q}} \hat{b}_{(k_+-q_+,k_--q_-),\sigma} (-\hat{s}^z_{(q_+,q_-)}) \bar{Z}_\text{I}(-k_+,k_-) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ & \qquad \left[ \begin{array}{rl} \bar{Z}_0(-k_+,k_-) &= Z_0^{C_{2v}}(\ve{k}) e^{-i k_-} \\ \bar{Z}_\text{I}(-k_+,k_-) &= -Z_\text{Yuk}^{C_{2v}}(\ve{k}) e^{-i k_-} \end{array} \right] \\ &= -\sum_{\ve{k}} \hat{a}_{\ve{k},\sigma}^\dag \Bigg( \hat{b}_{\ve{k},\sigma} Z_0^{C_{2v}}(\ve{k}) + \xi\sum_{\ve{q}} \hat{b}_{\ve{k}-\ve{q}} \hat{s}^z_{\ve{q}} Z_\text{Yuk}^{C_{2v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &=\mc{H}^{C_{2v}} \end{align*}\]

In real space, \(\hat{T}_-\) translates to:

\[\begin{split}\hat{T}_-^{-1} \begin{pmatrix} \hat{a}^\dag_{\ve{R},\sigma} \\ \hat{b}^\dag_{\ve{R},\sigma} \end{pmatrix} \hat{T}_- &= \begin{pmatrix} \hat{b}^\dag_{(-R_+, R_-),\sigma} \\ \hat{a}^\dag_{(-R_+, R_--1),\sigma} \end{pmatrix} \\ \hat{T}_-^{-1} s^z_{\ve{R}} \hat{T}_- &= s^z_{(-R_+,R_-)}\end{split}\]

2.2.2.2. The \(C_{4v}\) model

The Second model has a \(C_{4v}\) symmetry, consisting of a rotation by \(\tfrac{\pi}{2}\) and reflections on the x- and y-axis.

The corresponding operators in momentum space are:

\[\begin{align*} \hat{T}_{\pi/2}^{-1} \begin{pmatrix} \hat{a}_{\ve{k},\sigma}^{\dag} \\ \hat{b}_{\ve{k},\sigma}^{\dag} \end{pmatrix} \hat{T}_{\pi/2} &= \begin{pmatrix} \hat{b}_{(-k_y,k_x),\sigma}^{\dag} \\ \hat{a}_{(-k_y,k_x),\sigma}^{\dag} \end{pmatrix} & % \hat{T}_x^{-1} \begin{pmatrix} \hat{a}_{\ve{k},\sigma}^{\dag} \\ \hat{b}_{\ve{k},\sigma}^{\dag} \end{pmatrix} \hat{T}_x &= \begin{pmatrix} \hat{a}_{(k_x,-k_y),\sigma}^{\dag} \\ \hat{b}_{(k_x,-k_y),\sigma}^{\dag} \end{pmatrix} & % \hat{T}_{y}^{-1} \begin{pmatrix} \hat{a}_{\ve{k},\sigma}^{\dag} \\ \hat{b}_{\ve{k},\sigma}^{\dag} \end{pmatrix} \hat{T}_{y} &= \begin{pmatrix} \hat{a}_{(-k_x,k_y),\sigma}^{\dag} \\ \hat{b}_{(-k_x,k_y),\sigma}^{\dag} \end{pmatrix} \\ % % \hat{T}_{\pi/2}^{-1} \hat{s}^{z}_{\ve{q}} \hat{T}_{\pi/2} &= -\hat{s}^{z}_{ (-q_y,q_x)} & \hat{T}_{x}^{-1} \hat{s}^{z}_{\ve{q}} \hat{T}_{x} &= \hat{s}^{z}_{ (q_x,-q_y)} & \hat{T}_{y}^{-1} \hat{s}^{z}_{\ve{q}} \hat{T}_{y} &= \hat{s}^{z}_{ (-q_x,q_y)} \end{align*}\]

And in real space:

\[\begin{align*} \hat{T}_{\pi/2}^{-1} \begin{pmatrix} \hat{a}_{\ve{R},\sigma}^{\dag} \\ \hat{b}_{\ve{R},\sigma}^{\dag} \end{pmatrix} \hat{T}_{\pi/2} &= \begin{pmatrix} \hat{b}_{(-R_y,R_x),\sigma}^{\dag} \\ \hat{a}_{(-R_y,R_x),\sigma}^{\dag} \end{pmatrix} & % \hat{T}_x^{-1} \begin{pmatrix} \hat{a}_{\ve{R},\sigma}^{\dag} \\ \hat{b}_{\ve{R},\sigma}^{\dag} \end{pmatrix} \hat{T}_x &= \begin{pmatrix} \hat{a}_{(R_x,-R_y),\sigma}^{\dag} \\ \hat{b}_{(R_x,-R_y),\sigma}^{\dag} \end{pmatrix} & % \hat{T}_{y}^{-1} \begin{pmatrix} \hat{a}_{\ve{R},\sigma}^{\dag} \\ \hat{b}_{\ve{R},\sigma}^{\dag} \end{pmatrix} \hat{T}_{y} &= \begin{pmatrix} \hat{a}_{(-R_x,R_y),\sigma}^{\dag} \\ \hat{b}_{(-R_x,R_y),\sigma}^{\dag} \end{pmatrix} \\ \hat{T}_{\pi/2}^{-1} \hat{s}^{z}_{\ve{R}} \hat{T}_{\pi/2} &= -\hat{s}^{z}_{ (-R_y,R_x)} & \hat{T}_{x}^{-1} \hat{s}^{z}_{\ve{R}} \hat{T}_{x} &= \hat{s}^{z}_{ (R_x,-R_y)} & \hat{T}_{y}^{-1} \hat{s}^{z}_{\ve{R}} \hat{T}_{y} &= \hat{s}^{z}_{ (-R_x,R_y)} \end{align*}\]

In the Ising ordered phase, the Ising symmetry \(\hat{s}^{z}_{\ve{R}} \rightarrow -\hat{s}^{z}_{\ve{R}}\) is broken, which reduces \(\hat{T}_{\pi/2}\) to \(\hat{T}_{\pi}\), such that the \(C_{4v}\) symmetry is reduced to \(C_{2v}\). This reduced symmetry allows the cones to meander.

Particle-hole symmetry:

The particle-hole symmetry \(\hat{T}_{\text{ph}}\) implies that energy eigenstates satisfy \(E(\ve{k}) = - E(-\ve{k} + \ve{Q})\), \(\ve{Q} = (\pi, \pi)\).

(2.10)\[\begin{split}\hat{T}_\text{ph}^{-1} \alpha \begin{pmatrix} \hat{a}_{\ve{R},\sigma}^\dag \\ \hat{b}_{\ve{R},\sigma}^\dag \end{pmatrix} \hat{T}_\text{ph} &= \bar{\alpha} (-1)^{R_x + R_y} \begin{pmatrix} \hat{a}_{\ve{R},\sigma} \\ \hat{b}_{\ve{R},\sigma} \end{pmatrix} \\ \hat{T}_\text{ph}^{-1} \alpha s^z_{\ve{R}} \hat{T}_\text{ph} &= - \bar{\alpha} s^z_{\ve{R}}\end{split}\]

(2.11)\[\begin{split}\hat{T}_\text{ph}^{-1} \alpha \begin{pmatrix} \hat{a}^\dag_{\ve{k},\sigma} \\ \hat{b}^\dag_{\ve{k},\sigma} \end{pmatrix} \hat{T}_\text{ph} &= \bar{\alpha} \begin{pmatrix} \hat{a}_{\ve{Q}-\ve{k},\sigma} \\ \hat{b}_{\ve{Q}-\ve{k},\sigma} \end{pmatrix} \quad \ve{Q} = \begin{pmatrix} \pi \\ \pi \end{pmatrix} \\ \hat{T}_\text{ph}^{-1} \alpha s^z_{\ve{q}} \hat{T}_\text{ph} &= -\bar{\alpha} s^z_{-\ve{q}}\end{split}\]

Inserting this in Eq. (2.7), we obtain:

\[\begin{align*} \hat{T}_\text{ph}^{-1} \mc{H}^{C_{4v}} \hat{T}_\text{ph} &= -\sum_{\ve{k}} \hat{a}_{\ve{Q}-\ve{k}} \Bigg( \hat{b}^\dag_{\ve{Q}-\ve{k}} \bar{Z}_0^{C_{4v}}(\ve{k}) +\frac{\xi}{N} \sum_{\ve{q}} \hat{b}^\dag_{\ve{Q}-\ve{k}+\ve{q}} \left( -s^z_{-\ve{q}} \right) \bar{Z}_\text{Yuk}^{C_{4v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &= \phantom{-} \sum_{\ve{k}} \hat{a}^\dag_{\ve{Q}-\ve{k}} \Bigg( \hat{b}_{\ve{Q}-\ve{k}} Z_0^{C_{4v}}(\ve{k}) +\frac{\xi}{N} \sum_{\ve{q}} \hat{b}_{\ve{Q}-\ve{k}+\ve{q}} \left( -s^z_{-\ve{q}} \right) Z_\text{Yuk}^{C_{4v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &= \phantom{-} \sum_{\ve{k}} \hat{a}^\dag_{\ve{k}} \Bigg( \hat{b}_{\ve{k}} Z_0^{C_{4v}}(\ve{Q}-\ve{k}) +\frac{\xi}{N} \sum_{\ve{q}} \hat{b}_{\ve{k}-\ve{q}} \left( -s^z_{\ve{q}} \right) Z_\text{Yuk}^{C_{4v}}(\ve{Q}-\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ & \qquad \left[ \begin{array}{rl} Z_0^{C_{4v}}(\ve{Q}-\ve{k}) &= -Z_0^{C_{4v}}(\ve{k}) \\ Z_\text{Yuk}^{C_{4v}}(\ve{Q}-\ve{k}) &= Z_\text{Yuk}^{C_{4v}}(\ve{k}) \end{array} \right] \\ &= -\sum_{\ve{k}} \hat{a}^\dag_{\ve{k}} \Bigg( \hat{b}_{\ve{k}} Z_0^{C_{4v}}(\ve{k}) +\frac{1}{N} \sum_{\ve{q}} \hat{b}_{\ve{k}-\ve{q}} \left( s^z_{\ve{q}} \right) Z_\text{Yuk}^{C_{4v}}(\ve{k}) \Bigg) + \mathrm{H.c.} + H_\text{Ising} \\ &= \mc{H}^{C_{4v}} \end{align*}\]

As a result of this symmetry, the single particle spectral function satisfies

\[ A(\ve{k}, \omega) = A(-\ve{k} + \ve{Q}, -\omega), \text{ with } \ve{Q} = (\pi, \pi ). \]

2.3. Lattice mean-field theory

The key point of both models is that the point group and particle-hole symmetries are tied to the flipping of the Ising spin degree of freedom. In the large-\(h\) limit, the ground state has the full symmetry of the model Hamiltonian and at the mean-field level, we can set \( \langle \hat{s}_{\boldsymbol{R}}^z \rangle = 0\). In this limit, the Dirac cones are pinned by symmetry. In the opposite small-\(h\) limit, the Ising spins order, i.e. \( \langle \hat{s}_{\boldsymbol{R}}^z \rangle \neq 0\). Thereby, the \(C_{2v}\) (\(C_{4v}\)) symmetry is reduced to \(\hat{T}_+\) (\(C_{2v}\)). At the mean-field level, this induces a meandering of the Dirac points in the BZ, see Fig. 2.1(a), and an anisotropy in the Fermi velocities. In this section, we present the mean-field calculation. At this level of approximation, the transition turns out to be continuous, in agreement with the large-\(N\) analysis [53].

We expand Eq. (2.7) around \(\langle \hat{s}^z_{\ve{R}} \rangle \equiv \phi\). The resulting mean-field Hamiltonian reads:

(2.12)\[\begin{align} \mc{H}_\text{MF} &= \sum_{\sigma=1}^{N_\sigma} \sum_{\ve{k}} \hat{K}_{\ve{k},\sigma}(\phi\xi) + \sum_{\ve{R}} \hat{I}_{\ve{R}}(\phi). \end{align}\]

With:

\[\begin{align*} \hat{K}_{\ve{k},\sigma}(\phi\xi) &= -\hat{a}_{\ve{k},\sigma}^\dag \hat{b}_{\ve{k},\sigma} Z(\ve{k}, \phi\xi) & Z(\ve{k}, \phi\xi) &= Z_0(\ve{k}) + \phi \xi Z_\text{Yuk}(\ve{k}) & \hat{I}_{\ve{R}}(\phi) &= -4 J \left( \phi\hat{s}^z_{\ve{R}} - \tfrac{1}{2}\phi^2 \right) - h\hat{s}^x_{\ve{R}}. \end{align*}\]

The fermionic dispersion is

(2.13)\[\pm \left| Z(\ve{k}, \phi\xi) \right|.\]

To determine the nature of the zero-temperature phase transition, we determine the order parameter \(\phi\) for a given transverse field \(h\) by minimizing the ground state energy \(E_{0, \text{MF}} = \lim_{\beta \rightarrow \infty} E_\text{MF} = \lim_{\beta \rightarrow \infty} \Braket{\mc{H}_\text{MF}}_\text{MF}.\)

\[\begin{align*} E_\text{MF} &= \Braket{\mc{H}_\text{MF}}_\text{MF} = \frac {{\rm Tr}\big(\exp(-\beta \mc{H}_\text{MF}) \mc{H}_\text{MF} \big)} {{\rm Tr} \big(\exp(-\beta \mc{H}_\text{MF})\big)} \\ &= N_\sigma \sum_{\ve{k}} |Z(\ve{k}, \phi\xi)| \frac {1 - \exp\left(\beta |Z(\ve{k}, \phi\xi)| \right)} {1 + \exp\left(\beta |Z(\ve{k}, \phi\xi)| \right)} + L^2 \left( 2\phi^2 - \sqrt{h^2 + 16\phi^2} \tanh\left(\beta\sqrt{h^2 + 16\phi^2}\right) \right) \end{align*}\]

(2.14)\[\begin{split}\lim_{\beta \rightarrow \infty} \frac{E_\text{MF}}{L^2} &= \underbrace{- \frac{N_\sigma}{L^2} \sum_{\ve{k}} |Z(\ve{k}, \phi\xi)|}_{\epsilon_\text{F}(L, \phi \xi)} + \underbrace{\left( 2\phi^2 - \sqrt{h^2 + 16\phi^2} \right)}_{\epsilon_\text{I}(\phi, h)} \\ &\partial_{\phi} \epsilon_\text{F}(L, \phi) = -\frac{N_\sigma}{L^2} \sum_{\ve{k}} \frac{1}{|Z(\ve{k}, \phi)|} \Big( \Re\left(Z_0(\ve{k}) \bar{Z}_\text{Yuk}(\ve{k})\right) + \phi |Z_\text{Yuk}(\ve{k})|^2 \Big) \\ &\partial^2_{\phi} \epsilon_\text{F}(L, \phi) = -\frac{N_\sigma}{L^2} \sum_{\ve{k}} \left( \frac{|Z_\text{Yuk}(\ve{k})|^2}{|Z(\ve{k}, \phi)|} - \frac{\Re\left(Z_0(\ve{k}) \bar{Z}_\text{Yuk}(\ve{k})\right)} {|Z(\ve{k}, \phi)|^3} \Big( \Re\left(Z_0(\ve{k}) \bar{Z}_\text{Yuk}(\ve{k})\right) + \phi |Z_\text{Yuk}(\ve{k})| \Big) \right)\end{split}\]

Equation (2.14) separates into a fermionic and an Ising part, \(\epsilon_\text{F}\) and \(\epsilon_\text{I}\). While \(\epsilon_\text{I}\) has a well-behaved, closed form, \(\epsilon_\text{F}\) has some non-analytic points for finite lattices (cf. Fig. 2.4). Namely, \(\partial^2_\phi \epsilon_\text{F}\) diverges, if \(Z(\ve{k}, \phi\xi)\) vanishes.

This corresponds to a finite-size artifact which can be qualitatively understood with the help of Fig. 2.3. Essentially, \(\phi\) shifts the single particle energy and produces level crossing reminiscent of those produced when twisting boundary conditions [58, 59]. Fig. 2.3 shows the valence band of a one-dimensional Dirac cone on a lattice with five \(\ve{k}\) points at two different twists. As a function of the twist the \(\ve{k}\)-point will cross the Fermi surface and at this crossing point, a singularity in the kinetic energy – corresponding to a level crossing – will appear. This is explicitly shown in Fig. 2.5. This observation means that the thermodynamic limit and the derivative \(\partial_{\phi} \) do not commute: one should first take the thermodynamic limit prior to carrying out the derivative.

To avoid this artifact, we consider two different approaches:

1. Choose a weaker coupling \(\xi\), such that the Dirac points do not cross the Fermi surface in proximity to the critical point (see Fig. 2.5 vs Fig. 2.6). However, choosing a small \(\xi\) may result in a slow flow from the 3d Ising fixed point of the unperturbed Ising model to nematic criticality.

2. Choose antiperiodic boundary conditions in space for the fermions, so as to shift the \(\ve{k}\) points away from the Fermi surface: Fig. 2.7. However, this choice results in large finite size effects presumably due to the boundary-condition-induced finite size gap.

It turns out the first option is the best choice and that \(\xi\) can be chosen large enough so as to minimize the aforementioned crossover effects.

The \(C_{4v}\) model (Fig. 2.8) shows different behaviors between systems with linear sizes \(L = 2 + 4\mathbb{N}\) and \(L = 4\mathbb{N}\). At \(L = 4\mathbb{N}\) with periodic boundary conditions, the Dirac points in the disordered phase are located at \(\ve{k}\)-points resolved by the finite lattice. This is not the case at \(L = 2 + 4\mathbb{N}\) (cf. Fig. 2.9). As a result, the \(L = 4\mathbb{N}\) sizes have a smoothed-out phase transition. Nevertheless, both \(L = 4\mathbb{N}\) and \(L = 2 + 4\mathbb{N}\) converge to the same result in the thermodynamic limit. The Monte-Carlo simulations also have odd-even effects, as elaborated in Section 2.8.3. It turns out that even system sizes produce nicer numerical results for the phase transition.

Fig. 2.3 Sketch for understanding finite-size artifacts of Dirac systems as a function of \(\ve{k}\)-quantization. Shown is the valence band of a one-dimensional Dirac cone on a lattice of size \(L=5\). We can see that the choice of the momenta quantization (i.e. boundary conditions) affects the energy. Left: Dirac point belongs to the set of finite-size \(\ve{k}\)-points. Right: Dirac point is between two finite-size \(\ve{k}\)-points. In the quantizations shown on the left, the ground state energy is higher, and a level crossing occurs when a \(\ve{k}\) point crosses the Fermi surface. In nematic transitions, the translation symmetry is not broken such that momenta are well defined and the \(\ve{k}\) quantization for a given lattice size remains unchanged. However, the position of the Dirac point meanders. The energy level crossing, that originates, is reminiscent of that obtained when twisting boundary conditions [58, 59].

Fig. 2.4 Fermionic part of Mean-field ground state energy. The second derivative with respect to the Ising field \(\phi\) diverges at level crossings as described in the main text.

Fig. 2.5 Mean-field results for the \(C_{2v}\) model with coupling strength \(\xi = 0.75\). There are many discontinuities in the order parameter \(\phi\) when tuning the transverse field \(h\) trough the phase transition. These finite size artifacts are disappearing at \(L \rightarrow \infty\).

Fig. 2.6 Mean-field results for the \(C_{2v}\) model with coupling strength \(\xi = 0.25\). Most discontinuities visible for \(\xi = 0.75\) (cf. Fig. 2.5) are avoided at this coupling strength.

Fig. 2.7 Brillouin zone of \(C_{2v}\) model with \(\ve{k}\) points of a \(8*8\) lattice. Left: With periodic boundary conditions. Right: With antiperiodic boundary conditions for movement parallel to \((1,-1)\). Also sketched: Dispersion in disordered phase and trajectory of Dirac cones.

Fig. 2.8 Mean-field results for the \(C_{4v}\) model with coupling strength \(\xi = 1\). There is an qualitative difference between system sizes \(L = 2 + 4\mathbb{N}\) and \(L = 4\mathbb{N}\): The former converges faster, but both are still converging to the same result for \(L \rightarrow \infty\).

Fig. 2.9 Brillouin zone of \(C_{4v}\) model with \(\ve{k}\) points of \(6*6\) and \(8*8\) lattice. Also sketched: Dispersion in disordered phase and trajectory of Dirac cones. Left: \(6*6\) lattice, the Dirac cones in the disordered phase are each centered between four \(\ve{k}\) points. Right: \(8*8\) lattice, the Dirac cones in the disordered phase are directly resolved by the \(\ve{k}\) points.

To summarize, the mean-field calculation finds a continuous transitions for both models and to avoid finite size artifacts in the QMC simulations, \(\xi\) should not be chosen too large.

2.4. Continuum field theory

In order to investigate whether the above remains true upon the inclusion of order-parameter fluctuations, we derive corresponding continuum field theories, which are amenable to RG analyses.

We derive the low energy model from Eq. (2.7) by expansion in \(\ve{k}\) around the nodal points \(\ve{K}_i\) and for a scalar Ising field \(\phi(\ve{q})\):

\[\begin{align*} \mc{H} &= -\sum_{\sigma=1}^{N_\sigma} \sum_{i} \int\!\text{d}\ve{\kappa}\; \hat{a}_{\ve{K}_i+\ve{\kappa},\sigma}^\dag \Bigg( \hat{b}_{\ve{K}_i+\ve{\kappa},\sigma} Z_0(\ve{K}_i+\ve{\kappa}) \nonumber\\ & + \frac{1}{2\pi}\int\!\text{d}\ve{q}\; \hat{b}_{\ve{K}_i+\ve{\kappa}-\ve{q},\sigma} \phi({\ve{q}}) Z_\text{Yuk}(\ve{K}_i+\ve{\kappa}) \Bigg) + \mathrm{H.c.} + H_\text{Ising}\{\phi\} \end{align*}\]

In leading order in \( \ve{\kappa} \) we obtain:

\[\begin{align*} \mc{H} &= -\sum_{\sigma=1}^{N_\sigma} \sum_{i} \int\!\text{d}\ve{\kappa}\; \hat{a}_{\ve{K}_i+\ve{\kappa},\sigma}^\dag \Bigg( \hat{b}_{\ve{K}_i+\ve{\kappa},\sigma} \ve{\kappa} \nabla Z_0(\ve{K}_i) \nonumber\\ & + \frac{1}{2\pi}\int\!\text{d}\ve{q}\; \hat{b}_{\ve{K}_i+\ve{\kappa}-\ve{q},\sigma} \phi({\ve{q}}) Z_\text{Yuk}(\ve{K}_i) \Bigg) + \mathrm{H.c.} + H_\text{Ising}\{\phi\}. \end{align*}\]

Introducing the Fourier transformations:

\[\begin{align*} \begin{pmatrix} \hat{a}^i_{\sigma}(\ve{r}) \\ \hat{b}^i_{\sigma}(\ve{r}) \end{pmatrix} &= \frac{1}{2\pi} \int\!\text{d}\ve{\kappa}\; e^{i\ve{\kappa} \ve{r}} \begin{pmatrix} \hat{a}_{\ve{K}_i+\ve{\kappa},\sigma} \\ \hat{b}_{\ve{K}_i+\ve{\kappa},\sigma} \end{pmatrix} \\ \phi(\ve{r}) &= \frac{1}{2\pi} \int\!\text{d}\ve{\kappa}\; e^{i\ve{q} \ve{r}} \phi(\ve{q}). \end{align*}\]

and defining:

\[\begin{align*} \ve{v}_i &\equiv \nabla Z_0(\ve{K}_i) & I_i &\equiv Z_\text{Yuk}(\ve{K}_i) \end{align*}\]

The Hamiltonian takes the form:

(2.15)\[\mc{H} = -\sum_{\sigma=1}^{N_\sigma} \sum_{i} \int\!\text{d}\ve{r}\; {\hat{a}^i_{\sigma}}{} (\ve{r})^\dag \Big( i\ve{v}_i \cdot \nabla_{\ve{r}} + \phi(\ve{r}) I_i \Big) \hat{b}^i_{\sigma}(\ve{r}) + \mathrm{H.c.} + H_\text{Ising}\{\phi\}.\]

2.4.1. The \(C_{2v}\) model

The \(C_{2v}\) model has the nodal points \(\ve{K}_\pm = \begin{pmatrix} \pi/2 \\ \pm\pi/2 \end{pmatrix}\). By defining the four-component Dirac spinor

\[ \Psi_\sigma(\ve{r}) = \begin{pmatrix} \hat{a}^+_{\sigma}(\ve{r}) & \hat{b}^+_{\sigma}(\ve{r}) & \hat{a}^-_{\sigma}(\ve{r}) & \hat{b}^-_{\sigma}(\ve{r}) \end{pmatrix}^{\rm T}, \]

where \(\hat a_{\sigma}^{\pm}\) and \(\hat b_{\sigma}^{\pm}\) corresponds to hole excitations near \(\boldsymbol{K}_\pm\) on the \(A\) and \(B\) sublattices, respectively, and

\[\begin{align*} \tau_1 &= \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & \tau_2 &= \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & \tau_3 &= \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}. \end{align*}\]

Eq. (2.15) can be written as

\[\begin{align*} \mc{H}^{C_{2v}} =& \sum_{\sigma=1}^{N_\sigma} \int\!\text{d}\ve{r}\; \Psi_\sigma^\dag(\ve{r}) \Bigg[ 2i t \begin{pmatrix} \tau_2 & 0 \\ 0 & -\tau_1 \end{pmatrix} \partial_{r_+} + 2i t \begin{pmatrix} \tau_1 & 0 \\ 0 & -\tau_2 \end{pmatrix} \partial_{r_-} + 2\sqrt{2}\xi \phi(\ve{r}) \begin{pmatrix} -\tau_2 & 0 \\ 0 & -\tau_1 \end{pmatrix} \Bigg]\Psi_\sigma(\ve{r}) \\ &+ H_\text{Ising}\{\phi\}. \end{align*}\]

Introducing the gamma matrices, that realize a four-dimensional representation of the Clifford algebra

\[\begin{align*} \gamma_0 &= \begin{pmatrix} -\tau_3 & 0 \\ 0 & -\tau_3 \end{pmatrix} & \gamma_1 &= \begin{pmatrix} \tau_1 & 0 \\ 0 & \tau_2 \end{pmatrix} & \gamma_2 &= \begin{pmatrix} -\tau_2 & 0 \\ 0 & -\tau_1 \end{pmatrix} & \{\gamma_\alpha, \gamma_\beta\} &= 2 \delta_{\alpha\beta}, \end{align*}\]

we can write the action in the form

(2.16)\[\begin{align} S^{C_{2v}} = &\int\!\text{d}^D x\, \sum_{\sigma=1}^{N_\sigma} \Big[ \Psi_{\sigma}^\dag(x) \big[ \mathbb{1} \partial_\tau +v \gamma_0\gamma_1 \partial_+ +v \gamma_0\gamma_2 \partial_- +g \phi(x) \gamma_2 \big] \Psi_\sigma(x) \Big]+ S_\text{Ising}(\{ \phi \}). \end{align}\]

2.4.2. The \(C_{4v}\) model

The \(C_{4v}\) model has the nodal points \(\ve{K}_{+\pm} = \begin{pmatrix} \pi/2 \\ \pm\pi/2 \end{pmatrix}\) and \(\ve{K}_{-\pm} = -\ve{K}_\pm\). By defining the eight-component Dirac spinor

\[ \Psi_\sigma(\ve{r}) = \begin{pmatrix} \hat{a}^{++}_{\sigma}(\ve{r}) & \hat{b}^{++}_{\sigma}(\ve{r}) & \hat{a}^{-+}_{\sigma}(\ve{r}) & \hat{b}^{-+}_{\sigma}(\ve{r}) & \hat{a}^{+-}_{\sigma}(\ve{r}) & \hat{b}^{+-}_{\sigma}(\ve{r}) & \hat{a}^{--}_{\sigma}(\ve{r}) & \hat{b}^{--}_{\sigma}(\ve{r}) \end{pmatrix}^{\rm T}, \]

where \(\hat a_{\sigma}^{+\pm}\) and \(\hat b_{\sigma}^{+\pm}\) (\(\hat a_{\sigma}^{-\pm}\) and \(\hat b_{\sigma}^{-\pm}\)) correspond to hole excitations near \(\boldsymbol{K}_\pm\) (\(-\boldsymbol{K}_\pm\)). Eq. (2.15) can be written as

\[\begin{align*} \mc{H}^{C_{4v}} &= \sum_{\sigma=1}^{N_\sigma} \int\!\text{d}\ve{r}\; \Psi_{\sigma}^\dag(\ve{r}) \Bigg[ \\ & 2i t \Bigg( \begin{pmatrix} \tau_1 & 0 \\ 0 & -\tau_1 \end{pmatrix} \oplus \begin{pmatrix} -\tau_2 & 0 \\ 0 & \tau_2 \end{pmatrix} \Bigg) \partial_{r_+} + 2i t \Bigg( \begin{pmatrix} -\tau_2 & 0 \\ 0 & \tau_2 \end{pmatrix} \oplus \begin{pmatrix} \tau_1 & 0 \\ 0 & -\tau_1 \end{pmatrix} \Bigg) \partial_{r_-} \\ &+ 2\sqrt{2}\xi \phi(\ve{r}) \Bigg( \begin{pmatrix} \tau_2 & 0 \\ 0 & \tau_2 \end{pmatrix} \oplus \begin{pmatrix} \tau_2 & 0 \\ 0 & \tau_2 \end{pmatrix} \Bigg) \\ &\qquad\Bigg] \Psi_{\sigma}(\ve{r}) + H_\text{Ising}\{\phi\}, \end{align*}\]

where \(\oplus\) denotes the matrix direct sum. Introducing another set of gamma matrices, that also realize a four-dimensional representation of the Clifford algebra

\[\begin{align*} \tilde{\gamma}_0 &= \begin{pmatrix} \tau_3 & 0 \\ 0 & -\tau_3 \end{pmatrix} & \tilde{\gamma}_1 &= \begin{pmatrix} \tau_1 & 0 \\ 0 & \tau_1 \end{pmatrix} & \tilde{\gamma}_2 &= \begin{pmatrix} \tau_2 & 0 \\ 0 & \tau_2 \end{pmatrix} & \{\tilde{\gamma}_\alpha, \tilde{\gamma}_\beta\} &= 2 \delta_{\alpha\beta}, \end{align*}\]

we can write the action in the compact form

(2.17)\[\begin{align} \begin{split} S^{C_{4v}} = &\int\!\text{d}^D x\, \sum_{\sigma=1}^{N_\sigma} \Psi_{\sigma}^\dag(x) \Big[ \mathbb{1} \partial_\tau +v \left( \tilde{\gamma}_0\tilde{\gamma}_1 \oplus \tilde{\gamma}_0\tilde{\gamma}_2 \right) \partial_+ +v \left( \tilde{\gamma}_0\tilde{\gamma}_2 \oplus \tilde{\gamma}_0\tilde{\gamma}_1 \right) \partial_- +g \phi(x) \left( \tilde{\gamma}_2 \oplus \tilde{\gamma}_2\right) \Big] \Psi_{\sigma}(x) \\ &+ S_\text{Ising}(\{ \phi \}). \end{split} \end{align}\]

2.4.3. Finalized field theory

So, all in all, to leading order in the gradient expansion around the nodal points, we obtain the Euclidean action \(S = \int \du^2x \du \tau (\mathcal L_\Psi + \mathcal L_\phi)\) with

(2.18)\[ \mathcal{L}_{\Psi}^{C_{2v}} = \Psi_{\sigma}^\dagger \left(\partial_\tau + \gamma_0\gamma_1 v_\parallel \partial_+ + \gamma_0\gamma_2 v_\perp \partial_- + g \phi \gamma_2 \right) \Psi_\sigma. \]

and

(2.19)\[ \mathcal{L}_\Psi^{C_{4v}} = \Psi_{\sigma}^\dagger \bigl[\partial_\tau + \tilde{\gamma}_0 (\tilde{\gamma}_1v_\parallel \oplus \tilde{\gamma}_2v_\perp) \partial_+ + \tilde{\gamma}_0(\tilde{\gamma}_2v_\perp \oplus \tilde{\gamma}_1 v_\parallel) \partial_- + g \phi (\tilde{\gamma}_2 \oplus \tilde{\gamma}_2)\bigr] \Psi_{\sigma} \]

In the above Lagrangians, we have assumed the summation convention over repeated indices. Furthermore, we allow for anisotropic Fermi velocities \(v_\parallel\) and \(v_\perp\), corresponding to the directions parallel and perpendicular to the shift of the Dirac cones in the ordered phase, with \(v_\parallel = v_\perp \sim t\) at the UV cutoff scale \(\Lambda\). \(\partial_\pm\) denotes the spatial derivative in the direction along \(\boldsymbol{K}_\pm\). The fermions couple via \(g \sim \xi\) to the Ising order-parameter field \(\phi\). The dynamics of the Ising field is governed by the usual \(\phi^4\) Lagrangian,

\[ \mathcal L_\phi = \frac{1}{2} \phi(r - \partial_\tau^2 - c_+^2 \partial_+^2 - c_-^2 \partial_-^2)\phi + \lambda \phi^4, \]

with the tuning parameter \(r\), the boson velocities \(c_\pm\), and the bosonic self-interaction \(\lambda\).

2.5. \(\epsilon\) expansion

The presence of a unique upper critical spatial dimension of three allows an \(\epsilon = 3-d\) expansion, with \(\epsilon = 1\) corresponding to the physical case. Because of the lack of Lorentz and continuous spatial rotational symmetries in the low-energy models, it is useful to employ a regularization in the frequency only, which allows us to rescale the different momentum components independently, and evaluate the loop integrals analytically (See Appendix A.1 for details). Two central properties of nematic quantum phase transitions in Dirac systems are revealed by the one-loop RG analysis: First, both models admit a stable fixed point featuring anisotropic power laws of the fermion and order parameter correlation functions.

In the \(C_{2v}\) model, both components of the Fermi velocity remain finite at the stable fixed point with \(0 < v_\parallel^* < v_\perp^*\). At the critical point, a unique timescale \(\tau\) emerges for both fields \(\Psi\) and \(\phi\) [60, 61], which scales with the two characteristic length scales \(\ell_+\) and \(\ell_-\) as \(\tau \sim \ell_+^{z_+} \sim \ell_-^{z_-}\), with associated dynamical critical exponents \(z_\pm = [1-\frac12 \eta_\phi + \frac12 \eta_\pm]^{-1}\) as \((z_+, z_-) = (1 + 0.3695\epsilon, 1 + 0.1086\epsilon) + \mathcal O(\epsilon^2)\), reflecting the absence of Lorentz and rotational symmetries at criticality.

By contrast, in the \(C_{4v}\) model, the fixed point is characterized by a maximal velocity anisotropy with \((v_\parallel^*,v_\perp^*) = (0,1)\) in units of fixed boson velocities \(c \equiv c_+ = c_- = 1\). This result is consistent with the large-\(N\) RG analysis in fixed \(d=2\) [52]. The fact that \(v_\parallel^*\) vanishes leads to the interesting behavior that the fixed-point couplings \(g^2_*\) and \(\lambda_*\) are bound to vanish in this case as well. This happens in a way that the ratio \((g^2/v_\parallel)_*\) remains finite, such that the boson anomalous dimensions become \(\eta_\phi = \eta_+ = \eta_- = \epsilon\). Importantly, as the fixed-point couplings \(g^2_*\) and \(\lambda^*\) vanish, we expect the one-loop result for the critical exponents to hold at all loop orders in the \(C_{4v}\) model. For the correlation-length exponent, we find \(1/\nu = 2 - \epsilon\). The remaining exponents can then be computed by assuming the usual hyperscaling relations [62]. The susceptibility exponent, for instance, becomes \(\gamma = 1\), independent of \(\epsilon\). This result is again consistent with the large-\(N\) calculation and has previously already been argued to hold exactly [52].

We note that the values of the exponents in the \(C_{4v}\) model are independent of the number of spinor components, in contrast to the situation in the \(C_{2v}\) model, as well as to the usual Gross-Neveu universality classes [43, 44, 45, 46, 57, 63]. The unique dynamical critical exponent in the \(C_{4v}\) model becomes \(z = 1\). We emphasize, however, that the critical point still does not feature emergent Lorentz symmetry [64] due to the anisotropic fermion spectral function. The second important property revealed by the RG analysis is that the stable fixed points in both models are approached only extremely slowly as function of RG scale, Fig. 2.10. This is universally true for the \(C_{4v}\) model, in which case \(v_\parallel\) corresponds to a marginally irrelevant parameter, hence scaling only logarithmically to zero while other irrelevant operators rapidly die out. This defines a quasiuniversal flow [65, 66] in which only the velocity anisotropy and not the initial ultraviolet values of other parameters determine the slow drift of the exponents. The RG suggests that this regime emerges at scales \(1/b \lesssim 0.05\) (cf. Appendix A.1), such that it will dominate numerical as well as experimental realizations of this critical phenomena. For a reasonable set of ultraviolet starting values and \(\epsilon = 1\), we find that the effective correlation-length exponent \(1/\nu_\text{eff}\) (anomalous dimension \(\eta_\phi^\text{eff}\)) approaches one from above (below), with sizable deviations at intermediate RG scales. Moreover, we also observe that the initial flows at high energy in the two models resemble each other, despite the fact that they substantially deviate from each other at low energy. This suggests that the flow is generically slow in the \(C_{2v}\) model as well.

Fig. 2.10 Ratio of Fermi velocities \(v_\perp/v_\parallel\) as function of RG scale \(1/b\) for both models. We assume ultraviolet initial values of \(v_\parallel(b=1) = v_\perp(b=1) = 0.25\), and set \(g^2/(v_\parallel v_\perp)(b=1)\) to the value at the respective stable fixed point. (a) Semilogarithmic, (b) linear plots. Starting at a temperature scale representative of the ultraviolet initial parameters, one has to cool the system by 2 orders of magnitude to start observing the differences between both models.

2.6. QMC setup

For the numerical simulations, we used the ALF program package [11, 12] that provides a general implementation of the finite-temperature auxiliary field QMC algorithm [8, 9, 10]. To formulate the path integral, we use a Trotter decomposition with time step \(\Delta_\tau t=0.1\) and choose a basis where \(\hat{s}_{\boldsymbol{R}}^{z} | s_{\boldsymbol{R}} \rangle = s_{\boldsymbol{R}} | s_{\boldsymbol{R}} \rangle\). The configuration space is that of a \((2+1)\)-dimensional Ising model and we use a single-spin-flip update to sample it. As we will show in Section 2.6.1 both models are negative-sign-problem free for all values of \(N_{\sigma}\) [67]. For our simulations, we have used an inverse temperature \(\beta = 4 L\) for \(L \times L\) lattices, and have checked that this choice of \(\beta\) reflects ground-state properties. For the results shown in the main text, we have fixed the parameters as \(J = t = 1\) and \(N_\sigma = 2\). In the \(C_{2v}\) model, we choose \(\xi = 0.25\), since larger values of \(\xi\) lead to spurious size effects that could falsely be interpreted as first-order transitions, as discussed in Section 2.3, see also Ref. [57] for a detailed discussion. In the \(C_{4v}\) model, we set \(\xi = 1\). As shown in Appendix A.5, other values of \(\xi\) and \(N_\sigma\) do not alter the continuous nature of the transition.

2.6.1. Absence of negative sign problem

Here we use the Majorana representation to demonstrate, using the results of Ref. [67], the absence of negative sign problem for all values of \(N_{\sigma}\). Both models have \(\text{SU}(N_{\sigma})\) symmetry. Since the Ising spins couple symmetrically to the fermion spins, \(\text{SU}(N_{\sigma})\) symmetry is present for all Ising spin configurations. Thereby, the fermion determinant of the \(\text{SU}(N_{\sigma})\) model corresponds to that of the \(\text{U}(1)\) model (\(N_{\sigma} =1\)) elevated to the power \(N_{\sigma}\). It hence suffices to demonstrate the absence of negative sign problem at \(N_{\sigma} =1 \). In this section, we will hence omit the spin index. Additionally we include a chemical potential term \(\mathcal{H}_{\mu}\), to show that there is also no sign problem under doping for even values of \(N_{\sigma}\).

2.6.1.1. The \(C_{2v}\) model

Consider the canonical transformation,

\[\begin{equation*} \begin{pmatrix} \hat{a}_{\ve{k}} \\ \hat{b}_{\ve{k}} \end{pmatrix} \rightarrow \begin{pmatrix} 1 & 0 \\ 0 & e^{-i 3 \pi / 4 } \end{pmatrix} \begin{pmatrix} \hat{a}_{\ve{k} - \Delta \ve{K}} \\ \hat{b}_{\ve{k} - \Delta \ve{K}} \end{pmatrix} \end{equation*}\]

with

\[\begin{equation*} \hat{b}_{\ve{R}} = \frac{1}{\sqrt{N}} \sum_{\ve{k} \in BZ } e^{i \ve{k} \cdot \ve{R}} \hat{b}_{\ve{k}} \end{equation*}\]

and \(\Delta \ve{K} = \frac{1}{4} \left( \ve{b}_- - \ve{b}_+ \right) \). Here, \(\ve{e}_i \cdot \ve{b}_j = 2 \pi \delta_{i,j}\). This canonical transformation renders the Hamiltonian real: the \(\pi\)-flux, is realized by changing the sign of the intra unit-cell hopping with respect to the other hoppings. More precisely after the transformation, the fermionic part of the Hamiltonian takes the form:

(2.20)\[\begin{split} &\mathcal{H}_0^{C_{2v}} = -t\sum_{\ve{R}} \hat{a}_{\ve{R}}^\dag \Big( - \hat{b}_{\ve{R}} + \hat{b}_{\ve{R}+\ve{e}_-} + \hat{b}_{\ve{R}+\ve{e}_--\ve{e}_+} + \hat{b}_{\ve{R}-\ve{e}_+} \Big) + \mathrm{H.c.} \\ &\mathcal{H}_\text{Yuk}^{C_{2v}} = -\xi \sum_{\ve{R}} s_{\ve{R}} \hat{a}^\dag_{\ve{R}} \Big( - \hat{b}_{\ve{R}} - \hat{b}_{\ve{R}+\ve{e}_-} - \hat{b}_{\ve{R}+\ve{e}_--\ve{e}_+} + \hat{b}_{\ve{R}-\ve{e}_+} \Big) + \mathrm{H.c.} \\ &\mathcal{H}_\mu^{C_{2v}} = \mu \sum_{\ve{R}} \Big( \hat{a}^\dag_{\ve{R}} \hat{a}_{\ve{R}} + \hat{b}^\dag_{\ve{R}} \hat{b}_{\ve{R}} \Big)\end{split}\]

In the above, we have considered an arbitrary set of Ising spins \(s_{\ve{R}} = \pm 1 \).

Since equation (2.20) is real, the corresponding fermion determinant for \(N_\sigma = 1\) is also real and therefore positive for even \(N_\sigma\).

For the sign to remain positive with odd \(N_\sigma\), we have to dismiss \(H_\mu\) and introduce Majorana fermions:

\[\begin{align*} \hat{a}^{\dagger}_{\ve{R}} &= \tfrac{1}{2} \left( \hat{\gamma}_{\ve{R},1} - i \hat{\ve{\gamma}}_{\ve{R},2} \right) & \hat{b}^{\dagger}_{\ve{R}} = \tfrac{1}{2} \left( i \hat{\eta}_{\ve{R},1} + \hat{\eta}_{\ve{R},2} \right). \end{align*}\]

In the Majorana basis, the Fermionic part of the Hamiltonian reads:

\[\begin{align*} \begin{split} \mathcal{H}_0^{C_{2v}} &= \frac{it}{2}\sum_{\ve{R}} \hat{\ve{\gamma}}^{\rm T}_{\ve{R}} \Big( - \hat{\ve{\eta}}_{\ve{R}} + \hat{\ve{\eta}}_{\ve{R}+\ve{e}_-} + \hat{\ve{\eta}}_{\ve{R}+\ve{e}_--\ve{e}_+} + \hat{\ve{\eta}}_{\ve{R}-\ve{e}_+} \Big) \\ \mathcal{H}_\text{Yuk}^{C_{2v}} &= \frac{i \xi}{2} \sum_{\ve{R}} s_{\ve{R}} \hat{\ve{\gamma}}^{\rm T}_{\ve{R}} \Big( - \hat{\ve{\eta}}_{\ve{R}} - \hat{\ve{\eta}}_{\ve{R}+\ve{e}_-} - \hat{\ve{\eta}}_{\ve{R}+\ve{e}_--\ve{e}_+} + \hat{\ve{\eta}}_{\ve{R}-\ve{e}_+} \Big) \end{split} \end{align*}\]

In the above, \(\ve{\hat{\gamma}}^{T}_{\ve{R}} = (\hat{\gamma}_{\ve{R},1} , \hat{\gamma}_{\ve{R},2} )\) and a similar form holds for \(\ve{\hat{\eta}}_{\ve{R}}\). The fact that the Hamiltonian is diagonal in the Majorana index shows that it has a higher \(\text{O}(2N_{\sigma})\) as opposed to the apparent \(\text{SU}(N_{\sigma})\) one in the fermion representation. It also has as consequence that, for the \(N_\sigma = 1\) case, the fermion determinant is nothing but the square of a Pfaffian that takes real values. Hence the negative sign problem is absent [68, 69].

We close this subsection by making contact with the work of Ref. [67], showing the absence of the sign problem with a alternative approach. Let \(\ve{\mu}\) be a vector of Pauli matrices acting on the Majorana index. Adopting the notation of Ref. [67], we can define:

\[\begin{eqnarray*} & & \hat{T}^{-}_{1} \alpha \ve{\hat{\gamma}}_{\ve{R}} \left( \hat{T}^{-}_{1} \right)^{-1} = \bar{\alpha} i \ve{\mu}_y \ve{\hat{\gamma}}_{\ve{R}} \nonumber \\ && \hat{T}^{-}_{1} \alpha \ve{\hat{\eta}}_{\ve{R}} \left( \hat{T}^{-}_{1} \right)^{-1} = - \bar{\alpha} i \ve{\mu}_y \ve{\hat{\eta}}_{\ve{R}} \end{eqnarray*}\]

and

\[\begin{eqnarray*} & & \hat{T}^{+}_{2} \alpha \ve{\hat{\gamma}}_{\ve{R}} \left( \hat{T}^{+}_{2} \right)^{-1} = \bar{\alpha} \ve{\mu}_x \ve{\hat{\gamma}}_{\ve{R}} \nonumber \\ && \hat{T}^{+}_{2} \alpha \ve{\hat{\eta}}_{\ve{R}} \left( \hat{T}^{+}_{2} \right)^{-1} = - \bar{\alpha} \ve{\mu}_x \ve{\hat{\eta}}_{\ve{R}}. \end{eqnarray*}\]

Since \(\left[ \hat{T}^{+}_{2} , \mathcal{H}^{C_{2v}} \right] = \left[ \hat{T}^{-}_{1} , \mathcal{H}^{C_{2v}} \right] =0\) and \(\hat{T}^{-}_{1} \) and \(\hat{T}^{+}_{2} \) anti-commute, our Hamiltonian belongs to the so-called Majorana class, and is hence free of the negative sign problem.

2.6.1.2. The \(C_{4v}\) model

Consider the spinor \( \hat{\ve{c}}^{\dagger}_{\ve{R}} = \left( \hat{a}^{\dagger}_{\ve{R}}, \hat{b}^{\dagger}_{\ve{R}} \right)\). With this notation, the fermionic part of the \(C_{4v}\) model takes the form:

\[\begin{align*} \begin{split} \mathcal{H}_0^{C_{4v}} &= -t \sum_{\ve{R} \in A, \delta = \pm} \left( \hat{\ve{c}}_{\ve{R}}^\dag \, \ve{e}_{+} \cdot {\ve{\tau}} \, \hat{\ve{c}}_{\ve{R} + \delta \ve{e}_x} + \hat{\ve{c}}_{\ve{R}}^\dag \, \ve{e}_{-} \cdot {\ve{\tau}} \, \hat{\ve{c}}_{\ve{R} + \delta \ve{e}_y} + \mathrm{H.c.} \right) \\ \mathcal{H}_\text{Yuk}^{C_{4v}} &= \xi \sum_{\ve{R}} s_{\ve{R}} \hat{\ve{c}}^\dag_{\ve{R}} \ve{e}_y \cdot \ve{\tau} \hat{\ve{c}}_{\ve{R}} \\ \mathcal{H}_\mu^{C_{4v}} &= \mu \sum_{\ve{R}} \hat{\ve{c}}^\dag_{\ve{R}} \hat{\ve{c}}_{\ve{R}} \end{split} \end{align*}\]

In the above, \(\ve{\tau}\) denotes a vector of Pauli matrices that act on the orbital space, \(\ve{e}_{\pm} = \frac{1}{\sqrt{2}} \left( \ve{e}_x \pm \ve{e}_y\right)\). Further, \(\ve{R} \in A\) denotes the sum of the A sub-lattice, \((-1)^{R_x+R_y } =1\). Furthermore, we have to consider an arbitrary set of Ising spins \( s_{\ve{R}} = \pm 1 \). Consider the relation

\[\begin{equation*} U^{\dagger} (\ve{e},\theta) \ve{\tau} U (\ve{e},\theta) = R(\ve{e},\theta) \ve{\tau}. \end{equation*}\]

\(U (\ve{e},\theta) = e^{- i \theta \ve{e} \cdot \ve{\tau}/2}\) is an SU(2) rotation of angle \(\theta \) around axis \(\ve{e}\) (\( |\ve{e}| = 1\) ) and \(R(\ve{e},\theta) \) an SO(3) with the same angle and axis. We can hence carry out a canonical transformation,

\[\begin{equation*} \hat{d}_{\ve{R}} = U \hat{c}_{\ve{R}}, \end{equation*}\]

that rotates \(\ve{e}_{+} \rightarrow \ve{e}_z \), \(\ve{e}_{-} \rightarrow \ve{e}_x\), and \(\ve{e}_{y} \rightarrow - \frac{1}{\sqrt{2}} \left( \ve{e}_x - \ve{e}_z \right)\) by combing a \(\pi/4\) rotation around the z-axis and subsequently a \(\pi/2\) rotation around the x-axis. After this canonical transformation, the Hamiltonian is real, and takes the form:

\[\begin{align*} \begin{split} \mathcal{H}_0^{C_{4v}} &= -t \sum_{\ve{R} \in A, \delta = \pm} \left( \hat{\ve{d}}_{\ve{R}}^\dag \, {\tau}_z \, \hat{\ve{d}}_{\ve{R} + \delta \ve{e}_x} + \hat{\ve{d}}_{\ve{R}}^\dag \, {\tau}_x \, \hat{\ve{d}}_{\ve{R} + \delta \ve{e}_y} + \mathrm{H.c.} \right) %\nonumber \\ \mathcal{H}_\text{Yuk}^{C_{4v}} &= - \frac{\xi}{\sqrt{2}} \sum_{\ve{R}} s_{\ve{R}} \hat{\ve{d}}^\dag_{\ve{R}} \left( \tau_x - \tau_z \right) \hat{\ve{d}}_{\ve{R}} \\ %\nonumber \mathcal{H}_\mu^{C_{4v}} &= \mu \sum_{\ve{R}} \hat{\ve{d}}^\dag_{\ve{R}} \hat{\ve{d}}_{\ve{R}} \end{split} \end{align*}\]

We can now express the model in terms of Majorana fermions and choose

\[\begin{eqnarray*} \hat{\ve{d}}^{\dagger}_{\ve{R}} = \frac{1}{2} \left( \hat{\ve{\gamma}}_{\ve{R},1} - i \hat{\ve{\gamma}}_{\ve{R},2} \right) \end{eqnarray*}\]

as representation for \((-1)^{R_x+R_y } =1\) and

\[\begin{eqnarray*} \hat{\ve{d}}^{\dagger}_{\ve{R}} = \frac{1}{2} \left( i \hat{\ve{\gamma}}_{\ve{R},1} + \hat{\ve{\gamma}}_{\ve{R},2} \right) \end{eqnarray*}\]

as representation for \((-1)^{R_x+R_y } =-1\). Let \(\ve{\mu}\) be a vector of Pauli spin matrices that acts on the Majorana index. With this choice, the Hamiltonian then takes the form:

\[\begin{align*} \begin{split} \mathcal{H}_0^{C_{4v}} &= \frac{i t}{2} \sum_{\ve{R} \in A, \delta = \pm} \left( \hat{\ve{\gamma}}_{\ve{R}}^{\rm T} \, {\ve{\tau}_z} \, \hat{\ve{\gamma}}_{\ve{R} + \delta \ve{e}_x} + \hat{\ve{\gamma}}_{\ve{R}}^{T} \, {\ve{\tau}}_x \, \hat{\ve{\gamma}}_{\ve{R} + \delta \ve{e}_y} \right) \\ \mathcal{H}_\text{Yuk}^{C_{4v}} &= \frac{\xi}{4\sqrt{2}} \sum_{\ve{R}} s_{\ve{R}} \hat{\ve{\gamma}}^{\rm T}_{\ve{R}} \left( \ve{\tau}_x \ve{\mu}_{y}- \ve{\tau}_z \ve{\mu}_y \right) \hat{\ve{\gamma}}_{\ve{R}} \\ \mathcal{H}_\mu^{C_{4v}} &= \mu \sum_{\ve{R}} \left( 2 - \hat{\ve{\gamma}}^{\rm T}_{\ve{R}} \mu_y \hat{\ve{\gamma}}_{\ve{R}} \right) \end{split} \end{align*}\]

Using the notation of Ref. [67] we define:

\[\begin{equation*} \hat{T}^{-}_{1} \alpha \hat{\ve{\gamma}}_{\ve{R}} \left( \hat{T}^{-}_{1} \right)^{-1} = \bar{\alpha} i \ve{\tau}_y \ve{\mu}_x \hat{\ve{\gamma}}_{\ve{R}} \end{equation*}\]

and

\[\begin{equation*} \hat{T}^{-}_{2} \alpha \hat{\ve{\gamma}}_{\ve{R}} \left( \hat{T}^{-}_{2} \right)^{-1} = \bar{\alpha} i \ve{\tau}_y \ve{\mu}_z \hat{\ve{\gamma}}_{\ve{R}} \end{equation*}\]

that satisfy

(2.21)\[\begin{equation} \left[ \mathcal{H}^{C_{4v}}, \hat{T}^{-}_{1} \right] = \left[ \mathcal{H}^{C_{4v}}, \hat{T}^{-}_{2} \right] = 0. \end{equation}\]

Both above symmetries square to (-1) and anti-commute with each other. This hence places us in the Kramers class, see Ref. [67, 70], and no negative sign problem occurs.

2.7. QMC Observables

In the following, we define the QMC observables used throughout this project to study the quantum phase transition. We have used quantities based on both bosonic and fermionic degrees of freedom.

2.7.1. Bosonic degrees of freedom

2.7.1.1. Order parameters

The structure factor \(S(\ve{k})\) and susceptibility \(\chi(\ve{k})\) are defined as

(2.22)\[S(\ve{k}) = \frac{1}{N^2}\sum_{\langle \ve{R}, \ve{R}'\rangle} \left( \langle s^z_{\ve{R}} s^z_{\ve{R}'}\rangle - \langle s^z_{\ve{R}}\rangle \langle s^z_{\ve{R}'}\rangle \right) e^{i\ve{k} (\ve{R} - \ve{R}')},\]

and

(2.23)\[ \chi(\ve{k}) = \int\!\text{d}\tau\, \frac{1}{N^2}\sum_{\langle \ve{R},\ve{R}'\rangle} \left( \Braket{ s^z_{\ve{R}}(0) s^z_{\ve{R}'}(\tau) } - \Braket{ s^z_{\ve{R}}(0) } \Braket{ s^z_{\ve{R}'}(\tau) } \right) e^{i\ve{k} (\ve{R} - \ve{R}')}.\]

Both \(S(\ve{k}=0)\) and \(\chi(\ve{k}=0)\) are suitable order parameters to probe the paramagnetic-ferromagnetic phase transition.

2.7.1.2. RG-invariant quantities

Next, we define a set of renormalization group (RG) invariant observables. These are quantities with vanishing scaling dimension, that in the thermodynamic limit either converge to \(1\) if the state is ordered (ferromagnetically, in the case of this project), or to \(0\) in the absence of order.

At a quantum critical point, RG-invariant quantities follow the form \(f[L^z/\beta, (h-h_\mathrm{c})L^{1/\nu}, L^{-\Delta z}, L^{-\omega}]\). Here, we have taken into account the possibility of two characteristic length scales: \(\Delta z = 1 - z_-/z_+\). Since our temperature is representative of the ground state, we can neglect the dependence on \(L^{z}/\beta\). The term includes two corrections to scaling, the regular term —governed by an exponent \(\omega\)— and another correction that is present if \(z_- \neq z_+\). Up to these corrections, the data for different lattice sizes cross at the critical field \(h_\mathrm{c}\).

The first RG-invariant quantity is the well-known Binder ratio, \(B\) [34], defined as

(2.24)\[ B = \frac{1}{2}\left( 3 - \frac{ \Braket{(s^z)^4} }{ \Braket{(s^z)^2}^2 } \right).\]

Furthermore, we define two more RG-invariant quantities \(R_S\) and \(R_\chi\) through the correlation ratio: Given a local order parameter \(O\) at momentum \(\ve{p}\), one can define the correlation ratio \(R_O\) as

(2.25)\[ R_O \equiv 1 - \frac{O(\ve{p}+\delta\ve{p})}{O(\ve{p})} \quad \text{with} \quad O \in \{S, \chi\},\]

where \(O(\ve{p})\) is the two-point function of the order parameter in Fourier space, and \(\delta\ve{p}\) is the minimum nonzero momentum on a finite lattice. For the \(C_{2v}\) model \(\delta\ve{p}=(\pi/L,\pi/L)\) or \(\delta\ve{p}=(\pi/L,-\pi/L)\), while for the \(C_{4v}\) model \(\delta\ve{p}=(2\pi/L,0)\) or \(\delta\ve{p}=(0,2\pi/L)\); as usual, one can average over the two minimum displacements to obtain an improved estimator. In principle \(R_O\) has the same asymptotical behavior for bigger values of \(\delta\ve{p}\), as long as they scale with \(1/L\), however, we have found that using the minimal versions works best for us.

2.7.1.3. Derivative of the free energy

To provide further information on the nature of the transition, we use the derivative of the free energy,

(2.26)\[\begin{align} \frac{1}{N}\frac{\partial F}{\partial h} = \Braket{ \frac{1}{N} \sum_{\ve{R}} s^x_{\ve{R}} } \equiv X. \end{align}\]

2.7.2. Fermionic degrees of freedom

The fermionic observables consist of the momentum-resolved single-particle gap \(\Delta_\text{sp}(\ve{k})\) which we use to image the meandering of Dirac points. We furthermore use this quantity to determine the Fermi velocity anisotropy \(v_\perp / v_\parallel\).

2.7.2.1. Fermionic single-particle gap

To properly define \(\Delta_\text{sp}(\ve{k})\), we first introduce an energy basis:

\[\begin{align*} \mc{H} \Ket{\Psi^N_n(\ve{k})} = E^N_n(\ve{k}) \Ket{\Psi^N_n(\ve{k})}, \end{align*}\]

where \(\Ket{\Psi^N_n(\ve{k})}\) are also eigenstates of particle number \(\hat{N}\) and momentum \(\hat{\ve{k}}\) operators:

\[\begin{align*} \hat{N} \Ket{\Psi^N_n(\ve{k})} &= N \Ket{\Psi^N_n(\ve{k})} & \hat{\ve{k}} \Ket{\Psi^N_n(\ve{k})} &= \ve{k} \Ket{\Psi^N_n(\ve{k})} \end{align*}\]

In this basis, the gap is:

\[\begin{align*} \Delta_\text{sp}(\ve{k}) = E_0^{N_0+1}(\ve{k}) - E_0^{N_0}, \end{align*}\]

where \(N_0\) is the particle number of the half-filled system.

Now consider the time-displaced Green function

\[\begin{align*} G(\ve{k},\tau) = \Braket{\hat{c}_{\ve{k}}(\tau) \hat{c}_{\ve{k}}^\dag} \quad \text{with} \quad \hat{c}_{\ve{k}}(\tau) = e^{\tau \mc{H}} \hat{c}_{\ve{k}} e^{-\tau \mc{H}}. \end{align*}\]

Assuming a unique ground state, the \(T=0\) Green function reads:

\[\begin{align*} \lim_{\beta \rightarrow \infty} G(\ve{k},\tau) = \Braket{\Psi^{N_0}_0 | \hat{c}_{\ve{k}}(\tau) \hat{c}_{\ve{k}}^\dag | \Psi^{N_0}_0} = \sum_n e^{-\tau \left(E_n^{N_0+1}(\ve{k}) - E_0^{N_0}\right)} \left|\Braket{\Psi^{N_0+1}_n(\ve{k}) | \hat{c}_{\ve{k}}^\dag | \Psi^{N_0}_0} \right|^2. \end{align*}\]

Provided that the wave function renormalization, \(\left|\Braket{\Psi^{N_0+1}_n(\ve{k}) | \hat{c}_{\ve{k}}^\dag | \Psi^{N_0}_0} \right|^2\) is finite and that \(\ket{\Psi^{N_0+1}_0(\ve{k})}\) is non-degenerate, then

(2.27)\[ \lim_{\tau \rightarrow \infty} \lim_{\beta \rightarrow \infty} G(\ve{k},\tau) = e^{-\tau \left(E_0^{N_0+1}(\ve{k}) - E_0^{N_0}\right)} \left|\Braket{\Psi^{N_0+1}_0(\ve{k}) | \hat{c}_{\ve{k}}^\dag | \Psi^{N_0}_0} \right|^2\]

and we can extract \(\Delta_\text{sp}(\ve{k}) = E_0^{N_0+1}(\ve{k}) - E_0^{N_0}\) by fitting the tail of \(G(\ve{k},\tau)\) to an exponential form.

The source code used for this fit is displayed in Appendix A.4 and Appendix A.2.2.10 demonstrates how to apply the fitting function.

In Fig. 2.11 we show that this approach works, by comparing the dispersions deep in the disordered and ordered phases to mean field results. Note that in a fully ergodic QMC simulation we would sample both options for breaking the Ising \(\mathbb{Z}_2\) symmetry. To produce the results of Figs. 2.11(a2,c2), we have omitted the global move that flips all the spins and comes with a unit acceptance.

We observe a slight systematic derivation between mean field results and QMC data in the disordered phase. This stems from fluctuations of the order parameter in the vicinity of the critical point.

Fig. 2.11 Testing the dispersion of the fermionic single particle gap obtained from QMC data against mean field results. The first column shows the mean field results according to Eq. (2.13), while the second column shows numerical results obtained by employing Eq. (2.27) and the last column shows mean field minus QMC results. (a) \(C_{2v}\) model in ordered phase, at \(h=1.0\). The value of the mean field parameter \(\phi\) is set to \(\langle s^z \rangle\) from the QMC simulation. (b) \(C_{2v}\) model in disordered phase, at \(h=5.0\). The value of the mean field parameter is set to \(\phi=0\). (c) Same as (a), but for the \(C_{4v}\) model. (d) Same as (b), but for the \(C_{4v}\) model.

2.7.2.2. Fermi velocity anisotropy \(v_\perp / v_\parallel\)

With \(\Delta_\text{sp}(\ve{k})\) we can extract the anisotropy at the nodal points \(\ve{K}\), via,

(2.28)\[\frac{v_\perp}{v_\parallel} = \lim_{\delta \rightarrow 0} \frac{\Delta_\text{sp}(\ve{K}+\delta \ve{e}_\perp) - \Delta_\text{sp}(\ve{K})} {\Delta_\text{sp}(\ve{K}+\delta \ve{e}_\parallel) - \Delta_\text{sp}(\ve{K})}.\]

Where \(\ve{e}_\perp\) and \(\ve{e}_\parallel\) are unit vectors perpendicular and parallel to the meandering direction of the Dirac cones. We have considered three different strategies for approaching this limit on the finite size lattices, that are all equivalent in the thermodynamic limit:

1. The direct approach. The most straightforward implementation of Eq. (2.28) on a finite lattice would be

(2.29)\[\frac{v_\perp}{v_\parallel} = \frac{\Delta_\text{sp}(\ve{K}+\ve{\delta}_\perp) - \Delta_\text{sp}(\ve{K})} {\Delta_\text{sp}(\ve{K}+\ve{\delta}_\parallel) - \Delta_\text{sp}(\ve{K})},\]

where \(\ve{\delta}_\perp\), \(\ve{\delta}_\parallel\) are the shortest distances from the nodal point on the finite-size \(\ve{k}\)-Lattice.

2. Manually setting the finite size gap \(\Delta_\text{sp}(\ve{K}) = 0\). This approach makes sense, since we know that in the thermodynamic limit the gap vanishes. With this strategy, Eq. (2.28) takes the form

(2.30)\[\frac{v_\perp}{v_\parallel} = \frac{\Delta_\text{sp}(\ve{K}+\ve{\delta}_\perp)} {\Delta_\text{sp}(\ve{K}+\ve{\delta}_\parallel)}.\]

3. Avoid the nodal points. Another approach for avoiding the finite size gap is to measure one step away from it:

(2.31)\[\frac{v_\perp}{v_\parallel} = \frac{\Delta_\text{sp}(\ve{K}+2\ve{\delta}_\perp) - \Delta_\text{sp}(\ve{K}+\ve{\delta}_\perp)} {\Delta_\text{sp}(\ve{K}+2\ve{\delta}_\parallel) - \Delta_\text{sp}(\ve{K}+\ve{\delta}_\parallel)}\]

The results for these different approaches are shown in Fig. 2.12. The third strategy results in velocity anisotropies \(<1\), while Fig. 2.1(c) clearly shows that \(v_\perp / v_\parallel > 1\) at the critical point. This implies that the approach strongly underestimates the anisotropy due to the fact that the considered lattices sizes are too small for not measuring directly at the nodal point.

The other two approaches, while not giving quantitatively the same results, are qualitatively equivalent. We have opted to use the second strategy, corresponding to Eq. (2.30).

Fig. 2.12 Fermi velocity anisotropy at the critical point as function of \(1/L\), determined with (a): Eq. (2.29), (b): Eq. (2.30), (c): Eq. (2.31). Power law and logarithmic fits are shown, except for (c), where only a logarithmic fit is performed.

2.8. QMC results

In this section, we first give an overview of the QMC results, showing that there is indeed a continuous phase transition, after which we determine the critical exponents and finally discuss odd-even effects found in the \(C_{4v}\) model.

2.8.1. Overview

Fig. 2.13 displays the previously defined bosonic observables close to the phase transition. The RG-invariant quantities defined in Eqs. (2.24) and (2.25) shown in Figs. 2.13(a-c) cross for different system sizes at \(h_\mathrm{c} \approx 3.25\) (\(h_\mathrm{c} \approx 3.65\)) for the \(C_{2v}\) (\(C_{4v}\)) model, signifying a phase transition at that point. Both the order parameter \(S(\ve{k}=0)\) (Fig. 2.13(d)) and derivative of the free energy with respect to the tuning parameter, \(\partial F/\partial h\) (Fig. 2.13(e)) do not show any discontinuity, supporting a continuous phase transition.

Due to finite size effects, the RG-invariant quantities do not cross at the same point, but display a drift with system size. Fig. 2.14(c,d) show the crossing points between \(L\) and \(L+\Delta L\) lattices, with \(\Delta L=2\) \((4)\) for the \(C_{2v}\) (\(C_{4v}\)) model. As apparent, we obtain consistent results for \(h_\mathrm{c}\) when considering different RG-invariant quantities. We estimate the correlation-length exponents \(1/\nu\) by data collapse for the two models in Fig. 2.14(a,b). Considering values of \(L \geq L_{\text{min}} =12\) we obtain \(1/\nu = 1.376(6)\) [\(1/\nu = 1.38(1)\)] for the \(C_{2v}\) [\(C_{4v}\)] model. These values are in the ballpark of the \(\epsilon\)-expansion results in the quasiuniversal regime (cf. Appendix A.1). Section 2.8.2 will give data for various values of \(L_{\text{min}}\), that are standing in agreement with the above values. Although seemingly converged, the fact that the velocity anisotropy is expected to flow extremely slowly suggest that the exponents are subject to considerable size effects, see Section 2.5. The impact of critical fluctuations on the fermion spectrum is displayed in Figs. 2.1(b,c). In the disordered phase, Fig. 2.1(b), the dispersion relation exhibits rotational symmetry around the Dirac points. On the other hand, at criticality, Fig. 2.1(c), the dispersion relation suggests a velocity anisotropy, \(v_{\parallel} < v_{\perp}\) at the Dirac point. Figure 2.14(e) demonstrates that this anisotropy grows as a function of system size, in qualitative agreement with the RG predictions. Although our system sizes are too small to detect convergence or divergence of the velocity ratio, we find it reassuring that its dependence on system size qualitatively resembles the scale dependence predicted from the integrated RG flow; cf. Fig. 2.14(e) with Fig. 2.10.

Fig. 2.13 Bosonic observables as function of transverse Ising field \(h\) close to the critical point \(h=h_\mathrm{c}\). (a-c) The RG-invariant quantities: Structure factor correlation ration \(R_S\), Binder ratio \(B\) and Susceptibility correlation ratio \(R_{\chi}\). The point where lines for different system sizes cross indicates a phase transition. (d) Order parameter \(S(\ve{k}=0)\) [Eq. (2.22)]. (e) Derivative of free energy, exhibiting no discontinuities.

Fig. 2.14 (a) \(R_S\) as function of \((h-h_\mathrm{c}) L^{1/\nu}\) for the \(C_{2v}\) model, revealing a data collapse for \(L \gtrsim 12\), assuming \(1/\nu = 1.376\). (b) Same as (a), but for \(C_{4v}\) model, assuming \(1/\nu = 1.38\). (c) Crossing points of different RG-invariant quantities as function of \(1/L\) with \(\Delta L = 2\) in \(C_{2v}\) model, indicating a unique critical point \(h_\mathrm{c} = 3.27\) for \(L \to \infty\). (d) Same as (c), but for \(C_{4v}\) model and \(\Delta L = 4\), extrapolating to \(h_\mathrm{c} = 3.65\). (e) Ratio of Fermi velocities \(v_\perp/v_\parallel\) as function of \(1/L\) at \(h_\mathrm{c}\), revealing that the velocity anisotropy increases with increasing system size. The solid lines show power law fits for \(L \geq 8\) and logarithmic fits for \(L \geq 12\).

2.8.2. Critical exponents

2.8.2.1. Correlation length exponent \(\nu\) from RG invariant quantities.

A renormalization-group invariant quantity, \(R\), has by definition a vanishing scaling dimension. Consider a system at temperature \(\beta\), of size \(L_{+} \times L_{-}\) with a single relevant coupling \(h\). Under a renormalization group transformation that rescales \(L_{+} \rightarrow L_{+}/b\) with \(b>1\), we expect [71]:

(2.32)\[\begin{equation} R( (h - h_\mathrm{c}) , \beta, L_{+}, L_{-}) = R( (h-h_\mathrm{c})', \beta/b^{z}, L_{+}/b, L_{-}/b^{1+ \Delta z }). \end{equation}\]

In the above \(\Delta z \neq 0 \) encodes the difference in scaling between the \(L_{+} \) and \(L_{-}\) directions. Linearization of the RG transformation, \((h-h_\mathrm{c})' = b^{1/\nu} (h-h_\mathrm{c})\) and setting the scale \(b= L\) as well as \(L_{-} = L_{+} = L \), in accordance to our simulations, yields:

(2.33)\[ R \left( (h - h_\mathrm{c}) , \beta, L \right) = f \left( L^{1/\nu} (h -h_\mathrm{c}), L^z/\beta , L^{ - \Delta z } , L^{-\omega} \right).\]

In the above we have accounted for possible corrections to scaling \(L^{-\omega}\). In the presence of a single length scale \(\Delta z = 0\), such that the generic finite size scaling form is recovered.

Since in our simulations the temperature is representative of the ground state, we can neglect the dependence on \(L^{z}/\beta\). Up to corrections to scaling, \(\omega\), and the possibility of \(\Delta z \neq 0\), which would result in another correction to scaling term, the data for different lattice sizes cross at the critical field \(h_\mathrm{c}\) and should collapse when plotted as function of \((h-h_\mathrm{c})L^{1/\nu}\). The results for such data collapses are shown in Table 2.1 and Table 2.2 (cf. Appendix A.2.2.8 for a demonstration on how such a data collapse can be carried out). Furthermore, Fig. 2.15 shows \(1/\nu\) for the \(C_{2v}\) and \(C_{4v}\) models from pairwise data collapse of RG-invariant quantities, using system sizes \(L\) and \(L+2\) (\(L+4\)). Both suggest a relatively well converged result for \(L \geq 12\). Although seemingly converged, our system sizes are too small to detect the logarithmic drift in exponents suggested by the RG analysis.

Table 2.1 Data collapse results for RG-invariant quantities of \(C_{2v}\) model.

Observables	Used system sizes	\(h_\mathrm{c}\)	\(1/\nu\)	\(\chi^2\)
\(R_S\)	8, 10, 12, 14, 16, 18, 20	\(3.272715 \pm 0.000074\)	\(1.358934 \pm 0.001824\)	2.4
\(R_S\)	10, 12, 14, 16, 18, 20	\(3.272222 \pm 0.000065\)	\(1.373499 \pm 0.002887\)	1.8
\(R_S\)	12, 14, 16, 18, 20	\(3.272304 \pm 0.000099\)	\(1.376110 \pm 0.006104\)	1.9
\(R_S\)	14, 16, 18, 20	\(3.272617 \pm 0.000138\)	\(1.375085 \pm 0.007038\)	1.8
\(R_S\)	16, 18, 20	\(3.272521 \pm 0.000274\)	\(1.368184 \pm 0.014756\)	1.9
\(R_S\)	18, 20	\(3.273322 \pm 0.000449\)	\(1.373454 \pm 0.022000\)	1.9
\(B\)	8, 10, 12, 14, 16, 18, 20	\(3.270955 \pm 0.000089\)	\(1.317011 \pm 0.002273\)	13.4
\(B\)	10, 12, 14, 16, 18, 20	\(3.271532 \pm 0.000141\)	\(1.344485 \pm 0.004371\)	4.6
\(B\)	12, 14, 16, 18, 20	\(3.272535 \pm 0.000175\)	\(1.352409 \pm 0.006242\)	3.0
\(B\)	14, 16, 18, 20	\(3.273181 \pm 0.000152\)	\(1.361056 \pm 0.007298\)	2.6
\(B\)	16, 18, 20	\(3.273783 \pm 0.000210\)	\(1.370788 \pm 0.013022\)	2.2
\(B\)	18, 20	\(3.274325 \pm 0.000464\)	\(1.340845 \pm 0.028414\)	1.8
\(R_\chi\)	8, 10, 12, 14, 16, 18, 20	\(3.281138 \pm 0.000030\)	\(1.421387 \pm 0.000002\)	22.0
\(R_\chi\)	10, 12, 14, 16, 18, 20	\(3.279752 \pm 0.000073\)	\(1.381046 \pm 0.000002\)	7.3
\(R_\chi\)	12, 14, 16, 18, 20	\(3.275155 \pm 0.000500\)	\(1.369774 \pm 0.001609\)	5.8
\(R_\chi\)	14, 16, 18, 20	\(3.277021 \pm 0.000233\)	\(1.338662 \pm 0.010808\)	2.2
\(R_\chi\)	16, 18, 20	\(3.276434 \pm 0.000176\)	\(1.342788 \pm 0.004183\)	2.2
\(R_\chi\)	18, 20	\(3.275856 \pm 0.000693\)	\(1.369095 \pm 0.034679\)	2.4

Table 2.2 Data collapse results for RG-invariant quantities of \(C_{4v}\) model.

Observables	Used system sizes	\(h_\mathrm{c}\)	\(1/\nu\)	\(\chi^2\)
\(R_S\)	8, 12, 16, 20	\(3.64606 \pm 0.00007\)	\(1.328 \pm 0.006\)	18.7
\(R_S\)	12, 16, 20	\(3.64886 \pm 0.00011\)	\(1.381 \pm 0.011\)	3.2
\(R_S\)	16, 20	\(3.65108 \pm 0.00022\)	\(1.402 \pm 0.023\)	1.7
\(B\)	8, 12, 16, 20	\(3.64108 \pm 0.00009\)	\(1.254 \pm 0.007\)	73.2
\(B\)	12, 16, 20	\(3.64818 \pm 0.00012\)	\(1.340 \pm 0.014\)	5.1
\(B\)	16, 20	\(3.65138 \pm 0.00025\)	\(1.362 \pm 0.026\)	1.8
\(R_\chi\)	8, 12, 16, 20	\(3.66319 \pm 0.00027\)	\(1.309 \pm 0.018\)	16.5
\(R_\chi\)	12, 16, 20	\(3.65708 \pm 0.00031\)	\(1.368 \pm 0.028\)	3.1
\(R_\chi\)	16, 20	\(3.65537 \pm 0.00070\)	\(1.428 \pm 0.092\)	2.4

Fig. 2.15 Critical exponent \(1/\nu\) of \(C_{2v}\) (\(C_{4v}\)) model from pairwise data collapse of RG-invariant quantities, using linear system sizes \(L\) and \(L+2\) (\(L+4\)).

2.8.2.2. Scaling dimensions and scaling anisotropy

Next, we examine the scaling dimension of the bosonic field from the Ising spin correlations:

(2.34)\[\begin{align} S(\ve{x}) = \Braket{\hat{s}^z_{\bf R}(\tau) \hat{s}^z_{{\bf 0}}(0)} \end{align}\]

where \(\ve{x} = (\ve{R}, \tau)\) is a space-time coordinate. The models considered in this research are not Lorentz invariant such that the scaling dimension acquires a direction dependence. Following Eq. (2.33) we expect:

(2.35)\[\begin{align} S\left( r \hat{\ve{d}}_*, h \right) \propto \frac{1}{ | r \hat{\ve{d}}_*|^{2 \Delta_{s,*}}} f\left( L^z/\beta, (h-h_\mathrm{c})L^{1/\nu}, L^{-\Delta z} , L^{-\omega} \right) \end{align}\]

where \(\hat{\ve{d}}_*\) defines the direction.

To determine the scaling dimensions, we consider \(S(L \hat{\ve{d}}_*, h)\), for different system sizes \(L\) and use an RG-invariant quantity \(R\) to replace in leading order \(f( L^z/\beta, (h-h_\mathrm{c})L^{1/\nu}, L^{- \Delta z }, L^{-\omega} ) = \tilde{f}(R)\). Using this form, we perform data collapses using system sizes \(L\) and \(L+2\) (\(L\) and \(L+4\)), where the only free parameter is \(\Delta_{s,*}\). The considered directions are defined in Table 2.3, the \(C_{4v}\) symmetry of the second model enforces \(\Delta_{s,x}=\Delta_{s,y}\) and \(\Delta_{s,+}=\Delta_{s,-}\). As the results in Fig. 2.16 and Fig. 2.17 show, we cannot resolve a scaling anisotropy between the chosen directions. We conjecture that anisotropies in the exponents will emerge in the infrared limit. Given the very slow flow we believe that our numerical simulations are not in a position to probe these energy scales.

Table 2.3 Considered directions for the scaling dimension.

\(*\)	\(\hat{\ve{d}}_*\)
\(x\)	\((\hat{\ve{e}}_x, 0)\)
\(y\)	\((\hat{\ve{e}}_y, 0)\)
\(+\)	\(\frac{1}{2}(\hat{\ve{e}}_x+\hat{\ve{e}}_y, 0)\)
\(-\)	\(\frac{1}{2}(\hat{\ve{e}}_x-\hat{\ve{e}}_y, 0)\)
\(\tau\)	\((\ve{0}, 0.3)\)

Fig. 2.16 Scaling dimension of Ising field of \(C_{2v}\) model. For (a) \(\xi=0.25\), (b) \(\xi=0.4\). Note: \(\Delta_{s,y}\), \(R=B\) is indistinguishable from \(\Delta_{s,x}\), \(R=B\) and \(\Delta_{s,y}\), \(R=R_S\) is identical to \(\Delta_{s,x}\), \(R=R_S\).

Fig. 2.17 Scaling dimension of Ising field of \(C_{4v}\) model.

2.8.2.3. Dynamical exponent \(z\)

To determine the dynamical exponent of the \(C_{4v}\) model, we assume isotropic scaling in space, as suggested by the RG analysis. Then the RG-invariant quantities follow the form

(2.36)\[ R = f( L^z/\beta, (h-h_\mathrm{c})L^{1/\nu}, L^{-\omega} )\]

at the critical point.

At the crossing points \(h_*(L)\), with \(R(h_*(L), L) = R(h_*(L), L+\Delta_L)\) and \(\Delta_L = 4\), we measure \(R(\beta)\). Omitting corrections to scaling leads to \(R(L, \beta) = f( L^z/\beta)\). From this we derive

(2.37)\[\begin{align} z = \frac {\log\left( \frac{\partial_\beta R(L)}{\partial_\beta R(L+\Delta_L)} \right)} {\log\left(\frac{L+\Delta_L}{L} \right))}. \end{align}\]

The results are shown in Fig. 2.18, and are consistent with \(z=1\) as suggested in the RG analysis.

Fig. 2.18 Dynamical exponent \(z\) of \(C_{4v}\) model.

2.8.3. Odd-even effects

The \(C_{4v}\) model has strong odd-even effects. For linear system sizes \(L \in 4\mathbb{N}\) (\(\equiv\)even) and periodic boundary conditions, the Dirac points are included in the discrete set of \(\ve{k}\) vectors. This is not the case for odd lattices, \(L \in 4\mathbb{N}+2\). Interestingly, the value of the Binder and correlation ratios depend on this choice of the boundary, see Fig. 2.19(b),(c),(d). We believe that this stems form the fact that both quantities do not have a well defined thermodynamic limit at \(h = h_\mathrm{c}\). i.e. \(\lim_{L\rightarrow \infty} R_O( h = h_\mathrm{c})\) is mathematically not defined. However, the free energy (Fig. 2.19(a)), the critical field (Fig. 2.20(a)) and the exponents (Figs. 2.20(b-d)), should ultimately converge to the same value. On odd lattices, finite size effects seem to be larger.

The critical exponents \(2\beta/\nu\) and \(\eta\) in Figs. 2.20(c,d), stem from the scaling assumptions

(2.38)\[\begin{align} S(\ve{k}=0, h=h_\mathrm{c}, L) &\propto L^{2\beta/\nu} & \chi(\ve{k}=0, h=h_\mathrm{c}, L) &\propto L^{2-\eta_\phi}, \end{align}\]

where we omitted, as before, the dependence on the inverse temperature \(\beta\) and on corrections to scaling. Replacing \(h_\mathrm{c}\) by the crossing point \(h_*(L)\) of an RG-invariant quantity \(R\), meaning \(R(h_*(L), L) = R(h_*(L), L+4)\) with \(R \in \{R_S, R_\chi, B\}\), we obtain:

(2.39)\[ 2\beta/\nu = \log\left( \frac{S(\ve{k}=0, L+4, h=h_*(L)} {S(\ve{k}=0, L, h=h_*(L)} \right) / \log\left( \frac{L+4}{L} \right),\]

and

(2.40)\[ \eta_\phi = 2 - \log\left( \frac{\chi(\ve{k}=0, L+4, h=h_*(L)} {\chi(\ve{k}=0, L, h=h_*(L)} \right) / \log\left( \frac{L+4}{L} \right).\]

As apparent in Fig. 2.20, \(h_\mathrm{c}\) has the smallest corrections to scaling when determined from \(R_S\). However, the smallest corrections to scaling are when determining the critical exponents \(2\beta/\nu\), \(\eta_\phi\) and \(z\), are obtained by using \(h_\mathrm{c}\) as determined from \(R_\chi\). Finally, the velocity anisotropy at the critical point grows in both cases, but is much smaller for odd system sizes, cf. Fig. 2.21.

Fig. 2.19 Derivative of free energy and three RG-invariant quantities, showing a continuous transition around \(h\approx 3.65\). Notable is an odd-even effect between linear system sizes \(L \in 4\mathbb{N}\) (=even) and \(L \in 4\mathbb{N}+2\) (=odd).

Fig. 2.20 Demonstration of odd-even effects for the \(C_{4v}\) model. (a): Critical field \(h_\mathrm{c}\), extracted from the three RG-invariant quantities as determined by the crossing points between linear system sizes \(L\) and \(L+4\). Odd and even system sizes show different behavior. Even system size shows better convergence. (b): Critical exponent \(1/\nu\) as determined by data collapse of the three RG-invariant quantities, correlation ratio \(R\), Binder ratio \(B\) and susceptibility ratio \(R_\chi\), for linear system sizes \(L\) and \(L+4\). Odd and even system sizes show different behavior. Even system size shows better convergence. (c): Critical exponent \(2\beta/\nu\) as determined with Eq. (2.39). (d): Critical exponent \(\eta_\phi\) as determined with Eq. (2.40).

Fig. 2.21 Odd-even effects for the \(C_{4v}\) model on the anisotropy velocity of Dirac cones at the critical point.

2.9. Summary

Both the \(\epsilon\)-expansion analysis and the QMC simulations show that our two symmetry distinct models of Dirac fermions support continuous nematic transitions. In both cases, the key feature of the quantum critical point is a velocity anisotropy that is best seen in the QMC data of Fig. 2.1(c). For the \(C_{4v}\) model, the \(\epsilon\)-expansion shows that it diverges logarithmically with system size, in agreement with previous large-\(N\) results [52]. This law is supported by finite-size analysis based on QMC data up to linear system size \(L = 20\), which is close to the upper bound allowed by current computational approaches. Since the effective exponents flow with the velocity anisotropy, we foresee that lattice sizes beyond the reach of our numerical approach and experiments at ultralow temperatures will be required to obtain converged values. The QMC data captures a quasiuniversal regime [65, 66], in which irrelevant operators aside from the velocity anisotropy die out. In fact, the RG prediction for exponents in this intermediate-energy regime is roughly consistent with the finite-size QMC measurements, Fig. 2(c). Furthermore, for a reasonable set of starting values, the integrated RG flows of the two models are initially very similar and deviate from each other only at very low energy scales. A similar behavior of the two models is also observed in the QMC data.

An advantage of our models is that the Dirac points are pinned by symmetry, such that QMC approaches that take momentum-space patches around these points into account [72] represent an attractive direction for future work. Our models equally allow for large-\(N\) generalizations, such that QMC and analytical large-\(N\) calculations can be compared as a function of increasing \(N\). Finally, we can make contact to nematic transitions in \((2+1)\)-dimensional Fermi liquids [48, 49], since our models do not suffer from the negative-sign problem under doping.