Nonlinear optical selection rules of excitons in monolayer transition metal dichalcogenides

TL;DR

This paper introduces an analytical method to calculate linear and nonlinear optical responses of monolayer TMDs, accurately capturing excitonic effects and selection rules, offering physical insights and efficiency over traditional BSE calculations.

Contribution

It presents a new analytical approach incorporating trigonal warping to determine excitonic optical responses, simplifying calculations while maintaining accuracy compared to Bethe-Salpeter equation methods.

Findings

Accurately reproduces optical spectra of MoS2

Provides detailed excitonic selection rules

Enables physical insight into resonant transitions

Abstract

We propose an analytical approach for calculating the linear and nonlinear optical (NLO) responses of monolayer transition metal dichalcogenides (TMDs) including excitonic effects. An effective Hamiltonian reproducing the trigonal warping (TW) of the energy dispersion, is used to derive analytical expressions for excitonic matrix elements. Based on this approach, we provide an informative diagram, which encompasses all excitonic selection rules. The diagram enables us to identify main transitions for the first-, second- and third-order optical responses. As a case study, we calculate the optical conductivity and second/third harmonic generation responses of monolayer MoS2 and demonstrate that the analytical approach accurately reproduces the spectra obtained using the Bethe-Salpeter equation (BSE). Therefore, the analytical approach provides an efficient and reliable method for calculating the NLO spectra as accurate as the full BSE method. Moreover, it enables us to obtain valuable physical insight into the fundamental transitions responsible for individual resonances, which is not straightforward in the full BSE method.

Tables

Table 1. Table 1: Energies of the first few excitons computed using the full BSE or analytical (or Dirac) approach for monolayer MoS 2 .

Exciton	$|1s⟩$	$|2p^{+}⟩$	$|2p^{-}⟩$	$|2s⟩$	$|3d^{+}⟩$
Dirac, A	1.909	2.114	2.128	2.179	2.220
BSE, A	1.872	2.088	2.103	2.151	2.197
Dirac, B	2.050	2.258	2.271	2.323	2.364
BSE, B	2.017	2.228	2.242	2.291	2.345

Equations

$${\mathbf{J}}^{(N)}(t)={\mathbf{\sigma }}^{(N)}\underset{N\mathrm{times}}{\underbrace{{\mathcal{E}}_{\omega }...{\mathcal{E}}_{\omega }}}\exp (-iN\omega t)+c.c.,$$

$$\displaystyle\sigma_{\alpha\beta}^{(1)}=C_{1}\sum_{m}\bigg{[}\frac{\Pi_{0m}^{\alpha}X_{m0}^{\beta}}{\hbar\omega-E_{m}}-\frac{\Pi_{m0}^{\alpha}X_{0m}^{\beta}}{\hbar\omega+E_{m}}\bigg{]}\,,$$

$$\displaystyle\sigma_{\eta\alpha\beta}^{(2)}=C_{2}\sum_{n,m}\bigg{[}\frac{\Pi_{0n}^{\eta}Q_{nm}^{\alpha}X_{m0}^{\beta}}{(2\hbar\omega-E_{n})(\hbar\omega-E_{m})}+$$

$$\displaystyle\frac{X_{0l}^{\alpha}X_{m0}^{\beta}}{(\hbar\omega+E_{l})(\hbar\omega-E_{m})}\Big{(}\frac{Q_{ln}^{\eta}\Pi_{nm}^{\mu}}{2\hbar\omega+E_{n}}-\frac{\Pi_{ln}^{\mu}Q_{nm}^{\eta}}{2\hbar\omega-E_{n}}\Big{)}\bigg{]}\,.$$

$$X_{0n}^{\alpha }=(X_{n0}^{\alpha }{)}^{\ast }\equiv -\sum _{\mathbf{k}}{\psi }_{\mathbf{k}}^{(n)}\frac{p_{vc\mathbf{k}}^{\alpha }}{i{\epsilon }_{cv\mathbf{k}}},$$

$$\displaystyle Q_{nm}^{\alpha}\equiv i\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{(n)*}\big{[}\psi_{\mathbf{k}}^{(m)}\big{]}_{;k_{\alpha}}\,,$$

$${\epsilon }_{cv\mathbf{k}}{\psi }_{\mathbf{k}}^{(n)}-\sum _{{\mathbf{k}}^{\prime }}{V}_{cv}^{\mathbf{k},{\mathbf{k}}^{\prime }}{\psi }_{{\mathbf{k}}^{\prime }}^{(n)}=E_n{\psi }_{\mathbf{k}}^{(n)},$$

$$H_{0}(\mathbf{k})=\left[\begin{array}[]{cc}\Delta+\lambda g(\mathbf{k})s&-\gamma f^{*}(\mathbf{k})\\ -\gamma f(\mathbf{k})&-\Delta-\lambda g(\mathbf{k})s\end{array}\right]\,,$$

$$f(\mathbf{k})\equiv e^{i{k}_{x}a/\sqrt{3}}+2e^{-i{k}_{x}a/(2\sqrt{3})}\cos (k_{y}a/2),$$

$$g(\mathbf{k})\equiv 2\sin (k_{y}a)-4\sin (k_{y}a/2)\cos (k_{x}a/\sqrt{3}),$$

$$f(\mathbf{k})\approx\frac{\sqrt{3}}{2}e^{i\pi/3}\big{[}-i(\kappa_{x}-i\kappa_{y}\tau)+\zeta(\kappa_{x}+i\kappa_{y}\tau)^{2}\big{]}\,.$$

$$X_{0n}^x(s)=X_0{\int }_0^{\infty }{\varphi }_{-1}^{(n)}\frac{{\eta }_{\kappa }+1}{i{\epsilon }_{\kappa }}\kappa \mathrm{d}\kappa ,$$

$$X_{0n}^x(p^{+})=0,$$

$$X_{0n}^x(p^{-})=X_0\zeta {\int }_0^{\infty }{\varphi }_{-2}^{(n)}\frac{{\eta }_{\kappa }^3+4{\eta }_{\kappa }+3}{2{\epsilon }_{\kappa }}{\kappa }^2\mathrm{d}\kappa ,$$

$$X_{0n}^x(d^{+})=X_0{\int }_0^{\infty }{\varphi }_{+1}^{(n)}\frac{{\eta }_{\kappa }-1}{i{\epsilon }_{\kappa }}\kappa \mathrm{d}\kappa ,$$

$$X_{0n}^x(d^{-})=0.$$

$$\frac{Q_{nm}^x(l\to l\pm 1)}{\pi a}=i{\int }_0^{\infty }{\varphi }_l^{(n)}\frac{\partial {\varphi }_{l\pm 1}^{(m)}}{\partial \kappa }\kappa \mathrm{d}\kappa$$

$$\pm i{\int }_0^{\infty }{\varphi }_l^{(n)}{\varphi }_{l\pm 1}^{(m)}({\eta }_{\kappa }+l\pm 1)\mathrm{d}\kappa ,$$

$$\frac{Q_{nm}^x(l\to l\pm 2)}{\pi a}=\frac{\zeta }{2}{\int }_0^{\infty }{\varphi }_l^{(n)}{\varphi }_{l\pm 2}^{(m)}(5{\eta }_{\kappa }-{\eta }_{\kappa }^3)\kappa \mathrm{d}\kappa ,$$

