Electronic Ground State in Bilayer Graphene with Realistic Coulomb Interactions

TL;DR

This study uses quantum Monte Carlo simulations with realistic Coulomb interactions to determine the ground state of bilayer graphene, confirming the necessity of a finite onsite interaction for insulating behavior and analyzing the influence of interlayer coupling.

Contribution

It provides non-perturbative evidence that a finite onsite interaction induces an insulating ground state in bilayer graphene with realistic Coulomb interactions.

Findings

A finite critical onsite interaction is required to open a gap in bilayer graphene.

The critical onsite interaction decreases as interlayer coupling increases.

Results align with renormalization-group predictions for Hubbard models.

Abstract

Both insulating and conducting electronic behaviors have been experimentally seen in clean bilayer graphene samples at low temperature, and there is still no consensus on the nature of the interacting ground state at half-filling and in the absence of a magnetic field. Theoretically, several possibilities for the insulating ground states have been predicted for weak interaction strength. However, a recent renormalization-group calculation on a Hubbard model for charge-neutral bilayer graphene with short-range interactions suggests the emergence of low-energy Dirac fermions that would stabilize the metallic phase for weak interactions. Using a non-perturbative projective quantum Monte Carlo, we calculate the ground state for bilayer graphene using a realistic model for the Coulomb interaction that includes both short-range and long-range contributions. We find that a finite critical onsite interaction is needed to gap bilayer graphene, thereby confirming the Hubbard model expectations even in the presence of a long-range Coulomb potential, in agreement with our theoretical renormalization group analysis. In addition, we also find that the critical onsite interactions necessary to destabilize the metallic ground state decreases with increasing interlayer coupling.

Equations

$$H$$

$$R=1-\frac{S\left(\mathbf{\Gamma }+\mathbf{b}/L\right)}{S\left(\mathbf{\Gamma }\right)}$$

$${\mathcal{L}}_0=\left(\begin{array}{cccc}-i\omega & 0& 0& -t{\gamma }_{\vec{k}}^{\ast }\\ 0& -i\omega & -t{\gamma }_{\vec{k}}& 0\\ 0& -t{\gamma }_{\vec{k}}^{\ast }& -i\omega & t_{\perp }\\ -t{\gamma }_{\vec{k}}& 0& t_{\perp }& -i\omega \end{array}\right),$$

$$\Sigma (\vec{k})=\int \frac{d\omega d^2\vec{q}}{(2\pi {)}^3}\frac{2\pi e^2}{\left|\vec{k}-\vec{q}\right|}{G}_0\left(\omega ,\vec{q}\right),$$

$$\Sigma \left(\vec{k}\right)\approx$$

$$\times \left(\begin{array}{cccc}0& -\frac{{\left({\gamma }_{\vec{q}}^{\ast }\right)}^2}{{|{\gamma }_{\vec{q}}|}^2}& 0& -2{\gamma }_{\vec{q}}^{\ast }\frac{t}{t_{\perp }}\\ -\frac{{\gamma }_{\vec{q}}^2}{{|{\gamma }_{\vec{q}}|}^2}& 0& -2{\gamma }_{\vec{q}}\frac{t}{t_{\perp }}& 0\\ 0& -2{\gamma }_{\vec{q}}^{\ast }\frac{t}{t_{\perp }}& 0& 1\\ -2{\gamma }_{\vec{q}}\frac{t}{t_{\perp }}& 0& 1& 0\end{array}\right),$$

$$\frac{q}{\sqrt{k^2+q^2-2kq\cos ({\theta }_{qk})}}$$

$$=1+\frac{\cos ({\theta }_{qk})}{q}k+\frac{3{\cos }^2({\theta }_{qk})-1}{2q^2}{k}^2+\mathcal{O}\left(k^3\right)$$

$$\frac{{\left({\gamma }_{\vec{q}}^{\ast }\right)}^2}{{\left|{\gamma }_{\vec{q}}\right|}^2}=$$

$$+q^2{e}^{2i{\theta }_q}\frac{e^{6i{\theta }_q}-1}{48}+\mathcal{O}\left(q^3\right),$$

