Pinning the order: the nature of quantum criticality in the Hubbard model on honeycomb lattice

In numerical simulations, spontaneously broken symmetry is often detected by computing two-point correlation functions of the appropriate local order parameter. This approach, however, computes the square of the local order parameter, and so when it is {\it small}, very large system sizes at high precisions are required to obtain reliable results. Alternatively, one can pin the order by introducing a local symmetry breaking field, and then measure the induced local order parameter infinitely far from the pinning center. The method is tested here at length for the Hubbard model on honeycomb lattice, within the realm of the projective auxiliary field quantum Monte Carlo algorithm. With our enhanced resolution we find a direct and continuous quantum phase transition between the semi-metallic and the insulating antiferromagnetic states with increase of the interaction. The single particle gap in units of the Hubbard $U$ tracks the staggered magnetization. An excellent data collapse is obtained by finite size scaling, with the values of the critical exponents in accord with the Gross-Neveu universality class of the transition.


I. INTRODUCTION
Detecting spontaneous symmetry broken phases in numerical simulations often relies on the measure of correlation function. For instance, the magnetically ordered phase is characterized by long ranged spin-spin correlations, whereas the superconducting state exhibits long ranged pair correlations in the appropriate symmetry channel. A fundamental caveat with such an approach is that one measures the square of the order parameter. If the later is small, the quantity one attempts to obtain by extrapolating numerical data to the thermodynamic limit is quadratically smaller. As a consequence, very large system sizes at high precision are required in addition to an appropriate finite size extrapolation formula.
The aim of this article is two-fold. We will document a simple and very efficient alternative method to detect magnetically ordered phases in SU (2) invariant Hubbard type models in the realm of projective quantum Monte Carlo methods. With the enhanced resolution, we will revisit the semi-metal to insulator transition on graphene's honeycomb lattice, which has recently been under considerable debate [1,2].
Honeycomb lattice is a bipartite, non-frustrated lattice, which at half filling and small Hubbard repulsion U hosts the semi-metallic state of electrons, as in graphene. When the repulsion is increased one expects eventually a phase transition into an insulating state with antiferromagnetic order [3][4][5], which due to gapless Dirac fermionic excitations being present on the semimetallic side, should belong to a particular, Gross-Neveu universality class [5,6].
Starting from the strong coupling limit and noting that the insulator to metal transition occurs at values of the Hubbard interaction lesser than the bandwidth, allows for the proliferation of higher order ring exchange terms in an effective spin model aimed at describing the magnetic insulating state in the vicinity of the transition [7]. This point of view opens the possibility that the melting of the magnetic order is unrelated to the the metalinsulator transition. Recent quantum Monte Carlo calculations [1] suggested that there is an intermediate spin liquid phase with a single-particle gap but no magnetic ordering, separating the semi-metal and magnetic insulator. Similar results have been put forward for the related π-flux model on the square lattice [8]. The results of Ref. [1] have been challenged by recent studies. Entropy calculations do not favor ground state degeneracy as expected for the Z 2 spin liquid [9]. Moreover, Ref. [2] shows that extrapolating from significantly larger system sizes would suggest almost complete disappearance of the spinliquid from the phase diagram. The latter conclusion is reinforced here, where we find excellent data collapse and identical finite-size scaling of both single-particle gap and staggered magnetization, with the distinct values of critical exponents, in accord with the Gross-Neveu universality class [5,6].
Form the technical point of view, our approach is very similar in spirit to an approach considered in Ref. 10. By introducing a local magnetic field at say the origin, we explicitly break the SU(2) spin symmetry. In the presence of long range order and in the thermodynamic limit, any field will pin the order along the direction of the external field. Thereby, order can be detected by computing directly the magnetization infinitely far from the pinning field. The upside of such an approach is that one measures directly the order parameter rather than its square. This amounts to evaluating a single particle quantity which is often much more stable than correlation functions. The downsides are three-fold. One explicitly breaks SU(2) spin symmetry such that spin sectors mix and it becomes computationally more expensive to reach the ground state. Since the computational cost scales linearly with the projection parameter, this problem is tractable. The second difficulty lies in the ordering of limits. To obtain results which are independent on the magnitude of the pinning field, it is important to first take the thermodynamic limit and then the limit of infinite distance from the pinning field. In a practical implementation, this ordering of limits has as a consequence some leftover dependence of the magnetization on the magnitude of the pinning field. This is particularly visible when the pinning field is small. The final drawback is that it is not always possible to introduce a pinning field without generating a negative sign problem. For instance, in the Kane-Mele Hubbard model [11][12][13], the spin order lies in the x-y plane. Adding a magnetic field along this quantization axis introduces a sign problem. On the other hand, the method is applicable to SU(N) symmetric Hubbard-Heisenberg models [14].
The organization and main results of the article are the following. We will focus on the Hubbard model on honeycomb lattice at the filling one half, for which the presence of an intermediate spin-liquid phase has been controversial [1,2]. After introducing and testing the approach in the next section, we will provide a phase diagram of the Hubbard model in Sec. III. The data points to the fact that the staggered moment follows rather precisely the single particle gap, when the latter is measured in the natural units of the Hubbard U , suggesting a direct quantum phase transition between the semi-metallic and the insulating antiferromagnetic phases. Furthermore, an excellent finite size scaling of the data for both the staggered magnetization and the single particle gap is found by assuming the values of the critical exponents β = 0.79 and ν = 0.88. These are the values found in the first order expansion for the Gross-Neveu-Yukawa field theory of this quantum phase transition [6], around its upper critical (spatial) dimension of three. Altogether, the data strongly supports the existence of a single quantum critical point separating the semi-metallic and the insulating antiferromagnetic phases of the Hubbard model, with the quantum criticality belonging to the Gross-Neveu universality class [5].

