Externally controlled band gap in twisted bilayer graphene
A. O. Sboychakov A. V. Rozhkov A. L. Rakhmanov Franco Nori
cond-mat.str-el





Center for Emergent Matter Science, RIKEN, Wako-shi, Saitama,
351-0198, Japan
Institute for Theoretical and Applied Electrodynamics, Russian
Academy of Sciences, Moscow, 125412 Russia


Center for Emergent Matter Science, RIKEN, Wako-shi, Saitama,
351-0198, Japan
Institute for Theoretical and Applied Electrodynamics, Russian Academy of Sciences, Moscow, 125412 Russia
Moscow Institute for Physics and Technology (State
University), Dolgoprudnyi, 141700 Russia


Center for Emergent Matter Science, RIKEN, Wako-shi, Saitama,
351-0198, Japan
Moscow Institute for Physics and Technology (State
University), Dolgoprudnyi, 141700 Russia
Institute for Theoretical and Applied Electrodynamics, Russian Academy of Sciences, Moscow, 125412 Russia
Dukhov Research Institute of Automatics, Moscow, 127055 Russia


Center for Emergent Matter Science, RIKEN, Wako-shi, Saitama,
351-0198, Japan
Department of Physics, University of Michigan, Ann Arbor, MI
48109-1040, USA



We theoretically study the effects of electron-electron interaction in twisted bilayer graphene in applied transverse dc electric field. When the twist angle is not very small, the electronic spectrum of the bilayer consists of four Dirac cones inherited from each graphene layer. Applied bias voltage leads to the appearance of two hole-like and two electron-like Fermi surface sheets with perfect nesting among electron and hole components. Such a band structure is unstable with respect to exciton band gap opening due to the screened Coulomb interaction. The exciton order parameter is accompanied by the spin-density-wave order. The value of the gap depends on the twist angle. More importantly, it can be controlled by applied bias voltage which opens new directions in manufacturing of different nanoscale devices.


73.22.Pr, 73.21.Ac




























Externally controlled band gap in twisted bilayer graphene
    Franco Nori
    December 30, 2023
==========================================================



Introduction—
It is known that application of the bias voltage to AB stacked bilayer graphene opens a gap in its electronic spectrum <cit.>. This feature makes bilayer graphene promising for applications in electronics. Experiment shows, however, that, in many cases, the structure of the bilayer graphene samples is different from a simple AB stacking, and is characterized by a non-zero twist angle between layers <cit.>. Twisting makes the physics of the bilayer graphene more complicated and rich. For example, it leads to the appearance of the moiré patterns – alternating dark and bright regions seen in STM images <cit.>. For a countable set of twist angles the system has a superstructure with the period, which may coincide or be a commensurate with the moiré period <cit.>.

Rotational misorientation affects also the electronic properties of the bilayer.  Analysis in the single-particle approximation shows that one can distinguish three qualitatively different types of behavior of the spectrum at low energies. When the twist angle θ is close to commensurate value corresponding to the superstructure with considerably small size of the supercell, the spectrum has a gap at Fermi level <cit.>. This gap, however, is very sensitive to small variations of the twist angle, and is non-negligible only for a limited number of superstructures <cit.>. With the exception of those values, one can assume that when θ is greater  than critical value θ_c∼1–2^∘, the electron spectrum has a linear dispersion and consists of four Dirac cones inherited from two graphene layers. For a commensurate structure, which we only study in this paper, the four Dirac points of two graphene layers are distributed between two non-equivalent corners of the superlattice Brillouin zone, forming the band structure with two doubly-degenerate Dirac cones in the corners of the superlattice Brillouin zone. The Fermi velocity of these Dirac cones, however, turns out to be θ-dependent: it decreases monotonically down to zero when θ goes to critical value <cit.>. Finally, when θ<θ_c, the spectrum at low energies is characterized by flat bands and the density of states has a peak at the Fermi level <cit.>.

Despite the progress in understanding the electronic properties of twisted bilayer graphene, many issues still remain unclear. First of all this concerns many-body effects, since in the majority of cases, theorists are limited to a single-electron approximation <cit.>. In our work, we consider the effects of electron-electron interaction for superstructures with θ>θ_c and not too small size of the supercell, when the single-electron gap can be neglected. We also assume that the transverse electric field is applied to the bilayer. This lifts the double degeneracy of four Dirac cones of the system; two of them are shifted upwards in energy, and the other two – downwards. As a result, two hole-like and two electron-like Fermi surface sheets appear in the system with a perfect nesting between hole-like and electron-like components. This leads to Fermi surface instability with respect to exciton band gap formation due to electron-electron interaction. We show also that the exciton order parameter is accompanied by the spin-density-wave (SDW) order. The dependence of the gap on the twist angle and on the bias voltage is analyzed. It is shown that the gap can be effectively controlled by external field, which may be useful for applications.

Geometry of twisted bilayer graphene—
The geometry of the twisted bilayer graphene is described in details, in several papers <cit.>. Here we only present the basic properties necessary for further considerations. Graphene monolayer has a hexagonal crystal structure consisting of two triangular sublattices A and B. Carbon atoms in the layer 1 are located in positions 𝐫_𝐧^1A=𝐫_𝐧^1≡ n𝐚_1+m𝐚_2 and 𝐫_𝐧^1B=𝐫_𝐧^1+δ with 𝐧=(n, m) (n, m are integers), where 𝐚_1,2=a(√(3),∓1)/2 are the graphene lattice vectors (a=2.46Å) and δ=a(1/√(3),0). The positions of atoms in layer 2 are  𝐫_𝐧^2B=𝐫_𝐧^2≡ d𝐞_z+n𝐚_1'+m𝐚_2' and 𝐫_𝐧^2A=𝐫_𝐧^2-δ', where 𝐚_1,2', δ'  are vectors  𝐚_1,2, δ rotated by angle θ, symbol 𝐞_z is the unit vector along the z-axis, and d=3.35Å  is the interlayer distance. We also introduce vectors δ^iα≡𝐫_𝐧^iα-𝐫_𝐧^i, which  are independent on 𝐧. Limiting case θ=0 corresponds to the AB stacking. The superstructure exists for twist angles equal

    cosθ=3m_0^2+3m_0r+r^2/2/3m_0^2+3m_0r+r^2 ,

