Strong Coupling of Two-Dimensional Excitons and Plasmonic Photonic Crystals: Microscopic Theory Reveals Triplet Spectra

Monolayers of transition metal dichalcogenides (TMDCs) are direct-gap semiconductors with strong light–matter interactions featuring tightly bound excitons, while plasmonic crystals (PCs), consisting of metal nanoparticles that act as meta-atoms, exhibit collective plasmon modes and allow one to tailor electric fields on the nanoscale. Recent experiments show that TMDC-PC hybrids can reach the strong-coupling limit between excitons and plasmons, forming new quasiparticles, so-called plexcitons. To describe this coupling theoretically, we develop a self-consistent Maxwell-Bloch theory for TMDC-PC hybrid structures, which allows us to compute the scattered light in the near- and far-fields explicitly and provide guidance for experimental studies. One of the key findings of the developed theory is the necessity to differentiate between bright and originally momentum-dark excitons. Our calculations reveal a spectral splitting signature of strong coupling of more than 100 meV in gold-MoSe2 structures with 30 nm nanoparticles, manifesting in a hybridization of the plasmon mode with momentum-dark excitons into two effective plexcitonic bands. The semianalytical theory allows us to directly infer the characteristic asymmetric line shape of the hybrid spectra in the strong coupling regime from the energy distribution of the momentum-dark excitons. In addition to the hybridized states, we find a remaining excitonic mode with significantly smaller coupling to the plasmonic near-field, emitting directly into the far-field. Thus, hybrid spectra in the strong coupling regime can contain three emission peaks.


