Berry curvature contribution towards $1s-2p_{\pm}$ interlayer exciton ultrafast transition within a $R-WSe_{2}/MoSe_{2}$ heterobilayer.

In addition to the intrinsic properties of the excitons in heterobilayers, the latter are also affected by a possible stacking mismatch of the TMD layers. In our case, when $WSe_{2}$ and $MoSe_{2}$ are stacked at a weak twist angle $\theta_{t}$, a valley mismatch arises due to the rotational mismatch of the associated Brillouin zones. This mismatch is characterized by the vector $\bm{\Delta K}=\textbf{K}_{Mo}-\textbf{K}_{W}$, which measures the reciprocal-space distance between equivalent $K$ points associated with the two layers. Additionally, moiré patterns emerge, exhibiting a long moiré period $L_{m}>100$ nm that characterizes a new unit cell much larger than those in the individual layers. In reciprocal space, this larger unit cell yields a mini Brillouin zone (mBZ) featuring novel high-symmetry points denoted as $K_{m}$, $\Gamma_{m}$ and $K^{\prime}_{m}$, as depicted in figure 2(a). For a weak twist angle, the moiré period is approximated as

$$L_{m}=\frac{a_{0<}}{\sqrt{\theta^{2}_{t}+\delta_{0}^{2}}},$$

where $\delta_{0}=1-a_{0<}/a_{0>}$ is lattice mismatch and $a_{0<(>)}$ is the inferior (superior) lattice constant of the TMD monolayers. The moiré superlattice is composed of three distinct atomic registries denoted by $R^{h}_{h}$, $R^{X}_{h}$, and $R^{M}_{h}$, interconnected by the $C_{3}$ group symmetry [24]. Here, $R^{h}_{h}$ represents a hexagonal registry where the transition metal atoms of the $WSe_{2}$ layer align directly on top of the transition metal atoms of the $MoSe_{2}$ layer. In contrast, $R^{X}_{h}$ and $R^{M}_{h}$ stackings are characterized by a lateral shift between the two layers, where the transition metal atoms of one layer align directly with the chalcogen atoms of the other layer. These two atomic configurations exhibit $C^{1}_{3}$ and $C^{2}_{3}$ rotation symmetries with respect to the $R^{h}_{h}$ stacking, respectively [see Fig. 2(b)].

This moiré network can be described by a periodic moiré potential [28, 24] $\Gamma(\textbf{r})=2V_{0}\sum_{j=1,2,3}\exp(i\textbf{G}_{j}\textbf{r})$ that respects the underlying $C_{3}$ symmetry of the moiré lattice [see figures 2(b) and (c)]. Here $V_{0}\sim 9$ meV is the amplitude of the modulation the value of which has been determined by DFT calculations [24]. Notice that this potential is not an electrostatic one that acts differently on the two types of charges but it represents a spatial modulation of the minimal gap between conduction and valence band. Both electrons and holes are therefore attracted by the minima of this potential so that it can be directly viewed as an exciton potential. The reciprocal lattice vectors $\textbf{G}_{j}$ are associated with the moiré supercell and can be obtained from those, $\textbf{b}_{j}$, of the (untwisted) monolayer and the twist angle $\theta_{t}$ between the two layers, $\textbf{G}_{j}=\theta_{t}\textbf{b}_{j}\times\textbf{e}_{z}$. The $\textbf{G}_{j}$ vectors are given by $\textbf{G}_{1}=(4\pi/\sqrt{3}L_{m})(1,0)$, $\textbf{G}_{2}=C^{1}_{3}\textbf{G}_{1}$ and $\textbf{G}_{3}=C^{2}_{3}\textbf{G}_{1}$, in terms of the rotations $C_{3}^{j}$ of the $C_{3}$ group. Near the minima of the moiré potential [as shown in the Figs. 2(b) and 2(c)], we might consider the latter to be parabolic within a series expansion, $V_{M}(\textbf{r}_{\alpha})=m_{\alpha}\omega_{\alpha}^{2}\textbf{r}_{\alpha}^{2}/2$ [24]. Here, the subscript $\alpha=e,h$ designates the electron band as ’e’ and the hole band as ’h’, respectively. The term $\omega_{\alpha}$ signifies the frequency of charge carriers, while their corresponding moiré potential size is denoted by $R_{\alpha}=\sqrt{\hbar/m_{\alpha}\omega_{\alpha}}$. The frequency of charge carriers is anticipated to be expressed in terms of the moiré period as: $\omega_{\alpha}=4\pi\sqrt{V_{0}/m_{\alpha}}/L_{m}$. Consequently, we derive $R_{\alpha}=\sqrt{\ell_{\alpha}L_{m}}$, where the characteristic length $\ell_{\alpha}$ is defined as $\ell_{\alpha}=\hbar/4\pi\sqrt{m_{\alpha}V_{0}}$. The estimated values for the characteristic lengths $\ell_{e}$ and $\ell_{h}$ are 2.58 Å and 3.45 Å, respectively, based on charge carriers’ masses $m_{e}=0.8$ and $m_{h}=0.45$. As we show below, this approximation is justified since the moiré period is much larger than our other characteristic length scales of the exciton, namely the Bohr radius $a_{B}=\hbar^{2}\epsilon/e^{2}\mu$ and the Compton length (the square root of the Berry curvature) $\lambda_{C}\sim\sqrt{\Omega}\sim\hbar/\sqrt{\mu\triangle}$, in terms of the reduced mass $\mu=m_{e}m_{h}/M=m_{e}m_{h}/(m_{e}+m_{h})$. Both are on the order of some nm.

