Absence of dynamical gap generation in suspended graphene

There is an interesting proposal that the long-range Coulomb interaction in suspended graphene can generate a dynamical gap, which leads to a semimetal-insulator phase transition. We revisit this problem by solving the self-consistent Dyson-Schwinger equations of wave function renormalization and fermion gap. In order to satisfy the Ward identity, a suitable vertex function is introduced. The impacts of singular velocity renormalization and dynamical screening on gap generation are both included in this formalism, and prove to be very important. We obtain a critical interaction strength, $3.2<\alpha_{c}<3.3$, which is larger than the physical value $\alpha = 2.16$ for suspended graphene. It therefore turns out that suspended graphene is a semimetal, rather than insulator, at zero temperature.


Introduction
After the successful fabrication of monolayer graphene in laboratory [1], tremendous experimental and theoretical efforts have been devoted to exploring its novel and fascinating properties [2,3,4,5]. Without any interactions, the low energy quasiparticle excitations of graphene are massless Dirac fermions with linear dispersion [6]. The long-range Coulomb interaction between Dirac fermions is poorly screened due to the vanishing of density of states vanishes at Dirac points, and is therefore anticipated to be responsible for many unusual behaviors of graphene [2,4,5]. Interestingly, it can lead to singular renormalization of fermion velocity [7,8,9,10], which is able to induce unconventional properties in several observable quantities, such as specific heat [11,12], compressibility [12,13], and electrical conductivity [14,15,16,17]. It is remarkable that the predicted singular velocity renormalization has already been observed in recent experiments [18].
If the long-range Coulomb interaction is sufficiently strong, a finite fermion gap may be dynamically generated through the formation of excitonic particle-hole pairs [19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41]. The dynamical gap fundamentally changes the ground state of graphene [42]. As a consequence, there will be a quantum phase transition from semimetal to excitonic insulator. This gap generating mechanism is a concrete realization of the nonperturbative phenomenon of dynamical chiral symmetry breaking that has been extensively investigated in particle physics for five decades [43,44,45,46,47]. From a technological point of view, a finite gap would make graphene more promising as a candidate material for manipulating novel electronic devices [42,48]. For these reasons, the dynamical gap generation and the resultant semimetal-insulator transition have attracted intense theoretical interest in recent years. Many analytical and computational tools, including Dyson-Schwinger(DS) gap equation [19,20,21,22,23,24,25,26,27,28,29], Bethe-Salpeter equation [30,31], renormalization group (RG) [32,33,34], and large scale Monte Carlo simulation [35,36,37,38,39,40], have been applied to study this problem. A critical interaction strength α c and a critical fermion flavor N c are found in almost all these investigations: a dynamical gap is opened only when α > α c and N < N c . If one fixes the physical flavor N = 2, the semimetal-insulator transition is turned solely by the parameter α, with α c defining the quantum critical point. We list some existing values of α c for N = 2 in Table I. With a few exceptions, α c obtained in most calculations is smaller than the physical value α = 2.16 of graphene suspended in vacuum. Therefore, it is suggested by many that the suspended graphene should be an insulator at zero temperature.
However, there is little experimental evidence for the predicted insulating ground state. Actually, a recent experiment put an upper limit, as small as 0.1meV, on the possible gap in suspended graphene [18]. It is known that a finite gap is observed in graphene which is placed on specific substrate [49] or has finite-size configurations [50]. Nevertheless, this gap appears to be induced by the very particular environments, and can not be regarded as a true dynamical gap generated by chiral vacuum condensation. Generally speaking, there may be several reasons for the lack of clear experimental evidence for dynamical gap. For instance, the dynamical gap can be strongly suppressed by thermal fluctuation [19,20,21,22], doping [22], and/or disorders [22,27], which makes it technically hard to measure the dynamical gap in realistic experiments. Another possibility is that the Coulomb interaction is simply not strong enough to open a dynamical gap even in the clean and zero-temperature limits. Motivated by the recent theoretical and experimental progress, we revisit this problem in order to specify the genuine ground state of suspended graphene. Here we only consider clean suspended graphene at zero doping and zero temperature. The gap generation will be examined by means of DS equation, which is known to be a very powerful tool of analyzing dynamical gap generation in a number of strongly interacting models [43,44,45], such as Nambu-Jona-Lasinio model [43,44], quantum chromodynamics (QCD) [44,45], three-dimensional quantum electrodynamics (QED 3 ) [46,47,51], and graphene.
In Ref. [19], Khveshchenko applied the DS equation to examine the possibility of dynamical gap generation in graphene. In order to simplify the very complicated nonlinear gap equation, a number of approximations are introduced. First of all, the energy dependence of the polarization function is neglected, which is usually called instantaneous approximation in literature [19,20,21,22,23,24,26,27,28]. In addition, the renormalizations of wave function and fermion velocity are both ignored. In the subsequent years, the existence of dynamical gap generation in suspended graphene was rechecked after improving some of these approximations [22,23,24,25], and it became clear that these neglected effects can significantly increase or decrease α c . Unfortunately, due to the formal complexity of DS equation(s), it often happens that one specific approximation is improved at the cost of introducing another. In Refs. [23,25], the effect of fermion velocity renormalization on gap generation is investigated and showed to increase α c . However, the dynamical part of polarization function is not well addressed in their calculations. In Ref. [22], a full dynamical polarization function is utilized, but the velocity renormalization is not included. Ref. [24] performs a detailed analysis of the influence of dynamical polarization function without including the velocity renormalization and the energy-dependence of dynamical gap. Another potentially important effect is the strong fermion damping caused by Coulomb interaction. It is predicted that the long-ranged Coulomb interaction, even though being too weak to open a gap, produces marginal Fermi liquid behavior [9,19,52,53]. The strong fermion damping may affect dynamical gap generation, but is rarely considered in the existing literature. Since a precise value of α c is crucial to answer the question regarding the true ground state of graphene, it is necessary to calculate α c after going beyond all these approximations.
We will show that all the aforementioned effects can be incorporated selfconsistently by constructing a set of DS equations of wave function renormalization A 0,1 (p 0 , p) and dynamical fermion gap m(p 0 , p). To satisfy the Ward identity, we also include the vertex correction to the fermion self-energy. The dynamical fermion gap and fermion velocity renormalization can be obtained simultaneously after solving these equations. Our numerical calculations lead to a critical coupling 3.2 < α c < 3.3 for physical flavor N = 2. Apparently, α c is larger than α = 2.16, so the ground state of suspended graphene seems to be a semimetal, rather than an insulator. Our α c is quite different from those obtained in previous DS equation studies, which implies that an appropriate treatment of velocity renormalization and dynamical screening is very important. Moreover, from A 0,1 (p 0 , p), we recover the singular renormalization of fermion velocity.
The rest of the paper is organized as follows. In Sec.2, we construct and numerically solve the self-consistent DS equations. A critical interaction strength α c is obtained from the solutions. In Sec.3, we address the feedback effect of fermion velocity renormalization on the effective Coulomb interaction. We briefly summarize our results and discuss the implications in Sec.4.

