Combination of Wannier Functions and Continuum Electrostat-

In this section we explain our approach to derive appropriately screened Coulomb-interaction matrix elements for electrons in two-dimensional materials (e.g. graphene) and their heterostructures on the basis of Coulomb-interaction matrix elements from parent three-dimensional bulk systems (e.g. graphite). To this end, we first recall the continuum electrodynamic description of layered materials, free standing monolayers as well as heterostructures. Afterwards, the continuum formulation is embedded into a quantum lattice description in terms of localized Wannier functions.

Continuum Electrostatic Description

The problem of screening in terms of continuum electrostatics in a layered material is illustrated in Fig. 3.14 and has been considered in Refs. [199, 200, 201]. We aim to generalize these works to include non-local screening effects within the dielectric layer. Therefore, we will replace the dielectric constant ε1 by an appropriate dielectric function ε1pq_kq. To this end we have to relate the bulk dielectric function εBpqq to ε₁pq_kq. We assume that the former has been determined from first principles. As we will explicate later, we only modify the leading eigenvalue of the microscopic dielectric matrix derived from the dielectric kernel εBpr,r¹q so that we may assume εBpqqto be a scalar function. The embedding is not as easily done as in Refs. [201, 200, 199], because we face the problem of non-localities, in particular the qz-dependence which describes the periodicity of the bulk material in z-direction. Clearly, an assumption on how non-localities translate from bulk to monolayers or heterostructures has to be made. Here, we continue with the simplest possible approximation and neglect all non-localities inz-direction, i.e. we replace εBpqq by

ε1pq_kq “ h 2π

