{'fraction_non_alphanumeric': 0.035784563807302674, 'fraction_numerical': 0.014539756577508929, 'mean_word_length': 4.588067603339442, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microscopic theory of exciton and trion polaritons in doped monolayers of transition metal dichalcogenides

We study a doped transition metal dichalcogenide (TMDC) monolayer in an optical microcavity. Using the microscopic theory, we simulate spectra of quasiparticles emerging due to the interaction of material excitations and a high-finesse optical mode, providing a comprehensive analysis of optical spectra as a function of Fermi energy and predicting several modes in the strong light-matter coupling regime. In addition to exciton-polaritons and trion-polaritons, we report polaritonic modes that become bright due to the interaction of excitons with free carriers. At large doping, we reveal strongly coupled modes corresponding to excited trions that hybridize with a cavity mode. We also demonstrate that rising the carrier concentration enables to change the nature of the system's ground state from the dark to the bright one. Our results offer a unified description of polaritonic modes in a wide range of free electron densities.

Monolayers of transition metal dichalcogenides (TMDC) represent a class of two-dimensional (2D) materials with remarkable optical properties. They are direct bandgap semiconductors with hexagonal Brillouin zone, with two nonequivalent valleys at the K and K points 1-3 . Due to large spin-orbit interaction, they experience spin-valley locking, which opens the way for spin and valleytronic applications 4 . Moreover, the direct bandgap type and relatively large effective masses of electrons and holes lead to the formation of robust bright excitons 5 . Due to large exciton binding energy, they dominate an optical response of TMDC materials even at room temperature 6,7 . Also, thanks to peculiar Coulomb interaction screening in 2D 8-10 TMDC excitons have a non-Rydberg energy spectrum 11 . Combined with the remarkable compatibility of these materials with various semiconductor/dielectric platforms, this makes them promising candidates for the development of various nanophotonic components 12 , including logic circuits 13,14 , phototransistors 15 , and light sensors 16 as well as lightproducing and harvesting devices 17,18 . Studies of defects in TMDC monolayers and increased confinement also led to efficient single-photon emitters 19,20 .
The interaction of excitons with free electron gas, which can be present in a TMDC monolayer due to external doping, substantially modifies its optical response. The position of the excitonic line is shifted due to the effects of the bandgap renormalization [46][47][48][49][50] . In addition, complimentary screening and effects of Pauli blocking change the structure of the excitons, making them less compact and reducing the corresponding oscillator strength and Rabi splitting 51 . Other absorption peaks appear due to the brightening of an intervalley exciton 52 . The interaction between excitons and an electron gas leads to the appearance of charged quasiparticles being trions 4,53-57 or exciton polarons [58][59][60][61][62] . The presence of such quasiparticles can qualitatively modify the light-matter coupling and completely reshape the polaritonic spectrum 63,64 . Moreover, various nonlinear effects 62,64-73 arise from the phase-space filling, hybridization with light exciton-exciton and exciton-electron scattering, which offers a route towards polariton condensation 74 and the emergence of spontaneous coherence in monolayers strongly coupled to light 75 .
Here, we develop the quantitative microscopic theory for exciton-polariton states in TMDC monolayers with different doping levels. We consider a TMDC monolayer placed in an optical microcavity (Fig. 1). We solve the three-body problem by exact diagonalization of the Hamiltonian in a three-particle basis set 52,[76][77][78][79][80] , where an eigenmode solution is available at small and high free electron concentrations. We observe that an additional dressed-exciton mode appears at small doping, apart from the dominant exciton and trion polariton contributions. Qualitatively, its formation can be explained as a result of the brightening of an intervalley exciton in the presence of a Fermi sea, where free electrons contribute an additional momentum necessary for the optical recombination of the intervalley exciton hole with the doping electron. This couples strongly to the cavity mode and leads to a dipolariton-like spectrum. We observe a strongly modified spectrum at large doping at high frequencies, attributing it to the exotic excited-trion state. Our work lays the solid background for understanding the spectrum of TMDC monolayer excitations at the strong light-matter coupling.

A. Single-particle states
We start by describing the microscopic theory for polaritons in the electron-doped transition metal dichalcogenide monolayers. We consider a generic description valid for different TMDC materials. In our calculations, we use MoS 2 as an example. First, we calculate the single-particle states of a single TMDC monolayer. For this, we use the massive Dirac model 81 that is well suited to describe the single-particle band structure in the vicinity of K and K = −K points of the Brillouin zone. The The purple curve shows the absorption spectrum in the limit of zero doping level, where a single peak corresponds to 1s exciton. The blue curve corresponds to the absorption spectrum of three-particle states with τ sz = −1/2 at the level of doping equal to EF = 10.56 meV. The red curve shows the same dependence for τ sz = +1/2 three-particle states.
massive Dirac Hamiltonian reads as (h = 1) where ∆ is the bandgap, v F is the Fermi velocity, k x,y are the electron crystal momentum component, τ takes into account for the valley degrees of freedom and λ c,v are the SOC (spin-orbit coupling) constants. In addition, σ x,y,z are pseudospin Pauli matrices acting in the band subspace,σ ± = (σ x ± iσ y )/2 are raising and lowering operators, whileŝ 0 ,ŝ z are identity and spin Pauli Z matrix in the electron spin subspace, respectively. The role played by the SOC is included in the last term of Eq. (1). We note that the Hamiltonian (1) matches the density functional theory calculations at low energies 2,82 . Specifically, trions are formed by the states near the K and K points and the contribution from states away from the corners of the Brillouin zone are negligible 52 . Despite the relative simplicity, the model captures all relevant phenomena, including the interplay between the spin and the valley degrees of freedom.
To account for the effects of the dielectric environment, we introduce the bandgap dependence on the sub- strate and monolayer properties. The bandgap renormalization enters as the scissor shift δ for the bandgap ∆ = ∆ 0 + δ that modifies the bare monolayer bandgap ∆ 0 . The Fermi velocity is then adjusted as v F = ∆/2m to preserve the effective mass amplitude 83 . For the case of MoS 2 we consider the effective masses of electrons and holes to be 0.52m 0 (m 0 is the free electron mass), bandgap ∆ 0 of 2.087 eV, λ c of 1.41 meV, λ v of 74.6 meV, the bulk dielectric constant ε 0 of 14, the dielectric thickness of 6Å, and scissor operator δ equal to 0.64 eV.

B. Many-particle states
We introduce the interaction between electrons and holes and diagonalize the many-body three-particle Hamiltonian. In TMDC monolayers, the intravalley interaction is described by the screened potential in the Rytova-Keldysh form. 8,84,85 The Fourier transform of the screened Coulomb between particles in the same valley reads W (q) = V (q)/[ε + (1 + r 0 q)], where V (q) = 2πe 2 /q is the bare Coulomb potential in 2D, and q is the magnitude of exchanged wavevector. We use the average substrate permittivity ε + = (ε 2 + ε 1 )/2 and the screening length r 0 = ε 0 d/2. 86,87 The screened Coulomb potential for the intervalley processes reads W (q) = V (q)/ε 0 . 88 We neglect the wetting layer, vacuum spacing between the monolayer and substrate, which otherwise will modify the potential. 88,89 Excitons and trions in the presence of doping. To calculate the spectrum of three-body quasiparticles in the monolayer and find the corresponding eigenstates, we define operatorsT † ν (T ν ) that create (annihilate) states ν with two electrons in the conduction band (labeled as c 1,2 ) and a hole in the valence band (labeled as v). The corresponding basis is spanned by linear superpositionŝ whereâ † b is the fermion creation operator for a state in band b, and A ν,c1c2v denote the amplitudes. Note that c 1,2 are referring to the conduction band, and we use the usual anticommutation relations {â † c1 ,â † c2 } = δ c1,c2 . The valence band is completely separated, and we assume the relation {â c1,2 ,â v } = 0 to hold. At this point we stress that the introduced three-body operator can describe both a tightly bound state of two electrons and holes, as well as an electron-hole pair (exciton) in the presence of another electron.
Using the generalization of the Tamm-Dancoff approach 52,76-80 , we decompose the full many-body HamiltonianĤ in the three-particle state basis and solve the corresponding eigenvalue problem (see Methods for the details). We obtain the set of low-energy eigenstates and corresponding energies to high precision, which allows describing the optical properties of the system. To classify the states in trionic basis, we look at the product τ s z , with τ = τ c1 + τ c2 − τ v being the total valley index and s z = s c1 z +s c2 z −s v z being the total spin index. The total spin s z , in general, can have magnitudes 1/2 and 3/2. We restrict our analysis to s z = ±1/2 spin projections as states that participate in optical recombination processes and neglect spin-3/2 states. The case |s z τ | = 1/2 is a necessary condition for the trion state to be bright. In the numerical procedure, we also account for the symmetries stemming from the conservation of total spin s z and total momentum k = k c1 + k c2 − k v .
We model the presence of doping by introducing Pauli blocking through k-space discretization. The k-space grid is arranged as a mesh of N × N points, with the area of the primitive cell being Ω 0 . The doping density can be introduced as n = g v g s /(Ω 0 N 2 ), where g s and g v are the spin and valley degeneracies, respectively. The corresponding Fermi energy is equal to is the interval between k-points, and a is a lattice constant for the hexagonal lattice. The described approach allows obtaining both trion and exciton states in the presence of doping from the three-particle Hamiltonian, and in the zero doping limit (N → ∞) recovers the results from the Bethe-Salpeter equation for the excitonic states 52 . Therefore, we treat trionic and dressed excitonic states on equal footing and study polaritonic modes in the broad range of doping levels.
Together with eigenenergies of qus siparticles in the TMDC monolayer, we obtain the eigenstates and access various system observables. Specifically, we calculate dipole matrix elements D ν for interband transitions for the three-body state ν. Dipole moments are labeled by the corresponding band and read Absorption is shown as a function of detuning for increasing Fermi energy (left-to-right). We consider σ+ polarization for the cavity mode. Color bars define the intensity in arbitrary units. In all the spectra, from low (3.26 meV) to high (29.34 meV) Fermi energy, we see the emergence of trion modes at ∼ 2.1 eV with considerable Rabi splitting (∼ 10 meV), being significantly larger than the broadening γ = 1 meV . The exciton-polariton features the appearance of the additional mode inside the polaritonic gap. At high free electron concentration, the additional splitting is attributed to the higher-order trion state.
where d vc is the single-particle dipole matrix element for bands v and c, and δ c,c1 is a Kronecker delta. It should be noted that from the above expression, it follows that trions with a finite moment also have coupling with optical cavity mode. However, in our work, based on the assumption of zero temperature, we have calculated trions only with a moment equivalent to the moment of the centers of the K and K = −K valleys 80 .
We present the results for a TMDC monolayer (MoS 2 ) in the case of far-detuned optical mode (Fig. 2). The three-particle states are obtained by solving Eq. (9) at different levels of doping, i.e. Fermi energy E F . In Fig. 2(a), we plot the first 64 eigenenergies as a function of E F . We observe two dominant sets of modes at low energies for both spin-valley indices τ s z = +1/2 and τ s z = −1/2. The lowest set is attributed to trion-bound states. Trion energies, at around 2.15 eV at small free electron concentrations (E F = 2.93 meV), decrease sublinearly with increasing E F to 2.1 eV at large concentrations. With the increasing doping level, the redshift of the trion line can be attributed to the increase of the excitonic polarizability, which is necessary for the bind-ing of an additional electron. The size of the markers denotes the magnitude of dipole moments and shows the gradual growth of the trion oscillator strength. As expected, the excitonic modes, which are located ∼ 25 meV above trions, reveal opposite behavior with respect to trions, as they have the largest optical dipole moment at low concentrations. At the same time, with increasing E F we observe gradual redistribution of the oscillator strength between different excitonic modes as well as their blueshifts [ Fig. 2(a)], which can be attributed to the effects of screening and Pauli blocking 51,59 . We also note the appearance of a higher-lying band of states that gain significant optical dipole moments.
To plot the photoluminescence (PL) in the weak coupling regime, we assign the final lifetime for each quasiparticle and plot the PL spectrum as a set of Lorentzians [ Fig. 2(b)]. Choosing a large doping limit, we observe the strong trion mode, the shifted exciton peak, and an additional brightened quasiparticle at 2.155 eV 80 .

C. Light-matter coupling
Next, we consider a planar photonic microcavity system with a doped monolayer of TMDC placed in the antinode (Fig. 1). Confined cavity photons tuned in resonance with TMDC modes will couple strongly to excitonic and trionic optical transitions, resulting in a polaritonic spectrum with substantial modifications. Recently, the question of non-perturbative ab initio simulations of exciton-polariton in TMDCs was addressed in Ref. [94]. To describe TMDC polaritons at increasing doping, we use the diagrammatic approach. Starting from the bare Green's function of cavity photons G 0 (ω, q) with wavevector q, we obtain polaritonic from the photonic Green's function dressed by excitation in the system. This process is shown diagrammatically in Fig. 3, where multiple absorption and re-emission of a cavity photon lead to the mode hybridization 98,99 . Notably, this is a non-perturbative approach where the dressed propagator can be found exactly by solving the corresponding Dyson equation. In the algebraic form the solution reads where and ω cav (q) = ω 2 cav + c 2 q 2 with ω cav being the resonant frequency of a cavity at normal incidence, γ cav the broadening of the photonic linewidth due to a finite transmission coefficient of the cavity mirrors. The self-energy Σ(ω, q) describes material excitations and accounts for the interaction vertices defined by the strength of lightmatter interaction. It reads where Ω ν are Rabi frequencies of transitions, ν (q) are their resonant frequencies, and γ ν are non-radiative broadenings, depending on the quality of a sample. In total, we account for d = 64 transitions coming from the exact diagonalization and presented in Fig. 2(a). The broadenings are chosen based on typical experimentally achievable parameters. The light-matter interaction is defined by the optical dipole moments of transitions and the amplitude of the electric field (or vector potential) for the cavity mode. The associated Rabi frequencies are where A is the cavity area, and D ν are dipole matrix elements of the transitions given by Eq.  that the confinement of quasiparticles, as well as their physical origin, are already accounted for in D ν . Thus, this leads to the final result does not depend on the area of the sample A, but instead depends on the carrier concentration (see discussion in Ref. [? ]). For the dominant coupling to the exciton mode, the Rabi frequency is enhanced due to reduced exciton volume (area), independent of the number of particles. For the case of trion-dominant modes, the Rabi frequency is enhanced at increasing electron concentration (being area independent), which will analyze directly from the microscopic calculations.
An important consequence of cavity photon dressing by TMDC excitations corresponds to hybridization of these modes, as can be seen from the modification of poles of G(ω, q) for couplings |Ω ν,q | that dominate broadenings γ ν . The imaginary part of G(ω, q) then represents the absorption spectrum of the polaritonic system (up to scaling constant). 99 In the next section, we proceed to calculate the polaritonic spectrum and discuss the results.

RESULTS AND DISCUSSION
We consider a cavity polariton setup (see the sketch in Fig. 1) where excitons and trions couple to the cavity mode with a tunable frequency. The system can be conveniently realized in open cavities, where 0D cavity resonance can be shifted by cavity length change with a piezodrive, which allows coupling to excitations at energies ∼ 2 eV and resulting anticrossing of various modes. Using the developed theory and extensive numerical calculations, we present polaritonic spectra for varying TMDC doping. We concentrate on the gated MoS 2 monolayer case and cover a broad range of densities. It is convenient to present the doping level using the Fermi energy E F for the free electron gas. We start with an intrinsic value of E F = 2.92 meV and move up to larger doping levels up to E F = 29.3 meV. We plot the absorption spectra of the polaritonic system in Fig. 4, showing its modification at increasing E F . Here we assume mode broadening of γ = 1 meV being characteristic to modern setups. 63 Due to light-matter coupling with several modes simultaneously, we observe several distinct polaritonic modes with separate anticrossings.
At low doping in [ Fig. 4(a)], we observe two clear anti- crossings between the cavity mode and three-body excitations, corresponding to bound trions and 1s exciton states renormalized by the interaction with free electrons. 69,95 While the excitonic response dominates at higher energies, we note the significant oscillator strength for the low-energy trion modes and a clear sign of anti-crossing. Taking the lowest bright trion (attractive polaron) as an example, we check that the effective Rabi energy for this mode [Eq. 7] is equal to 4.2 meV at E F = 3.26 meV, being larger than 1 meV broadening. At the same time, we note that due to the multi-mode structure of our system the splitting between modes has to be estimated directly from the spectrum, as we show later in Fig. 5. At intermediate doping in [Fig. 4(b,c)], the most pronounced effects are: (1) the increase of trion-polariton splitting; (2) qualitative change of exciton-polariton coupling. For point (1), we track the energy splitting at resonant detuning and plot it as a function of the square root of the Fermi energy (and electron concentration), There is a known linear scaling of the trion-polariton splitting 64,91,102,103 , and we confirm the expected scaling at small and intermediate doping levels [see details in Fig. 5(a)]. Specifically, we show how the trion energy splitting grows, comparing numerical results [red circles and blue curve in Fig. 5(a)] with the linear scaling [ Fig. 5(a), yellow line]. For point (2), we note that additional mode that brightens with the increase of E F was noted in Ref. [52]. As it is located at the middle of the polariton gap, the anticrossing takes shape characteristic to dipolaritons, and features flat dependence on the cavity mode energy [ Fig. 4(c,d)]. This feature enriches the polaritonic spectrum, which was overlooked in the previous theoretical analysis 62,95 .
As doping level grows to the point when Fermi energy becomes of the order of trion binding energy, i.e. E F ∼ 20 meV, additional polaritonic modes emerge from excitonic modes dressed by the Fermi sea. We observe several modes that gain the oscillator strength at E F = 10.56 meV [in Fig. 4(e), see spectral features at around 2.18 eV]. These are attributed to the excited trion states, as follows from the charge distribution analysis presented below and in Fig. 6. Results obtained from our microscopic theory in the large doping limit have several remarkable features [ Fig. 4(f)]. The trion splitting is large in absolute values and reaches 18 meV value at E F = 29.34 meV [see spectrum cross-section in Fig. 5(b)]. At resonance with the trion mode, the crosssection shows two dominant peaks and weakly coupled excitonic modes at higher energies. The splitting shows a sublinear dependence with √ E F , which can be attributed to the increase of the role of the effects of the Pauli blocking and redistribution of the oscillator strength. However, we stress that strong coupling to trions remains at all doping levels.
To understand the nature of quasiparticles, we plot the charge distribution for different states (Fig. 6). This plot shows the contributions of single-particle states with different k into the full three-particle wavefunction for four many-body states. All states in Fig. 6 are bright because all four quasiparticles have direct electron-hole pairs with equal spins. The first state in Fig. 6(a) represents the trion state, which can be seen from the symmetrical conduction bands occupation. The latter manifests the equivalence between two electrons in the threeparticle state. The second state in Fig. 6(b) represents the momentum-indirect exciton state. Indeed, we observe a fully localized exciton-state for electron and hole being in K and K valleys, respectively, while the quasifree electron is located in the K-valley (described by small charge density spreading in the reciprocal space). We also note a signature of Pauli-blocking in the valence band occupation at K point. The third state, shown in Fig. 6(c), is the direct exciton with the Pauli-blocking signature in the valence band occupation, plus quasi-free electron in the same valley. Finally, the fourth state in Fig. 6(d) energetically is located above three previous states but does not exhibit Pauli blocking signature in the valence band occupation and has an s orbital-like occupation. It features the charge distribution that corresponds to trion modes. We thus conclude that it is an optically active 2s trion state, which was experimentally observed 92 . We have shown only a few states, which correspond only to the τ s z = +1/2 states. However, there are much more optically active three-particle states. For example, there are at least three optically active fully symmetrical trion states 93 . We anticipate similar behavior for the indirect/direct excitons and the excited trion states.
An essential consequence of strong coupling in MoS 2 monolayers is the qualitative change of the ground state properties. The three-body quasiparticles of a bare monolayer include four trionic modes, with the dark mode being energetically favorable [see Fig. 7(a)]. Qualitatively, this ground state corresponds to the suppressed photoluminesce from the uncoupled system. Similar physics is observed in carbon nanotubes, where tightly bound excitons have zero optical dipole moments 107 . The situation changes dramatically when the TMDC monolayer is placed in the optical cavity, where strong coupling modifies energies and redistributed dipole moments [ Fig. 7(b)]. In this case, a pair of polariton modes formed from the bright trion, and the lower one can go beyond the energy of a dark trion uncoupled to light, as shown in Fig. 1. As a consequence, the ground state of the system now becomes optically active.

CONCLUSIONS
In conclusion, we have developed a microscopic theory that describes strong light-matter coupling in doped TMDC monolayers in a wide range of free electron concentrations. The theory involves a numerical calculation of a set of three-body excitations, which are then coupled non-perturbatively to the cavity mode. We have calculated polariton spectra using this approach, which reveals the rich structure of the emerging composite light-matter quasiparticles. Our results confirm the robust nature of trion polaritons that holds in a large range of electron concentrations. At intermediate doping levels, we observe the three-mode structure of exciton-polaritons, where the presence of the additional mode in the middle of the polariton gap changes the anticrossing behavior qualitatively. Also, we predict that in polaritonic samples with large doping, the oscillator strength of the excited trion state becomes significant and may drive the system into a strong coupling regime with exotic charge distribution. Finally, we have demonstrated that increasing the electron concentration changes the nature of the system's ground state from dark to bright. Our results are crucial for building electrically tunable nonlinear polaritonic devices based on 2D materials.

METHODS
To calculate the energy spectrum and wavefunctions of many-body state we use the generalization of the Tamm-Dancoff approach 52,76-80 . We start by decomposing the full many-body HamiltonianĤ in the three-particle state basis [Eq. (2)]. Next, we solve the corresponding eigenvalue problem. In the trion basis, the full many-body Hamiltonian consists of several parts,Ĥ =Ĥ0 +Ĥc,c +Ĥc,v, withĤ0 being the free-electron Hamiltonian,Ĥc,c the Coulomb coupling between electrons, andĤc,v the Coulomb coupling between electrons and holes. Explicitly, these terms read H = H0 + Hcc + Hcv, (8) H0 = (εc 1 + εc 2 − εv)δ c 1 c 1 δ c 2 c 2 δ vv , where W ab a b = W (ka − k a ) a |a b|b and V ab a b = V (ka − k a ) a |a b|b are screened and bare Coulomb matrix elements, respectively. εc,v are free particle energy terms for the conduction (c) and the valence (v) band.
Being eigenstates of the full many-body Hamiltonian, one can find the elements of A by solving the eigenvalue problem where matrix elements H c 1 c 2 v c 1 c 2 v are described in Ref. [52. A direct solution of Eq. (9) is an extremely challenging task due to the high dimensionality of the trion Hamiltonian, which is proportional to NvN 2 c N 2 k , where Nv is a number of valence bands, Nc is a number of conduction bands, and N k is a number of k-points. For the converged calculations, the linear dimensions of the Hamiltonian matrix are more than ∼ 10 6 states. In the case of two-body exciton Hamiltonian, the dimensionality scales significantly slow as NvNcN k . However, the trion Hamiltonian matrix is sparse, with about 1% of the nonzero matrix elements. Therefore, we use the Arnoldi algorithm 109 , which is implemented in ARPACK 110 , for obtaining low-energy eigenstates.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
All the codes were implemented using home-built programs, which use standard open-source software packages.