Hamiltonian

The orthogonal tight-binding Hamiltonian, \(H(\mathbf{k})\), obtained via the Wannier90 framework [Arash2008] can be written as follows,

(1)\[\begin{equation} H(\mathbf{k}) = H_{0}+ \sum_{i=1}^{N} e^{i \mathbf{k} \cdot \mathbf{R_{i}}} H_{\mathbf{R_{i}}}~, \end{equation}\]

where \(H_{0}\) corresponds to the Hamiltonian in unit cell, which contains the on-site energies and hopping parameters inside the cell. \(H_{\mathbf{R_{i}}}\) corresponds to the hopping matrices, representing the interaction between unit cell and neighbors cells, while the matrix elements of \(H_{0}\) and \(H_{\mathbf{R_{i}}}\) are the output from the Wannier90 package [Wannier90] and the number of neighbors cell proportional to the \(\textbf{k}-mesh\) used for Wannier interpolation. In our code the values of \(H_{0}\) and \(H_{\mathbf{R_{i}}}\) are declared in the \(\texttt{PARAMS\_FILE}\).

The electronic energy levels are obtained as the following:

(2)\[\begin{equation} H(\mathbf{k}) |n,\mathbf{k} \rangle = E_{n,\mathbf{k}} |n,\mathbf{k} \rangle \end{equation}\]

where \(E_{n,\mathbf{k}}\) and \(|n,\mathbf{k} \rangle\) are the eigenvalues and eigenvectors respectively, \(n\) corresponds to the band index, that could be addressed as \(v\) for the valence (occupied) states and \(c\) for the conduction (unoccupied) states, \(\mathbf{k}\) are the \(\textbf{k}-point\) in Brillouin zone.

Optical properties

To obtain the optical properties, we calculated the real and imaginary parts of the frequency dependent dielectric tensor, \(\epsilon_{1,\alpha,\beta}(\omega)\) and \(\epsilon_{2,\alpha,\beta}(\omega)\), respectively, through the following expressions:

(3)\[\begin{align} \epsilon_{1,\alpha,\beta}(\omega) = \delta_{\alpha,\beta} + \frac{e^{2} S_{p}}{\epsilon_{0}\Omega N_{\mathbf{k}}} \sum_{\mathbf{k},c,v} F_{\alpha,\beta}^{c,v,\mathbf{k}} \frac{(E_{c,\mathbf{k}}-E_{v,\mathbf{k}})-\hbar \omega}{\left(\hbar \omega - (E_{c,\mathbf{k}}-E_{v,\mathbf{k}})\right)^{2}+\eta^{2}}~, \end{align}\]
(4)\[\begin{align} \epsilon_{2,\alpha,\beta}(\omega)= \frac{ e^{2} S_{p}}{\epsilon_{0}\Omega N_{\mathbf{k}}} \sum_{\mathbf{k},c,v} F_{\alpha,\beta}^{c,v,\mathbf{k}} \frac{\eta}{\left(\hbar \omega - (E_{c,\mathbf{k}}-E_{v,\mathbf{k}})\right)^{2}+\eta^{2}}~, \end{align}\]

where \(\delta_{\alpha,\beta}\) is a delta kroenecker, \(S_{p}\) is the spin factor, being \(1\) for SOC and spin polarized calculations and \(2\) for non-polarized calculations, \(F_{\alpha,\beta}^{c,v,\mathbf{k}}\), in the scope of independent particle approximation (IPA), corresponds to the oscillator force, defined by:

(5)\[\begin{equation} F_{\alpha,\beta}^{c,v,\mathbf{k}} = \frac{\langle c,\mathbf{k} |P_{\alpha}|v,\mathbf{k} \rangle \langle v,\mathbf{k} |P_{\beta}|c,\mathbf{k} \rangle}{\left(E_{c,\mathbf{k}}-E_{v,\mathbf{k}}-i\eta_{1}\right) \left(E_{c,\mathbf{k}}-E_{v,\mathbf{k}}+i\eta_{1}\right)}~, \end{equation}\]

where \(\Omega\) is the volume of unit cell, \(\epsilon_{0}\) is the vacuum permittivity constant, \(N_{\mathbf{k}}\) is the number of \(\textbf{k}-points\) employed for the BZ integration, \(\omega\) is the incident photon frequency, \(c(v)\) corresponds to the conduction(valence) states. \(\eta\) is a parameter to smooth the dielectric function, where \(\eta= \eta_{1}+\eta_{2}\), being \(\eta_{1}\) defined by the input flag \(\texttt{CSHIFT}\) and \(\eta_{2}\) automatically defined by the code through the separation energy between the transition states. \(\alpha\) and \(\beta\) corresponds to the components \(x, y\) and \(z\) in the dielectric tensor. \(P_{\alpha}\) corresponds to the light-matter interaction operator, which is given by,

(6)\[\begin{equation} P_{\alpha} = \frac{\partial H(\mathbf{k})}{\partial k_{\alpha}}~, \end{equation}\]