The Coulomb-type electron-hole interactions in the van-der-Waals heterostructure need to take into account non-local screening due to possible excitations from the valence to the conduction band over the effective gap $\Delta$. Notice that in the MoSe₂/WeSe₂ heterobilayer this gap is the energy difference between the top of the WeSe₂ (hole) valence band and the bottom of the MoSe₂ (electron) conduction band. To leading order in $q$, these excitations yield a dielectric function

$$\epsilon(q)\simeq 1+\xi q,$$

(4)

where the characteristic length

$$\xi=\#\frac{e^{2}}{\epsilon\Delta}$$

(5)

is the polarizability of the 2D material embedded in a dielectric environment characterized by the dielectric constant $\epsilon$, where $\#$ is a numerical prefactor [37]. This polarizability yield precisely the screening length

$$\ell_{s}=2\pi\xi,$$

(6)

which enters into the expression of the Rytova-Keldysh potential [43]

$$V_{RK}(r=|{\bf r}|)=-\pi R_{y}\frac{a_{B}}{\ell_{s}}K_{0}\left(\dfrac{\sqrt{r^{2}+d^{2}}}{\ell_{s}}\right),$$

(7)

which is the 2D Fourier transform of the screened Coulomb interaction $V(q)=V_{0}(q)/\epsilon(q)$, where $V_{0}=2\pi e^{2}\exp(-qd)/q$ represents the bare Coulomb potential that takes into account the layer separation $d$, which is the minimal distance between the electron in the MoSe₂ layer and the hole in the WSe₂ layer. Moreover, $|{\bf r}|={\bf r}_{e}-{\bf r}_{h}$ is the relative 2D distance between the electron and the hole that constitute the (neutral) exciton, and $K_{0}(x)=H_{0}(x)-Y_{0}(x)$, in terms of the Struve and Bessel functions of the first kind, $H_{0}(x)$ and $Y_{0}(x)$, respectively. This potential incorporates the dielectric environment [24, 27, 12] via the average dielectric permittivity $\epsilon$ (due to the substrates) in two quantities: first, it occurs in the expression of the effective Rydberg energy $R_{y}=e^{2}/2\epsilon a_{B}$, in terms of the Bohr radius $a_{B}=\hbar^{2}\epsilon/e^{2}\mu$ and the reduced mass $\mu=m_{e}m_{h}/(m_{e}+m_{h})$, and second it enters in the screening length $\ell_{s}$ [see Eqs. (5) and (6)]. If we neglect for the moment corrections due to the Berry curvature that are discussed in the following section, the interlayer exciton Hamiltonian therefore consists of six terms [31],

$$H_{ix}=\frac{p^{2}}{2\mu}+V_{RK}(r)+\frac{m_{e}^{3}+m_{h}^{3}}{2M^{2}}\omega_{cm}^{2}r^{2}+\frac{P^{2}}{2M}+\frac{M\omega_{cm}^{2}}{2}R^{2}+|m_{h}-m_{e}|\omega_{cm}^{2}\textbf{R}\cdot\textbf{r},$$

(8)

where

	$\displaystyle{\bf r}={\bf r}_{e}-{\bf r}_{h}\qquad\text{and}\qquad\frac{{\bf p}}{\mu}=\frac{{\bf p}_{e}}{m_{e}}-\frac{{\bf p}_{h}}{m_{h}},$
	$\displaystyle{\bf R}=\sum_{\alpha=e,h}\frac{m_{\alpha}\textbf{r}_{\alpha}}{M}\qquad\text{and}\qquad{\bf P}=\sum_{\alpha=e,h}\textbf{p}_{\alpha}.$		(9)

are relative and center-of-mass coordinates, respectively.

In the absence of a moiré potential, the Hamiltonian would only consist of the first two terms that describe the relative motion of the electron-hole compound, to which we will add the corrective terms due to the Berry curvature in the next section. Furthermore, the fourth term represents the motion of the center of mass of the exciton. The plane-wave character is modified mainly by the fifth term, which has a tendency to localize the exciton in the minima of the moiré potential that are approximated by parabolic potentials (harmonic approximation). The latter approximation is clearly justified as long as the moiré lattice spacing $L_{m}$ is much larger than all other length scales, namely the effective size of the exciton that is given by the Bohr radius $a_{B}$ that is on the order of some nm.

Within the above decomposition, we have assumed that both the electron and the hole exhibit the same potential size, leading to $m_{e}\omega_{e}=m_{h}\omega_{h}$. Furthermore, we employ the approximation $\sum_{\alpha=e,h}m_{\alpha}\omega_{\alpha}^{2}=M\omega_{cm}^{2}$, which characterizes the center-of-mass localization in terms of energies. This yields $\sum_{\alpha=e,h}\hbar\omega_{\alpha}/R_{\alpha}^{2}=\hbar\omega_{cm}/R_{cm}^{2}$. Here, we define $\omega_{cm}$ and $R_{cm}=\sqrt{\hbar/M\omega_{cm}}$ as the interlayer exciton’s center-of-mass frequency and localization radius, respectively. These components can be expressed in terms of the moiré period as $\omega_{cm}=4\pi\sqrt{2V_{0}/M}/L_{M}$ and $R_{cm}=\sqrt{\ell_{cm}L_{m}}$, with a characteristic length $\ell_{cm}=\hbar/4\pi\sqrt{2MV_{0}}\simeq 1.46$Å.