where m_0 and r are co-prime positive integers. Superlattice vectors 𝐑_1,2 are linear combinations of 𝐚_1 and 𝐚_2 with integer-valued coefficients <cit.>. The magnitude of these vectors is <cit.> |𝐑_1,2|=rL/√(g), where L=a/[2sin(θ/2)] is the moiré period and g=1 if r≠3n, or g=3 if r=3n (n is integer). Thus, only superstructures with r=1 coincide with the moiré lattice. The number of graphene unit cells of each layer inside the supercell is N_sc=(3m_0^2+3m_0r+r^2)/g. Thus, the number of carbon atoms in the superlattice cell is equal to 4N_sc.

We introduce 𝐛_1,2=2π(1/√(3),∓1)/a, which are the reciprocal lattice vectors of the layer 1, and 𝐛_1,2' for layer 2 (𝐛_1,2' are connected to 𝐛_1,2 by rotation on angle θ). Vectors G_1,2 are the reciprocal vectors for superlattice. All these vectors are related to each other according to the following formulas: 𝐛_1'=𝐛_1-r G_1 and 𝐛_2'=𝐛_2+r( G_1+ G_2) if r≠3n, or 𝐛_1'=𝐛_1+r( G_1+2 G_2)/3 and 𝐛_2'=𝐛_2-r(2 G_1+ G_2)/3, otherwise. Each graphene layer has two non-equivalent Dirac points located at corners of its Brillouin zone. Thus, the total number of Dirac points for the bilayer is four. The Brillouin zone of the superlattice has a shape of hexagon. It can be obtained by N_sc times folding of the Brillouin zone of the layer 1 or 2. As a result of this folding, Dirac points of each layer are translated to two non-equivalent Dirac points of the superlattice, 𝐊_1 and  𝐊_2, located at corners of the reduced Brillouin zone. Thus, one can say that the Dirac points 𝐊_1,2 are doubly degenerate since each of them corresponds to two non-equivalent Dirac points of constituent layers. Points 𝐊_1,2 can be expressed via vectors G_1,2 as 𝐊_1=( G_1+2 G_2)/3 and 𝐊_2=(2 G_1+ G_2)/3.

Model Hamiltonian—
We start from the tight-binding model for the p_z electrons in twisted bilayer graphene. Electrons are assumed to interact via the screened Coulomb interaction. Formally, we write

    H   =   ∑_i𝐧j𝐦αβσt(𝐫_𝐧^iα;𝐫_𝐦^jβ)
    d^†_𝐧iασd^†_𝐦jβσ+
    V_b/2∑_𝐧(n_𝐧1-n_𝐧2)+
          1/2∑_i𝐧j𝐦αβσσ'd^†_𝐧iασd^†_𝐧iασ
    V(𝐫_𝐧^iα-𝐫_𝐦^jβ)d^†_𝐦jβσ'd^†_𝐦jβσ' ,

where d^†_𝐧iασ and d^†_𝐧iασ are the creation and annihilation operators of the electron with spin projection σ, located at site 𝐧 in the layer i(=1,2) in the sublattice α(=A,B), and n_𝐧i=∑_ασd^†_𝐧iασd^†_𝐧iασ. The first term describes the intersite hopping. For intralayer hopping we consider only the nearest-neighbor term with amplitude -t, where t=2.57 eV. The interlayer hopping amplitudes are parameterized as described in Refs. <cit.>, with the largest interlayer hopping amplitude equal to t_0=0.4 eV. Second term describes the potential energy difference between layers due to the applied bias voltage V_b. Third term corresponds to the Coulomb interaction between electrons. The precise form of the function V(𝐫) will be discussed below.

Let us first consider single-particle part of the Hamiltonian (<ref>) (first and second terms). We proceed to the momentum representation introducing electronic operators

    d^†_𝐩𝐆iασ=1/√( N)∑_𝐧
    e^-i(𝐩+𝐆)𝐫_𝐧^id_𝐧iασ ,

where N is the number of graphene unit cells in the sample in one layer, momentum 𝐩 lies in the first Brillouin zone of the superlattice, while 𝐆=n G_1+m G_2 is the reciprocal vector of the superlattice, lying in the first Brillouin zone of the ith layer. The number of such vectors 𝐆 is equal to N_sc for each graphene layer. In this representation the single-particle part of Hamiltonian (<ref>) can be written as

    H_0   =   ∑_𝐩[∑_𝐆_1𝐆_2 ijαβσt̃_ij^αβ(𝐩+𝐆_1;𝐆_1-𝐆_2)
    d^†_𝐩𝐆_1iασd^†_𝐩𝐆_2jβσ.+
          .V_b/2∑_𝐆ασ(
    d^†_𝐩𝐆1ασd^†_𝐩𝐆1ασ-
    d^†_𝐩𝐆2ασd^†_𝐩𝐆2ασ)],

where

    t̃_ij^αβ(𝐤;𝐆)=1/N_sc∑'_𝐧𝐦e^-i𝐤(𝐫_𝐧^i-𝐫_𝐦^j)
    e^-i𝐆𝐫_𝐦^jt(𝐫_𝐧^iα;𝐫_𝐦^jβ) .

In the last formula, the summation over 𝐦 is performed over sites inside the zeroth supercell, while summation over 𝐧 is performed over all sites in the sample (or vice versa).