$$Q_{nm}^x(l\to l\pm 3)=0,$$

$$\frac{Q_{nm}^x(l\to l\pm 4)}{\pi a}=\frac{\zeta }{2}{\int }_0^{\infty }{\varphi }_l^{(n)}{\varphi }_{l\pm 4}^{(m)}({\eta }_{\kappa }^3-3{\eta }_{\kappa })\kappa \mathrm{d}\kappa .$$

$$i\hslash \frac{\partial {\rho }_{cv\mathbf{k}}}{\partial t}-{\epsilon }_{cv\mathbf{k}}{\rho }_{cv\mathbf{k}}+\sum _{{\mathbf{k}}^{\prime }}{V}_{cv}^{\mathbf{k},{\mathbf{k}}^{\prime }}({\rho }_{vv\mathbf{k}}-{\rho }_{cc\mathbf{k}}){\rho }_{cv{\mathbf{k}}^{\prime }}$$

$$\displaystyle\qquad-\sum_{\mathbf{k}^{\prime}}\big{[}V_{vv}^{\mathbf{k},\mathbf{k}^{\prime}}(\rho_{vv\mathbf{k}^{\prime}}-1)-V_{cc}^{\mathbf{k},\mathbf{k}^{\prime}}\rho_{cc\mathbf{k}^{\prime}}\big{]}\rho_{cv\mathbf{k}}$$

$$=(u_{cc\mathbf{k}}-u_{vv\mathbf{k}}){\rho }_{cv\mathbf{k}}+u_{cv\mathbf{k}}({\rho }_{vv\mathbf{k}}-{\rho }_{cc\mathbf{k}}),$$

$$\displaystyle i\hbar\frac{\partial{\rho_{cc\mathbf{k}}}}{\partial{t}}-\sum_{\mathbf{k}^{\prime}}\big{[}V_{vc}^{\mathbf{k},\mathbf{k}^{\prime}}\rho_{vc\mathbf{k}^{\prime}}\rho_{cv\mathbf{k}}-V_{cv}^{\mathbf{k},\mathbf{k}^{\prime}}\rho_{cv\mathbf{k}^{\prime}}\rho_{vc\mathbf{k}}\big{]}$$

$$=u_{cv\mathbf{k}}{\rho }_{vc\mathbf{k}}-u_{vc\mathbf{k}}{\rho }_{cv\mathbf{k}}.$$

$$V_{nm}^{\mathbf{k},{\mathbf{k}}^{\prime }}\approx \frac{e^2}{2{ϵ}_0}\frac{⟨n\mathbf{k}|n{\mathbf{k}}^{\prime }⟩⟨m{\mathbf{k}}^{\prime }|m\mathbf{k}⟩}{\left|\mathbf{k}-{\mathbf{k}}^{\prime }\right|({ϵ}_s+r_0\left|\mathbf{k}-{\mathbf{k}}^{\prime }\right|)},$$

$${\rho }_{cv\mathbf{k}}^{(\omega )}=\sum _n\frac{{\psi }_{\mathbf{k}}^{(n)}{U}_{n0}}{\hslash \omega -E_n},$$

$${\rho }_{cv\mathbf{k}}^{(2\omega )}=\sum _{n,m}\frac{{\psi }_{\mathbf{k}}^{(n)}{U}_{nm}{U}_{m0}}{(2\hslash \omega -E_n)(\hslash \omega -E_m)},$$

$${\rho }_{cv\mathbf{k}}^{(3\omega )}=\sum _{l,n,m}\frac{{\psi }_{\mathbf{k}}^{(l)}{U}_{ln}{U}_{nm}{U}_{m0}}{(3\hslash \omega -E_l)(2\hslash \omega -E_n)(\hslash \omega -E_m)},$$

$${\rho }_{cc\mathbf{k}}^{(\omega )}=0,$$

$${\rho }_{cc\mathbf{k}}^{(2\omega )}={\rho }_{cv\mathbf{k}}^{(\omega )}{\rho }_{vc\mathbf{k}}^{(\omega )},$$

$${\rho }_{cc\mathbf{k}}^{(3\omega )}={\rho }_{cv\mathbf{k}}^{(2\omega )}{\rho }_{vc\mathbf{k}}^{(\omega )}+{\rho }_{cv\mathbf{k}}^{(\omega )}{\rho }_{vc\mathbf{k}}^{(2\omega )}.$$

$$\displaystyle\varepsilon_{cv\mathbf{k}}\approx 2\varepsilon_{\kappa}\big{[}1-\zeta\tau\kappa(1-\eta_{\kappa}^{2})\sin(3\theta)\big{]}\,,$$

$$\displaystyle\frac{p_{vc\mathbf{k}}^{x}}{p_{0}}\approx\eta_{\kappa}\cos(\theta)+i\tau\sin(\theta)+\frac{\zeta\kappa}{2}\bigg{[}3i\cos(2\theta)-i\cos(4\theta)$$

$$\displaystyle\quad-\tau\eta_{\kappa}(4+\eta_{\kappa}^{2})\sin(2\theta)+\tau\eta_{\kappa}(2-\eta_{\kappa}^{2})\sin(4\theta)\bigg{]}\,,$$

$$\displaystyle\frac{p_{vc\mathbf{k}}^{y}}{p_{0}}\approx\eta_{\kappa}\sin(\theta)-i\tau\cos(\theta)-\frac{\zeta\kappa}{2}\bigg{[}3i\sin(2\theta)+i\sin(4\theta)$$

$$\displaystyle\quad+\tau\eta_{\kappa}(4+\eta_{\kappa}^{2})\cos(2\theta)+\tau\eta_{\kappa}(2-\eta_{\kappa}^{2})\cos(4\theta)\bigg{]}\,,$$

$$\displaystyle\frac{\Omega_{cc\mathbf{k}}^{x}-\Omega_{vv\mathbf{k}}^{x}}{a}\approx\eta_{\kappa}\bigg{[}\tau\frac{\sin(\theta)}{\kappa}+\frac{\zeta}{2}(5-\eta_{\kappa}^{2})\cos(2\theta)$$

$$\displaystyle\qquad\qquad\qquad-\frac{\zeta}{2}(3-\eta_{\kappa}^{2})\cos(4\theta)\bigg{]}\,,$$

$$\displaystyle\frac{\Omega_{cc\mathbf{k}}^{y}-\Omega_{vv\mathbf{k}}^{y}}{a}\approx-\eta_{\kappa}\bigg{[}\tau\frac{\cos(\theta)}{\kappa}+\frac{\zeta}{2}(5-\eta_{\kappa}^{2})\sin(2\theta)$$

$$\displaystyle\qquad\qquad\qquad+\frac{\zeta}{2}(3-\eta_{\kappa}^{2})\sin(4\theta)\bigg{]}\,,$$

$$E_n{\varphi }_l^{(n)}(\kappa )=2{\epsilon }_{\kappa }{\varphi }_l^{(n)}(\kappa )-\frac{e^2}{8(2\pi {)}^2{ϵ}_{0}a}{\int }_0^{\infty }\mathrm{d}{\kappa }^{\prime }$$