Notice finally that the moiré potential has also a small effect on the relative exciton motion, as it may be seen from the third term, which we explicitly take into account for a quantitative calculation of the exciton spectra for twisted bilayer samples. The last term indicates that relative and center-of-mass degrees of freedom are eventually coupled via the moiré potential. However, we anticipate that this term can be neglected in the case of a weak twist angle due to its proportionality to $\hbar\omega_{cm}(a_{B}/L_{m})$.

In the remainder of this paper, we are interested mainly in the spectrum of the interlayer excitons, taking into account the first two terms of Hamiltonian (8) plus the corrections to the relative motion due to the moiré potential (third) term, while the effect of the center-of-mass Hamiltonian is mainly to trap the excitons in the potential minima, and the last term (coupling between the relative and center-of-mass degrees of freedom) is neglected. Most saliently, this separation into relative and center-of-mass dynamics, which cannot be achieved within the more appropriate description of the exciton in terms of a Dirac exciton [14], remains possible if Berry-curvature corrections are taken into account. The latter take into account approximately the coupling between the electron and hole branches, which is natural in a description in terms of the Dirac equation and that is recalled in the following subsection.

As discussed in the preceding subsection, our focus shifts towards the role of Berry curvature in the IX relative motion part. The non-commutativity of the position operator, stemming from the presence of non-zero Berry curvature, prompts us to adopt a generalization of the Peierls substitution [19, 20, 21]. This substitution, defined by canonical coordinates, is expressed as follows: $\hat{\textbf{p}}_{\alpha}=\textbf{p}_{\alpha}$ and $\hat{\textbf{r}}_{\alpha}=\textbf{r}_{\alpha}+\bm{\Omega}_{\alpha}\times\textbf{p}_{\alpha}/2\hbar$. This change in coordinates results in non-commutative relations satisfying: $[\hat{\textbf{x}}_{\alpha},\hat{\textbf{y}}_{\alpha}]=i\bm{\Omega}_{\alpha}$. Consequently, new relative canonical variables are obtained: $\hat{\textbf{p}}=\textbf{p}$ and $\hat{\textbf{r}}=\textbf{r}+\left[\bm{\Omega}_{e}(\textbf{p}_{e})\times\textbf{p}_{e}-\bm{\Omega}_{h}(\textbf{p}_{h})\times\textbf{p}_{h}\right]/2\hbar$. To assess the impact of Berry curvature on the IX relative motion Hamiltonian which includes the moiré term effect, denoted as

$$\hat{H}^{ix}_{RM}(\hat{\textbf{r}},\hat{\textbf{p}})=\frac{\hat{\textbf{p}}^{2}}{2\mu}+V_{RK}(|\hat{\textbf{r}}|)+\frac{m_{e}^{3}+m_{h}^{3}}{2M^{2}}\omega_{cm}^{2}\hat{\textbf{r}}^{2},$$

(10)

we employ a second-order Taylor expansion of the charge carriers’ Berry curvature, aided by the Foldy-Wouthuysen transformation [56]. This yields the following modified IX relative motion Hamiltonian:

$$\hat{H}^{ix}_{RM}=\dfrac{p^{2}}{2\mu}+V_{RK}(\textbf{r})+\frac{m_{e}^{3}+m_{h}^{3}}{2M^{2}}\omega_{cm}^{2}r^{2}+H^{ix}_{B}+H^{ix}_{D},$$

(11)

where the two corrective terms are given by

$$H^{ix}_{B}=\frac{1}{2\hbar}\bm{\nabla}_{r}V_{q}(r)\left[\bm{\Omega}_{ix}(\textbf{p},\textbf{P})\times\textbf{p}+\bm{\beta}_{+}(\textbf{p},\textbf{P})\times\textbf{P}\right],$$

(12)

and

$$H^{ix}_{D}=\frac{|\Omega_{ix}|}{4}\bm{\nabla}^{2}_{r}V_{q}(r).$$

(13)

Here, the potential $V_{q}(r)=V_{RK}(r)+\frac{m_{e}^{3}+m_{h}^{3}}{2M^{2}}\omega_{cm}^{2}r^{2}$ takes into account both the effect of the direct electron-hole interaction given in terms of the Rytova-Keldysh potential (7) and the correction due to the moiré potential. We anticipate, here, that the second term is negligible with respect to the first one. Indeed, its relative weight may be estimated as

$$\frac{\hbar\omega_{cm}}{R_{y}}\sim\sqrt{\frac{V_{0}}{R_{y}}}\frac{a_{B}}{L_{M}},$$

(14)

