Energy spectrum of two-dimensional excitons in a non-uniform dielectric medium

TL;DR

This paper reveals that excitonic energy levels in monolayer transition metal dichalcogenides follow a simple Rydberg-like energy ladder, with the energy spectrum accurately modeled by an analytical Schrödinger equation using a modified potential.

Contribution

It introduces a simple analytical model for excitonic energy levels in 2D materials, validated by experimental data, and estimates the quantum defect parameter for various monolayers.

Findings

Excitonic states follow a Rydberg-like energy ladder with a small quantum defect.

The model accurately reproduces experimental magneto-optical data.

The quantum defect parameter is close to zero for WSe2 and MoS2 monolayers.

Abstract

We demonstrate that, in monolayers (MLs) of semiconducting transition metal dichalcogenides, the $s$-type Rydberg series of excitonic states follows a simple energy ladder: $\epsilon_n=-Ry^*/(n+\delta)^2$, $n$=1,2,\ldots, in which $Ry^*$ is very close to the Rydberg energy scaled by the dielectric constant of the medium surrounding the ML and by the reduced effective electron-hole mass, whereas the ML polarizability is only accounted for by $\delta$. This is justified by the analysis of experimental data on excitonic resonances, as extracted from magneto-optical measurements of a high-quality WSe$_2$ ML encapsulated in hexagonal boron nitride (hBN), and well reproduced with an analytically solvable Schr\"odinger equation when approximating the electron-hole potential in the form of a modified Kratzer potential. Applying our convention to other, MoSe$_2$, WS$_2$, MoS$_2$ MLs encapsulated in hBN, we estimate an apparent magnitude of $\delta$ for each of the studied structures. Intriguingly, $\delta$ is found to be close to zero for WSe$_2$ as well as for MoS$_2$ monolayers, what implies that the energy ladder of excitonic states in these two-dimensional structures resembles that of Rydberg states of a three-dimensional hydrogen atom.

Tables

Table 1. Table 1: Series of parameters ( E b exp subscript superscript 𝐸 exp 𝑏 E^{\text{exp}}_{b} , Δ ​ E 2 ​ s − 1 ​ s exp Δ subscript superscript 𝐸 exp 2 𝑠 1 𝑠 \Delta E^{\text{exp}}_{2s-1s} , δ exp superscript 𝛿 exp \delta^{\text{exp}} , μ exp superscript 𝜇 exp \mu^{\text{exp}} ) obtained from the analysis of PL spectra shown in Fig. 4 , compared with results of DFT calculations ( μ DFT superscript 𝜇 DFT \mu^{\text{DFT}} ( m 0 ) subscript 𝑚 0 (m_{0}) )  Kormányos  et al. ( 2015 ) .

Monolayer	$E_b^{\text{exp}}$ (meV)	$\Delta E_{2s-1s}^{\text{exp}}$ (meV)	${\delta }^{\text{exp}}$	${\mu }^{\text{exp}}$ $\left(m_0\right)$	${\mu }^{\text{DFT}}$ $\left(m_0\right)$
WSe2	167	130	-0.083	0.21	0.16
MoSe2	216	153	0.174	0.44	0.27
WS2	174	141	-0.229	0.15	0.15
MoS2	217	168	-0.095	0.26	0.24

Table 2. Table 2: Series of fitting parameters ( E g model subscript superscript 𝐸 model 𝑔 E^{\text{model}}_{g} , E b model subscript superscript 𝐸 model 𝑏 E^{\text{model}}_{b} and δ model superscript 𝛿 model \delta^{\text{model}} ) obtained from the analysis of data available in the literature  chernikovSM ; gorycaSM compared with the reported ones ( E g data subscript superscript 𝐸 data 𝑔 E^{\text{data}}_{g} and E b data subscript superscript 𝐸 data 𝑏 E^{\text{data}}_{b} ).

Publication	S-TMD ML	top medium	bottom medium	$E_g^{\text{data}}$ (eV)	$E_g^{\text{model}}$ (eV)	$E_b^{\text{data}}$ (meV)	$E_b^{\text{model}}$ (meV)	${\delta }^{\text{model}}$
Ref. 20	MoS2	hBN	hBN	2.16	2.16	221	223	-0.063
Ref. 20	WS2	hBN	hBN	2.238	2.24	180	180	-0.11
Ref. 20	MoSe2	hBN	hBN	1.874	1.87	231	232	0.044
Ref. 21	WS2	air	Si/SiO2	2.41	2.42	320	328	1.58

Equations

$$E_{ns}=E_g-\frac{R{y}^{\ast }}{(n+\delta {)}^2},$$

$$U_{ext}(r)=-\frac{e^{2}}{r_{0}}\Big{[}\frac{r_{0}^{*}}{r}-\frac{g^{2}r_{0}^{*2}}{r^{2}}\Big{]},$$

$${ϵ}_n=-R{y}^{\ast }/(n+g\kappa -1/2{)}^2,$$

$$U_{app}(\xi)=-U_{0}\Big{\{}\Big{[}\frac{1}{\xi}-\frac{0.21}{\xi^{2}}\Big{]}\theta(\xi-\xi_{0})+v_{0}\,\theta(\xi_{0}-\xi)\Big{\}},$$

$${ϵ}_n=-134\text{meV}/(n-0.099{)}^2.$$

$$\Big{\{}-\frac{\hbar^{2}}{2\mu}\Big{[}\frac{\partial^{2}}{\partial r^{2}}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\varphi^{2}}\Big{]}+U_{ext}(r)-\epsilon\Big{\}}\psi(r,\varphi)=0,$$

$$\Big{\{}\frac{d^{2}}{d\xi^{2}}+\frac{1}{\xi}\frac{d}{d\xi}+\frac{-k^{2}\xi^{2}+\kappa^{2}\xi-g^{2}\kappa^{2}-m^{2}}{\xi^{2}}\Big{\}}\phi_{m}(\xi)=0$$

$${\varphi }_{n,m}(r)=\frac{{\beta }_{n,m}}{\sqrt{2n+2{\delta }_m-1}}\sqrt{\frac{(n-|m|-1)!}{\Gamma (n+|m|+2{\delta }_m)}}\times ({\beta }_{n,m}r{)}^M{e}^{-{\beta }_{n,m}r/2}{L}_{n-|m|-1}^{2M}({\beta }_{n,m}r),$$

$${ϵ}_{n,m}=-\frac{\mu e^4}{2{\hslash }^2{\epsilon }^2}\frac{1}{(n+{\delta }_m-\frac{1}{2}{)}^2}.$$

$${ϵ}_n=-\frac{\mu e^4}{2{\hslash }^2{\epsilon }^2}\frac{1}{(n+g\kappa -\frac{1}{2}{)}^2},$$

$$⟨r^2{⟩}_{n,m}=\frac{2}{({\beta }_{n,m}{)}^2}[3-3m^2+5n(n-1)-5{\delta }_m-6|m|{\delta }_m+2{\delta }_m(5n+{\delta }_m)].$$

$$⟨r^2{⟩}_{n,0}=\frac{2}{({\beta }_{n,0}{)}^2}(2g^2{\kappa }^2+10g\kappa n-5g\kappa +5n^2-5n+3).$$