$$\displaystyle\qquad\times\phi_{l}^{(n)}(\kappa^{\prime})\kappa^{\prime}\Big{\{}2\sin(b_{\kappa})\sin(b_{\kappa^{\prime}})D_{l}(\kappa,\kappa^{\prime})$$

$$+[1+\cos (b_{\kappa })][1+\cos (b_{{\kappa }^{\prime }})]D_{l+\tau }(\kappa ,{\kappa }^{\prime })$$

$$\displaystyle\qquad+[1-\cos(b_{\kappa})][1-\cos(b_{\kappa^{\prime}})]D_{l-\tau}(\kappa,\kappa^{\prime})\Big{\}}\,.$$

Full text

††: \TDM

\ioptwocol

1 Introduction

The interest in the NLO response of 2D materials has recently grown drastically due to their large nonlinear coefficients, which facilitate a vast range of applications in optoelectronic devices [1, 2] as well as fundamental physics, e.g. nonlinear charge and spin Hall effects [3, 4, 5]. Various NLO processes in semi-metallic graphene have been extensively studied both theoretically and experimentally [6, 7, 8, 9, 10, 11]. In centrosymmetric graphene, even-order NLO processes vanish in the dipole approximation [7]. In contrast, the intrinsic broken inversion symmetry in monolayer TMDs leads to non-vanishing even-order NLO responses. Furthermore, strong excitonic effects due to reduced screening as well as huge spin-orbit coupling (SOC) in monolayer TMDs make their optical and electronic properties unique [12].

The linear optical response of monolayer TMDs has been extensively studied theoretically using the BSE in combination with single-particle band structures obtained by ab-initio techniques [13, 14, 15] or tight-binding (TB) methods [16, 17, 18]. Beside these rigorous calculations, analytical approaches have also been developed by employing the so-called massive Dirac Hamiltonian to obtain physical insight [16, 17, 19, 20]. Employing the analytical approach not only reduces the computational burden of excitonic optical response calculations but also leads to the identification of excitonic selection rules for the linear optical response in monolayer TMDs [21, 22]. Optical selection rules are fundamental principles, which determine allowed and forbidden transitions mediated by light. These rules are imposed by various symmetries of the system such as time-reversal, spatial or rotational symmetries. For instance, the brightness of $s$ and darkness of $p$ excitons in the linear optical spectrum of monolayer TMDs have been explained by selection rules [21, 22]. Note that although excitons in monolayer TMDs are of Wannier-Mott type, which are also found in conventional semiconductors, they possess different selection rules due to the chirality of the band structure at the bandgap [22]. Regarding the NLO response, a few calculations for monolayer TMDs can be found, in which the BSE solution is employed for computing the second harmonic susceptibilities [23, 24, 25]. While the linear optical response originates only from the transition between the ground state and excited states, the NLO responses depend also on the transition between excited states, i.e. inter-excitonic transitions. Therefore, more complicated selection rules are expected to govern NLO processes. The NLO selection rules have been studied qualitatively in Ref. [26] and quantitatively using a Wannier model (with $1/r$ Coulomb potential) in a recent work [27]. Nonetheless, the qualitative work can not provide any estimation for the excitonic matrix elements, while the Wannier model fails to accurately capture the nonlinear excitonic spectra due to the missing Bloch band curvature [22]. Therefore, a systematic study of excitonic selection rules and their effects on the NLO response including the Bloch band signature monolayer TMDs is missing.

In this work, we develop an analytical approach for systematic study of the NLO response in monolayer TMDs. By using an effective Hamiltonian, which includes the TW term, allowed and forbidden transitions between excitonic states are identified in the dipole approximation. A diagrammatic method is then suggested, as a means for predicting the dominant transitions and their origin for NLO processes up to any required order. We apply the proposed analytical approach to monolayer MoS2 as a case study and confirm the validity and accuracy of the analytical approach by comparing with the full BSE results. This paper is organized as follows. We begin by reviewing the theoretical framework for the NLO response in monolayer TMDs using the effective TW Hamiltonian. Based on the developed framework, optical selection rules for the first-, second- and third-order processes are discussed in section 3. In section 4, the linear and NLO conductivities computed for monolayer MoS2 are reported. Finally, a summary of the main findings is provided in section 5. A set of appendices explain the details of the derived expressions.

2 Theoretical framework

In this section, we present the theoretical framework for calculating the NLO response of monolayer TMDs using our analytical approach. First, general expressions for the optical conductivity (OC) and second/third harmonic generation (SHG/THG) are reviewed. Then, the TW Hamiltonian for monolayer TMDs is introduced. Finally, analytical expressions for the required excitonic matrix elements are presented. It should noted that throughout the text, all vectors and tensors are indicated by bold letters, and the single-particle/many-body operators and matrix elements are denoted by lowercase/uppercase letters, respectively.

2.1 Excitonic nonlinear optical response

The linear and NLO responses of a material can be characterized by linear and NLO conductivities (or equivalently susceptibilities) [28]. In the present work, we assume that the material is excited by a normal-incident monochromatic electric field written as $\boldsymbol{\mathcal{E}}(t)=\boldsymbol{\mathcal{E}}_{\omega}\exp(-i\omega t)+\mathrm{c.c.}$, where $\boldsymbol{\mathcal{E}}_{\omega}=\mathcal{E}_{\omega}^{x}\bm{e}_{x}+\mathcal{E}_{\omega}^{y}\bm{e}_{y}$ is the electric field phasor ($\bm{e}_{x}$ and $\bm{e}_{y}$ are unit vectors along the $x$- and $y$-directions, respectively). The light-matter interaction is then studied using a density matrix approach, which has been explained in detail in our previous works, c.f. Refs. [29, 30]. In this approach, a set of coupled equations for the reduced density matrix $\rho_{ij\mathbf{k}}=\langle{0}|{\hat{c}_{i\mathbf{k}}^{\dagger}\hat{c}_{j\mathbf{k}}}|{0}\rangle$ is derived, where $|{0}\rangle$ denotes the many-body ground state. Furthermore, $i$ and $j$ are band indices, $\mathbf{k}$ is the wavevector, and $\hat{c}^{\dagger}$ ($\hat{c}$) is the fermionic creation (annihilation) operator. Neglecting exchange terms, the dynamical equations for $\rho_{ij\mathbf{k}}$ for a two-band semiconductor within the mean-field approximation are presented in Eqs. (2cfhij) and (2cfhij) of A. These equations are then solved perturbatively up to any required order of the field at vanishing temperature, and the $N^{th}$-order density matrix elements $\rho_{ij\mathbf{k}}^{(N)}$ are determined. In the present work, we only focus on the NLO processes of second/third harmonic generation for simplicity, and the required density matrices are found in Eqs. (2cfhijla)-(2cfhijlf).

Upon determining the density matrix, the $N^{th}$ harmonic current density is evaluated as $\mathbf{J}^{(N)}(t)=\mathrm{Tr}[\hat{\rho}^{(N)}\hat{\mathbf{J}}]$, where the charge current density operator reads $\hat{\mathbf{J}}=-e\hat{\mathbf{\Pi}}/(mA)$ (for a crystal area $A$) [30]. Here, $\hat{\mathbf{\Pi}}$ is the Heisenberg momentum defined as $\hbar\hat{\mathbf{\Pi}}=im[\hat{\mathcal{H}}_{0},\hat{\mathbf{R}}]$, where $\hat{\mathcal{H}}_{0}$ and $\hat{\mathbf{R}}$ are the total unperturbed Hamiltonian and many-body position operators, respectively [30]. In turn, the $N^{th}$ harmonic generated current density reads

