The Quantum Electrostatic Heterostructure (QEH) model

Brief overview

In this tuturial we show how to calculate the linear response of a general van der Waals Heterostructure (vdWH) by means of the quantum electrostatic heterostructure (QEH) model. This method allows to overcome the computational limitation of the standard ab-initio approaches by combining quantum accuracy at the monolayer level and macroscopic electrostatic coupling between the layers. More specifically, the first step consists in calculating the response function of each of the layers composing the vdWHs in the isolated condition and encoding it into a so called dielectric building block. Then, in the next step the dielectric building blocks are coupled together through a macroscopic Dyson equation. The validity of such an approach strongly relies on the absence of hybridization among the layers, condition which is usually satisfied by vdWHs.

A thourough description of the QEH model can be found in [1] and [2]:

Constructing a dielectric building block

First, we need a ground state calculation for each of the layers in the vdWH. We will consider a MoS2/WSe2 bilayer. In the following script we show how to get the ground state gpw-file for the MoS2 layer:

from ase.build import mx2
from gpaw import GPAW, PW, FermiDirac


structure = mx2(formula='MoS2', kind='2H', a=3.184, thickness=3.127,
                size=(1, 1, 1), vacuum=5)

ecut = 600

calc = GPAW(mode=PW(ecut),
            xc='LDA',
            kpts={'size': (30, 30, 1), 'gamma': True},
            occupations=FermiDirac(0.01),
            parallel={'domain': 1},
            txt='MoS2_gs_out.txt')

structure.calc = calc
structure.get_potential_energy()

calc.diagonalize_full_hamiltonian(nbands=500, expert=True)
calc.write('MoS2_gs_fulldiag.gpw', 'all')

The gpw-file for WSe2 can be obtained in a similar way. The gpw-files stores all the necessary eigenvalues and eigenfunctions for the response calculation.

Next, we need the \(z\)-dependence of the full response function \(\chi(q,\omega, G_z, G'_z)\). Essentially, a Dyson equation for the isolated layer is solved to obtain the full response function \(\chi(q,\omega,G,G')\), and then the desired \(G_z\)-elements are chosen from its basis. This is all handled by the ChiHandler object, and saved to the file *-chi.npz. Subsequently, the method interpolate_chi_to_bb will interpolate the full response function onto a dense \(q\)-grid, and calculate the monopole and dipole component of the response function and the corresponding induced densities. These quantities are then condensed into the dielectric building block *-bb.npz. Here is how to get the MoS2 and building block:

from pathlib import Path
from gpaw.mpi import world
from gpaw.response.df import DielectricFunction
from gpaw.response.qeh import QEHChiCalc

from qeh.bb_calculator.chicalc import ChiHandler
from qeh.bb_calculator.bb_builder import interpolate_chi_to_bb

df = DielectricFunction(calc='MoS2_gs_fulldiag.gpw',
                        eta=0.001,
                        frequencies={'type': 'nonlinear',
                                     'domega0': 0.05,
                                     'omega2': 10.0},
                        nblocks=8,
                        ecut=150,
                        truncation='2D')

chicalc = QEHChiCalc(df)
chihandler = ChiHandler('MoS2', chicalc, q_max=3.0)
chihandler.calculate_chi_2d()

if world.rank == 0:
    interpolate_chi_to_bb('MoS2', aN=2)
    Path('MoS2_gs_fulldiag.gpw').unlink()

The same goes for the WSe2 building block. Once the building blocks have been created, they need to be interpolated to the same kpoint and frequency grid. This is done as follows:

from qeh.bb_calculator.bb_interpolator import interpolate_building_blocks
interpolate_building_blocks(BBfiles=['MoS2', 'WSe2'])

Specifically, this interpolates both building block to lowest \(q\)-grid of the two.

Finally the building blocks are ready to be combined electrostatically.

Interlayer excitons in MoS2/WSe2

As shown in [3] the MoS2_WSe2 can host excitonic excitations where the electron and the hole are spatially separated, with the electron sitting in the MoS2 layer and the hole in the WSe2 one.

Because of the finite separation distance we expect that the electron-hole screened interaction between the two particles will not diverge when the in-plane separation goes to 0. To illustrate this we show how to calculate the screened electron-hole interaction using the QEH model and the building blocks created above:

# web-page: W_r.png
import numpy as np
import matplotlib.pyplot as plt
from qeh import QEH


thick_MoS2 = 6.2926
thick_WSe2 = 6.718

d_MoS2_WSe2 = (thick_MoS2 + thick_WSe2) / 2
inter_mass = 0.244

HS = QEH.heterostructure(BBfiles=['MoS2-int', 'WSe2-int'],
                         layerwidth_n=[thick_MoS2, thick_WSe2],
                         wmax=0,
                         amax=2)

hl_array = np.array([0., 0., 1., 0.])
el_array = np.array([1., 0., 0., 0.])


# Getting the interlayer exciton screened interaction on a real grid
r, W_r = HS.get_exciton_screened_potential_r(
    r_array=np.linspace(1e-1, 30, 1000),
    e_distr=el_array,
    h_distr=hl_array)

plt.plot(r, W_r, '-g')
plt.title(r'Screened Interaction Energy')
plt.xlabel(r'${\bf r}$ (Ang)', fontsize=20)
plt.ylabel(r'$W({\bf r})$ (Ha)', fontsize=20)
plt.savefig('W_r.png')
plt.show()

Here is the generated plot:

As expected the interaction does not diverge.

If one is also interested in the interlayer exciton binding energy, it can be readily calculated by appending the following lines in the script above:

ee, ev = HS.get_exciton_binding_energies(eff_mass=inter_mass,
                                         e_distr=el_array,
                                         h_distr=hl_array)

print('The interlayer exciton binding energy is:', -ee[0])

We find an interlayer exciton binding energy of \(\sim 0.3\) eV.