$$⟨r^2{⟩}_{n,0}{|}_{g\kappa =1/2}=(a_0^{\ast }{)}^2{n}^2(5n^2+1)/2\approx 5(a_0^{\ast }{)}^2{n}^4/2,$$

$$\Big{\{}-\frac{\hbar^{2}}{2\mu}\frac{1}{r}\frac{d}{dr}\Big{[}r\frac{d}{dr}\Big{]}-\frac{\pi e^{2}}{2r_{0}}\Big{[}\text{H}_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}-Y_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}\Big{]}-\epsilon\Big{\}}\phi(r)=0,$$

$$\Big{\{}-b^{2}\frac{1}{\xi}\frac{d}{d\xi}\Big{[}\xi\frac{d}{d\xi}\Big{]}-\pi b\Big{[}\text{H}_{0}(\xi)-Y_{0}(\xi)\Big{]}-W\Big{\}}\phi(\xi)=0$$

$${ϵ}_n=R{y}^{\ast }{W}_n=-\frac{R{y}^{\ast }}{(\alpha n+\beta {)}^2}=-R{y}^{\ast }\frac{\gamma }{(n+\delta {)}^2},$$

$$H(\mu,\varepsilon,r_{0})=-\frac{\hbar^{2}}{2\mu}\Delta_{2D}-\frac{\pi e^{2}}{2r_{0}}\Big{[}\text{H}_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}-Y_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}\Big{]},$$

$$U_{app}(\xi)=\left\{\begin{array}[]{cc}-U_{0}\Big{[}\frac{1}{\xi}-\frac{0.21}{\xi^{2}}\Big{]},&\text{for}\quad\xi>\xi_{0};\\ -1.71134\,U_{0},&\text{for}\quad\xi<\xi_{0}.\end{array}\right.$$

$$\phi_{1}(\xi)\sim J_{0}\big{(}\kappa\sqrt{v_{0}-|\mathcal{E}|}\,\xi\big{)},$$

$$\phi_{2}(\xi)\sim e^{-k\xi}\xi^{g\kappa}\Psi\Big{(}-\frac{\kappa^{2}}{2k}+g\kappa+\frac{1}{2},1+2g\kappa;2k\xi\Big{)}$$

$$\Big{\{}z\frac{d^{2}}{dz^{2}}+(c-z)\frac{d}{dz}-a\Big{\}}\Psi(a,c;z)=0.$$

$$f_1(\mathcal{E})=-\sqrt{v_0-|\mathcal{E}|}\frac{J_1(2g^2\kappa \sqrt{v_0-|\mathcal{E}|})}{J_0(2g^2\kappa \sqrt{v_0-|\mathcal{E}|})},$$

$$f_{2}(\mathcal{E})=-\sqrt{|\mathcal{E}|}+\frac{1}{2g}+\big{[}\kappa-(2g\kappa+1)\sqrt{|\mathcal{E}|}\big{]}\times\frac{\Psi\Big{(}-\frac{\kappa}{2\sqrt{|\mathcal{E}|}}+g\kappa+\frac{3}{2},2+2g\kappa;4g^{2}\kappa\sqrt{|\mathcal{E}|}\Big{)}}{\Psi\Big{(}-\frac{\kappa}{2\sqrt{|\mathcal{E}|}}+g\kappa+\frac{1}{2},1+2g\kappa;4g^{2}\kappa\sqrt{|\mathcal{E}|}\Big{)}},$$

$${ϵ}_n=-\frac{134\text{meV}}{(n-0.099{)}^2}.$$

$$f_{3}(\mathcal{E})=-\sqrt{|\mathcal{E}|}+\frac{1}{2g}-\frac{\big{[}\kappa-(2g\kappa+1)\sqrt{|\mathcal{E}|}\big{]}}{1+2g\kappa}\times\frac{{}_{1}F_{1}\Big{(}-\frac{\kappa}{2\sqrt{|\mathcal{E}|}}+g\kappa+\frac{3}{2},2+2g\kappa;4g^{2}\kappa\sqrt{|\mathcal{E}|}\Big{)}}{{}_{1}F_{1}\Big{(}-\frac{\kappa}{2\sqrt{|\mathcal{E}|}}+g\kappa+\frac{1}{2},1+2g\kappa;4g^{2}\kappa\sqrt{|\mathcal{E}|}\Big{)}},$$

$$H=-\frac{\hbar^{2}}{2\mu}\Delta_{2D}-\frac{\pi e^{2}}{2r_{0}}\Big{[}\text{H}_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}-Y_{0}\Big{(}\frac{r\varepsilon}{r_{0}}\Big{)}\Big{]}.$$

$$T_{0}(\beta)=-\frac{\hbar^{2}}{2\mu}\int_{0}^{\infty}\!\!dr\,\psi_{0}(\mathbf{r})\frac{d}{dr}\Big{[}r\frac{d\psi_{0}(\mathbf{r})}{dr}\Big{]}=\frac{\hbar^{2}\beta^{2}}{8\mu}.$$

$$U(r)=-\frac{e^2}{2\pi \epsilon }{\int }_0^{2\pi }d\theta {\int }_0^{\infty }dk\frac{e^{ikr\cos \theta }}{1+k{r}_0^{\ast }},$$

$$U_0(\beta )=-\frac{e^2{\beta }^2}{2\pi \epsilon }{\int }_0^{\infty }\frac{dk}{1+k{r}_0^{\ast }}{\int }_0^{2\pi }d\theta {\int }_0^{\infty }drr{e}^{(ik\cos \theta -\beta )r}.$$

$$\epsilon(a)=\frac{e^{2}}{r_{0}}\Big{[}\frac{\hbar^{2}\varepsilon^{2}}{8\mu r_{0}e^{2}}a^{2}+af(a)\Big{]},$$

Full text

Coulomb interaction in a non-uniform dielectric medium  Rytova (1967); Keldysh (1979), is one of the central points in investigations of large classes of nanoscale materials, such as, for example, graphene Kotov et al. (2012); Faugeras et al. (2015) and other atomically thin crystals including their heterostructures Geim and Grigorieva (2013), as well as colloidal nanoplatelets Ithurria et al. (2011), and two-dimensional perovskites Stoumpos et al. (2016); Blancon et al. (2018). This problem has been, in recent years, particularly largely discussed in reference to a vast amount of investigations of excitons in monolayers (MLs) of semiconducting transition metal dichalcogenides (S-TMDs) Chernikov et al. (2014); Raja et al. (2017); Koperski et al. (2017); Wang et al. (2018); Trushin et al. (2018). Surprisingly, at first sight, the Rydberg series of $s$-type excitonic states in these archetypes of two-dimensional (2D) semiconductors, does not follow the model system of a 2D hydrogen atom MacDonald and Ritchie (1986); Koteles and Chi (1988); Christol et al. (1993), with its characteristic energy sequence, $\sim 1/(n-1/2)^{2}$, of states with a principal quantum number $n$. The main reason for that is a dielectric inhomogeneity of the 2D S-TMD structures, $i.e.$, MLs surrounded by (deposited/encapsulated on/in) alien dielectrics. At large electron-hole ($e$-$h$) distances, the Coulomb interaction scales with the dielectric response of the surrounding medium whereas it appears to be significantly weakened at short $e$-$h$ distances by the usually stronger dielectric screening in the 2D plane. A common approach to account for the excitonic spectra of S-TMD MLs refers to the numerical solutions of the Schrödinger equation, in which the $e$-$h$ attraction is approximated by the Rytova-Keldysh (RK) potential Rytova (1967); Keldysh (1979). The RK approach has been used to explain a number of excitonic features in S-TMD MLs Stier et al. (2018). However, it is only solvable numerically. A more phenomenological and intuitive approach, presented below, might be an optional solution to this problem.