I. INTRODUCTION
The light-matter interaction strength in transition metal dichalcogenide (TMDC) monolayers has been reported to be extremely strong [1], e.g., as demonstrated by absorption rates of up to 10% in the visible spectrum [2,3].Such a high absorption is particularly noteworthy given the two-dimensional nature of these materials, which possess a thickness of less than 1 nm.In addition to featuring a direct bandgap [2], TMDC monolayers support in-plane exciton formation due to their two-dimensional structure [4].Excitons (bound electronhole pairs) therefore dominate the optical spectrum below the band edge [5].In addition, the remarkably thin nature of TMDC monolayers results in an increased sensitivity to surrounding materials.Consequently, the atomically thin materials can easily be influenced by various factors, such as the choice of the substrate material, defects [6,7], and functionalization [8], e.g., with molecules [9,10] and heterostructure configurations [11,12].
In contrast, the optical response of metal nanoparticles (MNPs) is dominated by localized plasmons which are collective electron oscillations formed within the metal conduction band [13].A special feature of MNPs is a significant amplification of the electric near-field, which additionally allows for manipulating the electric field on dimensions far below the diffraction limit [14,15].Arranging MNPs as meta-atoms in a crystal structure yields a plasmonic crystal (PC) with extraordinary strong light-matter interaction.The localized plasmons couple with the electric field and form plasmon-polaritons * lara.greten@tu-berlin.de which can propagate within the crystal [16,17] and sharpen the single particle plasmon resonance.The optical properties of PCs strongly depend on a variety of parameters like lattice structure and nanoparticle shape.By manipulating these, it is possible to tune the optical properties of the crystal over a wide range [16,[18][19][20].The strong tunability and enhancement of the electric field make periodic plasmonic structures appealing for light harvesting and non-linear optics, yielding, e.g., applications for nanoscale lasing [19,21] and advanced optical spectroscopy [22,23].
Figure 1.Sketch of the hybrid system: the 2D semiconductor (TMDC) is covered by a square-structured 2D PC of metal nanodisks with the example of gold.The structure is periodic and infinite in the xy-plane and embedded in a surrounding medium with constant permittivity ε.

arXiv:2309.09673v1 [cond-mat.mes-hall] 18 Sep 2023
It has also been shown that the absorption of graphene is significantly enhanced by depositing plasmonic nanostructures near the graphene layer [24,25].An interaction between localized surface plasmons of a single MNP and excitons in the semiconductor plane can reach the strong coupling regime [26][27][28][29][30][31][32][33] and absorption rates up to 90% [34].This motivates that combining the particular properties of PCs with the environment-sensitive semiconductor monolayer promises a further increase of the light-matter-interaction of the excitons in the TMDC.
In the context of our study, we adopt a definition of "strong coupling" where the interaction strength between the two systems [35,36] exceeds the total losses of the system, where the former is quantified by an effective Rabi splitting Ω eff ≈ 2g eff > γ ex + γ pl [31,32,37] with g eff an effective excitonplasmon interaction strength.In a classical theory, this is usually termed normal mode splitting, though the underlying dissipative modes are quasinormal modes with complex eigenfrequencies, which are well characterized by Maxwell's equations [27].This definition is particularly relevant for excitonplasmon coupling as the losses in the plasmonic components are significantly larger than the losses in the excitonic component [27].In this work, we study the impact of a twodimensional (2D) PC on the excitonic dynamics in TMDCs.A sketch of the hybrid system is depicted in Fig. 1.The semiconductor is located parallel to the xy-plane.The 2D PC of gold disks (square array) is placed on top of the TMDC.Experimentally, a similar configuration for relatively large nanoparticles with an in-plane radius of more than 50 nm was already realized in Refs.[38][39][40][41], where TMDCs are coupled to 2D PCs with different lattice structures.
To theoretically study the exciton-plasmon coupling, we develop a self-consistent Maxwell-Bloch theory for the hybrid structure in Sec.II.We first give a short review on the solution of the Maxwell's equations (Sec.II A), the description of the TMDC monolayer with the excitonic Bloch equation (Sec.II B) and the 2D PC using Mie theory (Sec.II C).This allows us in Sec.II D to use the plasmonic dipole density within the PC as a source for the excitonic dipole density of the TMDC, by inserting it into the excitonic Bloch equations.The resulting formulas are Bloch equations that describe the dynamics of excitons in the presence of an electric field and the field-mediated interaction.We highlight the qualitatively different behavior of momentum-dark excitons, that couple to the plasmon-enhanced electric near-field, and originally bright excitons with vanishing center-of-mass momentum.Finally, an expression for the emitted electric field is derived from which we can deduce transmission, reflection, and absorption of the hybrid structure when using a plane wave excitation.
Using this theoretical framework, we perform numerical calculations to obtain the absorption spectra of the hybrid system in Sec.III.Our analysis resembles the broad range of coupling regimes, varying from weak to strong.In the strong coupling case, we find an additional peak at the unperturbed exciton energy, which exhibits a strong dependence with the temperature and coupling strength.Our theory suggests that this additional mode stems from the originally bright excitons that do not participate in the strong coupling to the plasmon-enhanced electric near-field.We find that this excitonic mode is uniformly distributed in real space, thus also located in the immediate vicinity of the MNP.The observed strong coupling between TMDCs and MNPs or PCs is based on the coupling between plasmon and momentum-dark excitons.We note that there is ongoing discussion regarding exciton-plasmon coupling in nanoshells, with some studies also suggesting the existence of an undisturbed excitonic mode [42].To the best of our knowledge, this mode has not been observed in experiments for nanoshells [43].However, our results agree with recent experiments [40] which observe the presence of such an additional excitonic peak for TMDC-PC hybrids.

II. THEORETICAL MODEL
In this section, we develop a theoretical framework that describes the interaction between TMDC excitons and PC plasmons via the radiation field.

A. Maxwell's Equations
Starting from Maxwell's equations, assuming a nonmagnetic and isotropic medium, the wave equation for the electric field reads with the source which is valid for a freestanding sample with dipole density P(r, t) embedded in a homogeneous, isotropic and nondispersive dielectric environment with permittivity ε.A solution of Eq. (1) can be formally obtained via the scalar Green's function G(r, r ′ , t − t ′ ), where the scalar Green's function depends on the boundary conditions.Both layers of the hybrid structure, cf.Fig. 1, are assumed to be aligned parallel to the xy-plane with a discrete translation invariance for the PC.We treat the z-direction separately and apply a Fourier transform in the in-plane coordinates x, y and time t, where we introduced the new Green's dyadic G q ∥ (z, z ′ , ω) that converts the dipole density at the position z ′ into an electric field and describes its propagation from z ′ to the observation position z.The term G q ∥ can be derived analyti-cally [44,45], Equation ( 5) gives the Green's dyadic in Cartesian coordinates regarding the z-direction whereas the in-plane components are expressed independent of a special basis for q ∥ , which will be specified in Sec.II D. The scalar Green's function for a constant surrounding permittivity is where we abbreviate the absolute value of a vector by, e.g., Next, using the thin film approximation for the individual layers [46,47], i.e., with TMDC (ex) and PC (pl), then allows to find an algebraic solution of Maxwell's equations in Fourier space.Adding the incident field E 0 q ∥ (z, ω) as a solution of the homogeneous wave equation, Eq. (1), yields which allows us to compute the transmission and reflection.Thus, we can connect the dipole density to experimentally measurable far-field signals.
We assume a perpendicular plane wave excitation that is propagating in the positive z-direction, with an amplitude perpendicular to the z-axis E 0 ⊥ e z .In the far-field limit, only a q ∥ = 0 Fourier component occurs, see Eq. (36) .The transmission is given by and the reflection is The absorption is then easily obtained from

B. Excitonic Dipole Density
To describe the response of the TMDC excitons to the electric field, we define the macroscopic 2D dipole density [46,48] with the dipole moment d ξ , ξ = +/− corresponding to the K + and K − valley, carrying the circular dichroism [49].The strength of the dipole moment d ξ is taken from DFT calculations [49]; the term φ r ∥ =0 accounts for the value of the 1s excitonic wave function in real space at r ∥ = 0 and is obtained from the solution of the Wannier equation [50,51], incorporating the dielectric environment in the Rytova-Keldysh approximation [52].We concentrate on the lowest 1s excitonic state, since it is energetically separated from higher transitions [5].A table with all parameters for MoSe 2 can be found in Appendix A. The excitonic transition amplitude is denoted as p ξ q ∥ and obtained via its Bloch equation in the Fourier domain [53], (14) with the dipole tensor and the excitonic transition in vector notation, The zero in the z component is added for convenience, to only allow for the appearance of square matrices.A detailed derivation is provided in Ref. [54], where the rotating wave approximation (RWA) is utilized.This approximation imposes a constraint on the light-matter interaction strength, which must be significantly smaller than the system energies.Therefore, the range of exciton-plasmon coupling strength that can be accurately described in this work is limited to g eff /E pl < 0.1 which is below the ultra-strong coupling regime [17,55,56].
The temperature-dependent dephasing rate γ(T ) accounts for non-radiative decay, which typically results from excitonphonon interactions and is calculated microscopically [51], where c 1 , c 2 , and the averaged phonon-energy Ω are given in Appendix A. All radiative corrections such as lineshift/splitting and broadening are incorporated via the self-consistently calculated electric field E q ∥ (z ex , ω) given in Eq. ( 8).The left-hand side of Eq. ( 14) accounts for the dispersion of excitons with center-of-mass wave number q ∥ (also referred to as inplane momentum), the excitonic mass M , the excitonic transition energy E ex of 1s-excitons with vanishing momentum.The parabolic exciton dispersion is illustrated in Fig. 2. In Eq. ( 14), an excitonic transition p q ∥ can be excited by an electric field E q ∥ with the same in-plane momentum.The excitonic transition p q ∥ is called bright if it can be excited from the far-field.However, propagating solutions of the wave equation for the electric field, Eq. (1), are only possible for real values of k q ∥ , cf.Eq. ( 6), meaning for small momenta q ∥ .We refer to excitons with higher momenta as momentum-dark, as they are inaccessible via far-field illumination.
To connect the solution of Eq. ( 14) with the Green's dyadic formalism solving Maxwell's equations in the circular basis (+, −), we adjust the notation of Eq. ( 13) to Due to the circular dichroism of TMDCs, the dipole tensor is diagonal in a circular polarized basis, d σξ = dδ σξ .

C. Plasmonic Polarizability
To describe the response of MNPs, we use Mie-Gans theory [57,58], providing a parametric frequency dependent polarizability α(ω).The optical response of every MNP is approximated by a point dipole p j at the lattice position R j , with where E(R j , ω) is the electric field, excluding the field generated by the MNP itself which is already incorporated in α(ω).
The quasi-static polarizability α qs for spheroids is given by, e.g., Mie-Gans theory [58,59].It can be used for a quasistatic description, if the extensions of the MNPs are significantly smaller than the wavelength of the exciting light-field r i ≪ λ.In Cartesian coordinates, with the MNP axes along the corresponding semi-axes (i = {x, y, z}), the quasi-static polarizability is diagonal with the entries where r i is the half-axis of the MNP along the corresponding direction, ε Au (ω) is the permittivity of gold, ε the surrounding permittivity and L i the shape factor for oblate spheroids, with the eccentricity, To construct the optical response of the MNP-based 2D PC depicted in Fig. 1, we need to take the radiative coupling between the nanoparticles into account.Therefore, we consider corrections to the quasi-static Mie-Gans solution of the single particle polarizability in the modified long-wavelength approximation (MLWA) including correction terms for dynamic depolarization and radiation damping [60][61][62], with the radiative wavevector k q ∥ =0 = √ ϵ ω c , compare for Eq.(6).A full self-consistent solution would inherently include the radiation damping.However, in Eq. ( 21) the electric field generated by the MNP itself is excluded, making this correction necessary.Whether dynamic depolarization and radiation damping are important for a single MNP depends on the size of the MNP: for particle diameters significantly smaller than λ/2π, the static Mie-Gans theory, Eq. ( 22), would be a sufficient description.In Sec.III, we evaluate the theory for MNPs with in-plane radii of 30 nm, where corrections due to the finite extent become essential for the description of their optical properties.
For the gold-MNP permittivity ε Au (ω), as occurring in Eq. ( 22), a heuristic, analytical model that fits experimental data for bulk [63] is parameterized by [64], In the high-frequency limit, the permittivity ε ∞ of gold differs from unity since d-bands in noble metals are filled and provide a high residual polarization [13].Equation (27) incorporates a Drude-like intraband response of the conduction band and two interband transitions with the possibility to fit asymmetric line shapes in linear bulk spectra (last two terms in Eq. ( 27)).We combine the room-temperature permittivity given in Ref. [64] with a temperature and spectrally dependent linewidth [65][66][67][68] via the damping term Γ(ω, T ) as a sum of electron-electron and electron-phonon scattering using a Debye model, with [65][66][67][68]] A table with all parameters, including the Debye temperature Θ, for the permittivity of gold is given in Appendix A. We consider the temperature dependence of the linewidth only for the Drude contribution in Eq. ( 27), since it is the main contribution in the well-studied energy regime below the energy range where interband transitions become relevant, ℏω < 2.4 eV [68].For simplicity, we disregard a subtle redshift of ω p for low temperatures.Equation ( 28) does not consider radiative damping, since it is included by self-consistently solving Maxwell's equations.Equations ( 21)-( 27) fully describe the response of single nanoparticles to a self-consistently calculated electric field E, Eq. ( 8).
In contrast to the single MNP properties, in a PC, the field generated by all other nanoparticles has to be considered self-consistently.The discrete translational invariance in the xy-plane allows to collect the interactions with the other MNPs in a modified polarizability α * q ∥ (ω) in Fourier space.A self-consistent solution for the two-dimensional dipole density of the 2D PC in coupled dipole approximation is given by Ref. [69] and experimentally confirmed for gold nanodisks with r x = r y ≈ 60 nm and a lattice constant a ≈ 500 nm in Ref. [61].The corresponding dipole density reads where the sum includes all reciprocal lattice vectors g ∥ corresponding to Umklapp processes.The effective polarizability tensor α * q ∥ contains corrections due to interactions between the nanoparticles in the coupled dipole approximation restricted to MNP center-to-center distances a ≥ 3r xy .In this limit [70], the effective polarizability is given by normalized by the unit cell area A pl UC of the PC.The form factor accounts for particle interactions, with the 3 × 3 unity matrix 1.The integer j indexes the nanoparticles at the in-plane lattice vector positions of the 2D PC.Therefore, the form factor depends on the lattice structure and lattice constant a.If the momentum q ∥ equals 0 or a reciprocal lattice vector g ∥ , the form factor becomes diagonal with the entries [71], .
The angle dependency of the summands is given by X j i = (cos 2 θ j , sin 2 θ j , 0) for the entries i = x, y, z, respectively, with the polar angle θ j of R j ∥ .In the numerical evaluation of the form factor, Eq. ( 34), the sum is only slowly converging, leading to spurious oscillations.This unphysical result can be circumvented using Ewald's onefold integral transform [72], yielding a fast converging expression for the form factor.
The Umklapp processes in Eq. ( 31) are considered in the near-field of the sample, but not in the far-field (|z − z pl | → ∞).This is appropriate since dipole densities with in-plane momentum q ∥ provide purely evanescent electric fields if the wave number k q ∥ is imaginary, cf.Eq. ( 6).It follows that the g ∥ = 0 summand is the only propagating contribution if The absolute value of the minimal non-trivial reciprocal lattice vector is |g min ∥ | = 2π/a with the lattice constant a, as depicted in Fig. 1.Consequently, Eq. ( 35) can be rewritten as a condition for the lattice constant a compared to the wavelength λ of the incoming light, If we consider the sum over reciprocal lattice vectors in Eq. ( 31) as a diffraction phenomenon, Eq. ( 36) states that the lattice constant a is small enough that the diffraction pattern in the far-field only exhibits the main maximum and not any additional side lobes.However, for the self-consistent solution, we have to consider all Umklapp processes in Eq. ( 31), since both propagating and evanescent electric fields couple to the TMDC excitons which are located in the near-field of the PC.

