Fermionic Quantum Critical Point of Spinless Fermions on a Honeycomb Lattice

Spinless fermions on a honeycomb lattice provide a minimal realization of lattice Dirac fermions. Repulsive interactions between nearest neighbors drive a quantum phase transition from a Dirac semimetal to a charge-density-wave state through a fermionic quantum critical point, where the coupling of Ising order parameter to the Dirac fermions at low energy drastically affects the quantum critical behavior. Encouraged by a recently discovery of absence of the fermion sign problem in this model, we study the fermionic quantum critical point using the continuous time quantum Monte Carlo method with worm sampling technique. We estimate the transition point $V/t= 1.356(1)$ with the critical exponents $\nu =0.80(3)$ and $\eta =0.302(7)$. Compatible results for the transition point are also obtained with infinite projected entangled-pair states.


Introduction
Interaction induced quantum phase transitions of Dirac fermions are of general interests in graphene [2], d-wave superconductors [3], topological insulators [4], ultracold atoms [5] and high energy physics [6]. One of the prototypical examples consists of halffilled spinless fermions on a honeycomb lattice interacting through nearest neighbor repulsionsĤ Eq. (1) is arguably the simplest model exhibiting a quantum phase transition of Dirac fermions in two dimension. However, despite its deceptively simple form, the model exhibits an unconventional quantum critical point which deserves detailed study because it may lay the foundation of understanding rich phenomena when other degrees of freedom or intertwined phases are involved. The phase diagram of Eq.(1) is easy to anticipate, see Fig.1. The system behaves like a classical lattice gas in the strong coupling limit (V t) and favors a staggered chargedensity-wave (CDW) state. The CDW state breaks the discrete sublattice symmetry and melts through a 2D Ising phase transition at finite temperature. In the weak coupling limit, quantum fluctuations due to fermion hopping destroy the CDW long range oder and restore the Dirac semimetal state. Since Dirac fermions are perturbatively stable against short range interactions, the quantum critical point separating the Dirac semimetal and the CDW state lies at a finite interaction strength V /t.
The topology of the phase diagram Fig.1 resembles that of the familiar 2D transverse field Ising model [3] where quantum fluctuations induced by the transverse field destroy the Ising long range order. However, in model Eq.(1) the coupling of the Ising order parameter to the Dirac fermions at low energy strongly affects its quantum critical behavior. It cannot be treated by the familiar scalar φ 4 -theory since integrating out the Dirac fermions will lead to a singular action for the Ising fields. The low energy physics is described by the Gross-Neveu-Yukawa theory [7,8,9,10] which features a fermonic quantum critical point. The Gross-Neveu-Yukawa theory has been studied intensively in the context of high energy physics [11,12,13] and quantum critical point scenario of the d-wave superconductors [14,15,16], however, there is no consensus concerning the critical exponents, partially due to uncontrolled approximations involved in various theoretic approaches.
Quantum Monte Carlo (QMC) simulations are valuable unbiased approach to study the quantum critical behavior if the notorious fermion sign problem is absent [17]. A recent example is the study of Dirac semimetal to antiferromagnetic insulator transition Figure 1. Schematic phase diagram of model Eq.(1). In strong coupling limit the system is in the charge-density-wave state. The long range order melts through a 2D Ising transition upon increase of temperature or a quantum phase transition upon decrease of V /t. The red dot represents a fermonic quantum critical point which is the focus of this paper.
in the half-filled repulsive Hubbard model on the honeycomb lattice [18,19,20,21]. Unfortunately, the seemingly simpler model Eq.(1) has a severe sign problem in the conventional auxiliary field QMC method [22]. Early QMC studies have thus been limited to high temperatures or small system sizes [23,24]. The meron-cluster algorithm [25] solves the sign problem for V ≥ 2t and simulations using it confirm the finite temperature Ising transition of several staggered fermion models [26,27]. However, the quantum critical point of the model Eq.(1) lies at V < 2t and is not accessible by the meron-cluster algorithm. The Fermi bag approach [28] has been used to study the 3D lattice massless Thirring model [29] and the Gross-Neveu model [30] with two flavors of four component Dirac fermions.
Recently, Ref. [1] discovered that the sign problem of the model Eq. (1) is absent in the continuous-time quantum Monte Carlo (CTQMC) formalism [31,32]. This allows us to access the quantum critical point in the CTQMC simulation. Using a standard finite size scaling analysis we estimate the critical point V /t ≈ 1.356 and the critical exponents ν ≈ 0.8, η ≈ 0.3. Our results are summarized in Table 1. We believe these results do not only apply to the specific microscopic model Eq.(1), but also hold for many intriguing problems including the chiral symmetry breaking of Dirac fermions [6] and the quantum critical point in the d-wave superconductors [3]. Future theoretical or experimental advances in either field [2,3,4,5,6] will be able to test our predictions.
The paper is organized as follows: In Sec.2.1 we briefly review the CTQMC method with focuses on the absence of the sign problem [1]. In Sec. 2.2 we introduce the worm sampling technique in the Monte Carlo calculation [33,34,32] to improve the efficiency of the simulation. Section 3 contains our main results and discussions as well as comparisons with results obtained by infinite projected entangled-pair states (iPEPS) calculations. In Sec.4 we briefly anticipate several future research directions based on this work. In the appendix we provide technique details of our numerical calculation (Appendix A) and additional results obtained on the π-flux square lattice (Appendix B).

