Trion states and quantum criticality of attractive SU(3) Dirac fermions

We perform the projector quantum Monte Carlo (QMC) simulation to study the trion formation and quantum phase transition in the half-filled attractive SU(3) Hubbard model on a honeycomb lattice. With increasing attractive Hubbard interaction, our simulations demonstrate a continuous quantum phase transition from the semimetal to charge density wave (CDW) at the critical coupling $U_c/t=-1.52(2)$. The critical exponents $\nu=0.82(3)$ and $\eta=0.58(4)$ determined by the QMC simulation remarkably disagree with those of the $N=3$ chiral Ising universality class suggested by the effective Gross-Neveu-Yukawa (GNY) theory, but coincide with the $N=1$ chiral Ising universality class. In the CDW phase, we show that on-site and off-site trions coexist and the off-site trion forms a local bond state. Our work not only illustrates the formation of off-site trions in two-dimensional Hubbard model, but also raises doubts about the extent of applicability of GNY model on the attractive SU(3) Dirac fermions.

The SU(3) 6 Li Fermi gas with tunable interactions is exceptionally appropriate for studying the three-body bound states [27][28][29]. In few-body spin-1/2 fermion systems, forming three-body bound states by on-site attractions is rather challenging due to the Pauli exclusion [30][31][32], although threebody bound states can be realized with finite-radius interactions [33]. However, in attractive few-body three-color fermions, due to an additional internal degree of freedom, three-body bound states can be formed by on-site attractions. Specifically, in the attractive SU(3) Hubbard model with only three fermions on the honeycomb lattice, the three-body bound * yu.wang@whu.edu.cn states have two configurations under different conditions [34]: when on-site triple occupancy is energetically favorable, the three-body bound state is an on-site trion, which consists of three fermions at one site; when on-site triple occupancy takes great energy penalty, the three-body bound state is an offsite trion, which consists of two fermions at one site and one fermion at the nearest-neighbor site. In many-body systems the three-body bound states can develop long-range order. The variational [35,36], self-energy functional [37,38] and dynamical mean-field theory (DMFT) [39,40] studies of the SU(3) Hubbard models on two-dimensional lattices have predicted a phase transition between the color superfluid and the on-site trion phase, which is reminiscent of the transition between quark superfluid and baryonic phase [20][21][22]. In one-dimensional lattices, the density matrix renormalization group (DMRG) studies have observed off-site trion phase either when on-site triple occupancy is prohibited in an attractive SU(3) Hubbard model [41] or when the attractions are color-dependent and thus anisotropic in a three-color Hubbard model with SU(3) symmetry breaking [42].
Quantum Monte Carlo (QMC) simulations of the SU(3) Hubbard model have long been absent due to the notorious sign problem. Owning to recent advances in the QMC algorithm [43][44][45][46], the QMC simulation of the attractive SU(3) Hubbard model can be proved to be sign-problem-free in bipartite lattices at half filling. In this work, we propose, for the first time, to conduct a projector determinant QMC simulations of the half-filled attractive SU(3) Hubbard model on the honeycomb lattice. The purpose of this work is twofold: firstly, QMC simulations of the attractive SU(3) Hubbard model can demonstrate the formation of trionic states, and in particular we shall show that the off-site and on-site trions coexist in this two-dimensional Hubbard model while previously off-site trions have only been found in one-dimensional Hubbard model; and secondly, the QMC simulations of the quantum phase transition of SU(3) Dirac fermions can be compared with the Gross-Neveu-Yukawa (GNY) models which is thought to provide a general description for the criticality of Dirac fermions in two spatial dimensions [47]. Given the symmetry breaking patterns, the GNY descriptions are irrelevant to the details of the microscopic Hamiltonian, but solely dependent on the number of fermion colors, N . So far various QMC calculations are consistent with the GNY models. In the spinless Dirac fermions, the transition between semimetal and chargedensity-wave (CDW) phases belongs to the N = 1 chiral Ising universality class [48,49]. In SU(2) Dirac fermions, the semimetal-CDW, semimetal-antiferromagnet and semimetalsuperconductor transitions fall into the N = 2 chiral Ising [50], chiral Heisenberg [51,52] and chiral XY [53] universality classes, respectively. In SU(4) Dirac fermions, the transition between semimetal and valence bond solid (VBS) belongs to the N = 4 chiral XY universality class, due to an emergent U(1) symmetry at the critical point [54][55][56]. Moreover, in SU(N ) Dirac fermions with singlet-bond interactions, the semimetal-VBS transitions fall into the chiral XY universality classes of N fermion colors [57]. However, we shall show that our QMC simulation of the quantum phase transition in the SU(3) Dirac fermions with attractive Hubbard interactions surprisingly conflicts with the N = 3 GNY model.