Dyson-Schwinger equations for gap and wave functions
The Hamiltonian for interacting Dirac fermions is where the density operator ρ i (r) =ψ i (r)γ 0 ψ i (r) and g = 2πe 2 /v F ε. The Dirac fermions have totally eight indices: two sublattices, two spins, two valleys. As usual, we adopt a four-component spinor field ψ to describe the Dirac fermions [19,20], and define its conjugate field asψ = ψ † γ 0 . Now the physical flavor of fermions is N = 2. The 4 × 4 γ-matrices, γ 0,1,2 , satisfy the Clifford algebra [19,20]. It is easy to check that the total Hamiltonian possesses a continuous chiral symmetry, ψ → e iθγ 5 ψ, where the matrix γ 5 anticommutes with γ 0,1,2 . If the massless Dirac fermions acquire a finite mass gap due to excitonic condensation, ψ ψ = 0, then this continuous chiral symmetry will be dynamically broken.
It is convenient to characterize the interaction strength by an effective fine structure constant, α = e 2 /v F ε. We need to consider the strong-coupling regime since the dynamical gap can only be opened by strong interaction. When α is close to or larger than unity, the conventional perturbation expansion in terms of α breaks down and one has to use 1/N as the expansion parameter. The free propagator for massless Dirac fermion is Due to the Coulomb interaction, this propagator will receive self-energy corrections and become where A 0,1 (p 0 , p) are the wave function renormalization and m(p 0 , p) is the fermion gap function. Here we would like to make a comparison with the gap generation problem of QED 3 . In QED 3 , the explicit Lorentz invariance ensures that A 0 (p 0 , p) = A 1 (p 0 , p) [46,47]. On the contrary, the Lorentz invariance is explicitly broken in the present system, so A 0 (p 0 , p) = A 1 (p 0 , p) and the fermion velocity is renormalized. The renormalized velocity is given by A 1 (p 0 , p)/A 0 (p 0 , p). The function A 0 (p 0 , p) itself is also important since it determines the fermion damping rate caused by the Coulomb interaction. The dynamical fermion gap is described by m(p 0 , p), which acquires a finite value when excitonic particle-hole bound states are formed due to sufficiently strong interaction. In previous DS equation analysis, a common approximation is to simply assume that, A 0 (p 0 , p) = A 1 (p 0 , p) = 1, which drastically simplify the complicated DS equations. However, the effects of velocity renormalization and fermion damping can not be well addressed by doing so. In order to include these effects in a selfconsistent manner, it is better to analyze the coupled equations of A 0 (p 0 , p), A 1 (p 0 , p), and m(p 0 , p). According to the Feynman diagram shown in Fig.1(b), the relationship between the free and renormalized fermion propagators is formally determined by the following DS equation, where Γ 0 (p 0 , p; k 0 , k) is the full vertex function and V (p 0 − k 0 , p − k) is the effective Coulomb interaction function. Regarding the Coulomb potential as the time component of a gauge potential, one can obtain a Ward identity [9] that connects the fermion propagator and the vertex function, Apparently, Γ 0 = γ 0 only when A 0 (p 0 , p) = 1. As we are willing to examine the effects of wave function renormalization on dynamical gap generation, the vertex function Γ 0 can not be simply replaced by the bare vertex γ 0 . In principle, one could build an equation of Γ 0 (p 0 , p; k 0 , k) and couple it self-consistently to Eq.(4). However, this will make the problem intractable. A more practical way is to assume a proper Ansatz.
Choosing a suitable vertex function is a highly nontrivial problem, and have been investigated extensively in the contexts of various quantum field theories [45]. The research experience that has been accumulated in QED 3 [51] is especially helpful since this model is very similar in structure to our present model. There are several frequently used Ansatze for the vertex function [51]. Among these Ansatze, for computational convenience we choose the following one, Other possible Ansatze of Γ 0 (p 0 , p; k 0 , k) can be analyzed by the same procedure. Substituting the full fermion propagator G(p 0 , p) into the above DS equation and then taking trace on both sides, we can obtain the equation for dynamical gap m(p 0 , p). In order to derive the equations of A 0 (p 0 , p) and A 1 (p 0 , p), we multiply both sides of Eq.(4) by γ 0 p 0 and γ · p respectively, and then take trace on both sides. After these manipulations, we finally arrive in a set of self-consistently coupled integral equations, During the derivation of these equations, we have used the following symmetries which originate from the particle-hole symmetry of undoped graphene.
We now consider the effective interaction function V (p 0 − k 0 , p − k). The bare Coulomb interaction is with effective fine structure constant, α = e 2 v F ǫ . The value of dielectric constant ǫ depends on the substrate on which the graphene is placed. For graphene suspended in vacuum, the corresponding value is roughly α = 2.16; for graphene placed on SiO 2 substrate, α = 0.79 [42]. Besides the static screening due to substrate, the Coulomb interaction is also dynamically screened by the collective particle-hole excitations. After including the dynamical screening, as displayed in Fig.(1a), the effective interaction function becomes ., The dynamical polarization function Π(q) is defined as to the leading order, which amounts to random phase approximation (RPA). It is straightforward to obtain Apparently, although the bare Coulomb interaction V 0 (q) is independent of energy, the effective interaction function V (q 0 , q) depend on both energy q 0 and momentum q due to dynamical screening. This V (q 0 , q) should be substituted to the self-consistent equations (7-9). Before doing numerical calculations, it is interesting to compare the present problem with QED 3 . In QED 3 , the inverse of bare photon propagator is ∝ q 2 , with q being three-dimensional energy-momentum. The corresponding polarization is known to be ∝ |q|, which is larger than q 2 at low energies [46]. Therefore, the effective photon propagator is dominated by the polarization and the bare term is relatively unimportant. In the present problem, the inverse of bare Coulomb interaction is V −1 0 (q) ∝ |q|. The polarization Π(q 0 , q) is more complex than that in QED 3 , but is evidently linear in momentum |q|. Therefore, the bare term V −1 0 (q) and the polarization Π(q 0 , q) appearing in the denominator of effective interaction V (q 0 , q) are equally important, which is different from QED 3 .
After time-consuming but straightforward numerical calculations, we obtain the relationship between m(0, 0) and α, shown in Fig.(2). It is difficult to get a precise α c because the necessary computational time becomes increasingly long as α is approaching its critical value. Therefore, here we can only give a narrow range of α c . For physical fermion flavor N = 2, the critical strength is 3.2 < α c < 3.3. This critical value is obviously larger than α = 2.16, so it turns out that the zero-temperature ground state of suspended graphene is a semimetal in the clean limit. This result is consistent with  the fact that so far no experimental evidence for insulating ground state of suspended graphene has been reported.
Our critical value α c is very different from the values presented in earlier publications. Compared with the previous gap equation computations, the key improvements in our analysis are that we have included the wave function renormalization A 0,1 (p 0 , p) self-consistently and maintained the energy-dependence of all the functions of A 0,1 (p 0 , p), m(p 0 , p), and Π(q 0 , q). We also have adopted an Ansatz for the vertex function in order not to violate the Ward identity. Our results indicate that fermion velocity renormalization, fermion damping, and dynamical screening of Coulomb interaction all play important roles in determining the critical interaction strength α c .
It is now useful to make a more concrete comparison between our α c and those presented in some recent literature. For example, our α c is much larger than the α c obtained by Liu et al. [22] and Gamayun et al. [24], who considered a fully dynamical polarization but ignored wave functions and velocity renormalization. Sabio et. al. [25] studied dynamical gap generation by means of a variational method that naturally includes velocity renormalization, and found a critical value α c = 8.1, which is much greater than our α c . Such a large quantitative difference is presumably owing to the instantaneous approximation of Coulomb interaction adopted in their analysis. Once the instantaneous approximation is assumed, the effective Coulomb interaction no longer depends on energy, and V (p 0 + k 0 , p − k) will be identically equal to V (p 0 − k 0 , p − k). As shown in Eq. (7), the time component of wave function becomes A 0 (p 0 , p) = 1. Now the fermion damping effect is automatically neglected and the velocity renormalization is solely determined by A 1 (p 0 , p). Moreover, A 0 (p 0 , p) = 1 implies that vertex function is simply Γ 0 = γ 0 . After these simplifications, our DS equations become similar to those presented in Ref. [25]. By solving these simplified DS equations, we find no evidence for dynamical gap even for infinite coupling α → ∞. This indicates that the instantaneous approximation misses very important fluctuation effects. However, despite the large difference in the magnitude of α c , our conclusion do agree with Sabio et al. that no dynamical gap is generated in suspended graphene.  Although the Coulomb interaction in suspended graphene is too weak to generate a dynamical gap, it still leads to unusual properties: strong fermion damping and singular velocity renormalization. Usually, the fermion damping rate and the fermion velocity renormalization are calculated separately. Our DS equation approach includes the mutual influence of these two quantities self-consistently, so the damping rate and velocity renormalization can be simultaneously obtained from the solutions of A 0 (p 0 , p) and A 1 (p 0 , p). The energy and momentum dependence of A 0 (p 0 , p) and A 1 (p 0 , p) are shown in Fig.(3), with (a-b) for suspended graphene with α = 2.16 and (d-e) for graphene placed on SiO 2 substrate with α = 0.79. It is clear that both A 0 (p 0 , p) and A 1 (p 0 , p) increase as interaction strength α is increasing. A 0 (p 0 , p) decreases with growing energy and momentum. Different from A 0 (p 0 , p), A 1 (p 0 , p) is a decreasing function of momentum |p| but an increasing function of energy p 0 . When |p| decreases gradually, A 1 (p 0 , p) keeps increasing monotonically; when p 0 → 0 or p 0 → ∞, A 1 (p 0 , |p|) saturates to a finite value A 1 (0, p) or A 1 (∞, p), respectively.
The fermion damping is determined by the time component of wave function renormalization, A 0 (p 0 , p). For any given energy p 0 , A 0 (p 0 , p) increases as momentum |p| decreases in the region |p| > p 0 , but saturates to a finite value A 0 (p 0 , 0) as |p| decreases in the region |p| < p 0 . For given momentum |p|, when p 0 decreases gradually, A 0 (p 0 , p) increases with decreasing p 0 in the region p 0 > |p|, but it saturates to a finite value A 0 (0, p) with decreasing p 0 in the region p 0 < |p|. In the low energy regime, is an increasing function of x. In principle, the fermion damping rate can be obtained from A 0 (p 0 , p) by making analytic continuation, ip 0 → ω + iδ. Unfortunately, the function F (Λ/ p 2 0 + p 2 ) is formally quite complicated and can not be written in terms of a simple analytical formula. Our numerical results suggest that F (Λ/ p 2 0 + p 2 ) grows a little more slowly than log(Λ/ p 2 0 + p 2 ) as p 2 0 + p 2 decreases continuously, hence the damping rate does not display an exact linear dependence on energy. By ignoring the singular velocity renormalization, some pervious perturbative [9,52] and self-consistent [19,53] calculations predict a marginal Fermi liquid behavior, namely Γ(ω) ∼ ω. Such linear-in-energy behavior disappears once the singular velocity renormalization is selfconsistently considered. However, the renormalization factor Z extracted from our numerical results is found to vanish, Z = 0, so the suspended graphene does not have well-defined quasiparticles. In this sense, our results are qualitatively consistent with the previous conclusions [19,52,53].
The fermion velocity renormalization is related to many properties of graphene [5,7,8,10,11,12,13,14,15,16,17]. Such effect is usually studied by RG method [7,8,9,10]. It can be naturally obtained by solving the self-consistent DS equations. The ratio between renormalized and bare fermion velocities is defined as We plot v R F (p 0 , p) for α = 2.16 and α = 0.79 in (c) and (f) of Fig.(3), respectively. The renormalized fermion velocity v R F (p 0 , p) increases singularly as |p| → 0, which reproduces the previous results of RG analysis [7,8,9,10]. We also obtain the energy dependence of renormalzied velocity v R F (p 0 , p), which is not well captured in RG analysis. For given momentum, v R F (p 0 , p) is an increasing function of energy. When p 0 → 0 or , as can be seen from Fig.(3c) and Fig.(3f) as well as Fig.(4). It is also interesting to consider the gapped phase with α > α c , though no known graphene material manifests such a large α. When α = 4.0, a dynamical gap is generated by the Coulomb interaction. The corresponding m(p 0 , p), A 0,1 (p 0 , p), and v R F (p 0 , q) are shown in Fig.(5). The gap m(p 0 , p) decreases with growing momentum but increases with growing energy. These behaviors are very different from those in QED 3 where the dynamical gap is a decreasing function of both energy and momentum. This difference reflects the different energy-dependence of Coulomb interaction and gauge interaction. From Eq.(13), we know that the effective Coulomb interaction decreases with growing momentum but increases with growing energy. However, the gauge interaction in QED 3 always decreases as energy or momentum grows, which is in accordance with the fact that QED 3 exhibits asymptotic freedom [46,47].
In the low-energy regime, p 2 0 + |p| 2 < m(0, 0), A(p 0 , p) saturates to a finite value and can be approximate as A 0 (p 0 , p) ∼ 1+F (Λ/ p 2 0 + |p| 2 + m 2 (0, 0)) with F (x) being certain increasing function of x. After analytic continuation, ip 0 → ω + iδ, the fermion damping rate is found to vanish when |ω| < |p| 2 + m 2 (0, 0). It implies that the fermions can not be excited below the scale of m(0, 0). Unlike the case of semimetal phase, v R F (p 0 , p) does not keep increasing as momentum decreases. For any given energy, v R F (p 0 , p) saturates to v R F (p 0 , 0) below the scale |p| ∼ m(0, 0), so the singular velocity renormalization is suppressed by the dynamical gap generation in the insulating phase, as shown in Fig.(4). These results are consistent with those obtained by perturbative calculations [54] and by functional RG calculations [55].