Method
Two properties of the model (1) are essential for the absence of the sign problem in the CTQMC simulation. First, the filling is fixed at 1/2 per site because of the particle-hole symmetry of the model. Second, the hopping matrix defined in Eq.(2) satisfies where the "parity index" η i = 1(−1) for a site i ∈ A(B) sublattice.

Interaction expansion CTQMC and the sign problem
We expand the partition function of the system in terms of the interaction vertices of Eq.(3) [31,32] wheren i (τ ) = eĤ 0 τn i e −Ĥ 0 τ and Z 0 is the partition function of the noninteracting system.
. . . 0 = T Tr(e −βĤ 0 . . .)/Z 0 denotes the average over the noninteracting Hamiltonian Eq.(2) and T is the time ordering operator. The interaction vertices {i 1 , i 2 }, . . . , {i 2k−1 , i 2k } consist of k pairs of neighboring sites. The delta functions in the first line of Eq. (5) indicates that the interactions are instantaneous. G k is a 2k × 2k matrix where G 0 ij (τ ) = ĉ i (τ )ĉ † j 0 is the noninteracting Green's function. The particle-hole symmetry ensures that G 0 ii (0 + ) = 1/2 and therefore the diagonal element of G k vanishes. In addition, one has while using the anti-periodicity of the Green's function and Eq.(4) one has Equations (6)(7)(8) show that the Green's function matrix satisfies G k qp = −η ip G k pq η iq . Introducing a diagonal matrix D k = diag(η i 1 , η i 2 , . . . , η i 2k ), it can be written as In other word, the matrix G k D k is skew-symmetric and its determinant is non-negative because it equals the square of the Pfaffian of the matrix. We will see in a moment that this ensures the absence of a sign problem in the CTQMC simulation. We write Eq.(5) in a form suitable for Monte Carlo sampling where In the second line we have used det The absence of a sign problem allows us to simulate fairly large systems at low temperatures to access the quantum critical point. In this paper, we simulate clusters with L × L unit cells with periodic boundary conditions. The number of sites is N s = 2L 2 . Close to the quantum critical point, nonrelativistic corrections are irrelevant and the dynamical critical exponent z = 1 [10]. We thus scale the inverse temperature linearly with the system length β = 4L/3. Because of the β 3 scaling of the CTQMC algorithm [31,32], the largest system size L = 15 considered in this paper is smaller than the one used in the projective auxiliary field QMC studies of the Hubbard model [19,20,21].
To detect the onset of the CDW order, we measure the density-density correlation function where . . . = Tr(e −βĤ . . .)/Z denotes the average over the full Hamiltonian Eq.(1). The second summation in Eq.(12) runs over all sites j (in total N R of them) whose graph distance to the site i is R §. Two sites with even (odd) graph distance have the same (different) parities. The other two important observables are the square and quartic of the CDW order parameter ‡ The 1/k! factor has been canceled by the k! permutations of the vertices. § Strictly speaking, these sites may not be symmetrically related and may have slightly different correlation functions.
The Binder ratio [35] is calculated as:

Worm Algorithm
Measuring M 2 and M 4 using the conventional approach [31] requires explicit loops over the i, j, (k, l) indices and the measurements will dominate the runtime of Monte Carlo simulations. This is especially inefficient when noticing that each term in Eqs. (13)(14) may differ by orders of magnitude. To overcome this difficulty we extend the configuration space and use the worm algorithm [33,34,32] to sample them efficiently.
Notice the similarity between the partition function Eq.(5) and the observables Eqs.(13-14) where G k;ijτ extends G k to the following (2k + 2) × (2k + 2) matrix: Similar to Eq.(9), G k;ijτ satisfies the following equation s ZM 4 and enlarge the configuration space into Now the configurations C may contain a two site worm {i, j; τ } or a four site worm {i, j, k, l; τ } in addition to the vertices described in Eqs. (10). By sampling the extended configuration space we can treat the summation over i, j, k, l in Eq. (16)(17) and the summations over the vertices on an equal footing. Here ξ 2 and ξ 4 are two positive numbers we can choose freely to balance the configurations in different sectors. We have devised several Monte Carlo updates and describe them in Appendix A. We use the following notation to denote the relative time spend in the each sector [34] δ The observables (13-15) then are The density correlation function is measured when the configuration is in the W 2 space and the distance between the two worm sites i, j is equal to R, The weight of a configuration C ∈ W 2 with worm at {i, j; τ } is where we have used det(D k;ij ) = ( k =1 η i 2 −1 η i 2 )η i η j =(−1) k η i η j . One can similarly show that the weight of C ∈ W 4 sector is positive. Therefore, there is no sign problem in the extended configuration space with worms.  , Figure 4. Extrapolation of the density correlations at the largest distance C(R max ) and the CDW structure factor M 2 using a/L+b/L 2 . The error bar of the extrapolation to the thermodynamic limit (1/L = 0) is evaluated using a jackknife analysis. Figure 3 shows the density-density correlations (12), which develop a staggered pattern as the interaction strength V increases. The density correlation at the farthest distance C(R max ) and the CDW structure factor M 2 approach the square of the CDW order parameter as the system size increases. Figure 4 shows the extrapolation of C(R max ) and M 2 to the thermodynamic limit using a/L + b/L 2 . The extrapolation suggests that the quantum critical point lies between V = 1.3 and V = 1.4.  To better estimate the the critical point we perform a finite size scaling (FSS) analysis based on the scaling ansatz

Quantum Monte Carlo Results
where F and G are universal functions and ν, η are the critical exponents. This scaling ansatz holds close to the critical point. The Binder ratios of different system sizes cross at the transition point. This provides a rough estimate of the transition point V c 1.36, as shown in Fig. 5. We next collapse the data of M 2 and M 4 to determine the transition point V c and the critical exponents ν, η. The results are summarized in Table 1 where we also list the estimates of the order parameter critical exponentβ = ν 2 (1 + η), which we will compare with iPEPS results in Sec.3.2. We have performed the data collapse using all available system sizes (L = 6, 9, 12, 15) and excluding the smallest system size (L = 6). They both give satisfactory data collapse where the χ 2 per degree of freedom (d.o.f) is close to one. To visually examine the quality of the data collapse, Fig. 6(a-b) shows the scaled M 2 and M 4 using η = 0.3 where all the curves intersect around V = 1.36. Further scale the horizontal axis using V c = 1.356 and ν = 0.8 collapse all the data onto a single curve, Fig.6(c-d).
Our estimation of the correlation length exponent ν agrees with earlier -expansion result ν = 0.797 [11] and functional renormalization group results ν = 0.738 ∼ The system sizes up to L = 15 do not allow us to reliably pin down the critical point based on 1/L extrapolation [19,20]. Table 1. The critical point and critical exponents determined using data collapses of M 2 and M 4 for all systems sizes (L = 6, 9, 12, 15) and excluding the smallest one (L = 9, 12, 15). The critical exponentβ (to avoid confusion with the inverse temperature β) is calculated usingβ = ν 2 (z + η). The estimated uncertainty [36] of the last digit is shown in the bracket. The χ 2 /d.o.f listed in the last row shows the quality of the data collapse. L = 6, 9, 12, 15 L = 9, 12, 15 0.927 [12,13]. However, our estimated anomalous dimension η ≈ 0.3 is smaller than the previous estimates η = 0.502 ∼ 0.635 [11,12,13]. We have checked that these values of η are not consistent with our QMC data. These field theory calculations [11,12,13] treated Dirac fermions with the same chiralities but our lattice model contains two Dirac fermions with opposite chirality. This difference might explain the discrepancies with our QMC data. On the other hand, since we observed subleading corrections in the Binder ratio crossing Fig. 5, further research using larger systems may be needed to determined the critical behavior more accurately.