For a given momentum 𝐩, Hamiltonian (<ref>) can be represented in the form of 4N_sc×4N_sc matrix. We diagonalize this matrix numerically, and calculate both the spectrum E_𝐩^(S) (S=1, 2, …, 4N_sc) and eigenvectors U^(S)_𝐩𝐆iα. Figure <ref> shows the band structure calculated for the sample with m_0=5, r=1 (θ≅6.009^∘), and for bias voltage V_b=0.15t. Only bands within the energy window |E|/t<0.5 are shown. Low-energy spectrum consists of four Dirac cones located in pairs near two Dirac points of the superlattice, 𝐊_1 and 𝐊_2. Bias voltage shifts apexes of two of these cones to positive energies, and two – to negative energies. The electron density corresponding to upper (lower) Dirac cones is concentrated mostly in layer 1 (layer 2), even though the interlayer hybridization tends to distribute it uniformly between the layers. When bias voltage is applied, the system acquires a Fermi surface, which consists of two closed curves (valleys) located near two Dirac points. These curves are nearly circular when the bias voltage is small enough, while trigonal warping reveals itself at larger V_b (see Fig. <ref>). The trigonal warping becomes more pronounced for superstructures with smaller θ.



The most important feature of the Fermi surface shown in Fig. <ref> is its double degeneracy; in each valley, the Fermi surface curve is created by both the electron-like band corresponding to the lower Dirac cone, and the hole-like band belonging to the upper Dirac cone. In other words, we have a situation with perfect nesting of Fermi surface. This leads to the instability of a Fermi liquid state with respect to the formation of some type of ordering due to the electron-electron interaction.