where \(H(\mathbf{k})\) corresponds to the electronic Hamiltonian, and \(\alpha=`{:math:`x,y,z\)}. For circularly polarized light, selected by the flag \(\texttt{CPOL= T}\), the light-matter interaction operator for \(\sigma_{\pm}\) are written by the following expression:

(7)\[\begin{equation} P_{\sigma_{\pm}} = \frac{1}{\sqrt{2}} \left( \frac{\partial H(\mathbf{k})}{\partial k_{x}}\pm i\frac{\partial H(\mathbf{k})}{\partial k_{y}}\right)~. \end{equation}\]

It’s well known the relation between the complex dielectric constant \(\epsilon(\omega)\) and the complex refractive index \(\tilde{n}(\omega)=n(\omega)+i\kappa(\omega)\), where \(n(\omega)\) is the refractive index and \(\kappa(\omega)\) the extinction coefficient. They are related by:

(8)\[\begin{equation} \epsilon(\omega) =\tilde{n}^{2}(\omega)~, \end{equation}\]

separating \(\epsilon(\omega)\) in real and imaginary part, \(\tilde{n}(\omega)\) in refractive index and extinction coefficient we have:

(9)\[\begin{equation} \epsilon_{1,\alpha,\beta}(\omega)+i\epsilon_{2,\alpha,\beta}(\omega) = n_{\alpha,\beta}^{2}(\omega) + 2 in_{\alpha,\beta}(\omega)\kappa(\omega)-\kappa_{\alpha,\beta}^{2}(\omega)~. \end{equation}\]

Then, we obtain the real and imaginary part of dielectric function, i.e.,

(10)\[\begin{equation} \epsilon_{1,\alpha,\beta}(\omega)= n_{\alpha,\beta}^{2}(\omega) - \kappa_{\alpha,\beta}^{2}(\omega)~, \end{equation}\]
(11)\[\begin{equation} \epsilon_{2,\alpha,\beta}(\omega)= 2 n_{\alpha,\beta}(\omega)\kappa_{\alpha,\beta}(\omega)~. \end{equation}\]

From these two equations we derive the expressions of the extinction coefficient and of the refractive index, as follows,

(12)\[\begin{equation} \kappa_{\alpha,\beta}(\omega)= \left[\frac{\sqrt{\epsilon_{1,\alpha,\beta}^{2}(\omega)+\epsilon_{2,\alpha,\beta}^{2}(\omega)}-\epsilon_{1,\alpha,\beta}(\omega)}{2}\right]^{\frac{1}{2}} \end{equation}\]
(13)\[\begin{equation} n_{\alpha,\beta}(\omega) = \left[\frac{\sqrt{\epsilon_{1,\alpha,\beta}^{2}(\omega)+\epsilon_{2,\alpha,\beta}^{2}(\omega)}+\epsilon_{1,\alpha,\beta}(\omega)}{2}\right]^{\frac{1}{2}}~. \end{equation}\]

The absorption coefficient \(A_{\alpha,\beta}(\omega)\) is defined as:

where \(c\) is the light speed. The reflectivity \(R_{\alpha,\beta}(\omega)\) is defined by the refractive index and extinction coefficient by the following expression:

(14)\[\begin{equation} R_{\alpha,\beta}(\omega) = \frac{(n_{\alpha,\beta}(\omega)-1)^{2}+\kappa_{\alpha,\beta}^{2}(\omega)}{(n_{\alpha,\beta}(\omega)+1)^{2}+\kappa_{\alpha,\beta}^{2}(\omega)}~. \end{equation}\]

The energy loss function \(L_{\alpha,\beta}(\omega)\) is defined as:

(15)\[\begin{equation} L_{\alpha,\beta}(\omega)=\frac{\epsilon_{2,\alpha,\beta}}{\epsilon_{1,\alpha,\beta}^{2}+\epsilon_{2,\alpha,\beta}^{2}}~. \end{equation}\]

Berry Curvature

In solid state physics, the crystal symmetry rules the electronic band structure, as well as the nature of Bloch states. For the crystal with inversion symmetry, the Berry curvatures of the Bloch states are zero [Xiao2010]. Nevertheless, in crystals with inversion asymmetry such as the \(\ce{MoS2}\) monolayer and the \(\ce{WS2}/\ce{MoS2}\) vdWs heterostructure, they are nonzero [Zhang2011], [Kormanyos2018] . Moreover, time reversal symmetry requires all Berry-phase related physical quantities to be valley-contrasting. To gain further insight into the electronic band structure and optical properties of the condensed-matter systems, let us explore the Berry curvature \(\Omega_{n}(\textbf{k})\), dictated by, [Cai2013].

(16)\[\begin{split}\begin{equation} \begin{split} & \Omega_{\rm n}(\bf k)= i \langle \nabla_{k} u_{n}\vert \times \vert \nabla_{k} u_{n}\rangle~, \\ \label{Berry-1} \end{split} \end{equation}\end{split}\]

where \(\vert u_{n} \rangle`is the :math:`n\)-th Bloch state.