żπ{h

´π{h

dqz εBpq_k, qzq, (3.12) wherehplays the role of an effective layer thickness. This definition is plausible as the two-dimensional embedding breaks the periodicity in z-direction and only the z-local

(a) dielectric interface (b) dielectric layer (c) screening within a di-electric layer

Figure 3.15.:

Dielectric problems, which can be solved with the help of image charges. (a)A single dielectric interface. (b) Two dielectric interfaces separated byh and (c)the arising effective dielectric function within theε₁ area for a constant and finiteε₁pq_kq.

term ofεB should remain relevant. In this sense, Eq. (3.12) gives this local term as the Fourier transformation to the center of the monolayer, i.e., z “0. A similar formula was used in Ref. [197] to define a “two-dimensional macroscopic dielectric function”.

We now assume thatε1pq_kqis constant on the whole width (or height) of the mono-layer, i.e., for |z| ď h{2, and consider two dielectric materials on both sides with dielectric constants ǫ2 and ǫ3 according to Fig. 3.14, which gives

εpq_k, zq “

$&

%

ε2 z ą ^h2

ε1pq_kq |z| ď ^h₂ ε3 z ă ^´h2

. (3.13)

To find the appropriately screened interactionU^2Dpq_kq between two electrons in the central layer, we consider the electrostatic problem of an oscillating two-dimensional charge densityρpq_k, zq “ρ^2Dpq_kqδpzqin the center of the monolayer, which is atz “0.

The resulting electrostatic potential Φpq_k, zq is not the same as in the bulk due to modified screening and the fact that we are now considering a two-dimensional charge density confined to z “0.

In the dielectric continuum model, the modification of the screening is caused by the formation of image charges at the interfaces between the dielectrics. Let us first consider a single interface like it is illustrated in Fig. 3.15 (a), where the dielectric constant changes from ǫ1 to ǫ2, and a test charge q1 is positioned in region 1. For an observer in region 1, the induced polarization has the form of an image charge of magnitude q2 “ q1pǫ1 ´ǫ2q{pǫ1 `ǫ2q that is located in region 2 at an equal distance from the interface as the test charge [89]. If one has more than one interface, as in the case of our dielectric model and like depicted in Fig. 3.15 (b), the charge is reflected infinitely many times (as light between two parallel optical mirrors), giving rise to an infinite number of image charges located at

zplq “l¨h with l“ ˘0,˘1,˘2, ... (3.14)

perpendicular to the interfaces with ever decreasing magnitudes qjplq “ q1

ˆǫ₁ ´ǫj

ǫ1 `ǫj

˙|l|

(3.15) with j “ 2 and j “ 3 for the upper (right) and lower (left) interface, respectively.

Several works have treated this situation, a dielectric monolayer sandwiched between two semi-infinite dielectric materials [199, 200, 201]. We use here a special case of a formula derived in Ref. [201], giving the effective two-dimensional dielectric function in momentum space

ε^2D_effpq_kq “ ε1pq_kq“

1´ε˜2pq_kqε˜3pq_kqe^´2qkh‰ 1`“

ε˜2pq_kq `ε˜3pq_kq‰

e^´qkh`ε˜2pq_kqε˜3pq_kqe^´2qkh, (3.16) where the ratio

ε˜jpq_kq “ ε1pq_kq ´εj

ε1pq_kq `εj

(3.17) has been introduced. With this equation, the two-dimensional screened interaction (in the long-wavelength limit, see below) can be written as

U^2Dpq_kq “ v^2Dpq_kq

ε^2D_effpq_kq, (3.18)

where v^2Dpq_kq “v^3Dpq_k, z “0qis the bare interaction in the 2D system.

By evaluating Eq. (3.16) with the help of Eq. (3.12) we are now able to cal-culate the screened interaction between electrons in the free standing or embedded two-dimensional monolayer directly from the three-dimensional layered bulk proper-ties with the main approximation being the neglect of all non-localiproper-ties of dielectric response in the vertical direction. We will assess the quality of this approximation in section 3.3.3.

Wannier Function Based Formulation

The generalized Hubbard model from Eq. (3.11) involves matrix elements in terms of Wannier functions which has to be linked to the continuum electrostatic description developed in the previous section. For real materials like graphene or transition metal dichalcogenides these models involve multiple Wannier orbitals per unit cell and the corresponding Coulomb matrix elements must in general not be restricted to density-density interaction terms. Then we have to deal with full tensorial representations of the interactions U_αβγδ^kk1q which depend in general on two initial momenta k, k¹ and the momentum transfer q and four orbitals as defined in Eq. (2.133). In this case

further approximations are helpful to derive from the bulk Coulomb interaction the corresponding monolayer terms. First, we will neglect the non-vanishing overlap be-tween Wannier functions in different unit cells, which means tracing out the k and k¹ dependencies. Then, the Coulomb interaction depends on momentum transfer q but not on the initial momenta k and k¹.

To combine the macroscopic electrostatic description of the previous section with the representation in a Wannier basis, we represent the bare and screened Coulomb interaction as well as the dielectric function as quadratic matrices using generalized indices α˜ “ tα, δu and β˜“ tβ, γu: Uαβγδpqq ” U_α˜β˜pqq. The resulting matrix elements are used within

Hee“ ÿ

˜ αβq˜

U_α˜β˜pqq ÿ

kσ

c^:_ασpk´qqcδσpkq ÿ

k¹σ¹

c^:_βσ1pk¹ `qqcγσ¹pk¹q (3.19)

“ ÿ

˜ αβq˜

U_α˜β˜pqqnˆα˜pqqˆnβ˜p´qq. (3.20) and might be interpreted as interaction energies between generalized charge-density waves

nα˜pqq “ÿ

kσ

@c^:_ασpk´qqcδσpkqD

and nβ˜p´qq “ ÿ

k¹σ¹

A

c^:_βσ1pk¹ `qqcγσ¹pk¹qE

. (3.21) In Fourier space the leading eigenvalue and eigenvector of U_α˜β˜pqq correspond to the charge-density modulations with the longest wavelength, i.e. those charge-density waves where screening effects due to the environment are supposed to be strongest.

This statement can most easily be understood in real space. To this end we utilize the real-space representation of the Coulomb matrix

U_α˜β˜pRq “ 1 Nq

ÿ

q

U_α˜β˜pqqe^´iqR (3.22) for a simplified situation in which solely density-density interactions are taken into account in a basis consisting of two Wannier functions (α P tA, Bu) only. In such a situation the Coulomb-interaction tensor is reduced to a matrix

Uij “

ˆUAApR_ijq UABpR_ijq UBApR_ijq UBBpR_ijq

˙

, (3.23)

which might describe graphene’s pz orbitals on the sub-lattices A and B. Taking into account that Uij has to be symmetric (UBA“UAB) and setting UAA“UBB (which is a reasonable assumption in the case of graphene) the Coulomb matrix can readily be diagonalized

Uij “

ˆUAApR_ijq UABpR_ijq UABpR_ijq UAApR_ijq

˙

“ 1 2

ˆ1 ´1 1 1

˙ ˆU1pR_ijq 0 0 U2pR_ijq

˙ ˆ 1 1

´1 1

˙

(3.24)

using the leading (U1) and the second eigenvalue (U2)

U1pR_ijq “ UAApR_ijq `UABpR_ijq (3.25) U2pR_ijq “ UAApR_ijq ´UABpR_ijq. (3.26) By plugging Eq. (3.24) into the non-local interaction terms (i‰j) of the Hamiltonian from Eq. (3.11) we find

H_ee^i‰j “ 1 4

ÿ

i‰j σ,σ¹

rU1pR_ijq pnˆiAσ `nˆiBσq pnˆjAσ<sup>1</sup>`nˆjBσ¹q (3.27)

`U₂pRijq pnˆiAσ´nˆiBσq pnˆjAσ¹ ´ˆnjBσ¹qs

“ 1 4

ÿ

i‰j

rU1pR_ijqnˆinˆj `U2pR_ijqδˆniδˆnjs (3.28) with

nˆi “ÿ

σα

nˆiασ and δˆni “ÿ

σ

pnˆiAσ´ˆniBσq. (3.29) Here, nˆi measures thetotal charge of thei-th unit cell, whileδˆni accounts for its inter-nal charge variations. The former can be linked to electric monopoles and the latter to electric dipoles. Hence, the leading eigenvalue U1 describes unavoidable Coulomb penalties of monopole-like interactions between different unit cells. In contrast,U2 cor-responds to Coulomb penalties due to charge variationswithin the unit cells. Therefore, the Coulomb eigenvalues are connected to different length scales: U1 renders effects on macroscopic (inter-cell) distances, while U2 is linked to microscopic (intra-cell) details. Correspondingly, the leading eigenvalue is indeed inseparably connected to charge fluctuations with the longest wavelength, which holds for arbitrary situations due to the symmetric and purely real form of the Coulomb matrix⁶. By definition, continuum medium electrostatics describes the macroscopic response of a medium to long-wavelength electric field / potential variations. We thus correct the leading eigen-value ofU_α˜β˜pqqaccording to the algorithm from the previous section, while we assume for all other eigenvalues the same screening as in the bulk.

The full algorithm which we refer to as “Wannier Function Continuum Electro-static” (WFCE) approach can be divided into several steps, which are illustrated in the flowchart of Fig. 3.16, to obtain the partially screened Coulomb interaction of a two-dimensional sub-system directly from its three-dimensional host. We start with the bare,v^3Dpqq, and partially screened Coulomb-interactionmatrices U^3Dpqqobtained

6 Being symmetric and real leads to the fact that the leading eigenvector will consists of positive entries only. Thus, the leading eigenvalue will always be connected to the sum of orbital and spin resolved occupations ˆniασ and not to their (partial) differences.

v^3Dpqq,U^3Dpqq ñ ε^3Dpqq ““

v^3Dpqq‰1{2“

U^3Dpqq‰´1“

v^3Dpqq‰1{2

qz Fourier transform:

v^3Dpq, z “0q “N_q^´1_z ř

qzv^3Dpqqe^iqz¨0 ε^3Dpq, z “0q “ N_q^´1_z ř

qzε^3Dpqqe^iqz¨0

v^2Dpq_kq “v^3Dpq_k, z“0q, ε^3Dpq_k, z“0q

Decomposition:

v^2Dpq_kq “v^2D₁ pq_kqˇ

ˇv₁^2Dpq_kqD @

v₁^2Dpq_kqˇ

ˇ`v_Rest^2D pq_kq ε^3D₁ pq_k, z “0q “ @

v₁^2Dpq_kqˇ

ˇε^3Dpq_k, z “0qˇ

ˇv₁^2Dpq_kqD ε^3Dpq_k, z “0q “ ε^3D₁ pq_k, z “0qˇ

ˇv^2D₁ pq_kqD @

v₁^2Dpq_kqˇ

ˇ`ε^3D_Restpq_k, z “0q

2D model dielectric matrix:

ε^2D₁ pq_kq “ε^2D_eff “

ε^3D₁ pq_k, z “0q, h‰ ε^2D_Restpq_kq “ε^3D_Restpq_k, z“0q ε^2Dpq_kq “ε^2D₁ pq_kqˇ

ˇv₁^2Dpq_kqD @

v₁^2Dpq_kqˇ

ˇ`ε^2D_Restpq_kq

ε^2Dpq_kq ñ U^2Dpq_kq ““

v^2Dpq_kq‰1{2“

ε^2Dpq_kq‰´1“

v^2Dpq_kq‰1{2

q_k Fourier transform:

U^2Dpr_kq “ _N¹_q

k

ř

q_kU^2Dpq_kqe^iqkrk

U^2Dpr_kq

Figure 3.16.:

Flowchart of the Wannier function continuum electrostatics (WFCE) algorithm to obtain the screened Coulomb matrix elements of a freestanding monolayer directly from the Coulomb interaction of the corresponding layered bulk material.

from the ab initio calculations for the three-dimensional bulk of the layered material.

With these quantities we define the symmetric 3D dielectric function ε^3Dpqq ““

v^3Dpqq‰1{2“

U^3Dpqq‰´1“

v^3Dpqq‰1{2

. (3.30)

We now aim to link this dielectric matrix to the interactions taking place between electrons within one monolayer and connect it to model dielectric functions derived in the context of continuum electrostatics. To this end, we consider the bare intra-layer electronic interaction

v^3Dpq_k, z“0q “ 1 Nqz

ÿ

qz

v^3Dpq_k, qzqe^iqz¨0, (3.31) which is the same in the bulk layered material and in the monolayer, bilayer, etc. [i.e.

v^3Dpq_k, z “ 0q “ v^2Dpq_kq]. Nqz is the number of points used in the qz-summation.

v^2Dpq_kqcan be decomposed exactly v^2Dpq_kq “v₁^2Dpq_kqˇ

ˇv₁^2Dpq_kqD @

v₁^2Dpq_kqˇ

ˇ`v_Rest^2D pq_kq, (3.32) where v₁^2Dpq_kq is the leading eigenvalue of v^2Dpq_kq, ˇ

ˇv₁^2Dpq_kqD

is the corresponding eigenvector and v^2D_Restpq_kq contains the rest. Using the leading eigenvector of the bare Coulomb interaction, we are able to perform an analogous exact decomposition of the intra-layer 3D dielectric matrix

ε^3Dpq_k, z “0q “ ε^3D₁ pq_k, z “0qˇ

ˇv^2D₁ pq_kqD @

v₁^2Dpq_kqˇ

ˇ`ε^3D_Restpq_k, z “0q, (3.33) where the “head element” ε^3D₁ pq_k, z“0qof the dielectric matrix is defined as

ε^3D₁ pq_k, z “0q “ @

v₁^2Dpq_kqˇ

ˇε^3Dpq_k, z “0qˇ

ˇv₁^2Dpq_kqD

. (3.34)

To obtain the dielectric matrix of the 2D system ε^2Dpq_kq we replace ε^3D₁ pq_k, z “0q in Eq. (3.33) by the model dielectric function ε^2D_effpq_kqfrom Eq. (3.16):

ε^2Dpq_kq “ε^2D_effpq_kqˇ

ˇv₁^2Dpq_kqD @

v₁^2Dpq_kqˇ

ˇ`ε^2D_Restpq_kq (3.35) with ε1pq_kq from Eq. (3.16) being set to ε^3D₁ pq_k, z “ 0q from Eq. (3.33), while the rest matrix remains unchanged. Here, we use the fact that the microscopic (short-range) screening properties of the monolayer should be very similar to that of the bulk (ε^2D_Rest “ε^3D_Rest). We will later see that this is a good approximation.

Together with the bare intra-layer interaction the screened intra-layer interaction is given by

U^2Dpq_kq ““

v^2Dpq_kq‰1{2“

ε^2Dpq_kq‰´1“

v^2Dpq_kq‰1{2

. (3.36)

dielectric function bare and screened Coulomb interaction

ε3D 1pqq

v3D 1pqq,U3D 1pqq(eV) _qz“0.13

qz“0.25 qz“0.38 qz“0.50 qz“0

qk“0 vab initio

Uab initio

vanalytical

100 200 300 400 500

q(Å^´1) q(Å^´1)

00.1 0.2 0.3 0.4 0.5 0.6 0.7 1.60 0.5 1 1.5

1.8 2 2.2 2.4 2.6 2.8 3 3.2

Figure 3.17.:

Left: Leading eigenvalues of the bare and screened Coulomb matrix elements of graphite for q_z “ 0 obtained from ab initio calculations together with the analytical description of the unscreened interaction (dashed blue line). The (light) gray area indicates the interval of all other eigenvalues of the (bare) screened Coulomb matrix. Right: Momentum-dependent leading eigenvalue of dielectric function of graphite obtained from cRPA calculations for different values ofq_z. The red markers forq “0indicate the parallel (circle) and perpendicular (square) limits of the screening.

From U^2Dpq_kq a Fourier transformation with respect to q_k finally leads to screened Coulomb matrices in real space

U^2Dpr_kq “ 1 Nq_k

ÿ

q_k

U^2Dpq_kqe^iqkrk, (3.37) which can be used in extended Hubbard models like given in Eq. (3.11). Here, Nq_k is the number of points used in the qk-summation.

Im Dokument Electronic Structure of Novel Two-dimensional Materials and Graphene Heterostructures (Seite 103-110)