[TABLE]

where $\boldsymbol{\sigma}^{(N)}$ is the $N^{th}$-order conductivity tensor ($N\geq 1$). Note that linear and NLO susceptibilities, $\boldsymbol{\chi}^{(N)}$ are related to their corresponding conductivities by $\boldsymbol{\sigma}^{(N)}=-iN\omega\epsilon_{0}\boldsymbol{\chi}^{(N)}$. The expressions of the first-, second-, and dominant part of the third-order conductivity tensors at vanishing temperature are given by [30]

[TABLE]

Here, $C_{3}\equiv eC_{2}\equiv e^{2}C_{1}\equiv-e^{4}\hbar/(m^{2}A)$, $\mu,\eta,\alpha,\beta=\{x,y\}$ are the Cartesian coordinates, and $E_{n}$ is the exciton energy. In addition, the excitonic matrix elements are defined as $\Pi_{0n}^{\alpha}=(\Pi_{n0}^{\alpha})^{*}\equiv-iE_{n}X_{0n}^{\alpha}$, $\Pi_{nm}^{\alpha}\equiv i(E_{n}-E_{m})Q_{nm}^{\alpha}$, and

[TABLE]

where $\varepsilon_{cv\mathbf{k}}\equiv\varepsilon_{c\mathbf{k}}-\varepsilon_{v\mathbf{k}}$ (indices $c$ and $v$ imply conduction and valence bands, respectively), and $p_{vc\mathbf{k}}^{\alpha}\equiv\langle{v\mathbf{k}}|{\hat{p}_{\alpha}}|{c\mathbf{k}}\rangle$ are the transition energy and momentum matrix element between single-particle states, respectively. Moreover, $\psi_{\mathbf{k}}^{(n)}$ is the exciton wavefunction projected onto a band-to-band transition, i.e. the exciton state is expressed as $|{\psi^{(n)}}\rangle=\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{(n)}\hat{c}_{c\mathbf{k}}^{\dagger}\hat{c}_{v\mathbf{k}}|{0}\rangle$ [31] (the exciton center-of-mass motion is neglected due to the negligible photon momentum). The summations over $\mathbf{k}$ should be understood as integrals over the Brillouin zone (BZ), i.e. $(2\pi)^{2}\sum_{\mathbf{k}}\rightarrow A\int_{\mathrm{BZ}}\mathrm{d}^{2}{\mathbf{\mathbf{k}}}$. Finally, $[\psi_{\mathbf{k}}^{(n)}]_{;k_{\alpha}}\equiv\partial\psi_{\mathbf{k}}^{(n)}/\partial k_{\alpha}-i(\Omega_{cc\mathbf{k}}^{\alpha}-\Omega_{vv\mathbf{k}}^{\alpha})\psi_{\mathbf{k}}^{(n)}$ denotes the generalized derivative of the exciton wavefunction with the Berry connections $\Omega_{nm\mathbf{k}}^{\alpha}\equiv\langle{n\mathbf{k}}|{i\partial/\partial k_{\alpha}}|{m\mathbf{k}}\rangle$. The generalized derivative stems from the intraband part of the position operator [32, 33]. Equations (2a) and (2) have been reported in our previous work [30] while Eq. (2) can be obtained in a similar manner by extending the general approach to the third order in the electric field. Although these expressions are obtained for a two-band semiconductor, they can be used for multi-band semiconductors by redefining the $Q_{nm}^{\alpha}$ matrix elements (see Ref. [30] for more details). Note that, with a minor modification, these expressions can also be employed for calculating the nonlinear spin conductivities [5].

To determine excitonic energies and wavefunctions, we solve the BSE [31, 34], given by

[TABLE]

where $V_{cv}^{\mathbf{k},\mathbf{k}^{\prime}}$ is the direct Coulomb matrix elements defined in Eq. (2cfhijk) of A. The Coulomb matrix elements include band signatures through the Bloch overlaps, $\langle n\mathbf{k}|n\mathbf{k}^{\prime}\rangle$. For 2D materials, a screened Coulomb interaction in the Keldysh form is employed, which accurately describes the excitonic spectrum [17, 30, 35, 36]. The BSE solutions contain both the bound and unbounded excitons.

2.2 Trigonal warping Hamiltonian for monolayer TMDs

It is well-known that the valence and conduction bands of monolayer TMDs in the vicinity of the Dirac points are mainly formed by $(|d_{x^{2}-y^{2}}\rangle+|d_{xy}\rangle)/\sqrt{2}$ and $|d_{z^{2}}\rangle$ orbitals of the metal atom, respectively [19, 37, 38]. Hence, we can construct an effective 2-by-2 tight-binding (TB) Hamiltonian [5, 16], which reproduces the low-energy band structure of monolayer TMDs. Including the SOC, the TB Hamiltonian in $\mathbf{k}$-space reads [5, 39, 40]

[TABLE]

where $s=\pm 1$, $\Delta$, $\gamma$ and $\lambda$ are the spin index, on-site energy, effective hopping and SOC strength, respectively. The wavevector-dependent functions read

[TABLE]

with the lattice constant $a$. A typical band structure obtained for this Hamiltonian is illustrated in Fig. 1. The SOC lifts the degeneracy of the bands at each valley, while the time-reversal symmetry (TRS) relates the bands in separate valleys, i.e. spin-up states at the K valley are degenerate with spin-down states at the K*′* valley.

The massive Dirac Hamiltonian can be derived from Eq. (2ce) by a Taylor expansion of $f(\mathbf{k})$ and $g(\mathbf{k})$ around the Dirac points to linear order in $\mathbf{k}$. Although this approximate Hamiltonian has proven useful for characterizing the linear optical properties of TMD monolayers, it fails to account for any even-order nonlinear response (in the dipole approximation) due to the presence of full rotation symmetry, $C_{\infty}$ [5]. This issue can be resolved by including terms up to second order in $\mathbf{k}$ for $f(\mathbf{k})$, which leads to restored threefold rotation symmetry, $C_{3}$, and the TW of isoenergy contours [41]. Hence, we form the TW Hamiltonian by expanding $g(\mathbf{k})$ and $f(\mathbf{k})$ at the K/K*′* valleys as $g(\mathbf{k})\approx 3\sqrt{3}\tau$ and

[TABLE]