In this Letter, we demonstrate that the energy spectrum, $\epsilon_{n}$ ($n$=1, 2,…), of Rydberg series of $s$-type excitonic states in S-TMD MLs may follow a simple energy ladder: $\epsilon_{n}=-Ry^{*}/(n+\delta)^{2}$. From magneto-optical investigations of a WSe2 ML encapsulated in hexagonal boron nitride (hBN), we accurately establish that $Ry^{*}$=140.5 meV and $\delta$=-0.083 in this particular S-TMD system. The $\epsilon_{n}$ spectrum, with $\delta\sim 0$, turns out to closely reflect the characteristic spectrum of a three-dimensional (3D) hydrogen atom. The $\epsilon_{n}=-Ry^{*}/(n+\delta)^{2}$ ansatz is well reproduced with an analytical theoretical approach in which the $e$-$h$ potential is assumed to have the form of a modified Kratzer potential Kratzer (1920). $Ry^{*}$ is identified with the effective (3D) Rydberg energy, $Ry\times\mu/(\varepsilon^{2}m_{0})$, scaled by the dielectric constant $\varepsilon$ of the surrounding hBN medium and the reduced $e$-$h$ mass $\mu=(m_{e}m_{h})/(m_{e}+m_{h})$, with $Ry$=13.6 eV; $m_{e}$ and $m_{h}$ are, correspondingly, the electron and hole effective masses, and $m_{0}$ is the free electron mass. Dispersion of $Ry^{*}$ and $\delta$ parameters in different studied samples, WSe2, MoSe2, MoS2, and WS2 MLs encapsulated in hBN, is discussed and the reduced $e$-$h$ masses in these ML structures are estimated.

To accurately determine the characteristic ladder of $s$-type excitonic resonances in the experiment, we profited of a particularly suitable for this purpose method of magneto-optical spectroscopy Akimoto and Hasegawa (1967); Potemski et al. (1991). The active part of the structure used for these experiments was a WSe2 ML embedded in between hBN layers. More details on samples’ preparation and on the experimental techniques can be found in the Supplemental Materials (SM). We measured the (circular) polarization resolved magneto-photoluminescence (PL) at low temperatures (4.2 K) and in magnetic fields up to 14 T, applied in the direction perpendicular to the monolayer plane. Here we focus on magneto-PL spectra of our WSe2 ML, observed in the spectral range from $\sim$1.7 to $\sim$1.9 eV. As shown in Fig. 1(a) and (b), these spectra are composed of up to five PL peaks, which are clearly resolved in the range of high magnetic fields. Following a number of previous investigations Manca et al. (2017); Chow et al. (2017); Chen et al. (2018a); Stier et al. (2018); Liu et al. (2019) on similar structures, the observed PL peaks are identified with a series of excitonic resonances forming the 1$s$, 2$s$, …, 5$s$ Rydberg series of the so-called A exciton Koperski et al. (2017); Wang et al. (2018). Each $ns$ PL peak, $n$=1, 2, …, 5, demonstrates the valley Zeeman effect. This is illustrated in Fig. 1(c) in which the energies of the $\sigma^{+}$ and $\sigma^{-}$—polarized PL peaks are plotted as a function of the magnetic field. In accordance with previous reports we extract $g$=-4.1 for valley $g$-factor of the 1$s$ resonance, but read a significantly stronger valley Zeeman effect for all excited states ($g\sim$-4.8). The later observation is intriguing and should be investigated in more details, which is, however, beyond the scope of the present paper. We have also carried out the magneto-PL experiments on MoS2 and WS2 monolayers encapsulated in hBN, but only the 1$s$ and 2$s$ resonances could be observed in these structures in the range of magnetic fields applied (see SM for details).

The magnetic field evolution of the mean energies of $\sigma^{+}$ and $\sigma^{-}$ PL peaks is illustrated in Figs 1(d) and (e). These energies, $E_{ns}$, are plotted as a function of the magnetic field $B$ in Fig. 1(d), and as a function of $B^{2}$ in Fig. 1(e), which illustrates the characteristic but distinct behavior of $ns$ states in the so-called low- and high-field regime MacDonald and Ritchie (1986); Potemski et al. (1991). The high field limit, for a given $ns$ resonance, appears when $l_{B}\ll r_{ns}$, or conversely when binding energy of the $ns$ state $E^{ns}_{b}\ll\hbar\omega^{*}_{c}/2$. Here $r_{ns}$ and $E^{ns}_{b}$ denote, correspondingly, the mean lateral extension and the binding energy $E^{ns}_{b}$=$E_{g}$-$E_{ns}$ of a given $ns$ state at $B$=0, $\hbar\omega^{*}_{c}=\hbar eB/\mu$, $l_{B}=\sqrt{\hbar/eB}$ is the magnetic length and other symbols have their conventional meaning. In the high-field limit, the energies of $E_{ns}$ resonances approach a linear dependence upon $B$, with a slope given by $(n-1/2)\hbar\omega^{*}_{c}$. In the low field limit ($l_{B}\gg r_{ns}$, $E^{ns}_{b}\gg\hbar\omega^{*}_{c}/2$), the $ns$ resonances display the diamagnetic shifts: $E_{ns}(B)=E_{ns}(B=0)+\sigma B^{2}$, where $\sigma=(er_{ns})^{2}/8\mu$ is the diamagnetic coefficient. The 1$s$ and 2$s$ resonances follow the low-field regime in the entire range of the magnetic fields investigated due to their small exciton’s radii and/or large binding energies, see Fig. 1(e). The high field regime is approached for the 5$s$ resonances with an approximate linear dependence of $E_{5s}$ with $B$, in the range above $\sim$8 T. This linear dependence, marked with a solid line in Fig. 1(d), displays a slope of 2.1 meV/T, which if compared to $(9/2)\hbar\omega^{*}_{c}$ dependence, provides an estimate of 0.25 $m_{0}$ for the reduced mass in the WSe2 ML. However, one may also argue that working with magnetic fields up to 14 T only, the high field limit is still barely developed even for the 5$s$ state. In this context, our estimation of the reduced effective mass should be seen as its upper bound and, in the following we assume $\mu$=0.2 $m_{0}$ for the WSe2 ML, following the results of experiments performed in fields up to 60 T Stier et al. (2018).

In the following we focus on the energy sequence $E_{ns}$ of 1$s$, 2$s$, …, 5$s$ excitonic resonances as they appear in the absence of magnetic field. As shown in Fig. 1(e), the apparent $E_{ns}$ values are accurately determined with linear extrapolations of $E_{ns}$ versus $B^{2}$ dependences to B=0. Next, we put forward a hypothesis that the energy sequence $E_{ns}$ obeys the following rule:

[TABLE]