Comparison with iPEPS Results
As an independent check of the results we have studied the model (1) with infinite projected entangled-pair states (iPEPS) -a variational tensor-network ansatz for twodimensional ground-state wave functions in the thermodynamic limit [45,46,47,48]. This ansatz is a natural extension of matrix product states (the underlying ansatz of the density-matrix renormalization group method) to two dimensions, and has been previously applied to the same model for attractive interactions [49,50]. Twodimensional tensor networks have first been introduced for spin models and later extended to fermionic systems [51,52,53,54,55,56,48,57]. The iPEPS ansatz consists of a cell of tensors with one tensor per lattice site, which is periodically repeated on the lattice. Each tensor has a physical index of dimension d which carries the local Hilbert space of a lattice site and z auxiliary indices which connect to the z nearest-neighboring tensors. The number of variational parameters (i.e. the accuracy of the ansatz) can be controlled by the bond dimension D of the auxiliary indices, where each tensor contains dD z variational parameters with d = 2 and z = 3 for the present model. For a general introduction to fermionic iPEPS we refer to Ref. [48].
We simulate the honeycomb model (1) by mapping it onto a brick-wall square lattice, as done in Ref. [58]. The variational parameters of the iPEPS ansatz are    optimized by performing an imaginary time evolution using a second order Trotter-Suzuki decomposition and the full-update scheme for the truncation of a bond index (see Ref. [48] for details). To evaluate the iPEPS wave function (e.g. for the computation of expectation values) we use a variant of the corner-transfer-matrix method [59,60] described in Refs. [61,62]. The U(1) symmetry of the present model is exploited [63,64] to increase the efficiency of the simulations. Since iPEPS represent a wave function in the thermodynamic limit, symmetries of the Hamiltonian can be spontaneously broken, and the order can be measured by a local order parameter. In Fig. 7 the iPEPS results for the CDW order parameter OP CDW = | n A −n B | as a function of V is shown, wheren A andn B correspond to the particle density on sublattices A and B, respectively. Since iPEPS is an ansatz in the thermodynamic limit, there are no finite size effects. However, the finite bond dimension D has a similar effect on the order parameter as a finite size system, i.e. there is no sharp transition but the order parameter is overestimated around the critical coupling V c . To obtain an estimate of the order parameter in the infinite D limit we extrapolate the data linearly in 1/D, shown by the black diamonds in Fig. 7. The error bar indicates the range of extrapolated values by taking into account different sets of data points. Although the analytical dependence of the order parameter on D is not known, empirically, one can get a reasonable estimate by such type of extrapolations (see e.g. Ref [65]). Based on these extrapolations of the iPEPS data up to D = 9 we obtain a value of the critical coupling of V c = 1.36 (3), in agreement with the CTQMC result.
The green crosses in Fig. 7 show the CTQMC data for the order parameter in the thermodynamic limit, computed as OP CDW = lim L→∞ 2 M 2 (L), which is in agreement with the iPEPS data. Similar results are obtained by estimating the order parameter from C(R max ), i.e. OP CDW = lim L→∞ 2 C(R max ). The extrapolation of QMC data close to the critical point is more difficult because the intersections at 1/L = 0 may become negative.
We also tried to extract the critical exponentβ by fitting the extrapolated iPEPS data to k(V − V c )β in the range [V c , 1.6]. However, due to the error bars and sensitivity of the exponent on the fitting range we can only give a crude estimate ofβ = 0.7 (15), which is somewhat larger than the CTQMC resultβ = 0.52 (6), but both are smaller than the mean-field resultβ M F = 1 [18] and consistent with the concave shape of the order parameter versus V curve.

