Charge Smoothening and Band Flattening due to Hartree corrections in Twisted Bilayer Graphene

TL;DR

This paper demonstrates that Coulomb interactions in doped twisted bilayer graphene smooth out charge inhomogeneities, flatten the electronic bands, and influence the material's correlated phases, using real-space Hartree calculations.

Contribution

It introduces a self-consistent Hartree approach on a relaxed lattice to show how Coulomb interactions modify charge distribution and band structure in doped twisted bilayer graphene.

Findings

Charge inhomogeneities are smoothed by Coulomb interactions.

Flat bands become further flattened with doping.

Van Hove singularity remains pinned at the Fermi level.

Abstract

Doping twisted bilayer graphene away from charge neutrality leads to an enormous buildup of charge inhomogeneities within each Moir\'{e} unit cell. Here we show, using unbiased real-space self-consistent Hartree calculations on a relaxed lattice, that Coulomb interactions smoothen this charge imbalance by changing the occupation of earlier identified `ring' orbitals in the AB/BA region and `center' orbitals at the AA region. For hole doping, this implies an increase of the energy of the states at the ${\bf \Gamma}$ point, leading to a further flattening of the flat bands and a pinning of the Van Hove singularity at the Fermi level. The charge smoothening will affect the subtle competition between different possible correlated phases.

Tables

Table 1. Table 1: The total accumulated charge in the AA region, δ ​ n AA 𝛿 subscript 𝑛 AA \delta n_{\mathrm{AA}} , and the maximum electric potential difference Δ ​ ϕ Δ italic-ϕ \Delta\phi , as a function of the filling ν 𝜈 \nu . We compare the results before, and after the self-consistent Hartree calculation. Both the electric potential and the charge inhomogeneity are severely reduced by the Hartree potential. The reduction factor (before divided by after) is shown in the last column.

$\nu$	$\delta n_{\mathrm{AA}}$		$\Delta \varphi$ (meV)		Reduction
	before	after	before	after	factor
0	-0.026962	-0.0049907	28.8934	8.60985	3 - 5
-1	-0.318072	-0.0285706	339.343	30.5505	11
-2	-0.61295	-0.0515452	660.412	56.6246	12
-3	-0.87059	-0.0641868	935.592	67.5796	14
-4	-1.03138	-0.0850411	1088.83	91.2731	12

Equations

$$-t(\vec{d})=V_{\pi }^0\left[1-{\left(\frac{\vec{d}\cdot {\vec{e}}_z}{d}\right)}^2\right]e^{-\frac{d-a_0}{r_0}}+V_{\sigma }^0{\left(\frac{\vec{d}\cdot {\vec{e}}_z}{d}\right)}^2{e}^{-\frac{d-d_0}{r_0}}$$

$$H_V=\frac{1}{2}\sum _{ij}\delta n({\mathbf{r}}_i)V({\mathbf{r}}_i-{\mathbf{r}}_j)\delta n({\mathbf{r}}_j)$$

$$V({\mathbf{r}}_i-{\mathbf{r}}_j)=\frac{1.438}{0.116+|{\mathbf{r}}_i-{\mathbf{r}}_j|}\mathrm{eV}$$

$$H_H$$

$${\varphi }_i$$

Full text

Studies of twisted bilayer graphene have revealed a number of striking electronic properties,Trambly de Laissardière et al. (2010); Luican et al. (2011); Brihuega et al. (2012); Wong et al. (2015); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011); Lopes dos Santos et al. (2012); Slotman et al. (2015); van Wijk et al. (2015); Fang and Kaxiras (2016); Nam and Koshino (2017) and recently culminated in the discovery of a correlated insulator and superconducting phase.Cao et al. (2018a, b) Several new experiments performed since thenCao et al. (2019); Yankowitz et al. (2019); Choi et al. (2019); Xie et al. (2019); Jiang et al. (2019); Lu et al. (2019); Sharpe et al. (2019) have revealed new correlated and topological states, and even more theoretical works have appearedXu and Balents (2018); Po et al. (2019); Kang and Vafek (2018, 2019); Zou et al. (2018); Po et al. (2018); Xie and MacDonald (2018); Rademaker and Mellado (2018); Pizarro et al. (2018); Wu et al. (2018); Kennes et al. (2018); Yuan and Fu (2018); Efimkin and MacDonald (2018); Liu et al. (2019); Choi and Choi (2018); Fu et al. (2018); Koshino et al. (2018); Venderbos and Fernandes (2018); Tarnopolsky et al. (2019); Song et al. (2018); Sherkunov and Betouras (2018); Thomson et al. (2018); Ochi et al. (2018); Zhang et al. (2019); Wu et al. (2019); Seo et al. (2019); Laksono et al. (2018); Guinea and Walet (2018); Cea et al. (2019); Guinea and Walet (2019); Haule et al. (2019); Bultinck et al. (2019); Leaw et al. (2019); Yan et al. (2019); Walet and Guinea (2019); Chou et al. (2019); González and Stauber (2019); Li et al. (2019); Wu and Das Sarma (2019); Gu et al. (2019); Yudhistira et al. (2019). However, there is no consensus yet on the correct low-energy effective band-structure, let alone the relevant many-body interactions.

It is well known that at fillings away from charge neutrality the electronic charge piles up in the AA regions of the Moiré unit cell.Trambly de Laissardière et al. (2010); Luican et al. (2011); Brihuega et al. (2012); Wong et al. (2015); Rademaker and Mellado (2018); Guinea and Walet (2018, 2019) The Coulomb interaction will counterbalance the formation of such large charge inhomogeneities, and is expected to significantly modify the non-interacting band dispersion. Here we study the effect of the electronic interactions on AA localization via an unbiased, self-consistent Hartree calculation.

Starting from a microscopic tight-binding model on a relaxed twisted bilayer at the magic angle $\theta=1.08^{\circ}$ at fillings at and below charge neutrality, we show that the Hartree corrections indeed cause the redistribution of charges within the unit cell. The resulting charge distribution is much smoother than that found with noninteracting theories. The smoothening itself is made possible by increasing the energy of the states at the ${\bf\Gamma}$ point in the Brillouin zone, which causes a further substantial flattening of the (already nearly) flat bands. Our results are qualitatively similar to earlier results based on a continuum model.Cea et al. (2019); Guinea and Walet (2018) Given the tight competition between different possible correlated phases,Xu and Balents (2018); Ochi et al. (2018); Kang and Vafek (2019) the charge transfer reported here will almost certainly affect this subtle competition.

In the remainder of this paper, we first introduce our model of twisted bilayer graphene and the details of our calculation. We then present and compare the charge inhomogeneities before and after the Hartree corrections. Finally, we discuss the ramifications of our result for future theoretical modelling.

The model – Our starting point is a tight-binding model of twisted bilayer graphene, Rademaker and Mellado (2018) which builds on earlier continuum models.Bistritzer and MacDonald (2011) We construct twisted bilayer graphene by starting with an AB stacked bilayer and rotating the top layer around an AB site. For a commensurate twist angle, we choose two integers $m_{1},m_{2}$ such that the large Moiré unit cell has unit vectors ${\bf G}_{1}=m_{1}{\bf a}_{1}+m_{2}{\bf a}_{2}$ and ${\bf G}_{2}=-m_{2}{\bf a}_{1}+(m_{1}+m_{2}){\bf a}_{2}$, where the graphene lattice unit vectors are ${\bf a}_{i}=\tfrac{a}{2}\left((-1)^{i}\hat{\bf x}+\sqrt{3}\hat{\bf y}\right)$ with $a=0.246$ nmReich et al. (2002); Castro Neto et al. (2009). The twist angle is now given by $\cos\theta=\frac{m_{1}^{2}+4m_{1}m_{2}+m_{2}^{2}}{2(m_{1}^{2}+m_{1}m_{2}+m_{2}^{2})}$, with $4(m_{1}^{2}+m_{1}m_{2}+m_{2}^{2})$ being the number of atoms per unit cell. In the remainder of this paper we focus on $m_{1}=30$, $m_{2}=31$, which corresponds to the twist angle $\theta=1.08^{\circ}$ and $N=11,164$ atoms in the Moiré unit cell. A schematic picture of the unit cell is shown in Fig. 1, left.

At such small twist angles the atomic positions are substantially changed when one allows for lattice relaxation. To include this effect, we use the relaxed atomic positions from Ref. [Gargiulo and Yazyev, 2018].111We thank Prof. Yazyev for kindly providing us with a list of atomic positions in the relaxed lattice. Using these atomic positions, we define a tight-binding hopping model with hopping between two atoms given by

[TABLE]

Here $a_{0}=0.142$ nm is the intralayer nearest neighbor atomic distance; $V^{0}_{\pi}=2.7$ eV is the nearest neighbor intralayer hopping strength; $d_{0}=0.335$ nm is the unrelaxed interlayer distance; $d=|\vec{d}|$ is the length of vector that connects the two carbon atoms; $r_{0}=0.045$ nm is the $p$-orbital decay length; and $V^{0}_{\sigma}=0.48$ eV is the interlayer hopping when two carbon atoms are exactly above each other. Nam and Koshino (2017); Trambly de Laissardière et al. (2010).

The electrons are subject to Coulomb interactions, which depend on the deviation of the electron density $n(\mathbf{r})$ from the average density $\bar{n}$. With $\delta n(\mathbf{r})\equiv n(\mathbf{r})-\bar{n}$ the interaction Hamiltonian reads

[TABLE]

where the sum $i,j$ runs over all pairs of atomic positions $\mathbf{r}_{i}$, $\mathbf{r}_{j}$. We follow estimates for single-layer graphene and interpolate between an onsite repulsion of $V(0)=9.3$ eV and an unscreened $\frac{e^{2}}{|\mathbf{r}|}$ at large distances. This leads to the following effective form of the screened Coulomb interaction in the tight-binding model,Wehling et al. (2011)

[TABLE]

where the distance between two carbon atoms $|\mathbf{r}_{i}-\mathbf{r}_{j}|$ should be measured in nm.222Eqn. (3) is an interpolation between the short-range coupling in Table I of Ref. [Wehling et al., 2011], and the long-range result $\epsilon\rightarrow 1$ for $r\rightarrow\infty$ given in the same article.

In the Hartree approximation, the product of two density operators in Eqn. (2) is replaced by the product of an operator and its expectation value. The interaction is now replaced by a site-dependent electric potential $\phi_{i}$, that is to be determined self-consistently,

[TABLE]

This provides us with a self-consistent iterative scheme: using the hopping Hamiltonian – including the Hartree potential of Eqn. (4) – we compute the band structure, from which we deduce the electronic density on each atom, which allows us to compute the electric potentials $\phi_{i}$ using Eqn. (5). We repeat these steps until we reach convergence, and the fields no longer change.

The band-structure and Bloch wavefunctions were computed on a $8\times 8$ Monkhorst-Pack momentum grid enlarged to reflect the sixfold rotational symmetryMonkhorst and Pack (1976), which means in total we have 192 momentum points. At each momentum point we compute the full band-structure, and then adjust the chemical potential to find the desired filling. We use mixing of the old and new electric potential to speed up convergence. After 20 iterations the electric potential had converged, in all cases, to within an accuracy of $10^{-4}$ eV. Note that we explicitly exclude lattice symmetry breaking in our computations, in contrast to other works focusing on, for example, $C_{3}$ breaking.Choi et al. (2019)

We introduce the filling parameter $\nu$, which gives the number of electrons per Moiré unit cell, relative to the charge neutrality (where the average charge density per atom is $\bar{n}=1$). Note that $\nu=\pm 4$ corresponds to either completely empty or completely filled flat bands, since there are $8$ of them, taking into account spin and valley indices. The average density per atom $\bar{n}$ is related to the filling parameter $\nu$ via $\bar{n}=1+\nu/N$, where $N$ is the number of atoms in the unit cell. In the remainder, we will focus on the fillings $\nu=0,-1,-2,-3$ and $-4$.

Charge distribution before self-consistency – We will first demonstrate the strong inhomogeneity of the charge distribution and electric potential arising in the non-interacting model, without self-consistent Hartree corrections. Interestingly, we find that the inhomogeneity in AB/BA regions is present even at charge neutrality.333Earlier we assumed that the charge distribution is smooth at charge neutrality,Rademaker and Mellado (2018) however, this turns out to be wrong The distribution of charge at $\nu=0$, illustrated in Fig. 1, right, reveals a clear charge deficiency $|\delta n|\lesssim 0.0026$ on the atoms that lie on top of an atom from the other layer. Note that such inhomogeneities exist in AB stacked bilayers only away from charge neutrality.444A quick calculation of AB stacked bilayers shows that at charge neutrality, the A and B sublattice have the same charge density. Away from charge neutrality, the charge on the sublattices is different. Furthermore, we find that the AB/BA regions have no net accumulated charge.

Upon hole doping, a significant charge inhomogeneity arises between the AA and the AB/BA regions of the unit cell. This is because of the well-documented fact that the local density of states is highest at the AA region.Trambly de Laissardière et al. (2010); Luican et al. (2011); Wong et al. (2015); Guinea and Walet (2018) The relative charge density is shown in Fig. 2, top row. Even though the charge density per atom is $|\delta n|<0.001$, the combined effect of thousands of atoms leads to a charge accumulation. We can express this by a total accumulated charge in the AA region, $\delta n_{\mathrm{AA}}$, obtained by summing the charge over the $N/2=5582$ atoms closest to the AA center. The values for $\delta n_{\mathrm{AA}}$ are shown in Table 1.

This charge inhomogeneity leads to a very strong electric potential, illustrated in Fig. 3, left. Small potential fluctuations arise because of the charge oscillations in the AB/BA region as shown in Fig. 1. At fillings away from charge neutrality, however, a major potential difference arises between the AA center and the AB/BA regions of the unit cell. For filling $\nu=-4$ this difference exceeds 1 eV. A full list of the potential difference $\Delta\phi$ as a function of filling is given in Table 1.

Self-consistent results – As shown in Fig. 3, the electric potential is negative around the AA regions. The potential therefore tends to pull more electrons there, thus negating the charge inhomogeneity. Consequently, the Hartree corrections will smoothen the charge distribution within each unit cell. Earlier we qualitatively discussed this phenomenon, which can be viewed as a charge transfer between AA and AB/BA regions.Rademaker and Mellado (2018) Here we provide a quantitative analysis and discuss its effect on the band-structure.

We find that the Hartree corrections reduce both the charge inhomogeneities and the electric potential by an order of magnitude at fillings other than charge neutrality, as shown in Fig. 2, Fig. 3 and Table 1.

The charge redistribution affects dispersion of the bands, as well as the local density of states (LDOS). In Fig. 4 we show the LDOS, measured at the AA region of the unit cell, for various fillings. At charge neutrality there are symmetric peaks representing the Van Hove singularities (VHS). When partially filling the flat bands (here the fillings $\nu=-1,-2$ and $\nu=-3$) the VHS remains close to the Fermi level, consistent with recent findings a continuum model.Cea et al. (2019)

The change of position of the VHS is a direct consequence of band structure changes. Recall that in the non-interacting picture the Bloch wavefunctions at the ${\bf K}$ point are highly centered at the AA regions of the unit cell (center orbitals),Guinea and Walet (2018); Rademaker and Mellado (2018). On the other hand, the real-space structure of the flat bands at the ${\bf\Gamma}$ point has a vanishing weight at the AA centers and instead resembles a ring occupying the AB/BA regions. Charge smoothening can be achieved by transferring charge from the center to the ring orbitals. Consequently, for hole doping away from charge neutrality ($\nu<0$), the charge smoothening will push up the energy of the band at the ${\bf\Gamma}$ point. This can be clearly seen in Fig. 5, where we show the band structure after including the Hartree corrections.

The states at the ${\bf\Gamma}$ point are clearly pushed upwards, to the point that at $\nu=-4$ they are higher in energy than the states at the ${\bf K}$ point. This is true both for the flat bands as well as the bands lower in energy. However, for the flat bands this upwards push directly causes an even further flattening, in the form of a reduction of the bandwidth from $18\,{\rm meV}$ to roughly $6\,{\rm meV}$ at its extreme. Especially noteworthy is the filling $\nu=-2$, where the charge smoothening leads to the appearance of an almost completely flat band segment between the ${\bf\Gamma}$ and ${\bf M}$ points. It will be interesting to investigate the effect of this significant flattening on the correlated insulator state observed at this filling.Cao et al. (2019, 2018a); Yankowitz et al. (2019)

Discussion – We have shown that Coulomb interactions, as captured by a fully self-consistent Hartree calculation, smoothen the charge inhomogeneities in the unit cell. This charge smoothening is achieved by pushing up the energy of the states at the ${\bf\Gamma}$ point, which in turn forces the VHS to remain close the Fermi level.

Note that in this work we assumed an interaction strength similar to those expected for single-layer graphene. However, both screening within the bilayer as well as the dielectric environmentPizarro et al. (2019) can alter the precise electron-electron interactions. Nevertheless, a different shape of the Coulomb interaction will only quantitatively change our results as its tendency to smooth out charge inhomogeneities remains. Similar quantitative changes are expected when one includes the lattice relaxation response to the Hartree potential, which is not included in this work.

Both the redistribution of charge and the shift in VHS position upon doping can, in principle, be observed using scanning tunnelling microscopy (STM). Indeed, a sharpening of the VHS closest to the Fermi level has been observed,Choi et al. (2019); Xie et al. (2019); Jiang et al. (2019) but further many-body effects, not considered here, make a quantitative comparison between our results and the observed STM spectra difficult.

Indeed, to quantitatively understand the plethora of interesting observed phases, ranging from superconductivity, correlated insulator states, quantum anomalous Hall effect, ferromagnetism, requires the inclusion of many-body interactions beyond Hartree corrections.Cao et al. (2018a, b, 2019); Yankowitz et al. (2019); Choi et al. (2019); Xie et al. (2019); Jiang et al. (2019); Lu et al. (2019); Sharpe et al. (2019); Ochi et al. (2018); Kang and Vafek (2019) This, in turn, demands the computation of localized Wannier orbitals and their Hubbard, exchange and Hund couplings. Given the strong band flattening and charge smoothening the Hartree correction gives, we argue that any effective low-energy model should start from a Hartree-renormalized band structure such as the one shown in Fig. 5.

Note added. Upon finalization of this manuscript, a related paper appeared.Cea et al. (2019) This work starts from a continuum model without lattice relaxation, while we start from a tight-binding model with lattice relaxation. Otherwise, the resulting band-structure and LDOS spectra are qualitatively similar. A major difference with Ref. [Cea et al., 2019] is that we discuss the renormalization of the band structure in terms of the charge inhomogeneities and the resulting electric fields, while authors of Ref. [Cea et al., 2019] interpret their results in terms of VHS pinning to the Fermi level.