Strong Interaction Effects at a Fermi Surface in a Model for Voltage-Biased Bilayer Graphene

Monte Carlo simulation of a 2+1 dimensional model of voltage-biased bilayer graphene, consisting of relativistic fermions with chemical potential mu coupled to charged excitations with opposite sign on each layer, has exposed non-canonical scaling of bulk observables near a quantum critical point found at strong coupling. We present a calculation of the quasiparticle dispersion relation E(k) as a function of exciton source j in the same system, employing partially twisted boundary conditions to boost the number of available momentum modes. The Fermi momentum k_F and superfluid gap Delta are extracted in the limit j tends to zero for three different values of mu, and support a strongly interacting scenario at the Fermi surface with Delta of order O(mu). We propose an explanation for the observation mu


Introduction
There are very few many-body systems permitting Monte Carlo simulation without the need to confront a Sign Problem. In Ref. [1] we introduced a new member to this class, based on an effective theory of bilayer graphene in which charge carrying excitations are modelled as N f = 4 relativistic fermions moving in a 2d plane. The introduction of chemical potential µ is via a bias voltage in the perpendicular direction which induces equal densities of electrons on one layer and holes on the other; this is analogous to isospin chemical potential in QCD and yields a real positive fermion determinant amenable to orthodox Monte Carlo methods. The simulation was performed in the vicinity of a quantum critical point (QCP) found at strong coupling, and the main result was the demonstration that the ground state is a superfluid formed by condensation of electronhole exciton pairs, and that the response to chemical potential is governed by the noncanonical scaling forms (7) and (8) given in Sec. 2.
The model [1] as originally devised for bilayer graphene [2] is artificial in a few respects. Firstly, the description in terms of N f = 4 relativistic species (ie. N f = 2 electrons and N f = 2 holes) is only justified by the band structure of the tight-binding model in the presence of an inter-layer "skew" coupling breaking the trigonal symmetry of the underlying lattice [3]. Secondly, the interaction between charge densities is simplified to be a local four-fermi contact, although a more realistic unscreened Coulomb interaction can be modelled with the introduction of a third spatial lattice direction to capture the electrodynamics [4]. Finally, intra-and inter-layer interactions have the same coupling strength, as a necessary condition of keeping the fermion determinant real. Nonetheless, it shares the essential features of a more general model for double-layer graphene systems in which there is some hybridization permitting interlayer tunnelling [5]. With N f = 1 this approach is also applicable to surface states of topological insulators, motivating study with variable N f [6].
The results of [1] were interpreted in terms of strong interaction effects at a Fermi surface; if exciton pairs within a shell of thickness ∆ condense around the Fermi surface centred at k F , then the anomalous scaling (7,8) is consistent with a BCS mechanism with ∆ ∼ O(µ). Everything is to be viewed in the context of an effective field theory valid near the QCP. In order to put this picture on firmer footing, and also to expose the Fermi surface, in this paper we use Monte Carlo simulation to calculate the quasiparticle dispersion relation E(k), identifying k F with the location of the minumum and ∆ with E(k F ). The main results are summarised in Fig. 5 below. In Sec. 2 we present the model and review the main findings of [1], then in Sec. 3 present the calculation of E(k). Our results, summarised in Sec. 4, indeed support the picture of a Fermi surface disrupted by strong interactions, leading to the formation of a gap ∆ increasing monotonically with µ. We also discuss our observation of the striking inequality µ < k F , and propose an explanation in terms of an estimate for the dynamical critical exponent z < 1.