II. MODEL AND METHOD
As in Ref. [1], we will consider the half-filled Hubbard model on the Honeycomb lattice (1) The hopping is restricted to nearest neighbors so that the bipartite nature of the lattice allows us to avoid the negative sign problem.
Generically, to detect anti-ferromagnetic ordering we compute spin-spin correlations: Here N = 2L 2 corresponds to the number of orbitals, and L is the linear length of the lattice. A finite value of m signalizes long range order and is equivalent to spontaneous symmetry breaking. In particular, including a magnetic field term with appropriate Fourier component, gives The ordering of limits is crucial. One first has to take the thermodynamic limit to allow for the collapse of Anderson's tower of states, and then the limit of vanishing magnetic field h. Such an approach was for instant used in Ref. 15. It is more convenient to consider a local field since, as we will see below, this lifts the burden of taking the limit h → 0 numerically. The local pinning field is given by the term in the Hamiltonian. Using the representation δ i,0 = 1 L 2 q e iq·i of the Kronecker symbol shows that each Fourier component comes with an amplitude h 0 /L 2 so that taking the thermodynamic limit is equivalent to taking the amplitude of the relevant Fourier component to zero. With the local field construction, the appropriate ordering of limits for an L × L lattice reads: That is, one first has to take the thermodynamic limit -again to guarantee the collapse of the tower of states in the presence of long range order -and only then can one take the distance from the pinning center to infinity [16]. In other words, the distance from the pinning center sets an energy scale which has to be larger than the finite size spin gap. As an efficient estimator for the evaluation of the ordered moment we thus propose: We have tested the above approach for the Hubbard model on the honeycomb lattice. Ground state calculations were carried out with the projective auxiliary field Here θ is a projection parameter, and the trial wave function is required to be non-orthogonal to the ground state.
For H = H tU + H loc the inclusion of the magnetic field does not generate a negative sign problem. We have chosen the trial wave function to be the ground state of the non-interacting Hamiltonian H T = H t + H loc in the S z = 0 sector. The implementation of the algorithm follows closely Refs. 1 and 17. The major difference is the use of a symmetric Trotter breakup which ensures hermiticity of the imaginary time propagator for any value of the time discretization ∆τ . It also leads to smaller systematic errors. Figure 1 (a) plots the local moment at U/t = 5 using different methods. In this case, magnetic ordering is robust such that various approaches can be compared. The data set at h 0 = 0 corresponds to the correlation functions of Eq. 2. For this set of runs we used a spinsinglet trial wave function and the projection parameter θt = 40 suffices to guarantee convergence to the ground state. This quick convergence stems from the fact that the trial wave function is orthogonal to the low lying spin excitations [18]. The runs at finite values of the pinning field correspond to the quantity of Eq. 7. In the presence of a finite pinning field, SU(2) spin symmetry is broken and the trial wave function overlaps with all spin sectors. Consequently, a large value of the projection parameter θt = 320 is required to guarantee convergence to the ground state within the quoted accuracy. Note that the CPU time scales only linearly with the projection parameter, so that such large projection parameters are still numerically tractable. It is also worth pointing out that the observable of Eq. 7 corresponds to a single particle quantity and shows very little fluctuations. As evident in Fig. 1(a), finite size effects are strongly dependent on the specific choice of the pinning field. Nevertheless, convergence to values consistent with the generic approach based on Eq. 2 is obtained for relatively large values of the pinning field. If the pinning field is chosen too small, larger lattices are apparently required to insure that the finite size spin gap, set by v/L with v the spin wave velocity, is smaller than the energy scale set by the pinning field. This expectation is confirmed by the data, which shows a systematic upturn as a function of the system size for smaller values of the pinning field. Such a non-monotonic finite size behavior complicates a finite size scaling analysis. For this reason, we propose to use a relatively large value of the pinning field [19]. at U/t=5 (a) and U/t = 4 (b). At U/t = 4, the local moment is smaller, and hard to detect. The h 0 = 0 data set of Fig. 1(b) compares our results for the correlation function of Eq. 2 to those of Ref. 2. As apparent, the agreement up to our largest lattice size, L = 18, is remarkable. Without the largest lattice sizes, L = 24 and L = 36, extrapolation to the thermodynamic limit is hard due to the downward turn present in the finite size results. The data sets stemming from the pinning field approach provide an alternative perspective, and, on the whole, confirm the result of Ref. 2. For the considered field range there is considerable scatter in the finite size results, but nevertheless the extrapolation to the thermodynamic limit seems to be field independent, as expected from the above considerations.
We conclude this section by mentioning that we have tested the approach for the non-interacting case and the method successfully demonstrates the absence of long range magnetic order. Hence both for critical states as well as for magnetically ordered phases the pinning field approach does provide an efficient tool. Further testing of Here Z corresponds to the single particle residue and ∆sp to the single particle gap. (b) Size dependence and extrapolation of the single particle gap.
this approach for Heisenberg bilayers is presently under progress [20] .