Here, $\kappa_{\alpha}\equiv a(k_{\alpha}-K_{\alpha})$ is the dimensionless wavevector measured with respect to K or K*′, $\zeta\equiv\sqrt{3}/12$ multiplies the trigonal warping term, and $\tau=\pm 1$ is the valley index. Note that the second-order term in $\kappa$ in the expansion of $g(\mathbf{k})$ is neglected here, since its contribution is much smaller than the equivalent term in $f(\mathbf{k})$ due to $\lambda\ll\gamma$. The band structure obtained from the TW Hamiltonian is compared with the TB band structure in Fig. 1. Although the TW band structure deviates from the TB band structure away from K or K′*, it accurately reproduces the optical response at low photon frequency [5]. Note that the band structure in Fig. 1 is slightly shifted compared to the typical band structure of monolayer TMDs, where the SOC mainly lifts the valence band degeneracy. However, this vertical shift of the band structure is not important in our calculations, since only energy differences (transitions) matter. Using the TW Hamiltonian, we derive analytical expressions for the momentum matrix elements and Berry connections to first order in $\zeta$ as presented in Eqs. (2cfhijlma)-(2cfhijlm) of B. These equations have previously been employed to obtain analytical expressions for the linear and NLO responses in the independent-particle approximation, which are in excellent agreement with the full TB results [5].

2.3 Analytical approach for excitonic optical response

The impact of excitons on the optical response can be accurately studied by incorporating the BSE solution. Although solving the BSE for a given single-particle band structure is possible, in principle, it is a numerically demanding practice. Furthermore, the physical interpretation of the optical response is not straightforward. Fortunately, the approximate solution of the BSE using the massive Dirac Hamiltonian provides not only accurate estimates of the excitonic energies and wavefunctions but also valuable physical insight [17, 42]. In the massive Dirac approximation, the BSE is transformed to a modified Schrödinger equation for a 2D Hydrogen-like atom [42] hereafter referred to as the analytical BSE, see Eq. (C) in C. Due to the full rotational symmetry, the solution of the analytical BSE yields an exciton wavefunction with a distinct angular momentum, given by $\psi_{\mathbf{k}}^{(n)}=\exp(il\theta)\phi_{l}^{(n)}(\kappa)$. Here, $\theta$ and $\kappa$ are the phase and magnitude of dimensionless wavevector, respectively, $l$ denotes the angular momentum of the exciton, and $\phi_{l}^{(n)}(\kappa)$ is the radial part of exciton wavefunction. Note that the angular momentum $l$ is defined in reciprocal space rather than real coordinate space.

In the BSE, the band overlaps $\langle n\mathbf{k}|n\mathbf{k}^{\prime}\rangle$, appearing in the integrand, are sensitive to the chosen phase of Bloch states. Although the choice of phase does not influence any measurable quantity such as the optical response, it can change the interpretation of angular momenta in the solutions of the analytical BSE. In analogy to the Hydrogen atom orbitals, Bloch phases can be chosen in such a way that the excitons with $l=-2,-1,0,+1,+2$ correspond to orbitals $d^{-},p^{-},s,p^{+},d^{+}$, respectively [21, 22, 43]. In the present work, however, we have chosen a different phase convention for the Bloch states, which leads to a $-1$ shift of the exciton angular momenta at the K valley, i.e. the $l=-3,-2,-1,0,+1$ states are labeled as $d^{-},p^{-},s,p^{+},d^{+}$ orbitals, respectively. At the K*′* valley, the angular momentum shift is $+1$, while the plus and minus orbitals interchange their roles due to the TRS, e.g. $p^{-}$ ($l=-2$) at the K valley is degenerate with $p^{+}$ ($l=+2$) at the K*′* valley. Transforming from our phase convention to the conventional Hydrogen atom phase convention is straightforward by relabeling excitons, yet it simplifies our analytical expressions considerably. Note that, analogous to the Hydrogen orbitals, excitons are labeled as $ns$, $np^{\pm}$, $nd^{\pm}$, etc, where $n$ is the principal quantum number starting from 1, 2 and 3 for $s$, $p^{\pm}$ and $d^{\pm}$ orbitals, respectively.

Besides the issue of exciton labeling, the calculated excitonic spectrum of the analytical BSE deviates from the conventional Rydberg series known for Hydrogen atoms due to the presence of Bloch band signatures and screened Coulomb interaction. In particular, the degeneracy between orbitals with the same principal quantum number is lifted. This has several implications for the energetic ordering of exciton states. Firstly, for the same principal quantum number, orbitals with lower angular momenta have higher energies (or smaller binding energies), e.g. $E_{3d^{\pm}}<E_{3p^{\pm}}<E_{3s}$ [44]. This can be explained by considering the screening property of the Keldysh potential. In fact, the 2D Coulomb interaction is only effectively screened for charges with small in-plane separation [35, 36]. The Keldysh potential has a weak logarithmic divergence at short distances, while it approaches the $1/r$ dependence of the conventional Coulomb potential at large distances [35]. For excitons with lower angular momenta, the wavefunction is confined closer to the origin, i.e. the hole and electron are closer to each other. Therefore, excitons with lower angular momenta experience less Coulomb interaction due to the pronounced short-range screening of the Keldysh potential and, hence, possess smaller binding energies. Secondly, orbitals with positive and negative angular momenta, e.g. $3p^{+}$ and $3p^{-}$, are not degenerate in a single valley, due to the non-zero Berry curvature [43]. Nonetheless, due to the TRS, orbitals with positive angular momenta at the K valley are degenerate with the negative ones at the K*′* valley and vice versa. Note that, if the TW Hamiltonian is used for the analytical BSE, a group of orbitals with angular momenta $l,l\pm 3,l\pm 6,\cdots$ are coupled [22], e.g. $s$ states are coupled to $f^{\pm}$ states. This coupling is typically weak ($\zeta<1$) and, hence, ignored here, which enables us to derive analytical expressions for the excitonic matrix elements.

Combining the analytical BSE solution with the matrix elements obtained using the TW Hamiltonian, i.e. Eqs. (2cfhijlma)-(2cfhijlm), both linear and NLO responses can be determined. This is done by performing the angular integral ($\theta$) for the exciton matrix elements between the ground state $|{0}\rangle$ and the excited states $|{n}\rangle$, i.e. $X_{0n}^{\alpha}$, or between the different excited states, i.e. $Q_{nm}^{\alpha}$. Allowed and forbidden transitions are identified based on the properties of these matrix elements. The analytical expressions for $X_{0n}^{x}$ at the K valley for $s$, $p^{\pm}$ and $d^{\pm}$ excitons are given by

[TABLE]

Here, $X_{0}=mv_{F}\pi/(2a^{2})$ and the Fermi velocity is given by $v_{F}\equiv\sqrt{3}a\gamma/(2\hbar)$. The wavevector dependent functions $\eta_{\kappa}$ and $\varepsilon_{\kappa}$ are defined as $\eta_{\kappa}\equiv\Delta_{s\tau}/\varepsilon_{\kappa}$ and $2\varepsilon_{\kappa}\equiv\sqrt{4\Delta_{s\tau}+3\gamma^{2}\kappa^{2}}$ with $\Delta_{s\tau}\equiv\Delta+3\sqrt{3}\lambda s\tau$. Moreover, if the initial excitonic state possesses an angular momentum of $l$, $Q_{nm}^{x}$ reads

[TABLE]

Similar expressions can be obtained for the $y$-direction, which differ from the above equations by factors of $\pm i$. In addition, excitonic matrix elements at the K*′* valley are found simply by interchanging plus and minus orbitals. These expressions clearly exhibit the excitonic selection rules as discussed in the next section.

3 Optical selection rules