D. Plexcitons
In the following, the electric field, Eq. 8, and the dipole densities, Eq. 13 and 31, are solved self-consistently to provide the near-field and far-field response of 2D TMDC-PC hybrids to an initially incident plane wave.To solve the set of equations, we insert the dipole densities of the individual layers, Eqs. ( 19) and (31), into the solution of the electric field, Eq. ( 8), and altogether in the excitonic Bloch equation, Eq. ( 14).From this procedure, we obtain a Bloch equation for the excitonic transition p q ∥ , where the left-hand side describes the free propagation of the excitons, cf.Eq. ( 14), and on the right-hand side the selfconsistent electric field occurs as a function of the dipole densities.
The first term appearing on the right-hand side of (37), after the left parentheses, E 0 q ∥ (z ex , ω), accounts for the undisturbed incoming electric field contribution at the TMDC position; the second term carries the renormalization of the exciton line due to radiative self-interactions; the third term accounts for the incoming field E 0 q ∥ (z pl , ω) scattered at the PC; and the fourth (final) term, is the electric field generated by the excitons and back-scattered from the PC.This term is a radiative near-field interaction that couples excitonic transitions with momenta q ∥ to those q ′ ∥ which are shifted by reciprocal lattice vectors g ∥ of the PC compared to the excitonic momentum on the lefthand side.This coupling results from the discretized translational invariance and will later be interpreted as a binding potential for originally spatially extended excitons.The terms two and four give rise to an intervalley coupling in Eq. ( 37), since their prefactor matrices are non-diagonal in the valley index.
Since for the solution of Eq. ( 37), we need to consider the observables near the sample, i.e., in a near-field coupling, we apply a quasi-static approximation in Eq. ( 5).We neglect propagation effects by setting ω to zero and otherwise only consider a parametric ω-dependency, e.g., in the permittivity ε Au (ω), the polarization densities or the electric field.This approximation is valid as long as [73] where the wavelength λ of the incoming light is compared to the distance between the two lattices.Thus, we obtain the simplified (quasi-static) scalar Green's function which shows an exponential decrease that depends on the distance between the source of the electric field at z ′ and the observer with position z.We indicate the quasi-static Green's function by dropping the dependence on the frequency ω.However, the Green's dyadic, Eq. ( 5), shows that the quasistatic approximation corresponds to neglecting ω 2 compared to the in-plane momentum q 2 ∥ which is only reasonable if In the following, we refer to electric fields and dipole densities that fulfill the condition in Eq. ( 40) as outside the lightcone.No propagating solution of the wave equations is possible, since the electric field is damped exponentially with the distance, cf.Eqs. ( 6) and (39).In contrast, in-plane momenta that do not fulfill the condition are inside the light-cone.The condition (40) is only valid for excitonic momenta q ∥ belonging to originally momentum-dark excitonic transitions.However, an inspection of Eqs. ( 8) and (31) shows that, for a perpendicular plane wave excitation, cf.Eq. ( 9), the plasmonic dipole density is only non-zero for in-plane momenta equal to a reciprocal lattice vector g ∥ of the PC.Since only the zeroth order scattering momentum g ∥ = 0 is a propagating solution (inside the light cone), it is allowed to apply the quasistatic approximation for all momenta q ∥ of the excitonic transition, except for q ∥ = 0. Therefore, we split Eq. ( 37) into a near-field Bloch equation (q ∥ > 0) and a radiative contribution (q ∥ = 0) and discuss them separately in the following (Sec.II D 1, II D 2, respectively).In particular, by inserting the Green's function for a homogeneous environment permittivity in quasi-static approximation, Eq. ( 39), one observes that the right-hand side of Eq. ( 37) is proportional to the in-plane momentum.The momentum-dark excitonic transition p q ∥ ̸ =0 (ω) (Sec.II D 1) does not couple to the radiative excitonic transition p q ∥ =0 (ω) (Sec.II D 2) and can be solved independently.