II. MODEL
At half filling, the attractive SU(3) Hubbard model is defined by the lattice Hamiltonian: (1) where ij denotes the nearest-neighbor sites on a honeycomb lattice and α, β are the color indices running from 1 to 3. U < 0 describes the on-site attractive Hubbard interaction, and n iα = c † iα c iα is the particle number operator for color α at site i. The nearest-neighbor hopping amplitude t = 1 is set as the energy unit in our simulations.
In the atomic limit t/U = 0, on-site trions are randomly formed by three fermions on a lattice point in the half-filled model. When the fermion hopping is turned on, the charge fluctuations induce the CDW order [15]. As shown in Fig. 1, the energy penalty for adding an extra fermion (hole) to the half-filled system is −U in the atomic limit. At strong coupling the extra fermion (hole) can hop on a triangular lattice via second-order perturbation, which expands the energy level into an energy band. Then the energy gap is ∆ at ≈ −U − W where W = − 3t 2 U . At the Mott transition point, the energy gap vanishes, and then the critical coupling strength is estimated as U at c /t ≈ − √ 3 (see Appendix A). The determinant formalism of the projector QMC method [58] will be employed to simulate the attractive SU(3) Hubbard model, which is sign-problem-free at half filling when decomposing the Hubbard interaction into the on-site colorflip channel. Implementation details of the algorithm are elaborated in Appendix B.

