Theory of magnetotrion-polaritons in transition metal dichalcogenide monolayers

Introduction
Monolayers of transition metal dichalcogenides (TMD) are representatives of the family of two-dimensional (2D) materials which reveal robust excitonic response. The corresponding bright exciton binding energies reach hundreds of meV1, and the presence of tightly bound excitons with large optical oscillator strength determine the optical spectra of TMD monolayers up to room temperatures2. Moreover, in contrast to the case of conventional semiconductor quantum wells, excited excitonic states are clearly seen in the reflection spectrum3. They are characterized by non-Rydberg energy scaling resulting from the reduced screening of Coulomb interaction peculiar for a truly 2D system4,5,6,7.
TMD monolayers have a hexagonal lattice structure with band edge minima at ± K points of the Brillouin zone forming two non-equivalent valleys characterized by opposite spin polarizations8. This gives rise to emergent spin-valley physics9, which in the domain of excitonics is characterized by a large variety of exciton species10, and determines the great potential of the considered materials for nanoscale device applications11,12,13.
A direct consequence of tightly bound character of excitons in TMD monolayers is an ability to reach the regime of high density of excitons with n ~ 1012–1013 cm−2 14,15, which is orders of magnitude larger than in conventional 2D systems based on GaAs quantum wells16. In this regime the nonlinear optical effects associated with exciton-exciton Coulomb interactions become prominent17,18,19,20,21,22. Such nonlinear behavior is further enhanced in the presence of an optical microcavity, where the regime of strong light-matter coupling occurs and exciton-polaritons are formed23.
Excitonic and polaritonic properties of TMD monolayers are strongly affected by the presence of free electron gas. The unintentional doping routinely appears during the fabrication process, and the density of excess electrons can be controlled via external gating field24, which thus becomes a control parameter for tuning of the exciton resonance. Moreover, the presence of free charge carriers results in a formation of three-body charged excitonic complexes—trions, characterized by a new emergent peak below the excitonic line in the optical spectrum25,26,27,28. It should be noted that at high densities of free carriers the underlying physics of exciton-electron interaction becomes more complex, and a crossover from trions to repulsive (exciton) and attractive (trion) polarons occurs29,30,31,32. Yet, at low densities of itinerant electrons, the exciton-trion picture is sufficient for description of the optical response33.
Similar to excitons, trions are efficiently coupled to light. The regime of strong light-matter coupling in microcavities with doped TMD monolayers is well elaborated34,35. Interestingly, the nonlinear response associated with trion polaritons is larger as compared to excitons36,37, and can be further enhanced by means of an external magnetic field38, which was also used to reveal the Rydberg series, and determine effective masses and g-factors of excitons39,40,41,42,43,44. However, the existing studies of the magnetic field impact on trions in TMD monolayers focus on a valley Zeeman splitting only45,46 with orbital effects being typically neglected33. In this paper we thus develop a microscopic formalism where orbital and spin effects of a magnetic field are treated on equal footing, demonstrating that in the case of trions they give comparable contributions. We also demonstrate the vital role played by the mixing of different orbital channels.
Results and discussion
Trion states in the presence of a magnetic field
We consider the magnetic field impact on trion-polariton resonances in monolayer WSe2, as schematically shown in Fig. 1a. The detailed description of the microscopic model is presented in the Methods section. The model takes into account band structure of the valleys near the ± K points of the Brillouin zone, as shown in Fig. 1b.

