Excitonic effects in the optical properties of 2D materials: An equation of motion approach

TL;DR

This paper introduces an efficient equation of motion method to analyze excitonic effects in 2D transition-metal dichalcogenides, accurately capturing binding energies and spectral features with less computational effort.

Contribution

It presents a unified, less computationally demanding approach to study excitonic properties in TMDC monolayers, incorporating exchange energy effects and providing semi-analytical formulas.

Findings

Accurately models exciton binding energies and spectra in TMDCs.

Shows the importance of exchange energy at the $b3$-point for spectral accuracy.

Achieves good agreement with experimental Rydberg series and absorption spectra.

Abstract

We present a unified description of the excitonic properties of four monolayer transition-metal dichalcogenides (TMDC's) using an equation of motion method for deriving the Bethe-Salpeter equation in momentum space. Our method is able to cope with both continuous and tight-binding Hamiltonians, and is less computational demanding than the traditional first-principles approach. We show that the role of the exchange energy is essential to obtain a good description of the binding energy of the excitons. The exchange energy at the $\Gamma-$point is also essential to obtain the correct position of the C-exciton peak. Using our model we obtain a good agreement between the Rydberg series measured for WS$_2$. We discuss how the absorption and the Rydberg series depend on the doping. Choosing $r_0$ and the doping we obtain a good qualitative agreement between the experimental absorption and our calculations for WS$_2$. We also derive a semi-analytical version of Ellitot's formula for TMDC's.

Tables

Table 1. Table 1: The first five columns give the material parameters used in all calculations in this paper. The parameters in the first four columns come from Ref. [ 27 ] and in the fifth column from Ref. [ 12 ] . The sixth column is the direct gap at the K -point calculated with the exchange self-energy obtained from our model. The last columns are experimental data (numbers in square braces refer to references). See also [ 41 ] for a different set of parameters and, in particular, the prediction of SOC induced splitting in the conduction band at the K -point.

TMDC	$\Delta$(eV)	$\hslash v_F$(eV/Å)	${\Lambda }_1$(eV)	${\Lambda }_2$(eV)	$r_0$(Å)	${\Delta }_K$(eV)		Experimental gap
MoS2	0.797	2.76	0.076	0.073	31.4	2.82	2.5 [22]	2.14 [34]	1.86 [35]
MoSe2	0.648	2.53	0.104	0.082	51.7	2.37	2.18 [36]	2.02, 2.22 [37]	1.58 [38]
WS2	0.685	3.34	0.164	0.230	37.9	2.78	2.14 [35]	2.41 [39]
WSe2	0.524	3.17	0.215	0.252	45.1	2.31	2.51 [40]	2.0, 2.18 [37]

Table 2. Table 2: Binding energy (in eV) for four different excitonic states of MoS 2 corresponding to the A-series. The first six data rows show the results for the numerical procedure based on a constant grid discretization; convergence is obtained only for very large grids of the order of 3000. GL accounts for Gauss-Legendre/Laguerre, where the integral in q 𝑞 q is divided in three intervals [ 0 , 0.3 ​ Δ ] 0 0.3 Δ [0,0.3\Delta] , [ 0.3 ​ Δ , 0.6 ​ Δ ] 0.3 Δ 0.6 Δ [0.3\Delta,0.6\Delta] , [ 0.6 ​ Δ , ∞ ] 0.6 Δ [0.6\Delta,\infty] . For the first interval we use 100 points, and for the second and third 50 points are used. In the first two intervals we use Gauss-Legendre quadrature, and the last one Gauss-Laguerre with a rescale of m / 150 𝑚 150 m/150 . The last row (LM) is the data from reference [ 18 ] , with excitonic binding energy computed using a lattice (tight-binding) model. We want to stress that the discretization procedure in a linear mesh of 3000 points takes several hours to run in a conventional laptop, whereas the GL method takes only few minutes in the same computer. At a given stage the calculation requires performing an angular integral of the Keldysh potential V ​ ( 𝐤 − 𝐪 ) 𝑉 𝐤 𝐪 V(\mathbf{k}-\mathbf{q}) . Although the integral can be proven finite, numerically the integral is ill behaved when k = q 𝑘 𝑞 k=q . For avoiding this pathology a finite constant of the order of the grid spacing is added to p k , q ​ ( θ ) subscript 𝑝 𝑘 𝑞 𝜃 p_{k,q}(\theta) in Eq. ( 97 ); this renders the integral finite, as it should.

	50	100	300	600	1200	3000	GL	LM Ref.[18]
1s	0.224	0.264	0.304	0.318	0.327	0.333	0.331	0.301
2s	0.033	0.055	0.081	0.092	0.099	0.104	0.103	0.099
2p+	0.051	0.077	0.107	0.119	0.126	0.132	0.132	0.125
2p-	0.062	0.090	0.122	0.134	0.142	0.148	0.147	0.150

Table 3. Table 3: Absorbance (in percentage) of A and B peaks for MoS 2 . The table gives a comparison between our theoretical model and the experimental results for samples deposited on silica ( [ 5 ] , [ 11 ] , [ 10 ] ), and a FET device [ 48 ] , where the MoS 2 is deposited on top of silicon and under a voltage gate of − 10 10 -10 V. The “Model” refers to the theoretical approach developed in this paper and we have considered MoS 2 on top of silica ( ϵ silica = 1.46 subscript italic-ϵ silica 1.46 \epsilon_{\rm{silica}}=1.46 ). This value of ϵ silica subscript italic-ϵ silica \epsilon_{\rm{silica}} translates into an ϵ m = ( 1 + 2.13 ) / 2 subscript italic-ϵ 𝑚 1 2.13 2 \epsilon_{m}=(1+2.13)/2 , which is the value we use in our equations. Remember that in our model both the peak intensity and the peak width are dependent on the choice of the relaxation rate γ K subscript 𝛾 𝐾 \gamma_{K} . We can artificially reduce the height of the peak at the expenses of increasing its width. It is worth noting the variation of the experimental values for the absorbance among themselves.

excitonic peak	Model	[5]	[11]	[10]	[48]
A	14	10.8	7.4	3.8	7.5
B	15	10.5	8.6	5.0	8.0

Equations

$$H_0^{s\tau }(\mathbf{k})=\hslash v_F\left(\tau {\sigma }_1{k}_x+{\sigma }_2{k}_y\right)+{\sigma }_3{m}_{s\tau }{v}_F^2-{\mu }_{s\tau }I,$$

$$m_{s\tau }=\Delta -\frac{s\tau }{2}\frac{{\Lambda }_1}{v_F^2},$$

$${\mu }_{s\tau }=\frac{1}{2}s\tau {\Lambda }_2.$$

$$V(q)=-\frac{e}{2{\epsilon }_0}\frac{1}{q(r_{0}q+{\epsilon }_m)},$$

$${\hat{H}}_I(t)=-e\mathcal{E}(t)\hat{x},$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_I(t)+{\hat{H}}_{\mathrm{ee}},$$

$${\hat{H}}_{\mathrm{ee}}=-\frac{e}{2}\int d{\mathbf{r}}_{1}d{\mathbf{r}}_2{\hat{\psi }}^{†}({\mathbf{r}}_1){\hat{\psi }}^{†}({\mathbf{r}}_2)V({\mathbf{r}}_1-{\mathbf{r}}_2)\hat{\psi }({\mathbf{r}}_2)\hat{\psi }({\mathbf{r}}_1),$$

$$\hat{\psi }(\mathbf{r},t)=\frac{1}{L}\sum _{\mathbf{k},\lambda ,s,\tau }{\varphi }_{\lambda }^{s\tau }(\mathbf{k}){\hat{a}}_{\mathbf{k}\lambda s\tau }(t){\text{rme}}^{-\mathrm{i}\mathbf{k}\cdot \mathbf{r}},$$

$$\phi^{s\tau}_{\lambda}(\mathbf{k})=\sqrt{\frac{E^{s\tau}_{k}+\lambda m_{s\tau}}{2E^{s\tau}_{k}}}\left(\begin{array}[]{cc}1\\ \frac{\tau k_{x}-\mathrm{i}k_{y}}{\lambda E^{s\tau}_{k}+m_{s\tau}}\end{array}\right),$$

$$E_k^{s\tau }=\sqrt{k^2+m_{s\tau }^2}.$$

$$\hat{P}(t)=\int d\mathbf{r}{\hat{\psi }}^{†}(\mathbf{r},t)(-ex)\hat{\psi }(\mathbf{r},t),$$

$$\int d\mathbf{r}{\varphi }_{{\mathbf{k}}^{\prime }{\lambda }^{\prime }{s}^{\prime }{\tau }^{\prime }}^{†}(\mathbf{r})$$

$$⟨s^{\prime }{\tau }^{\prime }{\mathbf{k}}^{\prime },{\lambda }^{\prime }\left|x\right|s\tau \mathbf{k},\lambda ⟩={\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}{\delta }_{s{s}^{\prime }}{\delta }_{\tau {\tau }^{\prime }}\frac{\mathrm{i}{v}_{{\lambda }^{\prime }}^{s\tau }(\mathbf{k})}{2{\lambda }^{\prime }{E}_k^{s\tau }},$$

$$v_{\lambda }^{s\tau }(\mathbf{k})=⟨s\tau \mathbf{k},\lambda \left|{\sigma }_1\right|s\tau \mathbf{k},-\lambda ⟩.$$

$$\hat{P}(t)=-\frac{\mathrm{i}e}{2}\sum _{s\tau \mathbf{k}\lambda }\frac{v_{\lambda }^{s\tau }(\mathbf{k})}{\lambda E_k^{s\tau }}{\hat{\rho }}_{\mathbf{k}\lambda ,-\lambda }^{s\tau }(t),$$

$$p_{\lambda }^{s\tau }(\mathbf{k},t)=⟨{\hat{\rho }}_{\mathbf{k}\lambda ,-\lambda }^{s\tau }(t)⟩.$$

$$n_{\lambda }^{s\tau }(\mathbf{k},t)=⟨{\hat{\rho }}_{\mathbf{k}\lambda \lambda }^{s\tau }(t)⟩.$$

$$-\mathrm{i}\frac{d}{dt}{\hat{\rho }}_{\lambda {\lambda }^{\prime }}^{s\tau }(t)=[\hat{H}$$

$$-\mathrm{i}{\partial }_t{p}_{\lambda }^{s\tau }(\mathbf{k},t)=$$

$$-{\partial }_t{n}_{\lambda }^{s\tau }(\mathbf{k},t)=2\Im \left[\left({\tilde{\Omega }}_{\mathbf{k}-\lambda }^{s\tau }(t)+{\mathcal{D}}_{\mathbf{k}\lambda }^{s\tau }(t)\right)p_{\lambda }^{s\tau }\right],$$

$${\tilde{\Omega }}_{\mathbf{k}\lambda }^{s\tau }(\mathbf{k},t)=-\frac{\mathrm{i}\mathcal{E}(t)v_{-\lambda }^{s\tau }(\mathbf{k})}{2\lambda E_{\lambda k}^{s\tau }}+{\mathcal{B}}_{\mathbf{k}\lambda }^{s\tau }(\mathbf{k},t),$$

$$\displaystyle{\cal B}^{s\tau}_{\mathbf{k}\lambda}(t)=\frac{1}{S}\sum_{\mathbf{q}}V(|\mathbf{k}-\mathbf{q}|)\Big{[}p^{s\tau}_{\lambda}(\mathbf{q},t)F^{s\tau}_{\lambda^{\prime}\lambda^{\prime}\lambda\lambda}(\mathbf{k},\mathbf{q})+$$

$$\displaystyle+p^{s\tau}_{\lambda^{\prime}}(\mathbf{q},t)F^{s\tau}_{\lambda^{\prime}\lambda\lambda^{\prime}\lambda}(\mathbf{k},\mathbf{q})\Big{]},$$

$${\tilde{\omega }}_{\lambda \mathbf{k}}^{s\tau }=2\lambda E_k^{s\tau }+\lambda {\Sigma }_{\mathbf{k},\lambda }^{s\tau ,\mathrm{xc}},$$

$${\Sigma }_{\mathbf{k},\lambda }^{s\tau ,\mathrm{xc}}(t)=\frac{\lambda }{S}\sum _{\mathbf{q}}V(q)\Delta n_{\lambda }^{s\tau }(\mathbf{k}-\mathbf{q},t)\times$$

$$\displaystyle\times\big{[}F^{s\tau}_{\lambda^{\prime}\lambda\lambda\lambda^{\prime}}(\mathbf{k},\mathbf{k}-\mathbf{q})-F^{s\tau}_{\lambda\lambda\lambda\lambda}(\mathbf{k},\mathbf{k}-\mathbf{q})\big{]},$$

$${\mathcal{W}}_{\mathbf{k}\lambda }^{s\tau }(t)$$

$${\mathcal{D}}_{\mathbf{k}\lambda }^{s\tau }(t)=\frac{1}{S}\sum _{\mathbf{q}{\lambda }_1}V(q)$$

$$F_{{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4}^{s\tau }({\mathbf{k}}_1,$$

$$n_{\lambda }^{s\tau }(\mathbf{k},t)=f_{\lambda }^{s\tau }$$

$${\Sigma }^{s\tau ,\mathrm{xc}}(\mathbf{k})=-\int \frac{d\mathbf{q}}{4{\pi }^2}V(q)\Delta f_{\mathbf{k}-\mathbf{q}}^{s\tau }\frac{\mathbf{k}\cdot \mathbf{q}+m_{s\tau }^2}{E_k^{s\tau }{E}_q^{s\tau }},$$

Full text

1 Introduction

The study of excitons in bulk transition-metal dichalcogenides (TMDC’s) is a research topic in condensed matter physics that dates back to 1960’s [1, 2]. With the advent of two-dimensional materials [3], this topic regained interest since it became possible to study single- and few-layers of TMDC’s [4, 5]. Together with its two-dimensional nature, this new class of materials also has a hexagonal lattice structure as does graphene. On the other hand, while the low-energy electronic excitations are in graphene described by a massless Dirac equation, in monolayer TMDC’s the same excitations can be described by a massive (with a gap) Dirac equation. The absence of a gap in graphene prohibits the existence of bound-states of excitons (but not of excitonic resonances [6]). On the contrary, we find in TMDC’s absorption spectrum fingerprints of both excitonic bound states (including the presence of a Rydberg series) and of excitonic resonances, due to electron-hole scattering processes, with energies above the non-interacting gap.

As a consequence of optical experimental studies in few-layers TMDC’s, the study of a new type of excitons in these novel 2D materials became possible. This has attracted a wealth of scientific research [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23]. The signature of excitons appeared first in the optical measurements of monolayer MoS2 [5], where two peaks in the absorbance, with energies $\sim 1.9$ eV and $\sim 2.1$ eV, were identified. These two peaks correspond approximately to the same results found in several layers of MoS2 [1].

The optical studies of other monolayers of TMDC’s soon followed at the pace of their synthesis. The optical properties of MX2, M={Mo,W}, X={S,Se}, in the range $1.5-3$ eV were experimentally studied by Li et al. [11], with reflectance and transmittance measurements followed by a Kramers-Kronig analyses, and by Morozov and Kuno [10], with differential transmission and reflectance measurements. It should be noted that all these four materials have similar optical properties. Their optical absorbance spectra show signatures of the spin-orbit splitting for excitons at the K(K’) points in the Brillouin zone, as well as signatures of excitonic resonances at the $\Gamma$-point.

The properties of the excitons at the K-point were extensively studied in the framework of tight-binding and Bethe-Salpeter equation (BSE) [17, 18], DFT+GW+BSE [15, 16], and gapped 2D Dirac-equation [18]. One of the most prominent features in the optical spectra of these materials is its dependence on the Berry phase, which generates a modified Rydberg series [19, 20]. The form of the electron-electron interaction potential, which deviates form the Coulomb one, also contributes to a modified Rydberg series [21]. The C-excitonic resonance in monolayer MoS2, due to transitions at the $\Gamma$-point, was first calculated by Qiu et al. [16], and it was associated with a minimum in the optical band structure around the $\Gamma$-point by Klots et al. [22]. The effects of temperature and carrier density in MoS2 were studied either solving the semiconductor Bloch equation (SBE) with a tight-binding Hamiltonian, whose parameters were obtained from a $G_{0}W$ calculation [23], or by combining a LDA+BSE approach with the inclusion of electron-phonon coupling [7].

In the present work we use the polarization concept formalism [24, 6] for describing the excitonic properties of monolayer TMDC’s. This formalism is easily applied to any system, both using low-energy effective models or tight-binding ones. The development of the formalism boils down to the solution of an eigenvalue problem for determining the excitonic bound states and to the solution of a linear system of equations for computing the optical conductivity of the system. We apply the resulting equations to a two-band gapped Dirac equation for describing the physics around the K-point; this originates the physics of the A and B excitons in the TMDC’s and of a modified Rydberg series. On the other hand, using as a starting point the three-band model for TMDC’s [25] we describe the formation of an exitonic resonance near the $\Gamma$-point. This approach allows us to make much analytical progress and clearly identify the origin of different bound-states and resonances in the absorption spectrum. In this regards, our approach is distinct from previous ones that consider the full band structure as a starting point. The advantage of our approach lies in the possibility of clearly identify the origin of the different peaks in the optical conductivity, or absorbance for the same matters, of TMDC’s.

We show that the optical properties have a strong dependence on external parameters, namely, temperature and dielectric function of the environment. This dependence on external parameters opens the possibility of engineering at will nano-materials showing strong optical response in the spectral range from the IR to visible. The application of these systems to opto-electronics, including photo-detectors, will launch a new set of devices in this area. Another possibility that these 2D materials may provide is the engineering Bose-Einstein condensation of excitons when a TMDC is put inside an optical cavity [26].

The paper is organized as follows: in section 2 we introduce the second quantized form of the Hamiltonian, which is composed of three pieces: the non-interacting part, the light-matter interaction term, and the Coulomb interaction. Using this Hamiltonian the excitonic properties at the Dirac point are worked out. In section 3 we derive the optical properties of four TMDC’s around the K-point in the Brillouin zone. Using a three-band tight-binding model we describe the excitonic properties of four TMDC’s in section 4. The excitonic effects around the $\Gamma-$point in the Brillouin zone are actually resonances, as they are above the continuum. In section 5 the optical properties of four TMDC’s are discussed in detail and compared with the existent experimental data for the absorption. We find a good agreement with the experimental data, although, since we do not include electron-phonon interaction, the agreement is not quantitative. We also note that the experimental values for the absorption present discrepancies among different experiments. This led us to conclude that there is a clear sample-dependence in the absorption measurements. In two of the TMDC’s we stress the absence of the B-excitonic series in the experimental data, which is a noticeable discrepancy with our theoretical calculations. This led us to believe that the experiments need to be repeated for encapsulated TMDC’s in h-BN at low temperatures. This approach screens away the effect of extrinsic disorder, and reduces the impact of phonons in the absorption spectrum, due to low temperatures. Finally, in section 6 we provide a summary of the main conclusions of the paper. A set of appendices give details of the calculations.

2 Formalism and $\mathbf{K}$-point excitons

In this section we introduce the effective model for electronic properties of TMDC’s around the $\mathbf{K}$-point in the Brillouin zone. Since we are dealing with a many-body problem, the second quantization formalism is used throughout the paper. Using the full interacting Hamiltonian, the equations of motion for the density matrix are obtained and from it the total polarization is derived. The electron-electron interaction generates a hierarchy of correlation functions that are truncated at the random-phase approximation (RPA) level. This procedure is equivalent, in a diagrammatic approach, to the inclusion diagrams considering only the interaction in the electron-hole propagator. The diagrams relevant to our calculation are given in figure 1, where $K$ represents the BSE kernel (see ahead).

Although the diagrammatic approach is a possible route to solve the problem of excitonic effects in TMDC’s, it is also possible to address it using an equation-of-motion approach. The latter formalism enables treating at the same level of approximation the exchange-energy correction and the excitonic effects, and that is the path we will follow in this paper.

2.1 Many-body Hamiltonian

The low-energy single-particle electronic-excitations of TMDC materials can be described, in the $\mathbf{k}\cdot\mathbf{p}$ approximation, by a 2D gapped Dirac equation. When we consider spin-orbit coupling (SOC), the effective mass and chemical potential become valley and spin dependent. With these aspects in mind, we can write the single-particle Hamiltonian for a single combination of valley($\tau$)/spin($s$) index as:

[TABLE]

with $\sigma_{i}$ the usual Pauli matrices and $I$ the identity matrix, $\tau=\pm$ the valley index, and $s=\pm$ the spin index. The effective mass, $m_{s\tau}$, and the on-site energy, $\mu_{s\tau}$, can be written in terms of the SOC parameters, $\Lambda_{1}$ and $\Lambda_{2}$ [27], and of the mass $\Delta$ as:

[TABLE]

[TABLE]

with $\Lambda_{1}=\Delta_{\rm{VB}}-\Delta_{\rm{CB}}$ and $\Lambda_{2}=\Delta_{\rm{VB}}+\Delta_{\rm{CB}}$, with $\Delta_{\rm{VB}}$ ($\Delta_{\rm{CB}}$) the spin-splitting of the valence (conduction) band.

The band structure implied by Hamiltonian (1) around the $\mathbf{K}-$point in the Brillouin zone is depicted in figure 2 (different colors correspond to opposite spin projections).

In three dimensions (3D), the electron-electron interaction in a dielectric medium is given by the Coulomb potential in vacuum but with the permittivity of free space $\varepsilon_{0}$ replaced by the medium permittivity $\varepsilon_{m}\varepsilon_{0}$. Contrary to 3D, the same procedure does not hold in 2D materials. In contrast, the electron-electron interaction is described by the Keldysh potential [28, 29]. This takes into account the surface charge polarization from a dielectric thin film and reads in momentum space:

[TABLE]

where $q$ is the 2D transfered momentum, $\varepsilon_{m}$ and $r_{0}$ are the capping dielectric function of the environment and a material-dependent constant, respectively, the latter measuring the deviation from the 2D Coulomb potential. Note that we recover the 2D Coulomb potential making $r_{0}\rightarrow 0$. The potential (4) is written in a slightly different manner than in reference [29] for removing the dependence of the parameter $r_{0}$ on the external dielectric constant.

To calculate the optical properties of TMDC’s, we consider the interaction of the electron gas with a time-dependent electric field ${\cal E}(t)$ polarized along the $x$ axis. For describing the light-matter interaction in this problem we use the dipole-coupling Hamiltonian

[TABLE]

with $\hat{x}$ the position operator.

From now on, we consider the full many-body Hamiltonian as:

[TABLE]

where $\hat{H}_{0}$ is built from the single-particle Hamiltonian (1) and $\hat{H}_{\mathrm{ee}}$ is the electron-electron interaction:

[TABLE]

where $\hat{\psi}(\mathbf{r})$ is the field operator (8) defined below and $V(\mathbf{q})$ is the Fourier transform of the Keldysh potential given by equation (4). For simplicity, from here on we choose units such that $v_{F}=\hbar=e=1$; the usual units are reintroduced at the end of the calculations. The field operator is given by:

[TABLE]

with $L$ the square-box side length, $\hat{a}_{\mathbf{k}\lambda s\tau}$ the usual annihilation operator that obeys anticommutation relations and $\phi_{\mathbf{k}\lambda s\tau}$ the eigenfunctions of $H_{0}^{\mathrm{s\tau}}$, $H_{0}^{\mathrm{s\tau}}\phi^{s\lambda}_{\tau}(\mathbf{k})=(\lambda E_{k}^{s\tau}+\mu_{s\tau})\phi^{s\lambda}_{\tau}(\mathbf{k}),$ with $\lambda=-$ ($+$) for the valence (conduction) band. The eigenfunctions and positive eigenvalues are given by:

[TABLE]

[TABLE]

The eigenfunctions (9) will be used extensively in this work for determining the four-body structure factor.

2.2 Polarization operator

For obtaining the optical conductivity and the absorbance, we have to compute the expectation value of the polarization operator. As noted above, we consider an electric field along the $x-$axis direction and define the polarization operator along the same spatial orientation. Using the field operators (8) the polarization operator reads:

[TABLE]

The integral can be explicitly computed with the help of the eigenfunctions of $H_{0}$ in position space: $\phi^{s\tau}_{\mathbf{k}\lambda}(\mathbf{r})=\phi^{s\tau}_{\lambda}(\mathbf{k})\rme^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}$. Using these eigenfunctions it follows that

[TABLE]

Noting that $[x,H_{0}(x)]=-\mathrm{i}\sigma_{1}$ we obtain for the dipole matrix element the result:

[TABLE]

for $\lambda\neq\lambda^{\prime}$ (inter-band transitions). We defined the matrix element of the velocity operator $\sigma_{1}$ as

[TABLE]

Finally, we can express the polarization operator as:

[TABLE]

where we have introduced the following matrix in the Heisenberg picture: $\hat{\rho}^{s\tau}_{\mathbf{k}\lambda\lambda^{\prime}}(t)=\hat{a}^{\dagger}_{\mathbf{k},\lambda s\tau}(t)\hat{a}_{\mathbf{k},\lambda^{\prime}s\tau}(t)$. This matrix is written in the basis that diagonalizes $\hat{H}_{0}$ and is equivalent to the density matrix for the states of $\hat{H}_{0}$. To determine the expectation value of the density matrix, and therefore the polarization, we use the Heisenberg’s equation-of-motion for $\hat{\rho}^{s\tau}_{\mathbf{k}\lambda\lambda^{\prime}}(t)$.

2.3 Equations of motion for the density matrix

As noted, for calculating the expectation value of the density matrix in the right-hand-side of equation (15), we use Heisenberg’s equation-of-motion method. We define the expectation value for the off-diagonal elements of the density matrix, corresponding to the transition probabilities $p_{\lambda}^{s\tau}(\mathbf{k},t)$, as:

[TABLE]

The diagonal elements of the matrix correspond to a new electronic distribution $n_{\lambda}^{s\tau}(\mathbf{k},t)$ defined as:

[TABLE]

Both the diagonal and off-diagonal matrix elements of the density matrix are time-dependent. Explicitly, Heisenberg’s equation-of-motion for the density matrix reads:

[TABLE]

where we split the full Hamiltonian into the three components introduced in equation (6). The commutators in right-hand-side of equation (18) are explicitly calculated in A. The resulting equations-of-motion for the expectation value of equation (18) are:

[TABLE]

[TABLE]

where the transition energy (also denoted bare optical band ahead) is $\omega_{\lambda\mathbf{k}}^{s\tau}=2\lambda E_{\mathbf{k}}^{s\tau}$, the difference in occupations reads $\Delta n^{s\tau}_{\lambda}(\mathbf{k},t)=n^{s\tau}_{\lambda}(\mathbf{k},t)-n^{s\tau}_{-\lambda}(\mathbf{k},t)$, and the renormalized Rabi frequency $\tilde{\Omega}^{s\tau}_{\mathbf{k}\lambda}(t)$ reads:

[TABLE]

and depends on the dipole moment, the electric field, and the transition energy. The different terms that appear in equations (19) and (22) are classified as:

• •

Excitonic Rabi frequency renormalization:

[TABLE]

and where in this expression we have the constraint $\lambda^{\prime}=-\lambda$, $S=L^{2}$ is the area of the system, and $F^{s\tau}_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}}(\mathbf{k_{1}},\mathbf{k_{2}})$ is defined in equation (27)

• •

Renormalized transition energy (or interacting optical band):

[TABLE]

where the exchange self-energy is:

[TABLE]

and where in this expression we have the constraint $\lambda^{\prime}=-\lambda$. This form of the exchange self-energy is the same we find in the jelium model, except for the non-trivial four-body structure factor $F^{s\tau}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}(\mathbf{k},\mathbf{k}-\mathbf{q})$.

