Excitons in Atomically Thin TMD in Electric and Magnetic Fields

Monolayer (1L) Transition Metal Dichalcogenides (TMDs) are a class of two dimensional (2D) semiconducting materials that feature a direct band gap at the K points of the Brillouin zone, quadratic band structure and large band gaps Wang et al. (2018); Scharf et al. (2019). The strong Coulomb interaction in 1L TMDs results in huge exciton binding energies which causes excitons to dominate the absorption spectrum even at room temperature Wang et al. (2018); Scharf et al. (2019); Chen et al. (2019); Meckbach et al. (2018); Qiu et al. (2016, 2019). 1L TMD semiconductors are 2D systems with hexagonal Brillouin zones, however they lack inversion symmetry. Both the conduction and valence bands are spin split by the spin-orbit interaction and the spin and valley degrees of freedom are coupled Kormányos et al. (2015); Stier et al. (2016); Rostami et al. (2013); Wang et al. (2018); Scharf et al. (2019). The unique properties of 1L TMDs render them an ideal platform for the study of exciton physics  Chen et al. (2019); Wu et al. (2015); Liu et al. (2021); Goldstein et al. (2020); Liu et al. (2020); Bange et al. (2023); Zhu et al. (2023) and dynamical screening effects, particularly in regards to quasiparticle band gap renormalization Gao et al. (2016); Liang and Yang (2015); Qiu et al. (2019, 2016); Zibouche et al. (2021); Liu et al. (2019). Further, 1L TMDs are ideal for the development of optoelectronic devices Wang et al. (2012); Mueller and Malic (2018); Taffelli et al. (2021), including photosensor devices  Yin et al. (2012); Velusamy et al. (2015); Chang et al. (2014); Zhang et al. (2013); Lopez-Sanchez et al. (2013); Gonzalez Marin et al. (2019) and the fabrication of vertical and lateral heterostructures with advanced optical performance, with such devices providing an excellent environment for light-matter coupling in the form of exciton-polaritons.  Dufferwiel et al. (2018); Zhao et al. (2023). As such, accurate determination of the effective Hamiltonian parameters of 1L TMDs would be beneficial for modelling these technologies, and several others, including single photon emitters Sortino et al. (2021); Schuler et al. (2020); Srivastava et al. (2015); Tonndorf et al. (2015).

Exciton physics in 2D TMD systems is highly sensitive to reduced (band) mass, characteristic screening length and the external dielectric constant and hence the appropriate analysis of exciton resonances can be used to determine these quantities. Further to this, we will show how the electron and hole band mass asymmetry can also be extracted from excitonic spectroscopy, providing a near complete picture of the underlying band structure (effective Hamiltonian).

In the present work we consider excitons in magnetic and electric fields. While our conclusions are generic, to be specific we concentrate on WSe₂. There are two major messages of our work. (i) At realistic out-of plane magnetic fields a simple perturbation theory approach in magnetic field is not sufficient and one needs an exact solution to analyse existing data. Developing the solution and performing analysis of the data we determine parameters of the system which differ significantly from those obtained with simple perturbation theory. (ii) While s-wave TMD excitons have been observed in a number of photoexcitation experiments, p-wave excitons are invisible because of the extremely small spectral weight. Application of a moderate in-plane electric field makes p-wave exciton states visible in photoexcitation. We predict properties of the p-wave states: energies, photoexcitation probabilities, valley-orbital level splitting and orbital Zeeman effect due to the electron-hole mass-asymmetry. Experimental studies of p-wave states are absolutely feasible and they will allow the study of these effects and to establish rather precisely the mass-asymmetry.

(i) Measurements of diamagnetic energy shifts of exciton s-states provides a route to determine the effective parameters of 1L TMDs. The perturbation theory formula for the diamagnetic shift reads Landau and Lifshitz (2007),

with $e=|e|$ the elementary charge, $\mu$ the reduced mass, $\langle r^{2}\rangle$ the square of the radius of the exciton in a given quantum state, and $B$ the out-of-plane magnetic field strength. This formula is used to fit experimental data on TMD excitons Walck and Reinecke (1998); Chen et al. (2019); Stier et al. (2016). We will show that this formula, obtained from a perturbation theory in $B$, is insufficient at strong $B$-fields. And instead, we will derive the appropriate strong field behaviour of the diamagnetic shift. This was previously shown for a 2D Coulomb problem in magnetic field in the strong field and weak field limits MacDonald and Ritchie (1986); Laird et al. (2022) but in this work we consider a Keldysh potential and present a solution for arbitrary magnetic field, similar to Ref. Kezerashvili and Spiridonova (2021) however we extend the method to account for arbitrary orbital angular momentum and particle-hole mass asymmetry which is essential for real systems. Using our non-perturbative approach, we can provide accurate fitting to available experimental data, and in doing so, establish the key parameters for the low-energy TMD Hamiltonian.