$$\Sigma \left(\vec{k}\right)=\left(\begin{array}{cccc}0& \alpha & 0& 0\\ {\alpha }^{\ast }& 0& 0& 0\\ 0& 0& 0& \beta \\ 0& 0& {\beta }^{\ast }& 0\end{array}\right)$$

$${\alpha }_1$$

$${\alpha }_2$$

$${\beta }_0$$

$${\beta }_2$$

$${\mathcal{H}}_0\left(\vec{k}\right)+\Sigma \left(\vec{k}\right)\approx \left(\begin{array}{cccc}0& {\alpha }_{1}k{e}^{-i\theta }+{\alpha }_2{k}^2{e}^{2i\theta }& 0& \left(\frac{1}{2}\sqrt{3}{e}^{i\theta }k-\frac{1}{8}{e}^{-2i\theta }{k}^2\right)t\\ {\alpha }_{1}k{e}^{i\theta }+{\alpha }_2{k}^2{e}^{-2i\theta }& 0& \left(\frac{1}{2}\sqrt{3}{e}^{-i\theta }k-\frac{1}{8}{e}^{2i\theta }{k}^2\right)t& 0\\ 0& \left(\frac{1}{2}\sqrt{3}{e}^{i\theta }k-\frac{1}{8}{e}^{-2i\theta }{k}^2\right)t& 0& t^{\prime }+{\beta }_0+{\beta }_2{k}^2\\ \left(\frac{1}{2}\sqrt{3}{e}^{-i\theta }k-\frac{1}{8}{e}^{2i\theta }{k}^2\right)t& 0& t^{\prime }+{\beta }_0+{\beta }_2{k}^2& 0\end{array}\right).$$

Full text

I Introduction

Stacking one layer of graphene on top of another dramatically changes the dispersion. In the simplest consideration where, in addition to the intra-plane hopping, only the hopping between the carbon atoms that are placed directly on top of the other is allowed, the non-interacting low-energy dispersion changes from linear to quadratic. Although the stacking does not directly gap out the bands, it does enhance the density of states at the band-touching point, which increases quantum fluctuations and the likelihood of correlated ground states. Moreover, onsite electron interactions are marginally relevant at half-filling, i.e. arbitrarily weak interactions will lead to an instability of the semi-metallic phase Sun et al. (2009). Bilayer graphene has additional degrees of freedom, which include spins, sublattices, layers and valleys, which has led to several theoretical predictions for different competing insulating phases including symmetry breaking of either real spin or pseudospin degrees of freedom (Castro et al., 2008; Min et al., 2008; Zhang et al., 2010, 2011; Jung et al., 2011; Zhang and MacDonald, 2012), quantum anomalous Hall states (Nandkishore and Levitov, 2010), nematic states (Vafek, 2010; Vafek and Yang, 2010; Lemonik et al., 2010, 2012) and canted antiferromagnetic states (Kharitonov, 2012). All of these theoretical works focused on the low-energy circular symmetric limit, where the dispersion is parabolic.

However, experiments remain inconclusive regarding what the ground state of bilayer graphene at low temperature is, with zero electric and magnetic field, and whether this system is conducting or gapped (Feldman et al., 2009; Weitz et al., 2010; Martin et al., 2010; Mayorov et al., 2011; Velasco et al., 2012; Freitag et al., 2012; Ulstrup et al., 2014). In particular, in the works by Bao et al. and Freitag et al. (Bao et al., 2012; Freitag et al., 2012), both insulating and conducting phases were observed in similarly prepared samples. The metallic samples were found to have low temperature conductivity of around $2$–$4\,e^{2}/h$, while that of the insulating samples was one order of magnitude lower. This sample-dependent differences in the experimental data remains unresolved.

In a recent theoretical development, Pujari et al. (Pujari et al., 2016) considered higher energies away from the band-touching point where the Hamiltonian is no longer circular symmetric, and showed that a linear term can be generated by contact interactions. Here, dimensional considerations show that the linear term once generated renders the onsite interaction irrelevant, and as such the system flows to a stable fixed point with Dirac cones. Qualitatively, we can understand the stability of this Fermi liquid as resulting from the vanishing density of states of the emergent Dirac Fermions. This conclusion was supported by quantum Monte Carlo simulations (Pujari et al., 2016) which show within this model, that for weak interactions, the system remains metallic, contradicting the previous expectation that interactions are marginally relevant. Pujari et al. expected that their findings would not hold for realistic long-range Coulomb interactions since the same dimensional analysis shows that long-range interactions are relevant for this system, and as a result lead to an instability of the Fermi liquid for vanishing interactions.