• •

Non-linear contribution:

[TABLE]

where in this expression $\lambda^{\prime}=-\lambda$, $\lambda_{1}^{\prime}=-\lambda_{1}$.

• •

Density term:

[TABLE]

Note that above we have defined the four-body spinor product (or structure factor) as:

[TABLE]

that follows from writing the electron-electron interaction (7) in the basis that diagonalizes $H_{0}$ (9). Explicit equations for the above terms are given in A. In the next sections we neglect the density ${\cal D}^{s\tau}_{\mathbf{k}\lambda}(t)$ and the non-linear ${\cal W}^{s\tau}_{\mathbf{k}\lambda}(t)$ contributions, as we are essentially focused on the exchange self-energy and excitonic effects to linear order. We also consider that the system is in thermodynamic equilibrium, where the electronic distribution $n^{s\tau}_{\lambda}(\mathbf{k},t)$ is given by the Fermi-Dirac distribution function $f^{s\tau}_{\lambda}(\mathbf{k})$:

[TABLE]

The latter approximation is valid in the linear regime (weak external electric fields).

2.4 Calculation of the exchange energy

The exchange self-energy (24) reshapes the electronic bands and, as will be shown later, it is essential to correctly describe the optical properties of TMDC’s, specially the value of the independent-particle energy gap.