Conclusion and Outlook
We presented a sign problem free CTQMC study of the Dirac semi-metal to charge-density-wave transition on the honeycomb lattice and compare it with theory and iPEPS results. Our main results about the transition point and the critical exponents are summarized in the Table 1. The present study uses the static densitydensity correlations as the diagnosis tool for the quantum critical point, it is however interesting to further study the transport and entanglement properties across the phase transition. Future studies may map out the finite temperature phase diagram and especially the crossover [38] from the Gross-Neveu to the 2D Ising universality class. The CDW transition of the spinful Dirac fermions [39,40] can also be studied using a similar method. Generalization of the model to include hopping and interactions beyond the nearest neighbors may allow us to address the intriguing question about the emergence [41] and stability of the topological insulating states [42,43] in the presence of interactions.

Appendix A. Monte Carlo Updates
A Monte Carlo update consists of proposing a move from a configuration C to a new configuration C with a priori probability A(C → C ). The acceptance probability R(C → C ) satisfies the detailed balance condition The Metropolis-Hasting solution of the detailed balance equation Eq.(A.1) is There are three classes of configurations shown in Fig. A1. We devised several updates to sample the configuration space. Most updates are in complementary pairs. Within each pair one can still fine tune the propose probability to enhance the acceptance rate.

Appendix A.1. Vertex add/remove
We add n vertices to a configuration with k vertices. The acceptance ratio is where N b = 3L 2 is the number of bonds of the honeycomb lattice. The move Eq.(A.3) is balanced by removing n vertices from a k-vertices configuration with the acceptance probability Replace G k by G k;ij(kl)τ one gets the formulas for adding/removing vertices in the worm space.
A cheaper way to go between the partition function and the W 2 space is to randomly select a vertex and interpret it as a worm. We call this process an open update and the reverse process a close update. These two updates change the perturbation order by one. However they are cheaper than creating/destroying worms Eq.(A.5-A.6) because no matrix operation is involved. The acceptance ratios are The factor 2 accounts for the fact that {i, j; τ } and {j, i; τ } are counted as two distinct worm configurations.
In the W 2 sector we insert another worm at {k, l; τ } choosing a random site k and a nearly site l (out of m sites). The time τ is the same as the imaginary time of the existing worm {i, j; τ }. Acceptance ratios are We create a worm at {i, j, k, l; τ } in the partition function sector. To improve the acceptation ratio, we choose the sites j, k, l in the neighborhood of a randomly chosen site i. The ratios are

. Worm shift
We shift the worm to a new space-time point. To enhance the acceptance probability, we randomly choose one site in the worm and shift it to one of its neighbors. The imaginary time τ is updated to τ by randomly by adding a random number in the range of [−0.05β, 0.05β). The matrix is updated to G k;ij τ and the acceptance probability is This process is self-balanced. The acceptance rate of worm shift in the W 4 space has a similar expression. To further confirm the critical exponent found in the main text, we simulated the model Eq.(1) on a square lattice with π-flux inserted in each plaquette, Fig.B1(a). The lattice also features two Dirac points in the Brillouin zone, Fig.B1(b). The Dirac semimetal to CDW transition should belong to the same universality class in the honeycomb lattice. In the simulation we use the Landau gauge for the flux and choose system sizes L to be divisible by 4. The inverse temperature scales linearly with length β = L.  Figure B2. The Binder ratio Eq.(15) versus V for different system sizes of the π-flux lattice. Lines are linear interpolations of the data. Fig. B2 shows the Binder ratios, from which we infer the transition point V c ≈ 1.3. A data collapse analysis of M 2 ( Figure B3) gives V c = 1.304(2), ν = 0.80(6) and η = 0.318 (8). These results indicate that the critical exponents we found for the honeycomb lattice (Table 1)