where, at this point, $E_{g}$, $Ry^{*}$, and $\delta$ should be regarded as unknown adjustable parameters. To test the above formula against experimental data, we note that Eq. 1 implies that (for example) the ratio ($E_{3s}-E_{1s})/(E_{2s}-E_{1s})$ only depends on $\delta$, and, reading this ratio from the experiment, we extract $\delta$=-0.083. With this value we find (see Fig. 2) that our experimental $E_{ns}$ series perfectly matches Eq. 1, and, at the same time, we determine two other parameters, $E_{g}$=1.873 eV and $Ry^{*}$=140.5 meV (or conversely, exciton binding energy $E_{b}=E_{g}-E_{1s}=Ry^{*}/(1-0.083)^{2}$=167 meV). The above $E_{g}$ and $E_{b}$ values are in very good agreement with those already reported in the literature Stier et al. (2018). Relevant for our further analysis, is the observation that the derived value for $Ry^{*}$ coincides well with the effective (3D) Rydberg energy $Ry^{*}$=13.6 eV$\cdot\mu/(\varepsilon^{2}m_{0}$)=134.3 meV, scaled by the dielectric constant of the surrounding hBN material $\varepsilon=\varepsilon_{hBN}$=4.5 Geick et al. (1966) and the reduced effective mass $\mu$=0.2 $m_{0}$ Stier et al. (2018) of the WSe2 ML. Intriguingly, the extracted $\delta$-parameter is close to zero which implies that the $\epsilon_{n}$=$E_{ns}$-$E_{g}$ Rydberg series found in a 2D system resembles that of a 3D hydrogen atom (\epsilon_{n}$$\sim$$-1/n^{2}).

On the theoretical ground, the problem of excitonic spectrum in S-TMD MLs is commonly solved by invoking the Rytova-Keldysh potential Rytova (1967); Keldysh (1979) $U_{RK}(r)$ (see purple curve in Fig. 3) to account for a specific character of the $e$-$h$ attraction in these systems. At large $e$-$h$ distances $r$, $U_{RK}(r)$ coincides with a usual Coulomb potential $U_{RK}(r)\sim-e^{2}/\varepsilon r$ (see blue curve in Fig.3), which scales with the dielectric constant $\varepsilon$ of the material surrounding the monolayer. On the other hand, $U_{RK}(r)\sim\log(r\varepsilon/r_{0})$ when $r$ is small, what accounts for the effective dielectric screening length $r_{0}=2\pi\chi_{2D}$ in the system, where $\chi_{2D}$ is the 2D polarizability of S-TMD ML. Distinctly, the apparent excitonic spectra and the related exciton binding energies critically depend on the efficiency of dielectric screening of the electron-hole attraction in the medium surrounding the monolayer.

Whereas previous efforts have been largely focused on the numerical study of such problem, we show that our model provides the analytical solution, which is in a good agreement with the experimental results discussed above. We propose to replace $U_{RK}(r)$ with the approximate potential $U_{app}(r)$, taken in the form of piecewise function. Namely, the sub-function $U_{cor}(r)$ defines $U_{app}(r)$ at small distances $r$ (the core domain), while the external potential $U_{ext}(r)$ corresponds to $U_{app}(r)$ in the region outside of the core.

We choose the external potential in the form of the modified Kratzer potential Kratzer (1920) (given in CGS units)

[TABLE]

where $r_{0}^{*}=r_{0}/\varepsilon$ is the reduced screening length and $g$ is a tunable parameter. For the case of $g^{2}=0.21$, $U_{ext}(r)$ fits $U_{RK}(r)$ in the region $r>r_{min}=0.46$ $r_{0}^{*}$ with the relative deviation less than $5\%$. For the WSe2 ML encapsulated in hBN, the minimal distance $r_{min}=4.6\,\mbox{{\AA}}$ is comparable with the lattice constant $a=3.28\,\mbox{{\AA}}$ Kormányos et al. (2015) of WSe2 (see Fig 3 for comparison).

The Schrödinger equation with the Kratzer potential (2) is exactly solvable providing the excitonic spectrum of the $s$-type states (see SM for details):

[TABLE]

in which $\kappa^{2}=2r_{0}^{*}/a_{B}^{*}$ and $a_{B}^{*}=\hbar^{2}\varepsilon/\mu e^{2}$ is the effective Bohr radius. The effective Rydberg constant $Ry^{*}=e^{2}/2\varepsilon a_{B}^{*}$ sets the energy scale in the system, while $\delta=g\kappa-1/2$ defines the relative positions of the energy levels in the spectrum. Since $g\kappa\propto\sqrt{\mu r_{0}}/\varepsilon$, the parameter $\delta$ is system dependent and its value can be tuned, in particular, by modifying the dielectric constant $\varepsilon$. We note that for a given material, $Ry^{*}\propto 1/\varepsilon^{2}$ and $\delta+1/2\propto 1/\varepsilon$. Such scaling laws as well as the energy sequence (3) can be derived from numerical simulations with $U_{RK}(r)$ potential at relatively large $\varepsilon$ (see SM for details). Note that Eq. 3 is an analogous of our experimentally found relation given by Eq. 1.

In the following, we introduce $U_{cor}(r)$ which replaces the Kratzer potential at small distances $r$, comparable with the lattice constant of WSe2 in our particular case. We choose the constant attractive potential $U_{cor}(r)=V_{0}$. Below we demonstrate, that $U_{cor}(r)$ does not change $\propto(n+\delta)^{-2}$ behaviour of the spectrum and modifies only $\delta$ parameter.

We consider the Kratzer and constant potentials as external and core ones, respectively. We choose the parameter $g^{2}=0.21$ and the region of validity of the Kratzer potential up to its minimum $\xi_{0}=2g^{2}$, where $\xi=r/r_{0}^{*}$. The parameter $V_{0}$ of the core potential is chosen as an average value of $U_{RK}(\xi)$ in the domain $\xi\in[0,\xi_{0}]$: $V_{0}=2\xi_{0}^{-2}\int_{0}^{\xi_{0}}d\xi\xi U_{RK}(\xi)$. Finally the approximate potential is

[TABLE]

where $\theta(x)$ is the step-function and $v_{0}=1.71134$. Please, note that the truncated Kratzer potential is applicable only if the radius of the core potential $r_{cor}$ is less or comparable with the lattice constant $a$. We estimate that this requirement is well satisfied for all monolayers encapsulated in hBN, and in particular for our WSe2 structure, for which $r_{cor}=2g^{2}r_{0}^{*}=4.2\,\mbox{\AA}$ and $a=3.28\,\mbox{\AA}$. Considering the $s$-type excitonic states in this later system, we derive the following formula (a detailed description is given in SM)

[TABLE]

Both found values: $Ry^{*}=134\,\text{meV}$ and $\delta=-0.099$ match their experimentally obtained counterparts (with the aid of Eq. 1) 140.5 meV and -0.083, respectively.

The applicability range of the formula given by Eq. 1 can be also considered from a different angle, i.e., when it is directly compared/fitted to numerical solutions obtained within the Rytova-Keldysh formalism. As demonstrated in the SM, the validity range of Eq. 1 can be defined with respect to a single, dimensionless parameter of a monolayer structure: $b=a_{B}^{*}/r_{0}^{*}$, and we find that our simple approach is valid when $b>0.3$, and estimate that this condition is well satisfied for all monolayers encapsulated in hBN. Nevertheless, even if $b$ is as small as $b\approx 0.1$, what may correspond to the case a monolayer deposited on Si/SiO2 substrate, the spectrum given by Eq. 1 coincides with that derived with the Rytova-Keldysh potential within the accuracy of $5\%$.