For graphene, described by a massless Dirac equation ($m^{s\tau}=0$), the exchange self-energy was calculated in [32] and [33] with the use of the Couloumb potential. Here we calculate the exchange self-energy using the Keldysh interaction for the gapped Dirac equation ($m\neq 0$).

The self-energy for the optical band structure, $\Sigma^{s\tau,\mathrm{xc}}(\mathbf{k})$, is calculated from equation (24) using the expressions in A. Its calculation boils down to an integral over all possible momentum values:

[TABLE]

where the difference between the electronic valence and conduction distribution functions is defined as $\Delta f^{s\tau}_{\mathbf{q}}=f^{s\tau}_{+}(\mathbf{q})-f^{s\tau}_{-}(\mathbf{q}).$

The direct gap renormalization, for each pair spin/valley, is given by the difference of the spin/valley top valence band and bottom conduction band energies, that is:

[TABLE]

and for $T=0$, can be calculated analytically from equation (29), resulting in

[TABLE]

with $\alpha\approx 1/137$ the fine structure constant, $\beta=v_{F}/c$ the ratio between Fermi-velocity and the speed of light, and the valley-spin dependent Fermi momentum $k_{F}^{s\tau}$ is given by:

[TABLE]

We depict the dependence of gap-renormalization in figure (3). The temperature dependence of the same quantity is given in figure (4), and the renormalized band in figure (5). We can see the strong dependence of the exchange energy on the external parameters (temperature and dielectric constant of the medium surrounding the TMDC), showing that the environment plays a key role on the optical properties of these materials.