a Sketch of the structure. A hBN-encapsulated TMD monolayer hosting excess electrons is embedded in an optical microcavity in the presence of a normal-to-plane uniform magnetic field. The microcavity is filled with dielectric media represented by PMMA. A pair of electrodes consisting of a metallic gate and graphene layer is used to control the density of excess electrons. A cavity photon creates an electron-hole pair, which is associated with a free electron forming a three-particle bound state (a trion). The radiative recombination of a trion generates a cavity photon and a free electron. In high-Q optical cavity, a successive repetition of these processes leads to the onset of strong light-matter coupling regime and formation of hybrid states, trion-polaritons. b Band structure in the vicinity of ±K points of Brillouin zone and optical excitonic transition (thick arrow) in WSe2 monolayer in the presence of free electron gas. The arrows indicate the subband’s spin projection. The lowest-energy trion state in the -K valley is formed by an electron from the upper conduction subband and a hole from the upper valence subband, coupled with an excess electron, located in the +K valley of the lower conduction subband.
We start with addressing the exciton energy EX defined as the ground state eigenvalue of electron-hole Hamiltonian in the absence of excess charge carriers. The exciton problem in a weak magnetic field is well studied47, including the case of TMD monolayers40,41,42,43,44. The magnetic field results in a diamagnetic energy shift, which in the lowest order of perturbation theory amounts to ΔEdia ≈ e2〈r2〉B2/8μ, where μ = μe1μh/(μe1 + μh) is the exciton reduced mass, with μe1, μh denoting the effective masses of upper subband electron and a hole, respectively. The quantity (langle {r}^{2}rangle =leftlangle {psi }_{{rm{X}}}rightvert {r}^{2}leftvert {psi }_{{rm{X}}}rightrangle) is calculated in the absence of the magnetic field, where ψX is the wave function of exciton relative dynamics. To treat the problem more rigorously, we invoke the stochastic variational method (SVM). The numerical computation efficiency of the method greatly depends on the chosen trial function. Specifically, we use an ansatz of fully correlated two-dimensional Gaussians. The trial wave function is formed as a linear combination of these basis functions, optimized through random adjustments of parameters to minimize the ground state energy. Since the basis functions are not mutually orthogonal, we solve a generalized eigenvalue problem to find the energy levels and corresponding eigenvectors, providing variational upper bounds for the energies of the respective states. Further details of the model and computation scheme are presented in Methods section.
We employ two sets of basis functions. The first one includes a single orbital channel — (me, mh) = (0, 0), as in the case considered in ref. 48. The results of calculation are shown in Fig. 2a, b for monolayers WSe2, and MoSe2, respectively. The material parameters are given in Table 1. While in the absence of the magnetic field, the exciton binding energy is well reproduced, this approach demonstrates a linear (instead of quadratic) scaling with respect to the magnetic field, and severely overestimates the diamagnetic shift. On the contrary, by using the multi-orbital basis with the orbital channels (0, 0), (1, −1), (−1, 1), (−2, 2), and (2, −2), one obtains the results which agree well with the experimental data and with the simple perturbative estimate discussed above.

Exciton energy versus the magnetic field calculated for atomic monolayer a WSe2, and b MoSe2 by means of three different approaches. The dotted lines correspond to the single channel (SC) calculations within the SVM method, where one orbital channel (me, mh) = (0, 0) is taken into account as discussed in ref. 48. The resulting diamagnetic shift scales linearly with the magnetic field and is essentially larger compared to the experimental measurements. The dashed lines correspond to the calculation within perturbation theory (PT), where the diamagnetic shift is estimated via ΔEdia ≈ e2〈r2〉B2/8μ, in accordance with the experimental results40,41,42,43,44. The solid lines show the SVM calculations of the exciton energy with a multi-orbital basis incorporating the following orbital channels: (0, 0), (1, −1), (−1, 1), (−2, 2), and (2, −2).
We next numerically analyze the generalized eigenvalue problem defined by Eq. (11) of Methods for the trion states.
We adopt the multi-orbital basis with all possible combinations −2 ≤ mi ≤ 2 such that ∑imi = 0. Such truncation of basis channels is enough for obtaining adequate precision of the calculation. We check the convergence of results expanding the basis to the case −3 ≤ mi ≤ 3 [see Supplementary Note 1]. The results of the calculations are summarized in Fig. 3. In Fig. 3a we display the trion binding energy as a function of the magnetic field for various TMD monolayers. First of all, we note that the value of the trion binding energy in monolayer WSe2 is above 20 meV in the absence of the magnetic field. Together with the exciton binding energy of ∣EX(B = 0)∣ = 146.3 meV, this provides essentially better agreement with experimental data than the results reported in ref. 28. The larger value of the trion binding energy was obtained due to different effective masses of electrons in the spin subbands. The presence of a magnetic field leads to the increase of the trion binding energy, defined as
where ET is the ground-state eigenvalue of Eq. (11), ωe = e∣B∣/μe2 is the cyclotron frequency of free electrons. The sublinear increase of trion binding energy with increasing magnetic field is due to the competition of linear in magnetic field shift of free electron Landau level ℏωe, and the diamagnetic shift of trion eigenenergy ET. We highlight that the suggested model adequately describes trion states at moderate density of free electron gas ({E}_{{rm{F}}},ll, {E}_{{rm{T}}}^{{rm{b}}})33. Here EF is the Fermi energy of free electron gas, defined in Methods.