III. TRION FORMATION
For studying the trion formation, we consider the correlation function [41,59], which measures the correlation between the color-1 fermion at site i and the color-2 and color-3 fermions at site j. Then the probability of triple occupancy can be defined as where L is the lattice size of a honeycomb lattice. At strong coupling, P (3) corresponds to the occupancy probability of on-site trion.
In the noninteracting limit, the density-density correlation can be decoupled directly, and thus T (i, i) follows the binomial distribution, lim U/t→0 P (3) = 1 8 . In the atomic limit, since only on-site trions exist, each site is either fully occupied or empty, and then lim t/U→0 P (3) = 1 2 . In Fig. 2(a), P (3) is plotted as a function of U for various lattice sizes L. For L = 9 and 12, the P (3) versus U curves are almost indistinguishable, and therefore the lattice size L 9 is sufficiently large to estimate the L → ∞ limit of P (3). As expected, P (3) increases monotonically with U .
As shown in Fig. 2(a), at strong coupling triple occupancy probability P (3) is observably smaller than the large-U limit 1 2 , which infers that there may exist another type of threebody bound state besides the on-site trion in the strong coupling regime. To explore the possible trion states, let us first consider a two-site half-filled attractive SU(3) Hubbard model. In this two-site model, there are obviously two possible trion states: on-site trion |Φ t = |123 and off-site trion |Φ ot = 1 √ 3 P ǫ P P|12, 3 , where the permutation P ∈ {(1), (13), (23)} and ǫ P = (−1) P is the parity. By the first-order perturbation theory, we obtain the ground-state wave function, This suggests that the ground state is the superposition of the on-site trion state |Φ t and the off-site trion state |Φ ot . We may use this result to estimate the probability of triple occupancy on the honeycomb lattice: where the coordination number z = 3. For |U |/t 3, the values of P (3) obtained by QMC simulations can be fitted into the equation P (3) = 1 2 − a t 2 U 2 . Surprisingly, the fitting coefficient a = 1.10(3) quantitatively agrees with Eq. (5).
In Fig. 2(b), the correlation function T (i, j) is plotted as a function of r ij . Here r ij ≡ |r i − r j | is the distance between sites i and j and the distance between nearest-neighbor sites is set as the length unit. For various U , the maximum of each T (i, j) appears at r ij = 0, due to the attractive interaction between fermions, while the minimum of each T (i, j) appears at r ij = 1, optimizing the kinetic energy gain that results in an effective nearest-neighbor repulsive interaction between on-site trions at strong coupling [39,60]. The minimum of T (i, j) decreases with the increase of |U | and vanishes in the strong coupling limit, which is consistent with the off-site trion term in Eq. (4). Additionally, the interaction-induced CDW phase transition can be illustrated in terms of the behavior of T (i, j). For |U |/t 1, T (i, j) converges to a constant when r ij is large, so the correlation function has the equal value for the color-1 fermion occupying the two sublattices, reflecting the lattice inversion symmetry of the semimetal phase. In contrast, for |U |/t 2, T (i, j) defined on the same sublattice is much larger than that defined on different sublattices, which corresponds to the lattice inversion symmetry breaking of the CDW phase.
The off-site trion at site i can be demonstrated via the prob- is the kinetic bond operator; s represents the simulated value during a QMC run;ê a denotes three nearest-neighbor bond orientations [61][62][63][64]. For a reference lattice point, each of the local bond vectors arising over the QMC simulation is tracked and collected. The entire collection of simulated bond vectors is plotted as a histogram which visualizes the probability distribution of the bond vectors in terms of density of data points. In the semimetal phase, the probability distribution P (N x , N y ) are symmetric around the origin since the bond vector 3 a=1 d i,êa sêa is homogenous. In the CDW phase, the off-site trion causes a nonzero d i,êa s . Since there is only one off-site trion on each site, a d i,êa sêa is polarized along oneê a direction. Figure 3(a) shows the symmetric probability distribution P (N x , N y ) at U/t = −1, which indicates the absence of an off-site trion. As shown in Fig. 3 is polarized along the threê e a directions, manifesting the bond state of an off-site trion. The histogram only illustrates the formation of a local off-site trion bond state, and does not imply the emergence of a longrange order with off-site trions due to the low density of off-site trions.