III. PHASE DIAGRAM OF THE HUBBARD MODEL ON HONEYCOMB LATTICE.
We have used the above approach to revisit the magnetic phase diagram of the Hubbard model on the honeycomb lattice. At weak couplings the model is known to have a stable semi-metallic state. In the strong coupling limit, and due to the absence of frustration, an antiferromagnetic Mott insulator is present. The nature of the transition between these two states has been studied in the past, [3][4][5] and is presently controversial [1,2]. A. Single particle gap.
To pin down the coupling strength beyond which the single particle gap opens, we have repeated calculations for the time displaced single particle imaginary time Green function at the nodal point: G(K, τ ) = σ c † K,σ (τ )c K,σ (τ = 0) . As evident in Fig. 2, and with the symmetric Trotter decomposition, the Trotter systematic error is negligible within our accuracy. Fitting the data to an exponential form allows us to extract the single particle gap, which we plot as a function of system size in Fig. 2. Assuming a polynomial form for the extrapolation to the thermodynamic limit, we find a small but finite single particle gap for U/t ≥ 3.7. This finding is particularly interesting when compared to the results of Sorella et al [2] which at U/t = 3.8 point to the absence of long range magnetic order. This could be taken as an indication for a possible intermediate phase.
However, the analysis that we present below suggests a different interpretation.