A more accurate approach to the calculation of renormalization of the band gap due to electron-electron interactions requires a self-consistent approach, where the unperturbed Hamiltonian $H_{0}^{s\tau}$ (1) is defined including the self-energy from electron-electron interactions. For graphene this procedure is used to study the possibility of dynamical generation of a gap in the spectrum [30, 31].

2.5 Excitonic effects in the Rabi frequency

The electron-electron interaction induces the creation of electron-hole bound states below the non-interacting energy gap (corrected by the exchange self-energy), that corresponds to excitonic bound states. Also excitonic resonances appear above the energy gap. These two effects are routinely measured in optical experiments in semiconductors [24]. In the equation-of-motion description, the electron-electron interaction renormalizes the Rabi frequency, and, as will be shown later, this corresponds to solve the Bethe-Salpeter equation in the ladder approximation and in the center-of-mass reference frame. This procedure allows the calculation of the renormalized expectation value of the $\hat{x}$ operator.

The renormalized Rabi frequency (22), term ${\cal B}^{s\tau}_{\mathbf{k}\lambda}(t)$, can be split into two parts (addition of two terms):

[TABLE]

When $\lambda=+$ the first term corresponds to the generation of an electron-hole pair before scattering from the Coulomb potential. The second term is not included in the usual approach to the Wannier equation, but has a non-negligible contribution to the optical conductivity.

3 Approximated equation-of-motion for the polarization operator: Linear response

In this section we derive the frequency response, in the linear regime, for equation (19) using the thermodynamic equilibrium density function given by the Fermi-Dirac distribution (28). The external electrical field is written as ${\cal E}(t)={\cal E}_{0}\rme^{\mathrm{i}\omega t}$, and the linear-response is obtained from the terms proportional to $\rme^{\mathrm{i}\omega t}$, $p^{s\tau}_{\lambda}(\mathbf{k},t)=p^{s\tau}_{\lambda}(\mathbf{k},\omega)e^{\mathrm{i}\omega t}$. In this regime we neglect the non-linear term ${\cal W}^{s\tau}_{\mathbf{k}\lambda}(t)$ (25) that only contributes to second-harmonic and higher frequency terms, and the density term ${\cal D}^{s\tau}_{\mathbf{k}\lambda}(t)$, that contributes indirectly to the linear response through equation (20). With these considerations, the equation-of-motion (19) becomes:

[TABLE]

with ${\cal B}^{s\tau}_{\mathbf{k}\lambda}(\omega)$ obtained replacing $p^{s\tau}_{\lambda}(\mathbf{k},t)$ by $p^{s\tau}_{\lambda}(\mathbf{k},\omega)$ in equation (34).

Once equation (35) is solved for the transition probabilities $p^{s\tau}_{\lambda}(\mathbf{k},\omega)$, the total polarization $P(t)=P(\omega)e^{\mathrm{i}\omega t}$ can be obtained from the expectation value of equation (15):

[TABLE]

and the optical conductivity follows from the macroscopic relation between the polarization current density and the polarization density $\mathbf{J}(t)=\partial_{t}\mathbf{P}(t)$:

[TABLE]

The presence of the excitonic term, equation (34) in ${\Omega^{1}}^{s\tau}_{\mathbf{k}\lambda}(\omega)$, makes the equation of motion (35) a system of two coupled Fredholm integral equations of the second kind for the transitions probabilities $p^{s\tau}_{\lambda}(\mathbf{k},\omega)$, $\lambda=\pm$, that has to be solved for each spin/valley pair $s\tau$. The correspondent homogeneous equation [that can be obtained making ${\cal E}_{0}=0$ in equation (35)], corresponds to a Fredholm integral equation of the first kind. The solution of the latter will be explored in the next section.

3.1 Study of homogeneous Bethe-Salpeter Equation: Excitonic states

The set of equations (35), in the homogeneous case, corresponds to the Bethe-Salpeter equation for the exciton wave function $\psi_{nl}^{s\tau}$ with energy $E^{s\tau,\mathrm{exc}}_{nl}$. In this limit this equation is also known as the Wannier equation [24]. Once the excitonic wave functions are known we can calculate the absorbance coefficient ${\cal A}(\omega)$ using Elliot’s formula [24], in a form appropriate for TMDC’s; this formula is derived in C following a procedure described in [17]; this approach leads to

[TABLE]

where $M^{s\tau}_{n\ell}$ is the oscillator strength, given by:

[TABLE]

with $\ell$ and $n$ the angular and radial quantum numbers, and we have explicitly reintroduced the Fermi velocity $v_{F}$ for defining the oscillator strength as a dimensionless quantity.

In the system of equations (35) we neglect the non-resonant term $p^{s\tau}_{-}(\mathbf{k},\omega)$, consider the system at zero temperature, and at the charge neutrality point, $\Delta f^{s\tau}_{\mathbf{k}}=-1$. Thus, we have for the homogeneous problem an integral equation for $p^{s\tau}_{+}(\mathbf{k},\omega)$:

[TABLE]

If we also neglect the exchange self-energy term in the previous result, equation (40) is formally equal to the Bethe-Salpeter equation for the two-body electron-hole wave function obtained in the center-of-mass reference frame in a gapped Dirac system. We can write equation (40) in the following matrix form:

[TABLE]

with $E_{\mathrm{exc}}=\omega$ the Exciton energy and $K^{BS}$ the integral operator of equation (40). We can use the cylindrical symmetry to write the eigenfunctions of equation (40) as:

[TABLE]

with $\ell$ the angular quantum number (note that we have omitted the dependence of $\psi^{s\tau}_{n\ell}(k)$ in $\omega$). The extra phase in definition (42) allows to write equation (40) in the heavy mass limit $m\rightarrow\infty$ as a hydrogen-atom equation in momentum space with a screened potential [13]. The introduction of a factor $\rme^{\mathrm{i}\theta}$ to classify the excitons can be seen as arbitrary. Since the spinors, given by equation (9), also have an arbitrary global phase that propagates to the four-spinor product (27), the best way to define the $s$-wave is by a limit condition, that is, taking the limit $m\rightarrow\infty$ we should recover the spectrum of the hydrogen atom.

Substituting (42) into (40), and using the orthogonal relations for the wave function $\psi^{s\tau}_{n\ell}(k)$, equation (40) becomes the following Wannier formula:

[TABLE]

with the kernel $\tilde{T}_{\ell}^{s\tau}(k,q)$ given in D. Equation (43) was solved before in references [17, 18] without the exchange self-energy contribution in $\tilde{\omega}^{s\tau}_{+k}$. Here we include the effect of the exchange term.

The results for convergence of the binding energies, $E^{\mathrm{binding}}_{n\ell}=\Delta_{s\tau}-E_{n\ell}$, with $\Delta_{s\tau}$ given by equation (30), of the integral equation for the $s$-wave ($\ell=0$) and the $p$-waves ($\ell=\pm 1$) are presented in columns labeled “3000” and “GL” of table 2 (in it we also discuss convergence issues of our numerical methods). In this table we compare our results with those of reference [18], and because of this we have neglected the exchange correction as those authors also did. The wave functions for the $s$-wave are presented in Fig. (6), where we can see the usual increase of nodes for higher modes. The breakdown of degeneracy in $p$-waves comes from the $\ell$ dependence in the matrix element of the four body spinor product, that appears inside the integral (97), and is a consequence of the Berry curvature of the Dirac Hamiltonian [19, 20].

The Rydberg series for excitons including exchange corrections is presented in figure 7 for all four TMDC’s considered in this work with the parameters of table 2. There are four combinations of spin/valley, but the time-reversal symmetry reduces this number to two independent combinations. Each of these will generate a distinct Rydberg series, that we call A for the lowest fundamental energy and B for the highest one. Note that the A and B series are split due to SOC. For WS2 we compare the Rydberg series with the experimental work of ref. [39]; our results show a good agreement with the experimental data, with a deviation smaller than $5$ meV.

3.2 Integral equation for the vertex function

With the knowledge of the solution of the homogeneous equation, we come back to the integral equation (35), where we add a phenomenological interband relaxation rate $\gamma_{p}$ to include disorder effects.

We write $p^{s\tau}_{\lambda}(k,\theta,\omega)={\cal E}_{0}\Psi^{s\tau}_{\lambda}(k,\theta)$ and proceed as we did in equation (42) expanding $\Psi^{s\tau}_{\lambda}(k,\theta)$ in the eigenstates of the homogeneous Bethe-Salpeter equation and angular momentum states:

[TABLE]

and as a consequence:

[TABLE]

where $\delta_{\ell,s}$ is the Kronecker-delta, and the kernels $T^{1,s\tau}_{\lambda,\ell}(k,q)$ and $T^{2,s\tau}_{\lambda,\ell}(k,q)$, and the velocity angular decomposition $v^{s\tau}_{\pm,\lambda}(k)$ are given in D.

The diagrammatic representation of equation (45) is shown in figure 1, where the internal electron-hole legs are understood to be dressed by the exchange interaction.

3.3 Bright and Dark excitons: optical selection rules

We are now in position to discuss the conditions that an exciton can absorb a photon. From the coupled set of equations (45), only $\ell=0$ ($s$-excitons) and $\ell=-2$ ($d$-excitons) contribute to the optical conductivity, but from our numerical calculations, the $\ell=-2$ mode contribution is negligible for real TMDC’s parameters. After performing the $\theta$-integration in equation (36), we obtain for the polarization:

[TABLE]

Finally, the numerical procedure to calculate the optical conductivity is the following: we solve the integral equation (45) to determine the eigenfunctions $\psi^{s\tau}_{\lambda\ell}(k)$, calculate the polarization from the integral in equation (46), and lastly the optical conductivity is calculated from relation (37). The absorbance for a MoS2 suspended sheet for a TEM wave with normal incidence can be obtained from the optical conductivity as [42]:

[TABLE]

with $f(\omega)=\sigma(\omega)/\sigma_{0}$ and $\sigma_{0}=e^{2}/4\hbar$.

4 Excitons at the $\Gamma$-point

To describe accurately the optical absorption in the frequency domain after the two first excitonic peaks (A and B), that is the region roughly located in the interval $2.3-3.5$ eV, we need to describe the excitonic effects due to electronic transitions at the $\Gamma$-point. This implies going beyond the $\mathbf{k}\cdot\mathbf{p}$ model at $\mathbf{K}$-point. To accomplish this, we use the three-band model of Liu et al. [25], which was shown to describe accurately the GGA band structure. We use the same equation of motion approach introduced in previous section. The dipole matrix element is calculated using a Peierls approximation [43, 23]:

[TABLE]

which takes in account vertical interband transitions only.

To proceed with the discussion about excitons (actually excitonic resonances) at the $\Gamma$-point, we have to look in detail into the TMDC band structure depicted in figure 9; this band-structure was calculated using a full relativistic method (see figure caption for details), that is necessary to correctly account for the spin-orbit coupling [44] in TMDC’s. Very close to the $\Gamma$-point, the top of the valence band (that from now on we label band [math]) and the last four conduction bands (that, from here on, are labeled $1$, $1^{\prime}$, $2$ and $2^{\prime}$, in increasing energy order) are essentially due to contributions from electrons belonging to the transition-metal $d$-orbitals. We note in passing that the bands [math], $1$, and $2$ are also used in the effective Hamiltonian valid near the $K$-point. The four lowest conduction bands at the $\Gamma$ point (labeled $3$, $3^{\prime}$, $4$ and $4^{\prime}$, in increasing energy order) are mostly composed of $p$-orbitals from the chalcogenides atoms. From here on, we do not consider SOC effects in the three band model as they are very small at the $\Gamma-$point (see right panels in Fig. 9). Therefore, we drop the prime notation of the bands, that is, the bands $i$ and $i^{\prime}$ are treated as spin-degenerated (at the computation level this amounts to a multiplicative factor of two affecting the absorbance curves).

The mirror symmetry, with respect to the plane formed by the metallic atoms, is important to discuss the optical properties of TMDC’s. The bands [math], $1$ and $2$ all have even symmetry, while the bands $3$ and $4$ both have odd symmetry. Therefore, the dipole matrix elements between bands with different mirror symmetry vanish (note that the mirror symmetry refers here to the $z-$coordinate). As a consequence, the existence of excitons composed of holes from the [math] band and electrons from the bands $3$ and $4$, they are optically dark under a single-photon experiment. This implies that we restrict the calculation of excitonic effects near the $\Gamma-$point considering only the optical properties of excitons composed of holes from band [math] and electrons from bands $1$ and $2$. These optical transitions generate a modified Rydberg series, which is a consequence of the non-parabolic dispersion relation of the optical band $[$see equation (52)$]$ and finite Berry curvature at the $\Gamma$-point and as well as from the form of the Keldysh potential.

To calculate the optical properties of the excitons at the $\Gamma$-point, we use the same equation of motion method developed in previous section. We define transition probabilities $p_{0i}(\mathbf{k},t)=\langle\hat{a}^{\dagger}_{i\mathbf{k}}(t)\hat{a}_{0\mathbf{k}}(t)\rangle$, that represent the annihilation of an electron at the valence band [math] and a creation of an electron at the conduction band $i=1,2$. The equation of motion is again given by (18), making $\lambda=0$ and $\lambda^{\prime}=1,2$, the latter two values represent the two possible bright excitons at the $\Gamma$-point. On the other hand, the $H_{0}=H_{0}^{3\mathrm{b}}$ Hamiltonian is the three-band model given by Liu et al. [25], that describes, up to the next-nearest-neighbor order, the effective tight-binding model between the metal atoms M={Mo,W}.

To calculate the equation of motion (18), we again need the commutators $[\hat{\rho}_{0i},\hat{H}^{3\mathrm{b}}_{0}]$, $[\hat{\rho}_{0i},\hat{H}_{I}]$, and $[\hat{\rho}_{0i},\hat{H}_{ee}]$, $i=1,2$. The first two commutators give analogous results to those of the previous section, and the last one is given in equation (66). We are interested only in the term equivalent to the Bethe-Salpeter equation (encoded in the Rabi frequency renormalization term):

[TABLE]

which is obtained from equation (66) from the terms with $\lambda_{1}=0$, $\lambda_{2}=0,i$ and $\lambda_{3}=i$. For the model we are considering, the four-body spinor is given by:

[TABLE]

and is obtained numerically from the diagonalization of the $3\times 3$ matrix defining the three-band model, with $\phi_{i}$ the wave function of the $i$-band.