Effects of velocity renormalization on polarization
In the calculations presented above, we have used the RPA expression of the dynamical polarization function Π(q 0 , q), in which the fermion velocity is a constant. However, this function may receive sizeable feedback effects from the velocity renormalization. Since the effective Coulomb interaction plays a crucial role in determining α c , it deserves to examine these effects. In a fully self-consistent analysis, one should consider the following defination where G(k 0 , k) is the full fermion propagator given by Eq.(3) and Γ 0 is the vertex function. By doing so, the polarization Π(q 0 , q) couples self-consistently to the equations of A 0,1 (p 0 , p) and m(p 0 , p). Unfortunately, this would make the computational time unacceptably long since the integrations over energy and momenta must be performed separately in the present problem. To simplify the calculations, we assume that the RPA expression of Π(q 0 , q) provides a reliable description of the dynamical screening effect, which then allows us to replace the constant velocity v F appearing in Π(q 0 , q) by the renormalized velocity v R F (p 0 , p). This strategy is analogous to that employed in Ref. [23]. After this substitution, the effective interaction becomes It is larger than 1 since v R F (q 0 , q) > v F , so the velocity renormalization tends to enhance the effective interaction strength. In order to examine the extend to which this feedback effect changes α c , we re-solve the coupled gap equations after substituting Eq.(18) into Eqs. (7 -9). We find the fermion gap is always zero when α ≤ 2.9. It implies that such feedback effect does not change α c substantially. In particular, α c obtained after considering this effect is still larger than the physical value α = 2.16 in suspended graphene.