However, as shown in Ref. Tang et al., 2018, in the realistic system where the long-range interactions are weaker than the contact interactions, the long-range interactions will stabilize the Fermi liquid instead of destabilizing it. One way to understand this is that the long-range interactions favor the charge-density-wave phase, which competes with the antiferromagnetic tendency favored by the contact interaction Hohenadler et al. (2014). From the point of view of renormalization group, the Coulomb interaction in 2D systems has a non-analytical form of $2\pi e^{2}/\left|\mathbf{q}\right|$. Since the Wilson’s renormalization group cannot generate non-analytic terms, the Coulomb interaction vertex does not get directly renormalized (González et al., 1994); due to Ward identity, however, the renormalization to the Coulomb interaction is tied to the renormalization of the fermionic propagator. In the case of Dirac fermions where the fermionic propagator is described solely by the Fermi velocity $v_{F}$, the renormalization of the long-range interactions is captured by a single parameter $r_{s}=e^{2}/v_{F}$. Similarly, it has also been shown for systems with parabolic dispersions that the renormalization group flow leads to a modification of the dynamical exponent, which prevents the long-range interaction coupling constant to run away (Janssen and Herbut, 2017). This can physically be understood as follows: Within the Hartree-Fock approximation, the effective mass of the bilayer graphene is renormalized by long-range interactions to a reduced level (Borghi et al., 2009; Viola Kusminskiy et al., 2009). Since the strength of the long-range interactions is determined by the parameter $r_{s}=2me^{2}/\sqrt{\pi n}$ in bilayer graphene (Das Sarma et al., 2011), the reduced effective mass is equivalent to a weaker interaction. Therefore, similarly to monolayer graphene where the enhanced Fermi velocity reduces the long-range interaction strength, the downward renormalization of the effective mass prevents the long-range interaction to escape to strong coupling.

The renormalization group arguments presented above are confirmed by large scale unbiased qauntum Monte Carlo simulation of electrons on Bernal stacked bilayer honeycomb lattice with on-site and long-range Coulomb interactions. Nevertheless, one needs to be careful in interpreting the QMC results. For non-interacting electrons on the bilayer honeycomb lattice with finite interlayer hopping $t_{\perp}$, the electronic bands interpolate between a parabolic dispersion close to the band touching point, and a linear dispersion at high momenta. Such a crossover happens at the momentum $ka=t_{\perp}/\sqrt{3}t$ (McCann and Fal’ko, 2006). However, limited by finite system sizes, quantum Monte Carlo simulations cannot probe momenta arbitrarily close to the band touching point. For a system with $t_{\perp}=t$, the crossover between a linear and parabolic dispersion occurs at $ka=1/\sqrt{3}$, which is not much larger than the smallest momentum that can be accessed in simulations. The quantum Monte Carlo simulations with nominally realistic values of inter-layer hopping are hence at the scales corresponding to the linear part of the band, where a large finite critical on-site interaction is expected. Instead, one needs to use unrealistically large inter-layer hopping to probe the effects of on-site and long range interactions on the low-energy quadratic dispersion.

In this work, using numerically exact, projective quantum Monte Carlo Sugiyama and Koonin (1986); Sorella and Tosatti (1992); Hohenadler et al. (2014); Assaad and Evertz (2008), we examine closely the argument that a linear term is dynamically generated in bilayer graphene. Our numerics suggests that the linear term is visibly generated only in the system with $t_{\perp}\gg t$. In addition, our numerics support the conclusion that the phase transition occurs at a finite onsite interaction strength even in the presence of long-range interactions. The rest of the paper is organized as follows: In Sec. II, we define the model and the parameters, while in Sec. III we present the quantum Monte Carlo results on the phase transition of bilayer graphene for multiple values of the inter-layer hopping, with $t\leq t_{\perp}\leq 10t$. The renormalized spectrum is examined in Sec. IV. Finally, we discuss the relevance of our results to the realistic system in Sec. V.