(ii) Physics of p-wave exciton states is more rich than that of s-wave states. It includes the valley-orbital effect (the splitting of $l=\pm 1$ states without magnetic field) and the orbital Zeeman effect (the splitting of $l=\pm 1$ states in magnetic field). These phenomena are intimately related to the valley anomalous Hall effect in electron transport. Unfortunately the p-wave states are invisible in photoexcitation; these are “dark” states. We propose to make the states visible by applying an in-plane electric field. The electric field should be sufficuently strong to brighten the p-wave states relative to the s-wave states, but still remain weak enough such that bound states are not destroyed. We find the appropriate electric field strength is $E\sim 1-3$V/$\mu$m. We calculate values of the aforementioned effects in combined electric and magnetic fields.

The rest of the paper is organised as follows: Section II discusses the diamagnetic shift and existing photoluminescence experimental data. Section III walks through the derivation of the effective Hamiltonian for exciton in 1L TMD. In Section IV we present our results for s-wave excitons. We compare our theory to experimental results and extract reduced mass, dielectric constant and characteristic screening length. In Section V we consider excitons in combined in-plane electric and out-of-plane magnetic fields and show how such measurements can allow the study of anomalous effects in 1L TMDs. We summarise our conclusions in Section VI.

To stress the shortcomings of Eq.(1), we present Fig. 1, taken from experimental paper Ref. Chen et al. (2019), which plots the excitonic energies against magnetic field in 1L WSe₂ encapsulated in hBN. The data shows a combination of the diamagnetic shift, which is quadratic in B, and a linear in B valley Zeeman effect that we discuss in Appendix A. The valley Zeeman effect is not exciton state specific, it is common for all exciton states. To exclude the valley Zeeman effect the data should be symmetrised in B. It is clear by inspection of Fig. 1 that the 3s and 4s diamagnetic shifts are not quadratic at high $B$ and hence cannot be described by Eq.(1); the perturbation theory is insufficient. A strong field approach is necessary.

There is an interesting theoretical aspect of the problem that we would like to stress. In the weak B-field limit the expression for the diamagnetic energy shift, Eq.(1), is independent of the dimensionality of the exciton. However, in strong $B$-field the exciton dynamics strongly depends on dimensionality. In the three-dimensional (3D) case the exciton wave function is like a needle aligned with the magnetic field, and the Larmor circle oscillates along the needle resulting in a series of 1D Coulomb levels built on each Landau level  Landau and Lifshitz (2007). In the case of a 2D exciton the $z$-confinement does not allow this oscillation. The exciton remains confined to the $x$-$y$ plane at any field MacDonald and Ritchie (1986). Hence, diamagnetic shifts in strong B are very different in 3D and 2D cases.

In the presence of a magnetic field, the formation of an exciton can be described by the following Hamiltonian

Here $\hat{\mathbf{p}}_{\mathrm{e}(\mathrm{h})}$ is the momentum operator for the electron (hole), and $m_{\mathrm{e}(\mathrm{h})}$ is its effective mass. In the cylindrical gauge, the magnetic field B can be described by the vector potential $\mathbf{A}_{\mathrm{e}(\mathrm{h})}=[\mathbf{B}\times\mathbf{r}_{\mathrm{e}(\mathrm{h})}]/2$. If TMD monolayer is deposited at the insulating substrate, the Coulomb attraction is accurately described by the Rytova-Keldysh potential given by

Here $H_{0}$ and $Y_{0}$ are the Struve and Bessel functions of the second kind respectively and $r_{0}$ is the characteristic screening length. The length $r_{0}=2\pi\alpha/\epsilon$ is determined by the polarizability of the TMD monolayer $\alpha$ and effective dielectric constant $\epsilon$ of the surrounding media.