Formulation and Simulation of the Model
The model we use to describe voltage-biased bilayer graphene in terms of N f = 4 relativistic fermions is described by the following Lagrangian [1]: Here ψ and φ are 4×2-component spinors each describing two Dirac flavors, with ψ,ψ corresponding to electron degrees of freedom on one layer and φ,φ holes on the other; V is an auxiliary field defined on the timelike links of the lattice which approximates an "instantaneous" Coulomb potential governed by 3+1d Maxwell electrodynamics; and j a symmetry-breaking gap parameter due to interlayer pairing. Because some weak interlayer hybridization is likely to be present in double-layer systems, in general j = 0 [5]; however we will attempt to extrapolate j → 0 so that exciton condensation can be viewed as a spontaneous symmetry breaking U(4)⊗U(4)→U(4) [1]. The bias voltage is given by 2µ where µ is formally equivalent to the isospin chemical potential in QCD.
In continuum notation the covariant derivative operator is where α, β run over N f = 2 Dirac flavors. The minimal coupling to V implies that ψψ, φφ and φψ interactions are all of equal strength, which is required for the action to be real following integration over the fermions. This corresponds to the interlayer separation d → 0 in the double-layer model [5].
In terms of staggered fermion fields living on the sites x, y of a 2+1d cubic lattice D is written (3) where the sign factors η νx ≡ (−1) x 0 +···x ν−1 ensure a covariant weak-coupling continuum limit. Eq. (3) was also used to model electron excitations in monolayer graphene [7]. Note that in 2+1d a single staggered fermion automatically describes N f = 2 continuum flavors [8]. We work in units in which both spatial a s and temporal a t lattice spacings are set to unity (equivalent to setting the bare Fermi velocity v F = 1); note that for a non-covariant action there is no reason a priori to assume a s = a t , though since the value of the ratio a t /a s is driven by UV physics it does not depend on µ and only weakly on j, so may be assumed constant throughout this paper. On the assumption that the dimension of D latt is even, it is straightforward to show and hence there is no obstruction to Monte Carlo simulation using orthodox numerical techniques.
In [1] we presented results from numerical simulation of the model (1,3) using a hybrid Monte Carlo algorithm. The coupling g −2 a ′ = 0.4, where the factor a ′ = a 2 s a −1 t follows because the interaction couples charge densities, was chosen in the vicinity of the quantum critical point (QCP), although since N f = 4 < ∼ N f c = 4.8(2) [6] it is hard to ascertain with confidence on which side of the phase boundary at µ = 0 we are sitting. Two principal observables were monitored as a function of µ and j, namely the carrier density and the exciton condensate Following extrapolation to the limit j → 0 (so long as µa t < 0.3 at which point saturation artifacts set in and the continuum approximation fails), both observables showed behaviour consistent with rising smoothly from zero as µ is increased, with This is to be contrasted with the expectations from weak-coupling. In free field theory the carrier density depends on the volume contained within the Fermi surface; for relativistic fermions k F ≈ µ and hence in 2+1d n c ∝ µ 2 . Similarly, in weak coupling we expect the exciton condensate to arise from electron-hole pairing with equal and opposite momenta from within a shell of thickness 2∆ centred on k F ; hence ΨΨ ∝ ∆µ. The non-canonical scaling is taken to be a symptom of strong field fluctuations near the QCP. Moreover, since according to Luttinger's theorem n c ∝ k 2 F even in the presence of interactions, we can adapt this argument to estimate the scaling of the gap ∆(µ): While the numerical value for the exponent should probably not be taken too seriously, the qualitative difference between the scaling predicted in weak coupling and observed near the QCP is striking; indeed, Fig. 13 of [1] shows the ratio ΨΨ /n 1 2 c extrapolated to j → 0 almost linearly proportional to µ. Since near a QCP µ is the only energy scale in the problem, naively ∆ ∝ µ is expected. As both arguments contain assumptions, it is clear that a direct calculation of ∆ from the quasiparticle propagator Ψ(0)Ψ(x) is needed.

