eduzhai > Physical Sciences > Materials Sciences >

Dynamic band structure tuning of graphene moire superlattice under pressure

  • sky
  • (0) Download
  • 20211014
  • Save S1. First principle calculations and parametrizations The atomic and electronic structures of graphene and hexagonal boron nitride (G/BN) interfaces are captured using input from ab-initio density functional theory (DFT) calculations. We rely on exact exchange and random phase approximation (EXX+RPA) input for atomic structure, while our treatment of the electronic structure is based on local-density approximation (LDA). Calculation of the dielectric constant from DFT One of the striking features of the experiment is an increase of the dielectric constant as a function of pressure. Our DFT calculations for bulk confirm this behavior showing increases on the order of 3% for the applied pressures up to 2.5 GPa as shown in Fig. S1. When the pressure is modified from P = 0 to P = 2.5 GPa, the EXX+RPA equilibrium interlayer distance is squeezed by ∼ 0.2 A˚ as we will discuss in the following subsection. Figure S1 Real part of the dielectric function as a function of energy ω obtained for different interlayer spacing in bulk AA stacked hexagonal boron nitride. The AA stacking have alternating B and N atoms in adjacent vertical layers. First principles calculations based on density functional theory (DFT) were performed using Quantum Espresso [1] under LDA within the plane wave basis set. Norm-conserving Vanderbilt WWW.NATURE.COM/NATURE | 1 pseudopotentials were used. The structures are fully optimized without any symmetry constraint by using conjugate gradient method. The convergence criterion for the force on each ion is taken to be less than 0.005 eV/A˚ while the energy is converged with a tolerance of 10−5 eV. In an insulating BN, the optical properties are dominated with the direct inter band contributions to the absorptive or imaginary part of dielectric function [2]. The matrix elements for a given interband transition β → α for a set of plane wave Bloch function: ∑ ψk,n = eiG·r uk,n = √1 V an,k,Gei(k+G)·r G (1) is given by: ∑ ∑ Mˆ α,β = a∗n,k,Gan ,k,GGα a∗n,k,Gan ,k,GGβ (2) G G These matrix elements account for the interband transitions with zero momentum transfer in the electric-dipole approximation, while the imaginary part of the dielectric tensor ε2α,β is a response function that can be obtained from a perturbation theory within the adiabatic approximation. All possible transitions from the occupied to the unoccupied states without local field effects are included in the imaginary dielectric function: ∑ ∑ ε2α,β = 4π e2 ΩNkm2 n,n k Mˆ α,β (Ek,n − Ek,n)2 f (Ek,n) + f (Ek,n) , Ek,n − Ek,n + h¯ ω + ih¯ Γ Ek,n − Ek,n − h¯ ω − ih¯ Γ (3) where Γ is the adiabatic parameter which must be zero due to conservation of energy. So, the equation can be rewritten in terms of Dirac delta functions: ∑ ∑ ε2α,β = 4π e2 ΩNkm2 n,n k Mˆ α,β f (Ek,n) (Ek,n − Ek,n)2 δ (Ek,n − Ek,n + h¯ ω) + δ (Ek,n − Ek,n − h¯ ω) , (4) However, interaction with the electromagnetic field which also take place even in the absence of photons, i.e spontaneous emission, brings intrinsic broadening that leads to a finite life time (Γ > 0). In the limit of non-vanishing Γ, the dielectric tensor takes the Drude-Lorentz form [2]: WWW.NATURE.COM/NATURE | 2 ∑ ∑ ∑ ε2α,β = 4π e2 ΩNkm2 n,k d f (Ek,n) dEk,n ηωMˆ α,β ω4 − η2ω2 + 8π e2 ΩNkm2 n=n k Mˆ α,β Ek,n − Ek,n [(ωk,n − Γω f (Ek,n) ωk,n)2 − ω2]2 + Γ2ω 2 , (5) and the real part of the dielectric function is derived from the Kramers-Kronig transformation: 2 ε1α,β = 1 + π ∞ ω ε2α,β (ω ω 2 −ω2 ) dω . (6) 0 The calculations of the dielectric function were carried out using the EPSILON package distributed with Quantum Espresso [1]. Our results were obtained based on local density approximation (LDA) electronic structure show that the enhancement of interlayer coupling due to pressure implies an overall redshift in the absorption that results in an increase in the dielectric constant. In the main text we have shown a qualitative agreement between theory and experiments pointing towards increase of dielectric constant with pressure but there is a quantitative disagreement of about a factor three. A better description of the many-body features missing in LDA calculations can be addressed [3] within the GW theory of Hedin for more accurate band gaps, often leading to two-fold enhancements, or the resolution of the Bethe-Salpeter equation for the excitonic effects. Other aspects such as the finite thickness of multilayer hBN used as the barrier material, the dielectric environment of the substrate, and details of the atomic structure as a function of pressure, the presence of multiple stacking arrangements, are other factors that should also be accounted for in a simulation that aims for a closer quantitative agreement with experiments. WWW.NATURE.COM/NATURE | 3 DFT deformation with pressure Figure S2 Left panel: Corrugation profiles as a function of pressure neglecting in-plane strains where the interaction with both the substrate and the top layers are accounted for. Right panel: The equilibrium interlayer distances are plotted against pressure for different stacking configurations in G/BN and for the most stable configuration in bulk hexagonal boron nitride (AA ) consistent with the EXX+RPA interlayer potentials of Ref. [4]. The mechanical behavior of nearly aligned graphene and boron nitride heterostructure (G/BN) is influenced by the encapsulating top BN layer. The interlayer interaction of nearly aligned G/BN is described using the assumption that they can be viewed as a collection of different stacking registries to be treated independently, from which the atomic structure of the whole sample can be obtained. We obtain the equilibrium spacing for each stacking point zs (s ∈ {AA, AB, BA}, see Fig. S2) by minimizing the total potential energy as a function of the separation distance between the top and bottom BN layer z. We also assume that the elastic resistance of the graphene sheet is negligible in the out-of-plane direction. The interlayer potential energy of the graphene sheet considering the effects of bottom and top BN layers is given by: ∑ 1 V (zAA, zAB, zBA, z) = 3 [Vs(zs) +V¯ (z − zs)] s (7) where V¯ (z − zs) = 1 3 ∑s Vs(z − zs) accounts for the averaging of the interlayer potential between the top BN and the graphene sheet due to their misalignment. We use for the interlayer inter- WWW.NATURE.COM/NATURE | 4 action potentials the exact exchange and random phase approximation (EXX+RPA) for G/BN and BN/BN interactions as implemented in VASP. [5] In equilibrium, the local pressure exerted by the encapsulating layer is equally balanced by that coming from the underlying substrate, i.e. Ps(zs) = P¯(z − zs). Defining the interlayer distance between graphene and the top layer zs = z − zs, we can write the equilibrium condition as 1 3 ∑s Ps(zs) = 1 3 ∑s P¯(zs) = P, where P is the external pressure applied on the system. Fig. S2 illustrates how the corrugation profiles of a graphene sheet change with pressure under the influence of two BN layers. Electronic coupling parameters and elastic properties of G/BN In our treatment of the moire´ bands in G/BN, we consider modifications in the low-energy regime of graphene electronic structure which are introduced by the interlayer tunneling. The modification is captured by two interlayer coupling terms V and V˜ whose values can be obtained from first-principle calculations. V and V˜ can be calculated from the following expressions [6]: √ V = 1 tN2C − tB2C , V˜ = 3 tN2C + tB2C , (8) 2 |εN| |εB| 2 |εN| |εB| where tNC (tBC) stands for the interlayer hopping term between C–N (C–B) atoms which highly depends on the interlayer separation, and εN (εB) are the on-site energies of N (B) atoms. These V and V˜ parameters assume values of 0.04 meV and 10 meV respectively when we use the interlayer coupling terms in Ref. [7] calculated for the interlayer separation of zr = 3.35 A˚ . In our semi-analytical model we adopt values of 0.01 meV and 18 meV respectively that account for the exchange enhancement of the primary Dirac cone gap [8]. The interlayer separation de- pendence of the electronic tunneling between the atomic orbitals of carbon and the underlying BN atoms is captured by three exponential factors {βAA, βAB, βBB} for each stacking configuration [8]. These parameters can be obtained by fitting ab-initio calculations data and yield values that range between 3.0–3.3 A˚ −1 [8]. As sample deformation is expected to play an important role in determining the band features, we model the structural relaxation using the Born-von Karman plate theory in which the elastic properties of graphene are fully characterized by two Lame´ constants: λg = 3.25 eV A˚ −2 WWW.NATURE.COM/NATURE | 5 and µg = 9.57 eV A˚ −2 [8]. One mechanism through which deformation leaves an imprint on the bands is the intralayer effects which are induced by the change of the carbon on-site energy and nearest-neighbor hoppings as each carbon atom is displaced with respect to its neighbors. They give rise to additional terms in the moire´ couplings which are equivalent to the electrostatic potentials and the pseudomagnetic fields, quantified by γV ≈ 4.0 eV and γB ≈ 4.5 eV respectively [9]. S2. Theoretical Model Graphene and BN are two-dimensional materials with hexagonal crystal structures, but with different lattice constants. When graphene is deposited on a BN substrate with near perfect alignment (small twist angles), the slight lattice mismatch (ε = aG−aBN aBN = −1.7%), leads to the formation of a moire´ structure which is characterized by a wavelength that is one or two orders of magnitude larger than graphene’s unit cell. As a consequence, the dynamics of the low energy Dirac electrons of graphene is governed by the corresponding long-wavelength component of the electronic couplings. Similarly, the interlayer potentials and deformations in G/BN are also largely characterized by this long-wavelength theory. It is by now established that the dominant contributions are captured by the first harmonics of the moire´ structure characterized by a set of vectors Gm, that is related with the original graphene’s lattice vectors gm: gm = Rˆ2π(m−1)/6 0, g , Gm = (1 + ε) − Rˆθ gm ≈ εgm − θ zˆ × gm, (9) where m ∈ {1, 2, ..., 6}, Rˆθ denotes a rotation by θ , g = 4π/3a is the length of graphene lattice vectors with a ≈ 1.42A˚ stands for the carbon-carbon distance, and the approximate sign indicates the approximation within small twist angle limit (θ 1). We also show in Fig. S3, that √ Gm defines the moire´ Brillouin Zone (mBZ), whose lateral dimension is scaled by ε˜ = ε2 + θ 2 with respect to graphene BZ that amounts to values 5% for twist angles θ 2◦. In G/BN system which possesses triangular symmetry, we can define two periodic functions within the first harmonics: (1) f1(r) = ∑m exp(iGm · r) which satisfies inversion and hexagonal symmetries, WWW.NATURE.COM/NATURE | 6 (a) (b) Figure S3 (a) Illustration of the moire´ Brillouin zone (mBZ) in relation to the reciprocal space of G/BN with slight misalignment θ . The presence of BN substrate introduces a characteristic long-wavelength potential which is one to two-order of magnitude larger than carbon-carbon distance, as defined by Gm. The Dirac cone is located on the black dot, at which BN sublattice asymmetry leads to a primary gap ∆p. The moire´ potentials lead to gapped secondary Dirac cones ∆s at the mBZ edges, marked by blue dots. (b) Deformation on graphene is characterized by the potential profiles on the special stacking points: AA (black), AB (blue) and BA (red). When pressure P is applied on the system, the carbon atoms increase their potential energies until a new equilibrium strain profile is reached. On the right panel we show for AB and BA stacking how the pressure can influence the in-plane and out-of-plane deformations, respectively proportional to variations in the potential (δU) and spacing (δ z). and (2) f2(r) = −i ∑m(−1)m exp(iGm · r) which is asymmetric under inversion. Definition of symmetric and antisymmetric deformations Deformation of the sample results from the minimization of the energy functional, leading to the following equation of motion for the in-plane components where we use notations similar to those in Ref. [10]. In the first Harmonic approximation we had used a scalar function defined in terms of three parameters U0, U1, U2 that captures the moire´ pattern in real space, WWW.NATURE.COM/NATURE | 7 Displacement vector AB BA AA B N Expansion Corrugation Figure S4 Plots of the symmetric (top) and asymmetric (bottom) components of the displacement vector, expansion and corrugation of graphene (defined in Eq. 11 and 12). To obtain the figures on the top row, we assume U1 = 3 meV, z1 = 0.01A˚ . In our convention, positive U1 places AA at the potential maximum with respect to AB and BA, leading to a shrinking AA region. Positive z1 also places AA at the maximum interlayer spacing from the underlying substrate. On the other hand, the bottom row corresponds to U2 = 3 meV, z2 = 0.01A˚ . Non-zero asymmetric terms account for the inequivalence between AB and BA sites. Positive U2 makes BA an energetically favored configuration, while positive z2 brings carbon atoms on AB to lie above those on BA. A typical deformation in graphene would be a superposition of both profiles with comparable strengths. the U0 for the average contribution, the U1 term for the symmetric and U2 for asymmetric contributions. This function can be used to describe any moire pattern including the interlayer coupling Hamiltonians, the interaction potentials, or the strain fields. Neglecting the effects from U0 the interlayer potentials can be written as: ε + θ (zˆ × ε˜ 2Ag 1) ∇[U1 f1 (r) + U2 f2(r)] ≈ 2µg ∇2u − 1 2 (zˆ × ∇)(∂xuy − ∂yux) + λg∇(∇ · u), (10) WWW.NATURE.COM/NATURE | 8 In the above formula we are treating the deformation within the first harmonics, which dominates the relaxation profile for small atomic displacements (|u| a) which accounts for the approximate sign in 10. Since the interlayer potential is completely defined by the first harmonics (Eq. 10, LHS), the solution for the displacement vector u(r) = ux xˆ + uy yˆ + h(r) zˆ is also constrained by moire´ periodicity, containing only the first harmonics components. The displacement and the corresponding in-plane expansion of carbon site at r is thus given by: u(r) = ks |ε 3 | G2ε˜ 2 ∇ + χRθ (zˆ × ∇) [U1 f1(r) +U2 f2(r)] + ε z1 f1(r) + z2 f2(r) zˆ, (11) 1 1 a/a0 ≈ 1 + 2 Tr[ui j(r)] = 1 + 2 ∂ ux + ∂ uy ∂x ∂y = 1 − ks |ε 3 | 2ε˜ 2 [U1 f1(r) + U2 f2(r)], (12) where U1 (U2) parametrize the inversion-symmetric (asymmetric) components of in-plane deformation, and z1 (z2) quantifies similar components for the deformation in the out-of-plane direction, i.e. corrugation. Meanwhile, the coefficient ks which relates u(r) with the potentials, and the factor χR which sets the magnitude of the curl-like terms are as follows ks = √ 2 3 3ε2a2(λg + 2µg) = 0.029 meV−1, χR = 2 + λg µg ≈ 2.34, (13) It can be seen in the the subsequent analysis that ks serves as the coefficient which encodes the contribution of in-plane deformation which is parametrized by U1 and U2, on the resulting gaps. We also assume the following definition for the symmetric and asymmetric components of the deformation:      U0 UAA 6(UAA +UAB +UBA)     1  U1 = M UAB =     18   2UAA − UAB − UBA √ ,   (14) U2 UBA 3(UAB −UBA) in which we take AA stacking point as the origin of our coordinate system, and the matrix elements of M are obtained in a way consistent with our definitions of f1(r) and f2(r), and in accordance with the convention in the literature [8, 11]. The symmetric and antisymmetric WWW.NATURE.COM/NATURE | 9 profiles are plotted in Fig. S4. We also use similar definition of symmetric and antisymmetric components for the corrugation and moire´ electronic couplings. Within this choice of coordinates, we can expect from symmetry consideration alone that a sizable primary gap only opens up when there is a relatively large symmetric component in the deformation (the asymmetric terms vanish after doing the mBZ average). In general, since in our choice of coordinate system the inversion-asymmetric electronic couplings dominate (V˜ V ), this implies that the dominant contribution to the global mass term results from symmetric deformations. Since the secondary Dirac cone gap is not a simple average over the mBZ, we expect both the symmetric and asymmetric contributions to be equally important. Moire´ couplings and gaps in G/BN The dynamics of low-energy electrons in G/BN can be understood as a relativistic Dirac particle subject to three types of perturbations: (1) H0(r), describing the periodic electrostatic potential, (2) Hz(r), acting as a local mass term which arises from the sublattice asymmetry, and (3) Hxy(r), quantifying the gauge field due to asymmetry in hoppings induced by the interaction with the BN substrate. For graphene electrons on valley K, the equation of motion can be expressed as follows HK(r) = υ p · σ + H0(r)1 + Hz(r)σ3 + ε(zˆ × 1) − θ ε˜ G ∇Hxy(r) · σ , (15) It is important to note that only the Hxy contribution is scaled by factors which depend on the twist angle. The term in Eq. 15 which is proportional to θ can be gauged away so that the influence of Hxy on the electronic structures is scaled by χθ = ε/ε˜. The three pseudospin components can be written in a compact manner using the following matrix notation: H0(r) Hxy(r) Hz(r) = 1 f1(r) f2(r) W, (16) where W stands for the electronic coupling matrix. Typically each pseudospin term for zero pressure has oscillating amplitudes on the order of ∼50 meV. Around AA, DFT parametriza- WWW.NATURE.COM/NATURE | 10

... pages left unread,continue reading

Document pages: 19 pages

Please select stars to rate!


0 comments Sign in to leave a comment.

    Data loading, please wait...