It is instructive to introduce the center of mass $\mathbf{R}=(m_{\mathrm{e}}\mathbf{r}_{\mathrm{e}}+m_{\mathrm{h}}\mathbf{r}_{\mathrm{h}})/(m_{\mathrm{e}}+m_{\mathrm{h}})$ and relative $\mathbf{r}=\mathbf{r}_{\mathrm{e}}-\mathbf{r}_{\mathrm{h}}$ coordinates. The corresponding momenta are ${\mathbf{P}}=-i\hbar\frac{\partial}{\partial{\mathbf{R}}}$ and ${\mathbf{p}}=-i\hbar\frac{\partial}{\partial{\mathbf{r}}}$. For the remainder of this work we shall use units $\hbar=c=1$. The Lamb transformation Lamb (1952) ${\hat{H}}\to U^{{\dagger}}{\hat{H}}U$ governed by the unitary operator $U=e^{i\frac{e}{2}\mathbf{B}\cdot[\mathbf{R}\times\mathbf{r}]}$ further simplifies Hamiltonian Eq.(2) as

with

We have also introduced the cyclotron frequency $\omega_{c}=eB/\mu$. We restrict ourselves only to the optically active excitons, i.e. excitons with zero center of mass momentum $\mathbf{P}=0$. Besides, the fourth term is determined by the orbital momentum for the relative motion $\hat{L}_{z}=[\mathbf{r}\times\mathbf{p}]_{z}$, which is a good quantum number $\hat{L}\to l$. It describes the orbital Zeeman shift of excitonic levels and is nonzero only for non-s-wave states.

The excitonic states are shaped by the interplay between the attractive Rytova-Keldysh potential and the harmonic confinement due to the magnetic filed. They smoothly evolve from the 2D hydrogenic-like series to the equidistant set of Landau levels. It is convenient to calculate the magnetic field dependence of excitonic states using the momentum space representation, i.e., $\psi(\mathbf{r})\to\psi_{\bf p}$. The eigenvalue problem with Hamiltonian (4) transforms to the integro-differential equation given by

where $V_{\bf p}=-2\pi e^{2}/\varepsilon p(1+r_{0}p)$ is the Fourier transform of the Rytova-Keldysh potential. It can be further simplified if we split the orbital part of the wave function $\psi_{\bf p}=\psi_{p}^{l}e^{il\theta_{p}}$ and rescale its radial part as $\psi_{p}^{l}=\chi_{p}^{l}/\sqrt{p}$, allowing the radial Hamiltonian to be written in a form that is strictly Hermitian. The resulting eigenvalue problem is given by

Here $V_{l}(p,p^{\prime})=\langle V_{\mathbf{p}-\mathbf{p}^{\prime}}e^{il\varphi}\rangle_{\varphi}$ is the multipole moment of interactions and involves averaging over the relative polar angle $\varphi$ between momenta $\mathbf{p}$ and $\mathbf{p}^{\prime}$. We solve Eq.(III) numerically using the linear algebra LAPACK package. Further details of these calculations are presented in Appendix B.

We wish to compare the measured magnetic field dependence of energies $E_{\mathrm{X}}$ for the s-wave excitonic resonances Chen et al. (2019) with the computed excitonic energies $\epsilon$. They are connected as $E_{\mathrm{X}}=\Delta+\epsilon$ where $\Delta\approx 1.9$ eV is the band gap in the TMD monolayer. Hence we have 3 fitting parameters: the reduced electron-hole mass $\mu$, the dielectric constant $\varepsilon$, and the screening length $r_{0}$. As presented in panel a of Fig. 2, our theoretical curves (solid lines) provide an excellent fit of the experimental data (asterisks) for all resolved excitonic states (1s - 4s)¹¹1We reiterate that the experimental points in Fig. 2 correspond to that of Fig. 1, yet averaged over $+B$ and $-B$. This averaging removes the linear valley Zeeman effect. The fitting parameters are

We make special mention of our fitted value $\varepsilon=3.9$, which is within the expected range for hBN Dean et al. (2010). With these parameters, the ground state binding energy at zero magnetic field is found to $|\epsilon_{1s}|=173$meV.

The original experimental work Chen et al. (2019) has used the simple parabolic approximation (1) for the diamagnetic shift. This has led to a different set of parameters: $\mu=0.22$ m₀, $\varepsilon=4.5$, and $r_{0}=4.5$ nm/$\varepsilon$. The theoretical curves (dashed lines in Fig. 2a) calculated with this set of parameters fit the data well only for the ground excitonic state, but the discrepancy for the excited excitonic states is evident. Our conclusion aligns well with a recent unpublished work ²²2David de la Fuente Pico, Jesper Levinson, Meera M. Parish and Francesca Maria Marchetti, private communication..