II Model

Our quantum Monte Carlo simulations assume single band interacting electrons on a bilayer honeycomb lattice at half-filling. In our model, the intralayer lattice vector $\mathbf{a}_{1}=a\hat{x}$, $\mathbf{a}_{2}=a\hat{x}/2+\sqrt{3}a\hat{y}/2$ and the interlayer vector $\mathbf{a}_{3}=1.38a\hat{z}$, where $a$ is the lattice constant $2.46\,\text{\AA}$ (Jung and MacDonald, 2014). Each unit cell consists of 4 sites $A_{1},B_{1},A_{2}$ and $B_{2}$, from the two sub-lattices of the two layers, with $B_{1}$ sites sitting on top of the $A_{2}$ sites. The Hamiltonian with both the onsite Hubbard interactions and the long-range Coulomb interactions takes the form

[TABLE]

The operator $a_{i\sigma l}^{\dagger}$ $\left(a_{i\sigma l}\right)$ acts on the atom in sub-lattice $A$ of the $l$-th layer, located in the unit cell positioned at $\mathbf{r}_{i}$, to create (annihilate) an electron of spin $\sigma=\uparrow\downarrow$. Similarly, $b_{i\sigma l}^{\dagger}$ and $b_{i\sigma l}$ act on sub-lattice $B$. The number operator $\hat{n}_{i}=\sum_{\sigma}\hat{n}_{i\sigma}=\sum_{\sigma}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}$ counts the number of electrons sitting at the atom $i$. For realistic bilayer graphene, while the intra-layer hopping integral is $t=2.6\text{--}3.1\text{ eV}$, the interlayer hopping $t_{\perp}$ is very small, usually about $10\%\text{--}15\%$ of $t$ (Jung and MacDonald, 2014). The interacting part of the Hamiltonian consists of the onsite Hubbard interaction $U$ and long-range Coulomb tail $V_{ij}=\gamma U/r_{ij}$, where $r_{ij}$ is the distance between atom $i$ and atom $j$. The coupling constant for long-range interactions $\alpha$ is related to the parameter in our model $\gamma$ as $\alpha=2U\gamma/3\sqrt{3}$. We ignore other hopping integrals such as interlayer hopping between atoms from the same sub-lattice which are negligible compared to $t$ and $t_{\perp}$.

III Phase transition of bilayer graphene

When we approach the strong coupling limit in the Hubbard model, the double occupancy at each site is suppressed, and the system maps to the $S=1/2$ Heisenberg model. Depending on the ratio of $t_{\perp}/t$, the system may develop a long-range Néel order or turn into a dimer phase. The interplay between the antiferromagnetic order and the dimer phase will be presented elsewhere Leaw et al. (2019). To study the phase transition, we measure the structure factor that measures both the antiferromagnetic order and the dimer order $S\left(\mathbf{q}\right)=\frac{1}{2L^{2}}\sum_{ij}e^{i\mathbf{q}\cdot\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}\left\langle M_{\text{AF}}\left(\mathbf{r}_{i}\right)\cdot M_{\text{AF}}\left(\mathbf{r}_{j}\right)\right\rangle$, where $M_{\text{AF}}\left(\mathbf{r}_{i}\right)=\sum_{l}\left[m_{Al}\left(\mathbf{r}_{i}\right)-m_{Bl}\left(\mathbf{r}_{i}\right)\right]$, and $m_{Cl}=\hat{c}_{i\uparrow l}^{\dagger}\hat{c}_{i\uparrow l}-\hat{c}_{i\downarrow l}^{\dagger}\hat{c}_{i\downarrow l}$ measures the magnetization of the $l$-th layer $C=A,B$ atom in unit cell positioned at $\mathbf{r}_{i}$. The correlation function $\left\langle M_{\text{AF}}\left(\mathbf{r}_{i}\right)\cdot M_{\text{AF}}\left(\mathbf{r}_{j}\right)\right\rangle$ is computed using the numerically exact zero-temperature projective quantum Monte Carlo method (Bercx et al., 2017; Tang, 2017), which projects the wave-function to zero temperature to obtain the ground state properties of the system. The resulting structure factor $S$ has an anomalous dimension, which would require an additional exponent when we perform the scaling analysis. We use instead the correlation ratio $R$ of the structure factor at the ordering momentum $\bm{\Gamma}$ and the momentum closest to it,