IV. THE CDW PHASE TRANSITION
The CDW ordering can be characterized by the CDW structure factor, which is defined at the Γ point via the densitydensity correlation function C(i, j) = α,β n iα n jβ [65], where ε i = +1 for sublattice A and ε i = −1 for sublattice B. Then the CDW order parameter is given by D = lim L→∞ 1 2L 2 S CDW (L, Γ). In the L → ∞ limit, the CDW order parameter for various U can be obtained by finite-size extrapolation. Figure 4(a) shows that the extrapolated CDW order parameters are substantially greater than zero when |U |/t > 1.5, suggesting that the critical point is at around U/t ≈ −1.5. The CDW order parameters near the critical point obey the power law D ∼ |U − U c | −ζ [53,66], which in turn can be used to fit the extrapolated CDW order parameters and then the critical exponent ζ = −0.67(3) and the critical coupling U c /t = −1.55(1) can be extracted, as shown in Fig. 4(b). Note that the few non-vanishing extrapolated values below critical point may originate from the simple finite-size extrapolation method and should be removed by the curve-fitting approach [53,67]. Compared to the mean-field critical point U MF c /t ≈ −1.11 (see Appendix C), it is reasonable that |U MF c | < |U c | , because quantum fluctuations are neglected in the mean-field approach.
To locate the transition point more accurately, we consider the dimensionless correlation ratio [50], where Γ is the CDW wave vector and Γ + δk represents a neighboring wave vector in the reciprocal lattice. In the semimetal phase, R CDW tends to zero as S CDW (L, Γ) ≈ S CDW (L, Γ + δk).
When the CDW order develops, S CDW (L, Γ) ≫ S CDW (L, Γ + δk), so R CDW approaches unity. For sufficiently large L, R CDW curves intersect at sizeindependent point corresponding to the critical point U c . For finite lattice sizes, the crossing point of R CDW curves defines a finite-size estimate of the critical value U c (L), which takes the form of U c (L) = U c + aL −b , when taking account of scaling corrections [51]. The critical point U c is then extracted in the L → ∞ limit. We use the resampling method [68,69] to extract the crossing points U c (L) between the R CDW (L) and R CDW (L + 3) curves. In Figs. 5(a) and 5(b), R CDW curves show a sizedependent crossing point inbetween −2.0 < U/t < −1.5 due to significant finite-size effects. Nevertheless, as shown in Fig. 5(c), the crossing points are fitted into the curve equation U c (L) = U c + aL −b where the fitting parameters are a = −3.7(11) and b = 2.2(3), and the critical point is found to be U c /t = −1.52 (2).
We shall derive the critical exponents of the semimetal-CDW transition. In the vicinity of the critical point, the CDW order parameter obeys the scaling equation [48,49], where δu = (U − U c )/U c and the exponent ζ = −ν(η + z)/2 (setting z = 1). Figures 6(a) and 6(b) show the scaling collapses of the CDW order parameters and correlation ratios, respectively. The exponent η = 0.57 (2) is extracted from the slope of the log-log plot of D versus L at the critical point U c /t = −1.52, as shown in Fig. 6(c). Then in Fig. 6(d), the exponents η and ν are extracted simultaneously by using the best-fitting analysis adapted from Refs. [70,71]. Typically, we randomly choose the value of U c /t at around −1.52 and randomly choose the initial guesses of the exponents ν and ζ, inside a small range. By the best-fitting procedure, the converged values of the critical exponents are found to be ν = 0.82 (3), ζ = −0.64(3) and η = 0.58 (4).

V. CONCLUSIONS AND OUTLOOK
We have performed the sign-problem-free QMC simulations to investigate the trion formation and quantum phase transition in the half-filled attractive SU(3) Hubbard model on a honeycomb lattice. With increasing attractive Hubbard interaction, the continuous semimetal-CDW transition occurs at the critical point U c /t = −1.52 (2) and the corresponding critical exponents are ν = 0.82(3) and η = 0.58 (4). In the CDW region, the off-site trions emerge due to density fluctuations, and therefore the on-site and off-site trions coexist in the deep CDW phase.
Our QMC simulations illustrate the formation of a local off-site trion bond state in two-dimensional Hubbard model, which extends the understanding of one-dimensional off-site trions suggested by previous DMRG study. It has been proposed that the trionic phase can be probed by shaking the optical lattice [35]. Moreover, the triple occupancy can be determined experimentally by measuring the loss of atoms governed by a three-body process [25]. Our work opens a new avenue for exploring the physical effects of the interplay between on-site and off-site trions in two-dimensional spatial models.
What is particularly interesting about our findings is that the critical exponents determined by QMC simulations remarkably disagree with those of the N = 3 chiral Ising universality class predicted by the effective GNY theory which has been believed to suggest a general description for the criticality of two-dimensional Dirac fermions and has been numerically verified in several SU(N ) models. Unexpectedly, our QMC results are in fair agreement with the N = 1 chiral Ising universality class. We argue that the formation of trions may affect the quantum criticality of attractive SU(3) Dirac fermions. At the critical point, two species of fermions (unbound fermions and trions) get involved in the ongoing development of CDW order. The on-site trion as a whole can be recognized as a spinless fermion, and dominate the long-range CDW ordering at strong coupling. Within the framework of GNY model, the criticality of spinless trions belongs to the N = 1 chiral Ising universality class. Therefore formation of trions completely deprive the attractive SU(3) Hubbard model of the N = 3 chiral Ising universality class. However, the critical point at which unbound fermions are still in the majority is far from the strong coupling regime. The reason that the QMC results are in a good coincidence with the N = 1 chiral Ising universality class is still uncertain. Our results evidently raise doubts about the extent of applicability of GNY model on attractive SU(3) Dirac fermions, which is definitely worth to think and need to solve in future study. . At half-filling, we define the Mott-insulating state by the nonzero energy penalty of adding (removing) a particle into (from) the system [63,75]. In the atomic limit the energy penalty of adding an extra hole into the Mott-insulating state is