Incorporation of Eqs. (2cfha)-(2cfhe) and (2cfhi)-(2cfhid) into the conductivity tensors of Eqs. (2a)-(2) can be performed by utilizing the schematic diagram presented in Fig. 2(a). In this diagram, the black dots represent the ground state $|{0}\rangle$ and first few excited states $|{s}\rangle$, $|{p^{\pm}}\rangle$ and $|{d^{\pm}}\rangle$. The lines represent transitions that are allowed by the optical selection rules, whereas the missing links indicate the forbidden transitions, for which excitonic matrix elements vanish. This diagram is obtained for the K valley, but a similar one can readily be provided for the K*′* valley by interchanging the plus and minus states. From a physical point of view, allowed and forbidden transitions are identified by examining the conservation of the angular momentum [26], i.e. preserving the total angular momentum including the photon, lattice and exciton momenta during a light-matter interaction. If the initial and final states carry angular momenta of $l_{i}$ and $l_{f}$, respectively, the transition between them is allowed in the Dirac approximation if $l_{f}-l_{i}=\pm 1$ (solid lines), while it is allowed by the TW if $l_{f}-l_{i}\pm 3=\pm 1$ (dashed lines). For the TW-enabled transitions, any combination of plus and minus signs is allowed, leading to a total of four possibilities. With our phase convention for Bloch eigenstates, the angular momentum of the ground state should be set to 0. Note that a $\pm$3 unit of angular momentum can be generated by the honeycomb lattice due to the $C_{3}$ symmetry. Moreover, a linearly-polarized photon can transfer $\pm$1 unit of angular momentum, since it is a linear combination of left- and right-handed circularly-polarized light. Using the momentum conservation rule, the transitions between $|{0}\rangle$ and $|{p^{+}}\rangle$ states are forbidden at the K valley (allowed by the TW at the K*′* valley), whereas the $|{s}\rangle\leftrightarrow|{p^{+}}\rangle$ transitions are allowed in both valleys. Among all allowed transitions, the TW-enabled transitions are typically weaker than the rest, because they are scaled by $\zeta$. In addition, since $\eta_{\kappa}\approx 1$ for small $\kappa$, the matrix elements of $|{0}\rangle\leftrightarrow|{d^{+}}\rangle$ transitions are considerably smaller than $|{0}\rangle\leftrightarrow|{s}\rangle$ transitions, i.e. $|X_{0n}^{x}(d^{+})|<|X_{0n}^{x}(s)|$. Therefore, the coupling between $|{0}\rangle$ and $|{s}\rangle$ states is generally stronger than the remaining transitions, which is illustrated by a thicker line for $|{0}\rangle\leftrightarrow|{s}\rangle$ transitions in Fig. 2.

Using the diagram in Fig. 2(a), we can systematically identify the dominant transitions for linear and NLO processes. The many-body system is initially in its ground state with all electrons occupying the valence band. Depending on the number of photons $N$, involved in the NLO process, precisely $N$ intermediate states are visited before the system eventually return to the ground state Therefore, two, three and four intermediate states are involved for the first-, second-, and third-order processes, respectively, with the successive transitions given by: First order $|{0}\rangle\rightarrow|{m}\rangle\rightarrow|{0}\rangle$, Second order $|{0}\rangle\rightarrow|{m}\rangle\rightarrow|{n}\rangle\rightarrow|{0}\rangle$, Third order $|{0}\rangle\rightarrow|{m}\rangle\rightarrow|{n}\rangle\rightarrow|{l}\rangle\rightarrow|{0}\rangle$. In each case, $|{l}\rangle$, $|{n}\rangle$ and $|{m}\rangle$ can be the ground or any of the excited states. Hereafter, any set of possible transitions is referred to as a (transition) path.

For a first-order process (such as linear absorption), the dominant transition path is $|{0}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$, and hence, mainly $s$ excitons are bight in the OC spectrum [22]. In addition to $s$ excitons, $d^{+}$ and $p^{-}$ excitons can also be excited, but their brightness is much lower than the $s$ ones, see Eqs. (2cfha)-(2cfhe). Figure 2(c) shows the dominant transition paths for the second-order processes (such as SHG and optical rectification), i.e. S1: $|{0}\rangle\rightarrow|{p^{-}}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$ and S2: $|{0}\rangle\rightarrow|{s}\rangle\rightarrow|{p^{-}}\rangle\rightarrow|{0}\rangle$. Note that the S1 and S2 paths lead to resonances at different frequencies, since the system is relaxed through distinct states, i.e. $|{s}\rangle$ and $|{p^{-}}\rangle$ for S1 and S2, respectively. The quadratic process requires at least one transition, which is enabled by the TW and hence, any quadratic response in monolayer TMDs vanishes if the TW is not included [5]. Besides these two dominant paths, two other paths are conceivable for the second-order process, $|{0}\rangle\rightarrow|{d^{+}}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$ and $|{0}\rangle\rightarrow|{s}\rangle\rightarrow|{d^{+}}\rangle\rightarrow|{0}\rangle$. However, their contributions are considerably smaller then the S1 and S2 paths due to weak coupling between the ground state and $d$ excitons, as well as between the $s$ and $d^{+}$ excitons [see Eq. (2cfhib)]. Moving to the third-order processes (such as THG or optical Kerr), four states are involved and, hence, the number of possible paths grows considerably. However, we can still identify the dominant paths as shown with labels in Fig. 2(c). In contrast to the quadratic response, the third-order response does not necessitate the TW to survive. The present algorithm for identifying the dominant transitions can readily be extended to higher order NLO processes in monolayer TMDs.

Combining the above-mentioned selection rules with the information about exciton energies, we can predict the spectral position and nature of resonance peaks appearing in OC, SHG and THG spectra. Starting with the linear optical process, it is expected that the first strong resonances in the OC spectrum correspond to $1s$ and $2s$ states as shown schematically in Fig. 3. This becomes more complicated for the NLO processes since contributions of all intermediate transitions should be summed to obtain the full spectrum, e.g. all transitions between $s$ and $p^{-}$ excitons should be considered for determining the SHG spectrum. Nonetheless, spectrally-neighboring states dominate, due to the energy-dependent denominators of Eq. (2) or (2). Therefore, the lowest frequency resonances in the SHG spectrum are anticipated to be at frequencies coinciding with half of $1s$ and $2p^{-}$ energies, and they are dominated by the transitions indicated in Fig. 3. Note that probing $p$ excitons in the linear optical spectrum of monolayer TMD is challenging due to the weak coupling between $|{0}\rangle$ and $|{p^{-}}\rangle$, whereas it can be explored experimentally using the SHG response as already demonstrated in Ref. [44]. Proceeding to the THG, we expect that the lowest frequency resonance occurs at one-third of the $1s$ state energy. In addition, multiple transitions contribute simultaneously to the resonance as illustrated in Fig. 3. In the next section, we numerically confirm the validity of these predictions for monolayer MoS2.

4 Numerical results for monolayer MoS2