[TABLE]

where $\mathbf{b}$ is the reciprocal lattice vector, and $\mathbf{b}/L$ is the momentum closest to the $\bm{\Gamma}$ point for a system with size $L\times L$. Since the correlation ratio $R$ is dimensionless, it scales to 1 in the antiferromagnetic phase and scales to 0 in the metallic phase (see Fig. 1). Exactly at the phase transition, the spin structure factor is scale invariant, which allows us to determine the critical interaction strength by pinpointing the Hubbard interaction that shows no change in the correlation ratio $R$ when the system size is increased.

The zero-temperature projective quantum Monte Carlo method has recently been adapted to include the long-range Coulomb interaction (Brower et al., 2012; Ulybyshev, 2013; Hohenadler et al., 2014). This allows us to study the phase transition in systems with long-range interactions. We perform the scaling analysis for various values of $\gamma$, and obtain a critical Hubbard interaction $U_{c}$ for each value of $\gamma$. Tracing out the $U_{c}$ in the phase space of onsite and long-range interaction, we obtain the phase boundaries as shown in Fig. 2a. We perform this analysis for various values of $t_{\perp}$, and compare with the monolayer results (Tang et al., 2018).

In the canonical Hubbard model, we find that when we increase $t_{\perp}$ from [math] to $10t$, the critical onsite interaction strength $U_{c}$ decreases gradually from that of the monolayer graphene. The critical onsite interaction strength $U_{c}$ is plotted against the interlayer hopping $t_{\perp}/t$ in Fig. 2b. For the case $t_{\perp}=t$, we find $U_{c}\approx 3t$ in agreement with the results in Ref. Pujari et al., 2016. However, the critical value for $t_{\perp}=t$ is close to that of monolayer graphene ($t_{\perp}=0$). This suggests that for $t_{\perp}=t$, the parabolic part of the low-energy dispersion is confined to momenta that are inaccessible to the finite system sizes studied, and instead, the QMC simulations yield results for the linear part of the spectrum. In other words, the system behaves more like a monolayer graphene than a bilayer graphene. On the other hand, for interlayer hopping as large as $t_{\perp}=10t$ where the electronic band is parabolic at scales accessible to the simulations (see Fig. 3), the critical value of the onsite interaction strength remains small but finite. In this limit, the QMC results describe the effects of on-site interactions on the parabolic band, and confirm that the a finite, non-zero critical $U_{c}$ is required to open up a gap at the Fermi surface. This supports the argument that a linear term is generated to stabilize the Fermi liquid.

When the long-range interaction is turned on, the critical onsite interaction strength increases. This can be understood in the following way. As the linear term is allowed by symmetry to be dynamically generated, the long-range interaction also contributes to the linear term during renormalization group flow. In addition, the mass of the parabolic term is also reduced by renormalization (Borghi et al., 2009; Viola Kusminskiy et al., 2009). Both of these amount to the increase in the renormalized energy spectrum, which decreases the density of states at the Fermi level. Therefore, a stronger critical onsite interaction is needed to gap out the system. To show that this is the case, we study the energy spectrum renormalization directly in the next section.

IV Spectrum renormalization

In addition to the phase transition, the zero-temperature projective quantum Monte Carlo method can also be used to study the spectrum renormalization, as demonstrated in Ref. Assaad and Herbut, 2014; Tang et al., 2018 for monolayer graphene. Using the projective quantum Monte Carlo method, we calculate the time-displaced single-particle imaginary time Green’s function $G_{\mathbf{k}}\left(\tau\right)=\sum_{\sigma}\left\langle c_{\mathbf{k},\sigma}^{\dagger}\left(\tau\right)c_{\mathbf{k},\sigma}\left(\tau=0\right)\right\rangle$. In the limit of $\tau\rightarrow\infty$, the Green’s function has only a single exponential decay $G_{\mathbf{k}}\left(\tau\rightarrow\infty\right)=Z_{\mathbf{k}}e^{-E_{\mathbf{k}}\tau}$, where $Z_{\mathbf{k}}$ corresponds to the single-particle residue and $E_{\mathbf{k}}$ is the renormalized single-particle energy. By fitting our quantum Monte Carlo results at large $\tau$ to an exponential form, we can extract the renormalized energy $E_{\mathbf{k}}$.