The behaviour of the eigenenergies of the Rydberg states strays from the parabolic diamagnetic approximation and it can be deduced from Eq.(III) that in the large $B$ limit, $B\to\infty$, the eigenergies approach those of 2D Landau levels. In this case the eigenenergy is quantized in radial quantum number $n_{r}$ and angular momentum $l$ as $\epsilon_{n_{r},l}=\omega_{c}\left(2n_{r}+|l|+\nu l+1\right)/2$ and is therefore linear in magnetic field $B$. The radial quantum number $n_{r}$ is related to the principle quantum number $n$ as $n_{r}=n-|l|-1$. In panel b of Fig. 2 we present the calculated first derivative of the excitonic energy calculated from Eq.(III) with respect to magnetic field. Solid lines in Fig. 2b correspond to solid lines in Fig. 2a. Dashed lines in Fig. 2b correspond to asymptotic values of the derivatives ($B\to\infty$) determined by Landau levels in 2D. While in the experimental region $0<B<30$T the derivative of the 1s state is practically linear in $B$ in accordance with perturbation theory, for higher states deviations from linearity are evident. Thus it is clear that the Rydberg states transition from the diamagnetic regime towards the Landau level regime at experimentally realistic $B$, rendering the quadratic diamagnetic approximation Eq.(1) invalid for the Rydberg excitons.

The region of validity of the diamagnetic approximation is also illustrated by considering the size of the exciton, or the root mean square (rms) radius $r_{rms}$. The size of the exciton may be calculated directly from the wavefunction.

The exciton size $r_{rms}=\sqrt{\langle\psi|r^{2}|\psi\rangle}$ is plotted against magnetic field in panel a of Fig. 3. The 1s exciton radius is $\sim 2$ nm and remains smaller than the magnetic length $l_{B}=\sqrt{e/B}$ over the entire range of magnetic field, however this is not true for the Rydberg excitons. Thus, to accurately model the Rydberg states over this range of $B$ the exact Hamiltonian (III) is required.

The physics of p-wave exciton states is richer than that of s-wave states. It includes the valley-orbital effect (the splitting of $l=\pm 1$ states without magnetic field) and the orbital Zeeman effect (the splitting of $l=\pm 1$ states in magnetic field). Unfortunately the spectral weight of direct photoexcitation of a p-wave exciton calculated in Appendix C is too small for direct detection, $\sim 10^{-3}$ in relative units compared with the s-wave state. However, one can mix s and p states by applying in-plane electric field and hence excite the p state due to admixture of the s state.

There is a recent experiment where the p-wave states have been observed using this method Zhu et al. (2023). However, in experiment Zhu et al. (2023) the electric field was relatively strong and an interpretation in terms of simple electron-hole bound states (excitons) is questionable. The major part of Ref. Zhu et al. (2023) data corresponds to events above the exciton ionisation limit where the observed resonances correspond to doorway scattering states ³³3O. P. Sushkov, J. N. Engdahl and D. K. Efimkin, to be published. Contrary to this in the present work we consider a relatively weak electric field when the description in terms of bound excition states is correct.

Energies of np-states with or without magnetic field can be calculated using Eq.(III). The np level is always slightly lower than the ns-level. However, Eq.(III) is missing an important effect: the valley-orbital level splitting between the states $p_{+}$ and $p_{-}$ that correspond to angular momentum $l=+1$ and $l=-1$, where this splitting exists even at B=0. Using the language of Dirac equation one can say that Eq.(III) corresponds to the “non-relativistic approximation” to the Dirac equation that at $B=0$ gives degenerate $p_{\pm}$ states, and the $p_{+}$ - $p_{-}$ splitting is the 1st “relativistic” correction to this equation, The splitting is discussed in Appendix D, it is small, from only a few to several meV, and if it is positive in one valley, $\Delta\epsilon_{\tau}=\epsilon_{-}-\epsilon_{+}>0$, it is negative in the other valley.

Having in mind that energies of $ns$, $np_{+}$ and $np_{-}$ states with the same n are close and assuming that the electric field $E=E_{x}$ is not too strong we can restrict our analysis to a 3-level approximation with the effective Hamiltonian

Here $\epsilon_{a,0}$ is the eigenenergy of the state with angular momentum $l_{a}$ calculated with Eq.(III). $\Delta\epsilon_{\tau}$ is the valley-orbital splitting, and $v_{a,b}=eE_{x}x_{a,b}$, where

In panel b of Fig. 3 we present plots of the electric dipole radial matrix element $r_{np,ns}$ for $n=2,3,4$ as a function of magnetic field.