In this code, one can obtain the total Berry curvature (semiconductors only) either along certain \(\textbf{k}\)-path \((\texttt{KPATH\_FILE})\), using the flag \(\texttt{BERRY= T}\), or the whole BZ, using the flag \(\texttt{BERRY\_BZ= T}\). Total Berry curvature \(\Omega_{\alpha,\beta}(\mathbf{k})\) are written by the following expression:

(17)\[\begin{equation} \Omega_{\alpha,\beta} (\mathbf{k}) = -2 \text{Im} \sum_{v,c} \frac{\langle c,\mathbf{k} |P_{\alpha}|v,\mathbf{k} \rangle \langle v,\mathbf{k} |P_{\beta}|c,\mathbf{k} \rangle}{\left(E_{c,\mathbf{k}}-E_{v,\mathbf{k}}\right)^{2}}~. \end{equation}\]

The code only outputs the non-zero \(\Omega_{\alpha,\beta}(\mathbf{k})\) components: \(\Omega_{x,y}(\mathbf{k})\), \(\Omega_{x,z}(\mathbf{k})\) and \(\Omega_{y,z}(\mathbf{k})\). More details of Berry curvature calculation in the scope of MLWF, can be found in Wang’s work [Wang2006].

Bethe–Salpeter Formalism

The Bethe–Salpeter equation (BSE) is a many-body equation that describes the interaction between an electron and a hole. The excitonic states are obtained through the solution of the two-particle problem via BSE [Dias2020], [Dias2021]. The excitonic Hamiltonian, \(H_{exc}\), is composed by the electron, \(H_{e}\), and hole, \(H_{h}\), single particle Hamiltonians plus Coulomb interaction potential, \(V_{eh}\), which binds the electron-hole pairs, i.e.,

(18)\[\begin{equation} H_{exc} = H_{e} + H_{h}+V_{eh}~. \end{equation}\]

The excitonic states with momentum center of mass \(\mathbf{Q}\) can be expanded in terms of the product of electron and hole pairs wave functions, as follows,

(19)\[\begin{equation} \Psi_{ex}^{n}(\mathbf{Q}) = \sum_{c,v,\mathbf{k}} A_{c,v,\mathbf{k},\mathbf{Q}}^{n} \ \left( |c,\mathbf{k}+\mathbf{Q} \rangle \otimes |v,\mathbf{k} \rangle \right)~, \label{eq:Exc_Basis} \end{equation}\]

where the index \(c\) and \(v\) corresponds to the conduction and valence bands states, with momentum \(\mathbf{k}+\mathbf{Q}\) and \(\mathbf{k}\), respectively. The problem of excitons eigenvalues, can be transformed into BSE, [Salpeter1951] using the expansion of Equation (19), it reads,

(20)\[\begin{eqnarray} \left( E_{c,\mathbf{k}+\mathbf{Q}}-E_{v,\mathbf{k}}\right) A_{c,v,\mathbf{k},\mathbf{Q}}^{n} + \frac{1}{N_{k}} \sum_{\mathbf{k'},v',c'} W_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}} \ A_{c',v',\mathbf{k'},\mathbf{Q}}^{n} = E^{n}_{\mathbf{Q}} A_{c,v,\mathbf{k},\mathbf{Q}}^{n} \label{bse}~, \end{eqnarray}\]

where \(E^{n}_{\mathbf{Q}}\) are the energy of n-th excitonic state with momentum \(\mathbf{Q}\), \(A_{c,v,\mathbf{k},\mathbf{Q}}^{n}\) are the exciton wave function (eigenvector), which are obtained by solving Equation (20). \(E_{c,\mathbf{k} + \mathbf{Q}} - E_{v,\mathbf{k}}\) are the single-particle energy difference between a conduction band state \(c\) with momentum \(\mathbf{k}+\mathbf{Q}\) and a valence band state \(v\) with momentum \(\mathbf{k}\), and \(W_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}}\) are the many-body Coulomb interaction matrix element, which can be divided into two parts, direct interaction, \(W^d\), and exchange interaction, \(W^x\), respectively, i.e.,

(21)\[\begin{equation} W_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}}=W^d_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}}+W^x_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}}~. \end{equation}\]

Since the Coulomb potential varies slightly inside unit cell in comparison with Bloch functions, we can approximate the orbital character of Coulomb term in the following way:

(22)\[\begin{equation} W^{d}_{(\mathbf{k},v,c),(\mathbf{k'},v',c'),\mathbf{Q}}=V(\mathbf{k}-\mathbf{k'}) \ \langle c,\mathbf{k}+\mathbf{Q}|c',\mathbf{k'}+\mathbf{Q} \rangle \ \langle v',\mathbf{k'}|v,\mathbf{k} \rangle \end{equation}\]

and

where, \(N_{k}\) is the number of \(\textbf{k}-points\) in the BZ.

In BSE formalism, the real and imaginary parts of the frequency dependent dielectric tensor \(\epsilon^{BSE}_{1,\alpha,\beta}(\omega)\) and \(\epsilon^{BSE}_{2,\alpha,\beta}(\omega)\) are obtained from the following expressions:

(23)\[\begin{align} \epsilon^{BSE}_{1,\alpha,\beta}(\omega) = \delta_{\alpha,\beta} + \frac{e^{2} S_{p}}{\epsilon_{0}\Omega N_{\mathbf{k}}} \sum_{n} F_{\alpha,\beta}^{n,BSE} &\frac{E_{0}^{n}-\hbar \omega}{\left(\hbar \omega - E_{0}^{n}\right)^{2}+\eta^{2}}~, \end{align}\]
(24)\[\begin{align} \epsilon^{BSE}_{2,\alpha,\beta}(\omega)= \frac{e^{2} S_{p}}{\epsilon_{0}\Omega N_{\mathbf{k}}} \sum_{n} F_{\alpha,\beta}^{n,BSE} \frac{\eta}{\left(\hbar \omega - E_{0}^{n}\right)^{2}+\eta^{2}}~, \end{align}\]

where \(E_{0}^{n}\) is the n-th direct \((\mathbf{Q}=0)\) excitonic state energy, \(F_{\alpha,\beta}^{n,BSE}\) is the excitonic modulated oscillator force, given by,

(25)\[\begin{align} F_{\alpha,\beta}^{n,BSE} = \left( \sum_{c,v,\mathbf{k}} \frac{A^{n}_{c,v,\mathbf{k},0}\langle c,\mathbf{k} |P_{\alpha}|v,\mathbf{k} \rangle}{\left(E_{c,\mathbf{k}}-E_{v,\mathbf{k}}+i\eta_{1}\right)} \right) \left(\sum_{c',v',\mathbf{k'}} \frac{A^{n*}_{c',v',\mathbf{k'},0}\langle v',\mathbf{k'} |P_{\beta}|c',\mathbf{k'} \rangle}{\left(E_{c',\mathbf{k'}}-E_{v',\mathbf{k'}}-i\eta_{1}\right)} \right)\label{eq:bse-fosc} ~. \end{align}\]

Finite temperature BSE approach

The code also calculate the exciton energy levels, dielectric functions and other optical properties at a finite temperature for the excitons with \(\mathbf{Q}=0\) by the flags \(\texttt{BSET= T}\), and \(\texttt{BSET\_BND= T}\) to calculate exciton band structure. These flags cannot be combined with their counterparts without thermal effects \(\texttt{BSE= T}\) and \(\texttt{BSE\_BND= T}\), respectively. To incoporate the temperature effects, we have to modify BSE as follows [Albrecht1998], [Marini2008], and [MolinaSanchez2016].

where \(f(E)\) is the Fermi-Dirac distribution, and \(\Xi_{c,\mathbf{k}+\mathbf{Q}}\) is the conduction band energy with a correction of temperature effects, defined as :

where \(T\) is the temperature in \(K\), given by the flag \(\texttt{TEMP}\), \(\alpha_{B}\) and \(\Theta\) are empirical parameters to be fitted to experimental data, given by the flags \(\texttt{ST}\) and \(\texttt{PHAVG}\), respectively, and the analytical expression of the correction \(\Delta(\alpha_{B},\Theta,T)\) depends on the the approximations made in the calculation. Usually, the following three formulations have been widely used in the literature. The first one is the frozen atom approximation in which we assume \(\Delta(\alpha_{B},\Theta,T)=0\), given by the flag \(\texttt{TA= FA}\). The second one is Bose-Einstein approximation [Liu2020] given by the flag \(\texttt{TA= BE}\), where

where \(\alpha_B\) represents the strength of electron-phonon interactions in si{electronvolt} and \(\Theta\) is the average phonon temperature in \(K\). The third one is given by the flag \(\texttt{TA= VE}\), in which the correction term is described by [ODonnell1991].

where \(\Theta\) is a dimensionless coupling constant and \(\alpha_B\) is an average of phonon energy in \(eV\).

These temperature effects also correct the expressions for the oscillator force \(F_{\alpha,\beta}^{c,v,\mathbf{k}}\) and excitonic modulated oscillator force \(F_{\alpha,\beta}^{n,BSE}\) by

Coulomb Potentials

Our code posses several options for the Coulomb interaction \(V_eh\) (reciprocal space) in the BSE formalism. In this section we’ll briefly discuss each of them. To avoid divergence in these potentials, a tolerance parameter \(\mathbf{Q}_{tol}\) was introduced in our implementation, its value given by the flag \(\texttt{KTOL}\).

Bare Coulomb Potential (V3D)

This is the traditional Coulomb Potential in reciprocal space, being identified by the input flag \(\texttt{COULOMB\_POT= V3D}\), and described by the following expression

(28)\[\begin{equation} V_{e h}(\mathbf{Q}) = \frac{4 \pi e^{2}}{\epsilon_{0} \mathbf{Q}^{2}}~, \end{equation}\]

where \(\epsilon_{0}\) is the vacuum permittivity.

(28)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2}}{2 V_{uc} \epsilon_{0} \mathbf{Q}^2}~, \end{equation}\]