ACKNOWLEDGMENTS
where N tot is the total number of fermions. When tuning on the hopping term, the extra hole (particle) can move in the system, which further modifies the energy penalty. Below, we derive an effective model for the description of the extra hole (particle).
Following the steps in Ref. [76], we define the ground state Let P 0 be the projection operator onto the subspace D spanned by |n (0) . The projection operator outside the subspace D is then defined as P 1 = 1 − P 0 . The degenerate Rayleigh-Schrödinger perturbation theory yields the effective Hamiltonian, H eff |n (0) = E n |n (0) , up to the second order, where the hopping operator H 0 P 0 moves one fermion from the local singlet bound state to the NN sites, as shown in Fig. 7. Consider the Pauli exclusion principle, the Hamiltonian can be written as where ij represents the next-nearest-neighbor (NNN) sites and the NNN hopping process is restricted to the fermion color of the extra hole. C (1) = P 0 H 0 P 0 is the first-order perturbation (see Fig. 7(a)), and C (2) is the second-order perturbation where one fermion hops to the NN sites and then hops back (see Fig. 7(b)). The last term represents the effective tightbinding model originated from the desired second-order process, as shown in Fig. 7(c). In the ground state, we assume that the extra hole (particle) occupies the lowest energy level −yt ′ , where t ′ = − t 2 U(N −1) and y > 0 is a constant relying on the lattice structure.
Similarly, the effective Hamiltonian of the half-filled system yields H eff = E Ntot + C (2) where the NNN hopping process is absent up to the second order. Combined with Eq. (A2), we obtain the energy penalty where δC (2) is the difference between the NN second-order correction terms. In fact, δC (2) is marginal and we may also safely discard C (1) because of the CDW ordering. At the Mott transition point, ∆ at = 0 and then an estimate of the transition point is For the attractive SU(N ) Hubbard model on the honeycomb lattice, the NNN hopping processes are on a triangular lattice as shown in Fig. 1, and thus ǫ . Thus, the lowest energy level is ǫ h (Γ) = −6t ′ . Substituting y = 6 into Eq. (A4), we obtain For N = 3, the critical point is estimated as U at c /t = − √ 3.

