Microscopic screening theory for excitons in two-dimensional materials: A bridge between effective models and ab initio descriptions

Abstract

We present a computational approach for exciton calculations in two-dimensional (2D) materials within the Bethe-Salpeter equation (BSE) framework, employing an atomistic description with point-like orbitals. Unlike widespread efficient calculations that rely on classical or effective interaction models, such as the Rytova-Keldysh model, our method incorporates quantum screened interactions. By explicitly computing the 2D dielectric function at the random-phase approximation level, we capture screening effects beyond such approximations with an accuracy akin to first-principles methods. Consequently, we can realistically estimate excitonic binding energies with a bearable computational cost. A detailed account of the various convergence parameters sheds light on a possible cause of the large dispersion of binding energies reported in the literature using first-principles GW/BSE implementations. This work thus provides an alternative pathway towards efficient and faithful dielectric screening and exciton computations in low-dimensional materials.

Figures (6)

• Figure 1: Band structures of monolayer hBN and monolayer MoS2.a/b Band structure of monolayer hBN/MoS2 obtained through DFT calculations, using the exchange--correlation functional HSE06. For both materials, the bandgap is located at the K point of the BZ, and it is $\Delta=6.08$ eV for hBN and $\Delta=2.08$ eV for MoS2. For a better visualization, only the first few bands around the bandgap of each material are displayed.
• Figure 2: Head matrix element of the polarizability. Convergence of the head matrix element of the non-interacting polarizability with the number of conduction bands included in the sum in Eq. \ref{['eq:Chi_0_static_results']}. Panel a is for hBN and $\mathbf{q} = (0.5,0)$ Å$^{-1}$, while panel b is for MoS2 and $\mathbf{q} = (0.2,0)$ Å$^{-1}$. Both polarizability elements are converged with the number of momentum points in the BZ and appear in eV$^{-1}$.
• Figure 3: Head matrix element of the polarizability as a function of momentum. Head matrix element of the irreducible polarizability $\chi^0_{\mathbf{0} \mathbf{0}}(\mathbf{q})$, as per Eq. \ref{['eq:Chi_0_static_results']}, as a function of momentum magnitude, for hBN in panel a (scaled by a factor of $10$) and for MoS2 in panel b. The polarizability was computed in both cases along the symmetry line $\Gamma - \mathrm{K}$. The momentum $q=|\mathbf{q}|$ is in Å$^{-1}$ and the polarizability in eV$^{-1}$. Note the different scales for the vertical and horizontal axes in each plot. In both panels, all the valence and conduction bands were included in the sum \ref{['eq:Chi_0_static_results']}.
• Figure 4: Macroscopic dielectric function. Comparison between the effective 2D dielectric function obtained in this work with ab initio methods. Panels a and b display the results for hBN and MoS2, respectively. The ab initio results were obtained using the QEH package for Python, which can be found in the literature Andersen_Latini_Thygesen_2015, and they are represented by the full black curves. The macroscopic dielectric functions obtained using the strictly 2D description developed in this work as per Eq. \ref{['eq:macroscopic_epsilon_def_results']} are represented in red. The blue curves represent the same quantity, but when obtained through the Q2D approach that we also introduce in this work, explained with detail in Section \ref{['subsubsection:Q2D_dielectric_function']}. For hBN, in panel a, we applied the cutoff $G_{\mathrm{c}} = 3$ Å$^{-1}$ for the size of the dielectric function. This accounts for seven reciprocal lattice vectors, including the null vector. For MoS2 we applied the same cutoff, which results in seven $\mathbf{G}$ vectors as well, but we intentionally discarded two reciprocal lattice vectors in the Q2D approach for the blue curve to better agree with the ab initio one. We found that including all vectors underestimated the screening with respect to the black curve. We compare as well the full numerical $q$-dependence with the Rytova--Keldysh model in dashed black, with the $r_0$ parameter estimated by our numerical curves (forcing $\varepsilon_{\mathrm{RK}}(q=0)=1$). For the Q2D results, we used for the monolayer thicknesses the values provided by the QEH package for each material.
• Figure 5: Convergence of the excitonic binding energy with the BZ mesh size and regularization radius. Panels a1-d1 present the results pertaining to hBN, while panels a2-d2 pertain to MoS2. Panels a1 and a2 show the ground state excitonic energy level $E_X$ as a function of $N$ for different realizations of $\varsigma$, which defines the radius of the regularization region as $q_0=\varsigma k_0$, $k_0$ being the length of the momentum vector closest to the origin in the BZ. Panels b1 and b2 show the dependence of $E_x$ on $\varsigma^{-1}$ for different realizations of $N$, which follows a linear model as $E_X = m_N \varsigma^{-1} + b_N$. The plot in panels c1 and c2 shows the behavior of the slope $m_N$ with $N$, whereas panels d1 and d2 plot the origin off-set $b_N$ versus $N$.