where \(V_{uc}\) is the unit cell volume. To avoid the divergence at \(\mathbf{Q} \rightarrow 0\), we define

(28)\[\begin{equation} V(\mathbf{|Q|<Q_{tol}}) = 0. \end{equation}\]

2D Keldysh Potential(V2DK)

A Keldysh type Coulomb potential normally be used to properly describe the 2D Coulomb interaction [Ridolfi2018]. It’s selected by the flag \(\texttt{COULOMB\_POT= V2DK}\) in main input file, with the following equation

(32)\[\begin{equation} V(\mathbf{Q})=-\frac{e^{2}}{2 A_{uc} \epsilon_{0} \epsilon_{d} |\mathbf{Q}| (1+r_{0} |\mathbf{Q}|)}~. \label{keldysh-k} \end{equation}\]

with the objective to avoid Coulomb potential singularity, we make the following approximation for \(|\mathbf{Q}| \le \mathbf{Q}_{tol}\),

(32)\[\begin{equation} V(|\mathbf{Q}| \le \mathbf{Q}_{tol} ) = \frac{-e^{2} a N_{k}}{4 \pi A_{uc} \epsilon_{0} \epsilon_{d}} \left[\alpha_{1}+\alpha_{2} \Delta+\alpha_{3} \Delta^{2} \right], \end{equation}\]

where \(\Delta=2 \pi r_{0} /a N_{k}\), and \(\alpha_{1}=1.76\), \(\alpha_{2}=1\), and \(\alpha_{3}=0\) [Dias2020]. \(N_{k}\) is the number of \(\textbf{k}\)-points used in \(\textbf{k}\)-mesh, \(A_{uc}\) is the unit cell area and \(a\) is the average of the first and second lattice vectors size. The Coulomb potential screening length \(r_0\) is given by the expression:

(32)\[\begin{equation} r_0 = l_c \frac{(\epsilon_{2}-1)}{(\epsilon_{1}+\epsilon_{3})} ~, \end{equation}\]

and the effective dielectric constant \(\epsilon_{d}\) is given by:

(32)\[\begin{equation} \epsilon_{d} = \frac{\epsilon_{1}+\epsilon_{3}}{2} \end{equation}\]

where \(l_{c}\) is a free parameter, given by the flag \(\texttt{LC}\), that corresponds to the third lattice vector size or the distance between two neighbor layers of the 2D systems, \(\epsilon_{2}\) are in plane macroscopic effective dielectric constant, given in the main input file by the flag \(\texttt{EDIEL}\), \(\epsilon_{1}\) and \(\epsilon_{3}\) corresponds to the dielectric constants of the substrate above and bellow, given by the flags \(\texttt{EDIEL\_T}\) and \(\texttt{EDIEL\_B}\) respectively.

2D Rytova–Keldysh Potential (V2DRK)

2D Rytova–Keldysh Potential is another version of a Keldysh type Coulomb potential used for 2D systems [VanTuan2018]. It’s selected by the flag \(\texttt{COULOMB\_POT= V2DRK}\) in main input file, being described,

(35)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2}}{2 A_{uc} \epsilon_{0} |\mathbf{Q}| F(\mathbf{Q})} e^{-w_{0} |\mathbf{Q}|} ~, \end{equation}\]

where

(35)\[\begin{equation} F(\mathbf{Q}) = \frac{\left[1-\left(p_{b} \ p_{t} \ e^{-2 |\mathbf{Q}| \eta l_{c} }\right)\right] \kappa}{\left[1-\left(p_{t} e^{-\eta |\mathbf{Q}| l_{c}}\right)\right] \left[1-\left(p_{b} e^{-\eta |\mathbf{Q}| l_{c}}\right)\right]} + r_{0} |\mathbf{Q}| e^{-w_{0} |\mathbf{Q}|} \end{equation}\]

and

where \(A_{uc}\) is the area of crystal unit cell, \(w_{0}\) and \(r_{0}\) are free parameters given by the flags \(\texttt{W\_COUL}\) and \(\texttt{R\_0}\) respectively. \(\epsilon_{2}\) is the in plane effective dielectric constant of the crystal, given in the main input file by the flag \(\texttt{EDIEL}\), \(\epsilon_{z}\) is the out of plane effective dielectric constant of the crystal, given by the flag \(\texttt{EDIEL\_Z}\), \(\epsilon_{b}\) is the dielectric constant of the environment bellow the 2D system, given by the flag \(\texttt{EDIEL\_B}\), \(\epsilon_{t}\) is the dielectric constant of the environment above the 2D system, given by the flag \(\texttt{EDIEL\_T}\). \(l_{c}\) refers to the crystal thickness (in \(\mathrm{\mathring{A}}\)), given by the flag \(\texttt{LC}\).

To avoid the divergence at \(\mathbf{Q} \rightarrow 0\), we consider the following

(35)\[\begin{equation} V(\mathbf{|Q|<Q_{tol}}) = 0~. \end{equation}\]