which is $\sim 0.03$ for a Bohr radius in the nm range, as compared to a characteristic moiré spacing of $L_{m}\sim 10$ nm, and where we have used $V_{0}\sim 10$ meV and $R_{y}\sim 100$ meV. The Berry cuvature’s quantities $\bm{\Omega}_{ix}(\textbf{p},\textbf{P})$ and $\bm{\beta}_{+}(\textbf{p},\textbf{P})$ are given by:

	$\displaystyle\bm{\Omega}_{ix}(\textbf{p},\textbf{P})$	$\displaystyle=$	$\displaystyle\bm{\Omega}_{e}(\textbf{p}+\frac{m_{e}}{M}\textbf{P})-\bm{\Omega}_{h}(\textbf{p}-\frac{m_{h}}{M}\textbf{P})\quad\text{and}\quad$
	$\displaystyle\bm{\beta}_{+}(\textbf{p},\textbf{P})$	$\displaystyle=$	$\displaystyle\frac{m_{e}}{M}\bm{\Omega}_{e}(\textbf{p}+\frac{m_{e}}{M}\textbf{P})+\frac{m_{h}}{M}\bm{\Omega}_{h}(\textbf{p}-\frac{m_{h}}{M}\textbf{P}).$

The quantity $\Omega_{ix}$ appearing in Eq. (12) and Eq. (13) represents the interlayer exciton Berry curvature in the two-band model case and can be calculated analytically within the massive-Dirac-fermion model. Its derivation is presented in detail in Appendix A. The component $\bm{\beta}_{+}(\textbf{p},\textbf{P})$ in Eq. (12) is defined as the sum of two distinct Berry curvature vectors with different signs corresponding to the electron and hole contributions. We anticipate that this term will have a weaker effect compared to the term $\bm{\Omega}_{ix}(\textbf{p},\textbf{P})$ in Eq. (12), due to its smaller amplitude $\beta_{+}(0,0)\sim 0.2$ Ų compared to $\Omega_{ix}(0,0)\sim 10$ Ų. Therefore, in the remainder of this paper, we will only consider the term $\bm{\Omega}_{ix}(\textbf{p},\textbf{P})$ in Eq. (12).

The term $H^{ix}_{B}$ appearing in Eq. (12), labeled as the Berry correction, can be understood as an effective spin-orbit coupling. This component can lift the degeneracy of states with angular momentum different from zero. It plays the role of an effective Zeeman effect in reciprocal space. It is important to point out that the interlayer excitons in the present heterobilayers are formed in the vicinity of $K_{m}$ and $K_{m}^{{}^{\prime}}$ valleys, i.e. at points in reciprocal space where the respective Berry curvature reaches its maximal absolute values, whence the importance of the two corrective terms (12) and (13). The Berry correction term affects exciton spectra, as it couples to the gradient of Keldysh potential which is generated by the attractive interaction between the electron of $MoSe_{2}$ and the hole of $WSe_{2}$ forming the interlayer exciton states. This is a result of the Dirac character of the low-energy charge carriers in these materials, which are typically described by a massive 2D Dirac equation [12, 50, 51, 52].

In addition to the Berry correction, there is another term known as the Darwin term Eq. (13), where the Berry curvature couples with the Laplacian of the Rytova-Keldysh potential. This occurs in the relativistic treatment of the Dirac Hamiltonian of the two adjacent bands within the heterobilayer.

As shown in appendix A, the charge carriers’ Berry curvature present two components: a zero order term, and a second order component in terms of $\textbf{p}^{2}$ and $\textbf{P}^{2}$, which can be neglected [see Eq. (34) and Eq. (35)]. The resulting interlayer exciton Berry curvature reads

$$\bm{\Omega}_{ix}\simeq\bm{\Omega}^{(0)}_{ix}(0,0)=\dfrac{\hbar^{2}}{4}\left(\dfrac{1}{m_{e}\triangle_{Mo}}+\dfrac{1}{m_{h}\triangle_{W}}\right)\vec{e}_{z}.$$

(16)

To understand the impact of Berry curvature on the interlayer exciton spectrum energy, we have simplified the expressions of the Berry correction and Darwin correction terms as:

	$\displaystyle H^{ix}_{B}$	$\displaystyle=$	$\displaystyle\frac{i\pi R_{y}}{2\ell^{2}_{s}}a_{B}\Omega_{ix}\frac{1}{\sqrt{r^{2}+d^{2}}}K_{1}\left(\frac{\sqrt{r^{2}+d^{2}}}{\ell_{s}}\right)\frac{\partial}{\partial\theta}$		(17)
			$\displaystyle+i\hbar\omega_{cm}\frac{\alpha_{0}}{2}\frac{\Omega_{ix}}{R_{cm}^{2}}\frac{\partial}{\partial\theta},$

and

$$H^{ix}_{D}=\frac{\Omega_{ix}\pi R_{y}}{4\ell^{2}_{s}}\Xi(r)+\hbar\omega_{cm}\frac{\alpha_{0}}{4}\frac{\Omega_{ix}}{R_{cm}^{2}}.$$

(18)

Here, the quantity $\alpha_{0}=(m_{e}^{3}+m_{h}^{3})/M^{3}$ is related to the (less relevant) contribution of the exciton moiré potential to the corrective terms, and the term $\Xi(r)$ appearing in Eq. (18) is given by

	$\displaystyle\Xi(r)$	$\displaystyle=$	$\displaystyle\frac{a_{B}r^{2}}{\ell_{s}(r^{2}+d^{2})}K_{0}\left(\frac{\sqrt{r^{2}+d^{2}}}{\ell_{s}}\right)$		(19)
			$\displaystyle+\frac{a_{B}(d^{2}-r^{2})}{(r^{2}+d^{2})^{\frac{3}{2}}}K_{1}\left(\frac{\sqrt{r^{2}+d^{2}}}{\ell_{s}}\right).$