Quasiparticle Dispersion
While the results outlined in the previous section are intriguing the conclusions, drawn from simulations of systems with both significant UV and IR artifacts on a restricted range of µ values, are necessarily provisional. In this paper we will present complementary information through analysis of the quasiparticle dispersion relation E( k), which will enable identification of the Fermi momentum k F and direct estimation of the superfluid gap ∆. The basic observables are the timeslice correlators in momentum space: where we distinguish between normal propagation of an electron or hole within a layer, and anomalous propagation in which an electron moving in one layer is absorbed by an exciton, transferring its momentum to an electron moving in the other. On a finite volume in the absence of explicit symmetry breaking j = 0, the anomalous component C A necessarily vanishes. Many momentum modes must be available in order to obtain a good resolution for E( k). Naively this would entail making at least one of the spatial dimensions of the lattice as large as possible [9]. It is much more efficient, however, to employ partially twisted boundary conditions [10] in which the constraint ψ(L x ) = e iθx ψ(0), with the angle θ x adjustable, is implemented in the calculation of the propagator (10) so that accessible modes have In practice the twist is implemented via a field redefinition so that each x-link in (3) with ε(x) ≡ (−1) x 0 +x 1 +x 2 , means that C N (k, t) vanishes for t even and C A (k, t) for t odd. Moreover only ReC N and ImC A survive the ensemble average. The resulting correlators are then fitted to the forms [11] (Lt−t) ). (14) to yield the k-dependent amplitudes A, B, |C| and the energy E, which is extracted independently from both C N and C A . Stable fits for C N (t) and C A (t − 1) were found for the windows ta t ∈ [7, 25] (32 3 ) and ta t ∈ [7, 41] (48 3 ). It is slightly unusual to have to deal with excited state contamination when fitting a fermion propagator, which may be a symptom of the strong fluctuations near a QCP; this was found to be an issue particularly for k > ∼ k F .  It is instructive first to consider the fitted amplitudes: Fig. 1 shows results from µa t = 0.1 on 48 3 , for two values of the exciton source j. In the normal channel the timeasymmetric form of C N is manifest, and changes in character as k increases. For small k the forwards propagating signal is stronger, but B/A grows with k and for ka s > ∼ 0.2 the backwards signal dominates. The interpretation is as follows [11]: for k < k F the dominant excitations are hole-like, and for k > k F particle-like. Increasing jã from 0.005 to 0.05 has the effect of smearing the Fermi surface so that quasiparticles tend to become an admixture of both, and the disparity between A and B diminishes. The effect of smearing is also seen in the anomalous channel, where C A grows steadily in magnitude with increasing j. In weakly-coupled models [11,9] |C(k)| is non-monotonic with a maximum near k F where particle-hole mixing is strongest, but here the behaviour is less clear-cut. The same trends with both k and j are observed at larger µ.  Fig. 2 shows results from fits to the dispersion E(k) for various µ at jã = 0.01 from (13,14). The common feature is that E(k) is non-monotonic with a minimum in the neighbourhood where the amplitude ratio A/B ≈ 1, which we have identified as the Fermi momentum k F . For k < k F quasiparticle excitations are hole-like, and the energy needed to excite them from the ground state decreases as k ր k F . For k > k F , excitations are particle-like and the opposite holds true: indeed, in this regime results from all three µ-values studied are plausibly consistent with being drawn from the same branch of the dispersion curve appropriate to the vacuum (ie. with zero bias voltage). Fig. 2 compares results from two volumes 32 3 and 48 3 , and also for µa t = 0.15 with fits in both normal (13) and anomlaous (14) channels. On the assumption that a smooth curve may be drawn through the admittedly noisy data, there is no evidence for any significant finite volume artifacts, or systematic difference between the two channels. Fig. 3 plots E(k) for various j at µa t = 0.1 on 48 3 . The non-monotonicity observed above becomes more pronounced as j → 0, and since there is little shift in the minimum with j, we adopt the pragmatic procedure of identifying the Fermi momentum k F with the value of k where E is minimum. The resulting estimates are shown in Table 1; the quoted error is half the mode spacing on 48 3 , except at µa t = 0.2 where the dispersion is flatter and a full mode spacing is taken. The first observation is that k F is systematically greater than µ, consistent with the precocious saturation of n c (µ) observed in [1]. This is discussed further in Sec. 4, but already we note a discrepancy with the expectation k F = µ for free massless fermions with a s = a t . Setting aside the issue of the smearing of the Fermi surface by an exciton gap ∆ > 0, we can at this stage test consistency with Luttinger's theorem, which states that n c depends solely on the geometry of the Fermi surface characterised by k F , independent of the nature of the interactions. The third   column of Table 1 shows the carrier density n free c evaluated for free massless fermions on the same 48 3 lattice at the reference value of µ, and the fourth with µ set equal to k F in the second column; the fifth column gives the value of n c in the interacting theory in the limit j → 0 obtained in [1]. In all cases n free c (µ) ≪ n free c (k F ), while the values of n free c (k F ) and n c (µ) are comparable, which is encouraging; however as µ increases n free c (k F ) < n c (µ) indicative of the difficulties in precisely locating k F due to both the limited momentum resolution and perhaps the absence of a sharp Fermi surface due to exciton condensation.
In Fig. 4 we plot the gap ∆ defined as the energy E(k F ). Data from all available fits on both volumes is shown, in both normal and anomalous channels. As before there is little evidence for a systematic effect with lattice volume and both channels yield consistent results. The dashed lines show an extrapolation to j = 0 based on the quadratic form ∆ = ∆ 0 + aj + bj 2 .
It is clear lim j→0 ∆(j) = 0 which is direct evidence for spontaneous gap formation via exciton condensation for all three values of µ. The resulting extrapolations, together with the estimates for k F (µ) from Table 1 are shown in Fig. 5, which is the main result of this paper. In the region of µ studied it is clear that ∆ varies strongly with and is of the same order of magnitude as the chemical potential µ, which are notable features of this particular model proposed in [1], and indicative of strong interactions at the Fermi surface. This conclusion holds in both normal and anomalous channels, appears to be independent of volume, and is striking enough to be robust against uncertainties introduced by the ad hoc nature of our analysis (an analytic form against which to fit the dispersion E(k) would be very valuable), the IR and UV artifacts inherent in studies of lattice models with µ = 0 [12], and lack of knowledge of the physical anisotropy a t /a s . On the basis of the three chemical potentials studied both ∆ and k F appear to scale superlinearly with µ, in qualitative agreement with (9).