The model proposed above accounts well for the experimental results obtained for the WSe2 monolayer and it is obviously interesting to test this model for other S-TMD materials, as well. Unfortunately, the observation of the rich Rydberg spectrum of excitonic states in S-TMD MLs seems to be, so far, uniquely reserved for WSe2 MLs. Nevertheless, for all other S-TMD MLs studied, $i.e.$, MoS2, WS2, and MoSe2 MLs encapsulated in hBN, we do experimentally observe the 2$s$ in addition to the 1$s$ excitonic resonance (PL and reflectance contrast spectra), see Fig. 4 and Fig. S7 in SM. The energy positions, $E_{1s}$ and $E_{2s}$, of, correspondingly, the 1$s$ and 2$s$ resonances (of A exciton) are directly read from the data shown in Fig. 4. Of interest is the energy difference ($E_{2s}$-$E_{1s}$)=$\Delta E^{\text{exp}}_{2s-1s}$ listed in Table 1, for all four MLs investigated. As shown in Fig. 4, the PL peaks associated with the excited excitonic states are followed by noticeable PL tails developed at higher energies. We believe that these PL tails penetrate above the band-gap energies which are, however, not spectacularly marked in the spectra. We note, however, that in the case of our exemplary WSe2 ML, the PL intensity at the band-gap energy (accurately estimated from magneto-PL data and marked with a pink arrow in Fig. 4) consists of 5$\%$ of the intensity of the 2$s$ exciton PL peak. Applying the same convention to all spectra presented in Fig. 4, we estimate the band gaps in the three other MLs, as illustrated with pink arrows in this figure. Most critical is estimation of the band gap in MoSe2 ML, which requires a convolution of the PL spectra due to an additional signal associated with the B-exciton resonance (see SM for details). With estimation of the band gap and reading the energies of 1$s$ excitonic resonances directly from the spectra (see Fig. 4), we extract exciton binding energies $E_{b}^{exp}=(E_{g}-E_{1s})$ and show their values in Table 1. Having estimated $\Delta E^{\text{exp}}_{2s-1s}$ and $E_{b}^{exp}$ parameters, and following our predictions that $E_{ns}=E_{g}-Ry^{*}/{(n+\delta)^{2}}$, where $Ry^{*}=Ry\times\mu/(\varepsilon_{hBN}^{2}m_{0})$, we derive the $\delta^{exp}$ and $\mu^{exp}$ parameters for all MLs studied, see Table 1. We found very good agreement between our estimations and results of DFT calculations Kormányos et al. (2015) for the reduced masses in WS2 and MoS2 MLs, while we note an apparent discrepancy for WSe2 and MoSe2 MLs. We also applied our model to estimate values of band gaps and binding energies to the experimental data available in the literature Chernikov et al. (2014); Goryca et al. (2019), which is discussed in SM.

Concluding, the presented experimental and theoretical study let us proposed that the $ns$ Rydberg series of excitonic states in S-TMD monolayers encapsulated in hBN follows a simple energy ladder: $\epsilon_{n}=-Ry^{*}/(n+\delta)^{2}$. $Ry^{*}=Ry\times\mu/(\varepsilon^{2}m_{0})$, where $Ry$ is the Rydberg energy, $\mu$ denotes the reduced $e$-$h$ mass, and $\varepsilon$ is the dielectric constant of the surrounding material. The dielectric polarizability $\chi_{2D}$ of a monolayer is only encoded in $\delta$ which, in the first approximation, is given by $\delta=0.21\kappa-1/2$, where $\kappa^{2}=2\mu r_{0}e^{2}/\hbar^{2}\varepsilon^{2}$ and $r_{0}=2\pi\chi_{2D}$ is a characteristic 2D screening length. Strikingly, $\delta$ is found to be close to zero for WSe2 (and MoS2) ML whose $\epsilon_{n}$ spectrum resembles that of a 3D hydrogen atom. The proposed model may be applicable to other Coulomb bound states ($e.g.$ donor and/or acceptor states), also to other systems, such as colloidal platelets Ithurria et al. (2011) or 2D perovskites Stoumpos et al. (2016). Finally, we note that our $\epsilon_{n}=-Ry^{*}/(n+\delta)^{2}$ solution coincides with that expected for a hypothetical hydrogen atom in fractional dimension $N$, ($N=2\delta+3$), which was indeed speculated Christol et al. (1993) to mimic the spectrum of Coulomb bound states in low-dimensional semiconductor structures.

The work has been supported by the EU Graphene Flagship project (no. 785219), the ATOMOPTO project (TEAM programme of the Foundation for Polish Science, co-financed by the EU within the ERD-Fund), the National Science Centre, Poland (grants no. 2013/10/M/ST3/00791, 2017/27/B/ST3/00205, and 2018/31/B/ST3/02111), and the Nanofab facility of the Institut Néel, CNRS UGA, and the LNCMI-CNRS, a member of the European Magnetic Field Laboratory (EMFL). M.B. acknowledges the financial support from the Ministry of Education, Youth and Sports of the Czech Republic under the project CEITEC 2020 (LQ1601). K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan, and the CREST (JPMJCR15F3), JST.

S1 Samples and experimental setups

The active parts of our samples consist of a monolayer (ML) of semiconducting transition metal dichalcogenides (S-TMD), $i.e.$ WSe2, MoS2, WS2, and MoSe2, which has been encapsulated in hexagonal boron nitride (hBN) and deposited on a bare Si substrate. They were fabricated by two-stage polydimethylsiloxane (PDMS)-based gomezSM mechanical exfoliation of S-TMD and hBN bulk crystals.

The encapsulating hBN layers were rather thick, the bottom hBN layers were 35-40 nm-thick in case of WS2 and WSe2 monolayers and 60-65 nm-thick in case of MoS2 and MoSe2 monolayers. The thickness of the top hBN was of about 20nm in all structures studied.

The $\mu$-photoluminescence ($\mu$-PL) and $\mu$-reflectance contrast ($\mu$-RC) experiments were performed using a $\lambda$=515 nm CW laser diode and a 100 W tungsten halogen lamp, respectively.

Micro-magneto-PL measurements were performed in the Faraday configuration using an optical-fiber-based insert placed in a superconducting magnetic coil producing magnetic fields up to 14 T. The sample was mounted on top of an $x-y-z$ piezo-stage kept in gaseous helium at $T$= 4.2 K. The excitation light was coupled to an optical fiber with a core of 5 $\mu$m diameter and focused on the sample by an aspheric lens (spot diameter around 1 $\mu$m). The signal was collected by the same lens, injected into a second optical fiber of 50 $\mu$m diameter, and analyzed by a 0.5 m long monochromator equipped with a CCD camera. A combination of a quarter wave plate and a polarizer are used to analyse the circular polarization of signals. The measurements were performed with a fixed circular polarization, whereas reversing the direction of magnetic field yields the information corresponding to the other polarization component due to time-reversal symmetry.