The function $K_{0}(x)$ has already been defined in the context of the Rytova-Keldysh potential (7), and we have used $K_{1}(x)=H_{1}(x)-Y_{1}(x)$, where $H_{1}(x)$ and $Y_{1}(x)$ are the first-order Struve and Bessel functions of the second kind, respectively.

The diagonalization of the Hamiltonian matrix, constructed from Eq. (11) with the Rytova-Keldysh potential Eq. (7) and with the corrective terms Eq. (17) and Eq. (18), in the 2D-hydrogenic basis state $|\varphi_{n,\ell}\rangle$, leads to the following eigenvalues $E_{\tilde{n},\tilde{\ell}}$ and their corresponding eigenvectors

$$|\psi_{\tilde{n},\tilde{\ell}}\rangle=\sum_{n,|\ell|<n}A^{\tilde{n},\tilde{\ell}}_{n,\ell}|\varphi_{n,\ell}\rangle.$$

(20)

Here, $A^{\tilde{n},\tilde{\ell}}_{n,\ell}$ are their coefficients obtained after numerical diagonalization. The corresponding basis state used in this work has the following expression

$$\varphi_{n,\ell}(r,\theta)=\dfrac{Y_{n,\ell}(r)e^{i\ell\theta}}{\sqrt{2\pi}},$$

(21)

where the radial part

$$Y_{n,\ell}(r)=\frac{C_{n,\ell}}{a_{B}}e^{-\frac{r\alpha_{n}}{2a_{B}}}\left(\alpha_{n}\frac{r}{a_{B}}\right)^{|\ell|}L_{n-|\ell|-1}^{2|\ell|}\left(\alpha_{n}\frac{r}{a_{B}}\right),$$

(22)

is given in terms of associated Laguerre polynomials $L^{a}_{n}(x)$, and the normalization constants $C_{n,\ell}$ read

$$C_{n,\ell}=\frac{4}{(2n-1)^{\frac{3}{2}}}\left(\frac{(n-|\ell|-1)!}{n+|\ell|-1)!}\right)^{\frac{1}{2}}.$$

(23)

The terms $n$ and $\ell$ are the principal quantum number and the angular momentum, respectively, with $n=1,2,3...$ and $-(n-1)\leq\ell\leq n-1$. The states are $(2n-1)$-fold degenerate, and we recall that they are currently labeled s for $\ell=0$, p for $\ell=\pm 1$ and d for $\ell=\pm 2$. The index numbers $\tilde{n}$ and $\tilde{\ell}$ labeling the eigenstates refer as well to the dominant contribution of the coefficients $A^{\tilde{n},\tilde{\ell}}_{\bar{n},\bar{\ell}}$.

To illustrate the present framework, we have identified several key figures that capture the main results. These will be presented in the following section.

Figure 3 shows the energy levels of the interlayer exciton Rydberg states 1s, 2p, 2s, and 3d, along with their corresponding eigenvectors $|\psi_{\tilde{n},\tilde{\ell}}\rangle$, within the $R-\text{MoSe}_{2}/\text{WSe}_{2}/\text{SiO}_{2}$ heterostructure. In the calculation of the spectra, we have considered perfect alignment between the layers so that the effect of the twist is negligible ($L_{m}\rightarrow\infty$ and thus $\omega_{cm}\rightarrow 0$). Furthermore, we have used physically relevant values for the interlayer distance $d=7$ Åand the dielectric environment $\epsilon=1.5$ corresponding to a SiO₂ substrate with a dielectric constant of $\epsilon_{\text{SiO}_{2}}=2.1$. The left column of Fig. 3 shows the energy levels without incorporating the Berry curvature terms, while the energies on the right hand side take into account the latter via the terms (17) and (18). Most saliently, we observe a splitting of the degenerate Rydberg states, such as the $2p_{\pm}$ and $3d_{\pm}$ dark states. We observe a splitting energy of $|\Delta E_{2p_{\pm}}|=|E_{2p_{+}}-E_{2p_{-}}|=5$ meV for interlayer exciton states, akin to the behavior of a spontaneous Zeeman effect, attributable to the Berry correction term (17). Notice that the Darwin term (18) has a much weaker effect on the exciton spectra. Indeed, it is extremely local – in the case of a pure Coulomb potential $\sim 1/r$, Eq. (13) yields a Dirac delta term, $\nabla_{r}^{2}V(r)\propto\delta(r)$ – and therefore mainly affects the $ns$ states, which have a non-zero amplitude at the origin. However, even in the case of the most prominent $1s$ state, its energy is only shifted by $1.5$ meV for the parameters used in our calculations, i.e. a shift of only $\sim 0.7\%$ as compared to its binding energy.

The relevant parameters used in our calculations are given in Table 1. Here, we also show the binding energies of the intralayer excitons for WeSe₂ and MoSe₂ monolayers, for comparison. The energy splitting between their $2p_{+}$ and $2p_{-}$ states is approximately 12 meV and 20 meV, respectively. The observed splitting energy for interlayer excitons is therefore lower compared to that for intralayer excitons, possibly due to the smaller spacing between the layers, leading to the excited states of interlayer excitons being closely situated, thereby creating a denser spectrum. These findings agree well with previously reported results [14, 23, 12, 22, 17].

