Numerical evidence of conformal phase transition in graphene with long-range interactions

TL;DR

This study uses advanced simulations to investigate phase transitions in graphene, finding that a scaled-up interaction strength induces a semimetal to spin-density wave transition, with evidence suggesting a conformal phase transition rather than a traditional universality class.

Contribution

The paper provides the first unbiased simulation evidence of a conformal phase transition in graphene with long-range interactions and introduces an improved HMC algorithm for such studies.

Findings

A 1.6x increase in interaction strength induces a SDW transition in graphene.

No evidence found for charge-density wave order.

Results favor a conformal phase transition over the predicted universality class.

Abstract

Using state of the art Hybrid-Monte-Carlo (HMC) simulations we carry out an unbiased study of the competition between spin-density wave (SDW) and charge-density wave (CDW) order in suspended graphene. We determine that the realistic inter-electron potential of graphene must be scaled up by a factor of roughly 1.6 to induce a semimetal-SDW phase transition and find no evidence for CDW order. A study of critical properties suggests that the universality class of the three-dimensional chiral Heisenberg Gross-Neveu model with two fermion flavors, predicted by renormalization group studies and strong-coupling expansion, is unlikely to apply to this transition. We propose that our results instead favor an interpretation in terms of a conformal phase transition. In addition, we describe a variant of the HMC algorithm which uses exact fermionic forces during molecular dynamics trajectories and avoids the use of pseudofermions. Compared to standard HMC this allows for a substantial increase of the integrator stepsize while achieving comparable Metropolis acceptance rates and leads to a sizable performance improvement.

Figures

Equations

$$\hat{\mathcal{H}}=-\kappa \sum _{⟨x,y⟩,\sigma }({\hat{c}}_{x,\sigma }^{†}{\hat{c}}_{y,\sigma }+\text{h.c.})+\frac{\lambda }{2}\sum _{x,y}{\hat{\rho }}_x{V}_{xy}{\hat{\rho }}_y.$$