APPENDIX B: THE AUXILIARY FIELDS COUPLED TO ON-SITE COLOR FLIPS
In this section, we adopt the theorem from Ref. [43] to prove that the QMC simulation of attractive SU(3) Hubbard model           is sign-problem-free at half-filling. We consider the finitetemperature formalism and apply the Suzuki-Trotter decomposition to separate the kinetic and interaction terms in the partition function, where ∆τ = β M is the Trotter decomposition step. For the attractive Hubbard interaction, The coupling matrix between colors α and β reads and its eigenvectors correspond to the complex fermion basis [44],c It can be verified thatc αβ andc βα obey fermionic anticommutation relations: {c † αβ ,c αβ } = 1, {c † αβ ,c βα } = 0. After the diagonalization of coupling matrix Eq. (B3), the Hubbard interaction term (c † α c β − c † β c α ) 2 is refined to −(c † αβc αβ − c † βαc βα ) 2 . In the subspace spanned by the complex fermion basis, we apply the discrete Hubbard-Stratonovich (HS) transformation to decouple the Hubbard interaction term: It is seen that the Hubbard interaction term in the SU(3) Hamiltonian Eq. (1) can be decoupled in the on-site color-flip channels via three auxiliary fields s 12 , s 13 and s 23 : with a systematic error O(∆τ 2 ). We expand Eq. (B7) in the on-site Fock space, and plot the numerical error of the diagonal element as a function of U , as shown in Fig. 8. According to this correctness test, the HS transformation used in our work is accurate in the intermediate coupling regime (U/t < 10) but debatable at very strong couplings (U/t > 10).
Next, we prove that QMC simulations using the HS transformation of Eq. (B7) can avoid the sign problem. In the bipartite lattice, we arrange the color-orbital operators in such order: and Λ = −Λ T . It is easy to verify that the decoupled Hamiltonian satisfies where η = diag(1, . . . , 1, −1, . . . , −1). This condition guarantees the sign-problem-free determinant QMC simulations [43]. However, when the system is away from half filling, D A(B) is not a real anti-symmetric matrix, which may cause the sign problem. Furthermore, we can use the Rodrigues formula to simplify the matrix exponential: e Λ = I 3 + sin θ θ Λ + (1−cos θ) θ 2 Λ 2 with θ = (s 2 12 + s 2 13 + s 2 23 )λ 2 = √ 3λ 2 . The auxiliary fields associated with on-site color-flip channels can be generalized to the attractive SU(2N +1) Hubbard model. This HS transformation, however, enlarges the size of matrix which leads to the time complexity O(β(2L 2 N ′ ) 3 ) with N ′ = 2N +1. The HS transformation of Eq.(B7) can also be applied to the projector QMC (PQMC) method in which the expectation values of observables with the ground state |Ψ G are calculated. A trial ground-state wave function |Ψ T is chosen based on the Hamiltonian Eq.(1) in the noninteracting limit under the anti-periodic boundary condition [43]. Given the inner product Ψ G |Ψ T = 0, the ground state |Ψ G can be reached by applying the projection operator e − βH 2 onto |Ψ T for a sufficiently long projection time β [58,77].
In our simulations, the honeycomb lattice in real space is subject to the periodic boundary condition for L = 3, 6, 9, 12. Under this boundary condition, the lattice vectors in the Brillouin zone can meet the Dirac points. The projection time β = 8 3 L and the Trotter decomposition step ∆τ 0.1 are sufficient for the accurate description of the ground-state properties of the attractive SU(3) Hubbard model.
In the half-filled attractive SU(3) Hubbard model, the results of PQMC are compared with those of exact diagonalization (ED) method on a 2 × 2 square lattice. In Fig. 9, the ED results are shown with the red square, while the PQMC results for β = 10, 20 are presented respectively by the blue circle and blue pentagon. For most of U , PQMC results agree well with ED calculations within the margin of numerical errors. However, when U/t = −0.05, we obtain inconsistent results except the ground state energy. The reason could be that the trial wave function becomes inappropriate due to Fermi surface nesting.