Investigations at zero magnetic field were carried out with the aid of a continuous flow cryostat mounted on $x-y$ motorized positioners. The sample was placed on a cold finger of the cryostat. The excitation light was focused by means of a 50x long-working distance objective with a 0.5 numerical aperture producing a spot of about 1 $\mu$m. The signal was collected via the same microscope objective, sent through a 0.5 m monochromator, and then detected by a CCD camera.

S2 Excitonic spectrum and eigenfunctions in the Kratzer potential

We solve two-dimensional (2D) Schrödinger equation with the Kratzer potential kratzerSM for wave-function $\psi(\mathbf{r})=\psi(r,\varphi)$

[TABLE]

in which $U_{ext}(r)=-e^{2}/r_{0}(r_{0}^{*}/r-g^{2}r_{0}^{*2}/r^{2})$ is the modified Kratzer potential. $r$ is in-plane electron-hole distance, $\mu$ denotes the reduced electron-hole mass, $\varepsilon$ represents the dielectric constant of the material surrounding the monolayer, $r_{0}^{*}=r_{0}/\varepsilon$ is the reduced screening length, and $g$ is a tunable parameter. Taking $\psi_{m}(r,\varphi)=e^{im\varphi}\phi_{m}(r)/\sqrt{2\pi}$ and introducing the new variable $\xi=r/r_{0}^{*}$, we obtain the equation

[TABLE]

with $k^{2}=-2\mu\epsilon r_{0}^{*2}/\hbar^{2}>0$ and $\kappa^{2}=2\mu r_{0}e^{2}/\hbar^{2}\varepsilon^{2}>0$. The solution to this eigenvalue problem is

[TABLE]

with $M=\sqrt{m^{2}+g^{2}\kappa^{2}}$, $\delta_{m}=M-|m|$ and $\beta_{n,m}=2\mu e^{2}/\hbar^{2}\varepsilon(n+\delta_{m}-1/2)$, respectively. $n$=1, 2… is a principal quantum number, $m$=0, $\pm 1$, $\pm 2$… is an angular momentum quantum number, and $L_{n}^{\alpha}(x)$ is the modified Laguerre polynomial. The energy spectrum for such system is described with:

[TABLE]

For $g$=0, our result coincides with 2D hydrogen model YangSM . In the case of $s$-type states ($n$=1, 2… and $m$=0), the excitonic spectrum simplifies to

[TABLE]

We mention the following consequences of this model: (i) the energy scale (prefactor in $\epsilon_{n}$) does not depend on the screening length $r_{0}$. It coincides with the Rydberg constant $Ry^{*}$ for an exciton with reduced mass $\mu$ in an environment with dielectric constant $\varepsilon$, $i.e.$ $Ry^{*}=\mu e^{4}/2\hbar^{2}\varepsilon^{2}$; (ii) the information about relative positions of the energy levels of the system is encoded in the denominators in Eq. 9 and 10; (iii) the Kratzer potential lifts the Coulomb degeneracy of the $s$- ($m$=0) and $p$-type ($m$=$\pm$1) states, as one can be noticed from Eq. (9); (iv) since $Ry^{*}\propto 1/\varepsilon^{2}$ and $\delta+1/2\propto 1/\varepsilon$, the energy ladder of the excitons can be progressively tuned by changing the dielectric constant $\varepsilon$ of the surrounding medium. Surprisingly, the results of numerical simulations performed in the Rytova-Keldysh potential rytovaSM ; keldyshSM , discussed in the next section, demonstrate the similar behavior. This fact can be interpreted as the indirect confirmation, that the Kratzer potential is a good approximation for the considered model.

Using the wave-functions obtained in Eq. (8), we calculate the mean value of $r^{2}$, which can be useful for analysis of diamagnetic shift of excitons, which reads

[TABLE]

For $s$-type states characterized by $m$=0, it takes the form

[TABLE]

Moreover, for the special case $g\kappa=1/2$, while Eq. 10 resembles the three-dimensional (3D) hydrogen model, the mean value of the $r^{2}$ parameter is given by

[TABLE]

where $a_{0}^{*}$=$\hbar^{2}\varepsilon/\mu e^{2}$ is the effective Bohr radius. It is interesting to note that the latter formula coincides with the mean value of $r^{2}$ for 3D hydrogen atom betheSM .

The eigenfunctions of $s$-states of the Schrödinger Hamiltonian with the Kratzer potential (8) tend to zero at $r\rightarrow 0$. This is the consequence of the repulsive part of the potential at short distances. Therefore, such solutions can not be good approximation for $s$-state exciton wave-functions at small distances, since the Rytova-Keldysh potential is attractive. In order to improve the current result, one needs to modify the Kratzer potential at small distances.

S3 Numerical analysis of the excitonic spectrum in the Rytova-Keldysh potential

We solve numerically the eigenvalue problem for 2D Schrödinger Hamiltonian with the Rytova-Keldysh potential rytovaSM ; keldyshSM . We analyse the scale laws for the spectrum both as a function of principal number $n$ and dielectric constant of surrounding medium $\varepsilon$. We start from the radial equation on wave function $\phi(r)$ for the $s$-type states characterized by zero angular momentum:

[TABLE]

where $\mathrm{H}_{0}(x)$ and $Y_{0}(x)$ are the zeroth order Struve and Neumann functions. Introducing new variables $\xi=r\varepsilon/r_{0}=r/r_{0}^{*}$ and $\epsilon=(\mu e^{4}/2\hbar^{2}\varepsilon^{2})W=WRy^{*}$, we rewrite the equation

[TABLE]

with $b=\hbar^{2}\varepsilon^{2}/(\mu e^{2}r_{0})=a_{B}^{*}/r_{0}^{*}$ – the ratio of the natural scales in the system. We derive the spectrum of this differential equation as a function of $\varepsilon$ and $n$.

According to our hypothesis described in the main text, the excitonic spectrum should be described with

[TABLE]

where $\gamma=\alpha^{-2}\simeq 1$, while $\beta$ (and hence $\delta=\beta/\alpha$) strongly depends on $\varepsilon$. Therefore, we estimate the linear behaviour of $(-W_{n})^{-1/2}$ with $n$.

The results of numerical simulations (performed in ”Mathematica”) for the case of WSe2 ($\mu$=0.2 $m_{0}$ stierSM , $r_{0}=45$ Å berkelbachSM ) are presented in Figs S1 and S2. Indeed, the linear growth of $(-W_{n})^{-1/2}$ as a function of $n$ for different values of $\varepsilon$ can be appreciated in Fig. S1. The Fig. S2 qualitatively confirms that $\alpha\simeq 1$. The precision of this result (relative deviation of the curve $(-W_{n})^{-1/2}-n$ from its average value) becomes higher for larger values of $\varepsilon$.

Let us unify the aforementioned result for the case of other S-TMD monolayers and determine the limits of the applicability of our model. First, we derive numerically the spectrum $W_{n}\,(n=1,2\dots 5)$ from Eq. 15 for different values of the parameter $b$. Then we fit the obtained data with the formula $\mathcal{W}_{n}=-\gamma/(n+\delta)^{2}$ and extract the parameters $\gamma$ and $\delta$ for each $b$.