a Trion binding energy versus the magnetic field, calculated using the expression (1) for various TMD monolayers. b Trion resonance energy shift as a function of the magnetic field for WSe2 and MoSe2 atomic monolayers. The dotted lines correspond to the diamagnetic shift. The solid and dashed lines demonstrate the total shift accounting for the Zeeman splitting for the two spin orientations.
In Fig. 3b the trion resonance energy shift versus magnetic field for monolayers of WSe2 and MoSe2 is shown. The pure diamagnetic shift defined as (Delta {E}_{{rm{T}}}^{{rm{dia}}}({boldsymbol{B}})=| {E}_{{rm{T}}}({boldsymbol{B}})-{E}_{{rm{T}}}({boldsymbol{B}}=0)|) is depicted by the dotted lines. We account for the magnetic field-induced Zeeman splitting of trion energy ΔETZ = gXμBB33, where μB is the Bohr magneton, gX = gc − gv is the exciton Lande factor defined via the Lande factors of upper conduction subband gc, and valence band gv49. The overall shift of trion resonances for the two spin orientations
is displayed in Fig. 3b by solid and dashed lines. Our direct calculations indicate that at the scale of 10 T of the applied magnetic field, the diamagnetic shift of the trion energy is comparable with the energy of Zeeman splitting.
Trion-photon coupling
We now consider the regime of strong light-matter coupling of trions with a near-resonant cavity eigenmode characterized by Eq. (29) of Methods. The negative trion resonance in hBN-encapsulated WSe2 in the absence of the magnetic field and at low density of free carriers corresponds to ({E}_{{rm{T}}}^{0}=1685) meV28. We use this value to choose the resonance energy of a cavity mode according to (hslash {omega }_{{rm{C}}}={E}_{{rm{T}}}^{0}). The dielectric constant of the media is set ε = 3, typical for PMMA50, and the cavity length is found from the resonance condition in the form ({L}_{{rm{C}}}=2pi c/(sqrt{varepsilon }{omega }_{{rm{C}}})), where c is the speed of light. We set the density of excess electrons at ne = 2 × 1011 cm−2, and also set T = 4 K. The dipole matrix element of the interband transition is chosen as dcv = 18 Deb, resulting in ΩT ≈ 5 meV in the absence of the magnetic field, typical for the trion-photon coupling36,38,51.
The polarization-resolved magnetic field dependence of the polariton energies
for monolayer WSe2 is shown in Fig. 4a. The σ+ and σ− polarized branches demonstrate the opposite behavior with increasing magnetic field: the Rabi splitting of the former reduces, while for the latter it increases. The reason is as follows. The magnetic field results in an imbalance of excess electrons, where a major part of them is transfered into +K valley, see Eq. (18) of Methods. As the trion photon-coupling rate scales as ({Omega }_{{rm{T}}}propto sqrt{{n}_{e}}), this gives rise to the reduction of trion-photon coupling for σ+ polarization, and enhancement for σ− polarization, as shown in Fig. 4b. The electron gas becomes fully polarized around 20 T of applied magnetic field, resulting in complete suppression of the trion-photon coupling for σ+ polarization. We note here that the trion-photon coupling ΩT depends on the magnetic field via the factor IT as well, defined via the trion wave function, see Eq. (23). This particularly results in a minor reduction of trion-photon coupling at elevated values of magnetic field, as depicted in Fig. 4b.