Near-Field Exciton-Plasmon Interaction
The resulting Bloch equations for q ∥ ̸ = 0 are identical to Eq. ( 37), with the quasi-static approximation applied to every Green's dyadic.To account for the circular dichroism of the TMDC excitons, we choose a circular polarized basis (e + , e − , e z ).Any quantity in Cartesian coordinates is transformed to the circular polarized basis by multiplying with the unitary (change-of-basis) matrix T circ , In the quasi-static limit, a diagonalization of Eq. ( 37) in circular polarized basis corresponding to the valleys K + and K − is given by a transformation T ϕ , Eq. ( 42), that depend on the polar angle ϕ of the in-plane momentum q ∥ ,    The transformation T ϕ corresponds to an in-plane rotation of the momentum space, orienting the momentum q ∥ along the V -axis.We find that the transformation into the (e U , e V , e z ) basis is the diagonalization of the in-plane contribution of the quasi-static Green's dyadic, Applying the rotation to the Bloch equations (37) finally allows a diagonalization process, ℏω where Eq. ( 43) possesses the undisturbed, parabolic (Ushaped) exciton dispersion for p U q ∥ .In contrast, we find with β = |d| 2 2ϵ 0 ϵ , yielding a conic (V -shaped) dispersion for p V q ∥ similar to results reported in Ref. [74].The excitonic self-interaction term, q 2 ∥ G ex−ex q ∥ , dominates the kinetic energy contribution for small q ∥ .The optical sources of p V q ∥ , cf.Eq. ( 44), occur on the right.Obviously, the near-field interaction of exciton and plasmon is fully encoded in p V q ∥ .The interaction with the PC is denoted by the contributions (ex ↔ pl) and provides a momentum and frequency dependent PC-induced excitonic self-interaction C ex↔pl q ∥ q ′ ∥ ω , namely, a plasmon mediated effective exciton-exciton interaction.For a concise notation without loss of information, we indicate the z-dependencies of the prefactor Green's functions in their superscripts.On the righthand side of Eq. ( 44), the existence of an incoming field E 0 q ∥ , that was scattered at the PC, results in the excitation of the conic excitonic transition.The interaction with the PC provides a coupling to in-plane momenta shifted by a reciprocal lattice vector g ∥ of the PC due to Umklapp processes depending on the PC lattice structure.The PC-induced excitonic selfinteraction reads and the source term that arises from the scattering of the incoming field at the PC becomes if we express E 0 q ′ ∥ (z pl , ω) in Cartesian coordinates.To achieve these expressions, derivatives of the Green's function with respect to z ′ have been carried out, cf.Eq. ( 5).
Since we were able to decouple the valleys of the excitonic transition {p + q ∥ , p − q ∥ } in Eq. ( 37) into a conic (optically driven p V q ∥ , Eq. ( 44)) and an undisturbed parabolic (p U q ∥ , Eq. ( 43)) contribution, we can now solve them independently.The parabolic Bloch equation for p U q ∥ immediately yields In contrast, in Eq. ( 44), the conic Bloch equation for p V q ∥ contains a coupling of excitons with different center-of-mass momenta, mediated by the PC.We project the excitonic transition p V q ∥ onto a set of eigenstates of the conic dispersion, The corresponding symmetric, non-Hermitian eigenvalue problem reads similar to Ref. [26], As a consequence of the non-Hermiticity we have to distinguish between left and right eigenstates v Lλ q ∥ , v Rλ q ∥ .To justify the projection onto these states, we verify the existence of solutions of Eq. ( 50) and their completeness numerically.By applying a suitable normalization, we ensure orthonormality using With this expression, we project the conic excitonic transition in the Bloch equation (44) on the new states according to Eq. ( 49).Substituting the momentum dependencies with the corresponding eigenvalue in Eq. ( 50) and taking advantage of the orthonormality condition between the eigenstates finally allows us to give a solution for the conic contribution of the excitonic transition To connect the excitonic transition with observables, we insert it in the dipole density of the TMDC layer.The contribution arising from the conic part becomes The full dipole density stemming from momentum-dark excitons trapped in the plasmonic potential is finally given by For momenta that do not correspond to the wave number of the incoming field shifted by a reciprocal lattice vector, the TMDC polarization vanishes due to the absence of a suitable excitation.