Fig. 4 shows the energy renormalization $E-E_{0}$ for the Hubbard model $\gamma=0$ and system with the long-range interaction $\gamma=1$, each for the case of $t_{\perp}=t$ and $t_{\perp}=10t$. In the system with large interlayer hopping $t_{\perp}=10t$, we cannot exclude that a linear term is dynamically generated. Close to the Dirac point the data are consistent with $E(\mathbf{k})=v_{f}|\mathbf{k}|+|\mathbf{k}|^{2}/2m$ with $v_{f}>0$ such that the density of states at small frequencies reads: $N(\omega)=\frac{2\pi}{v_{f}}|\omega|$. Together with the result that the phase transition for $t_{\perp}=10t$ occurs at a finite value, quantum Monte Carlo data support the renormalization group argument proposed by Pujari et al Pujari et al. (2016). On the other hand, we see a negative energy renormalization in the case of the Hubbard model with $t_{\perp}=t$, but a positive energy renormalization in the case with long-range interaction. A similar trend is shown in the energy renormalization in monolayer grapheneTang et al. (2018). This further supports the view that for $t_{\perp}=t$ the QMC simulations on accessible system sizes can only probe the linear part of the spectrum and the system behaves more like a monolayer than a bilayer graphene.

We note that for large interlayer hopping $t_{\perp}=10t$, the system shows a larger positive energy renormalization in the presence of long-range Coulomb interactions. This indicates that the long-range interactions contribute to generating the linear term, alongside with the contact interactions. Here we show that this is possible within the first-order perturbation theory. Considering the tight-binding model for Bernal-stacked bilayer honeycomb lattice with all the four bands, the non-interacting part of the Lagrangian density is $L_{0}=\psi^{\dagger}\mathcal{L}_{0}\psi$, with matrix structureMcCann and Fal’ko (2006)

[TABLE]

where $\gamma_{\vec{k}}=\sum_{i}e^{i\left(\vec{K}+\vec{k}\right)\cdot\vec{\delta}_{i}}$ is the sum of the phase factors, $\vec{K}=(4\pi/3,0)$ is one of the Dirac points, and $\vec{\delta}_{i}$ are the positions of the nearest neighbours. In the first-order approximation, the self-energy due to the long-range interaction is

[TABLE]

where the non-interacting Green’s function is given by the inverse of the Lagrangian density $G_{0}=\mathcal{L}_{0}^{-1}$. After integrating out the frequency $\omega$, we keep the terms that are dominating in the limit of large interlayer hopping $t/t_{\perp}\rightarrow 0$,

[TABLE]

where $\theta_{qk}\equiv\theta_{q}-\theta_{k}$. In the limit of small external momentum $k$,

[TABLE]

and approximating $\left(\gamma_{\vec{q}}^{*}\right)^{2}/\left|\gamma_{\vec{q}}\right|^{2}$ to second order in $q$,

[TABLE]

we integrate out the angle $\theta_{q}$ to obtain

[TABLE]

where $\alpha=\alpha_{1}ke^{-i\theta_{k}}+\alpha_{2}\left(ke^{i\theta_{k}}\right)^{2}+\mathcal{O}\left(k^{3}\right)$ and $\beta=\beta_{0}+\beta_{2}k^{2}+\mathcal{O}(k^{3})$. Note that only certain terms of $k$ remains, the other terms vanish upon integration over $\theta_{q}$. The subscripted $\alpha$ and $\beta$ are integrals of $q$,

[TABLE]