2D truncated Coulomb Potential (V2DT)

A Coulomb potential truncated for 2D systems, proposed by Rozzi textit{et. al.} [Rozzi2006] being selected by the flag texttt{COULOMB_POT= V2DT} in main input file. This potential has the following expression:

(37)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2}}{2 V_{uc} \epsilon_{0} \mathbf{Q}^2} \left[ 1 - e^{-C_{2}} \left( C_{1} \sin(C_3) - \cos(C_3) \right) \right]~, \end{equation}\]

where \(V_{uc}\) is volume of the unit cell and \(\mathbf{r}_{3}\) is the size of the third lattice vector (real space). C:math:_1, C:math:_2 and C:math:_3 are defined as,

(37)\[\begin{split}\begin{eqnarray} C_{1} = \frac{|Q_{z}|}{\sqrt{Q_{x}^{2}+Q_{y}^{2}}}~,\nonumber \\ C_{2} = \frac{\mathbf{r}_{3}}{2} \sqrt{Q_{x}^{2}+Q_{y}^{2}}~, \nonumber \\ C_{3} = \frac{\mathbf{r}_{3}}{2} |Q_{z}|~. \end{eqnarray}\end{split}\]

If \(\sqrt{Q_{x}^{2}+Q_{y}^{2}} < |\mathbf{Q}_{tol}|\) and \(|Q_{z}| \ge |\mathbf{Q}_{tol}|\), then the Coulomb potential is described by,

(39)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2} }{2 V_{uc} \epsilon_{0} \mathbf{Q}^2} \left[ 1 - \cos(0.5 |Q_{z}| \mathbf{r}_{3}) - (0.5 |Q_{z}| \mathbf{r}_{3}) \sin(0.5 |Q_{z}| \mathbf{r}_{3}) \right] . \end{equation}\]

Else if \(\sqrt{Q_{x}^{2}+Q_{y}^{2}} < |\mathbf{Q}_{tol}|\) and \(|Q_{z}| < |\mathbf{Q}_{tol}|\), then Coulomb potential becomes,

(39)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2} }{2 V_{uc} \epsilon_{0}} \left(\frac{\mathbf{r}_{3}^{2}}{8} \right)~, \end{equation}\]

2D truncated simplified Coulomb Potential (V2DT2)

A variation for the truncated Coulomb potential for 2D systems proposed by IsmailBeigi textit{et. al.} [IsmailBeigi2006] being selected by the flag \(\texttt{COULOMB\_POT= V2DT2}\) in main input file. It has the following expression:

(41)\[\begin{align} &V(\mathbf{Q}) = -\frac{e^{2}}{2 V_{uc} \epsilon_{0} \mathbf{Q}^2} \left[1-\text{e}^{-Q_{xy}Z_{c}}\cos(Q_{z}Z_{c})\right] ~, \end{align}\]

where \(Q_{xy}=\sqrt{(Q_{x}^{2}+Q_{y}^{2})}\), \(V_{uc}\) is the unit cell volume and \(Z_c\) is half of the lattice vector in \(\hat{z}\) direction, given by the flag \(\texttt{LC}\).

To avoid the divergence at \(\mathbf{Q} \rightarrow 0\), we suppose

(41)\[\begin{equation} V(\mathbf{|Q|<Q_{tol}}) = 0~. \end{equation}\]

2D Ohno Potential (V2DOH)

Ohno Coulomb potential for 2D systems [Trolle2014] being selected by the flag \(\texttt{COULOMB\_POT= V2DOH}\) in main input file. It is given by the following expression:

(43)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2}}{2 A_{uc} \epsilon_{0} \epsilon_{d} |\mathbf{Q}|} e^{-w_{0} |\mathbf{Q}|}~, \end{equation}\]

where \(A_{uc}\) is the unit cell area, \(\epsilon_{d}\) is an effective dielectric constant parameter, given by the flag \(\texttt{EDIEL}\), \(w_{0}\) is a free parameter given by the flag \(\texttt{W\_COUL}\).

To avoid the divergence at \(\mathbf{Q} \rightarrow 0\), we define

(43)\[\begin{equation} V(\mathbf{|Q|<Q_{tol}}) = 0~. \end{equation}\]

0D Potential (V0DT)

An adaption of the Coulomb potential for 0D systems, proposed by Rozzi \(\textit{et. al.}\) [Rozzi2006], being selected by the flag \(\texttt{COULOMB\_POT= V0DT}\) in main input file. It is described by the following expression:

(44)\[\begin{equation} V(\mathbf{Q}) = -\frac{e^{2}}{2 V_{uc} \epsilon_{0} \mathbf{Q}^{2}} \left[1- \cos(0.5 \mathbf{r}_{min} |\mathbf{Q}|) \right] \end{equation}\]

To avoid the divergence at \(\mathbf{Q} \rightarrow 0\), we consider the following condition

where \(V_{uc}\) is the unit cell volume and \(\mathbf{r}_{min}\) is the size of the smallest lattice vector (real space).

Exciton radiative lifetime