a Energy of trion-polariton branches versus the magnetic field for atomic monolayer WSe2. Red [blue] curves correspond σ+[−] circular polarizations. The increase of the magnetic field leads to a collapse of Rabi splitting for σ+-polarized trions, and the respective increase of Rabi splitting for σ−-polarized trions. b Trion-photon coupling energy versus the magnetic field calculated via Eq. (26). c Zeeman splitting of bare trions (dotted line) and trion-polaritons (solid line). The reduction [enhancement] of Rabi splitting with increasing magnetic field for σ+[−]-polarized results in a large increase of Zeeman splitting of trion-polaritons.
The valley-dependent nature of the trion-photon coupling ΩT in the presence of a magnetic field results in drastic modulation of the trion Zeeman splitting. As the lower trion-polariton energies of the two polarizations shift in opposite directions with increasing magnetic field, in the strong-coupling regime, the effective Zeeman splitting (Delta {E}_{{rm{T}}}^{{rm{LP}}}={E}_{{rm{T}}}^{{rm{LP+}}}-{E}_{{rm{T}}}^{{rm{LP-}}}) essentially increases, as shown in Fig. 4c. This result agrees well with a recent experimental observation of giant Zeeman splitting of trion-polaritons in MoSe2-based structure38. The saturating character of effective Zeeman splitting at large magnetic field is due to complete polarizaion of the free electron gas. We note that the enhanced trion-polariton Zeeman splitting will result in respective circular polarization of photoluminescence spectra52.
The spectra of WSe2 atomic monolayer is characterized by two trion peaks, attributed to intravalley and intervalley trions28. While this additional trion peak has considerably smaller oscillator strength53, its presence in the strong light-matter coupling Hamiltonian can modify the effective Zeeman splitting of lower trion-polariton branch (see Supplementary Note 2 for the details).
We next consider the impact of free electron gas density on trion-polariton effective Zeeman splitting. Generally the variation of free electron gas density affects the trion states via the additional screening of Coulomb interaction, and the Pauli blocking effects, resulting in slight modulation of trion wave function and spectral position32,35,54. Here we neglect this effect, assuming that the free electron gas density primarily affects the trion-photon coupling as ({Omega }_{{rm{T}}}propto sqrt{{n}_{e}}), see Eq. (26) in Methods. The trion-polariton effective Zeeman splitting for varying free electron gas density is depicted in Fig. 5a. The magnetic field-induced free electron valley polarization rate is defined via the ratio (Delta {E}_{{rm{TZ}}}(B)/{E}_{{rm{F}}}^{0}), see Eqs. (18), (17) in Methods. Hence, for smaller density ne the complete valley polarization occurs at smaller magnetic field, as shown in Fig. 5a. One should note, however, that the further reduction of free electron gas density will result in the collapse of strong coupling regime between trion and photon states.