Radiative Interaction
In the previous section, the excitonic dipole density outside the light-cone was derived.However, to determine the farfield response, we need an expression for the dipole density within the light cone which is a source for propagating solutions of the electric field.Therefore, we evaluate the Bloch equation (37) for q ∥ = 0, i.e., the radiative Bloch equations, For clarity, we suppress the index of the vanishing in-plane momentum q ∥ = 0 in the notation for the radiative Green's dyadic.The previous provides us with all the necessary tools to solve this equation.Simplifying the individual terms on the right-hand side is demanding but straight-forward by explicitly evaluating the matrix products between the dipole elements, Eq. ( 15), the Green's dyadic, Eq. ( 5) and the PC polarizability α * q ∥ (ω).Finally, we collect all summands containing p q ∥ =0 on the left-hand side and multiply with the inverse of its prefactor.This procedure yields with the unity vector in polar coordinates when we express the excitonic transition p q ∥ =0 and the PC polarizability α * q ∥ in a linear polarized basis.
The denominator of Eq. ( 56) corresponds to harmonic oscillators at frequency ω with renormalized eigenenergies Re e(α * q ∥ =0 (ω)), (57) and the damping The ex-ex term accounts for the radiative dephasing of the TMDC excitons whereas the last term governs radiative interference phenomena between the TMDC and the PC.For a resonant TMDC-PC interaction (E ex = E pl ), the real part of the PC polarizability can be neglected Re e(α * q ∥ (ℏ −1 E pl )) ≈ 0. This yields an unchanged excitonic eigenenergy and the coupling to the PC mainly modifies the damping.The numerator includes from left to right: the direct excitation of excitonic transitions via the incident electric field, the plasmon mediated excitation from the incoming field that was first scattered at the PC, and finally the plasmon-mediated influence of momentum-dark excitonic transitions.From the derived excitonic transition, we deduce the dipole density via Eq.( 19).