In this code the exciton radiative lifetime is calculated directly from the excitonic state oscillator force (\(F_{\alpha,\beta}^{n,BSE}\)), described in equation (25). Being calculated using the flag \(\texttt{BSE= T}\) combined with \(\texttt{SPEC= T}\), where the lifetime is obtained in seconds. This quantity is only calculated for direct excitons.

For 3D systems (\(\texttt{SYSDIM= 3D}\)), the lifetime, for the n-th excitonic states, is obtained by,

(45)\[\begin{equation} \tau_{3D,\alpha,\beta}^{n} = \frac{3 c^{2} \hbar^{3} N_{\mathbf{k}} }{4 \chi \left(E_{0}^{n}\right)^{3} F_{\alpha,\beta}^{n,BSE}} ~. \end{equation}\]

While for 2D systems (\(\texttt{SYSDIM= 2D}\)), a similar expression is given by [Palummo2015]

(46)\[\begin{equation} \tau_{2D,\alpha,\beta}^{n} = \frac{\hbar A_{uc} N_{\mathbf{k}} }{8 \pi \chi E_{0}^{n} F_{\alpha,\beta}^{n,BSE}}~. \end{equation}\]

Finally for 1D systems (\(\texttt{SYSDIM= 1D}\)), it is described by [Spataru2005]

(47)\[\begin{equation} \tau_{1D,\alpha,\beta}^{n} = \frac{c \hbar^{2} \mathbf{r}_{1} N_{\mathbf{k}} }{2 \pi \chi \left(E_{0}^{n}\right)^{2} F_{\alpha,\beta}^{n,BSE}} ~, \end{equation}\]

where \(c\) is the light speed, \(E_{0}^{n}\) is the n-th excitonic state energy, \(A_{uc}\) is the 2D crystal unit cell area, \(\mathbf{r}_{1}\) is the 1D unit cell size and \(\chi\) is the fine structure constant, given by:

(48)\[\begin{equation} \chi = \frac{e^{2}}{4 \pi \epsilon_{0} \hbar c} \approx \frac{1}{137}~. \end{equation}\]

k-mesh generator

The \(\textbf{k}\)-meshes are generated using the Monkhorst–Pack scheme [Monkhorst1976], where the \(\textbf{k}\)-points are defined by the following expression \(\texttt{NGY}\) and \(\texttt{NGZ}\) in the input file, it’s possible to define it through the flags \(\texttt{MESH\_TYPE= RK3D}\) or \(\texttt{MESH\_TYPE= RK2D}\), choosing a density of \(\textbf{k}\)-points defined by the flag \(\texttt{RK}\). The values of \(\texttt{NGX}\), \(\texttt{NGY}\) and \(\texttt{NGZ}\), are obtained for the option \(\texttt{MESH\_TYPE= RK3D}\) by the following expression:

(49)\[\begin{split}\begin{eqnarray} NGX = int \left(max\left(1,\texttt{RK} \frac{|\mathbf{b}_{1}|}{2 \pi}+0.5\right)\right)\nonumber \\ NGY = int \left(max\left(1,\texttt{RK} \frac{|\mathbf{b}_{2}|}{2 \pi}+0.5\right)\right) \nonumber \\ NGZ = int \left(max\left(1,\texttt{RK} \frac{|\mathbf{b}_{3}|}{2 \pi}+0.5\right)\right) \end{eqnarray}\end{split}\]

for the option \(\texttt{MESH\_TYPE= RK2D}\) the only difference is that \(\texttt{NGZ=1}\) independent of \(\texttt{RK}\).

Solar Cell Power Conversion Efficiency

To estimate the power conversion efficiency (PCE), of the studied system, in the scope of Shockley-Queisser Limit (SQ-limit), [Shockley1961] and spectroscopy limited maximum efficiency (SLME) scheme [Yu2012]. This simulation is available by the input flag \(\texttt{PCE= T}\).

Both approximations are based on the principle of detailed balance between absorbed photons and emitted photons. Below, we will summarize the mathematical formalism [Bercx2018] and the SQ-limit and SLME approximations.

The PCE is defined as the maximum output power density (\(P_{\text{PV}}\)) from the photovoltaic device divided by the total incoming power density from the solar spectrum, \(P_{solar}\), as follows,

(50)\[\begin{equation} \text{PCE} = \frac{P_{\text{PV}}}{P_{solar}}~, \end{equation}\]

where

(51)\[\begin{equation} P_{solar} = \int_{0}^{\infty}P(E)dE~, \end{equation}\]

and \(P(E)\) is the solar energy flux, which is selected by the flag \(\texttt{SES}\), and has the options: \(\texttt{AM15G}\), \(\texttt{AM15D}\) and \(\texttt{AM0G}\), which corresponds to the standard solar spectrum for non-concentrated photovoltaic conversion, taking light absorption and scattering in the atmosphere into account [Astm2012]. The output power density is described by the product \(J(V)V\), as the maximum output power density (\(P_{\text{PV}}\)) is obtained maximizing the \(J-V\) characteristic of an illuminated solar cell :

