Microscopic theory of nonlinear phase space filling in polaritonic lattices

We develop a full microscopic theory for a nonlinear phase space filling (NPSF) in strongly coupled two-dimensional polaritonic lattices. Ubiquitous in polaritonic experiments, the theoretical description of NPSF, remains limited to perturbative treatment and homogeneous samples. In this study, we go beyond the existing theoretical description and discover the broad scope of regimes where NPSF crucially modifies the optical response. Studying the quantum effects of non-bosonicity, cooperative light-matter coupling, and Coulomb blockade, we reveal several regimes for observing the nonlinear Rabi splitting quench due to the phase space filling. Unlike prior studies, we derive nonlinear Rabi frequency scaling all the way to the saturation limit and show that the presence of a lattice potential leads to qualitatively distinct nonlinearity. We concentrate on three regimes of NPSF: 1) planar; 2) fractured; and 3) ultralocalized. In planar saturation, the Rabi frequency decreases exponentially as a function of exciton density. For the fractured case, where excitons form a lattice with sites exceeding the exciton size, we discover fast NPSF at low occupation in the lattice. This is followed by slower NPSF as the medium becomes fully saturated. This behavior is particularly pronounced in the presence of Coulomb (or Rydberg) blockade, where regions of fast and slow NPSF depend on the strength of repulsion. For the ultralocalized NPSF, we observe the square-root saturation typical to the collection of two-level systems. Our findings can help describing recent observations of strong nonlinearity in heterobilayers of transition metal dichalcogenides where Moir{\'e} lattices emerge naturally [Nature \textbf{591}, 61 (2021)]. The theory also opens the prospects for engineering strongly nonlinear responses of polaritonic lattices with patterned samples, driving polaritonics into the quantum regime.

The second type of nonlinearity comes from statistical properties of matter excitations and corresponds to the nonlinear phase space filling (NPSF), which can be referred to interchangeably as a nonlinear optical saturation [56,59,84].Namely, creating two excitons in exactly the same state of electrons (e) and holes (h) is forbidden by their fermionic statistics.Similarly, a collection of two-level systems (TLS) can only be excited until further excitations are prevented by Pauli statistics [5].Thus, increased number of excitations leads to filling the available phase space (no room for creating new quasiparticles) [85,86], and effectively reduced light-matter coupling in remaining sites or area.At high powers NPSF can lead to the power-dependent quench (or collapse) of the associated Rabi splitting [24].For Wannier excitons (delocalized e-h pairs) this nonlinear mechanism was discussed already in the seminal paper of Tassone & Yamamoto [56], where first-order powerdependent contribution to the reduction of Rabi splitting was described.However, in III-V semiconductors this nonlinearity was not considered to be dominant [59].For Frenkel-type excitons (quasiparticles based on localized e-h pairs in molecular lattices) NPSF plays a major role [87,88], and is ubiquitous in various experiments [87,[89][90][91].Finally, recent results in transition 1. Sketch of an optical microcavity with polaritonic lattice as an active medium.The cavity is formed by a concave top mirror (e.g.fiber-based) and a bottom mirror as a distributed Bragg reflector (DBR) of high reflectivity.Polaritonic lattice is arranged as a patterned semiconductor with excitonic potential, as typically realized in heterobilayers of transition metal dichalcogenides (TMDC) with moiré potential.
In this work, we develop a unified treatment of nonlinear quantum optical effects based on phase space filling.Our theory is applicable to a wide range of excitonpolariton lattices, also in the presence of both Pauli and Coulomb blockade.We describe three distinct NPSF regimes being the planar, fractured, and ultralocalized regimes.For each case, we present an analysis and show that in the fractured case a sharp decrease of Rabi frequency in the low-density regime can be facilitated by the Coulomb blockade.Our theory can shed light onto recent experiments in moiré heterobilayers, and open a way for enhancing the nonlinear response.

II. THEORY OF EXCITONS COUPLED TO A CAVITY MODE
To develop the microscopic model of nonlinear phase space filling, we start by considering excitons (electronhole pairs) strongly coupled to a cavity mode.This is typically formed by distributed Bragg reflectors or fiberbased mirrors (see sketch in Fig. 1).
Specifically, we consider bound excitons being confined to minima of potential, thus forming a polaritonic lattice (Fig. 2).A system Hamiltonian with N s lattice sites reads where ĉ † (ĉ) is a field operator creating (annihilating) a cavity photon with energy ω 0 (we set ℏ = 1).The first term in Eq. ( 1) thus describes the energy of the cavity mode at normal incidence (small photon momenta).The second term

Ns i=1
ĤX i describes the available quantum states on each lattice sites, labeled by index i.This term may also be viewed as a lattice of identical quantum emitters which can host up to ℓ excitons, depending on the confinement potential and interaction between excitons.Therefore, the Hamiltonian of each separate sites can be written as where 1 is the identity operator in the Hilbert space of a single site, and In the above, |n⟩ denotes the n-exciton eigenstate with energies ε n in an emitter.The ground state of a emitter is |0⟩ (no exciton) with energy ε 0 = 0, see Fig. 3.
Intuitively, this model is analogous to a lattice of decoupled harmonic oscillators.However, we note that we assume a general form of the energy-level structure of ĤX i to account for exciton-exciton interaction and arbitrary confinement potential, introducing some degree of unharmonicity.With this, we make sure that the microscopic model can describe excitons of different types (Wannier and Frenkel), as well as remaining valid for generic quantum emitters with ℓ + 1 levels.
Finally, the third term in the system Hamiltonian (1) describes the strong light-matter coupling, where Ω 0 is the bare exciton-photon coupling strength.This term defined the transition between the emitter's energies levels by absorbing or emitting a photon.These transitions at site i is described by the excitonic ladder operator as where r n ≡ |⟨n + 1|x † |n⟩| 2 with x † being the field operator for creating an exciton.This exciton transition matrix element r n describes the effective rate for creating an additional exciton in a n-exciton state, |n⟩.It is determined by the microscopic details of a n-exciton state.We note that, although our polaritonic lattice model Eq. ( 1) is similar to an atomic system, the optical transition of each emitter [Eq.( 3)] only takes place between the adjacent energy levels, |n⟩ and |n + 1⟩.Namely, the optical multiexciton processes are forbidden.In contrast to atoms, the transitions between levels do not have such a constraint and they are determined by optical selection rules.Importantly, excitons are composite bosons with a nontrivial n-dependence of r n .In contrast to atomic level transition, each atomic level cannot host more than one electron.This corresponds to r n = δ n,0 in Eq. ( 3), where the only allowable transition is between states |0⟩ and |1⟩.This is the key difference between the two systems which leads to distinct optical saturation effects qualitatively.The transition matrix element r n depends on the underlying physics of the quantum emitter (exciton).The quantitative description for r n is presented in Section III.
In the lattice of identical emitters, the exciton created by a uniform cavity mode must preserve the translational

Moiré pattern in TMDC heterostructure
Three regimes of nonlinear phase space filling: FIG. 2. Regimes of NPSF.Here we sketch different regimes of nonlinear phase space filling (saturation) in polaritonic lattices.As a prominent example, we present moiré heterostructures with different twist angles, demonstrating the planar, fractured, and ultralocalized NPSF at progressively increased twist angle.

Quantum emitters:
3. Quantum emitter.Each moiré cell of a bilayer can be viewed as a quantum emitter that hosts excitons in its excited state.The quantum states of the emitter with energy ε0, . . ., εn are denoted by |0⟩, . . ., |n⟩ in the diagram.
symmetry in the lattice.This implies that the photon can only couple to the lattice collective mode (bright state) in the form [5,86] This operator creates a coherent lattice excitation where the probability of finding an exciton is uniformly distributed across the lattice instead of sitting on a particular site.It corresponds to a collective exciton mode that participates in SC.This allows us to rewrite the coupling Hamiltonian, meaning the coupling only acts in the subspace of bright states.The collective excitonic quasiparticle B † is exactly the mode that couples to light in optically active materials, with the corresponding √ N s enhancement for Frenkel excitons [1,87] and exciton area-based enhancement of light-matter coupling for Wannier excitons [3].Importantly, the collective exciton mode B † also acts like a quasiparticle with well-defined particle numbers, but it has a peculiar statistical property.In the low ratio of total number of excitons N to N s in the lattice (N/N s , exciton number per site), the B † quasiparticle excitations obey statistics similar to the bosonic one.However, in a large N/N s regime, this quasiparticle excitation may strongly deviate from the Bose statistics due to the composite nature of exciton and blockade effects arising from the Coulomb interaction between excitons.This non-bosonicity determines how the phase space is depleted by B † -excitations.This will eventually appear as a nonlinear correction to the light-matter coupling of the lattice.
To see the relation between the phase spacefilling effect and nonlinear light-matter coupling, we consider in the SC regime that the two states {( B † ) N +1 |∅⟩, ĉ † ( B † ) N |∅⟩} hybridize and form a quasiparticle -polariton.Here, we defined the global ground state |∅⟩ = |0⟩ c Ns i=1 |0⟩ as the product of emitter ground state (|0⟩) and cavity (|0⟩ c ) ground states.To investigate the nonlinearity of the polariton in the lattice of emitters, we construct the polaritonic Hamiltonian as block diagonal matrices that couple states with N + 1 excitons, containing B † -excitations in total. Here is the total energy of N collec-tive excitons and we have introduced the normalization factor of the N -exciton states [85,102] This ensures the states are properly normalized such that the matrix elements in Eq. ( 6) retain the correct physical meaning at increasing densities (as phase space filling increases).The off-diagonal elements then can be written as [103] Ω meaning that the light-matter coupling Ω N (effective Rabi frequency) depends on the number of excitons.This introduces the nonlinearity in the system in the form of nonlinear phase space filling.
To study the effect of NPSF in full generality, we develop a strategy for calculating the normalization factor Eq. ( 7) for a generic structure of onsite excitations and arbitrary N s .We do this by using a multinomial expansion.The compact expression of the normalization factor reads (9) where P (n i ) = r 0 . . .r ni−1 with P (0) ≡ 1 (see Appendix A for the derivation).The renormalization factor in Eq. ( 7) completely determines the NPSF of the Hamiltonian (1).We note that while F (Ns) N was evaluated before in some limiting cases, such as Wannier exciton [85,102] and Frenkel exciton (two-level emitter) [86], here we present the generic non-perturbative treatment.This is imperative for accessing the NPSF in the deep saturation regime, as demonstrated in the next sections.
To summarize, we have developed a general theoretical framework for the exciton-polariton on a lattice which is described by Eqs.( 6)-( 9).This reduces the study of the NPSF problem to specifying r n and |n⟩, which are ultimately determined by the properties of the n-excitons state of an emitter.As we will discuss later, r n is the key parameter that gives rise to a different F (Ns) N , Ndependence of the Rabi frequency, and thus enable various qualitative modes of NPSF.Before carrying out the full analysis, it is instructive to consider two special cases and recover known results.First, let us consider |n⟩ being the state of purely bosonic excitations.In this case, the transition matrix element is r n = n+1 (bosonic enhancement factor) and the number of available states is ℓ → ∞.Using this assumptions in Eq. ( 9), we recover the usual bosonic normalization factor F (Ns) N = N !for N -boson state.Next, let us then consider the opposite limit corresponding to two-level systems (fermionic limit).This is modeled by a quantum emitter with r 0 = 1, r n≥1 = 0.With this, Eq. ( 9) reduces to the quantum limit [86] with In the rest of the paper, we demonstrate that the polariton Hamiltonian in Eqs. ( 6)-( 9) allows us to go beyond these two limits.We consider transition matrix elements that are neither purely bosonic nor purely fermionic.This is controlled by the confinement, composite particle properties, and interactions.Our theory reveals the behavior of F (Ns) N for a polaritonic lattice with strong NPSF (deep saturation regime).This regime of SC has not been explored before.In the following, we investigate this regime and discuss the relevant system in moiré structure of 2D materials.

III. INTRODUCING LOCALIZATION OF EXCITONS AND POLARITONIC LATTICES
In the previous section, we presented the general theory of excitons in a cavity and introduced the transition matrix elements r n that characterize the polaritonic system [Eq.( 3)].As we mentioned previously, the form of r n has important implications for the nonlinear phase space filling of the lattice.As described in Eq. ( 3), the transition matrix element r n is a measure of the rate for creating an additional exciton-like quasiparticle in a nexcitons system.This quantity depends on the specific correlation effects between excitons.Calculating it for general cases is a highly non-trivial many-body problem.Therefore, in this section, we consider two microscopic mechanisms that induce NPSF.One is due to the correlation of the fermionic statistics (Pauli blockade) where the electron and hole from excitons cannot occupy the same state in phase space.The other mechanism is due to the strong exciton-exciton interactions (Coulomb blockade) where the high-energy states are effectively projected out in the phase space.In the Pauli blockade, we can evaluate r n exactly, and for the Coulomb blockade we use the phenomenological approach that allows to keep the system tractable in the presence of nonlinearity.

A. Pauli blockade
To investigate the Pauli blockade of a polaritonic lattice, we begin with the simplest case of N s = 1 and extend the consideration to excitons with a spatial shape.For a single site (or emitter), we construct the n-exciton ground state as where the exciton creation field operator is Here, Ψ ν (r e , r h ) is the exciton wavefunction, with ν being the state index (quantum number), r e and r h are the electron and hole position vectors in real space, and a † re and b r h are the electron and hole field operators, respectively.
Due to the finite sample size, there is a limited amount of quantum numbers to assign to each electron and hole composing the exciton.Moreover, two different fermions cannot be labeled by the same quantum numbers (Pauli exclusion principle), meaning that the single site cannot host an unlimited number of excitons.The site will saturate and exhibit a nonlinear optical response as the exciton density becomes higher and higher.In Eq. ( 10) the effect of Pauli blockade appears in the normalization factor F (1) n .We write the latter in the recursive form, where we introduced the Pauli scattering terms σ m (see Appendix B and Refs.[85,103]).The Pauli scattering is defined as an overlap between excitonic wavefunctions, Intuitively, this term may be understood as the degree of bosonicity of the exciton.For the pure boson limit we have σ 1 = 1 and σ m>1 = 0, leading to F (1) n = n! as a standard bosonic state normalization factor.In the Frenkel limit, we have σ m = 1 for all m, which yields F (1) 1

= 1 and F
(1) n>1 = 0.In between these two limits σ m is a monotonically decreasing function of m with σ m ≤ 1. Away from the bosonic limit, the effect of nonlinear phase space filling is always present, thus lowering the effective probability for creating more excitons in confined regions.That is, F n /n! is a monotonically decreasing function of the integer index n.The effect of the Pauli blockade can be seen as a modification of the transition matrix element r n defined through Once r n is set for a single site, we can construct the N ssite of the polaritonic lattice by using Eq. ( 9) with r n specified in Eq. ( 14).In this case, calculating matrix elements in Eq. ( 6) by using Eqs.( 9) and ( 14) is not convenient, since the number of available excited states (ℓ) may be large.Instead of using multinomial expansion, we derive the recursive formula for F (Ns) N (see Appendix B for full derivation): However, we note that when using Eq. ( 15) for numerical calculation, one needs to ensure high accuracy for each terms in the sum due to large cancellations (we achieve this by developing a semi-symbolic treatment of expressions when only the final evaluation is numerical).As compared to the single site case in Eq. ( 12), the multisite result has a reduced Pauli scattering term σ m /N m−1 s which takes into account two contributing effects.One is the Pauli blockade that we have already described, relating the effect to the statistical properties of the particles.Second, we have the real space enhancement factor, meaning there are more sites to place our emitters.

B. Coulomb blockade
Next, we note that placing two excitons together may be impossible due to the fact that their constituents interact with each other via Coulomb potential.The enhancement of onsite exciton-exciton (X-X) Coulomb repulsion due to weak screening in two-dimensional structures translates into a strong Kerr shift for the excitonic mode.As the exciton energy is shifted out of resonance, this effectively reduces the transition matrix element between excited states, resulting in similar implications as the Pauli blockade.The main difference with the Pauli blockade is that the number of available excitations per site is sharply depleted at finite ℓ.The cut-off value ℓ depends on the depth of localization potential and the cost of Coulomb energy by creating additional excitons at the same site.However, we note that the Coulomb blockade is a complex correlated problem since the onsite X-X Coulomb repulsion and the maximum number of excitons per site ℓ implicitly depend on the n-exciton rather than one-exciton wave function as in the Pauli blockade problem.Calculating the transition matrix element r n microscopically is beyond the scope of this paper, as it reduces to exclusively numerical modeling, while here we concentrate on predominantly analytical treatment.Nevertheless, our general theory in Eq. ( 1) is still applicable, and we introduce the Coulomb blockade effect as the phenomenological reduction of the transition matrix element at increasing n.This may be modeled by corresponding to the rate for driving the transition between energy levels with detuning δ n = (n − 1)U/Ω 0 (dimensionless) due to the Kerr shift.Here, U is the strength of X-X interaction and n is the number of excitons.In this case the Pauli statistics plays no role, where the exciton is considered as an elementary boson.However, we will see later that different blockade mechanisms eventually lead to similar qualitative results.
Although the form of Coulomb blockade in Eq. ( 16) is considered for this paper, we remark that the physics of saturation effects due to X-X interaction is very rich.Particularly, the interaction depends on both the exciton's wavefunction and the screened Coulomb potential which have important consequences for forming the nexciton state |n⟩.The details of |n⟩ ultimately impact the scaling of the saturation effect.For instance, as we compare s-wave and p-wave exciton (Rydberg state), in the case of dipole-dipole interactions, the long-range part of the effective X-X interaction can set a blockade radius, thus effectively enhancing the saturation (we have now seen this in Cu 2 O polaritons [104].For s-excitons, the exchange interaction is dominant and generally short-range.However, it depends heavily on the details of wavefunctions and this leads to a qualitatively different scaling from Rydberg state.

IV. DIFFERENT REGIMES OF NONLINEAR PHASE SPACE FILLING
Previously, we have derived a microscopic expression for light-matter coupling and nonlinear phase space filling in the general form, accounting for multiple localized excitonic sites, the spatial structure of excitons, and X-X interactions.Now, let us analyze the behavior of NPSF for some qualitatively distinct cases, in particular driven by a sample geometry.The size of excitons is described by the Bohr radius a X .This has to be compared with the exciton localization length L, defined by the lattice potential.Depending on the L/a X ratio as well the number of sites N s , the NPSF contribution has different scaling, both in high and low occupation limits.We suggest the three regimes corresponding to: 1) planar; 2) fractured; and 3) ultralocalized NPSF.
We visualize the three sample geometries in Fig. 2. The planar case corresponds to an exciton that is delocalized over the entire sample (N s = 1; Fig. 2, left).The saturation effect in this regime was studied before for Wannier excitons in III-V semiconductors [56,59] and TMDC monolayers [24,105].In this case, the prior analysis was performed in the perturbative limit, where only the first nonlinear correction to Ω N is derived as ), where the beta-factor depends on the ratio of exciton area to the total area [56,59], and the secondorder correction was derived in Ref. [24].
Another regime corresponds to the opposite limit, where excitons are ultralocalized on the lattice of many sites, and the localization length is comparable to the size of the exciton (Fig. 2, right).In this case, the exciton behaves like a Frenkel exciton [86], and we note that similar behavior can be attributed to trion polaritons [24,94].The Rabi frequency of the ultralocalized case scales as a square root of the deviation from single occupation case.
Finally, we reveal that for the intermediate localized length, the fractured regime can be realized.In this case the analytical form of Rabi frequency is not known, and we will show that the NPSF scaling behavior of this regime is rather nontrivial.Namely, it is not a simple interpolating result between planar and ultralocalized cases.This may shed light on various polaritonic experiments where lattice sites can host multiple but finite number of excitons.
Let us now proceed with calculating the scaling of NPSF for the three regimes.To account for the effect of Pauli blockade we consider an excitonic wavefunction in the Gaussian form where R = mc M r e + mv M r h is the center of mass, and r = r e − r h is the relative position of an electron-hole pair.Here, M = m c + m v is the mass of exciton, and m c and m v are the masses of electron and hole, respectively.We note that the similar ansatz was chosen for studying phase space filling of a quantum dot [103], and that the Gaussian exciton shape is particularly suitable for TMDC heterostructures [95,101,106] due to the quadratic scaling of an interlayer potential.We remind that L is the localization length of the lattice potential and a X defines the size of the exciton.The exciton localization length can be characterized by the wavefunction spread in the exciton's center-of-mass motion which depends on the confinement potential.This can be dictated by an excitonic disorder, or induced by twisting in a TMDC bilayer forming the moiré potential [107,108] (visualized in Fig. 2).The exciton Bohr radius a X depends on the electron-hole screened interacting potential (defined as intralayer and interlayer potentials for TMDC bilayers).While other choices of wavefunction are possible, Gaussian shape allows evaluating many scattering terms efficiently, and we do not expect the shape modification to induce qualitative changes.
We proceed by evaluating Pauli scatterings for defining r n .Substituting the wavefunction from Eq. ( 17) in Eq. ( 13), we obtain the analytic expression of σ n .The compact expression reads (see details in Appendix C) where γ c,v = (m c,v /M )(a X /L) [103].We report the effects of the Pauli blockade in Fig. 4. We run the numerical simulation with various L/a X ratios and m c = m v by solving Eq. ( 8).In Fig. 4(a), we demonstrate the effects of saturation in the single site limit.This case corresponds to the single parabolic potential well with excitons that couple strongly to the cavity mode.In the limit of delocalized particles (light blue line), it approaches the ideal bosonic limit (gray dashed line), where NPSF is not present due to Bose statistics.The different curves and bullet points are summarized in the legend in Fig. 4(b).
In the other limit of particles localized over the full sample size (dark blue curve with L = 0.5a X , closest to the bottom), particles exhibit "close-to-fermionic" behavior, 1 1 Herein, the term fermionic is used in a loose sense, as the quasiparticles (excitons) do not exhibit anticommutation.Technically correct definition, albeit unusual, is a spin-1/2 or paulino statistics.The black dashed line is the bosonic limit, where the exciton is treated as an elementary boson.(b) NPSF for Ns = 20 sites.We can see the transition from Frenkel (L/a X = 0.5) to bosonic limit by changing the localized length L. Three qualitative regimes of the NPSF are labeled.The planar regime smoothly transitions over to the ultralocalized regime via the fractured regime (blue-shaded region).
where the sample is saturated with more than one excitation.In the intermediate region, 0.5 < L/a X < 10, one can see a smooth transition from bosonic-to-fermionic behavior.Also, the renormalized Rabi frequency in this intermediate region may be well approximated by the exponential decay with the exciton packing fraction in the sample η = N (a X /L) 2 , where we have used v ) in Eq. ( 18).Note that this result follows from the definition (8) and the exponential form of the normalization factors [85,102], which is valid for very high orders in η.The exponential scaling in Eq. ( 19) can be also recovered as a continuation of the diagrammatic expansion [24].
Next, we consider a lattice with many sites.As light couples to the lattice with N s = 20, the saturation effect exhibits a rather nontrivial behavior shown in Fig. 4(b).In this case, we observe similar behavior for large localization length with L = 10a X and recover the expected bosonic limit for point-like particles.As the exciton size increases (L ≃ 0.6a X ), we observe a kink in the proximity of N ≈ N s due to the transition from the low occupation low-density regime with a fast saturation rate to the high occupation regime with a slower rate.This behavior may be understood as follows.As N < N s , the blockade of available phase-space is very effective which resembles Frenkel exciton (see the blue curve with L = 0.5a X ).At the point where N = N s , the N excitons almost occupies the entire phase space with less available space.Therefore, the saturation due to the newly created excitations in N > N s regime is less effective.We can intuitively describe it as "fill all sites by at least one exciton first, and only then fill the rest".
To see the nontrivial changes in the saturation rate, in Fig. 5 we plot a derivative of the Rabi frequency with respect to the number of particles, dΩ N /dN (assuming Ω N is a continuous function in N ).We observe that in the planar regime (large L/a X ) the saturation rate is relatively small and constant in N .Once the exciton becomes localized with L ≤ a X , the saturation rate exhibits strong N -dependent behavior.In particular, in the fractured regime (L/a X ≈ 0.6, yellow curve in Fig. 5), we can see clearly the transition of fast-to-slow saturation rate from low-occupation (N < N s ) to high-occupation (N > N s ) regime.In the dilute limit (N ≪ N s ), the saturation rate can be characterized by a constant β 1 .As demonstrated in the inset of Fig. 5, this β 1 constant grows exponentially as the exciton localization length decreases.
As the localization length becomes smaller (L = 0.5a X ), the lattice localization eventually reaches the Frenkel limit [86].In this limit, the Rabi frequency can be calculated as FIG. 6. Coulomb blockade.The NPSF under the Coulomb blockade shown as a function of exciton occupation, for varying number of exciton states ℓ and X-X interaction strength U .We model the transition matrix element using Eq. ( 16) with ℓ ≤ 7. Upper, middle and lower panels are plotted by using the different onsite X-X repulsion as U/Ω0 = 0.5, U/Ω0 = 1, and U/Ω0 = 2, respectively.
which exhibits the square root scaling in N .We note that this expression follows from the normalization prefactor that is valid for all N and N s , and includes nonperturbative corrections.We stress that the exponential and square root scalings are qualitatively distinct from the conjectured saturation dependence Ω N = Ω 0 / 1 + N/N s used in some studies [95,100].Now, let us add the Coulomb blockade.In this case, the number of excitons at each site is limited due to the finite confinement potential and Coulomb repulsion.Also, once many excitons occupy a single site, it becomes more difficult to drive the transition from |n⟩ to |n + 1⟩ due to the Kerr shift of the excitonic mode.Using the blockade model in Eq. ( 16), we study the nonlinear phase space filling effect on the Rabi frequency.The results are shown in Fig. 6.Similar to the case of Pauli blockade (Fig. 4), the saturation curves exhibit similar behavior in each distinct regime.In this case, the fractured saturation may be realized as ℓ = 2, depending on the strength of Coulomb repulsion.The kink-like feature may not appear near N ≈ N s if the Coulomb blockade is weak such that the r n does not decrease fast enough in large n, see Fig. 6.Moreover, the quench of Rabi frequency nearing the full saturation (Rabi collapse) behaves closer to the ultralocalized regime.
Summarizing the discussion in this section, the nonlinear phase space filling due to excitons localization can be divided qualitatively into three different regimes (Fig. 2).Particularly, as the excitons started to localize at a different region of the lattice, NPSF enters the fractured regime.In this regime, the Rabi frequency of the polaritonic lattice is strongly renormalized.The phase space is saturated at a higher rate in the low (N < N s ) density regime rather than in the high (N > N s ) density regime, exhibiting an exponential tail.