In our calculation, the exchange self-energy is included as an energy shift from the renormalization of the band gap [see equation (32)]. We also ignore spin-orbit effects in the calculations at the $\Gamma$-point, as these are very small. The equation that we need to solve for obtaining the transition probability $p_{0i}(k,\theta,\omega)$ reads:

[TABLE]

with $\omega_{i}(k,\theta)=E_{ik,\theta}-E_{0k,\theta}$, $E_{ik,\theta}$, and $E_{0k,\theta}$ the eigenvalues of the three-band Hamiltonian $H^{3\mathrm{b}}_{0}$, $\Delta f_{i}=f_{ik,\theta}-f_{0k,\theta}$ is the difference in occupation numbers between bands $i$ and [math] (given in terms of the Fermi-Dirac function), and the dipole element $d^{\Gamma}_{i}(k,\theta)=\langle\lambda\mathbf{k}|\hat{x}|\lambda^{\prime}\mathbf{k}\rangle$, with the expectation value calculated using equation (48).

Since we only need the band-structure near the $\Gamma-$point we approximate the optical band structure near this point by the Fourier series:

[TABLE]

with $h^{i}_{\ell}(k)$, $\ell=0,6,12$ are polynomials of degree six (this expression is valid up to momentum values of the order of $2\pi/(3a_{0})$, where $a_{0}$ is the lattice parameter). Expression (52) describes accurately the optical band structure near the $\Gamma$-point. For $k\rightarrow\infty$ the optical band approximation (52) diverges. Although the contributions for $k\rightarrow\infty$ becomes negligible to the excitons’ wavefunction, the approximation (52) makes the numerical convergence faster. Using the angular decomposition (10), $p_{i0}(k,\theta,\omega)={\cal E}_{0}\sum_{\ell}c_{\ell}^{i}(k)e^{\mathrm{i}\ell\theta}$ (note that we have omitted the dependence of $c_{\ell}^{i}(k)$ in $\omega$), we can write the Bethe-Salpeter equation (51) as (where we have made $\Delta f_{i}=-1$ since are interested in a neutral system):

[TABLE]

and $\tilde{\omega}^{i}_{\ell}(k)=\int_{0}^{2\pi}\frac{d\theta}{2\pi}\,e^{\mathrm{i}\ell\theta}\omega_{i}(k,\theta)$.

The previous equation couples coefficients $c_{\ell}^{i}$ with different angular momentum numbers $\ell$ through two terms: the kinetic term $\sum_{\ell^{\prime}}\tilde{\omega}^{i}_{\ell^{\prime}}(k)c^{i}_{\ell-\ell^{\prime}}(k)$ and the electron-electron interaction term (Rabi frequency renormalization term). The kinetic term couples only coefficients having $\ell^{\prime}=0,\pm 6,\pm 12$, a consequence of the sixfold symmetry of equation (52). The potential term also couples coefficients $c^{i}_{\ell}(k)$ with different angular momentum values, but this term gives a negligible contribution to the optical response when $\ell\neq\ell^{\prime}$. This is a consequence of the fast vanishing of the potential $V(k,q,\theta-\theta^{\prime})$ whenever $\theta-\theta^{\prime}\neq 0$ and $\theta-\theta^{\prime}$ is varied. Therefore, we can replace the four-body spinor function $F^{3b}_{0i}(k,q,\theta,\theta^{\prime})$ by its average angular value as follows:

[TABLE]

and the effective potential is given in this approximation by:

[TABLE]

The approximation (55) and (54) keeps the hermiticity of equation (53). Finally, we replace the potential term in equation (53) by (55), where we arrive at the coupled set of integral equations:

[TABLE]

To solve (56), the summation in $\ell^{\prime}$ gives five additional terms [$\ell^{\prime}=0,\pm 6,\pm 12$; see equation (52)] that are coupled together. This generates a hierarchy of equations for the coefficients $c^{i}_{\ell}(k)$. Therefore, the solution of equation (56) has to be truncated at some $\ell$ value. In this procedure we have assumed that the contributions above $c^{i}_{18}(k)$ are vanishing small. This is confirmed by figure 10, which shows that for $\textsection\ell=5$ the contribution is already small (note that for this curve we have all the coefficients $c^{i}_{\ell}(k)$, with $\ell=-7,-1,11,17$ entering the calculation of the conductivity, see below).

In terms of the coefficients $c^{i}_{\ell}(k)$ the conductivity is computed as follows. The expectation value of the polarization operator can be calculated as we did in section 2, and results in:

[TABLE]

where we account for the spin degeneracy introducing a factor of two. The conductivity can be obtained from equation (37), and we can separate the contribution for each band $i$. Performing the angular integral in the equation for $P(\omega)$ we obtain:

[TABLE]

Once the coefficients $c^{i}_{\ell}(k)$ are determined from the solution of (56) the conductivity follows from the previous equation.

The solution of (56) also give us the excitonic wave functions in momentum space. The results for the first excitonic energy, for each angular momentum mode, is shown in figure 11 for the exciton composed from an electron in band $1$, and in figure 12 for an electron in band $2$. From a careful inspection of figures 11 and 12, we can see that the nodes of exciton with band index $i=1,\ell=0$ lies along the $\Gamma-\mathbf{K}$ direction, while the nodes of exciton with band index $i=2,\ell=0$ lies along the $\Gamma$-M point.

5 Results

In this section we perform a thorough analysis of the absorption spectrum of four TMDC’s. For computing the absorbance, the optical conductivity is needed. Taking the example of MoS2, the decomposition of the real part of the optical conductivity, coming from different angular momentum contributions of the exciton at $\Gamma$-point, associated with the transition $0\rightarrow 1$, is shown in figure 10; remember that each contribution is composed of $\ell=0$, $\ell\pm 6$, and $\ell\pm 12$ angular momentum components.

It is important to introduce here a note on notation: the peak at lowest energy is denoted by A$=1s$ and the next Rydberg energy level in the A-series is denoted by A’$=2s$; this corresponds in a given valley and to a given spin projection. In the same valley, and for the other spin projection, the peaks belong to the B-series, with the lowest energy is denoted by B=$1s$ and the next one by B’$=2s$. For MoX2 TDMC’s the energy order is A, B, A’, and B’, whereas for WX2 TMDC’s the energy order is A, A’, B, and B’. This agrees with the notation introduced in figure 7.

The absorbance, and the real and the imaginary parts of the optical conductivity, for four TMDC’s considered in this work, are shown in figure 13, with the parameters of table 1. That is, in this figure we do not try to fit the data but simply use the parameters characterizing the potential and the band-structure of the TMDC’s given in other papers. In figure 14, on the contrary, we fit the A peak position changing $r_{0}$ and we also add a chemical potential, since, as noted in Refs. [48, 54] all TMDC’s samples have a certain and undetermined amount of negative doping. We note in passing that at the time of writing different experiments report distinct percentages for absorption of radiation for two, supposedly identical, TMDC’s. Table 3 gives, from four different references, the measured values of the absorbance of MoS2 samples; as it can be seen the values fluctuate among different experiments. Also, our model predicts, at low temperatures, larger absorption peaks than those measured at room temperature. This result makes sense, but when we increase the temperature we never obtain values as small as those reported in the experiments for MoS2. It is now known [55] that excitonic spectrum of TDMC’s samples in SiO2 are strongly influence by the disorder of the substrate. In this reference it is shown that encapsulated samples in h-BN have much narrower excitonic peaks. Therefore our results should agree with absorbance measurements in these encapsulated samples (measurements yet to be made).

Next, we analyze each aspect of the optical spectrum of each TMDC and compare our results with the experimental measurements available to date in a large frequency window. The parameters used in our calculations are: (i) at the $\mathbf{K}$-point we used the values in table 1 and a broadening $\gamma_{K}=50$ meV; (ii) at the $\Gamma$-point we used the GGA parameters of the three-band model given by Liu et al. [25], the same Keldysh parameters of table 1, and a broadening $\gamma_{\Gamma}=100$ meV. Note that exception made to the broadening parameters, all the other values were taken of the literature and no attempt was made to choose them in order to fit the data, with exception to the case reported in figure 14.

• •

MoS2 The two first peaks in the absorbance, A and B, correspond to the A($1s$) and B($1s$) excitons. The different position of the two peaks is a consequence of SOC splitting of the bands. The last two peaks in the absorbance spectrum, $C_{01}$ ans $C_{02}$ (having about the same intensity –see the conductivity curve), correspond to the sum of different angular momenta contributions from the $\Gamma$-excitons (actually excitonic resonances). The third ($C_{01}$) peak is associated with the transition from the top of the valence band to both the 1 and 1’ conduction bands at the $\Gamma$-point (conduction bands number 5 and 6 in figure 9); we have considered these two degenerated since SOC is small in this case. Finally, the fourth peak ($C_{02}$) comes from transitions connecting the top of the valence band and the 2 and 2’ conduction bands (also taken degenerated; conduction bands 7 and 8 in figure 9).