(52)\[\begin{equation} P_{\text{PV}} = J(V_{max}) V_{max}~, \end{equation}\]

where \(V_{max}\) is the voltage that results in the maximum output power density. The open circuit voltage \(V_{oc}\) is the voltage that minimizes \(J(V)\). In this method, the current density, \(J(V)\), are described by the following expression:

(53)\[\begin{equation} J(V) = J_{sc}-\frac{J_{0}}{fr}\left(exp\left(\frac{e V}{k_{\text{B}}T}\right)-1\right)~, \end{equation}\]

where \(k_{\text{B}}\) is Boltzmann’s constant, \(e\) is the elementary charge, \(fr\) is the radiative electron-hole recombination fraction and \(T\) is the temperature of the solar cell, defined by the flag \(\texttt{CTEMP}\). \(J_{sc}\) is the short-circuit current density, also known as the illuminated current or photogenerated current, calculated from the following expression,

\[\]

begin{equation} J_{sc} = eint_{0}^{infty}a(E)frac{P(E)}{E}dE~. end{equation} \(a(E)\) is the absorbance, which is defined as the ratio of power absorbed by the solar device over the power of incident sunlight.

\(J_{0}\) is the reverse saturation current density, calculated considering the detailed balance principle when in equilibrium conditions, the photon emission rate from radiative recombination is equal to the photon absorption from the surrounding environment. This was done considering our solar cell attached to an ideal heat sink, expecting the cell temperature to be the ambient temperature. Hence, the spectrum of the surrounding environment is that of a black body at cell temperature \(T\),

(54)\[\begin{equation} J_{0} = e \pi \int_{0}^{\infty} a(E) \Phi_{bb}(E) dE~, \end{equation}\]

where,

(55)\[\begin{equation} \Phi_{bb}(E) = \frac{2E^{2}}{h^{3}v_{c}^{2}} \left(e^{\frac{E}{k_{\text{B}}T}}-1\right)^{-1}~, \end{equation}\]

and \(h\) is Planck’s constant and \(v_{c}\) is the speed of light. Complementary we can also calculate the Fill-Factor (\(FF\)) of our solar device, by the expression:

(56)\[\begin{equation} FF = \frac{V_{max} J(V_{max})}{V_{oc} J_{sc}}~. \end{equation}\]

Shockley–Queisser Limit

In the SQ-limit approximation, we consider the absorbance, \(a(E)\), as a Heaviside step function, where all photons with energy higher or equal to the energy bandgap, \(E_g\), are absorbed, i.e., \(a(E) = 1\) for \(E{~}\ge{~}E_{g}\), and \(a(E) = 0\) for \(E{~}<{~}E_{g}\). Furthermore, this approximation assumes \(fr = 1\), i.e., the radiative recombination process is the only recombination process for all systems, neglecting the non-radiative recombinations (e.g., Auger recombination) due to indirect bandgap [Huldt1971], [Green_671_1984].

Spectroscopy Limited Maximum Efficiency

In contrast with the SQ-limit approach, the SLME approximation requires the total absorption coefficient,

(57)\[\begin{equation} \alpha(E) = \frac{4\pi}{\lambda}k(E)~, \end{equation}\]

where \(k(E)\) is the imaginary part of the dielectric function, \(\epsilon(E)\), and \(\lambda\)

(58)\[\begin{equation} A(\omega)=A_{x,x}(\omega)+A_{y,y}(\omega)+A_{z,z}(\omega)~, \end{equation}\]

material thickness, \(\Delta\), and the bandgaps. Its also considers the non-radiative recombination in the solar cell by modeling the fraction of radiative recombination, \(fr\), using a Boltzmann factor [Yu2012],

(59)\[\begin{equation} fr = e^{-\frac{\delta}{k_{B}T}} \end{equation}\]

for calculations in the scope of IPA (without excitons), we can define \(\delta = E_{op}-E_{g}\), where \(E_{op}\) is the optical bandgap, given by the flag \(\texttt{EGD}\) and \(E_{g}\) is the fundamental energy bandgap, given by the flag \(\texttt{EG}\). In the scope of BSE (with excitonic effects), we can define \(\delta = Ex_{br}-Ex_{gs}\), where \(Ex_{br}\) is the exciton bright ground state (optical bandgap considering excitonic effects), given by the flag \(\texttt{EBGS}\), while \(Ex_{gs}\) is the exciton ground state, given by the flag texttt{EGS}. The absorbance, \(a(E)\), is obtained considering the solar cell with the same conditions of SQ-Limit described previously, given by the following expression [Yu2012], [Duan2019].

(60)\[\begin{equation} a(E) = 1 - e^{-2 \ A(\omega) \ \Delta}~, \end{equation}\]

where \(E = \hbar\omega\), and \(\Delta\) is the material thickness, in this code we calculate the PCE, in the scope of SLME method from \(\Delta=0\) until \(\Delta= \Delta_{max}\), being \(\Delta_{max}\) defined by the flag \(\texttt{THMAX}\).

References