Figures 4(a) and 4(b) illustrate the dependence of interlayer exciton energies on the average dielectric environment and spacing separation, respectively, displaying 144 eigenvalues. In panel (a), where we vary the dielectric constant $\epsilon$, we have used a fixed value of $d=7$ Å, while the variation of the interlayer distance $d$ in panel (b) is carried out for a fixed value of $\epsilon=1.5$.

They clearly demonstrate a significant decrease in binding energies with an increase in either $\epsilon$ or $d$, as it is expected from the dependence of the effective potential on these parameters. Notice that the original $1/r$ potential is replaced by $1/\sqrt{r^{2}+d^{2}}$ due to the layer separation and thus the separation between the electron in one layer and the hole in the other one. Furthermore, the binding energies are controlled by the Rytova-Keldysh potential (7), which scales as $\sim\ell_{s}/r\sim 1/\epsilon$ in the limit $\ell_{s}\rightarrow 0$ [i.e. $\epsilon\rightarrow\infty$, see Eqs. (5) and (6)]. The dependence of IX splitting energies $|E_{1s}-E_{2p_{\pm}}|$ on $\epsilon$ and $d$ is also depicted in figures 4(c) and 4(d), respectively. These figures reveal a considerable decrease in spacing energy, tending toward the range of Terahertz frequencies as $\epsilon$ or the spacing separation between the two layers increases. These findings may have significant implications for the development of high-performance ultrafast devices.

Figures 5(a) and 5(b) represent the dependency of $|\Delta E_{2p_{\pm}}|$, the interlayer exciton splitting energy, on spacing separation and average dielectric environment for various values of $\epsilon$ and $d$, respectively. Figure 5(a) shows a decrease in spacing energy with increasing $d$, as one may have anticipated from the scaling arguments mentioned above. Indeed, the splitting energy increases when both $d$ and $\epsilon$ approach weak values. As shown in Fig. 5(a), $|\Delta E_{2p_{\pm}}|$ reaches a maximum value of approximately $8$ meV when $\epsilon=1$ and $d$ is close to zero, consistent with findings reported by Tang et al. [22].

It is noteworthy that adjusting the spacing $d$ between the two monolayers of TMDs can be achieved by incorporating layers of hBN, a typical experimental technique. This interplay can significantly impact the interlayer band-gap, although it does not notably affect the effective masses. Conversely, analyzing $|\Delta E_{2p_{\pm}}|$ as a function of the average dielectric environment, taken from different values of interlayer distance, reveals an unexpected non-monotonic behavior as depicted in Fig. 5(b). We identify three distinct regimes for the splitting $|\Delta E_{2p_{\pm}}|$: a decreasing trend within $\epsilon$ $\in$ [1, 3] when $d=7$ Å, followed by an increasing trend within $\epsilon$ $\in$ [3, 12.8], and a subsequent decrease within $\epsilon$ $\in$ $[12.8,\infty]$. The latter is expected again based on the scaling arguments.

While we do not provided a full-fletched theoretical explanation of this non-monotonic behavior on the dielectric function, we would emphasize that it is likely due to the hidden dependence of the effective interaction potential on the dielectric constant, which enters as a global scaling factor only in the large-$\ell_{s}$ limit, as mentioned above. In an intermediate regime, one needs to take into account that the Rytova-Keldysh potential is itself a (non-monotonic) function of the screening length $\ell_{s}\sim\epsilon$. Notice finally that this non-monotonic behavior cannot be obtained when using pure Coulomb interaction, which yields $|\Delta E_{2p_{\pm}}|=\frac{64}{81}\frac{\Omega_{ix}}{a_{B}^{2}}R_{y}$, as reported in [12, 16]. Consequently, this suggests a decrease in splitting energy as the dielectric environment increases, owing to its proportionality with $1/\epsilon^{4}$.

In the previous section, we have not taken into account the role of the twist between the two layers on the exciton energies, i.e. we have considered a situation with a twist angle near zero degree. Within this case, the quadratic potential which incorporates the moiré traps effects has not a significant effect on IX energy spectra, where it shifts the ground state by approximately 1 meV. Thus, to assess the impact of moiré effect on IX energies, Fig. 6(a) shows the dependence of the first interlayer exciton Rydberg states on the twisting angle effects, which reveals a decreasing in IX binding energies as increasing the twist angle from 0° to 4° with more significant impact on the s-type Rydberg states.

In counterpart, Fig. 6(b) showcases the dependence of twisting effect on the $|\Delta E_{2p_{\pm}}|$ splitting energy. This result reveals a considerably increasing in the splitting energy between $2p_{+}$ and $2p_{-}$ Rydberg states as increasing the twist angle from 0° to 4°.

The Berry curvature does not only modify the exciton binding energy spectrum but also the corresponding eigenvectors ($\psi_{2p_{\pm}}$ and $\psi_{3d_{\pm}}$) by mixing the $2p_{x}$ and $2p_{y}$ degenerate states in the linear combination $|\psi_{2p_{\pm}}\rangle=\frac{1}{\sqrt{2}}[|\psi_{2p_{x}\rangle}\pm i|\psi_{2p_{y}\rangle}]$. The exciton $nd_{\pm}$ states are also known to exhibit dark properties in linear optics, necessitating visualization through methods such as Third Harmonic Generation (THG) [3]. The $|\psi_{3d_{\pm}}\rangle$ eigenstate will be also superposed on $|\psi_{3d_{xy}\rangle}$ and $|\psi_{3d_{x^{2}-y^{2}}}\rangle$ states. The $2p_{+}$ and $2p_{-}$ exciton states have the same real part sign as $2p_{x}$ and opposite imaginary part signs of $2p_{y}$ (refer to Fig. 3). As a consequence, the two new eigenstates will exhibit opposing rotations due to symmetry, presenting an opportunity to explore the $1s-2p_{\pm}$ transition through the utilization of circularly polarized optical absorption originating from the $1s$ ground state. This will be the focal point of our discussion in the subsequent subsection.