B. Magnetization from pinning fields.
We have used the pinning field approach to compute the staggered moment as a function of U/t. The results at h 0 = 5t are reported in Fig. 3. The extrapolation to the thermodynamic limit is carried out using a polynomial scaling up to second order in 1/L. Fig. 4 plots the so obtained staggered moment for two choices of the pinning field as well as the single particle gap. Several comments are in order.
• Within our accuracy, and maybe most importantly, with the polynomial fit used in extrapolating the data to the thermodynamic limit, it appears that the single particle gap opens right when magnetic ordering sets in. The only mismatch is at U/t = 3.7 where we do not detect magnetic ordering but we do detect a small single particle gap.
• The QMC data in Fig. 4 shows that over a wide parameter range, the single particle gap measured in units of the Hubbard U, tracks the staggered magnetization. We take this as a strong indication, that the magnetization provides the only relevant scale in the problem, determining directly the single particle gap. We will see below, that this conclusion, based here on a simple, polynomial extrapolation of the finite size data, is also obtained, if a more refined data analysis is performed.
• The data in Fig. 4 exhibits an unusual inflection point at approximately U/t = 4.1. Such an inflection point is clearly absent at the mean-field level (see inset of Fig. 4). We will discuss the implications of this inflection point in the next section. Let us finally note, that in previous calculations [1] we were unable to resolve staggered moments lesser than m 0.03. We thereby missed this inflection point in the polynomially extrapolated magnetization curve and concluded the presence of an intermediate phase [21].

C. Finite size scaling
As mentioned above, one of the particularities of the data presented in Fig. 4 is the occurrence of an inflection point at U/t = 4.1. It is a natural question to ask if this rather peculiar feature may be an artifact of using a simple polynomial fitting procedure, which one would indeed expect to fail close to criticality. This could result in an overestimation of the magnetization in the vicinity of the critical point between the semi-metallic and the insulating phase of the Hubbard model. As we explain next, arguments in favor of this conjecture are provided by the large-N treatment of the Gross-Neveu model [5], and the -expansion around three spatial dimensions in the equivalent Gross-Neveu-Yukawa field theory, formulated in Ref. 6. Given the order parameter exponent, β, as well as the correlation length exponent, ν, the staggered magnetization scales as Using the standard scaling laws [22], the exponent β/ν may conveniently be expressed in terms of the anomalous dimension for the order parameter η, as where d + z is the effective dimensionality of the system. If we assume that the Lorentz invariance is emergent at the critical point, as it indeed is close to the upper critical dimension d up = 3 of the Gross-Neveu-Yukawa theory [6], and maybe even more generally [23], the dynamical critical exponent is z = 1. If then the anomalous dimension of the order parameter is such that η < 3−d, we find that the combination of the exponents β/ν < 1, and our polynomial fitting procedure in the previous section could very well overestimate the value of the staggered magnetization. In fact both the large-N approach [5] and the expansion around the upper critical dimension [6] suggest that this is indeed the case. Within the first order of the expansion in the parameter = 3 − d, for example, η = 4 /5, so that In two dimensions then, β/ν 0.9.
To look for the signs of the Gross-Neveu criticality in the Hubbard model we have carried out a finite size scaling analysis based on the usual scaling form Figure 5(a) plots mL β/ν versus U for the magnetization data at the fixed field h 0 = 5. (We will omit at this point the second scaling variable, h 0 L y−d , since the scaling dimension y − d = ( − η)/2 1. This second argument of the scaling function, present in principle, is therefore effectively constant at a fixed h 0 , and its inclusion does not visibly affect the quality of scaling. For further discussion of this point, see the Appendix.) As a guide, we have used the first-order -expansion value of β/ν = 0.9. Five curves then all cross at a single point U c /t 3.8, thereby providing a first non-trivial indication of the critical point. This value of U c is slightly larger than that obtained with the polynomial fit. The -expansion value of the correlation length exponent reads [6], With this value of ν, again at = 1, we obtain an excellent data collapse, as shown in Fig. 5(b).
In accord with the Gross-Neveu-Yukawa theory, the numerical data of Fig. 4 support the interpretation that the magnetization is the only scale in the problem. To further check this interpretation, we have scaled the single single particle gap to the form: again using the -expansion exponents in two dimensions. Fig. 6(a) shows that the crossing point of the ∆sp U L β/ν curves again occur at U c /t 3.8 and Fig. 6(b) shows the collapse. It is quite remarkable, that within our precision, the two scaling functions are equal,F = F .
Hence, the scaling analysis of our QMC within the Gross-Neveu scenario is consistent with a single continuous quantum phase transition between the semimetal and the antiferromagnetic insulator, and suggests that the first-order expansion around the upper critical dimension in the Gross-Neveu-Yukawa theory may already yield rather accurate values of the critical exponents.