APPENDIX C: MEAN-FIELD ANALYSIS
For the CDW order on a honeycomb lattice, the Hubbard interaction term can be decoupled in terms of the on-site parameters c † iα c iα = ε i ∆ at the mean-field level, where ∆ is the CDW order parameter; ε i is +1 on sublattice A and is −1 on sublattice B. The mean-field CDW order parameter can be defined as The mean-field Hamiltonian of the attractive SU(3) Hubbard model is then written in the reciprocal space, and the bipartite basis c † k ≡ (c † Ak , c † Bk ) is used. The offdiagonal term ǫ(k) = −t a e −ik·êa comes from the noninteracting part of Eq. (1), where a is the sum over three FIG. 9. Correctness test on a 2 × 2 square lattice for (a) ground state energy, (b) pairing order parameter, (c) CDW order parameter, and (d) triple occupancy probability.
. The distance between nearest-neighbor (NN) sites is set as the unit of length. At half filling, the self-consistent equation of ∆ reads where λ k = ± (N − 1) 2 U 2 ∆ 2 + ǫ(k) 2 are the eigenvalues of the matrix h k . The nonzero value of ∆ opens the energy gap 2(N − 1)|U |∆ in the energy spectrum λ k . Thus the ground state is an insulator and Eq. (C1) is often referred to as the gap function. Furthermore, according to Eq. (C3) the number of fermion colors actually rescales U by the factor N − 1 [40]. For each Hubbard U , one can solve Eq. (C3) selfconsistently by the root-finding method [78,79]. As shown in Fig. 10, the CDW phase transitions occur at the critical points U MF c /t ≈ −1.11 for N = 3 and U MF c,N =2 /t = 2U MF c /t ≈ −2.23 for N = 2.
On the other hand, the pairing gap function is defined as can be mapped onto (0, 0, ∆ 12 ) by a global gauge change, leaving only one pair gapped [16]. By solving the standard BCS problem, the energy spectrum has the gapped branches ± ǫ(k) 2 + U 2 ∆ 2 12 with the energy gap |2U ∆ 12 |. Note that Eq. (C3) also holds for ∆ 12 when N = 2. Thus, in the half-filled attractive SU(3) Hubbard model, the pairing gap (N = 2) is certainly not larger than the CDW gap (N = 3), and the ground state is associated with CDW order [15].
Next, we investigate the semimetal-CDW transition by analyzing the Ginzburg-Landau (GL) free-energy density f (∆).
Since σ x ∆σ x = −∆ in which σ x is the Pauli matrix in the basis of sublattices, the CDW order breaks the lattice inversion symmetry. However, f (∆) needs to maintain the lattice inversion symmetry. Hence, the analytic part of f (∆) can be written as f a = r 2 ∆ 2 + r 4 ∆ 4 . (C4) According to the GL theory, f a describes a second-order transition with the critical exponent ζ = − 1 2 . Due to the coupling between ∆ and the gapless Dirac fermions, the total freeenergy density f (∆) potentially contains a nonanalytic part as explained below. When taking account of the spin degeneracy, there are six Dirac cones in the first Brillouin zone. In the CDW phase, the mean-field energy spectrum around each Dirac point is E k = v 2 k 2 + (N − 1) 2 U 2 ∆ 2 , where k represents the deviation from the Dirac point. Denote the nonanalytic part by f non (∆, β) where β is the inverse temperature. At the mean-field level, we can estimate f non (∆, β) arising from the low-energy spectra around the Dirac points as where Λ is the momentum cutoff. By taking the limit of β → ∞ and then ∆ → 0, we solve the integral, where r 3 = (N −1) 3 U 3 πv 2 > 0. Combined with the analytic part f a , we obtain the total GL free-energy density f (∆) = f a + f non = r 2 ∆ 2 + r 3 |∆| 3 . (C7) According to the GL theory, f (∆) describes a second-order transition with the critical exponent ζ = −1.

APPENDIX D: ABSENCE OF THE PAIRING ORDER
The pairing structure factor can be defined as where 2 × L × L is the number of lattice sites and P (i, j) = α<β c † iα c † iβ c jβ c jα + c iβ c iα c † jα c † jβ is the equal-time pairpair correlation function. By extrapolating the structure factor to thermodynamic limit, the long-range pairing order parameter can be obtained: P s = lim L→∞ 1 2L 2 S pair (L). At half filling, the ground state of the attractive SU(3) Hubbard model on the honeycomb lattice is a semimetal in noninteracting limit U/t → 0. As the attractive interaction increases, the system may enter an ordered phase. At weak coupling, the pairing gap function on the honeycomb lattice is vanishingly small because of the zero density of states at Dirac points [16], and therefore the quantum fluctuations prevent pairings on the half-filled honeycomb lattice. At strong coupling, the attractive SU(2) Hubbard model enters the ordered phase where the pairing order and the CDW order are degenerate, and thus the (spin) SU(2) symmetry is preserved [65,[80][81][82].
As shown in Fig. 11, the finite-size extrapolation to the L → ∞ limit shows no sign of pairing order on the honeycomb lattice in a wide range of coupling strengthes. Hence, the pairing order is absent and therefore the color superfluid is not the ground state of our model. In Appendix C, the meanfield analysis also indicates that the energy of pairing order is higher than that of CDW order on a honeycomb lattice, which is consistent with our QMC results.