V. SATURATION IN TMDC MOIRÉ LATTICE
In this section, we concentrate on drawing a connection between the proposed microscopic modeling of polariton lattice saturation with the recent exciting results for bilayers of TMDC.In the case of TMDC bilayers, the lattice ordering appears naturally, and we consider several possible contributing effects.Specifically, by twisting a TMDC bilayer one can drastically modify its optical properties [107][108][109].One of the salient features of such a bilayer is the formation of modulating optical absorption with a moiré period.Additionally, the inhomogeneity due to the moiré pattern can also generate a potential landscape that results in the localization of excitons in a moiré cell [108,110].Moreover, the variation of twisting angles and domain relaxation can lead to the formation of a lattice structure that is composed of smaller moiré patches.Therefore, twisted TMDC bilayers form an ideal system to realize near-perfect exciton-polaritonic lattices.moiré excitons in a twisted heterobilayer TMDC is the subject of current intense research [106,108,[111][112][113][114].In particular, the recent experiment [95] for the TMDC heterobilayers has found strong nonlinear saturation effects in moiré exciton-polariton.Experimental results suggest that there are several regimes of saturation where the saturation rate changes qualitatively.This was conjectured to originate from the Coulomb blockade due to the strong dipole-dipole interactions of moiré excitons, but a full explanation is missing.
We approach this system with the developed theory of NPSF.To investigate the optical nonlinearity of this twisted bilayer we use the lattice model in Eq. ( 1).Assuming m c = m v and using M = 0.71 m 0 from Ref. [95], we plot the saturation curve for Rabi frequency due to Pauli blockade and compare it to Coulomb blockade with Ω 0 = 10 meV and U = 30 meV.Using the estimation of the moiré cells density in Ref. [95], being 6 × 10 4 µm −2 , this gives the estimated exciton density in the horizontal axis in Fig. 7.As we can see, the Pauli blockade shows a rather smooth saturation in this bilayer system.In contrast to the Coulomb blockade, we can see a clear kink near N ≈ N s .We see that this nontrivial Rabi frequency renormalization qualitatively resembles the saturation behavior observed in Ref. [95].Namely, the saturation rate is fast in the low-density exciton regime where N < N s , and the rate reduces in the N > N s case.
We note the kink-like feature in the saturation curve appears in the higher density as compared to the density estimated in the experimental study.One possible reason may be that our calculation assumes a perfect correlation of exciton density with the laser power.In reality, the relation between them may be complicated due to the presence of disorder such as residue strain and can be impacted by the lattice reconstruction [108].This results in a smaller effective area of light-matter interaction.Furthermore, these effects also lead to the light interacting with the dark collective modes.Therefore, deducing the actual excitonic density in the sample may not be an easy task, since the number of excitons cannot be counted directly in the experiments.Nevertheless, keeping this as a hypothesis, we suggest that the experimentally observed nontrivial renormalization of Rabi frequency may be the manifestation of fractured NPSF which is facilitated by the Coulomb blockade.

VI. CONCLUSIONS
In summary, we developed a non-perturbative microscopic theory for describing nonlinear optical effects arising from the phase space filling.While the developed quantum theory can be applied for many systems at strong coupling, we concentrate on the lattice geometries in the limit of strong NPSF (deep saturation).We unified different phase space filling mechanisms such as Pauli and Coulomb blockade by writing them in term of transition matrix elements between the excitations of an emitter, r n .First, we described the effect of Pauli blockade arising from the finite system size and accounting for the exciton shape.This allows describing the change of NPSF from the planar quasi-bosonic case all the way to the case of two-level system saturation.We generalized the considerations to multiple sites, and studied the NPSF behavior between planar and ultralocalized (Frenkel) limit.Intriguingly, we find a distinct regime of fractured NPSF, where the Rabi frequency decrease has a kink-like feature at the number of sites, while the full saturation at high occupation reveals an exponential tail.This behavior is particularly pronounced in the presence of Coulomb blockade effects.Looking into specific examples, we analyzed the NPSF effects for moiré lattices of TMDC heterobilayers, suggesting that fractured NPSF may be relevant to the recent experimental observations of strong nonlinearity [95].
Looking into the future, we note that the understanding of nonlinear phase space filling can help engineering moiré structures, or rather patterned samples, such that the nonlinearity is maximized.This will allow to push further the limits of quantum polaritonics [29,38,39,94].
As for potential open questions, in our polaritonic lattice model we considered only the NPSF of a uniform collective excitation (bright state) [Eq.( 4)] coupled to cavity photons.However, cavity photons can couple to other non-uniform collective excitations (dark states) in the presence of lattice inhomogeneity arising from disorder, strain [115], lattice reconstruction [116], and non-uniform cavity modes.This shall further enrich the NPSF behavior leading to unique nonlinear effects.Also, in this paper we did not take into account the tunneling between sites.As indicated in Ref. [95], the moiré band dispersion of the exciton center-of-mass motion is not very flat.This implies that excitons tunneling between localized regions in the moiré lattice may not be negligible and can be tunable.Another consideration in Ref. [64] is that the light can induce changes in the exciton radius.This effect may add an extra contribution to the nonlinear saturation, and affect the Pauli and Coulomb blockade that we consider in this paper.We believe it is an interesting topic for future studies.It is also interesting to generalize our theoretical approach to investigate the saturation effect of moiré trions [117][118][119].Studying the trion-polariton NPSF in the lattice geometry and accounting for particle attraction [120] is an important next step.
Assuming that the light-matter interacting Hamiltonian of a single emitter as where the exciton in i-emitter is labeled by the subscript, we obtain Eq. ( 1) by sandwiching (B3) by the identity operator in Eq. (B2).Namely, To derive the transition matrix element for the Pauli blockade, we calculate Here, ∆ 1 = [x 0 , x † 0 ] and we introduce the notation and rewrite Eq.(B7) into To proceed further, we derive the iteration formula for Π Using this iterating formula and setting F N <0 = 0, we obtain Eq. ( 12) in the main text.

Appendix C: Calculation of σn
In this section, we explore the property of σ n for a single emitter (e.g.quantum dot [103]) in Eq. (B6).Here, we only focus on one emitter.Therefore, we drop the emitter index i for simplicity.To calculate σ n Eq. (B6), we use the commutation relation where ν is the quantum number of the exciton (ν = 0 being the ground state).The non-bosonicity is with σ 2 = Λ 00 00 .Similarly, it is not difficult to show that where σ 1 = 1, σ 2 = Λ 00 00 , and To do the integration in σ n analytically, we estimate the exciton wavefunction by using a Gaussian function.Then, the integration reduces to The eigenvalues of A can be obtained exactly by using a tight-binding approach (nearest-neighbor hopping with two sublattices).The eigenvalues are with k = 0, . . ., n − 1, Finally, we integrate out r ei and r hi .This gives as a closed expression for the Pauli scatterings in the finite size system.
FIG.1.Sketch of an optical microcavity with polaritonic lattice as an active medium.The cavity is formed by a concave top mirror (e.g.fiber-based) and a bottom mirror as a distributed Bragg reflector (DBR) of high reflectivity.Polaritonic lattice is arranged as a patterned semiconductor with excitonic potential, as typically realized in heterobilayers of transition metal dichalcogenides (TMDC) with moiré potential.

FIG. 4 .
FIG. 4. Pauli blockade and the dependence.The dependence of the (renormalized) Rabi frequency ΩN (Ns, a X ) = ΩN / √ N on number of excitons N .(a) NPSF effect of a single localization site (Ns = 1) for different ratios of L/a X .L is the exciton localized length, and the values are shown in the legend in (b).The black dashed line is the bosonic limit, where the exciton is treated as an elementary boson.(b) NPSF for Ns = 20 sites.We can see the transition from Frenkel (L/a X = 0.5) to bosonic limit by changing the localized length L. Three qualitative regimes of the NPSF are labeled.The planar regime smoothly transitions over to the ultralocalized regime via the fractured regime (blue-shaded region).