Scattered Electric Field
To give explicit results for the electric field emitted by the TMDC layer E ex q ∥ , we multiply the dipole density, that was derived in the previous section, with the Green's dyadic.We find the momentum-dark dipole density to be an eigenstate of the quasi-static Green's dyadic, cf.Eq. ( 5).Therefore, we find In fact, also a z-component of the electric near-field E ex z q ∥ ̸ =0 (z, ω) occurs.However, it only causes a non-vanishing z-component of the PC polarization since α * q ∥ is assumed to be diagonal regarding its z-entries.This z-polarization of the PC is not transferred to the far-field for q ∥ = 0. To obtain the contribution of the PC to the far-field, E pl q ∥ =0 , we deduce the dipole density of the 2D PC, according to Eq. (31).It is As already justified above for dense PCs, cf.Eq. ( 36), there is no propagating solution of Maxwell's equations with an inplane wavenumber equal to a non-trivial reciprocal lattice vector.Therefore, we may drop the sum over g ∥ ̸ = 0 for the exciting field.As for the radiative TMDC contribution, the dyadic Green's function has to contain the full time dependency.Furthermore, since z ̸ = z ′ , the z-entries of the Green's dyadic vanish.The in-plane far-field plasmon contribution reads • Adding all contributions of the electric field according to Eq. ( 8), results in

Overview of the Theory
With the theoretical framework now fully established, our attention turns to shedding some light on the physical meaning of the derived equations.Figure 3 provides a schematic of the theory.
We excite the TMDC-PC hybrid via far-field illumination E 0 q ∥ , Eq. ( 9), represented by the sun-like schematic in Fig. 3.It directly acts on the PC plasmons, Eq. ( 60), and the bright excitons p q ∥ =0 , Eq. ( 56).The scattering respectively Umklapp processes of the incoming light at the PC provide access to momentum-dark excitonic transitions, Eq. ( 44), via the plasmon-enhanced electric near-field.These interactions are typical dipole-dipole near-field interactions, proportional to q ∥ e −q ∥ δz , which can be seen by explicitly inserting the quasi-static Green's function, Eq. ( 39), into the conic Bloch equation, Eq. (44).The coupling between the plasmon and the bright exciton behaves qualitatively different with the coupling strength significantly reduced, proportional to ω/c, cf.Eq. (56).The implications of these different coupling mechanisms are numerically evaluated in Sec.III.Finally, we note that, in the limit a → ∞, the developed theory is also applicable for a TMDC coupled to a single MNP without qualitatively changing the results.

III. NUMERICAL RESULTS AND DISCUSSION
We numerically evaluate the key equations ( 61), ( 59), ( 58), ( 54) and ( 50) to calculate the absorption, Eq. ( 12), of the hybrid structure.If not stated differently, we perform calculations for a MoSe 2 -PC stack with the geometry specified in Table I.The given parameters yield a PC with its resonance energy at E pl = 2.03 eV.We choose the lattice constant a in a way that the collective PC resonance sharpens the single MNP response that is governed by the bulk metal parameters and aspect ratio r z /r x .Material specific parameters for gold and MoSe 2 are given in Appendix A. To study different coupling regimes, the modification of the layer distance δz allows to directly access the strength of the near-field mediated interaction.Figure 4 shows the absorption of the hybrid over the detuning ∆ = ℏω − E pl between the electric field and the plasmon for the resonant case for exciton and plasmon (E ex = E pl ) for different exciton-plasmon distances δz.Modifying the distance δz qualitatively changes the hybrid's spectrum.All cases shown in Fig. 4 have been observed in experiments with TMDC excitons coupled to MNPs or PCs, yet lacking a microscopic theory, which we provide from the excitonic perspective.For MNP-TMDC distances where δz is approximately equal to or exceeds the effective extensions of the MNP, the exciton-plasmon coupling is weak.Thus, for δz ≥ 30 nm, the line shape is qualitatively preserved compared to the uncoupled case.However, the interaction induces an additional damping, increasing the linewidth and reducing the maximal absorption.This is because the interaction terms in the denominators of Eqs. ( 54) and ( 56) are imaginary, thus effectively acting like an additional damping to the phonon damping γ.This behavior is analogous to the classical coupled oscillator model [75,76] and has been experimentally observed between TMDC excitons and MNP plasmons in Refs.[29,77].
With further decreasing TMDC-PC distance, the absorption line shape changes.The absorption of the hybrid develops two maxima, while the absorption at the original, non-perturbed, resonance energy shrinks.The splitting between the formed maxima can be identified with an effective Rabi energy Ω eff .Its distance dependency is depicted in the inset of Fig. 4b, illustrating a (δz) −2 proportionality.To assign the hybrids to a particular coupling regime, the effective Rabi energy Ω eff ̸ = 0 is compared with the linewidth γ ex and γ pl of the individual constituents [76,78].
A coupling strength leading to a peak splitting that is small compared to the linewidths as e.g. in the δz = 17 nm case is still assigned to weak coupling and has been observed in experiments with TMDCs coupled to a MNP or a PC in Refs.[28,38,40,79].