In this section, we apply the proposed approach to compute OC, SHG, and THG spectra of suspended monolayer MoS2 as a representative member of the TMD family. Nonetheless, the main findings hold true for other monolayer TMDs due to similar physics. The required parameters, $\Delta$, $\gamma$ and $\lambda$, are determined by fitting to the experimental data or calculated quasi-particle band structures and, here, are set to 1.25 eV, 1.51 eV, 7.2 meV, respectively [5]. A lattice constant of $a=3.18$ $\mathrm{\AA}$ and screening length of $r_{0}=44.3$ $\mathrm{\AA}$ are assumed [45]. Due to the symmetry of the honeycomb lattice, at normal incidence, it is sufficient to consider the diagonal components of the conductivity tensors, i.e. $\sigma_{xx}^{(1)}$, $\sigma_{xxx}^{(1)}$, and $\sigma_{xxxx}^{(1)}$ [33, 46]. The line-shape broadening is accounted for by adding a small phenomenological imaginary part, $i\Gamma$, to the frequency, i.e. $\omega\rightarrow\omega+i\Gamma$. We set $\Gamma=10$ meV for our calculations here. For the full BSE calculations, a dense $\mathbf{k}$-mesh with more than 17000 points is employed for getting converged excitonic energies.

We begin by comparing the exciton energies obtained using the full BSE with the ones computed by the analytical approach (in the massive Dirac approximation) in Table 1. In the analytical approach, Eq. (C) should be solved multiple times for different angular momenta and using two different bandgaps, $\Delta\pm 3\sqrt{3}\lambda=\{\Delta_{A},\Delta_{B}\}$, for the A and B excitons, whereas, in the full BSE method, solving Eq. (2cd) once with the TB Hamiltonian provides all required exciton energies. Nonetheless, the one-dimensional integral equation of the analytical BSE approach is solved with much less numerical effort compared to the 2D integral equation of full BSE method. Results show that the binding energy of 1s exciton is approximately $\sim$ 0.5 eV, which is in agreement with previously reported experimental [14, 47] and theoretical works [15, 45, 48, 49]. Excitonic energies obtained in the analytical approach are overestimated slightly by approximately 30 meV, which is due to the underestimation of the Coulomb potential. Although this minor difference can be compensated merely by a slight reduction of the screening parameter, we use the same screening for both approaches here. Precise identification of the angular momentum characteristics of excitons obtained using the full BSE requires investigating the spatial distribution of the exciton wavefunction [44], whereas it is automatically accomplished in the analytical approach.

Upon determining the exciton energies, the OC spectrum is computed using Eq. (2a) and shown in Fig. 4. The excitonic OC spectrum displays several distinct resonances inside the bandgap due to the fundamental and higher order excitons and agrees with the previous results for monolayer MoS2 [16, 17]. Apart from the minor energy shift, the analytical approach provides a very accurate replication of the BSE spectrum. Furthermore, we can readily identify the contributions of different excitons in the spectrum. The results show that the resonances in the OC spectrum mainly stem from the $s$ excitons and, hence, the first two resonances in the spectrum appear at the frequencies of the $1s$ and $2s$ states at each valley. Although transitions mediated by $p^{-}$ and $d^{+}$ excitons are allowed by the selection rules, they are barely observable in the spectrum due to their negligible contribution. Indeed, we numerically confirm that the magnitude of the matrix element between the ground state and $1s$ (or $2s$) exciton is at least 10 times larger than that of the ground state and the $2p^{-}$ or $3d^{+}$ exciton.

For the second-order processes, we limit our analysis to SHG, but similar results can be obtained for other quadratic processes such as optical rectification [5]. Figure 5 illustrates the SHG spectrum computed for monolayer MoS2 using the full BSE and analytical approaches, which shows a good agreement between the two methods. As expected for a typical excitonic SHG spectrum, excitons appear as strong resonances at frequencies below half the bandgap. The spectral positions of the fundamental peaks agree well with the experimental reports [50, 23]. Using the analytical approach, we identify and label the resonance peaks according to the responsible exciton. The contributions of the two dominant paths in Fig. 2(c) are distinguished. These two paths generate resonances with similar strength but at different frequencies, corresponding to the $s$ and $p^{-}$ exciton energies. We also numerically examine the magnitude of the matrix element between the $1s$ and $2p^{-}$ or $3d^{+}$ states, and confirm that the former is approximately 10 times larger the latter. Combining this fact with the weak coupling of the ground state and $d^{+}$ excitons leads to a negligible contribution of the transition paths via the $d^{+}$ excitons in the total SHG response. Note that, for circularly polarized light, the resonances corresponding to $s$ and $p$ excitons possess different valley dichroism due to their distinct angular momenta [5].

In Fig. 6, we compare the THG conductivity computed using the full BSE with the analytical approach. The THG results follow the main trends observed in the OC and SHG responses. However, multiple transition paths should be included in order to obtain an accurate spectrum. In contrast to the OC and SHG responses, where each resonance is mainly formed by a single transition path, the resonances in THG spectrum have contributions from several distinct paths. For instance, the fundamental peaks, i.e. A1s and B1s, originate from three paths; T1: $|{0}\rangle\rightarrow|{s}\rangle\rightarrow|{p^{-}}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$, T2: $|{0}\rangle\rightarrow|{s}\rangle\rightarrow|{p^{+}}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$, and T3: $|{0}\rangle\rightarrow|{d^{+}}\rangle\rightarrow|{p^{+}}\rangle\rightarrow|{s}\rangle\rightarrow|{0}\rangle$. Note that, although the T3 path includes a weak transition between the ground and $d^{+}$ states, the coupling between $d^{+}$ and $p^{+}$ states can be sufficiently large so as to compensate this weak transition and, hence, contribute significantly to the THG response. Indeed, our calculations show that the magnitude of the coupling between the $2p^{+}$ and $3d^{+}$ states is almost 3 times larger than the coupling between the $1s$ and $2p^{+}$ states.

5 Conclusion

In summary, we have elucidated an analytical approach for calculating the linear and NLO responses of intrinsic monolayer TMDs, including excitonic effects. Analytical expressions for excitonic matrix elements were derived and used to determine the NLO response of monolayer MoS2. With a much less numerical effort, the analytical approach generates almost identical spectra to that of the numerically demanding BSE, while it also provides a transparent picture of the excitonic transitions in the NLO response. Our approach can readily be extended to identify the main transitions for any other NLO process in monolayer TMDs and, hence, provides a simple but accurate tool for analyzing experimental NLO spectra of monolayer TMDs.

The authors thank F. Hipolito for helpful discussions throughout the project. This work was supported by the QUSCOPE center sponsored by the Villum Foundation and the CNG center under the Danish National Research Foundation (project DNRF103).

Appendix A Equation of motion and its perturbative solution

Here, we review the derivation of the dynamical equation for the density matrix and present its perturbative solution up to third order in the external field. We follow the procedure outlined in Refs. [29, 30]. For a two-band semiconductor, the dynamical equations for the population densities and coherences in the mean-field and Tamm-Dancoff approximations read [29, 30]

[TABLE]

Here, $u_{nm\mathbf{k}}\equiv\langle n\mathbf{k}|\hat{u}(t)|m\mathbf{k}\rangle$ is the single-particle matrix element of the light-matter interaction Hamiltonian $\hat{u}(t)$. In the dipole approximation (long wavelength regime), the single-electron interaction Hamiltonian reads $\hat{u}(t)\equiv e\hat{\mathbf{r}}\cdot\boldsymbol{\mathcal{E}}(t)$. In Eqs. (2cfhij) and (2cfhij), we ignore the exchange terms due to their negligible influence [15]. Moreover, $V_{nm}^{\mathbf{k},\mathbf{k}^{\prime}}$ are the direct Coulomb matrix elements, given by