Before diagonalization of (10) the energies $\epsilon_{s,0}$, $\epsilon_{p_{+},0}$ and $\epsilon_{p_{-},0}$ must be calculated from Eq.(III). From the fit of the s-wave states we know all the parameters of this equation except for $\nu$, the particle-hole asymmetry term defined in (5). Here we take $\nu=-0.1$, which is obtained from the effective electron and hole masses determined from DFT calculations in Ref. Kormányos et al. (2015). However, armed with our formalism, we suggest that analysis of future experiments can establish $\nu$, providing vital information on the underlying band structure. The valley-orbital splitting $\Delta\epsilon_{\tau}$ is discussed and calculated in Appendix D. In principle it depends on magnetic field, but for n=2 this dependence is negligible over the range of magnetic field that we consider. The splitting depends on the velocity in the effective Dirac equation, see Appendix D. We approximate the velocity from the band gap and the effective mass as $v^{2}\approx\Delta/4\mu$, which results in a valley-orbital splitting $\Delta\epsilon_{\tau}\approx\pm 3.8$meV. Finally we reiterate the point that we already have made in Section II: the linear in B valley Zeeman effect that is not exciton specific and can be determined from the 1s state data is subtracted from all our results. After this preparatory procedure, diagonalization of Eq.(10) is straightforward. We use nomenclature $s^{\prime}$, $p_{+}^{\prime}$ and $p_{-}^{\prime}$ to denote the states that originate from $s$, $p_{+}$ and $p_{-}$ at $E=B=0$. The photoexcitation probability is proportional to the weight of the bare s-wave state in the corresponding wave function. In Fig.4 we plot the probabilities versus electric field for n=2 states. The probabilities are normalized such that $P_{2s}=1$ at $E=0$. As previosuly mentioned, the approach used in the present work fails at strong electric fields. This is due to the effective potential of the Keldysh potential combined with the strong electric field failing to permit exciton bound states. A separate analysis of strong electric fields ⁴⁴4O. P. Sushkov, J. N. Engdahl and D. K. Efimkin, to be published shows that for n=2 states in WSe₂ the approach employed in the present work is valid up to $E=3-4$V/$\mu$m, above this field bound n=2 exciton states do not exist. This is where we terminate our plots in Fig.4. Here the probability of the $p^{\prime}$-state excitation can be up to 20% of that of the s-wave. The relative probabilities in Fig.4 assume linear polarization of the excitation laser that does not separately select single valleys. Panel a in Fig.4 displayes probabilities at B=0. Plots for different valleys are identical, but what is the $p_{+}^{\prime}$ curve for one valley is the $p_{-}^{\prime}$ curve for the other valley. Panel b in Fig.4 displayes probabilities at B=30T.

In Fig.5 we plot energies of $2s^{\prime}$, $2p_{+}^{\prime}$, and $2p_{-}^{\prime}$ states versus external in-plane electric field E for two values of magnetic field, $B=0$ and $B=30$T. These plots depend on the valley-orbital energy splitting between the $p_{\pm}$ levels, $\Delta\epsilon_{\tau}$, and also on the mass asymmetry parameter $\nu$-term in Eq.(4). In Fig.5 we assume values that follow from DFT calculationsKormányos et al. (2015) and indirect analysis of ARPES dataNguyen et al. (2019). Measurements of these energy levels would shed light on true values of these parameters.

We have shown that the perturbation theory approximation for diamagnetic shift of exciton energy levels is not sufficient to describe existing data. We develop the theory for strong magnetic field, reanalyse the data for monolayer WSe₂, and arrive to the set of exciton parameters presented in Eq.(8). This set is different from what was known in literature, although the reduced mass is similar to that obtained from DFT calculations.

Only s-wave excitons are visible in photoluminescence, the spectral weight of p-wave states is too small at just $\sim 10^{-3}$ relative to the s-wave. We propose to study p-wave states by mixing them with s-wave states by application of external in-plane electric field and show that a moderate electric field $E\approx 3V/\mu m$ results in the intensity about 20% relative to s-wave. We calculate energy levels in combined electric and magnetic fields and demonstrate that our proposal opens a path to experimentally study the valley-orbital level splitting and the electron-hole mass asymmetry.

We acknowledge discussions with Zeb Krix, Francesca Marchetti and Jesper Levinson. This work was supported by the Australian Research Council Centre of Excellence in Future Low- Energy Electronics Technologies (CE170100039).