{'fraction_non_alphanumeric': 0.03165566802394598, 'fraction_numerical': 0.010517850988601406, 'mean_word_length': 4.678335870116692, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum theory of two-dimensional materials coupled to electromagnetic resonators

We present a microscopic quantum theory of light-matter interaction in pristine sheets of two-dimensional semiconductors coupled to localized electromagnetic resonators such as optical nanocavities or plasmonic particles. The light-matter interaction breaks the translation symmetry of excitons in the two-dimensional lattice, and we find that this symmetry-breaking interaction leads to the formation of a localized exciton state, which mimics the spatial distribution of the electromagnetic field of the resonator. The localized exciton state is in turn coupled to an environment of residual exciton states. We quantify the influence of the environment and find that it is most pronounced for small lateral confinement length scales of the electromagnetic field in the resonator, and that environmental effects can be neglected if this length scale is sufficiently large. The microscopic theory provides a physically appealing derivation of the coupled oscillator models widely used to model experiments on these types of systems, in which all observable quantities are directly derived from the material parameters and the properties of the resonant electromagnetic field. As a consistency check, we show that the theory recovers the results of semiclassical electromagnetic calculations and experimental measurements of the excitonic dielectric response in the linear excitation limit. The theory, however, is not limited to linear response, and in general describes nonlinear exciton-exciton interactions in the localized exciton state, thereby providing a powerful means of investigating the nonlinear optical response of such systems.


I. INTRODUCTION
Over the last decade, there has been a growing interest in excitonic properties of two-dimensional (2D) semiconductors, especially in monolayers of the transitionmetal dichalcogenide family 1 . Owing to a direct bandgap in the visible frequency range and large exciton binding energies, these materials are particularly interesting for polaritonic physics and technology 2 . Indeed, pristine sheets of these materials have been interfaced with optical nanocavities or plasmonic resonators, leading to coupling strengths of the order of 100 meV [3][4][5][6][7][8][9][10] . Different models have been used to describe the experiments and to account for the fact that the high interaction strengths have been reached even without the need for careful positioning of the nano-resonator to align it with a local defect in the 2D material. Phenomenological treatments originating in the quantum optics literature view the electronic excitations in the 2D-material that couple to the resonant electromagnetic field as a collection of N independent dipole particles, where N is typically fitted to match the experimentally observed light-matter coupling strength 3,[6][7][8]10 . Within such N -dipole theories, the effective light-matter coupling strength is proportional to N/V , where V is an effective electromagnetic mode volume. Although the phenomenological models can be well fitted to experimental data, there appears to be no convincing microscopic theory explaining the origin or nature of the dipolar particles in the seemingly pristine 2D materials. As a consequence, the usefulness FIG. 1. Illustration of an electromagnetic resonator in the form of a plasmonic nanorod situated above an infinite sheet of 2D semiconductor material. Light-matter interaction between the resonant electromagnetic field and excitons in the semiconductor leads to the formation of a localized collective exciton mode with a center-of-mass wave function matching the in-plane electric field profile indicated by the color coding.
of the quantum optical concepts has been questioned 11 , and it has been argued that a framework rooted in the condensed-matter theory of quantum wells is more appropriate. In particular, it has been recognised that such an approach can lead to models in which the electromagnetic resonator couples to a single effective exciton state 11,12 .
In this article we build on similar ideas as put forward in Ref. 11 and develop a microscopic quantum theory for excitons in 2D materials coupled to electromagnetic res-onators (see Fig. 1). The resonator may be realized using plasmonic resonances in metal [3][4][5][6][7][8][9][10] or dielectric nanocavities [13][14][15] , including a new generation of dielectric cavities with extreme confinement of light [16][17][18] . Our analysis shows that by breaking the translation symmetry, the interaction with the electromagnetic field of the resonator leads to the formation of a localised exciton state with a center-of-mass wave function exactly matching the electromagnetic field profile. We find that this localised exciton state can be formally described as an excitonic reaction coordinate -a very successful theoretical concept developed in the context of open quantum systems [19][20][21][22][23][24][25] . With this formalism, we derive analytical results for the coupling strength between the exciton and the resonant electromagnetic field as well as the nonlinear excitonexciton interaction. These results show explicitly that the coupling strength does not scale with the effective mode volume V or the number of excitons N . Rather, the coupling strength is independent of the lateral field confinement L and depends only on the confinement in the out-of-plane direction L z . The independence of L arises from a perfect spatial overlap between the exciton reaction coordinate and the resonant electromagnetic field within the plane.
The representation of the excitons in terms of a single reaction coordinate comes at the price of introducing a residual environment of exciton states that are coupled to the reaction coordinate, but not directly to the electromagnetic field. We show, however, that these residual exciton states influence the dynamics of the system only at very small lateral confinement length scales, typically below a few nanometers. Within the same reaction coordinate formalism, one can also conveniently account for nonlinear exciton-exciton interactions, and we find that the lateral field confinement here plays a crucial role: It determines the effective area of the exciton reaction coordinate, and the nonlinear interaction of excitons within the reaction coordinate therefore scales as 1/L 2 , reflecting the fact that the optical confinement dictates the multi-excitonic co-localisation.
We present three different approaches for calculating the time evolution of the system and use these methods to assess the influence of the residual exciton environment and the nonlinear response. In this way, we identify a range of lateral confinement lengths where L is sufficiently large that the residual excitons can be ignored but at the same time small enough that the nonlinear response significantly alters the dynamics.
The article is organised as follows. In Sec. II, we describe the Wannier-Mott exciton states of the 2D materials, the electromagnetic fields of the resonators, and their interaction. In Sec. III, we derive the exciton reaction coordinate formulation for coupling to a single mode of the electromagnetic resonator. In Sec. IV, we calculate the time evolution of the system using three different approaches, which we benchmark against each other to assess their regimes of validity. In Sec. V, we derive the effective linear dielectric function of the Wannier-Mott excitons and use this function to make a reference calculation of the excitation spectrum, which we compare to the microscopic quantum reaction coordinate approach in the linear-response limit. Finally, we summarise our findings in Sec. VI.

II. GENERAL FRAMEWORK
In this section, we present the fundamental structure of the theory, which is based on Wannier-Mott exciton states and their interactions with a resonant electromagnetic field described by a single quasi-normal mode (QNM). We will generally study excitons in directbandgap semiconductors with discrete in-plane translational symmetry.
In this sense, monolayer twodimensional semiconductors such as transition metal dichalcogenides share many physical features with semiconductor quantum wells, although the excitons of the former are often more strongly bound due to their reduced dielectric screening 1 . Monolayer transition-metal dichalcogenides have direct bandgaps at the K and K points 26,27 and feature a rich electronic band structure (cf. Fig. 2a,c). In the vicinity of these K and K points, however, the conduction and valence bands can be well approximated by parabolic bands, leading to an effective-mass approximation, which we shall use here. The optically bright excitons generated from these bands are Coulomb-bound electron-hole states. Comprehensive theoretical treatments of excitons resolve their composite fermionic electron-hole structure [28][29][30][31][32][33][34][35] . Here, we shall employ a simpler description of the excitons in terms of interacting bosons [36][37][38][39][40] .
The resonant electromagnetic fields in optical cavities and plasmonic particles share many characteristics with the bound states of electrons. It is a distinct feature of electromagnetic resonators, however, that the modes are not truly bound, and this gives rise to discrete peaks with finite widths in scattering spectra, for example. From a mathematical point of view, it is advantageous to treat these resonances as modes of the electromagnetic field with finite lifetimes, and the theory of QNMs provides a rigorous framework for doing this [41][42][43][44] . In this work, we will start from a quantum description of the electromagnetic field in electromagnetic resonators 45 and extend the theory to describe the interaction with Wannier-Mott excitons in 2D materials. At positions far from the resonator, a QNM description of the electromagnetic field is non-trivial, and this poses a challenge for coupling to excitons in the nominally infinite sheet of 2D material. Such a description is in principle influenced by retardation effects, but for the present purpose of describing interactions very close to the resonator we avoid these complications by treating the interaction in the quasistatic limit. . The band energies are plotted as function of in-plane wave vector and are given relative to the Fermi energy, EF. Within the effectivemass approximation, the conduction and valence band energies near the band gap, (K or K ), are approximated as parabolic (blue and green dashed lines). b. Conduction band (blue) and valence band (green) band electronic energies relative to valence band maximum, E vbm (left axis) as a function of the wavevector, k, given relative to K (or K ). Also shown is the lowest-lying 1s exciton energy (right axis) as a function of the center-of-mass momentum, k. The band gap, Eg and the exciton gap, E0 = Eg − E b , are indicated with vertical arrows. The effective masses of the valence and conduction bands are independent of the direction in k-space 46 and the exciton dispersion is thus direction-independent as well. c. The Brillouin zone of a 2D hexagonal lattice with special symmetry points marked as red dots and the path over which the band structure in panel a is calculated (grey line).

A. Exciton states
The description of the excitonic degrees of freedom in the 2D material will be based on a Wannier-Mott framework 47 , which provides a useful analytical description of the excitons and has proven successful in this context 48,49 . To this end, we consider the lowest-energy conduction band and the highest-energy valence band of a two-dimensional semiconductor, separated by a band gap energy E g . We denote by |0 the Fermi sea, i.e. the fermionic state in which the valence band is fully occupied by electrons and the conduction band is empty. The fermionic creation operator for a hole in the valence band with in-plane wavevector k is denoted by v † αk and, simi-larly, c † αk denotes the creation operator for an electron in the conduction band. Here, the index α labels the highsymmetry point with wavevector K α , when the band gap is degenerate, as is the case for e.g. monolayer transitionmetal dichalcogenides, where α ∈ {K, K } 1 . The singleparticle wave functions, ψ c,αk (r) and ψ v,αk (r) are taken to be of Bloch form, where i = c, v; N is the number of unit cells in the 2D sheet with surface area S, and u i,α (r) is a Bloch function, which has the periodicity of the crystal lattice and is normalised over a single unit cell, VUC d 3 r |u i,α (r)| 2 = 1. Within the Wannier-Mott framework, the conduction band and valence band energies are described through effective masses, m e and m h , respectively, which reflect the local parabolic approximation to the band structure near the K and K valleys, cf. Fig. 2a.
The Wannier-Mott exciton states can be written as a momentum-superposition of electron-hole pairs 50 where φ n (q) is the momentum-space exciton wave function with shell index n, and M = m e + m h is the total exciton mass. The label k is thus the center-of-mass momentum of the exciton. In Ref. 49, it was found that a very good approximation to the lowest-energy exciton (n = 1) for several two-dimensional transition-metal dichalcogenides is given by the simple hydrogenic form in which a B is the exciton Bohr radius and q = |q| . This form is often found in the context of semiconductor quantum wells 51,52 . In the present work, we restrict our discussion to the lowest-lying exciton and generally use Eq. (3) for the exciton wave function and denote its quantum state by |Φ α,k , dropping the index n in Eq. (2). The total exciton energy is given by the sum of the kinetic energy of the center-of-mass coordinate, the band gap E g and the exciton binding energy, E b : where E 0 := E g − E b is the exciton gap, as illustrated in Fig. 2b. We will often use the corresponding frequencies as well, ω k = E k /h, and ω 0 = E 0 /h. Being composed of electron-hole pairs, it is favourable to describe the excitons through a set of bosonic creation operators,b † αk , which generate single-and multi-excitonic quantum states. Noting, however, that operators of the form have commutation relations that are neither bosonic nor fermionic 53 , a description of the excitons as noninteracting bosons is infeasible. Nevertheless, it is possible to use an approximate description in terms of interacting bosons [36][37][38][39][40] . In this work, we include the leading interaction term, which arises due to exciton-exciton Coulomb interactions 51,52 , such that the total bosonic form of the exciton Hamiltonian (including the free exciton energy) iŝ wherehW kk q is the momentum-dependent excitonexciton interaction energy andb αk andb † αk are bosonic exciton annihilation and creation operators, obeying the commutation relation [b αk ,b † α ,k ] = δ kk δ αα . These operators can be understood as a bosonic approximation to the excitonic operatorsĈ αk andĈ † αk , such thatb † αk |0 is the bosonic representation of |Φ αk . The state |0 is the bosonic exciton vacuum state, which is the bosonic equivalent to the Fermi sea 36 . In the present context, the correspondence between excitonic fermion pairs and bosons is only of formal interest, since the practical calculations of the microscopic light-matter coupling strength can be phrased in terms of the fermionic operators. Similarly, the interaction strengths W kk q are directly related to four-particle scattering matrix elements evaluated in the original fermionic space 36,38,51,52 . For convenience, we can split the exciton Hamiltonian into a noninteracting part,Ĥ x,0 = αkh ω kb † αkb αk , and an interacting part,Ŵ = αkk qh W kk qb † αk+qb † αk −qb αk b αk . Detailed discussions and calculations of the interaction strengths W kk q can be found in Refs. 51,52,54,55 . In general, a momentum cutoff of the order a −1 B is reported, which means thathW kk q hW 000 for k, k , q a −1 B . In Ref. 55, variational calculations of the exciton wavefunctions and the corresponding Coulomb matrix elements showed that for monolayer WS 2 , the interaction strength can be approximated ashW 000 2.07E b a 2 B /S. It was also shown that the interaction has contributions from a direct part, where the constituent electron and hole within the two interacting excitons remain fixed, and an exchange part, where the constituent particles are exchanged. In particular, it was shown that the exchange part dominates at low momenta, and that the direct part vanishes identically at zero momentum, meaning that the matrix element W 000 is determined solely by the exchange contribution. We note, furthermore, that intervalley exchange Coulomb effects in transition-metal dichalcogenides 56 are not included in this work. Besides from an overall shift, the linear exchange coupling vanishes for vanishing exciton center-of-mass momentum [56][57][58] , and is therefore small for optically bright excitons, which have small center-of-mass momenta. The nonlinear intervalley exchange Coulomb interaction is small compared to the intravalley Coulomb interaction,Ŵ 34 . Thus, we expect that the intervalley exchange interaction will give rise to minor corrections to the overall physics of light-matter coupling and the weak nonlinear optical response studied here.
A collection of relevant material parameters for a selection of monolayer transition-metal dichalcogenides is presented in Table I.

B. Resonant electromagnetic fields
Optical cavities and plasmonic particles support a number of resonances, which show up as distinct peaks in scattering spectra, for example, and whose corresponding field distributions are commonly referred to as quasinormal modes (QNMs) [41][42][43][44] or resonant states 61,62 . They are defined as solutions to the wave equation subject to suitable radiation conditions, such as the Silver-Müller condition for resonators in free space 63,64 . For the present application, we limit the analysis to cases in which there is only a single QNM in the frequency range of interest, and we denote the associated electric-field distribution byf c (r); the corresponding eigenfrequencyω c = ω c − iγ c is complex with a negative imaginary part, where γ c accounts for the cavity decay rate. For positions close to or inside the resonator, the operator describing the quantized electric field can be expanded in terms of this QNM by the method presented in Ref. 45. In this work, we shall be interested in coupling to sheets of 2D materials, which extend to regions far away from the resonator, and where the formulation of Ref. 45 is not directly applicable. As discussed in Appendix A, however, it is possible to extend the general framework of Ref. 45 by use of the Lippmann-Schwinger equation and ideas originally put forward in Ref. 65 to write the electric-field operator at general positions r outside the resonator in terms of a convolution aŝ Here, the creation and annihilation operatorsâ † c andâ c obey the commutation relation [â c (t),â † c (t)] = 1 45 , and the memory kernelF c (r, t), which ensures a proper causal relation to locations far from the resonator, derives from analytical continuation of the electric field QNM onto the real frequency axis. In the present work, we shall focus on the local and non-retarded coupling dynamics by setting F c (r, t) =F c (r)δ(t) and calculatingF c (r) fromf c (r) by use of the quasistatic Green tensor, see Appendix A for details. Moreover, as we shall see below, the coupling will be phrased in terms of the electromagnetic vector potential, which we write in the single-QNM case aŝ Due to radiative loss and possibly absorption in the resonator material, the energy of the electromagnetic field  in the resonator is not conserved. Therefore, it is convenient to describe the state of the field through its density operator,ρ c , which in the absence of interactions is governed by the master equation 45 where is the Lindblad dissipator.

C. Light-matter interaction
We now proceed to consider interactions between the excitons and the resonant electromagnetic field. The exciton-field coupling is generated by the minimalcoupling Hamiltonian, of which the dominating part (written in first quantization) is 66 where the sum runs over the electrons in the system. The light-matter interaction Hamiltonian can be cast into the second-quantised form 67 where the interaction strengths are evaluated in the fermionic space as andhg αk is obtained by replacingF c (r) withF * c (r). Assuming that the light-matter coupling is weak compared to the relevant electromagnetic frequencies and exciton frequencies, we make the rotating-wave approximation, corresponding to setting g αk = 0 in Eq. (11). With these approximations, the light-matter interaction iŝ Using the Slater-Condon rules 47 , the matrix element entering g αk can be evaluated, again in the fermionic space, as where we used the expression for the single-particle wave functions of Eq. (1). The summation over q only involves the wave function φ(q) and can thus be evaluated independently as Assuming that the mode function does not vary appreciably over the unit cell, and that the relevant wavevectors are much smaller than the inverse lattice constant, we can approximate the integral in Eq. (14) as a sum over unit cells (indexed by j) and a Bloch matrix element where p α cv = VUC d 3 r u * c,α (r)pu v,α (r), r j is the lateral coordinate of the jth unit cell, and we have taken the 2D-sheet to be located at z = z 0 . For 2D materials in the transition-metal dichalcogenide family, the two degenerate exciton modes at the K and K valleys have the matrix elements 27,68 , p K cv = m 0 V (x + iy), p K cv = m 0 V (x − iy), where V is a material-dependent velocity parameter (see Table I for values). These matrix elements are circularly polarized due to spin-orbit coupling. The summation over j in Eq. (16) can then be rewritten as an integral as j → N S d 2 r. Combining everything, the coupling strength becomes where the integral is over the infinite extent of the 2D material. Whereas the QNMs diverge (exponentially) at sufficiently large distances from the resonator 44 , the func-tionsF c (r) behave as the free-space quasistatic Green tensor at sufficiently large distances and are therefore square integrable. We note that in addition to the contributions to the light-matter interactions derived until this point, a nonlinear saturation term of the form appears due to the non-bosonic nature of the excitons 51,52 . Within the hydrogenic Wannier-Mott exciton approximation, the nonlinear interaction due to saturation can be approximated in the zero-momentum limit, where it takes its maximal value, ashσ α000 (4π/7)(a 2 B /S)2G 0 , where G 0 is the coupling strength 51 . For a coupling strength of 50 meV and an exciton binding energy of 500 meV, this gives a ratio relative to the exciton-exciton Coulomb interaction of σ α000 /W 000 0.17. Thus, for the systems of interest in this article, the Coulomb exciton-exciton interactionŴ is considered to be the strongest nonlinear effect, and the saturation interaction is neglected. We do note, however, that such saturation effects have been shown to play an important role for the nonlinear dynamics of trions in transitionmetal dichalcogenides, owing to the three-fermion composite structure of trions 69,70 .

III. EXCITON REACTION COORDINATE
In Section II, we have seen that the excitons are described by a HamiltonianĤ x , which conserves the centerof-mass momentum k, thus reflecting the discrete translation symmetry of the 2D material sheet. This symmetry is then explicitly broken through the interaction Hamil-tonianĤ I , because the electric-field distributionF c is not translationally invariant. In this section, we introduce a basis transformation of the excitons that provides a natural starting point for analyzing this system. The transformation defines a localized exciton mode, which we denote the exciton reaction coordinate in the spirit of quantum chemistry, where the concept of reaction coordinates has been developed in a similar fashion to describe nuclear motion in molecules [19][20][21][22][23][24] . In the new basis, the light-matter interaction is greatly simplified, because the electromagnetic field couples only to the exciton reaction coordinate. The transformation comes at the cost of introducing a reservoir of residual exciton modes, which in turn are coupled to the reaction coordinate. However, as we shall see in Section IV, there are several successful approximate strategies for treating these residual modes.
A. Exciton reaction coordinate and residual excitonic spectral density We now define a new basis of exciton modes, {B i }, generated by a unitary transformation U as 25 in which the first row in the transformation matrix U is given by We denote this collective mode as the exciton reaction coordinate. The remaining rows in U are constructed via Gram-Schmidt orthogonalisation as orthonormal vectors, The light-matter interaction can now be written in the much simpler form where G 0 = αk |g αk | 2 . This illustrates the advantage of introducing the reaction coordinate: by construction, the resonant electromagnetic field now only interacts with the weighted sum of Wannier-Mott excitons defined byB 0 . The transformation U thus precisely captures the notion of a localized excitonic state, which is defined by the interaction with the electromagnetic field. This physically appealing reformulation of the dynamics comes at the computational price that the associated free-exciton HamiltonianĤ x,0 is no longer diagonal, where Ω ii = αk E k U iαk U * i αk and thus Ω ii = Ω * i i . In this way, the exciton reaction coordinate is coupled to a residual environment of other exciton modes,B i , i > 0. As we shall see later, it is useful to have access to the spectral density of this environment, and this requires a re-diagonalisation as follows. The procedure starts by separating out the terms inĤ x,0 that contain the reaction coordinateB 0 , This can be written in the compact matrix form where Ω is constructed from Ω by removing the first row and column, andB and λ are vectors with elementsB i and λ i = Ω i0 , respectively, with i > 0. We now wish to re-diagonalise the part of the Hamiltonian governing the modes with i > 0. To do this, we defineŨ as the transformation that diagonalises Ω , such thatΩ :=Ũ † Ω Ũ is diagonal. We also define a transformed set of modes, B =ŨB. Since Ω is a Hermitian matrix,Ũ is unitary, and the columns ofŨ are the eigenvectors of Ω . In this way, the last term in Eq. (25) can be written in a diagonal form, and the free-exciton Hamiltonian can be written aŝ whereλ =Ũ † λ and we have defined Ω 0 := Ω 00 ,Ω i := Ω ii . From Eq. (26) it is evident that the exciton states of the residual environment are mutually uncoupled and couple only to the exciton reaction coordinate. The spectral density of the residual environment can be calculated from Eq. (26) as This is in short called the residual spectral density of the exciton reaction coordinate. Although the procedure of constructing the reaction coordinate transformation matrix U and subsequently the re-diagonalisation matrixŨ can be carried out numerically, this scales poorly with the number of elements in the matrix. Fortunately, the procedure of extracting reaction coordinates and residual spectral densities has been extensively studied in the literature of open quantum systems 24,25,[71][72][73] , and there are direct ways of obtaining the residual spectral density without going through the intermediate steps as above.
The starting point for such analyses is the exciton spectral density of the full set of exciton states: from which we can derive the necessary quantities related to the exciton reaction coordinate. Starting with the reaction coordinate frequency, we have which is simply the first moment of J(ω). The coupling strength can be calculated similarly, The residual spectral density is more involved. Here we shall not prove the relation between J and J res , but state a result from Ref. 71, where is the so-called reducer and the interval [a, b] is the frequency support of J(ω) 71 ; in the present case, the support of the spectral density is [ω 0 , ∞). When calculating the reducer numerically, it is convenient to use the form 71 where the upper limit of the support, b, is chosen sufficiently large that J res can be considered independent thereof. To deepen the intuitive understanding of the reaction coordinate transformation for our 2D excitonic systems, and for practical reasons, we can rewrite the creation operator of the exciton reaction coordinate by introducing the real-space exciton operators where ψ α 0 (r) is the real-space exciton reaction coordinate wave function, From the form of ψ α 0 (r), it is clear that the exciton reaction coordinate inherits its spatial distribution from the vectorial projection of the electromagnetic field distribu-tionF c (r) onto the plane of the 2D material.
It is useful to establish the connection between the coupling strength G 0 , and the electromagnetic field distribution given byF c (r). This can be obtained by using the general exciton-field coupling coefficient, Eq. (17), along with the expression for G 0 , Eq. (30), which leads directly tō Importantly, the coupling strength is independent of the lateral confinement length of the optical mode. This can already be seen in Eq. (36), where the integral over the mode extends over the entire 2D surface. The independence of G 0 on the lateral confinement stems from the fact that the exciton reaction coordinate wave function, Eq. (35), is perfectly matched with the optical mode within the 2D plane. In Sec. III C, we show this analytically for the case of a mode profile that is separable in the in-plane and out-of-plane coordinates and see that the coupling strength only depends on the out-of-plane confinement length. In Appendix D, we have also performed numerical calculations with a QNM of a gold nanorod resonator as in Fig. 1, where the resonator length is varied. The calculations show that the contribution to G 0 from the integral over the normalised QNM profile in Eq. (36) varies only about 1% when the length of the nanorodand hence the lateral confinement length of the QNMis varied between 80 nm and 95 nm.

B. Exciton-exciton interactions within the reaction coordinate
The exciton-exciton interaction,Ŵ , can be rewritten in the transformed basis of exciton modes In the dynamical model derived below, we account only for the exciton-exciton interactions within the exciton reaction coordinate. This is justified by the assumption that this is the only region of the exciton Hilbert space where the exciton density is sufficiently large to give a significant contribution to the dynamics. Thus, keeping only the term i, i , j, j = 0 in the summation, we end up with the interaction term with the effective nonlinear interaction strength When the characteristic confinement length of the electromagnetic field is large compared with the exciton Bohr radius, the transformation elements U 0k decay on a momentum scale that is small compared to the momentum variation of W kk q . In this regime, we can approximate W kk q W 000 . Furthermore, using the spatial reaction coordinate wave function, ψ 0 , the effective nonlinear interaction strength can be vastly simplified as This result is consistent with Ref. 12, where it is assumed that there exists an excitonic eigenmode with the same spatial wave function as the resonator mode. As we shall see explicitly later, the integral of |ψ α 0 (r)| 4 is a measure of only the lateral confinement of the optical mode. Thus, the lateral optical confinement is inherited by the exciton reaction coordinate, thereby determining the interaction strength between excitons within the reaction coordinate.

C. Localized and separable mode profiles
In this subsection, we investigate the coupling dynamics in an idealized limit of an electromagnetic mode function that is localized and separable. We do this in order to illustrate the mechanism by which the perfect colocalization of the exciton reaction coordinate and the resonant field leads to a coupling strength that is independent of the lateral extent of the electromagnetic field.
The electromagnetic field distributionF c (r) is localized in space, as discussed in Section II B. It is advantageous, therefore, to think ofF c (r) as a hypothetical localized solution to the wave equation with a purely real frequency. Even if all resonant electromagnetic modes of optical cavities and plasmonic particles have finite Qvalues and are leaky in nature, the abstraction of perfect temporal and spatial confinement allows us to derive a number of interesting analytical results regarding the effect of the electromagnetic field distribution on the coupling strength. At positions close to the resonator,F c (r) closely resembles the QNMf c (r) since they are related by the analytical continuationω c → ω c . Therefore, we expect the general findings to apply qualitatively also to electromagnetic resonators with finite Q-values.
With this motivation, we consider now a localized mode functionF c (r) in an environment of constant and real permittivity , and we assume that the field is separable in the lateral and perpendicular coordinates, such thatF where n is the unit polarisation vector of the mode. Sincẽ F(r) is localized in space, it obeys the normalization requirement in which eff is an effective dielectric constant to account for the dielectric response of the surrounding material. While (r) can be considered constant within the 2D semiconductor sheet, it generally varies across the surrounding structure. In order to focus on the influence of the confinement length scales, however, we approximate the combined effect by use of an effective dielectric constant eff . We note that this assumption is not generally a requirement for using the theory, and for the numerically calculated QNM that we describe in Sec. V and App. D, we do not make any such simplifying assumptions about the dielectric environment.
Using the separable mode function in Eq. (41) and the normalization reguirement in Eq. (42), we can demonstrate that G 0 is entirely independent of the lateral mode distribution by rewriting it as where is the out-of-plane confinement length. It follows from Eq. (43) that the coupling strength is independent of the lateral mode distribution, and that it scales with the outof-plane confinement length as G 0 ∝ L −1/2 z . For ease of notation, we shall use an implicit summation over the valley index as |n · p cv | = ( α |n · p α cv | 2 ) 1/2 when stating the polarisation overlap used in specific calculations.
The reaction coordinate exciton-exciton interaction strength W 0 can be written using Eq. (40) as where is a polarisation-dependent prefactor. This prefactor takes values between 1/2 and 1. The maximal value, 1, is obtained when n is orthogonal to one of the momentum matrix elements, i.e. n·p α cv = 0. Conversely, the minimal value, 1/2, is obtained when |n · p α cv | is equal for the two polarisations.

Gaussian lateral field distribution
We now proceed by assuming that the lateral field distribution defined by the electromagnetic resonator is Gaussian with a confinement length scale L, We note that d 2 r |F (r)| 2 = 1, so that the normalisation requirement is eff dz |F z (z)| 2 = 1, and the exciton coupling strengths in Eq. (17) can be written as This relatively simple expression for the coupling strength allows us to evaluate the expression for the spectral density in Eq. (28) by writing the summation over k as an integral, k → S (2π) 2 d 2 k. In this way, we find that we can write the spectral density compactly in terms of the reaction coordinate coupling strength in Eq. (43) as where ξ =h/(2M L 2 ) is a cutoff frequency and Θ is the Heaviside function. Using this analytical expression for the spectral density, it is also possible to evaluate the residual spectral density in Eq. (31) analytically, where Ei(x) = x −∞ dz exp(z)/z is the exponential integral function. Fig. 3 shows an example of the spectral density and residual spectral density generated by the Gaussian field distribution in Eq. (47) coupled to WS 2 and MoS 2 , respectively. As can be seen from Eq. (49), the magnitude of the spectral density scales as the ratio G 2 0 /ξ. Here, G 0 depends on the material-specific parameters p α cv and a B along with the out-of-plane confinement length scale L z .
The cutoff frequency ξ depends on the material-specific total exciton mass M as well as the in-plane confinement length scale of the resonator, L. The magnitude of the residual spectral density, on the other hand, scales as the cutoff frequency ξ. Thus, the cutoff frequency determines not only the relevant frequency scale of the exciton spectral density, but also the relative strength of the interactions with the residual exciton environment.
The reaction coordinate exciton-exciton interaction strength W 0 can be obtained from Eq. (45) as As anticipated in Sec. I, this nonlinear interaction strength scales inversely with the confinement area, L 2 , reflecting that the co-localisation of multiple excitons (and thereby their interaction) is fully determined by the electromagnetic field distribution. For the calculations in this article that involve nonlinear interactions with Gaussian electromagnetic modes (Sec. IV C 3), we take the polarisation prefactor to be unity, corresponding to a circularly polarised optical mode.

IV. TIME EVOLUTION
Based on the general framework in Sec. II and the reformulation in terms of the reaction coordinate in Sec. III, we are now in a position to calculate the time evolution of the excitations in the coupled system comprising the electromagnetic resonator and the 2D material.
In this section, we present three strategies of increasing complexity for calculating the system dynamics. The first method, which is exact when excitonic broadening effects are ignored, is based on direct time evolution of the single-excitation product states formed by the Fock states of a single electromagnetic excitation and the continuum of single-exciton states with given momentum. The second method is based on a Markovian treatment of the residual exciton modes, and a master equation for the reduced density operator of the resonant electromagnetic field and the exciton reaction coordinate is derived. In this formulation, one can account for excitonic decay and dephasing within the reaction coordinate, and external driving of the system can be included as well. The Markovian master equation is benchmarked against the exact treatment when excitonic broadening effects and external driving are ignored, thereby providing a reference calculation to assess the Markov approximation for the interaction with the residual exciton modes. The third method is based on an iterative extension of the reaction coordinate mapping, allowing one to represent the residual exciton modes by a one-dimensional chain of bosonic modes. This, in turn, enables the model to account for non-Markovian features of the residual environment, although our implementation can only evolve the system up to a finite time. Furthermore, using this approach, it is also straightforward to include external driving with a time-dependent amplitude, as is the case when the system is excited by a laser pulse.
A. Exact evolution in the single-excitation sector As a relatively simple starting point, and to establish a benchmark for evaluating the precision of more complicated approaches for time-evolution calculations, we consider in this section the Fock state representation of a simplified model system. This approach is motivated by the fact that the dynamics can be solved exactly, when the coupled system is restricted to the one-excitation sector. In practice, we do this by initialising the resonant electromagnetic field in a single-photon Fock state, such that nonlinearities can be neglected. We can then expand the combined state as where |1; {0} denotes the state with a single electromagnetic excitation and zero excitons and |0; 1 αk denotes the state with no electromagnetic excitations combined with a single exciton with momentum k. The amplitude φ c then obeys the equation of motion 74 where Compared to Ref. 74, the roles of the electromagnetic and electronic degrees of freedom are interchanged, so that the excitons act as a continuum with which the single electromagnetic field of the resonator interacts. Since J(ω) does not suffer from an ultraviolet divergence, there is no need to introduce a cut-off in order to evaluate the integral. Although this method is in principle exact, it is limited to single-excitation problems and cannot be used to model problems with external driving or dephasing effects. For this reason, we use the method mostly for reference calculations for the master equation formulations to be described below.

B. Secular Markovian master equation
In this section we derive a secular Markovian master equation for the reduced density operatorρ of the system comprising the resonant electromagnetic field and the exciton reaction coordinate, where the residual exciton modes are traced out. We subsequently benchmark this master equation approach against the exact method from section IV A. We do this in the limit where dephasing and driving effects are turned off in order to assess the accuracy of the master equation and study the conditions for describing the interactions with the residual exciton modes within the Markov approximation.
In the Hamiltonian governing the evolution of the resonant electromagnetic field, we include external laser driving with frequency ω d and constant amplitude F in addition to the free evolution described in Section II B, such that the part of the Hamiltonian governing the field evolution isĤ c =hω câ † câc +hF (e iω d tâ c + e −iω d tâ † c ). Next, we divide the total Hamiltonian into contributions describing the exciton reaction coordinate and the resonant mode, being the system (S), the residual exciton modes, being the reservoir (R), and their mutual coupling (SR), asĤ =Ĥ S +Ĥ R +Ĥ SR . Putting it all together, we havê in which we have expressed the Hamiltonian in a reference frame rotating with the driving frequency ω d . The standard Markovian master equation obtained by tracing out the environment is given by 75 in which the Lindblad dissipator D(x,ρ) is specified in Eq. (9), and whereĤ SR (−τ ) = e −i(ĤS+ĤR)τ /hĤ SR e +i(ĤS+ĤR)τ /h is the interaction-picture time evolution ofĤ SR andρ 0 R is the initial state of the residual excitonic environment, which we take to be the exciton vacuum. As discussed in Appendix B, we can use the secular approximation to write the exciton dissipator K[ρ] in a simpler form as where the sum is over all system eigenfrequency differences ω, and the exciton dissipation rate Γ res (ω) can be written in terms of the residual spectral density as Γ res (ω) = 2πJ res (ω + ω d ).
The operatorsB 0 (ω) andB † 0 (ω) are eigenstate-projected exciton operators, which are described in detail in App. B, where the term ∆ res (ω) is also defined. In the weak-driving limit, F G 0 , γ c , only the low-energy states of the system are populated, and G 0 dominates the structure of the coupled system. In this limit, we can calculate the eigenstates ofĤ S in the basis {|0, 0 , |1, 0 , |0, 1 }, where |n c , n x denotes a Fock state with n c electromagnetic energy quanta in the resonant field mode and n x excitons in the reaction coordinate mode. We find that there are only two non-zero contributions to the summation over ω, corresponding to the upper and lower polariton modes, with frequencies The corresponding eigenstate-projected operators are the annihilation operators of the upper and lower polaritons, where δ cx = ω c − Ω 0 is the exciton-resonator detuning and η = 4G 2 0 + δ 2 cx is the polariton splitting. The bars overω ± signify that the polariton frequencies are given in the rotating frame. They are related to the corresponding lab frame frequencies ω ± as ω ± =ω ± + ω d . In the case where the lower-polariton energy ω − is below the exciton gap ω 0 , (at resonance, this amounts to G 0 > (Ω 0 − ω 0 )) the contribution to the dissipator fromB 0 (ω − ) vanishes, and we can write the master equation in the simplified whereB 0 :=B 0 (ω + ). These derivations put us in a position to understand the impact of the different energy scales in the system on a deeper level. In Fig. 4, the residual spectral density is shown along with an indication of the reaction coordinate frequency Ω 0 . The upper and lower polariton frequencies ω ± are also indicated for the resonant case ω c = Ω 0 , where ω ± = Ω 0 ± G 0 . From Eq. (59), we know that the dissipation rates from the polaritons into the residual environment are given by 2πJ res (ω ± ). Thus, when the coupling strength is increased far beyond the cutoff frequency ξ, the polaritonic peaks are pushed away from the peak of J res (ω), and the effective interaction strength with the residual excitons is consequently reduced. Correspondingly, if the cutoff frequency is decreased, for example by increasing the lateral confinement length L, the same reduction in effective interaction strength is observed. In the following section, we shall see that this phenomenon of effectively decoupling the residual excitons when G 0 ξ leads to a large parameter regime where the residual environment can be safely ignored.

Benchmarking
In Fig. 5a we show the temporal evolution of the resonator population for three different lateral confinement length scales. The initial state is a single excitation in the resonator. We compare the time evolution calculated by the exact equation of motion of the resonant electromagnetic field in Eq. (53) (black solid lines) with the secular Markovian master equation with all terms included in Eq. (58) (green dots), as well as the simplified Markovian master equation in Eq. (62) (dashed orange lines). For reference, we have also shown the time evolution generated when the residual excitons are entirely ignored (red dotted lines). In general, the dynamics show Rabi oscillations due to the interaction between the exciton reaction coordinate and the resonator field. These oscillations are damped due to Markovian losses of the resonator and interactions with the residual excitons. As discussed earlier, the effect of the residual excitons becomes more pronounced as the lateral size decreases. In Fig. 5b, the relative errors of the three approximate approaches are shown. The errors are calculated as where â † c (t)â c (t) is the population of the electromagnetic resonator obtained with one of the master equations and φ c (t) is obtained from the exact time evolution. Furthermore, the Markovian decay rate into the residual exciton modes Γ res is shown (blue solid line, right y-axis). The dependence of the errors on the lateral length scales confirms that the influence of the residual excitons is more pronounced at small length scales. For the chosen parameters, we find that the residual excitons can be ignored when L > ∼ 4 nm. Furthermore, we see that the Markovian master equation provides a useful description of the interactions with the residual environment, which improves the accuracy of the calculated time evolution. We also conclude that the error of the simplified Markovian master equation, Eq. (62), is comparable with the full Markovian dissipator, Eq. (58) with all secular terms included. In a significant portion of the parameter regime, the simplified master equation even performs slightly better.

Inclusion of excitonic broadening effects
Because of interactions with lattice phonons, the excitons experience population decay and dephasing, which result in a temperature-dependent broadening of the exciton line 76,77 . In Ref. 76, it was found that the dominant mechanisms generating excitonic decay are radiative recombination (which is temperature-independent) and scattering with phonons at the Λ point. In addition, an important contribution to the exciton linewidth was found from intra-valley scattering with phonons at the Γ point. Specifically, the total exciton linewidth of various monolayer transition-metal dichalcogenides was found to be well described by the temperature-dependent expression Γ x = γ 0 + c 1 T + c 2 /[exp{Ω/k B T } − 1], where c 1 , c 2 and Ω are material-dependent coefficients. For WS 2 and WSe 2 , the temperature independent term γ 0 contains contributions from radiative decay and spontaneous emission of Λ-phonons. The linear coefficient, c 1 , stems from intra-valley scattering with thermally excited phonons at the Γ point, and the last term accounts for interactions with thermally excited phonons at the Λ point. The first and last terms describe processes that lead to a decay of excitons. Whereas radiative decay is accounted for through the interaction with the electromagnetic field, H I , we include the possibility of scattering into dark exciton states as a decay term in the master equation of the form 2γ x D(B 0 ,ρ), where γ x = γ nr 0 +c 2 /[exp{Ω/k B T }−1] and γ nr 0 is the non-radiative contribution to γ 0 . Intravalley scattering with Γ-phonons, on the other hand, does not lead to a population decay but rather to a dephasing, similar to virtual phonon transitions to higherlying excited states seen in systems with localised exciton states [78][79][80] . This process can be included in the master equation by a dephasing term, 2γ x D(B † 0B 0 ,ρ), where γ x = c 1 T . Thus, although the total exciton linewidth, Γ x = γ x + γ x , depends only on the sum of the two contri-  76 for calculation of temperature-dependent exciton decay and dephasing rates. The parameters are the temperature-independent contribution to the non-radiative decay from phonon interactions, γ nr 0 , the contribution to dephasing from intravalley phonon scattering, c1, the contribution to decay (and for MoS2 also dephasing) due to intervalley phonon scattering, c2, and the typical phonon energy for intervalley processes, Ω. butions, it is important to note that there is a difference in nature between population decay and dephasing processes. In Ref. 76, monolayer MoS 2 was also studied, and it was found that the coefficient c 2 contains contributions both from intra-valley scattering, i.e. dephasing, and from interactions with Λ phonons, i.e. decay into dark exciton states. The parameters from Ref. 76 are presented in Table II.
We note that recent studies indicate strain as a possible way of energetically shifting the direct K-valley exciton below the indirect, momentum-dark K − Λ exciton in monolayer WSe 2 81 . Such an energetic cross-over would arguably lead to significant reduction of phononinduced decay of the bright exciton, even with a vanishing decay at zero temperature, as is the case for Mo-based transition-metal dichalcogenides, where the K − Λ exciton is below the direct K-exciton in the absence of strain.

C. Non-Markovian treatment of residual excitons using chain mapping
In some situations, e.g. for small lateral confinement scales, it may be necessary to account for non-Markovian effects in the interaction with the residual excitonic environment. This can be done by extracting additional reaction coordinates from the residual environment, thereby extending the system Hilbert space to include the most important environmental degrees of freedom. Since the exciton reaction coordinate Hamiltonian in Eq. (26) is structurally equivalent to the original Hamiltonian, we can iterate the procedure of extracting reaction coordinates, thereby generating a one-dimensional chain of coupled modes, as illustrated in Fig. 6. This strategy has been formally studied in the literature 71,73 ; here, we derive it iteratively, starting from the reaction coordinate transformation. First, we re-write Eq. (26) by adding the superscript '(0)' to the residual modes and quantities related to them, The superscript '(0)' indicates that the residual environment is coupled to the reaction coordinateB 0 . Correspondingly, we label the spectral density of theB We can now define a new reaction coordinatê By repeating the procedure of defining a new residual reservoir and re-diagonalising it, as in Sec. III, we obtain where res (ω). Thus, from the residual exciton spectral density, we can extract a new collective mode with frequency Ω 1 that couples toB 0 with strength G 1 and to a new residual exciton environment with spectral density generated by Eq. (31) (cf. Fig. 6b). This process can be iterated indefinitely, thereby generating a one-dimensional chain of coupled modes. We already know that the first residual spectral density J (ω) dω is the coupling strength between the modes n and n − 1, and The corresponding mode frequencies are Ω n = J (n−1) res (ω) dω, and the Hamiltonian at the n'th iteration of the prodedure is thuŝ wherê i . b. When the reaction coordinate mapping is repeated, a new reaction coordinate, B1, is extracted from the residual environment. This reaction coordinate is then coupled to a new residual environment with operatorsB (1) i . c. After n repetitions of the mapping, the excitons are represented as a chain with n + 1 bosonic modes and a residual excitonic environment.
A possible strategy for making the non-Markovian time evolution tractable is to neglect the nth residual environment, thereby truncating the chain at the nth level, This procedure converges towards the exact result as n → ∞, and the number of chain links required to obtain an error that is low enough for the method to be useful depends on the specific system studied. For example, it was previously found that structures with Fano interference effects are very challenging to capture with a one-dimensional chain representation 82 . Furthermore, the error of a particular n-truncation depends on the time scale over which one is interested in the dynamics: If the chain is initially unpopulated, population will flow from the resonator mode through the chain. This is the relevant situation in the present case because the initial state is the thermal state at or below room temperature where there are no excitons. At longer times, a larger portion of the chain is explored by non-vanishing populations and thus more chain links are necessary to resolve the evolution. In the time evolution of the truncated chain system we also include dissipation of the electromagnetic field and excitonic line broadening effects of the reaction coordinate as discussed in Sec. IV B 2, by evolving the density operator of the resonant field and n chain modes with the master equation Although beyond the scope of this article, we note that the one-dimensional chain is amenable to numerically efficient and exact renormalisation-group methods [83][84][85][86][87][88] .
Here we shall discuss how the time evolution within a limited time window can be accurately captured by using a truncated chain mapping with a finite number of sites.

Benchmarking
To test the precision of a given truncation n, we can consider the situation in Sec. IV A, where the drive is turned off and the electromagnetic resonator is initialised in a single-excitation state. In this case, we can benchmark the truncated chain expansion against the exact time evolution and evaluate the error. An example of such a comparison is presented in Fig. 7. The time evolution in Fig. 7a demonstrates that a chain with more links allows one to evolve the system further in time before the error becomes pronounced. This is further supported in Fig. 7b, which shows that the relative error of the truncated chain mapping, as defined in Eq. (63), decreases when more links are included. Here, the error of the Markovian master equation, Eq. (62) is indicated with a dotted grey line, demonstrating that it is possible to go below this error and thus resolve non-Markovian effects in the residual exciton environment. Additionally, Fig. 7c shows that the maximum population of the last chain link during the evolution time decreases monotonically with the number of chain links. To carry out this analysis, we have neglected excitonic line broadening effects by setting γ x and γ x to zero, since these effects are not compatible with the exact strategy described in Sec. IV A. We do, however, note that the influence of the residual exciton environment is expected to be less important when additional decay channels are present. As such, we should think of the error in Fig. 7b as an upper bound on the truncation error that is expected when excitonic line broadening effects are included.
For reference, the numerical benchmark calculations of the Markovian master equations and the chain-mapped master equation have also been performed with a larger resonator linewidth and at higher temperature of 300 K, which means that the phonon-induced broadening is more pronounced. These calculations can be found in Appendix E and show the same overall behaviour as the calculations presented here. However, the effect of the residual excitons is seen to be smaller, which is attributed to the increased dissipation into the other decay channels, i.e. through resonator losses and phonon-induced exciton decay.

Pulsed driving in the linear regime
One of the benefits of the chain-mapping technique is the simplicity of introducing driving with a timedependent amplitude, F (t). For the Markovian master equation derived in Sec. IV B, this would result in a decay term, K, which is explicitly time dependent. Here, we use the chain mapping to study the response of the exciton-resonator system to a short laser pulse that weakly perturbs the system. By making sure that the initially induced population of the resonant field is far below unity, the nonlinear exciton-exciton interaction can be neglected, while the interactions with the residual exciton modes are captured by the chain of bosonic modes. The Hamiltonian for the truncated chain in a frame ro- tating at the driving frequency iŝ and we shall focus on a Gaussian pulse amplitude of the form We also re-instate the excitonic line broadening effects described in Sec. IV B 2 with the rates γ x and γ x . In Fig. 8, the dynamics of the electromagnetic field (orange solid) and exciton reaction coordinate (blue solid) is plotted for two different lateral length scales of the electromagnetic resonator and compared to the corresponding evolution obtained by ignoring the residual excitons (dots). As the lateral length scale decreases, the exciton cutoff frequency ξ increases as 1/L 2 , effectively increasing the interaction between the reaction coordinate and the residual exciton modes. For L = 4 nm, we find a relative error of 0.8% in the resonator field evolution and 0.4% for the exciton reaction coordinate evolution, when the residual excitons are ignored. This means that we can justify the use of a Markovian theory with a single exciton mode when L > 4 nm for the given parameters. In the next section, we shall investigate the nonlinear response in this regime. As we shall see, there is a parameter regime where L is large enough (here above 4nm) to neglect the residual exciton modes, yet small enough that the nonlinear interactions become pronounced.

Pulsed driving in the nonlinear regime
To investigate the nonlinear response of the system, we include exciton interactions in the local exciton reaction coordinate through the termŴ 0 , as described in Sec. III B. To see the nonlinear response, we need to increase the amplitude A of the driving pulse. In practice, this means that we need to resolve a higher number of excitations in the system, whereby it becomes significantly more challenging to include a higher number of chain sites. At this point, therefore, we choose the lateral size of the mode to be sufficiently large that the residual modes can be neglected (L ≥ 4 nm for the parameters in Fig. 8).
As before, we show in Fig. 9 the time evolution of the system (solid lines) for different lateral optical confinement lengths, L. We compare the evolution with the linear response, as obtained by removing the interaction term W 0 from the Hamiltonian (dots, dashed lines), thereby uncovering the role of exciton-exciton interactions in the dynamical evolution. Due to the increase of the nonlinearity as L is decreased, a larger deviation between the full dynamics and the linear response is observed. The dominating nonlinear effect is a reduced energy transfer from the resonant field to the exciton reaction coordinate, which arises because the transition from one to two excitons is shifted away from resonance due to the nonlinear interactions. This is a signature of a polariton blockade process, i.e. the inhibition of multi-polariton excitations. The dynamics in Fig. 9 shows that there exists a parameter regime, where the lateral confinement length is sufficiently small that the nonlinear interactions influence the dynamics, yet large enough that the residual excitons can be ignored. In Ref. 89, we study the impact of the nonlinear interactions in terms of polariton blockade in detail and establish the conditions for reaching blockade. Most importantly, we find that polariton blockade is reached when the nonlinear interaction strength, W 0 , exceeds the polariton dephasing generated by γ x .

V. SEMICLASSICAL LIMIT
As an alternative to the reaction-coordinate approach to exciton-resonator interactions, we now develop a semiclassical description of the interaction of the electromagnetic field with the excitons in terms of the excitonic dielectric response. This allows us to connect the fundamental material parameters to the dielectric function, which has been experimentally measured for several materials 59,60 . Furthermore, it allows us to carry out an independent classical reference calculation, which we expect to agree with the result of our microscopic theory in the weak-excitation limit of linear response. This reference calculation thus serves as an important consistency check of the microscopic theory.

A. Exciton susceptibility
In a purely classical framework, we can model the electromagnetic response of a two-dimensional material in the plane z = z 0 as a thin polarizable sheet of thickness d with a relative permittivity distribution given by In cases where the sheet is illuminated by an incoming electromagnetic field, we can calculate the total electric field by use of the Lippmann-Schwinger equation. For the present analysis, we can confine the discussion to the case of normal incidence and consider incoming electric fields of the general form E in (r, ω) = E in (z, ω)n, where n is a unit polarization vector, taken to be linear for simplicity. The total field can then be calculated as the solution to the equation in which k 0 = ω/c is the ratio of the angular frequency to the speed of light, and G B (z, z , ω) = i exp[ik 0 |z − z |]/(2k 0 ) is the one-dimensional electric field Green function of the homogeneous background. For sufficiently thin materials, we can assume the integrand to be approximately constant, wherefore we can solve the equation to find In order to establish a link between the fundamental excitonic properties and the corresponding susceptibility, we consider the total field generated by weakly driving the excitons with a set of normal-incidence plane waves of the form where Z is the depth of the quantization volume. To calculate the total field, we use the full interaction Hamiltonian with non-rotating wave terms, Eq. (11), and follow Refs. 90-92 to derive a Lipmann-Schwinger equation for the vector potential operatorÂ(z, ω) = nÂ(z, ω), expressed in terms of the incoming vector potential opera-torÂ (0) (z, ω) = nÂ (0) (z, ω), of the form where V (z, ω) = V 0 (ω)δ(z − z 0 ), and we have ignored a possible source term for the vector potential, because we consider a scattering problem, where no excitons are initially excited. The Green's function entering Eq. (79) is related to G B as G(z, z , ω) = −G B (z, z , ω)/S. We can solve Eq. (79) to find and since the electric field is the derivative of the vector potential, we find that the same equation holds for the electric field, under the substitutionsÂ →Ê andÂ (0) → E (0) . By comparing Eqs. (77) and (81), we can directly see that the two expressions are equivalent, when k 2 0 dχ(ω) = −V 0 (ω)/S. This requirement yields the expression for the exciton susceptibility, To connect to susceptibility measurements in the literature, we note that the exciton line broadening due to non-radiative decay and phonon interactions should also be taken into account in the susceptibility. We include these effects by a complex shift of the poles of the susceptibility at ω = ±ω 0 into the lower part of the complex plane by the total line broadening Γ x = γ x + γ x . The resulting susceptibility, which serves as a model of what one would measure in linear-response measurements, then takes the form The poles of interest are both located in the lower part of the complex plane, as they should be due to causality of the electromagnetic response in the time domain, and one of the poles is at negative real frequencies, wherefore it can often be neglected. When more exciton transitions are present, a generalization gives an expression with multiple pole terms. Combining these pairwise and dropping small terms of order Γ 2 x , we find the familiar form where the oscillator strengths, f m , are given by With the assumption that the field polarisation is in the plane, such that n 2 x + n 2 y = 1, we can write the oscillator strength as This dielectric parameter f m has been measured for several 2D materials and provides a very useful means to determine the exciton Bohr radius from experiments 59,60 . For WS 2 , using the parameters in Table I

B. Reference calculation
We are now in a position to compare the microscopic theory to a semiclassical reference calculation. To this end, we consider a gold nanorod coupled to a monolayer of WS 2 and calculate the linear excitation spectrum when driving with an external laser at different frequencies. The nanorod is modeled as a cylinder with spherical end caps, as depicted in Fig. 1, a diameter of 30 nm, and a total length of 90 nm. For convenience, we define a coordinate system in which the nanorod is oriented in the x-direction, and its center is at the position (x, y, z) = (0, 0, 20 nm). For these calculations, we use a Drude permittivity model of the form withhω p = 6.9 eV andhγ = 0.2 eV. The numerical QNM calculations were carried out with the boundary-element method "MNPBEM" 93 (5)i. Here, 0 = 100 nm is a fixed length scale that is used to express the QNM frequency and norm in dimensionless units. In the single-QNM approximation, it can be shown that the normalization factor S c appearing in the derivations in Ref. 45 always takes the value S c = 1, which is the value that we used for the calculations in this work. The exciton spectral density, J(ω) and the residual spectral density, J res (ω) (green dashed), are shown in Fig. 10a along with the absolute value of the in-plane components of the field profileF c , as derived from the fundamental dipolar QNM and discussed in Appendix A. The dipolar nature of the field profile results in a non-trivial in-plane distribution of the electric field as shown in the inset of Fig. 10 for the case of z = z 0 . It follows from the analysis in Sec. III A that this is the field distribution defining the exciton reaction coordinate and the coupling constant G 0 , recall Eq. (36). Running a benchmark calculation as described in Sec. IV B 1, we find that the relative error obtained by neglecting the residual exciton modes is less than one in a thousand, wherefore we neglect the residual excitons in the analysis. We then calculate the excitation spectrum as the steady-state photon number, n ss (ω d ) := Tr â † câcρss (ω d ) whereρ ss (ω d ) is the steadystate density operator of Eq. (62) with constant driving (frequency ω d and driving strength F ). For these calculations, the driving strengthhF was set to a very low value of 36µeV to ensure that the system is in the linearresponse regime.
As a semiclassical reference calculation, we can also calculate the excitation spectrum through a solution of the classical electromagnetic problem, where the WS 2 monolayer is modelled as a sheet of thickness d = 0.618nm with a dielectric function corresponding to Eq. (84) with experimentally measured parameters 59 , h 2 f 0 = 1.9 eV 2 ,hω 0 = 2.014 eV, and the exciton linewidthhΓ x = 16.1 meV, calculated from Table II using T = 300 K. The reference calculations were done with MNPBEM and the same mesh for the nanorod as was used for the QNM calculations; additional scattering from the 2D material was included by the use of the appropriate background Green function 96 corresponding to a thin sheet with a single Drude-Lorentz pole. In practical experiments, the sample will rest on a substrate and there will be additional contributions to the optical response of the WS 2 -material and other similar corrections, which will lead to shifts in the resonance frequencies of the nanorod and possibly the 2D material. In Appendix C, we investigate the effect of a substrate on the electromagnetic response and show that the nanorod can always be tuned into resonance with the excitons by various means, for example by varying the length of the nanorods. In order to simplify the model and focus on the dominating physics, we have left out these effects here.
The driving laser field was modeled as an incoming plane wave with frequency ω d , and for each value of this frequency, the resonator excitation was measured as the field intensity in the middle of the nanorod, I(ω d ). In order to compare the field intensity I(ω d ) and the steadystate photon number n ss (ω d ), we normalise both to their maximum values and plot them together in Fig. 10b. The asymmetry in the semiclassical spectrum can be attributed to a frequency-dependent incoupling factor between the external driving field and the resonator field term, which can be derived from coupled-mode theory 44 , but which is not accounted for in the present approach. The remaining discrepancy is attributed to the nonretarded coupling and the approximation that only a single QNM is taken into account in the microscopic model. Importantly, we find that the calculated splitting of the spectrum in the two independent calculation methods differ only by 0.5%. In combination with the general qualitative agreement between the two spectra, we interpret this as a demonstration of consistency between the microscopic quantum model and the semiclassical theory based on measurements of the linear exciton susceptibility. In this limit of linear response, it is an interesting fact that one can also treat the problem from a purely electromagnetic point of view and model the response by use of two quasi-normal modes, as was recently presented in Ref. 97.
In closing, we emphasize that the general microscopic model is applicable also beyond the linear, semiclassical regime, when nonlinear effects and few-exciton statistics become important, as discussed in Secs. III B and IV C 3, as well as in Ref. 89.

VI. CONCLUSION
In conclusion, we have developed a microscopic quantum theory for the interaction between an electromagnetic resonator and excitons in a pristine sheet of 2D semiconductor material. In particular, by invoking a basis change of the exciton continuum, we have identified a collective exciton mode, termed the exciton reaction coordinate, that effectively accounts for the light-matter interaction. We have derived analytic expressions for the coupling strength between the resonant electromagnetic field and the reaction coordinate, thereby showing that it is independent of the lateral confinement of the field.
To calculate the dynamical evolution of the system, we have introduced and analyzed several Markovian and non-Markovian approaches and assessed their regimes of validity. Using these strategies, we have evaluated the importance of the residual exciton environment, which is coupled to the reaction coordinate. We find that the influence of the residual excitons becomes more pronounced when the lateral optical mode dimensions become smaller. In many cases, however, the residual excitons can be ignored altogether. For the extreme regime where the electromagnetic field is laterally confined to a characteristic length scale of a few nanometers, it becomes necessary to account for the residual excitons. We have developed an iterative chain-representation of the residual exciton environment, which is able to resolve non-Markovian effects and thus to go beyond the Markovian master equation.
We have also derived the linear dielectric response of the excitons, which allows one to connect the material parameters to the dielectric function and to consistently interface the microscopic theory with a semiclassical approach. Furthermore, we have calculated the nonlinear interaction strength of the excitons within the reaction coordinate and found that it scales as the inverse area of the electromagnetic field in the 2D material, meaning that laterally confined electromagnetic fields lead to stronger exciton-exciton interactions. In this context, we have found that there exists an interesting parameter regime, where the lateral confinement length scale is large enough that the residual excitons can be ignored, but small enough that nonlinear effects are significant.
this Appendix, we apply the approach to the present situation, where the resonator and exciton reaction coordinate are treated as an open quantum system, which is coupled to an environment consisting of the residual exciton modes. Starting from Eq. (57), we decompose the interaction Hamiltonian asĤ SR =h(Ŝ 1R1 +Ŝ 2R2 ), whereŜ To simplify the expression for the interaction-picture time evolution, we introduce the eigenstate-projected system operatorŝ with i ∈ {1, 2}, where Π(E) is the projector onto the system subspace with eigenenergy E with respect toĤ S , We distinguish betweenŜ † i (ω) and [Ŝ i (ω)] † , such that the former expression refers to E −E=ω Π(E)Ŝ † i Π(E ) and the latter to . The interaction-picture time evolution of the projected system operators is then Using the completeness of the eigenstates ofĤ S , we find and the dissipator due to the residual excitons can thus be written as where Λ ij (τ ) = Tr R R i e −iĤ R τ /hR j e +iĤ R τ /hρ0 R is a residual excitonic correlation function. Here,ρ 0 R is the initial density operator of the residual exciton environment, which is taken to be the vacuum state as a good approximation to the thermal state of a semiconductor. As described in Ref. 75, a so-called secular approximation is enforced by keeping only terms with ω = ω in the summation, which is justified by the fact that a factor of exp{i(ω − ω )t} appears in the sum for the interactionpicture time evolution of the reduced density operator; if ω = ω , the exponential is assumed to average out to zero. In addition to simplifying the expression, the secular approximation ensures that the dynamics generated by the master equation is completely positive and trace preserving 75,98 . Noting that, due toρ 0 R being the vacuum state, the only nonzero correlation function is Λ 12 (τ ) = i>0 |λ i | 2 exp{−i(Ω i − ω d )τ }, the secularised residual exciton dissipator becomes where Γ res (ω) = 2 Re The second term in Eq. (B8) amounts to a shift of the resonance energies and will be neglected here. The remaining part describes exciton dissipation with a rate that can be written in terms of the residual spectral density as in Eq. but the general anti-crossing trend in the curves is unchanged, since the dominating physics is still that of two strongly coupled harmonic oscillators. In order to substantiate the claim that the coupling strength G 0 is largely independent of the lateral confinement length scale, we explicitly compare the value of G 0 for gold nanorods of different lengths ranging from L = 75 nm to L = 95 nm, but otherwise identical to the one that was investigated in Sec. V B, see also App. C. Fig. 12a shows G 0 as a function of the nanorod length, L, and the relative difference compared to L = 90 nm is shown in Fig. 12b (blue data points and lines). While the coupling strength increases as a function of length, this effect is mainly due to the change in the resonance frequency ω c , which decreases with increasing length, as seen in Fig. 11. In Eq. (36), a factor of 1/ √ ω c appears in the coupling strength. Thus, to make a meaningful comparison of the effect of the spatial mode distribution on the coupling strength, we should multiply the coupling strength by √ ω c to correct for the shift in resonance frequency. This comparison is shown in Fig. 12b with orange data points and lines and reveals that the change in resonator length by 20% generates a vanishing shift in √ ω c G 0 of around 1%.
Appendix E: Comparison of time-evolution methods with increased dissipation and decoherence In Sec. IV, three different approaches for calculating the time evolution of the exciton-resonator system were presented and compared, in order to assess their validity. For completeness, we present the same comparison calculations with the only difference that the resonator decay rate has been increased to 2hγ c = 20 meV, and the temperature has been increased to 300 K; the resulting phonon-induced exciton decay rate ishγ x = 7.7 meV, and the dephasing ishγ x = 8.4 meV. Figure 13 corresponds to Fig. 5  culation. We remind that excitonic line broadening due to phonon interactions is not included in these benchmark calculations, since the exact memory-kernel equation, Eq. (53), is incompatible with these effects.  Fig. 8, but with the resonator decay rate increased to 2hγc = 20 meV and the temperature increased to 300 K, leading to the phononinduced exciton decay ratehγx = 7.7 meV and decoherence ratehγ x = 8.4 meV.