Strong Coupling (δz = 11 nm)
Further decreasing the distance of TMDC and PC allows one to enter the strong coupling regime as Ω eff > γ pl + γ ex with peak splittings up to Ω eff ≈ 140 meV.In this case, mathematically, in Eqs.(50)(51)(52)(53)(54), a real part is added to the excitonic eigenenergy, thus changing the root of the denominator in Eq. ( 54), and consequently the plexcitonic resonance energies.Experimentally, strong coupling has been observed up to room temperature [29-32, 41, 77, 80-82] for TMDC-MNP hybrids with effective Rabi energies similar to the values shown in Fig. 4.
It would be worth also considering the possibility of ultrastrong coupling in a similar platform.However, this limit cannot be addressed by the developed theory in this work due to the constraint g eff /E pl < 0.1, since the rotating wave approximation is applied (Sec.II).

B. Low Temperatures
To study low-temperature effects, Fig. 5a depicts the absorption spectrum in the strong coupling case for coinciding exciton and plasmon resonance E ex = E pl at liquid nitrogen temperature T = 77 K.In addition to the typical strong coupling spectrum, we find a third peak at the undisturbed exciton resonance.At room temperature (Fig. 4a) this peak is not well resolved due to the increased damping.To investigate the dispersion of the three-peak spectrum, we numerically vary the exciton energy around the plasmon resonance in Fig. 5b.The color coding breaks the full absorption of the hybrid down to exciton (blue), plasmon (yellow) and the hybridized plexciton contribution (green).The individual plasmon and exciton dispersions are indicated by the horizontal line and the main diagonal, respectively.The avoided-crossing behavior of the upper and lower branch substantiates that the system is indeed in the strong coupling regime.However, the spectral position of the additional middle peak coincides with the unperturbed exciton energy E ex .The inset depicts the energy separation of the hybrid branches, where the minimum is the effective Rabi energy of exciton and plasmon.For comparison, the linear gray plot indicates the energy separation between exciton and plasmon without any interaction.

Interpretation of the Bright Excitonic Mode
The additional peak in Figs.5a and 5b stems from the bright excitonic transition p q ∥ with in-plane momenta q ∥ inside the light cone, see Eq. ( 40) and Fig. 3.It is caused by the qualitative difference between the exciton-plasmon coupling strength for bright and momentum-dark excitons that was already discussed in Sec.II D 4. With our restriction to PCs that fulfill Eq. ( 36), the only non-vanishing contribution to the third peak is p q ∥ =0 , Eq. ( 56), meaning it constitutes a weakly coupled bright excitonic mode.
The appearance of the bright excitonic mode in Fig. 5a can be traced back to a qualitative difference of our description to previous models [27,33,75,83].It results from the geometry of the considered subsystems: the PC consists of discrete dipoles that feature scattering processes.This allows the plasmon to couple to excitonic transitions outside the light cone via the plasmon-enhanced electric near-field.In contrast, the TMDC facilitates a continuous translational invariance.Therefore, it is necessary to distinguish between excitonic transitions with in-plane momenta in-and outside the light cone, providing that p q ∥ =0 does not couple to the plasmonenhanced near-field but only to its q ∥ = 0 Fourier mode.
As plasmons are popular for the amplification of the electric field in their vicinity, the near-field exciton-plasmon interaction is much stronger than the radiative (far-field) contribution.However, the q ∥ = 0 Fourier component is the propagating mode, also transferring the plasmon response to the far-field.It is physically intuitive, due to the argument of energy conservation, that this mode cannot provide a field enhancement and thus yields a reduced coupling strength to the TMDC excitons.
In Fig. 5a, we find a subtle damping and broadening due to the radiative exciton-plasmon interaction for the excitonic mode compared to the unperturbed excitonic TMDC absorption (blue), similar to the weak exciton-plasmon coupling regime.

Visibility of the Bright Excitonic Mode
Comparing Fig. 4 and Fig. 5, we find that the visibility of the weakly coupled excitonic mode strongly depends on the temperature.This trend is analyzed in Fig. 6.It illustrates the hybrid absorption spectrum for several temperatures between T = 50 K and T = 300 K with the excitonic mode clearly visible for low temperatures.With increasing temperature, the absorption peak shrinks until the increasing exciton damping reaches values of approximately half the peak splitting at room temperature.The larger the peak splitting Ω eff , that in turn characterizes the coupling strength, the smaller is the plexcitonic absorption contribution at ∆ = 0, cf.Fig. 4, which improves the visibility of the excitonic mode.
The bright excitonic mode, beyond the usually found peak splitting [28-33, 38, 79, 84], was recently observed in Ref. [40] at cryogenic temperatures.Measured was the differential reflection of a WSe 2 -monolayer, encapsulated in hBN, on top of a PC with lattice constant a = 300 nm consisting of gold nanodisks with radii r x = r y ≈ 60 nm, r z ≈ 8 nm and TMDC-PC distance δz = 12 nm.To the best of our knowledge, this was the first observation of the additional excitonic mode, since most experimental works operate either at room temperature or with the interaction strength too small.