Summary and discussions
In conclusion, we have presented a detailed analysis of dynamic gap generation due to Coulomb interaction in graphene using the coupled DS equations of wave function renormalization and fermion gap. After including the effects of fermion velocity renormalization and dynamical screening of Coulomb interaction, we obtain a critical interaction strength, 3.2 < α c < 3.3. It is apparently greater than the physical strength α = 2.16 in suspended graphene. Therefore, the Coulomb interaction in suspended graphene is too weak to generate a dynamical fermion gap, and the semimetal ground state of suspended graphene is robust, which is consistent with the fact that so far no insulating phase has been observed in experiments.
Although our calculations suggest that the Coulomb interaction in suspended graphene is not strong enough to open a dynamical gap, the possibility of dynamical fermion gap generation can not be entirely precluded. Indeed, there do exist several possible mechanisms for the dynamical generation of fermion gap. For instance, the Dirac fermions can acquire a tiny bare gap for some reasons, such as Kekule distortion [56] and spin-orbit coupling [57]. When this happens, a large dynamical fermion gap can be generated by a relatively weak Coulomb interaction [26], which then produces properties analogous to the remarkable phenomena of QCD [26]. Moreover, the additional on-site repulsive interaction may help to generate a dynamical gap even if the Coulomb interaction itself is not sufficiently strong [22,24,27]. Finally, it is well known that an external magnetic field perpendicular to the graphene plane can lead to dynamical gap generation and insulating phase even at infinitesimal coupling [58,20].