$$V(r)=\left\{\begin{array}[]{lr}U_{00},U_{01},U_{02},U_{03}&,\;r\leq 2a\\[2.0pt] e^{2}\left(c\,\displaystyle\frac{\exp(-m_{D}r)}{a(r/a)^{\gamma}}+m_{0}\right)&,\;r>2a\end{array}\right.$$

$${\hat{c}}_{x,\uparrow },{\hat{c}}_{x,\uparrow }^{†}$$

$${\hat{c}}_{x,\downarrow },{\hat{c}}_{x,\downarrow }^{†}$$

$$\frac{V_{xx}}{2}{\hat{\rho }}_x^2=$$

$$+V_{xx}(1-\eta ){\hat{\rho }}_x^{\prime },$$

$$\mathcal{Z}$$

$$e^{-\frac{{\delta }_{\tau }}{2}\sum _{x,y}{V}_{xy}{\hat{\rho }}_x{\hat{\rho }}_y}\cong \int D\varphi e^{-\frac{1}{2{\delta }_{\tau }}\sum _{x,y}{\varphi }_x{V}_{xy}^{-1}{\varphi }_y}{e}^{i\sum _x{\varphi }_x{\hat{\rho }}_x}.$$

$$e^{\frac{{\delta }_{\tau }}{2}(1-\eta )\sum _x{V}_{xx}({\hat{\rho }}_x^{\prime }{)}^2}\cong \int D\chi e^{-\frac{1}{2{\delta }_{\tau }}\sum _x\frac{{\chi }_x^2}{(1-\eta )V_{xx}}}{e}^{\sum _x{\chi }_x{\hat{\rho }}_x^{\prime }}.$$

$$\chi \to \chi -{\delta }_{\tau }{V}_{xx}(1-\eta ).$$

$$\mathrm{Tr}\left(e^{-{\hat{A}}_1}{e}^{-{\hat{A}}_2}\dots e^{-{\hat{A}}_n}\right)=$$

$$\displaystyle=\det\left(\begin{array}[]{cccc}1&-e^{-A_{1}}&0&\ldots\\ 0&1&-e^{-A_{2}}&\ldots\\ \vdots&&\ddots&\\ e^{-A_{n}}&0&\ldots&1\\ \end{array}\right),$$

$$\mathcal{Z}=\int D\Phi |\det M(\Phi ){|}^2{e}^{-S_{\eta }(\Phi )},$$

$$S_{\eta }(\Phi )=$$

$$+\sum _{x,t}\frac{({\chi }_{x,t}-(1-\eta ){\delta }_{\tau }{V}_{xx}{)}^2}{2(1-\eta ){\delta }_{\tau }{V}_{xx}}.$$

$$M(\Phi )=$$

$$\displaystyle=\left(\begin{array}[]{cccccc}1&-e^{-\delta_{\tau}h}&0&0&0&\ldots\\ 0&1&-e^{i\Phi_{1}}&0&0&\ldots\\ 0&0&1&-e^{-\delta_{\tau}h}&0&\ldots\\ 0&0&0&1&-e^{i\Phi_{2}}&\ldots\\ \vdots&&&&\ddots&\\ e^{i\Phi_{N_{\tau}}}&0&0&&\ldots&1\\ \end{array}\right),$$

$$S_{\text{eff.}}(\Phi )=S_{\eta }(\Phi )-\ln (\det M(\Phi )\det M(\Phi {)}^{†})$$

$$|\det M(\Phi ){|}^2=\int D\Psi {\Psi }^{†}{e}^{-{\Psi }^{†}(M{M}^{†}{)}^{-1}\Psi },$$

$$\frac{\partial \ln \det M}{\partial {\varphi }_{x,t}}=\mathrm{Tr}\left(M^{-1}\frac{\partial M}{\partial {\varphi }_{x,t}}\right).$$

$$M(\Phi )=$$

$$\displaystyle=\left(\begin{array}[]{cccccc}1&D_{1}&0&0&0&\ldots\\ 0&1&D_{2}&0&0&\ldots\\ 0&0&1&D_{3}&0&\ldots\\ 0&0&0&1&D_{4}&\ldots\\ \vdots&&&&\ddots&\\ D_{2N_{\tau}}&0&0&&\ldots&1\\ \end{array}\right),$$

$$M^{-1}(\Phi )=$$

$$\displaystyle=\left(\begin{array}[]{cccccc}g_{1}&\ldots&\ldots&\ldots&\ldots&\bar{g}_{2N_{\tau}}\\ \bar{g}_{1}&g_{2}&\ldots&\ldots&\ldots&\ldots\\ \ldots&\bar{g}_{2}&g_{3}&\ldots&\ldots&\ldots\\ \ldots&\ldots&\bar{g}_{3}&g_{4}&\ldots&\ldots\\ \vdots&&&&\ddots&\\ \ldots&\ldots&\ldots&&\ldots&g_{2N_{\tau}}\\ \end{array}\right).$$

$$g_i={\left(I+D_i{D}_{i+1}...D_{2N_{\tau }}{D}_1...D_{i-1}\right)}^{-1},$$

$$g_{i+1}=D_i^{-1}{g}_i{D}_i.$$

$${\bar{g}}_i=D_{i+1}...D_{2N_{\tau }}{g}_i,$$

$${\bar{g}}_{i+1}=D_{i+1}^{-1}{\bar{g}}_i{D}_i.$$

$${\bar{g}}_i=D_i^{-1}(I-g_i),$$

$${\bar{M}}^{(l)}{X}^{(l)}=Y^{(l)},$$

$$D_k^{(l+1)}$$

$$Y_k^{(l+1)}$$

$$X_{2k}^{(l)}$$

Full text

I Introduction

Graphene, with its orders of magnitude higher charge carrier mobility, is considered silicon’s ideal replacement for semiconductor-based devices. However, clean suspended graphene, which features the maximal carrier mobility at the same time lacks an energy gap in its band structure, the existence of which is the prerequisite for building graphene-based transistors.

Hypothetically such a gap could exist, since the small Fermi velocity of $v_{F}\approx c/300$ (where $c$ is the speed of light) leads to strong interactions between electronic quasi-particles with an effective fine-structure constant of $\alpha_{\textrm{eff}}=e^{2}/(\hbar v_{F})\approx 2.2$. Thus one expects, based on the analogy to chiral symmetry breaking in quantum field theories, numerous theoretical arguments Gamayun et al. (2010); Khveshchenko and Leal (2004); Araki and Hatsuda (2010); Araki (2011, 2012) and numerical simulations Drut and Lähde (2009a, b); Hands and Strouthos (2008); Armour et al. (2010); Buividovich and Polikarpov (2012), that for $\alpha_{\textrm{eff}}$ larger than some critical value $\alpha_{c}\approx 1$, interactions destabilize the system towards spontaneous formation of gapped ordered phases. Besides the well-known anti-ferromagnetic spin-density wave (SDW) order favored by sufficiently strong on-site repulsion Herbut (2006a); Herbut et al. (2009); Juričić et al. (2009); Semenoff (2012), various combinations of nearest neighbor and other short-range couplings might also induce such phases as a charge-density wave (CDW) Sorella and Tosatti (1992); Semenoff (1984); Herbut (2006b); Araki and Semenoff (2012); Gracey (2018), topological insulators Raghu et al. (2008), spontaneous Kekulé distortions Hou et al. (2007); Classen et al. (2017) as well as coupled spin-charge-density-wave phases Makogon et al. (2012) and spin spirals Peres et al. (2004). Also the existence of triple or multicritical points at which semimetal, CDW and SDW phases meet has been discussed Herbut (2006b); Classen et al. (2015, 2016).

Experimentally on the other hand, it has been firmly established that suspended graphene is a semimetal Elias et al. (2011); Mayorov et al. (2012). This implies that electronic two-body interactions are still too weak to induce a semimetal-insulator quantum phase transition. The absence of an energy gap has been reproduced in first-principle numerical simulations Ulybyshev et al. (2013); Smith and von Smekal (2014); Boyda et al. (2016) which properly take into account the screening of the bare Coulomb potential by electrons in lower $\sigma$-orbitals Wehling et al. (2011). This screening increases the critical coupling for a semimetal-insulator transition up to roughly $\alpha_{c}\approx 3.1$, which is noticeably higher than the effective coupling strength $\alpha_{\textrm{eff}}\approx 2.2$ in suspended graphene. Similarly, numerical calculations of the conductivity yielded a finite result almost equal to that of non-interacting graphene, implying the absence of a band gap Boyda et al. (2016).

Despite being in the weak-coupling gapless regime, suspended graphene may still be quite close to a semimetal-insulator transition. The knowledge of how close real graphene is to a phase transition might still help to describe the strong-coupling aspects of the many-body physics of this material. An obvious example of the relevance of the position and order of the closest phase transition in the weak-coupling regime are the convergence radius and rate of the perturbative expansion. On the experimental side, applying mechanical strain can move suspended graphene closer to the phase transition Tang et al. (2015); Xiao et al. (2017).

Guided by the results obtained within renormalization group techniques Herbut (2006a); Herbut et al. (2009); Juričić et al. (2009) and strong-coupling expansion Semenoff (2012), most numerical studies Ulybyshev et al. (2013); Smith and von Smekal (2014); Tang et al. (2015) have focused on detecting the onset of spin-density wave (SDW) order, which is expected to be a second-order phase transition in the universality class of the $N_{f}=2$ chiral Heisenberg Gross-Neveu model in three space-time dimensions.111The chiral Gross-Neveu universality class has been verified for the hexagonal Hubbard model with purely on-site interactions through numerous numerical studies. See e.g. Refs. Assaad and Herbut (2013); Otsuka et al. (2016); Parisen Toldin et al. (2015); Hohenadler et al. (2014). Within the perturbative renormalization-group analysis the robustness of this scenario is supported by the observation that the long-range Coulomb potential is a marginally irrelevant interaction Herbut (2006a); Herbut et al. (2009); Juričić et al. (2009), and thus the semimetal-insulator phase transition should be driven by on-site interactions.

Since the screening of the bare Coulomb potential by electrons in $\sigma$-orbitals mostly suppresses short-distance interactions and the long-range potential is only weakly affected Wehling et al. (2011); Tang et al. (2015), the long-range Coulomb interaction might still dominate the near-critical behavior Juričić et al. (2009); Gamayun et al. (2010). This might favor ordered phases other than the anti-ferromagnetic SDW phase, with the charge-density wave being the most likely candidate Sorella and Tosatti (1992); Semenoff (1984); Herbut (2006b); Araki and Semenoff (2012); Gracey (2018). An unbiased study from first principles of the competition between different ordered phases in the vicinity of suspended graphene, considered as a point in the space of all possible inter-electron interactions, is thus desirable.

By extension, the universality class (and even the order) of the possible semimetal-insulator transition also remains unclear. At present, the only first-principles calculations of critical exponents were carried out in the Dirac cone approximation Drut and Lähde (2009a, b); Hands and Strouthos (2008); Armour et al. (2010) and have not unambiguously settled the issue. The prediction of Gross-Neveu scaling relies upon an identification with the Hubbard model with on-site interactions only, which is a rather drastic modification.

A transition to a phase other than SDW could certainly imply different critical properties: It has been argued, for instance, that a semimetal-CDW transition should fall into the chiral Ising universality class Janssen and Herbut (2014). Another interesting possibility is that the scale-invariant $\sim 1/r$ Coulomb interaction induces the so-called conformal phase transition (CPT) of infinite order Kaplan et al. (2009), at which physical observables exhibit exponential (“Miransky”) rather than powerlaw scaling Miransky and Yamawaki (1997). CPT generalizes the concept of the Berezinskii-Kosterlitz-Thouless transition Berezinskii (1971); Kosterlitz and Thouless (1973) to higher than two dimensions where long-range order is possible. It is a continuous transition characterized by an exponentially increasing correlation length and occurs when changes of some control parameter cause the merging of infrared and ultraviolet renormalization group fixed points which marks the transition from a conformal to a non-conformal phase Kaplan et al. (2009). For graphene modelled as $2+1$-dimensional Dirac fermions with bare Coulomb interaction, a CPT is predicted by the analysis based on Schwinger-Dyson equations Gorbar et al. (2002); Khveshchenko and Leal (2004).

In this work, we use our state of the art Hybrid Monte-Carlo simulation code with numerous improvements discussed recently in Ref. Buividovich et al. (2018) to check how close suspended graphene might be to a semimetal-insulator transition, and to address the properties of this transition. We use a realistic partially screened Coulomb potential Smith and von Smekal (2014) which accounts for screening by electrons in lower $\sigma$-orbitals Wehling et al. (2011) and drive the system towards the transition by uniformly rescaling this potential, as in Ref. Ulybyshev et al. (2013). We thus improve and re-check the results of previous studies Ulybyshev et al. (2013); Smith and von Smekal (2014); Tang et al. (2015) which might have significant systematic errors due to small lattice sizes, large discretization artifacts, and artificial mass terms. The most important improvements include:

• •

Unlike in Refs. Ulybyshev et al. (2013); Smith and von Smekal (2014) we study the competition of SDW and CDW phases in a completely unbiased way, without symmetry-breaking mass terms that favor a specific phase. To this end we use quadratic observables as order parameters, such as the squared spin or charge per sublattice Buividovich et al. (2018), which unlike the corresponding condensates are non-zero in finite volume even without external sources and allow for an unambiguous determination of the ground state.

• •

To avoid the loss of ergodicity of HMC simulations caused by zero modes of the tight-binding Hamiltonian without mass terms, we represent inter-electron interactions in terms of complex-valued Hubbard-Stratonovich fields Beyl et al. (2017); Ulybyshev and Valgushev (2017); Buividovich et al. (2018).

• •

An improved fermionic lattice action with exact sublattice (chiral) symmetry and strongly suppressed discretization errors Buividovich et al. (2018, 2016).

• •

An efficient non-iterative Schur complement solver which significantly speeds up the simulations Ulybyshev et al. (2018a).

• •

Molecular dynamics trajectories which use exact fermionic forces and avoid the use of pseudofermions, which leads to another sizable performance improvement. This is a very recent development and is described in detail in Section II.2.

• •

Several improvements of the simulation parameters such as lower electronic temperatures ($T=0.125\,{\rm eV}$ instead of $0.5\,{\rm eV}$), larger spatial lattice sizes (up to $L=24$ instead of $18$) and finer discretization of the Euclidean time axis ($\delta_{\tau}=0.0625\,{\rm eV}^{-1}$ instead of $0.1\,{\rm eV}^{-1}$).

Using infinite-volume extrapolations of order parameters, we are able to determine that SDW order spontaneously forms, without being favored by a source term, at a critical coupling of $\alpha_{c}\approx 3.5$. This is larger than the previous estimate $\alpha_{c}\approx 3.1$ and implies that the scenario of suspended graphene being in the semimetal phase remains stable under our improvements of the HMC method. With high confidence we also rule out the presence of CDW order for couplings of $\alpha_{\textrm{eff}}\lesssim 5$.

Furthermore, we address the question of the universal properties of the phase transition induced by rescaling the screened inter-electron interaction potential in suspended graphene Wehling et al. (2011). Quite intriguingly, we find indications that the ratio between on-site and non-local interactions in the screened potential might favour the infinite-order phase transition scenario predicted in Refs. Gorbar et al. (2002); Khveshchenko and Leal (2004).

By studying the finite-size scaling of the squared spin per sublattice, we find that the ratio of critical exponents $\beta/\nu$ is close to exactly one in good approximation for the semimetal-SDW transition. Furthermore, the collapse of data points from different lattice sizes onto a universal scaling function is rather insensitive toward the choice of correlation length exponent $\nu$, with the optimal choice drifting slightly towards larger values when smaller lattices are excluded from the analysis. We obtain evidence that a collapse may occur naturally in infinite volume, without the need for a rescaling factor $L^{1/\nu}$, which formally corresponds to the limit $\nu\to\infty$. We argue that such behavior is consistent with a phase transition governed by Miransky scaling, with finite-volume corrections which mimic a second order phase transition on small systems. We also compare our data with reference data obtained for the Hubbard model with purely on-site interactions, and conclude that the interpretation of the numerical results in terms of Gross-Neveu scaling is much less convincing for our non-local interaction potential than for purely on-site interactions.

II Simulation setup

II.1 The path-integral formulation of the partition function

The basic idea of first-principle Monte-Carlo simulations is to carry out a stochastic integration of the functional integral representation of the grand-canonical partition function ${\cal Z}={\rm Tr}\,e^{-\beta\hat{\mathcal{H}}}$, in which operators are replaced by fields, by using a Markov process which evolves the fields in computer time such that their time histories are consistent with the equilibrium distribution. Thermodynamic expectation values of observables $\langle\hat{O}\rangle=\frac{1}{{\cal Z}}{\rm Tr}\,(\hat{O}e^{-\beta\hat{\mathcal{H}}})$ are then obtained from measurements on a representative set of field configurations.

Our starting point is the interacting tight-binding Hamiltonian on the hexagonal lattice, in second quantized form:

[TABLE]

Here $\kappa$ is the hopping parameter, $\langle x,y\rangle$ denotes nearest-neighbor sites, $\sigma=\uparrow,\downarrow$ labels spin components and $\hat{\rho}_{x}=\hat{c}^{\dagger}_{x,\uparrow}\hat{c}_{x,\uparrow}+\hat{c}^{\dagger}_{x,\downarrow}\hat{c}_{x,\downarrow}-1$ is the electric charge operator. The creation- and annihilation operators satisfy the anticommutation relations $\{\hat{c}_{x,\sigma},\hat{c}^{\dagger}_{y,\sigma^{\prime}}\}=\delta_{x,y}\delta_{\sigma,\sigma^{\prime}}$. $V_{xy}$ is the partially screened Coulomb potential used in Refs. Smith and von Smekal (2014); Körner et al. (2017). To drive the system towards the semimetal-insulator phase transition, we rescale this potential by a factor $\lambda>1$, so that suspended graphene corresponds to $\lambda=1$. As the cRPA data of Tang et al. (2015) suggests, for not very large distances of order of few lattice spacings, the effect of strain can be roughly described in terms of such a rescaling.

The potential $V_{xy}$ contains the exact values obtained from calculations within a constrained random-phase approximation (cRPA) by Wehling et al. in Ref. Wehling et al. (2011) for the on-site $U_{00}$, nearest-neighbor $U_{01}$, next-nearest-neighbor $U_{02}$, and third-nearest-neighbor $U_{03}$ interaction parameters.222There is still some minor disagreement over the exact values of these parameters in graphene Tang et al. (2015). The uncertainties are most likely too small to have any significant effect on the results of this work however. At longer distances a momentum dependent phenomenological dielectric screening formula, derived also in Ref. Wehling et al. (2011) based on a thin-film model, is used for a smooth interpolation to an unscreened Coulomb tail. Both results are combined via a parametrization based on a distance dependent Debye mass $m_{D}$. The matrix elements $V_{xy}$ are then filled using

[TABLE]

where $a$ is the nearest-neighbor distance and $m_{D}$, $m_{0}$, $c$ and $\gamma$ are piecewise constant chosen such that $m_{D},m_{0}\to 0$ and $c,\gamma\to 1$ for $r\gg a$. Tables with precise values of these parameters can be found in Ref. Smith and von Smekal (2014). Fig. 1 shows this interaction potential in comparison to the unscreened Coulomb potential.

To avoid a fermion sign problem (where the measure of the functional integral becomes complex or of indefinite sign, which prevents importance sampling), we apply the following canonical transformation to the Hamiltonian (1): Hole creation and annihilation operators $\hat{b}^{\dagger}_{x},\hat{b}_{x}$ are introduced for the spin-down electrons and the sign of these is then flipped on one of the triagonal sublattices of the hexagonal lattice. The transformation law can be summarized as

[TABLE]

where the signs in the second line alternate between the two sublattices. This leads to $\hat{\rho}_{x}=\hat{a}^{\dagger}_{x}\hat{a}_{x}-\hat{b}^{\dagger}_{x}\hat{b}_{x}$. We also apply the following Fierz transformation to the on-site interaction term:

[TABLE]

Here $\hat{\rho}_{x}^{\prime}=\hat{a}^{\dagger}_{x}\hat{a}_{x}+\hat{b}^{\dagger}_{x}\hat{b}_{x}$ is the spin-density operator and the constant $\eta$ can be chosen in the range $(0,1)$. The purpose of this transformation is to extend the Hubbard fields (introduced below) to complex numbers. This is necessary when the Hamiltonian (1) contains no mass terms, as the configuration spaces of both purely real and purely imaginary auxiliary fields then form disconnected regions, separated by infinitely high potential barriers (extended manifolds where the fermion determinant vanishes). The additional degrees of freedom of complex fields allow our Monte-Carlo algorithm to circumvent the barriers and ensure ergodicity Beyl et al. (2017); Ulybyshev and Valgushev (2017); Buividovich et al. (2018). The constant $\eta$ interpolates between real and imaginary fields.

To derive the functional integral, we start with a symmetric Suzuki-Trotter decomposition which yields

[TABLE]

where the exponential is factorized into $N_{\tau}$ terms and the kinetic $\hat{\mathcal{H}}_{0}$ and interaction $\hat{\mathcal{H}}_{\text{int}}$ contributions are separated. This introduces a finite step size $\delta_{\tau}=\beta/N_{\tau}$ in Euclidean time and a discretization error $O(\delta^{2}_{\tau})$. The four-fermion terms appearing in $\hat{\mathcal{H}}_{\text{int}}$ are now converted into bilinears by Hubbard-Stratonovich (HS) transformation. We use two distinct variants: The first term $\sim\eta\hat{\rho}_{x}^{2}$ appearing on the right hand side of Eq. (II.1) is re-absored into the interaction matrix $V_{xy}$ appearing in (1) and the combined expression is then transformed using

[TABLE]

The $\sim(1-\eta)(\hat{\rho}_{x}^{\prime})^{2}$ term in Eq. (II.1) is transformed by its own, using

[TABLE]

In effect, we have introduced a complex bosonic auxiliary field $\Phi$ (“Hubbard field”) with real part $\chi$ and imaginary part $i\phi$. Note that the transformations are applied once to each timeslice, leading to $\phi\equiv\phi_{x,t}$ and $\chi\equiv\chi_{x,t}$. The third term in Eq. (II.1) is already a bilinear and doesn’t need to be transformed. Due to the translational invariance of the integration measure in Eq. (7) it can be absorbed into the real part of the Hubbard field through the transformation

[TABLE]

To compute the trace in the fermionic Fock space (with anti-periodic boundary conditions) appearing in Eq. (5) we use

[TABLE]

where $\hat{A}_{k}=\left(A_{k}\right)_{ij}\hat{c}^{{\dagger}}_{i}\hat{c}_{j}$ are the fermionic bilinear operators and $A_{k}$ (without hat) contain matrix elements in the single-particle Hilbert space. This identity is derived in Refs. Hirsch (1985); Blankenbecler et al. (1981); Montvay and Muenster (1994). Applying (II.1) to Eq. (5), we obtain

[TABLE]

with

[TABLE]

Here $\widetilde{V}$ denotes a modified interaction matrix wherein the diagonal elements have been rescaled by a factor of $\eta$ by the Fierz transformation (II.1). The constant shift of $\chi$ in the second sum results from Eq. (8). The fermion matrix is given by

[TABLE]

where we use the short-hand notation $e^{i\Phi_{t}}\equiv\textrm{diag}\left(e^{\chi_{x,t}+i\phi_{x,t}}\right)$ and $h$ denotes the single-particle tight-binding hopping matrix. $|\det M(\phi)|^{2}$ appears in (14) since the fermionic matrices for spin-up and spin-down electrons are $M$ and $M^{\dagger}$, respectively.

Eq. (14) is exactly of the form required for Monte-Carlo simulations: ${\cal Z}$ is expressed as a functional integral over classical field variables with a positive-definite measure. We point our here that the appearence of matrix exponentials $e^{-\delta_{\tau}h}$ in the fermion matrix $M(\Phi)$ leads to a non-local fermion action. This action has an exact sublattice-particle-hole symmetry, even at finite $\delta_{\tau}$ and in the presence of the fluctuating Hubbard fields Buividovich et al. (2018, 2016). In contrast, the linearized action used in the previous studies Smith and von Smekal (2014); Körner et al. (2017); Ulybyshev et al. (2013) corresponds to expanding the blocks $e^{-\delta_{\tau}h}$ in $M(\Phi)$ to linear order in $\delta_{\tau}$. The main disadvantage of this linearized formulation is that the leading discretization errors generate a strong explicit breaking of the spin rotational symmetry in this case, which is only suppressed at very large $N_{\tau}$ Buividovich et al. (2018, 2016). $M(\Phi)$ is a dense matrix here, which makes iterative inversion methods such as the standard conjugate-gradient solver rather inefficient. $M(\Phi)$ can be efficiently inverted however, using a recently developed solver based on Schur decomposition Ulybyshev et al. (2018a).

II.2 Hybrid-Monte-Carlo with exact fermionic forces

This work employs the Hybrid-Monte-Carlo (HMC) algorithm based on the formalism originally developed in Refs. Brower et al. (2011, 2012) to study the graphene tight-binding model with interactions. HMC has its origins in lattice QCD simulations DeGrand and DeTar (2006); Montvay and Muenster (1994); Buividovich and Ulybyshev (2016) but is increasingly being applied also in condensed matter physics Hands and Strouthos (2008); Armour et al. (2010, 2011); Del Debbio and Hands (1996); Drut and Lähde (2009a, b, c, 2011); Brower et al. (2012); Ulybyshev et al. (2013); Smith and von Smekal (2014); Beyl et al. (2017); Ulybyshev and Katsnelson (2015); Ulybyshev et al. (2018b); Boyda et al. (2016); DeTar et al. (2017); Yamamoto and Kimura (2016); Luu and Lähde (2016); Berkowitz et al. (2018); Buividovich et al. (2017); Ulybyshev et al. (2017) alongside determinantal Quantum-Monte-Carlo simulations following Blankenbecler, Scalapino and Sugar (BSS-QMC) Blankenbecler et al. (1981); Scalettar et al. (1986). As we have described the individual steps of HMC in detail in several publications Smith and von Smekal (2014); Körner et al. (2017); Buividovich et al. (2018) we will focus entirely on a recent development here, whereby the algorithm is implemented with exact fermionic forces rather than using pseudofermions.

The HMC algorithm includes molecular dynamics (MD) trajectories, during which the Hubbard field is evolved in computer time by an artificial Hamiltonian process. During these trajectories the effective action

[TABLE]

plays the role of potential energy for the Hubbard field $\Phi$. Obviously, one needs to compute the derivative of the effective action with respect to Hubbard field in order to solve Hamilton’s equations. The standard approach is to use a stochastic representation of the determinants in Eq. (23):

[TABLE]

which introduces an additional pseudofermionic field $\Psi$. Calculations of derivatives of $S_{\textrm{eff.}}(\Phi)$ with respect to $\Phi$ then require just one solution of the linear equation $MM^{\dagger}\Psi=X$, where $X$ is a Gaussian distributed field. This solution can be obtained using an iterative solver or a non-iterative solver Ulybyshev et al. (2018a). The latter strategy was used in Ref. Buividovich et al. (2018).

However, one can go even further and avoid pseudofermions entirely, by computing the derivatives of $S_{\textrm{eff.}}(\Phi)$ directly, starting from Eq. (23). Calculations of derivatives of $S_{\eta}(\Phi)$ are trivial, while derivatives of $\ln(\det M(\Phi)\det M(\Phi)^{\dagger})$ can be computed using:

[TABLE]

It turns out that this requires the knowledge of only a few elements of the fermion propagator $M^{-1}$. Due to the special band structure of the matrix $M$ given by Eq. (II.1), we need only elements of $M^{-1}$ which are located in blocks immediately off the main diagonal.

To proceed, let us write the fermionic operator (II.1) in the general form:

[TABLE]

where even blocks $D_{2k}$ with $k=1...N_{\tau}$ correspond to diagonal matrices containing the exponentials $\pm e^{i\phi_{x,t}+\chi_{x,t}}$, and odd blocks are equal to exponentials of the tight-binding Hamiltonian $-e^{-\delta_{\tau}h}$. The inverse fermionic matrix can also be written in terms of spatial blocks:

[TABLE]

The matrix $M^{-1}$ is dense, but here we explicitly show only those blocks which are needed for our calculations. In fact, in the trace in Eq. (25) only the even blocks $\bar{g}_{2k}$ for all $k=1...N_{\tau}$ will contribute to the exact derivatives for computing the fermionic force.

We can now use part of the BSS-QMC algorithm Scalettar et al. (1986) to compute the desired blocks of the propagator. Due to the structure of $M$, the diagonal blocks of $M^{-1}$ can be formally written as

[TABLE]

and the following iteration formula can be proven:

[TABLE]

Analogously, off-diagonal blocks of $M^{-1}$ can be written as

[TABLE]

which leads to the relation:

[TABLE]

We can now either directly use Eq. (43) to obtain the $\bar{g}_{i}$ or first obtain the $g_{i}$ and use the relation

[TABLE]

between diagonal and off-diagonal blocks. By iterating either (43) or (41) we can easily find all elements of $M^{-1}$ needed for the computation of the derivative, starting from just one block, which is computed from scratch using the Schur complement solver Ulybyshev et al. (2018a). This is done by applying the solver to point sources in the corresponding time slice.

An important point here is that the whole procedure scales as ${N_{S}}^{3}N_{\tau}$, where $N_{S}$ is the number of sites in one Euclidean time slice of the lattice, so the scaling is not worse than that of the Schur complement solver itself. In practice however, the iterations (43) and (41) suffer from the accumulation of round-off errors, which limits the number of times they can be applied (this number depends mostly on the condition number of $e^{-\delta_{\tau}h}$). Afterwards, the block of $M^{-1}$ in the subsequent time slice must be computed from scratch. This is the so-called stabilization which is routinely used in BSS-QMC Bercx et al. (2017).

Finally, an additional simplification comes from the fact that we do not even need the full Schur complement solver for the computation of the blocks $g_{i}$ or $\bar{g}_{i}$. In order to demonstrate this, we sketch the essential parts of the solver. A more detailed description can be found in Ref. Ulybyshev et al. (2018a).

Essentially, the solver consists of tree stages. In the first stage we decrease the size of the linear system in an iterative procedure. At each step, the system has the form

[TABLE]

where $l$ denotes the step number. We start from the initial system with the matrix $M^{(0)}=M$, the unknown vector $X^{(0)}$ containing elements of the fermionic propagator, and a point source vector $Y^{(0)}$. In the simplest case, when $N_{\tau}$ is some power of 2, the size of the system decreases as $\bar{N}_{\tau}^{(l)}=N_{\tau}/2^{l-1}$. The general case is only slightly more complicated and described in Ref. Ulybyshev et al. (2018a).

The matrix $M^{(l)}$ always has the same form, with unit matrices in the diagonal blocks and with off-diagonal blocks $D^{(l)}_{k}$ for $k=1...\bar{N}_{\tau}^{(l)}$. Iterations are described by the relations

[TABLE]

for matrices and

[TABLE]

for vectors. $Y^{(l)}_{k}$ denotes the $k$-th timeslice of the vector $Y^{(l)}$.

The second stage is LU decomposition and solution of the compactified system at $l=l_{max}$. Thus we know the vector $X^{(l_{max})}$. Finally, the third stage is the reversed iterative process of reconstruction of the initial solution starting from $X^{(l_{max})}$, using matrix blocks $D^{(l)}_{k}$ and vectors $Y^{(l)}_{k}$ computed during the first stage:

[TABLE]

In the end, we arrive at the initial vector $X^{(0)}$ which gives us the matrix elements of the fermionic propagator.

One should note that the initial vector $Y^{(0)}$ contains non-zero elements only in one time slice. Due to the structure of the iterations (47), this feature is preserved at each step, thus we actually do not need to make the full loop over $k$ in (47). The same is true for backward iterations (48), for a different reason: we need only one time slice in the final solution $X^{(0)}$, since we are interested either only in diagonal blocks $g_{k}$ or only in off-diagonal blocks $\bar{g}_{k}$. Due to this simplification we need only one matrix-vector operation for each of the few time slices in which we actually recompute the elements of fermionic propagator from scratch. Thus the numerical cost of the method is dominated by matrix-matrix operations (46) and (43). This means that the number of floating-point operations scales as ${N_{S}}^{3}N_{\tau}$ with possible logarithmic corrections $\sim\log{N_{\tau}}$ from the sparse LU decomposition. Such a mild scaling with $N_{\tau}$ allows us to enlarge the Euclidean time extent of the lattice and work in the regime where systematic errors produced by the Trotter decomposition are negligible. In terms of scaling with $N_{\tau}$ at fixed $N_{S}$, the Schur complement solver definitely outperforms the Conjugate Gradient solver, see e.g. Fig. 2 in Ulybyshev et al. (2018a).

Generally, as it was shown in Ref. Ulybyshev et al. (2018a), the Schur complement solver is faster than preconditioned conjugate gradient for moderate lattice sizes up to $N_{S}=10^{3}...10^{4}$, depending on the model. This is one source of speedup. But an even more important source of speedup is that we can typically increase the integrator stepsize in MD trajectories by at least a factor of $50$ without losing the acceptance rate, if exact fermionic forces are used. The reason is a much smoother profile of fermionic forces in this case. A comparison of algorithm with pseudofermions (24) and exact (25) force calculations is shown in Fig. 2. For these tests, it was possible to achieve an acceptance rate of $\sim 0.7$ with exact fermionic forces with an integrator stepsize of $0.2$. Conventional HMC using stochastic representation of determinant (24) could achieve the same acceptance rate only with the stepsize $0.0032$. In this case we could decrease the number of steps in MD trajectories by a factor of $60$.

The actual speedup in terms of computer time is approximately half as much, since the iterative computation of the fermionic propagators (43) makes each integrator step twice as expensive.

II.3 Observables

SDW and CDW phases are characterized respectively by the separation of spin and charge between the two triangular sublattices. To study the competition between them in an unbiased way, we introduce order parameters which develop a non-zero expectation value in a finite volume even without any external sources. We use the square of charge and square of spin per sublattice, which are given by

[TABLE]

and

[TABLE]

where $L$ is the linear lattice size and

[TABLE]

As the sublattices $A$ and $B$ are equivalent, contributions from both are added to improve the signal-to-noise ratio.

A non-zero value of $\langle S^{2}\rangle$ in infinite volume does not unambiguously signal SDW order, as the same observable becomes finite in a ferromagnetic phase. To rule out ferromagnetic order, we also compute the mean squared magnetization

[TABLE]

See Appendix B of Ref. Buividovich et al. (2018) for expressions for $\langle S^{2}\rangle$, $\langle q^{2}\rangle$ and $\langle m^{2}\rangle$ in terms of fermionic Green functions.

II.4 Simulation parameters and data analysis

Using HMC, we simulate graphene sheets with an equal number $L$ of unit cells along each of the crystallographic axes. We simulate lattices with $L=6,12,18,24$. We impose periodic boundary conditions across the borders of rectangular sectors. This choice of geometry corresponds to that used in Ref. Smith and von Smekal (2014) but differs from the Born-von Kármán boundary conditions used in Körner et al. (2017). All results were obtained at temperatures $T=0.125\textrm{eV}=0.046\kappa$ with $N_{\tau}=128$, which leads to a time discretization $\delta_{\tau}=0.168\kappa^{-1}$. This choice is justified by the study of discretization effects in the exponential fermion matrix (II.1) made in Buividovich et al. (2016), where it was shown that the value of squared spin per sublattice (II.3) already stabilizes at this $\delta_{\tau}$. Thus, we can skip a rather expensive study of the $\delta_{\tau}\rightarrow 0$ limit. We stress here that a similar conclusion does not necessarily follow for other observables, so a convergence with respect to $\delta_{\tau}$ should carefully be checked in all future work. We choose $\eta=0.9$ as the mixing factor between real and imaginary parts of the Hubbard field introduced in Eq. (II.1), which is sufficient to ensure ergodic trajectories Buividovich et al. (2018); Ulybyshev and Valgushev (2017). For each produced lattice configuration we compute the full fermionic equal-time Green function $g(x,y)=\langle\hat{a}_{x}\hat{a}^{\dagger}_{y}\rangle=M^{-1}_{x,t,y,t}$, which is then used to compute all observables. To account for possible autocorrelation effects in our data, we use binning to calculate statistical errors. Typical sample sizes are on the order of several hundreds of independent measurements for a fixed set of parameters.

III Results

We begin with an unbiased study of the competition between CDW and SDW order along the lines of what was recently done for the extended Hubbard model with onsite and nearest-neighbor interactions Buividovich et al. (2018). We compute $\sqrt{\langle S^{2}\rangle}$ and $\sqrt{\langle q^{2}\rangle}$ for values of $\lambda$ in the range $[1.45\ldots 1.8]$. We use the square roots of (49) and (II.3) here, as they are characterized by an approach to the infinite volume limit which is linear in $1/L$ to good approximation, which is convenient for extrapolations to the thermodynamic limit. We use linear fits of the form $f(1/L)=a+b\cdot(1/L)$ to lattice sizes $L=12,18,24$ to carry out an $L\to\infty$ extrapolation. The approach to infinite volume is demonstrated in Fig. 3, while the final results is shown in Fig. 4.

From these results we can immediately conclude that SDW order is favored over CDW order: while the extrapolation of $\sqrt{\langle q^{2}\rangle}$ is consistent with zero for any of the coupling strengths considered, $\sqrt{\langle S^{2}\rangle}$ develops a non-zero expectation value around $\lambda_{c}\approx 1.65$. A more precise estimate of $\lambda_{c}$ will be given below. To rule out a ferromagnetic phase we also measure $\sqrt{\langle m^{2}\rangle}$. We find that it is smaller than $\sqrt{\langle S^{2}\rangle}$ by an order of magnitude for each parameter set and extrapolates to zero within errors for all cases as well. See Fig. 5 for an illustration.

To investigate the universal properties of the semimetal-SDW transition, we study the critical scaling of $\langle S^{2}\rangle$, which under the assumption of a second-order phase transition should respect

[TABLE]

where $f(x)$ is a universal finite-size scaling function and $x=L^{1/\nu}\left(\lambda-\lambda_{c}\right)/\lambda_{c}$ is the finite-size scaling parameter. Assuming naively that the transition is indeed of second order, we will use Eq. (55) to estimate values for $\lambda_{c}$, the ratio $\beta/\nu$, as well as for the correlation length exponent $\nu$ itself.

A priori the most obvious candidate for the universality class is the $D=3$, $N_{f}=2$ chiral Heisenberg Gross-Neveu model, based on the hypothesis that the main driving force of the phase transition are onsite interactions (as suggested by RG studies Herbut (2006a); Herbut et al. (2009); Juričić et al. (2009); Semenoff (2012)). While critical exponents for this class have been obtained in various ways (see e.g. Table I in Ref. Buividovich et al. (2018) for a summary), it is useful to also examine the Hubbard model (with onsite potential $U$ only) directly here, and obtain a data set which can be used as a point of reference in a one-to-one comparison. This will be the first step of our analysis. In principle the required simulations could be carried out using our HMC code, but for practical reasons333For the pure Hubbard model BSS is still faster by an order of magnitude than HMC, and remains the method of choice for simulations of contact interactions. HMC is advantageous for long-range interactions, as the number of auxiliary fields is independent of the choice of potential, whereas each interaction term requires an additional field in BSS-QMC.; For a publicly available BSS-QMC code see https://git.physik.uni-wuerzburg.de/ALF. we choose to produce this data using a GPU implementation of BSS determinantal Quantum-Monte-Carlo Blankenbecler et al. (1981) instead, which we will not describe here as the method is widely known (interested readers are referred e.g. to Refs. Gerlach (2017); Bercx et al. (2017)).

Using BSS to simulate the Hubbard model, we obtain data for $\langle S^{2}\rangle$ for several values $U$ around the phase transition, which has been estimated to occur at $U_{c}/\kappa\approx 3.8$ Assaad and Herbut (2013), with all other external parameters matching those of our graphene simulations. To extract both $U_{c}$ and $\beta/\nu$ from the data we determine a choice of $\beta/\nu$ for which the $U$-dependent curves of $\langle S^{2}\rangle L^{2\beta/\nu}$ for different $L$ intersect in one point. This is done by fitting the data from $L=6,12,18,24$ with linear functions close to the presumed transition point and adjusting $\beta/\nu$ in steps of $0.001$ until the enclosed triangles of all subsets of three lattices are minimized. Note that this requires choosing smaller fit windows on larger lattices, as the curvature of the order paramter grows with system size. Through this method, we find $\beta/\nu\approx 0.8$ and $U_{c}/\kappa\approx 3.9$ consistent with the previous measurements. After fixing the value of $\beta/\nu$ we then obtain $\nu$ by optimizing the collapse of $\langle S^{2}\rangle L^{2\beta/\nu}$ onto a universal scaling function, by fitting all data points from $L=6,12,18,24$ with a single polynomial function of $x=L^{1/\nu}(U-U_{c})/U_{c}$ (we find that we must use a polynomial of third order) and adjusting $\nu$, also in steps of $0.001$, until the $\chi^{2}$ per degree of freedom becomes minimal. With $\beta/\nu=0.812$ and $U_{c}/\kappa=3.942$ this yields $\nu=0.928$.444Note that the three digits quoted for all results here reflect the resolution used in our optimization procedure and by no means imply a corresponding accuracy. As a cross-check $U_{c}$ is also allowed to shift, yielding an optimal value of $U_{c}/\kappa=3.944$ for the data collapse. Fig. 6 summarizes these results.

A few comments are in order here: In principle our predictions for $U_{c}$ and $\beta/\nu$ should be affected by a systematic uncertainty due to a sensitivity to the windows in which linear fits are applied. We find however that these values are remarkably stable under variations of the fit windows and, quite conservatively, estimate these errors to be $\sim 1\%$ for $U_{c}$ and $\sim 2\%$ for $\beta/\nu$, which is also most likely larger than our statistical errors. Our results are quite close to the values $\beta/\nu\approx 0.9$ and $U_{c}/\kappa\approx 3.8$ quoted in Ref. Assaad and Herbut (2013), with the difference likely being due to finite size and temperature effects, which we suspect are the leading source of errors in our case. Likewise, finite size is likely the leading source of uncertainty for $\nu$. If we exclude the $L=6$ lattice and both the $L=6,12$ lattices from the optimized collapse we obtain $\nu=1.037$ and $\nu=1.024$ respectively, suggesting a combined finite-size and statistical error of at least $\sim 6\%$ (note that $\nu$ is slightly larger when only $L=6$ is excluded, suggesting that the we are already close to the thermodynamic limit). We point out here that our values are very much in line with those typically seen in in Monte-Carlo simulations of Hubbard-type models believed to fall into the Gross-Neveu universality class Otsuka et al. (2016); Parisen Toldin et al. (2015); Hohenadler et al. (2014); Assaad and Herbut (2013). Renormalization group studies tend to observe slightly larger values, whereby results as large as $\beta/\nu\approx 1.0$ and $\nu\approx 1.2$ have been predicted Gracey (2018); Knorr (2018); Zerf et al. (2017).

We now turn to the data generated with the realistic potential of graphene. The next logical step is to determine whether we can characterize the critical properties using the same exponents as for the Hubbard model. Thus, we fix $\beta/\nu=0.812$ and $\nu=0.928$ and test both the intersection of $\langle S^{2}\rangle L^{2\beta/\nu}$ for different $L$ and the collapse of data onto a universal function $f(x)$. We find that we can clearly rule out this possibility: Not only do $\langle S^{2}\rangle L^{2\beta/\nu}$ intersect nowhere in the region $\lambda\leq 2.0$, contradicting the results shown in Fig. 4, the quality of collapse is also very poor and leads to an unreasonably large estimate of $\lambda_{c}$. The best possible collapse is shown in Fig. 7 while we refrain from even showing any figures for the intersection.

In order to obtain first-principles estimates of $\lambda_{c}$, $\beta/\nu$ and $\nu$ for graphene we now repeat the same steps as for the Hubbard model: To extract $\lambda_{c}$ and $\beta/\nu$ from the data we determine a choice of $\beta/\nu$ for which linear fits to $\langle S^{2}\rangle L^{2\beta/\nu}$ for different $L$ intersect in one point. Using the resulting $\beta/\nu$ we then determine $\nu$ by optimizing the collapse. In doing so we observe a somewhat odd behavior: It appears that the data, while exhibiting just as small statistical errors for each data point, don’t constrain $\beta/\nu$ and $\lambda_{c}$ nearly as strongly as for the Hubbard model. By choosing different fit windows we can obtain values of $\beta/\nu$ in the range $[0.95,1.0]$ (with the different lines crossing in one point in each case), while $\lambda_{c}$ shifts in the range $[1.62,1.7]$. The left panel of Fig. 8 represents our best possible choice of fit windows (with smallest $\chi^{2}$ for the fits) and leads to $\lambda_{c}\approx 1.70$ and $\beta/\nu=0.967$. Compared to Fig. 4 this estimate for $\lambda_{c}$ seems slightly too large.

When determining $\nu$ a similarly peculiar observation is made: While an optimized collapse of the $L=6,12,18,24$ data onto a universal function $f(x)$, $x=L^{1/\nu}\,(\lambda-\lambda_{c})/\lambda_{c}$ yields $\nu=1.473$ (right panel of Fig. 8), we find that the optimal value of $\nu$ increases slightly when the $L=6$ lattice is excluded from the fit, resulting in the optimal choice of $\nu=1.635$ for fits to $L=12,18,24$ and $\nu=1.602$ for $L=18,24$ (see Fig. 9). While the spread of these numbers suggests that the combined statistical and finite-size error can be as large as $\sim 10\%$, the observation that the optimal value moves further away from that of the Hubbard model on larger lattices is unexpected. We point out here that it was sufficient to model the universal scaling function with a quadratic polynomial for the case of graphene.

One could still be inclined to interpret these results in terms of Gross-Neveu scaling: As our value of $\beta/\nu$ for graphene is closer to the RG result ($\beta/\nu\approx 1.0$) than for the Hubbard model, one might speculate that non-universal corrections to scaling are more severe in the case of a pure onsite potential for some reason, that RG exponents are closer to the true values than MC predictions, and that our results with the long-range potential reflect the true universal behavior more closely. Without commenting on the plausiblity of this scenario we note that our smallest prediction for $\nu$ is still $\sim 20\%$ larger than the RG result ($\nu\approx 1.2$), however.

To get a sense of how significant this deviation is (in principle a visually only slightly less convincing collapse can be obtained with $\nu=1.2$, see Fig. 10) we do the following: After fixing $\beta/\nu=0.967$ and $\lambda_{c}=1.7$ we shift $\nu$ in the range $[0.5,5.0]$ and obtain the $\chi^{2}$ per degree of freedom resulting from fits of quadratic polynomials to $\langle S^{2}\rangle L^{2\beta/\nu}$ as a function of $x=L^{1/\nu}(\lambda-\lambda_{c})/\lambda$ for each choice of $\nu$. We do this for the full set of lattice sizes $L=6,12,18,24$ and for the cases where $L=6$ and $L=6,12$ are ignored. A similar procedure is repeated for the Hubbard model (with $\beta/\nu=0.812$, $U_{c}=3.944$ and third order polynomials) for comparison. The results are shown in Fig. 11.

The first thing we can say is that, on a quantitative level, $\nu=1.2$ is clearly disfavored for graphene. While for the $L=6,12,18,24$ lattices the $\chi^{2}/\textrm{d.o.f}$ is about $1.7$ times larger at $\nu=1.2$ compared the optimum, the ratio grows even larger as small lattices are excluded, up to about $2.0$ on lattice sizes $L=18,24$.

An interesting observation is that the $\chi^{2}$ curves for both the Hubbard model and graphene become flatter at large $\nu$ as small lattice sizes are excluded. To some degree this is not surprising, as a smaller number of data points places weaker constraints on the parameters of the polynomials used to model the universal scaling function. What is striking however is the substantial difference between graphene and the Hubbard model. For the Hubbard model, the scale of the y-axis is larger by an order of magnitude, despite the fact that higher order polynomials were used, reflecting the fact that critical exponents are constrained much stronger by the data. For graphene on the other hand, it appears as if the curves very quickly converge to a completely flat curve with $\chi^{2}\sim 1$ for $\nu\gtrsim 2.0$ on larger volumes, suggesting that for large systems any sufficiently large choice of $\nu$ will work equally well. This is not at all what one expects for a regular second order phase transition.

To conclude our analysis, we conservatively state that the numerical data for the long-range interaction potential is quite different from the data obtained for purely on-site interaction, and is hardly consistent with the Gross-Neveu universality class. We conjecture that instead our results could be interpreted as signatures of Miransky scaling.

A CPT is known to occur in quantum electrodynamics (QED) with massless fermions, both in $3+1$ and $2+1$ dimensions, with the number of fermion flavors being the control parameter in the later case Gusynin et al. (1998); Appelquist et al. (1995, 1988); Braun et al. (2014).555QED2+1 has been considered as a model for a similar transition believed to occur in $SU(N)$ gauge theories (see e.g. Refs. Banks and Zaks (1982); Kusafuka and Terao (2011); Appelquist et al. (2008); Braun and Gies (2010); Gies and Jaeckel (2006) and references therein). In QED2+1 one speaks of a “pseudo-conformal” transition as conformal symmetry is explicitly broken by a dimensionful coupling constant but the theory nevertheless exhibits an effective low-energy scale invariance. The proper low-energy effective theory of graphene (“reduced QED4”), combines features of both, as electron motion is restricted to a plane but photons can propagate in the three dimensional bulk (one also takes the smallness of the Fermi velocity into account, which leads to effectively instantaneous Coulomb interactions with an exact conformal symmetry). In this theory, a chiral phase transition exhibiting Miransky scaling666A sufficiently strong four-fermion coupling (corresponding roughly to onsite interactions in the microscopic theory) can change the transition to one of second order Gamayun et al. (2010). has been demonstrated by solving Dyson-Schwinger equations Gorbar et al. (2002); Khveshchenko and Leal (2004). In Ref. Gamayun et al. (2010) it was argued that this CPT formally corresponds to the limit $\beta,\nu\to\infty$, $\delta=1$ of a second order transition and that the usual hyperscaling relations may apply. In our case $d=2$ (where $d$ is the number of spatial dimensions) and with $\delta=1$ the relation

[TABLE]

would thus lead to $\beta/\nu=1$, which agrees with our estimate for the optimal value at a level of about $3\%$ and is thus well within the present errors.

For QED2+1Gusynin and Reenders (2003); Goecke et al. (2009) and reduced QED4 Liu et al. (2009) it is well known that that the CPT exhibits a strong sensitivity to an infrared cutoff. For many-flavor QCD it was shown that Miransky scaling receives powerlaw corrections from an infrared RG fixed point of the gauge coupling Braun et al. (2011). It is thus reasonable to assume that finite-size effects mimic a second order phase transition, and that $\beta,\nu\to\infty$ as the thermodynamic limit is approached, with their ratio being fixed by the hyperscaling relation. The slight drift of the exponent $\nu$ towards larger values observed above could be interpreted in these terms as can be the relative insensitivity of the quality of collapse towards further increases of $\nu$ on large lattices.

As a final test we therefore study the finite-size scaling properties of $\langle S^{2}\rangle$ for graphene with fixed $\beta/\nu=1.0$. We find that linear fits to $\langle S^{2}\rangle L^{2\beta/\nu}$ with $L=6,12,18,24$ cross in one point to good precision if we choose the fit windows as $\lambda=[1.525;1.725]$ (see left panel of Fig. 12) and thereby estimate $\lambda_{c}\approx 1.62$, which is substantially lower than the result $\lambda_{c}\approx 1.70$ obtained with $\beta/\nu=0.967$ but falls much closer to the location one would expect, judging by the $L\to\infty$ extrapolated results for $\sqrt{\langle S^{2}\rangle}$ shown in Fig. 4: By fitting the Miransky scaling function

[TABLE]

as appropriate for reduced QED4 (see e.g. Ref. Gamayun et al. (2010)) to $\sqrt{\langle S^{2}\rangle}|_{L=\infty}$ (which works extremely well), we can independently estimate $\lambda_{c}=1.61\pm 0.02$ which agrees with our prediction from finite-size scaling to very good precision. This agreement can be seen as evidence in favor of the CPT scenario.

Perhaps more importantly, however, we note that the slopes in the intersection plot in the left panel of Fig. 12 do not appear to increase towards infinity with increasing volumes as they should for a second-order phase transition. In fact, the solid black line in this plot represents our infinite-volume extrapolation which is clearly not vertical. If true, there would certainly be no way that this could ever happen in a second-order phase transition for any finite value of $\nu$, i.e. without rescaling the reduced control parameter by $L^{1/\nu}$ along the abscissa.

To investigate this somewhat more carefully, we plot the inverse slopes of our linear fits to $\langle S^{2}\rangle L^{2\beta/\nu}$ from the intersection plots over $1/L$ in the right panel of Fig. 12:

While the Hubbard model data once again shows the expected behavior, the inverse slopes for graphene obtained from the left panel in this figure (with $\beta/\nu=1.0$) are well described by a linear model fit to

[TABLE]

This non-zero intercept $b$ then provides our best numerical evidence of a finite slope in the infinite-volume limit and hence of $\beta,\nu\to\infty$ as CPT characteristics.

As mentioned, this also implies that on sufficiently large lattices a collapse onto a universal finite-size scaling function occurs for $x=(\lambda-\lambda_{c})/\lambda_{c}$ alone, i.e. without rescaling the reduced coupling by the factor $L^{1/\nu}$, which one expects if $\nu\to\infty$, but which also explains the difficulty in determining a stable value for $\beta/\nu$ from intersection points on larger lattices where the curves do not intersect anymore when the data collapses all by itself in an ever growing region around $\lambda_{c}$.

For completeness we close this section with adding that a behavior as expected for a second-order phase transition, modelling the inverse slopes with $f(1/L)=a\cdot(1/L)^{1/\nu}$ without intercept, can also be used to fit the graphene data, resulting in $\nu=1.6\pm 0.1$, but yielding a $\chi^{2}$/d.o.f. which is larger by about a factor of six as compared to the linear fit described above.777Note that $\chi^{2}$ values are $<1$ for both models here, which reflects the fact that the slopes and their error bars are themselves the results of linear fits. That such second-order fits work as well, on sufficiently small lattices, might rather be a manifestation of the difficulty in distinguishing CPT behavior from second-order scaling in finite volumes.

The direct comparison between Hubbard model and graphene data, however, provides quite compelling evidence that the second-order scaling, which works very well for the former, gets increasingly difficult to maintain with increasing volumes in the case of graphene with long-range interactions where the CPT scenario appears to provide the much more natural explanation.

IV Conclusion and Outlook

In this work we carried out a detailed unbiased study of the competition of SDW and CDW orders in graphene, with a realistic two-body potential that includes an unscreened long-range Coulomb tail. Using state of the art Hybrid-Monte-Carlo simulations we were able to determine that the potential must be uniformly rescaled by roughly a factor of $\sim 1.6$ to induce a semimetal-SDW phase transition, corresponding to a critical value of effective fine-structure constant $\alpha_{c}\approx 3.5$. This is substantially larger than the value $\alpha_{c}\approx 3.1$ predicted by previous studies Ulybyshev et al. (2013); Smith and von Smekal (2014); Tang et al. (2015) which were much more strongly affected by discretization artifacts, finite-size effects and larger temperatures. We find no evidence for CDW order, confirming that SDW is the preferred phase as predicted by renormalization group analysis Herbut (2006a); Herbut et al. (2009); Juričić et al. (2009) and strong-coupling expansion Semenoff (2012).

A careful study of the critical properties suggested that the expected $D=3$, $N_{f}=2$ chiral Heisenberg Gross-Neveu universality class, as expected for the corresponding Hubbard model, is unlikely to apply to the semimetal-insulator transition in graphene with long-range Coulomb interactions. Our lower-bound estimates for the correlation-length exponent $\nu$ are significantly larger than the largest values predicted by RG studies for this class. A direct comparison with a data set produced for the Hubbard model with on-site interactions only also ruled out with high confidence that both systems can be described by a common set of critical exponents, demonstrating clearly that the long-range part of the inter-electron potential plays an important role in non-perturbative many-body physics of graphene.

In studying the finite-size scaling of the squared spin per sublattice an unexpected property of $\nu$ was observed: The optimal choice, which produces the best possible collapse of data from different lattice sizes onto a universal finite-size scaling function, seems to drift slightly towards larger values as smaller lattices are disregarded during the analysis, instead of converging towards the RG predictions as one would expect if the difference was due to corrections to scaling. Further investigations revealed that constraints on $\nu$ become weaker on large lattices, such that $\nu$ can be increased without affecting the quality of collapse substantially. This stands in stark contrast to the Hubbard model, where critical exponents are constrained much tighter by the data. Furthermore, we obtained evidence that a collapse may occur naturally in infinite volume, without the need for rescaling the reduced coupling by a factor $L^{1/\nu}$, which formally corresponds to the limit $\nu\to\infty$.

We have proposed that the observed behavior can be explained in terms of a conformal phase transition, exhibiting exponential “Miransky” scaling, which is predicted for $2+1$-dimensional Dirac fermions with bare Coulomb interactions by Dyson-Schwinger studies Gorbar et al. (2002); Khveshchenko and Leal (2004), and power-law corrections that mimic a second order phase transition caused by finite size effects Gusynin and Reenders (2003); Goecke et al. (2009); Liu et al. (2009); Braun et al. (2011). A formal hyperscaling relation between the exponents $\beta$ and $\nu$ seems to be fulfilled in good approximation.

Let us note that the phase transitions to CDW and SDW ordered states and an infinite-order conformal phase transition are basically the only theoretical predictions for graphene with long-range interactions of which we are aware, and our numerical analysis is intended to find the most likely scenario out of these three. This turns out to be the conformal phase transition scenario. It can be that our data is also consistent with some other exotic scenario which has not been studied so far.

On the technical side, we have described a variant of HMC which achieves a substantial performance improvement by using exact fermionic forces and avoiding the use of pseudofermions.

An obvious direction for future studies is to push the simulations towards even larger system sizes. The $L=24$ lattices studied in this work represent the largest systems which are feasible with our current computational resources. Repeating the finite-size scaling analysis with $L=24,30,36$, to test whether the trend of growing correlation length exponent $\nu$ continues, would be extremely beneficial and should become feasible in the near future. The infinite volume extrapolations shown in Figs. 3, 4 and 5 as of now use only three points and would also be greatly improved by including additional lattice sizes. Another possibility is to study in detail how the order of the phase transition depends on the balance of short- and long-range parts of the potential. This could guide experimental efforts to induce a conformal phase transition in real graphene samples through techniques such as mechanical strain Tang et al. (2015); Xiao et al. (2017).

Acknowledgements.

P.B. is supported by a Heisenberg Fellowship from the the Deutsche Forschungsgemeinschaft (DFG), grant BU 2626/3-1. M.U. is also supported by the DFG under grant AS120/14-1. D.S. and L.v.S.  are supported by the Helmholtz International Center for FAIR within the LOEWE initiative of the State of Hesse. Calculations were carried out on GPU clusters at the Universities of Giessen and Regensburg and on the JUWELS system at the Jülich Supercomputing Centre. We thank F. Assaad, Ch. Fischer and B.-J. Schaefer for helpful discussions.