Comparison to Previous Models
Previous theoretical analysis was restricted to a phenomenological coupled oscillator model [75] and the Jaynes-Cummings model [76] applied on TMDC-MNP and PC hybrids [28-30, 32, 38-41, 77, 81, 82, 84] or to near-field excitation involving quasinormal modes [27,83,85].Numerical approaches using Maxwell solvers have been used to confirm theoretical and experimental results at room temperature [27-32, 39-41, 86, 87].The theory developed in this contribution provides a solution for the electric near-and far-field explicitly and exceeds the limits of phenomenological models by giving a microscopic description of the TMDC excitons by including the exciton center-of-mass momentum to derive the coupling strength.
When included in the coupled oscillator model, the additional excitonic mode was described as a third oscillator and assigned to excitons spatially separated from the electric field hot spots caused by a single MNP [30] or a PC [40]; this phenomenological explanation is consistent with the developed microscopic theory, since the strong coupling is restricted to q ∥ ̸ = 0 Fourier modes, corresponding to spatially localized states [26].However, the q ∥ = 0 mode, responsible for the bright excitonic mode, constitutes a constant spatial distribution in real space.Our analysis therefore shows that the additional excitonic mode P ex q ∥ =0 is uniformly distributed also in the vicinity of the MNPs.
Another interesting feature of the absorption spectra is the asymmetry of the spectral peaks line shape with respect to width and height.The asymmetry is also covered by phenomenological models [75,76] and is consistently observed in experiments [28-31, 39, 77, 82].The reason for the asymmetry are different line shapes, line widths and dipole moments of exciton and plasmon.Within our microscopic model, another factor contributes to the observed asymmetry.The momentum-dark excitons responsible for the peak-splitting exhibit energy distributions up to 50 meV on the exciton dispersion, Eq. 18, well above the bright exciton state.Consequently, these excitons, that facilitate the strong excitonplasmon coupling, are detuned from the joint plasmon and bright exciton resonance.The magnitude of this energy difference depends on the dielectric environment, cp.Eq. 18, and the distance δz, as discussed in sec.II D 4. The detuning leads to an incomplete hybridization between plasmons and excitons.Consequently, the peak at lower energies effectively contains more plasmonic contributions which results in a stronger absorption.Figure 6.Absorption for different temperatures T between 50 and 300 K in the case of coinciding exciton and plasmon resonance.For room temperature T = 300 K (red), the hybrid features a typical strong coupling spectrum.

IV. CONCLUSIONS
We have introduced a self-consistent theory that highlights the extraordinary optical properties that can be achieved by combining ultrathin semiconductors, such as transition metal dichalcogenides (TMDCs), with plasmonic crystals (PCs) composed of metal nanoparticles.Our approach yields new analytical insights into exciton-plasmon coupling and allows to significantly reduce the numerical costs compared to general Maxwell solvers, such as ANSYS Lumerical (finitedifference time-domain) or COMSOL (finite element).
We developed a theoretical formalism to describe the coupling between a TMDC and a 2D PC, mediated by the self-consistently solved electric field.More explicitly, we considered the optical response of collective plasmons of 2D PCs within the excitonic Bloch equation.The arising effective intra-and intervalley exciton-exciton interactions for momentum-dark excitons have been decoupled in the eigenbasis of the quasi-static Green's dyadic into excitonic branches with conic respectively parabolic dispersion.Any near-field effects appeared to be incorporated via the conic contribution.In contrast, momentum-bright excitons show a qualitatively different behavior, as the coupling strength is proportional to the negligible center-of-mass momenta.Nevertheless, the momentum-dark excitons are visible in the far-field signal due to a second scattering process at the PC.
By evaluating the absorption for different TMDC-MNP separations, we were able to compare different coupling strengths and observe weak and strong coupling line shapes, the latter obeying an avoided-crossing dispersion behavior.We also find a new feature in the hybrid spectrum at low temperatures beyond the phenomenological model of two coupled oscillators: the momentum-bright TMDC excitons do not participate in the strong coupling due to their reduced coupling strength but emit undisturbed into the far-field.

Figure 2 .
Figure 2. Parabolic exciton dispersion over center-of-mass momentum q ∥ .A bright and a momentum-dark exciton are indicated exemplary.The schematic "sun" represents a radiative excitation from the far-field.

Figure 3 .
Figure 3. Schematic of the theory: the PC couples to the bright and momentum-dark TMDC excitons which are illustrated on the parabolic dispersion curve.The sun represents the far-field excitation.The gray arrows account for the exciton-plasmon coupling, with the indicated coupling strengths.

Figure 4 .
Figure 4. Absorption for different distances δz = |zpl − zex| between TMDC and PC at room temperature, corresponding to a scan over an effective coupling strength g eff .(a) Absorption spectra for δz = [11, 17, 30] nm belonging to different interaction regimes, from weak (2g eff < γex + γ pl ) to strong (2g eff > γex + γ pl ) coupling.For comparison, the artificially uncoupled case (gray), is the sum of single TMDC (blue) and PC (yellow) absorption in the inset.(b) Color map of the absorption depending on the distance δz and the Detuning ∆.The vertical lines indicate the spectra shown in panel a).Inset: Effective Rabi energy Ωeff over δz.It reaches values significantly larger than 100 meV.The peak splitting extracted from the full theory is compared to a fit b/(δz) 2 with b = 17.24 nm 2 eV.

Figure 5 .
Figure 5. (a) Absorption spectrum at liquid nitrogen temperature T = 77 K of the hybrid structure (green) compared to the uncoupled absorption of the TMDC (blue) and the 2D PC (yellow) for E ex = E pl .The effective Rabi energy Ωeff ≈ 140 meV is identified as the energetic separation between the two outer plexcitonic resonances.(b) Dispersion of the plexciton branches depending on the excitonic resonance energy E ex .The upper and lower branches are highlighted with solid lines as a guide to the eye.The hybrid spectrum shows an exciton-plasmon splitting, where plasmon (yellow) and exciton (blue) hybridize to plexcitonic modes (green).The horizontal and main diagonal mark the undisturbed plasmon (yellow) respectively exciton (blue) resonances.Inset: Splitting Ωeff (green) between plexcitonic branches depending on the difference between the resonance energies E ex − E pl .The gray linear plot indicates the splitting without any interaction.