Consider now the interaction term of Hamiltonian (<ref>). In terms of electron operators d^†_𝐩𝐆iασ, Eq. (<ref>), the last term in Eq. (<ref>) takes the form

    1/2 N∑_iασ jβσ'∑_𝐩_1𝐩_2𝐩_1'𝐩_2'∑_𝐆_1𝐆_2𝐆_1'𝐆_2'∑_𝐆δ_𝐆, 𝐤_1+𝐤_2-𝐤_1'-𝐤_2'×
          d^†_𝐩_1𝐆_1iασd^†_𝐩_1'𝐆_1'iασV_iα;jβ(𝐤_1-𝐤_1';𝐆)
    d^†_𝐩_2𝐆_2jβσ'd^†_𝐩_2'𝐆_2'jβσ',

where 𝐤_1,2=𝐩_1,2+𝐆_1,2, 𝐤_1,2'=𝐩_1,2'+𝐆_1,2', and

    V_iα;jβ(𝐤;𝐆)=1/N_sc∑'_𝐧𝐦
    e^-i𝐤(𝐫_𝐧^i-𝐫_𝐦^j)
    e^-i𝐆𝐫_𝐦^jV(𝐫_𝐧^iα-𝐫_𝐦^jβ) .

As before, in the last equation the summation over 𝐦 is performed over sites inside the zeroth supercell, while summation over 𝐧 is performed over all sites in the sample. For intralayer (i=j) interaction one can separate the summation on 𝐧 and 𝐦 by substitution 𝐧→𝐧+𝐦. As a result, we obtain V_iα;iβ(𝐤;𝐆)=(∑_𝐛^iδ_𝐛^i, 𝐆)
∑_𝐧e^-i𝐤𝐫_𝐧^iV(𝐫_𝐧^i+δ^iα-δ^iβ), where the first and second summations are performed over all reciprocal lattice vectors (𝐛^i) and all lattice sites of the layer i, correspondingly. Let us denote by 𝐤̃^i the `vector 𝐤 modulo 𝐛^i', that is, the vector lying in the first Brillouin zone of the layer i and coinciding with 𝐤 upon the translation on some reciprocal vector 𝐛^i. By definition (<ref>), we have V_iα;iβ(𝐤;𝐆)=V_iα;iβ(𝐤̃^i;𝐆). Below we will use the continuum (low-𝐤) approximation for V_iα;iβ(𝐤;𝐆), when one can substitute the summation over lattice sites by the 2D integration. As a result, we obtain

    V_iα;iβ(𝐤;𝐆)=1/ V_cV_ii(𝐤̃^i)
    e^i𝐤̃^i(δ^iα-δ^iβ)∑_𝐛^iδ_𝐛^i, 𝐆 ,

where V_c=√(3)a^2/2 is the graphene's unit cell area, and V_ii(𝐤)=∫ d^2rV(𝐫)e^-i𝐤𝐫 is the Fourier transform of the function V(𝐫).

We introduce also the Fourier transform for interlayer interaction  as V_ij(𝐤)=∫ d^2rV(d𝐞_z+𝐫)e^-i𝐤𝐫 (i≠ j). Substituting this equation to Eq. (<ref>), one obtains

    V_iα;jβ(𝐤;𝐆)   =   1/ V_c∑_𝐛^iV_ij(𝐤+𝐛^i)
    e^i(𝐤+𝐛^i)(δ^iα-δ^iβ)×
          ∑_𝐛^jδ_𝐛^j, 𝐛^i+𝐆 .

For functions V_ij(𝐤) we use expressions for the screened Coulomb potential in the form <cit.>:

    V_ij(𝐪)=v_𝐪/1+Π_𝐪v_𝐪(
    [     1 e^-qd; e^-qd     1 ]) ,

where bare Coulomb potential is v_𝐪=2π e^2/ϵ q. The permittivity of the substrate is ϵ and Π_𝐪≡-P(ω=0,𝐪), where P(ω,𝐪) is the polarization operator of the bilayer. Since interlayer potential V_ij(𝐪) decays exponentially with q and e^-|𝐛_1,2|d≈5×10^-5, one can take only one term in the sum over 𝐛^i in Eq. (<ref>), such that  𝐤+𝐛^i=𝐤̃^i. Below we will use the long wave length approximation for the function Π_𝐪. In this case, it is independent on 𝐪 and equals to the density of states of the bilayer at the Fermi level. We calculate the latter quantity following procedure described below.

Exciton order parameter—
To go further, we introduce new electronic operators ψ_𝐩Sσ defined according to the relations:

    ψ_𝐩Sσ=∑_𝐆iαU^(S)*_𝐩𝐆iαd_𝐩𝐆iασ ,  
    d_𝐩𝐆iασ=∑_SU^(S)_𝐩𝐆iαψ_𝐩Sσ .

In this representation,  the single-particle part of the Hamiltonian (<ref>) is diagonal. The interaction term takes a form ∑ V(𝐩_1S_1…𝐩_4S_4)
ψ^†_𝐩_1S_1σψ^†_𝐩_2S_2σψ^†_𝐩_3S_3σ'ψ^†_𝐩_4S_4σ', where the summation is performed over all indices in this expression, and V(𝐩_1S_1;𝐩_2S_2;𝐩_3S_3;𝐩_4S_4) is the convolution of the product of the function V_iα;jβ(𝐤;𝐆) and four wave functions U^(S)_𝐩𝐆iα.

Let us denote the indices of two electron-like and two hole-like bands closest to the Fermi level by S^μ_+ and S^μ_-, correspondingly, where μ=1,2. The perfect nesting of hole-like and electron-like Fermi surfaces gives rise to a dielectric order parameter driven by the Coulomb interaction. In general, this order parameter is a superposition of the expectation values of the type ⟨ψ^†_𝐩Sσψ^†_𝐩S'σ'⟩, where S≠ S'. Here we assume that non-zero expectation values are only those, which couple the bands S^μ_+ and S^μ_- (see Fig. <ref>). Analysis shows also that corresponding order is of the SDW type. Assuming planar spin configuration, we obtain that non-zero expectation values are the following

    η_𝐩μσ=⟨ψ^†_𝐩S^μ_+σψ^†_𝐩S^μ_-σ̅⟩=
    ⟨ψ^†_𝐩S^μ_-σ̅ψ^†_𝐩S^μ_+σ⟩^*,

where σ̅ means `not σ'. Note that possible charge-density-wave order, corresponding to the non-zero η̃_𝐩μ=∑_σ⟨ψ^†_𝐩S^μ_+σψ^†_𝐩S^μ_-σ⟩ is energetically unfavorable in comparison to the SDW ones.

We will study the model (<ref>) in the mean-field approximation. This scheme involves the replacement of the product of two operators O_1O_2→ O_1⟨ O_2⟩+⟨ O_1⟩ O_2-⟨ O_1⟩⟨ O_2⟩, where ⟨ O_1,2⟩∝η_𝐩μσ. As a result, the interaction term of the Hamiltonian becomes quadratic in ψ^†_𝐩Sσ operators. In addition, we truncate the total mean-field Hamiltonian keeping only bands S^μ_±. Moreover, in the interaction part of this Hamiltonian we keep only terms coupling the electron band S^μ_+ with the hole band S^μ_- with the same μ. As a result, the effective mean-field Hamiltonian becomes

    H_MF=∑_𝐩μΨ^†_𝐩μĤ_𝐩μΨ^†_𝐩μ+U_c ,

where we introduce the following 4-component operator Ψ^†_𝐩μ=(ψ^†_𝐩S^μ_+↑,ψ^†_𝐩S^μ_-↑,
ψ^†_𝐩S^μ_+↓,ψ^†_𝐩S^μ_-↓)^T.
In equation (<ref>), Ĥ_𝐩μ is the 4×4 matrix

    Ĥ_𝐩μ=([ E^(S^μ_+)_𝐩           0           0    -Δ_𝐩μ↑^*;           0 E^(S^μ_-)_𝐩      -Δ_𝐩μ↓           0;           0    -Δ_𝐩μ↓^* E^(S^μ_+)_𝐩           0;      -Δ_𝐩μ↑           0           0 E^(S^μ_-)_𝐩 ]),

where Δ_𝐩μσ is the order parameter having the form

    Δ_𝐩μσ=1/ N∑_𝐪ν[A_μν(𝐩;𝐪)η_𝐪νσ+
    B_μν(𝐩;𝐪)η^*_𝐪νσ̅],

where  A_μν(𝐩;𝐪)≡ V(𝐩S^μ_-;𝐪S^ν_-;𝐪S^ν_+;𝐩S^μ_+) and B_μν(𝐩;𝐪)≡ V(𝐩S^μ_-;𝐪S^ν_+;𝐪S^ν_-;𝐩S^μ_+). The term  U_c in Eq. (<ref>) is the c-number

    U_c   =   1/2 N∑_𝐩𝐪μνσ[A_μν(𝐩;𝐪)η^*_𝐩μση^*_𝐪νσ+.
          .B_μν(𝐩;𝐪)η^*_𝐩μση^*_𝐪νσ̅+H.c.].

The precise form of the functions A_μν(𝐩;𝐪) and B_μν(𝐩;𝐪) is the following:

    A_μν(𝐩;𝐪)   =   ∑_iα jβ∑_𝐆_1𝐆_2𝐆_1'𝐆_2'
    U^(S^μ_-)*_𝐩𝐆_1iαU^(S^ν_-)_𝐪𝐆_1'iαU^(S^ν_+)*_𝐪𝐆_2jβU^(S^μ_+)_𝐩𝐆_2'jβ×
          V_iα;jβ(𝐩-𝐪+𝐆_1-𝐆_1';𝐆_1+𝐆_2-𝐆_1'-𝐆_2'),
    
    B_μν(𝐩;𝐪)   =   ∑_iα jβ∑_𝐆_1𝐆_2𝐆_1'𝐆_2'
    U^(S^μ_-)*_𝐩𝐆_1iαU^(S^ν_+)_𝐪𝐆_1'iαU^(S^ν_-)*_𝐪𝐆_2jβU^(S^μ_+)_𝐩𝐆_2'jβ×
          V_iα;jβ(𝐩-𝐪+𝐆_1-𝐆_1';𝐆_1+𝐆_2-𝐆_1'-𝐆_2').


Minimizing the total energy at zero temperature and at half-filling, we obtain the system of equations for the order parameters:

    Δ_𝐩μσ=1/2 N∑_𝐪ν[A_μν(𝐩;𝐪)Δ_𝐪νσ/√(|Δ_𝐪νσ|^2+E^2_𝐪ν)+
    B_μν(𝐩;𝐪)Δ^*_𝐪νσ̅/√(|Δ_𝐪νσ̅|^2+E^2_𝐪ν)],

where E_𝐪μ=[E^(S^μ_+)_𝐪-E^(S^μ_-)_𝐪]/2.

For a given superstructure and bias voltage we calculate the functions A_μν(𝐩;𝐪) and B_μν(𝐩;𝐪) numerically. Analysis shows, that with a good accuracy the following relations take place:

    A_μν(𝐩;𝐪)   ≈   δ_μνA(𝐩;𝐪) ,
    
    B_μν(𝐩;𝐪)   ≈   δ_μν̅B(𝐩;𝐪) ,

where ν̅ means `not ν'. The deviations from these equalities do not exceed 1% for any superstructures considered. Note that  in the limit of uncoupled graphene layers, t_0=0, Eqs. (<ref>) become exact. This follows from the symmetry of the wavefunctions presented in the definition of the functions  A_μν and B_μν, Eqs. (<ref>): for uncoupled layers the electrons are localized either in layer 1 or 2.  Further, by choosing appropriate phases of the wave functions U^(S^μ_±)_𝐩𝐆iα one can make A(𝐩;𝐪) and B(𝐩;𝐪) real-valued functions. In this case, one can choose Δ_𝐩μσ to be real-valued functions satisfying the relations Δ_𝐩μ↑=Δ_𝐩μ↓≡Δ_𝐩μ. The equation for the order parameters then becomes

    Δ_𝐩μ=1/2∫d^2q/v_BZ[A(𝐩;𝐪)Δ_𝐪μ/√(Δ_𝐪μ^2+E^2_𝐪μ)+
    B(𝐩;𝐪)Δ_𝐪μ̅/√(Δ_𝐪μ̅^2+E^2_𝐪μ̅)],

where v_BZ=8π^2/(a^2√(3)) is the area of the graphene's Brillouin zone, while the integration is performed over the Brillouin zone of the superlattice. In numerical calculations we take A(𝐩;𝐪)=[A_11(𝐩;𝐪)+A_22(𝐩;𝐪)]/2 and B(𝐩;𝐪)=[B_12(𝐩;𝐪)+B_21(𝐩;𝐪)]/2.

Approximate solution, weak interaction limit—
The main interest is the value of the function Δ_𝐩1 at the Fermi surface since it gives us the energy gap.  Numerical analysis shows that if momenta 𝐩 and 𝐪 belong to different valleys, that is, they are located near different Dirac points, we have A(𝐩;𝐪)≈0 and B(𝐩;𝐪)≈0. This property makes it possible to consider the functions Δ_𝐩μ near each Dirac point, 𝐊_1 and 𝐊_2, independently. Let's take for example the 𝐊_1 point. The Fermi surface near each Dirac point is a closed curve having near-circular shape. Below we neglect the trigonal warping and approximate the Fermi surface by a circle with radius q_F^* calculated numerically by averaging over the Fermi surface. We assume also that Δ_𝐩μ are the step-like functions, describing by the equations Δ_𝐊_1+𝐩μ=Δ_μΘ(q_Λ-|p-q_F^*|), where the cutoff momentum q_Λ will be specified below. Thus, the region of integration in Eq. (<ref>) becomes a ring centered at the Dirac point 𝐊_1 and having radii q_F^*-q_Λ and q_F^*+q_Λ if q_Λ<q_F^*, or the circle with radius q_Λ, otherwise. In further approximation, we replace the functions A(𝐩;𝐪) and B(𝐩;𝐪) by constants A̅ and B̅ obtained by averaging of  A(𝐩;𝐪) and B(𝐩;𝐪) over the Fermi surface. We also approximate energies E_𝐪μ by linear functions

    E_𝐊_1+𝐪1≈ v_F^*(q-q_F^*) ,  E_𝐊_1+𝐪2≈ v_F^*(q+q_F^*) ,

where the renormalized Fermi velocity is calculated numerically by averaging the function ∑_μ|∂ E_𝐊_1+𝐪μ∂𝐪|/2 over the Fermi surface. As a result, the system of equations for order parameters becomes

    Δ_1   =   ∫_q_1^q_2dq[1/2qλ_AΔ_1/√(Δ_1^2+v_F^*2(q-q_F^*)^2)+1/2qλ_BΔ_2/√(Δ_2^2+v_F^*2(q+q_F^*)^2)],
    Δ_2   =   ∫_q_1^q_2dq[1/2qλ_AΔ_2/√(Δ_2^2+v_F^*2(q+q_F^*)^2)+1/2qλ_BΔ_1/√(Δ_1^2+v_F^*2(q-q_F^*)^2)],

where q_1=q_F^*-q_Λ, q_2=q_F^*+q_Λ if q_Λ<q_F^*, or q_1=0, q_2=q_Λ, otherwise, λ_A=2πA̅/v_BZ, and λ_B=2πB̅/v_BZ. We solve this system numerically.

The magnitude of the order parameters Δ_μ depends on the values of A̅ and B̅, as well as on the cut-off momentum q_Λ. Our analysis shows that the main contribution to the functions A(𝐩;𝐪) and B(𝐩;𝐪) comes from the interlayer interaction. Following Refs. <cit.> we define q_Λ from the condition V_12(q_Λ)=V_12(0)/2.  Assuming that q_Λd≪1 (which is correct for V_b≲ t_0 and e^2/ϵ v_F≲1), from Eqs. (<ref>) we obtain the estimate q_Λ≈2πΠ_0v_Fα, where Π_0 is the density of states of the bilayer at the Fermi level, α=e^2/ϵ v_F is the graphene's fine structure constant, and v_F=ta√(3)/2 is the Fermi velocity of the single layer graphene. Neglecting the trigonal warping, the density of states is expressed as Π_0≈4q_F^*/(π v_F^*). Thus, we obtain for q_Λ

    q_Λ=8v_F/v_F^*α q_F^* .


Approximate solution, strong interaction limit—
The limit of weak coupling corresponds the case of α≪1. Simple estimates show, however, that for bilayer suspended in vacuum, ϵ=1, the parameter α≈2.6. When α increases, the cut-off momentum can exceed the size of the superlattice's Brillouin zone. In this case, we should increase the number of bands in our effective Hamiltonian. Simultaneously, we should increase, the number of order parameters Δ_𝐩μ, with μ now changing from 1 to some N>2. The rank of the matrix functions A_μν and B_μν becomes equal to N. This consideration can be substantially simplified if we consider  the added `high-energy' bands in the limit of decoupled (t_0=0) graphene layers. This approximation is justified, since for these bands we have |E_𝐩^(S)|≳ t_0. For decoupled layers, one can associate the band index S to the momentum 𝐩 lying in the Brillouin zone of the layer 1 or 2, that is, one can perform the band unfolding <cit.>. As a result, one can assume that the number of order parameters Δ_𝐩μ is still equal to 2, but the integration in Eqs. (<ref>) or (<ref>) is extended to the momenta exceeding the reciprocal  unit cell of the superlattice. Applying this procedure, one must keep in mind that the two valleys should be still considered independently. To understand why this is so, let us consider the situation in the unfolded Brillouin zone from the beginning. Layer 1 has two non-equivalent Dirac points, 𝐊 and 𝐊'(=-𝐊). Rotation by twist angle θ transforms them into Dirac points of the layer 2, 𝐊_θ and 𝐊'_θ. The considered ordering corresponds to the formation of the electron-hole pair consisting of the electron with momentum 𝐊+𝐩 and the hole with momentum 𝐊_θ'+𝐩 (valley 1), and the electron with momentum 𝐊'+𝐩 and the hole with momentum 𝐊_θ+𝐩 (valley 2). [For superstructures with r≠3n. Otherwise, the electron and the hole have momenta 𝐊+𝐩 and 𝐊_θ+𝐩 in valley 1, and 𝐊'+𝐩 and 𝐊_θ'+𝐩 in valley 2.] It is clear, that the annihilation of such a pair in the valley 1 with simultaneous creation of the pair in the valley 2 is prohibited by the momentum conservation law. Thus, in our model, the intervalley scattering can be neglected [This picture is valid, of course, for E_Λ≡ v_Fq_Λ≲ t. At larger energies, the spectrum cannot be described by separate Dirac cones.].

The procedure described above implies the knowledge of the energies E_𝐩μ at large momenta. We approximate E_𝐊_1+𝐩μ by Eqs. (<ref>), when p<q_0, where q_0=| G_1,2|/(2√(3)) is the radius of the circle centered at Dirac point 𝐊_1 and touching the edges of the reciprocal  unit cell of the superlattice (see Fig. <ref>).  At larger p we use the limit of decoupled layers

    E_𝐊_1+𝐩1≈ E_𝐊_1+𝐩2≈ v_Fp .

With this accuracy we neglect the effect of the bias voltage at high energies, since V_b≲ t_0. As a result, approximate equation for the order parameters Δ_μ become

    Δ_1   =   ∫_0^q_0dq1/2qλ_AΔ_1/√(Δ_1^2+v_F^*2(q-q_F^*)^2)+
    ∫_q_0^q_Λdq1/2qλ_AΔ_1/√(Δ_1^2+v_F^2q^2)+
          ∫_0^q_0dq1/2qλ_BΔ_2/√(Δ_2^2+v_F^*2(q+q_F^*)^2)+
    ∫_q_0^q_Λdq1/2qλ_BΔ_2/√(Δ_2^2+v_F^2q^2) ,
    Δ_2   =   ∫_0^q_0dq1/2qλ_AΔ_2/√(Δ_2^2+v_F^*2(q+q_F^*)^2)+
    ∫_q_0^q_Λdq1/2qλ_AΔ_2/√(Δ_2^2+v_F^2q^2)+
          ∫_0^q_0dq1/2qλ_BΔ_1/√(Δ_1^2+v_F^*2(q-q_F^*)^2)+
    ∫_q_0^q_Λdq1/2qλ_BΔ_1/√(Δ_1^2+v_F^2q^2) .




Results and Discussion—
We solve the system of equations (<ref>) [or (<ref>), if q_Λ<q_0] numerically for several superstructures with r=1 in a wide range of V_b and α. We found that the ratio Δ_2/Δ_1≳0.3 for any interaction strength considered and goes to 1 when α increases. Our main interest is the value of Δ_1 since it gives us the energy gap. The dependencies of the band gap Δ_1 on α for superstructures with m_0=2, 4, 6, 8, 10 (corresponding twist angles are θ≅17.7, 7.3, 5.1, 3.9, 3.2 degrees) are presented in Fig. <ref>(a). The gap strongly (exponentially) depends on the interaction strength for all superstructures. It is seen from this figure that the gap is considerably large only when α≳1, that is, when ϵ≲2.5. The ratio Δ_1/t∼10^-2 corresponds to the value of the gap about 300 K. Thus, according to our calculations, in order to observe the gap at room temperatures, the permittivity of the sample should not be large.

The important result demonstrating in Fig. <ref>(a) is that for any α, the band gap is larger for superstructures with smaller twist angles (larger m_0). This is illustrated in Fig. <ref>(b), where we plot the dependence of Δ_1 on θ calculated at α=1.044 and V_b/t=0.037. We see that the band gap increases by about 4 orders of magnitude, when twist angle changes from θ≅17.7^∘ (m_0=2) down to θ≅3.2^∘ (m_0=10). Such a strong enhancement can be explained by the reduction of the Fermi velocity due to the interlayer hybridization. To illustrate this, let us consider the weak interaction limit. Assuming that q_λ≪ q_F^* and Δ_μ≪ V_b, one can solve the system of equations (<ref>) analytically. This gives

    Δ_1≈2v_Fq_F^*α e^-1/Λ+4α^*,

where α^*=e^2/(ϵ v_F^*) is the `renormalized α' and

    Λ≈2πA̅q_F^*/v_BZv_F^* .

Calculations show that good approximation for A̅ and, consequently, for Λ can be obtained if in Eq. (<ref>) for A_μν(𝐩;𝐪)
we will use wave functions corresponding to the limit of decoupled layers, t_0=0. In this case, expressions for U^(S)_𝐩𝐆iα can be found analytically, and as a result, we obtain

    A̅≈⟨cos^2(φ/2)V_12(2q_F^*|sin(φ/2)|)⟩ ,

where φ is the polar angle parameterizing points on the Fermi surface, while averaging is performed over this angle. The factor cos^2(φ/2) before potential V_12 is inherited from the wave functions. Substituting this expression into Eq. (<ref>), we obtain

    Λ≈1/8⟨cos^2(φ/2)/1+1/4α^*|sin(φ/2)|⟩.

Deriving the latter equality we assumed that q_F^*d≪1, which is valid for all bias voltages considered. According to Eq. (<ref>), Λ  increases when α^* increases. When α^*≪1, we have Λ∝α^*ln(1/α^*). Since α^* inversely proportional to v_F^*, parameter Λ and, consequently, Δ_1 increases when v_F^* decreases. This can explain the dependence of Δ_1 on θ since v_F^* decreases with the twist angle <cit.>. Numerical calculations show also, that actual value of A̅ is greater than estimate (<ref>). This can be explained by the fact that at finite interlayer hybridization, the quasiparticles are no longer localized in one particular layer. As a result, intralayer potential contributes to A̅, and this effect is stronger for smaller twist angles.

It is seen from Eqs. (<ref>), that exponent in Eq. (<ref>) is independent on the bias voltage. This result correlates with that obtained in Refs. <cit.>. Such a behavior can be understood if we realize that parameter Λ is a product of the interaction strength A̅∝1/q_F^* and the density of states at the Fermi level ρ_0∝ q_F^*∝ V_b. Thus, only pre-exponential factor depends the bias voltage. Numerical analysis shows, that at small V_b, the band gap linearly depends on V_b for any α, while at strong interaction the function Δ_1(V_b) can show non-monotonous behavior at larger bias voltages [see the inset to Fig. <ref>(a)].

The order parameter goes to zero when V_b→0 for any superstructures and for any strength of the interaction considered. At the same time, it is seen from Fig. <ref>(b), that Δ_1 increases very fast with the decrease of the twist angle. We did not analyze what happens at twist angles very close and below critical value (θ_c≅1.89^∘ in our model), but we can expect that for these θs the Coulomb interaction can stabilize some type of ordering even at zero bias voltage. This suggestion is confirmed by recent studies performed in Ref. <cit.>, where authors predict the SDW ground state for bilayers with θ<θ_c in the framework of the Hubbard model. Note that our `exciton-plus-SDW' order parameter is stabilized mainly by the interlayer interaction, and it would not occur (or would be strongly suppressed) for Hubbard interaction. Similar type of order has been considered in Ref. <cit.> devoted to the study of the AA-stacked bilayer graphene in the applied electric field. The detailed investigation of the ordering type for bilayers with the twist angles close to the critical value is the subject for future study.

In conclusion, we studied the ground state of the twisted bilayer graphene when the bias voltage is applied to the system. We showed that the bias voltage forms two hole-like and two electron-like Fermi surfaces with perfect nesting. As a result of such a band structure, the screened Coulomb interaction stabilizes the exciton order parameter in the system. The exciton order parameter is accompanied by the spin-density-wave order. The value of the gap depends on the twist angle and on the applied voltage. The latter property is quite useful for different applications in electronics.

Acknowledgments.— This work is partially supported by the
Russian Foundation for Basic Research (Projects 17-02-00323, and
15-02-02128), JSPS-RFBR grant No 17-52-50023,




the RIKEN iTHES Project,
MURI Center for Dynamic Magneto-Optics via the AFOSR Award
No. FA9550-14-1-0040,
the Japan Society for the Promotion of Science (KAKENHI),
the IMPACT program of JST,
JSPS-RFBR grant No 17-52-50023,
CREST grant No. JPMJCR1676,
and the Sir John Templeton Foundation.

apsrevlong_no_issn_url


19
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1

[Rozhkov et al.(2016)Rozhkov, Sboychakov,
  Rakhmanov, and Nori]ourBLGreview2016
authorA. Rozhkov,
  authorA. Sboychakov,
  authorA. Rakhmanov,
  and authorF. Nori,
  “titleElectronic properties of graphene-based bilayer
  systems,” journalPhysics Reports
  volume648, pages1 (year2016).

[Mele(2012)]MeleReview
authorE. J. Mele,
  “titleInterlayer coupling in rotationally faulted multilayer
  graphenes,” journalJournal of Physics D: Applied Physics
  volume45, pages154004
  (year2012).

[Brihuega et al.(2012)Brihuega, Mallet,
  González-Herrero, Trambly de Laissardière, Ugeda, Magaud,
  Gómez-Rodríguez, Ynduráin, and Veuillen]STM1
authorI. Brihuega,
  authorP. Mallet,
  authorH. González-Herrero,
  authorG. Trambly de Laissardière,
  authorM. M. Ugeda,
  authorL. Magaud,
  authorJ. M. Gómez-Rodríguez,
  authorF. Ynduráin,
  and authorJ.-Y.
  Veuillen, “titleUnraveling the Intrinsic and
  Robust Nature of van Hove Singularities in Twisted Bilayer Graphene by
  Scanning Tunneling Microscopy and Theoretical Analysis,”
  journalPhys. Rev. Lett. volume109,
  pages196802 (year2012).

[Luican et al.(2011)Luican, Li, Reina,
  Kong, Nair, Novoselov, Geim, and Andrei]STM2
authorA. Luican,
  authorG. Li,
  authorA. Reina,
  authorJ. Kong,
  authorR. R. Nair,
  authorK. S. Novoselov,
  authorA. K. Geim, and
  authorE. Y. Andrei,
  “titleSingle-Layer Behavior and Its Breakdown in Twisted
  Graphene Layers,” journalPhys. Rev. Lett.
  volume106, pages126802
  (year2011).

[Ohta et al.(2012)Ohta, Beechem,
  Robinson, and Kellogg]SynthesysMoire
authorT. Ohta,
  authorT. E. Beechem,
  authorJ. T. Robinson,
  and authorG. L.
  Kellogg, “titleLong-range atomic ordering and
  variable interlayer interactions in two overlapping graphene lattices with
  stacking misorientations,” journalPhys. Rev. B
  volume85, pages075415
  (year2012).

[Lopes dos Santos et al.(2012)Lopes dos
  Santos, Peres, and Castro Neto]dSPRB
authorJ. M. B. Lopes dos Santos,
  authorN. M. R. Peres,
  and authorA. H.
  Castro Neto, “titleContinuum model of the
  twisted graphene bilayer,” journalPhys. Rev. B
  volume86, pages155449
  (year2012).

[Mele(2010)]MelePRB1
authorE. J. Mele,
  “titleCommensuration and interlayer coherence in twisted
  bilayer graphene,” journalPhys. Rev. B
  volume81, pages161405
  (year2010).

[Shallcross et al.(2010)Shallcross,
  Sharma, Kandelaki, and Pankratov]Pankratov1
authorS. Shallcross,
  authorS. Sharma,
  authorE. Kandelaki,
  and authorO. A.
  Pankratov, “titleElectronic structure of
  turbostratic graphene,” journalPhys. Rev. B
  volume81, pages165105
  (year2010).

[Sboychakov et al.(2015)Sboychakov,
  Rakhmanov, Rozhkov, and Nori]ourTBLG
authorA. O. Sboychakov,
  authorA. L. Rakhmanov,
  authorA. V. Rozhkov,
  and authorF. Nori,
  “titleElectronic spectrum of twisted bilayer graphene,”
  journalPhys. Rev. B volume92,
  pages075402 (year2015).

[Rozhkov et al.(2017)Rozhkov, Sboychakov,
  Rakhmanov, and Nori]ourTBLG2017
authorA. Rozhkov,
  authorA. Sboychakov,
  authorA. Rakhmanov,
  and authorF. Nori,
  “titleSingle-electron gap in the spectrum of twisted bilayer
  graphene,” journalPhysical Review B
  volume95, pages045119
  (year2017).

[Lopes dos Santos et al.(2007)Lopes dos
  Santos, Peres, and Castro Neto]dSPRL
authorJ. M. B. Lopes dos Santos,
  authorN. M. R. Peres,
  and authorA. H.
  Castro Neto, “titleGraphene Bilayer with a
  Twist: Electronic Structure,” journalPhys. Rev. Lett.
  volume99, pages256802
  (year2007).

[Bistritzer and MacDonald(2011)]PNAS
authorR. Bistritzer and
  authorA. H. MacDonald,
  “titleMoiré bands in twisted double-layer graphene,”
  journalProceedings of the National Academy of Sciences
  volume108, pages12233
  (year2011).

[Suárez Morell
  et al.(2010)Suárez Morell, Correa, Vargas, Pacheco, and
  Barticevic]Morell1
authorE. Suárez Morell,
  authorJ. D. Correa,
  authorP. Vargas,
  authorM. Pacheco, and
  authorZ. Barticevic,
  “titleFlat bands in slightly twisted bilayer graphene:
  Tight-binding calculations,” journalPhys. Rev. B
  volume82, pages121407
  (year2010).

[Trambly de Laissardière
  et al.(2010)Trambly de Laissardière, Mayou, and
  Magaud]NanoLettTB
authorG. Trambly de Laissardière,
  authorD. Mayou, and
  authorL. Magaud,
  “titleLocalization of Dirac Electrons in Rotated Graphene
  Bilayers,” journalNano Letters volume10,
  pages804 (year2010).

[Lozovik and Sokolik(2008)]LozovikSokolik
authorY. E. Lozovik and
  authorA. Sokolik,
  “titleElectron-hole pair condensation in a graphene
  bilayer,” journalJETP letters volume87,
  pages55 (year2008).

[Lozovik and Sokolik(2009)]LozovikPLA2009
authorY. E. Lozovik and
  authorA. Sokolik,
  “titleMulti-band pairing of ultrarelativistic electrons and
  holes in graphene bilayer,” journalPhysics Letters A
  volume374, pages326 (year2009).

[Nishi et al.(2017)Nishi, Matsushita, and
  Oshiyama]BandUnfoldingPRB2017
authorH. Nishi,
  authorY.-i. Matsushita,
  and authorA. Oshiyama,
  “titleBand-unfolding approach to moiré-induced band-gap
  opening and Fermi level velocity reduction in twisted bilayer graphene,”
  journalPhys. Rev. B volume95,
  pages085420 (year2017).

[Gonzalez-Arraga
  et al.(2017)Gonzalez-Arraga, Lado, Guine, and
  San-Jose]GuineaArXiv2017
authorL. A. Gonzalez-Arraga,
  authorJ. Lado,
  authorF. Guine, and
  authorP. San-Jose,
  “titleElectrically controllable magnetism in twisted bilayer
  graphene,” journalarXiv preprint arXiv:1702.08831
  (year2017).

[Akzyanov et al.(2014)Akzyanov,
  Sboychakov, Rozhkov, Rakhmanov, and Nori]AkzyanovAABLG2014
authorR. S. Akzyanov,
  authorA. O. Sboychakov,
  authorA. V. Rozhkov,
  authorA. L. Rakhmanov,
  and authorF. Nori,
  “titleAA-stacked bilayer graphene in an applied electric
  field: Tunable antiferromagnetism and coexisting exciton order
  parameter,” journalPhys. Rev. B
  volume90, pages155415
  (year2014).