In further, we focus mainly on the range of parameters $b\in[0.05,2.5]$. Such a domain contains all $b$, accessible in the experiment, for all S-TMD monolayers and/or average dielectric constant of surrounding medium (such as PDMS, sapphire, Si/SiO2 or hBN). The parameters $\gamma$ and $\delta$, as functions of $b$, are presented in Figs. S3 and S4 respectively.

Note that the values of $\gamma$ slightly deviate from $1$ at $b>0.3$ and approach the limit $\gamma_{\infty}=1$ at $b\rightarrow\infty$. Therefore, at the relatively large $b$ the energy ladder of our problem coincides with the experimentally observed spectrum of excitons $\epsilon_{n}\approx-Ry^{*}/(n+\delta)^{2}$.

The values of $\delta$ are represented by monotonically decreasing function, with the limit $\delta_{\infty}=-1/2$ at $b\rightarrow\infty$, which is nothing but the case of 2D hydrogen atom. It is interesting to note that for $b\approx 0.9$ the parameter $\delta$ becomes zero and the exciton spectrum reproduces 3D hydrogen atom energy ladder.

In order to define the limits of the applicability of our model, we calculate the relative deviations $[(\mathcal{W}_{n}-W_{n})/W_{n}]\times 100\%$ for $n=1,2,\dots 5$ as a function of parameter $b$. The corresponding plots are presented in Fig.S5. One can see that starting from $b>0.3$ all the relative deviations become smaller than $2\%$ and tend to zero value at $b\rightarrow\infty$. The variational procedure, described in the next chapter of Supplementary Materials, predicts the ground state energy of 2D exciton in the Rytova-Keldysh potential with the same precision. Therefore, in order to have the self-consistent picture of all our calculations we choose the “$2\%$ deviation rule” as a formal criterium, which defines the limits of the applicability of our model.

Finally, let us discuss the specific symmetry of the studied Hamiltonian

[TABLE]

where $\Delta_{2D}$ is a two-dimensional Laplacian. The Hamiltonian is a homogeneous function $H(\lambda\mu,\lambda\varepsilon,\lambda r_{0})=\lambda^{-1}H(\mu,\varepsilon,r_{0})$ of its parameters $\{\mu,\varepsilon,r_{0}\}$ for any $\lambda>0$. Therefore the spectrum of the problem should have the same property: $\epsilon_{n}(\lambda\mu,\lambda\varepsilon,\lambda r_{0})=\lambda^{-1}\epsilon_{n}(\mu,\varepsilon,r_{0})$. One can demonstrate that our expression for the energy ladder (Eq. 16) satisfies this scaling law too.

S4 Derivation of the excitonic spectrum in WSe2 monolayer encapsulated in hBN

We consider the solution of 2D Schrödinger equation in the potential, defined in the main text, which reads

[TABLE]

We restrict our consideration only to the $s$-type states characterized by zero angular momentum. The regular radial solution of the Schrödinger equation in the region $\xi<\xi_{0}$ is

[TABLE]

where $J_{0}(x)$ is the zero-order Bessel function of the first kind, $\mathcal{E}=\epsilon/U_{0}$ and $v_{0}=1.71134$. The solution for the region $\xi>\xi_{0}$ has the form

[TABLE]

with $g^{2}$=0.21, $k^{2}=-2\mu\epsilon r_{0}^{*2}/\hbar^{2}>0$ and $\Psi(a,c;z)$ is the Tricomi’s function, which is regular at $\xi\rightarrow\infty$ and solves the degenerate hypergeometric equation batemanSM

[TABLE]

Introducing the normalized logarithmic derivatives for both solutions $f_{n}(\mathcal{E})=\kappa^{-1}[d\ln\phi_{n}(\xi)/d\xi]_{\xi=\xi_{0}},\,n=1,2$

[TABLE]

[TABLE]

one derives the energy spectrum from continuity equation $f_{1}(\mathcal{E})=f_{2}(\mathcal{E})$. We solve the latter equation for WSe2 ML encapsulated in hBN with a set of parameters: $\mu$=0.2 $m_{0}$stierSM , $\varepsilon=4.5$ geickSM , and $r_{0}=45\,\mbox{{\AA}}$ berkelbachSM . The derivatives as a function of $|\mathcal{E}|$ are presented in Fig. S6 and the excitonic spectrum is defined by the intersection points of $f_{1}(\mathcal{E})$ and $f_{2}(\mathcal{E})$ curves.

The excitonic spectrum, obtained within the described above method, can be presented in the same form as for the Kratzer potential (compare with Eq. 5 in the main text) and is given by

[TABLE]

This result demonstrates the good coincidence with the excitonic spectrum reported in Ref. 7, with relative errors $8\%$ for $n=1$, $3.5\%$ for $n=2$ and less than $2\%$ for higher excited states. In order to check the precision of the graphical method, we applied it for the case, when the core potential is described by the Kratzer one with the same parameter $g^{2}=0.21$. For the Kratzer potential, graphical solution provides the excitonic spectrum, which follows Eq. 10 of the SM. In this case, the logarithmic derivative is

[TABLE]

where ${}_{1}F_{1}(a,c;z)$ is the confluent hypergeometric function of the first kind, and corresponds to the blue curve in Fig. S6.

One can mention, that the calculation with the potential $U_{app}(\xi)$ does not approximate perfectly the $1s$-exciton state. From the technical point of view such discrepancy can be the consequence of the modification of the Rytova-Keldysh potential at small distances. In order to check this hypothesis, we estimate the ground state energy of exciton using another method. Namely, we derive the ground state energy of S-TMD excitons using the Ritz variational procedure. We take the variational wave-function in the form $\psi_{0}(\mathbf{r})=\beta\exp(-\beta r/2)/\sqrt{2\pi}$ and evaluate the average of the Hamiltonian

[TABLE]

The kinetic energy can be calculated directly with

[TABLE]

To determine the potential energy, a few steps need to be performed. First, we present the Rytova-Keldysh potential in the integral form CudazzoSM

[TABLE]

Then, we substitute it in the expression for the average potential energy

[TABLE]

Consequently, after evaluation of integrals and adding the kinetic energy part, we determine formula

[TABLE]

where $a=\beta r^{*}_{0}$ is the dimensionless parameter and

[TABLE]

The minimum of $\epsilon(a)$ can be found straightforwardly (with the help of ”Mathematica”, for example). For the case of $\mu=0.2$ $m_{0}$ stierSM , we got the ground state energy of exciton, $\epsilon_{0}$=-157 meV, which is in good agreement with the results of numerical simulation reported in Ref. 7, obtained for the same values of parameters. Moreover, the relative deviation of the variational exciton ground-state energy, calculated for the different values of $\varepsilon$, deviates from the numerical results discussed in Ref. 7 less than $2\%$.

Surprisingly, the numerical simulations for WSe2 monolayer encapsulated in hBN with an effective mass $\mu$=0.21 $m_{0}$ are in better agreement with the experimentally obtained excitonic spectrum,$\epsilon_{n}=-140.5\,\mbox{meV}/(n-0.083)^{2}$ (see the main text), than for $\mu$=0.20 $m_{0}$ and leads to the energy ladder of excitons given by

[TABLE]

The excitonic binding energy, $E_{b}$=162 meV, calculated with the aid of the Ritz variational method also nicely matches to the experimental one, $E_{b}$=167 meV, obtained in the main text.