The significant contribution of the Berry curvature lies in its ability to leverage the lifting of the $2p_{\pm}$ degeneracy, offering two degrees of freedom generated by optical circular polarization. This opens up a novel avenue in quantum information and THz devices. Moreover, harnessing the $1s-2p_{\pm}$ interlayer exciton transition in TMD heterobilayers holds immense promise across various domains including valleytronics, and bio-sensing. These examples merely scratch the surface of the myriad potential applications stemming from the utilization of the $1s-2p_{\pm}$ interlayer exciton transition in TMDs heterobilayers.

Among the well-known techniques used to observe exciton dark states in TMD monolayers, pump-probe exciton spectroscopy stands out as a highly effective method in nonlinear optics, particularly for elucidating the brightening of exciton $np_{\pm}$ dark states [3, 6]. By pumping the exciton ground state with a near-infrared (NIR) frequency, the ground state can then be probed by a medium infrared (MIR) frequency. This enables the $1s-2p_{\pm}$ exciton transition which falls within the TeraHertz frequency range for interlayer excitons in TMD heterobilayers [6]. In this study, we aim to uncover the exciton transition between the state $1s-2p_{\pm}$ for both interlayer and intralayer excitons. We consider a circularly polarized electric field F in interaction with interlayer exciton ground state. We analyze the interaction Hamiltonian $H_{int}=e\hat{\textbf{r}}\cdot\textbf{F}$, where the exciton dipole moment operator is $\hat{\textbf{r}}=\textbf{r}+\bm{\Omega_{ix}}\times\textbf{p}/2\hbar$. The electromagnetic-field interaction can be re-expressed as $\hat{\textbf{r}}\cdot\textbf{F}=\hat{\textbf{r}}^{+}\cdot\textbf{F}^{+}+\hat{\textbf{r}}^{-}\cdot\textbf{F}^{-}$, where $\hat{\textbf{r}}^{\pm}=(\hat{\textbf{r}}_{x}\pm i\hat{\textbf{r}}_{y})/\sqrt{2}=\textbf{r}^{\pm}\pm\bm{\Omega}_{ix}\times\textbf{p}^{\pm}/2\hbar$ and $\textbf{F}^{\pm}=(\textbf{F}_{x}\pm i\textbf{F}_{y})/\sqrt{2}$ is the circular polarized electric field. The dipole moment then becomes

$$\hat{\textbf{r}}^{\pm}=\frac{e^{\pm i\theta}}{\sqrt{2}}\left[r\pm\frac{i\Omega_{ix}}{2}\frac{1}{r}\frac{\partial}{\partial\theta}\right].$$

(24)

To determine the exciton polarizability, we calculate the dipole moment matrix element, which read

$$\langle\psi_{\tilde{n}^{{}^{\prime}},\tilde{\ell}^{{}^{\prime}}}|\hat{\textbf{r}}^{\pm}|\psi_{\tilde{n},\tilde{\ell}}\rangle=\delta_{|\ell-\ell^{{}^{\prime}}|,1}\sum_{n,\ell,n^{{}^{\prime}},\ell^{{}^{\prime}}}A^{\tilde{n},\tilde{\ell}}_{n,\ell}A^{\tilde{n^{{}^{\prime}}},\tilde{\ell}}_{n^{{}^{\prime}},\ell^{{}^{\prime}}}\\ \times\int_{0}^{\infty}Y_{n^{\prime},\ell^{\prime}}(r)\left[\frac{r^{2}}{\sqrt{2}}\mp\dfrac{\ell\Omega_{ix}}{2\sqrt{2}}\right]Y_{n,\ell}(r)dr,$$

(25)

the dipole moment matrix element Eq. (25) only exists between states with different symmetries (such as $s$ and $p$ states) and its magnitude increases as the difference between $\tilde{n}$ and $\tilde{n^{\prime}}$ increases. Meanwhile, to understand the role of Berry curvature in the transition between the $1s$ and $2p_{\pm}$ states, it is important to note that around the $K_{m}$ valley point, the transition of the $1s$ to $2p_{\pm}$ interlayer exciton is exclusively influenced by coupling to $\sigma^{\pm}$ circularly polarized light. Conversely, the $K_{m}^{{}^{\prime}}$ valley, connected to the $K_{m}$ valley via time reversal symmetry (TRS), results in an equivalent transition of $1s$ to $2p_{\mp}$ interlayer exciton state, coupling solely to $\sigma^{\mp}$ polarized light. This relationship is illustrated in the accompanying figure 7(a).

The oscillator strength of the $1s-np_{\pm}$ transition is defined as [6]:

$$f_{1s-np_{\pm}}=\dfrac{2\mu}{\hbar^{2}}|E_{1s}-E_{np_{\pm}}|\langle\psi_{1s}|\hat{\textbf{r}}^{\pm}|\psi_{np_{\pm}}\rangle|^{2}$$