[TABLE]

where a Keldysh potential with substrate screening of $\epsilon_{s}$ and screening length of $r_{0}$ is assumed [5]. The approximation in Eq. (2cfhijk) is mainly because of the neglected terms for non-zero reciprocal vectors [30]. The equation of motion for $\rho_{cv\mathbf{k}}$ is simply obtained by complex conjugating the first equation, since the density matrix operator is Hermitian. In addition, the population conservation, i.e. normalization of the density matrix, leads to $\rho_{vv\mathbf{k}}=1-\rho_{cc\mathbf{k}}$.

For a time-harmonic electric field, i.e. $\hat{u}(t)=\hat{u}_{\omega}\exp(-i\omega t)+\mathrm{h.c.}$, the dynamical equation for the density matrix can be solved perturbatively up to any required order of the field using the Green’s functions [29]. Here, we only focus on the expressions for second/third harmonic generation for brevity, but similar results can be readily obtained for the optical rectification and Kerr processes. Hence, the $N^{th}$-order density matrix reads $\rho_{ij\mathbf{k}}^{(N)}(t)=\rho_{ij\mathbf{k}}^{(N\omega)}\exp(-iN\omega t)+\rho_{ij\mathbf{k}}^{(-N\omega)}\exp(iN\omega t)$, where $\rho_{ij\mathbf{k}}^{(N\omega)}$ for $N=1,2,3$ are given by

[TABLE]

Here, $U_{0n}=U_{n0}^{*}\equiv\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{(n)}u_{vc\mathbf{k}}^{(\omega)}$ and $U_{nm}\equiv\sum_{\mathbf{k}}\psi_{\mathbf{k}}^{(n)*}\psi_{\mathbf{k}}^{(m)}[u_{cc\mathbf{k}}^{(\omega)}-u_{vv\mathbf{k}}^{(\omega)}]$ are the excitonic matrix elements of the interaction Hamiltonian, where $u_{nm\mathbf{k}}^{(\omega)}\equiv\langle{n\mathbf{k}}|{\hat{u}_{\omega}}|{m\mathbf{k}}\rangle$. Indeed, $U_{0n}$ corresponds to the coupling between the ground state $|{0}\rangle$ and excited state $|{n}\rangle$, whereas $U_{nm}$ is the inter-excitonic matrix element. Note that \rho_{vc\mathbf{k}}^{(N\omega)}=\Big{[}\rho_{cv\mathbf{k}}^{(-N\omega)}\Big{]}^{*} and $\rho_{vv\mathbf{k}}^{(N\omega)}=-\rho_{cc\mathbf{k}}^{(N\omega)}$ due to the Hermiticity and normalization of the density matrix, respectively. Upon obtaining the density matrix, the expectation value of the $N^{th}$-order current density is determined straightforwardly as $\mathbf{J}^{(N)}(t)=\mathrm{Tr}[\hat{\rho}^{(N)}\hat{\mathbf{J}}]$ leading to Eqs. (2a)-(2) for the conductivities.

Appendix B Trigonal-warping Hamiltonian

For the 2-by-2 TB Hamiltonian, Eq. (2ce), the eigenenergies for a given spin read $\varepsilon_{c\mathbf{k}}=-\varepsilon_{v\mathbf{k}}=\sqrt{\Delta_{s}^{2}+\gamma^{2}F^{2}}$ with $\Delta_{s}\equiv\Delta+s\lambda g(\mathbf{k})$. The corresponding eigenvectors are then given by $|c\mathbf{k}\rangle=[\cos(b/2),-\sin(b/2)e^{i\phi}]$ and $|v\mathbf{k}\rangle=[\sin(b/2),\cos(b/2)e^{i\phi}]$, where $F$ and $\phi$ are defined using $f(\mathbf{k})=F\exp(i\phi)$, and $\cos(b)\equiv\Delta_{s}/\varepsilon_{c\mathbf{k}}$. Furthermore, the momentum operator in the TB method is given as the $\mathbf{k}$-derivative of the Hamiltonian, $\hbar p^{\alpha}=m\partial{H_{0}}/\partial{k_{\alpha}}$. Hence, the momentum matrix elements can readily be determined using $p_{nm\mathbf{k}}^{\alpha}=\langle{n\mathbf{k}}|{p^{\alpha}}|{m\mathbf{k}}\rangle$, e.g. $\hbar p_{vc\mathbf{k}}^{\alpha}=-m\gamma[\cos(b)F_{\alpha}+iF\phi_{\alpha}]$ with $F_{\alpha}\equiv\partial{F}/\partial{k_{\alpha}}$ and $\phi_{\alpha}\equiv\partial{\phi}/\partial{k_{\alpha}}$.

Taylor-expanding $f(\mathbf{k})$ and $g(\mathbf{k})$ functions in the vicinity of the Dirac points, K: $2\pi[3^{-1/2},3^{-1}]/a$ and K*′*: $2\pi[3^{-1/2},-3^{-1}]/a$, the massive Dirac and TW Hamiltonian are determined. Using Eq.(2cfg) for $f(\mathbf{k})$, we calculate any required parameter and keep only 0th and 1st order terms in $\zeta$. Hence, the transition energies, $\varepsilon_{cv\mathbf{k}}$, momentum matrix elements, $p_{vc\mathbf{k}}^{\alpha}$, and Berry connections, $\Omega_{cc\mathbf{k}}^{\alpha}-\Omega_{vv\mathbf{k}}^{\alpha}$, are given by

[TABLE]

with $p_{0}\equiv-mv_{F}$ . In these expressions, the terms proportional to $\zeta$ stem from the TW effect, and well-known expressions in the massive Dirac approximation are obtained by setting $\zeta=0$.

Appendix C Analytical BSE

In the massive Dirac approximation, a simple expression can be obtained for the Coulomb matrix element. Then, using an ansatz $\psi_{\mathbf{k}}^{(n)}=\exp(il\theta)\phi_{l}^{(n)}(\kappa)$ and performing the angular integral, the BSE is reduced to the following 1D eigenvalue problem for the radial part of exciton wavefunction [17]:

[TABLE]

Here, $\cos(b_{\kappa})\equiv\Delta_{s\tau}/\varepsilon_{\kappa}$ and $D_{l}(\kappa,\kappa^{\prime})$ is defined as

[TABLE]

with $q\equiv\sqrt{\kappa^{2}+\kappa^{\prime 2}-2\kappa\kappa^{\prime}\cos(\theta)}$. From Eq. (C), one can confirm that the $l$ state at one valley is degenerate with the $-l$ state at the other valley [note that $D_{l}(\kappa,\kappa^{\prime})=D_{-l}(\kappa,\kappa^{\prime})$]. For each angular momentum, Eq. (C) can be solved numerically by discretizing $\kappa$, using the technique explained in Ref. [51]. Upon determining the exciton eigenenergies and eigenvectors, the excitonic matrix elements are obtained using Eqs. (2cfha)-(2cfhe) and (2cfhi)-(2cfhid).