Discussion
In this paper we have used lattice Monte Carlo simulation techniques to explore the quasiparticle dispersion relation in an interacting field theory with non-zero charge density, and shown that for k ∼ k F the excitations are gapped with ∆ ∼ O(µ), and ∆ scaling faster than linearly with µ. This is in sharp contrast to results from comparable studies of other simulable models with µ = 0. In [9] the gap in the 3+1d Nambu Jona-Lasinio model (a relativistic analogue of the original BCS model) was shown to be approximately constant above onset, independent of and numerically much smaller than µ, consistent with the BCS result ∆ ∼ Λ U V exp(−cΛ 2 U V /µ 2 ). In QCD with gauge group SU(2) there is a so-called quarkyonic regime above onset where ΨΨ ∝ µ 2 [12]; this is consistent with degenerate fermions in 3+1d with a gap ∆ ∼ O(Λ QCD ) independent of µ. It is also very different from the result ∆/µ ∼ O(10 −7 ) obtained by self-consistent diagrammatic  Table 1. techniques [13], although comparable with the large values of ∆/µ obtained in [5], where it was found that ∆ depends sensitively on the treatment of screening effects, and in particular on the reduction of screening once the superfluid gap forms. One feature of our approach which does merit comparison with the treatment in [5] is that competition between inter-and intra-layer pairing condensates can be addressed; see Fig. 11 of [1].
We now return to Table 1 and the issue of why µa t < k F a s . In [1] it was suggested this is because in a strongly self-bound system the Fermi energy is necessarily less than the Fermi momentum in natural units. While this may be plausible for a system where µ ≫ ∆ is by far the largest scale, it is difficult to see how this picture can persist in the regime we have been focussing on. Another possibility, which we cannot dismiss, is that the disparity is a lattice artifact caused by a large induced anisotropy a t /a s ∼ O(0.5). Indeed, if we assume the Fermi velocity remains close to one even in the presence of interactions, then the dispersion data of Fig. 2 might suggest a t /a s ≈ 0.3. However, a more compelling possibility is that E(k) is not a linear relation, but rather a power law characteristic of a nearby QCP. Rewrite the scaling form (7) as n c ∝ µ 2 z ; the Luttinger scaling n c ∝ k 2 F then gives E ∝ k z , where guided by Fig. 2 we assume the relation between E F and k F completely characterises the quasiparticle dispersion. In this scenario the scaling (7) extracted from bulk observables in [1] thus yields an estimate for the dynamical critical exponent z ≈ 0.6.
Much greater numerical precision than achieved in Fig. 2 would be needed to distinguish these two possibilites unambiguously. Another route would be to perform a "biased bilayer" study for a related 2+1d theory, the Thirring model [14], whose be-haviour as a function of N f and µ is qualitatively similar to the model here [15], but whose continuum action is manifestly covariant implying a t ≡ a s throughout.
Finally, we note that the superlinear scaling of ∆(µ) is also suggestive of a power law ∆ ∝ µ σ , with σ > 1, modifying our naive expectation ∆ ∝ µ. Clearly, a much more extensive simulation campaign is required to verify this interesting possibility.