The real part of the conductivity follows closely the absorbance spectra, as expected. Usually, the imaginary part of the conductivity from a single excitonic contribution is negative for $\hbar\omega<E_{b}$ and positive for $\hbar\omega>E_{b}$, where $E_{b}$ is the binding energy, a result that can be obtained by inspection of Elliot’s formula for the optical conductivity (91).

Let us now discuss the differences between experimental data and our model. We note that the rigid shift to the left performed by Wu et al. [18], and Steinhoff et al. [23] is not necessary in our case. The difference in intensity of A and B peaks is probably a consequence of the phonons that exist at finite temperature. This effect was not considered in this work but was shown to be important for the peaks’s broadening [16, 7].

Lastly, we discuss the excitonic effects at the $\Gamma$-point. The optical measurements identify only one peak, which seems to correspond to the $C_{02}$ exciton. The work of Qiu et al. [16] obtain a rich structure of peaks in this region that was washed out when they include quasi-particle lifetimes from phonon terms.

• •

MoSe2

The aborbance spectrum of this TMDC share many similarities with MoS2: two peaks from the $\mathbf{K}$-point split (A and B) by the SOC and two wider peaks from the $\Gamma$-point are also present. The exciton at the $\Gamma$-point contributes with two peaks at $\sim 2.4$ eV and $\sim 2.7$ eV. The experimental data shows a single peak at $2.6$ eV. This discrepancy comes possibly from the phonons already discussed for the MoS2. Overall there is a good agreement between the data and the calculated curves, both in position of the peaks and in intensity.

The imaginary part of the conductivity is only positive for frequencies $\hbar\omega>2.6$ eV, meaning that exciton-polaritons can only be excited for energies in the visible.

• •

WS2

For this material we note the very good agreement of the position and magnitude of the calculated A peak in comparison with those in the experimental data. We also see that the second experimental peak coincides with a small computed peak from the 2s state (A’) associated with the series of first exciton A-peak (see figure 7). There is at least three reports [49, 50, 51] of measurement of the A’ peak in WS2 in the temperature range of 4-300 K. Unfortunately, in the literature the A’ peak has been dubbed B, using an analogy with the MoX2 case. However, looking at the central panel of figure 13 we clearly see that the A’$=2s$ peak appears at lower energy than the B$=1s$. Note that from our analysis we can separate each spin/valley contribution. Also note that the A’ peak has a similar absorbance to the experimental one (identified in the experimental literature on WX2 TMDC’s as B, because it is the second to appear in the energy scale). Studying the dependence of light absorption of different peaks on an external magnetic field, for breaking spin degeneracy, together with the use of strong circular polarized light to populate the two valleys differently[52], is a possible way of clarifying the microscopic origin of the different peaks.

The third theoretical (B) peak (which is the SOC counterpart of the first peak) is absent in the experimental data. Note that from figure 7, all but one (1$s$=B) contributions from the B family of peaks are excitonic resonances (above the interacting band gap). The proximity of the B-peak to the continuum may provide a scattering channel to transfers spectral weight from this peak to the resonances in the continuum. An additional and possible mechanism is based on extrinsic doping of these materials as shown in figure 14. It has been shown that doping has a strong effect in attenuating the excitonic peaks [14], specially the B-peak in MoS2. There is no reason to believe that the same mechanism would not work in WX2 TMDC’s. Indeed, from figure 2 we expect a strong attenuation of the high-energy excitonic peak whereas the low energy one should survive. This should happen since the doping with electrons tends to block first the higher energy transition whereas maintaining the low energy one. In figure 14 we see a comparison of the absorption spectrum of WS2 with the data taking into account the effect of doping; the agreement is excellent. The suppression of the B-peak is evident from our results, thus confirming doping by electrons as a possible mechanism for suppressing the B excitonic state.

The first excitonic resonance ($C_{01}$) at the $\Gamma$-point is in very good agreement with the experimental one, while the second excitonic ($C_{02}$) resonance at the $\Gamma$-point is at an energy range above the measured one (although its intensity is rather small). Therefore nothing can be said about the possible agreement with the experimental data, since this does not cover that spectral region. Lastly the imaginary part of the optical conductivity becomes positive above the energy $\sim 2.6$ eV nd therefore the system can support exciton-polaritons in that spectral region.

• •

WSe2

We end our analysis with a comparison between the calculated absorbance curve and the one measured for WSe2. For this material the disagreement between the calculated curves and the experimental data is the largest of the four TMDC’s studied in this work. Indeed, the data seems stretched relatively to the calculated curves. The first peak in the WSe2 absorbance spectrum is in very good agreement with the experimental data, with a difference in position less than $0.1$ eV.

As in the case of WS2, we see that the B-family peak is present in the data as a small shoulder. The third and fourth experimental peaks, when compared with our theoretical model, come from resonances at the $\Gamma$-point. The theoretical calculations show a red shift of about $0.2$ eV for these two peaks, indicating that higher order exchange corrections, which reshape the band structure around the $\Gamma$-point, might be important.

The imaginary part of the optical conductivity is positive above $\sim 2.8$ eV, thus allowing for excitons-polaritons.

Next we present an analysis of the effects associated with changing the Fermi energy and the Keldysh potential parameter $r_{0}$. Given a Fermi energy the parameter $r_{0}$ can be adjusted to fit the A peak. Results of this procedure are shown in figure (14) for WS2. We can see an excellent agreement between our results and the experimental curve. This highlights the importance of a finite Fermi energy in describing the experimental data. As noted before a finite Fermi energy comes from the spontaneous negative doping observed in TMDC’s samples. The better agreement with the data shown in figure (14) relatively to the results of figure (13) shows the non-negligible effect of the doping in the optical properties. On the other hand, the parameter $r_{0}$ should also be a function of the electronic density. At the moment of writing this dependence is unknown.

One aspect that our calculation does not take into account in an exact way is the self-consistent solution of the exchange energy. Since this calculation is outside the scope of this work, we can mimic it using a different value of $\Lambda_{2}$ [see Eq. (3)]. This leads to a narrow A peak and a broaden B peak in WS2, as seen in the experimental data. In this regime, the B peak is no longer an exciton but rather an excitonic resonance. The mechanism leading to the broaden of the B peak can be explained by the self-consistent solution of the exchange energy. For a given carrier density, the iterative calculation of the exchange energy reduces its value and therefore the importance of the doping increases for the lowest conduction band. In WS2 the effect is much stronger in the lowest band than in the next conduction band due to the large spin-orbit splitting. This mechanism due to exchange increases the splitting on the two conduction bands.

Another aspect of the doping is its influence on the decreasing of the band gap. We show in figure (15) the dependence of the band gap and the exciton energies on the Fermi energy. We can see that increasing Fermi energy makes the binding energy (difference between the thick blue curve and all the others) smaller. The energy of the first excited state (squares) increases with the doping while the energy the second (triangles) and third (circles) have the opposite behavior. We also see that it exists a critical doping that makes the exciton states collapsing into resonances when they merge with the band gap. For the energy of first and third excited states we see the same qualitative behavior as measured in Ref. [54].

This concludes the analysis of our theoretical results when compared with the experimental data. Globally, the agreement is good, but some points need further research. Measurements performed at low temperatures in encapsulated TMDC’s using hexagonal boron-nitride should reveal the fine structure of the excitonic spectrum predicted by our model.

6 Discussion and conclusions

In summary, we performed a study of excitons in TMDC’s monolayers including in the same foot both excitons at the $\mathbf{K}$- and $\Gamma$-points. The excitons at the $\mathbf{K}$-point were calculated with a gapped Dirac equation including electron-electron interactions and SOC. The excitonic resonances at the $\Gamma-$point were calculated with the tight-binding three-band model expanded around that point in the Brillouin zone. We compared our theoretical results with the experimental data available from reference [11]. We clarified the microscopic origin of each observed excitonic peak and discussed the reasons for some disagreement between our theoretical model and the experimental data. Note that the measurements where made at room temperature. Therefore, the effect of a self-energy, which will be energy dependent, due to electron-phonon interactions might play an important role in modeling the absorbance spectrum at room temperature. We note here that our equation of motion method also allows for treating electron-phonon interactions at the expenses of a more lengthly calculation.

Also, as noted by Mak et al. [48]: “Spontaneous negative doping, presumably from defects within the MoS2 layer and/or substrate interactions, has been commonly reported in mechanically exfoliated samples”. This seems be the reason [14] why the B-series is not visible in WX2 when the material is electron-doped (see figure 14). To conclude, given the uncertainties in the experimental data reported in table 3 we consider the agreement between our calculation and the data to be quite good.

We have also studied the variation of the $1s=$A peak with the dielectric function of the capping medium (results not shown). We found that the A-peak position varies little with $\epsilon_{m}$. This happens because the exchange energy correction compensates the binding energy coming from the BSE.

Although we have considered in this work the response to linearly polarized electromagnetic radiation, it is simple to generalized the calculation to circularly polarized one. This would allow us to discuss the additional appearance of more selection rules associated with spin.

Finally, we note that our method can be easily generalized to the calculation of excitonic effects in other 2D materials, such as phosphorene, a line of research we will pursue in a forthcoming paper. This case will be particularly interesting due to the strong anisotropy of the material.