(26)

Using the classical sum-over-states method, the frequency-dependent polarizability due to transitions from the $1s$ exciton state to $np_{\pm}$ states can be expressed as [45]:

$$\alpha(\hbar\omega)=\dfrac{e^{2}\hbar^{2}}{\mu}\sum_{n,\ell}\dfrac{f_{1s-np_{\pm}}}{|E_{1s}-E_{np_{\pm}}|^{2}-(\hbar\omega)^{2}-i\hbar\omega\gamma}$$

(27)

where $\gamma$ is a damping term that controls the line width of the resonances and is equal to 5 meV [3]. The imaginary and real parts of the exciton polarizability (in atomic units) provide information about absorption and refractivity for each transition, respectively. We have determined the exciton polarizability for both interlayer and intralayer excitons, as shown in figure 7(b), and found that the peak absorption of the $1s-2p_{\pm}$ interlayer excitons is more significant than in the intralayer excitons case. The resonant peak of the real and imaginary part of polarizability corresponding to the $|E_{\tilde{1s}}-E_{\tilde{2p}_{\pm}}|$ energy showed a high peak absorption for interlayer excitons compared to intralayer excitons. This can be especially explained by the short energy separation between $\tilde{1s}$ and $\tilde{2p}_{\pm}$ states.

Following to our previous results, we will now thoroughly investigate the interlayer exciton two-level system coupled with photons field. Let us first consider the initial state, denoted as $|1s,N\rangle$, and the final state, denoted as $|2p_{\pm},N-1\rangle$. Here, the integer number $N$ is defined as $N=aa^{\dagger}$, where $a$ ($a^{\dagger}$) represents the annihilation (creation) photon operator. This system can be described by the following Hamiltonian [23]:

$$H_{i}=\begin{pmatrix}E^{1s}_{IX}&\dfrac{V_{int}}{2}\\ \dfrac{V_{int}}{2}&E^{1s}_{IX}-\Delta_{0}\\ \end{pmatrix}.$$

(28)

Here $\Delta_{0}=E_{ph}-|E_{IX}^{1s}-E^{2p_{\pm}}_{IX}|$ is the detuning energy. $E_{IX}^{1s(2p_{\pm})}=2\Delta-E_{b}^{1s(2p_{\pm})}$ is the interlayer exciton ground state ($2p_{\pm}$ state) energy. $E_{ph}=\sum\hbar\omega(N+\frac{1}{2})$ is the photon excitation energy, describes the pump process. The coupling term $V_{int}$ is associated with the probe process, which is proportional to $1s-2p_{\pm}$ interlayer dipole moment $\hat{\textbf{r}}_{1s-2p_{\pm}}$, througth the expression $V_{int}=\bm{\hat{r}}_{1s-2p_{\pm}}\cdot\textbf{F}^{\pm}$. By directly diagonalized the Hamiltonian shown in Eq. (28), we can obtain the eigenenergies:

$$E_{\pm}=E_{IX}^{1s}-\frac{\Delta_{0}}{2}\pm\frac{\sqrt{V_{int}^{2}+\Delta_{0}^{2}}}{2}.$$

(29)

The corresonding eigenvectors are given by:

$$|\psi_{\pm}\rangle=C_{\pm}|1s,N\rangle+C_{\mp}|2p_{\pm},N-1\rangle$$

(30)

Here, $C_{\pm}$ is the Hopfield coefficient satisfying $|C_{+}|^{2}+|C_{-}|^{2}=1$, with

$$|C_{\pm}|^{2}=\frac{1}{2}\left(1\pm\frac{\Delta_{0}}{\sqrt{\Delta_{0}^{2}+V_{int}^{2}}}\right).$$

(31)

As anticipated, Eq. (29) illustrates the quantum mechanical coupling between matter and photon-dressed states, manifesting as energy-level repulsion. This interaction leads to the emergence of upper and lower exciton bands. Figure 8(i) delineates these two repulsive branches, while their respective Hopfield coefficients are elucidated in Figure 8(ii). Notably, the unmistakable anticrossing between the bands is evident around $\Delta_{0}=0$ detuning energy, indicating an anticipated self-hybridization between the $1s$ and $2p_{+}$ states influenced by photon fields.

To assess the hybridization of the $1s$ and $2p_{+}$ orbitals under the influence of an optical Stark effect, we scrutinize the dynamics of the upper and lower exciton bands. Utilizing the time-dependent eigenvectors: $\psi_{\pm}(r,t)=e^{-i(E_{+}+N\hbar\omega)\frac{t}{\hbar}}C_{\pm}|1s,N\rangle\pm e^{-i(E_{-}+(N-1)\hbar\omega)\frac{t}{\hbar}}C_{\mp}|2p_{+},(N-1)\rangle$, we present in figures 8(a-h) the evolving density $|\psi_{\pm}(r,t)|^{2}$ under strong regime coupling at $\Delta_{0}=0$, observed at various time points. These figures distinctly illustrate the mixing of orbitals, showcasing ultrafast transitions occurring within the femtosecond range. The two-level dressed model employed in this study to describe the optical Stark effect holds substantial significance in optics, offering valuable insights into the interaction between the $1s-2p_{\pm}$ interlayer exciton transition and photon fields. Moreover, the Berry curvature assumes a pivotal role by providing two degrees of freedom can be generated by circular optical polarization.