S5 Magneto-photoluminescence investigation of MoS2 and WS2 monolayers

To study the ladder of excitonic states in S-TMD monolayer, we additionally performed the helicity-resolved magneto-photoluminescence experiments on MoS2 and WS2 monolayers encapsulated in hBN, see Fig. S7. For these materials, however, only emissions related to the ground 1$s$ state and the first excited 2$s$ state of the A exciton is observed at whole range of applied magnetic fields. It is in contrast to the WSe2 monolayer studied in the main text, in which the emissions of higher excitonic states is better visible for bigger values of magnetic fields.

S6 Dependence of excitonic diamagnetic coefficients in WSe2 monolayer

To test our assumption that the excitonic spectrum in the WSe2 ML encapsulated in hBN resembles a 3D hydrogen atom ($\sim-1/n^{2}$), we investigate dependence of the obtained diamagnetic coefficients $\sigma$ of excitonic states, $ns$, in this system. We found theoretically (see Eq. 13 in the SM) that the mean value of $r^{2}$ calculated with the aid of our Kratzer potential approach approximate the one for 3D hydrogen atom, in which $r^{2}$ parameters of excitonic states scales with $n^{4}$. With the aid of Eq. 13 and $\sigma=(er)^{2}/8\mu$, we calculate theoretical $\sigma$ values of excitons for WSe2 ML encapsulated in hBN with parameters: $\mu$=0.2 $m_{0}$stierSM and $\varepsilon=4.5$ geickSM . The theoretical values are compared with the experimental diamagnetic coefficients in Fig. S8. The theoretical dependence fits very well the experimental data up to the 4$s$ state.This additionally confirms that the excitonic spectrum of the WSe2 monolayer encapsulated in hBN corresponds to that of a 3D hydrogen atom. The apparent discrepancy between the theory and the experiment for the 5$s$ state results in our opinion from the small range of the low-field limit, which affects the determined $\sigma$ value for this state.

S7 Low temperature reflectance contrast spectra of S-TMD monolayers

The low temperature RC spectra of WSe2, MoS2, WS2, and MoSe2 encapsulated in hBN are presented in Fig. S9. We define the RC spectrum as $RC(E)=\frac{R(E)-R_{0}(E)}{R(E)+R_{0}(E)}\times 100\%$, with $R(E)$ and $R_{0}(E)$, respectively, the reflectance of the dielectric stack composed of a monolayer encapsulated in hBN supported by a bare Si substrate and of the two alone layers of hBN on top of Si substrate. Note that the presented spectra correspond to the PL ones shown in Fig. 4 in the main text. First, the spectra display two pronounced resonances labelled 1$s_{\text{A,B}}$ which arise from the ground state of the so-called A and B excitons liSM ; aroraWSe2SM ; arorMoSe2SM; molasSM . In addition to them, less pronounced features, labelled 2$s_{\text{A,B}}$ and 3$s_{\text{A}}$, appear at about 200 meV higher in energy as compared to the 1$s_{\text{A,B}}$ ones. The assignment of the 2$s_{\text{A}}$ and 3$s_{\text{A}}$ features to the first and the second excited state of the A exciton is straightforward and corresponds to many other investigations on S-TMD MLs encapsulated in hBN stierSM ; robertSM ; hanSM ; slobodeniukSM ; gerberSM. The origin of the 2$s_{\text{B}}$ is less clear, as it has not been reported so far. Due to the similar energy separation between 1$s_{\text{B}}$ and 2$s_{\text{B}}$ as compared with the 1$s_{\text{A}}$ and 2$s_{\text{A}}$, we ascribed tentatively the 2$s_{\text{B}}$ features to the first excited states of the B exciton, which, however, requires further investigations.

S8 Estimation of the band-gap energy in MoSe2 monolayer

As it has been discussed in the main text, the estimation of band-gap energy is essential for our analysis of excitonic ladder in S-TMD MLs. The estimation of band-gap energy can be easily carried out for WSe2, MoS2, and WS2 MLs, however, this issue is more complex for the MoSe2 one. This is, because PL related to the ground state of B exciton, $1s_{\text{B}}$, appears in the spectral range of the emissions related to the $2s$ and higher $ns$ states of the A exciton, labelled as $1s_{\text{A}}$ and $ns_{\text{A}}$ in Fig. S10. To determine the band-gap energy of MoSe2 ML, we use the procedure described in the main text. The PL intensity at the band-gap energy of WSe2 ML equals 5$\%$ of the intensity of the 2$s$ exciton PL peak. In order to apply the same approach for MoSe2 ML, we deconvolute the spectrum shown in Fig. S10 with three Gaussian profiles (the PL due to the $ns_{\text{A}}$ and $1s_{\text{B}}$ are resolved spectrally). We set the linewidth of the $1s_{\text{B}}$ emission equals to 28 meV, as obtained from the upconversion PL spectrum of MoSe2 ML reported in Ref. 17. The result of the procedure is presented in Fig. S10. The band-gap energy in MoSe2 ML determined using the procedure equals 1.861 eV and is marked in Fig. S10 with pink arrows.

S9 Application of the proposed model to the data available in the literature

In the main text of this report, we have clearly demonstrated that the simple $E_{ns}=E_{g}-Ry^{*}/(n+\delta)^{2}$ ansatz can be well applied to reproduce the energy ladder of excitonic $s$-resonances in WSe2 monolayer encapsulated in hBN. Recently, such ladders has been also extrapolated from observation of a series of excitonic resonances in very high magnetic fields, for other, MoS2, WS2 and MoSe2, monolayers encapsulated in hBN gorycaSM . Moreover, the observation of spectral series of $s$-type excitons has been also inferred from the refined analysis of the reflectance spectra in early reports on monolayer WS2 deposited on Si/SiO2 substrate chernikovSM .

It is interesting to test the applicability of our formula against the data reported in the above references. Thus, the energies of $s$-excitonic resonances, extracted from Refs. 20 and 21, have been fitted to the $E_{ns}=E_{g}-Ry^{*}/(n+\delta)^{2}$ formula, considering $E_{g}$, $Ry^{*}$, and $\delta$ as independent, adjustable parameters. As shown in Figs. S11 and S12, our formula fits perfectly the data for monolayers encapsulated in hBN and it is also well applied to the data for WS2 on Si/SiO2 system. The extracted $E_{g}^{\text{model}}$, $\delta^{\text{model}}$ and $E_{b}^{\text{model}}$ parameters ($E_{b}=Ry^{*}/(1+\delta)^{2}$) are listed in the Table 2. The corresponding energy values $E_{g}^{\text{data}}$ and $E_{b}^{\text{data}}$, originally reported in Refs. 21 and 20, are in perfect agreement with our findings.

According to our preceding discussion (see section S3), the applicability of our simple approach to the case of monolayers in hBN is well understood. The fact that our formula is operational for the data for WS2 on Si/SiO2 structure is surprising. As a matter of fact, these data follow our formula but only in sort of the “effective” way. Description of these data, within our approach and notably within the Rytova-Keldysh formalism in Ref. 21 as well, implies the use of rather unrealistic structure parameters, e.g., the underestimated dielectric constant ($\varepsilon=1$) and overestimated screening radius ($r_{0}=75\,\mbox{{\AA}}$).