A future working direction is the inclusion of density and non-linear terms coming from the commutator of the density matrix with the electro-electron interaction. In particular the latter term will be important for the discussion of excitonic-assisted nonlinear optics in TMDC’s (see also [53]). Note that contrary to 3D semi-conductors the optical spectrum is dramatically affected by excitonic excitation even for frequency ranges well above the non-interacting gap. Thus, a description of the nonlinear optics properties of TMDC’s using a single electron picture is doomed to fail.

Acknowledgments

A. J. Chaves acknowledges a scholarship from the Brazilian agency CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico). N.M.R. Peres and R.M. Ribeiro acknowledge support from the European Commission through the project “Graphene-Driven Revolutions in ICT and Beyond” (Ref. No. 696656), project PTDC/FIS-NAN/3668/2014, and the Portuguese Foundation for Science and Technology (FCT) in the framework of the Strategic Financing UID/FIS/04650/2013. T. Frederico thanks for the support of the Brazilian Agencies CNPq and FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo). The authors acknowledge Fanyao Qu and Alan MacDonald for useful discussions about the numerical solution of the homogeneous Bethe-Salpeter equation, and Ermin Malic for sharing with us useful supplementary notes.

Appendix A Calculation of the commutators

The equation of motion for density matrix $\hat{\rho}_{\lambda\lambda^{\prime}}^{s\tau}$ will be determined in this appendix. We need do calculate the following commutators: $[\hat{H}_{0},\hat{\rho}^{s\tau}_{\lambda\lambda^{\prime}}((\mathbf{k},t)]$, $[\hat{H}_{\mathrm{ee}},\hat{\rho}^{s\tau}_{\lambda\lambda^{\prime}}(\mathbf{k},t)]$ and $[\hat{H}_{I},\hat{\rho}^{s\tau}_{\lambda\lambda^{\prime}}(\mathbf{k},t)]$. For the first of these we have:

[TABLE]

For the external field, we only consider the interband terms of the Hamiltonian $\hat{H}_{I}$ as we are discussing a neutral system. In this condition $\hat{H}_{I}$ reads:

[TABLE]

and therefore:

[TABLE]

with the corresponding expectation values:

[TABLE]

and

[TABLE]

Finally, the commutator with the electron-electron interaction reads:

[TABLE]

The expectation value of four body operators is truncated at the RPA level:

[TABLE]

and within this approximation, the expectation value of equation (64) is given by ($S=L^{2}$, the area of the system):

[TABLE]

where we have used the property $V(\mathbf{q})=V(-\mathbf{q})$. Equation (66) is valid for a system composed of any number of bands. For the particular case of two band systems, as in the case given by the Hamiltonian (6), we split equation (66) into various terms. For $\lambda=\lambda^{\prime}$ we have:

[TABLE]

and for $\lambda^{\prime}=-\lambda$ we find:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally we need to add (59), (62), and (67) to obtain the equation of motion for $n^{s\tau}_{\lambda}$, and (59), (63), and (68) to obtain the equation of motion for $p^{s\tau}_{\lambda}$.

Appendix B Overlap of the four-body wavefunctions

The four-body overlap functions are explicitly defined below for the massive Dirac Hamiltonian:

[TABLE]

For simplicity of writing, we omit in this appendix the superscript $s\tau$ in the $F$’s-functions and in the energy $E_{k}^{s\tau}$. For the case $\lambda_{1}=\lambda_{4}$ and $\lambda_{2}=\lambda_{3}$ the overlap function reads:

[TABLE]

whereas when $\lambda_{1}=\lambda_{4},\lambda_{2}=-\lambda_{3}$ we find:

[TABLE]

Finally, in the conditions $\lambda_{1}=-\lambda_{4},\lambda_{2}=-\lambda_{3}$ we have:

[TABLE]

When $\lambda_{1}=-\lambda_{4}$ and $\lambda_{2}=\lambda_{3}$ we have the following symmetry:

[TABLE]

that is, in expression (76) we have an identity upon the exchange of indexes $\lambda_{1}\leftrightarrow\lambda_{2}$, $\mathbf{k}_{1}\leftrightarrow\mathbf{k}_{2}$.

Appendix C Derivation of Elliot’s formula

The solution of the homogeneous problem presented in equation (43) can be used to calculate the optical conductivity of the system. We now detail the derivation of Elliot’s formula for TMDC’s

First we decompose the excitonic wave function into a complete set of eingenfunctions of (43):

[TABLE]

where we have separated the discrete and continuum states of the exciton spectrum, with $n$ and $q$ the radial quantum numbers, respectively. We further recall the orthogonality relations:

[TABLE]

The non-homogeneous equation (40), after we substitute the expansion for $\Psi^{s\tau}_{\ell}(k)$ into the eigenfunctions of the Kernel $K^{BS}$ (41), and integrating in $\theta$, becomes:

[TABLE]

where the angular decomposition of the velocity (14),

[TABLE]

is composed of two terms:

[TABLE]

remembering that the angular decomposition has an extra $\theta$ phase. The explicit expression for $v^{s\tau}_{\ell,\lambda}(k)$ reads:

[TABLE]

Using the eigenfunction orthogonality and neglecting the continuum part $\psi^{s\tau}_{\ell}(q,k)$, we arrive at:

[TABLE]

from where it follows the coefficients $c^{s\tau}_{\ell s}$:

[TABLE]

Using the last result, the exciton contribution to the polarization is, using equation (36), given by:

[TABLE]

remembering here that we are using units such that $v_{F}=\hbar=e=1$. Note that the $\theta$ integral has been performed, and we are summing over all the spin/valley indexes. If we use the definition of the weight function $M^{s\tau}_{n\ell}$ (39) we have:

[TABLE]

where we have introduced a phenomenological relaxation rate $\gamma$. The 2D optical susceptibility comes from $\mathbf{P}=\varepsilon_{0}\chi_{2D}{\cal\mathbf{E}}$:

[TABLE]

where we reintroduce the units, and the conductivity reads:

[TABLE]

Finally, the absorbance coefficient ${\cal A}(\omega)=1-{\cal T}(\omega)-{\cal R}(\omega)$, where ${\cal T}(\omega)$ [${\cal R}(\omega)$] is the electromagnetic transmission [reflection] for a TEM wave, is given by:

[TABLE]

where $\alpha\approx 1/137$ is the fine-structure constant.

Appendix D The Bethe-Salpeter kernel

In this appendix we give the explicit forms of the BSE kernel. Firstly, from equation (22) we have:

[TABLE]

For the homogeneous case we only consider the first term in the previous equation and choose with $\lambda=-\lambda^{\prime}=+$, which corresponds to the resonant term. Thus we have the BSE kernel reading:

[TABLE]

Using the expression (76) for the $F^{s\tau}_{--++}(\mathbf{k},\mathbf{q})$, and the Keldysh potential (4), and after the angular decomposition

[TABLE]

we have the corresponding kernel $\tilde{T}^{s\tau}_{\ell}(k,q)$:

[TABLE]

where we have defined:

[TABLE]

[TABLE]

and

[TABLE]

where $\alpha\approx 1/137$ is the fine structure constant and $c$ is the speed of light.

For the kernel in the non-homogeneous BSE, we write the renormalization of Rabi frequency as:

[TABLE]

Thus, we have two kernels to consider:

[TABLE]

[TABLE]

Note that, in this case and contrary to the homogeneous BSE, we have to keep both terms in the renormalization of the Rabi frequency, as otherwise the real part of the optical conductivity would not have the correct positive sign. That is because both the resonant and off-resonance terms contribute to the optical conductivity, as is well known in the non-interacting case. Proceeding as before, the angular decomposition (95) of the kernels leads to:

[TABLE]

with

[TABLE]

and

[TABLE]

[TABLE]

where $\lambda=\pm$ and $\lambda^{\prime}=-\lambda$, and $I_{\ell}(k,q)$ is defined in equation (97).

The matrix element of the velocity operator appearing as the source term in the non-homogeneous BSE, was calculated in the previous section and reads:

[TABLE]

This concludes the presentation of the mathematical steps leading to the BSE equation discussed in the main text.

Appendix E Exchange correction around the $\Gamma$ point

From equation (66), the exchange self-energy correction to the transition energy between bands $i$ and $j$ is:

[TABLE]

Neglecting temperature and doping effects, for the three band-model of reference [25], the only term that contributes to the exchange self-energy is the one with $\lambda=0$ (the valence band):

[TABLE]

and we make $i=1,2$ and $j=0$.

To remove the integrable divergence at $q=0$, that comes from the Keldysh potential, we use polar coordinates leading to the need of computing the following integral

[TABLE]

where $Q_{i}(\mathbf{k},q,\theta)=F_{i00i}(\mathbf{k},\mathbf{k}-\mathbf{q})-F_{0000}(\mathbf{k},\mathbf{k}-\mathbf{q})$ and $q_{0}$ reads

[TABLE]

In figure 16 we show the exchange energy in the Brillouin zone for the $j=0$ to $i=1,2$ transitions, for the three-band tight-binding model and for the parameters of MoS2. The magnitude of the correction to the bare bands is of the order of 2.1-2.6 eV, which is a substantial value.