a Zeeman splitting of bare trions (black dotted line) and trion-polaritons for different values of free carrier density indicated on the color bar (solid orange curves). b Dependence of effective ({g}_{{rm{T}}}^{{rm{eff}}}=Delta {E}_{{rm{T}}}^{{rm{LP}}}(B)/({mu }_{{rm{B}}}B)) at B → 0 on the free electron density ne. The orange dots of varying color correspond to lines on the (a), and black dotted line indicates the bare g-factor gX = 4.
The renormalization of trion Zeeman splitting in strong light-matter coupling regime can be further quantified via an effective g-factor, defined as ({g}_{{rm{T}}}^{{rm{eff}}}=Delta {E}_{{rm{T}}}^{{rm{LP}}}(B)/({mu }_{{rm{B}}}B)) at B → 0. At moderate density of free electrons ne = 1011 cm−2 one has ({g}_{{rm{T}}}^{{rm{eff}}}=10.78), about three times larger as compared to gX = 449. The modulation of ({g}_{{rm{T}}}^{{rm{eff}}}) with varying free electron gas density is shown in Fig. 5b.
Discussion
We developed a microscopic formalism for the description of magnetotrions and magnetotrion-polaritons in TMD monolayers. The spin and orbital effects of the magnetic field were treated on the same footing, and it was demonstrated that they give comparable contributions to the trion energy. We also showed that if a trion resonance is coupled to an optical cavity mode, the Zeeman splitting in the conductance band leads to a strong valley and polarization dependence of the trion-polariton energies, ultimately resulting in an effective giant Zeeman splitting for trion-polaritons.
Methods
Trion states in the presence of a magnetic field
We consider a trion state in the presence of a normal-to-plane uniform magnetic field B in monolayer WSe2 placed inside a planar optical cavity represented by a pair of distributed Bragg reflectors (DBR), as shown in Fig. 1a. The excess electrons are located in the lowest spin split subband in the vicinity of K point of Brillouin zone, whereas an optical interband transition occurs to the upper subband (see Fig. 1b). This separation allows one to discard the phase space filling effects associated with Pauli blocking33. The Hamiltonian of a three-particle system consisting of two electrons with coordinates r1,2 and a hole with coordinate rh is given by
where
Here the momenta operators read ({hat{{boldsymbol{pi }}}}_{j}=-ihslash {{boldsymbol{nabla }}}_{j}+e{boldsymbol{A}}({{boldsymbol{r}}}_{j})), ({hat{{boldsymbol{pi }}}}_{h}=-ihslash {{boldsymbol{nabla }}}_{h}-e{boldsymbol{A}}({{boldsymbol{r}}}_{h})), e is the positive unit charge, and A is the vector gauge potential, ∇ × A = B. The parameters μe1, μe2, μh are the effective masses of upper and lower subband electron, and a hole, respectively. The Coulomb interaction between charge carriers can be described by means of the Rytova-Keldysh potential4,5:
where ε0 is the static dielectric constant, r0 is the effective screening length, κ is the effective dielectric permittivity of surrounding media represented by a hexagonal boron nitride (hBN). H0 and Y0 are the Struve function and Bessel function of the second kind, respectively.
We employ a variational approach for the calculation of the ground state energy and wave function of the bound state. The numerical efficiency of the method strongly depends on the type of a trial function. The possible choice is the use of the anzats of fully correlated two-dimensional Gaussians55,56:
where ({mathcal{A}}) is the antisymmetrization operator, r = (r1, r2, r3), the superindex M = (m1, m2, m3), ({chi }_{S{M}_{S}}) are linear combinations of single-particle spin functions corresponding to a given spin S and spin projection MS. Throughout the paper we assume r3 ≡ rh.
The basis functions (10) are used to construct a trial wave function as a linear combination Ψ = ∑iciΦi, where the superindex i carries information about the orbital and spin quantum numbers. The functions Φi should be optimized within the procedure of random parameter modification in order to minimize the ground state energy. Given that these functions are not orthogonal to each other, one has to solve a generalized eigenvalue problem:
where as the output we obtain E1, E2, E3, … as variational upper bounds of energies of respective states and corresponding eigenvectors. In what follows we denote the trion ground state by ΨT(r1, r2, r3). We perform a similar procedure for the exciton problem, where the ground-state wave function will be denoted by ΨX(r1, r2).
The parameters of matrix ({A}_{ij}^{n}) can be found by means of the SVM method, which has proven its high efficiency in typical eigenvalue problems for few-body systems57,58,59,60. Detailed description of the numerical procedure is presented in refs. 48,55,56,61. Here we use the svm-2d package55, in which we recalculate the matrix elements of the interparticle interaction with the potential (9) by Gauss-Legendre quadrature. We first make sure that in the absence of the magnetic field, the resulting trion binding energy matches with the previous reports61. As we only consider the ground state, we set for an exciton S = 0, MS = 0, and for a trion S = 1/2, MS = 1/2.
Valley distribution of free electron gas
No magnetic field. In the absence of magnetic field the density of states for each valley is D±(E) = μe2/(2πℏ2), and the electron density
Here ({E}_{{rm{F}}}^{0}) is the being the Fermi energy in the absence of the magnetic field, kB denotes the Boltzmann constant, and T stands for the temperature. At the limit ({E}_{{rm{F}}}^{0},gg, {k}_{{rm{B}}}T) one has
yielding in ({E}_{{rm{F}}}^{0}=pi {hslash }^{2}{n}_{e}/{mu }_{e2}).
Magnetic field impact. We consider the impact of magnetic field accounting for the valley Zeeman splitting of the lower conduction subband ΔEZ,c2 = gc2μBB, where gc2 is the Lande factor of the lower conduction subband. For the sake of simplicity we discard the formation of Landau levels and the respective modulation of 2D electron gas density of states, which can lead to some quantitative modification of the results. To obtain the Fermi energy dependence on magnetic field, the following conditions must be satisfied: (i) ({n}_{e,+{rm{K}}}+{n}_{e,-{rm{K}}}={n}_{e}={rm{const}}), (ii) ({E}_{{rm{F}}}^{-}={E}_{{rm{F}}}^{+}={E}_{{rm{F}}}), (iii) ne,±K ≥ 0. One has
yielding in
Solving with respect to Fermi energy EF, we get
The electron densities in two valleys then read
Trion-photon coupling
We consider a direct interband optical transition within a two-band approximation. We consider the case of the normal incidence, so that the photon wave vector qph = 0. The corresponding Hamiltonian can be represented as
with
where ωC is the cavity eigenfrequency, ε is the dielectric constant of a media typically represented by an organic polymer, such as polymethyl methacrylate (PMMA)62,63,64. Then LC and S are the cavity length and sample area, respectively, dcv is the dipole matrix element of the interband transition, ({hat{a}}_{{{boldsymbol{k}}}_{e}}^{dagger }), ({hat{b}}_{{{boldsymbol{k}}}_{h}}^{dagger }) are the creation operators for electrons and holes, respectively. The trion state can be expressed in the form
where ({hat{tilde{a}}}_{{{boldsymbol{k}}}_{2}}^{dagger }) corresponds to excess electrons, and FT(k1, k2, kh) is the trion wave function in the momentum space. Accordingly, the excess electron state is (leftvert erightrangle ={hat{tilde{a}}}_{{boldsymbol{k}}}^{dagger }vert phi rangle).
The coupling rate is then evaluated via
In the presence of a magnetic field, the center-of-mass dynamics of trions cannot be straightforwardly decoupled from its internal dynamics, and one obtains
For the case of moderate dopings (| {boldsymbol{k}}| sim | {{boldsymbol{k}}}_{F}|, ll, {a}_{{rm{tr}}}^{-1}), where ({a}_{{rm{tr}}}=leftlangle {Psi }_{{rm{T}}}rightvert rleftvert {Psi }_{{rm{T}}}rightrangle) is the effective Bohr radius of trion, and ({{boldsymbol{k}}}_{F}propto sqrt{{n}_{e}}) is the Fermi wave vector defined by excess electron density ne. Hence, one can assume eikr ≈ 1, so that ({I}_{{rm{T}}}({boldsymbol{k}})approx {I}_{{rm{T}}}={rm{const}})33.
We further introduce quasi-bosonic operators36
where Ne,±K = ne,±KS is the total number of excess electrons in each valley. Then we find
where
In the absence of the magnetic field, the trion wave function can be factorized as ({Psi }_{{rm{T}}}({{boldsymbol{r}}}_{1},{{boldsymbol{r}}}_{2},{{boldsymbol{r}}}_{h}){| }_{{boldsymbol{B}} = 0}={psi }_{{rm{T}}}({{boldsymbol{r}}}_{1}-{{boldsymbol{r}}}_{h},{{boldsymbol{r}}}_{2}-{{boldsymbol{r}}}_{h})/sqrt{S}), where ψT is the wave function of relative dynamics, and the center-of-mass dynamics is represented by a plane wave with a zero wave-vector. Then the trion-photon coupling is simplified:
Analogously, the exciton-photon coupling reads33
Note that similarly to the case of trions, the dynamics of center of mass cannot be decoupled from the internal dynamics if the magnetic field is present. In the absence of the magnetic field, one has ({Psi }_{{rm{X}}}({{boldsymbol{r}}}_{e},{{boldsymbol{r}}}_{h})={psi }_{{rm{X}}}(| {{boldsymbol{r}}}_{e}-{{boldsymbol{r}}}_{h}| )/sqrt{S}), where ψX is the wave function of relative dynamics involved in the well-known expression ({Omega }_{{rm{X}}}=Omega sqrt{S}{psi }_{{rm{X}}}(0))65.
Given the large energy separation (above 30 meV) between the exciton and trion optical transitions28, we consider a strong coupling of trions and cavity modes only, i.e., we disregard the exciton resonance. The corresponding Hamiltonian reads
where the superscript ± denotes circular polarization of the cavity mode and the trion valley index. The eigenmodes of this Hamiltonian correspond to the upper and lower trion-polaritons.
Responses