The integrals are evaluated in the limit of $ka\ll 1$ and the ultraviolet cutoff $\Lambda a$ is of the order of 1. In this limit, $\alpha_{1}$, $\alpha_{2}$, $\beta_{0}$ and $\beta_{2}$ have some real finite values. We may now find the renormalized energies by taking the eigenvalues $\epsilon_{n}\left(\vec{k}\right)$ of $\mathcal{H}_{0}\left(\vec{k}\right)+\Sigma\left(\vec{k}\right)$, where $\mathcal{H}_{0}=\mathcal{L}_{0}+i\omega I_{4}$. Keeping only terms up to $k^{2}$, we have

[TABLE]

The contribution $\sim\beta$ of the self-energy would just renormalize the quadratic bands. However, the contribution $\sim\alpha$ to the self-energy implies that a linear term is generated by the long-range interactions.

V Discussion

Using complimentary RG analysis and large scale unbiased QMC simulations, we have shown unambiguously that the parabolic low-energy dispersion of non-interacting electrons on Bernal stacked bilayer honeycomb lattice is renormalized by the dynamic generation of a linear dispersion due to the effects of on-site interactions. This linear dispersion stabilizes the Fermi liquid phase and an interaction-driven transition to a Mott insulating state occurs at a finite, non-zero critical value that depends on the magnitude of the inter-layer hopping. Long range interactions further enhance the linear dispersion and drive the critical interaction strength to even larger values. Our results further demonstrate that QMC simulations on bilayer honeycomb lattice needs to be treated with caution. This is exemplified by the results for the interacting model at small to intermediate values of the inter-layer hopping ($t_{\perp}\leq t$). While the results yield a non-zero critical $U_{c}$, qualitatively in agreement with expectations from RG analysis, data for spectrum renormalization at low energies reveal that this is an artefact of finite system sizes accessible to QMC simulations. The parabolic dispersion of the non-interacting model is confined to small momenta $ka\lesssim t_{\perp}/\sqrt{3}t$. For $t_{\perp}\lesssim t$, to probe the relevant momenta, one needs to go to length scales larger than the system sizes that are accessible to QMC simulations. As a result, the QMC results in this limit correspond to the linear (at larger momenta) part of the dispersion and effectively reproduces the monolayer physics.

The effects of the dynamically generated linear term in the dispersion can be probed by QMC simulations at large values of the interlayer coupling $t_{\perp}\gg t$. The parabolic dispersion extends to higher momenta that are accessible to available system sizes in our study and the results are relevant to large length scale physics.

Despite the large discrepancy between the interlayer hopping in the simulations $t_{\perp}=10t$ and the experimental interlayer hopping $t_{\perp}\approx 0.1t$, we can still glean information from the simulations to understand the experiments. More specifically, the renormalization group flow filters out the high energy physics to focus on the low energy physics. During a renormalization group flow, the high frequency modes are integrated out and the resulting effective theory is rescaled to be compared with the original theory. Often in this renormalization group approach, a high energy cutoff is introduced to specify the validity regime of the theory, where the effects of renormalization beyond this cutoff are neglected. One can imagine to start a renormalization group flow from realistic bilayer graphene parameter $t_{\perp}\approx 0.1t$. During the flow only the interlayer hopping $t_{\perp}$ is rescaled, all other renormalization effects are neglected. At the point when the interlayer hopping is rescaled to $t_{\perp}=10t$, we arrive at a new theory for which is cut off at $t_{\perp}=10t$, and that we study using the quantum Monte Carlo simulations. Although the interaction will also be rescaled to a different value, the conclusion that the phase transition occurs at finite critical onsite interaction strength remains valid.

Additionally, in realistic bilayer graphene where the interlayer next-nearest-neighbour hopping is allowed, a linear dispersion appears in the momentum scale below the parabolic dispersion. Taken together, both these arguments suggest a finite critical interaction for the transition to a Mott insulating phase in realistic graphene bilayers.

The work was made possible by allocation of computational resources at the CA2DM (Singapore), National Supercomputing Centre in Singapore and the Gauss Centre for Supercomputing (SuperMUC at the Leibniz Supercomputing Center), and funding by the Singapore Ministry of Education (MOE2017-T2-1-130), Deutsche Forschungsgemeinschaft (SFB 1170 ToCoTronics, project C01) and NSERC of Canada.