IV. DISCUSSION AND CONCLUSIONS
We have introduced an alternative method to compute the staggered magnetization. By introducing a local magnetic field we pin the quantization axis of the ordered state. The staggered magnetic moment then corresponds to the local magnetization infinitely far away from the pinning center. The approach has the major advantage that we compute directly the magnetization as opposed to its square when measuring correlation functions. We can therefore expect improved resolution when the local moment is small. One advantage of the approach is an internal cross-check which requires the staggered moment to be independent on the numerical value of the pinning field. We have been able to reach this internal cross-check only in the case of large pinning fields. This is consistent with the approach proposed by [10] where the pinning field is set to infinity on the boundary of the lattice. If the pinning field is too small, the approach suffers from large and non-monotonic size effects, since the energy scale set by the local magnetic field is unable to overcome the finite size spin gap. Under these circumstances the extrapolated value of the magnetization has the tendency of underestimating the order. In this article we have only considered a very specific form of the pinning field. The number of different choices of pinning fields provides a playground for optimization of the approach, and for minimization of the size effects.
The application to the Hubbard model on the honeycomb lattice sheds new light on the phase diagram of this well known problem. The enhanced precision in comparison to Refs. 1 and 2 reveals that the staggered moment has the same functional form as the single particle gap (as measured in units of U ). Remarkably, an excellent data collapse onto a single universal curve is found in the finite size scaling of both quantities, with the values of the critical exponents characteristic of the Gross-Neveu criticality between the semimetallic and the magnetic insulating phases.

ACKNOWLEDGMENTS
FFA is very indebted to T. Lang, Z.Y. Meng, A. Muramatsu and S. Wessel for many valuable comments. We equally acknowledge discussions with A. Chernyshev, L. Fritz, M. Hohenadler and S. Sorella. This work was initiated at the KITP during the workshop on "Frustrated Magnetism and Quantum Spin Liquids From Theory and Models to Experiments" coordinated by K. Kanoda, P. Lee, A. Vishwanath and S. White and finalized at the MPI-PKS (Dresden). Funding from the DFG under the grant number AS 120/9-1 and AS120/10-1 (Forschergruppe FOR 1807) is acknowledged. IFH has been suported by the NSERC of Canada. We thank the Jülich Supercomputing Centre and the Leibniz-Rechenzentrum in Münich for generous allocation of CPU time.

V. APPENDIX
Strictly speaking, the finite-size scaling form for the magnetization in presence of the external field is where G(X, Y ) is the scaling function of two variables, and ξ is the (diverging) correlation length. Here we have used the fact that the relevant Fourier component of the local pinning field scales as h0 L d . The dimension y of the uniform external field is [22] y = d + z + 2 − η 2 .
(16) Away from criticality ξ is bounded, and in the thermodynamic limit, the magnetization scales to its field independent value. The data in Fig. 1 confirms this and shows that the extrapolated value of the magnetization is h 0 independent provided that, for the considered lattice sizes, h 0 is chosen to be large enough. It is also interesting to point out that for values of h 0 in the range 1t < h 0 < 5t the finite size value of the magnetization is next to independent on the value of the pinning field.
At criticality we can replace ξ by L to obtain: With the Lorentz symmetry at the critical point the value of the dynamical critical exponent is pinned to z = 1, and the scaling dimension of the local field h 0 is In the -expansion then, and rather small, presumably even for = 1 (d = 2). For a fixed value of h 0 , therefore, the second argument of the scaling function G is almost constant. To be more precise, we can use the asymptotic form G(1, Y ) ∝ Y 1/δ with 1 δ = β νy [22] such that in the large-L limit, m ∝ h Hence on the whole, the scaling function G depends rather weakly on the second argument and for practical purposes it suffices to neglect it.