Multistable excitonic Stark effect

The optical Stark effect is a tell-tale signature of coherent light-matter interaction in excitonic systems, wherein an irradiating light beam tunes exciton transition frequencies. Here we show that, when excitons are placed in a nanophotonic cavity, the excitonic Stark effect can become highly nonlinear, exhibiting multi-valued and hysteretic Stark shifts that depend on the history of the irradiating light. This multistable Stark effect (MSE) arises from feedback between the cavity mode occupation and excitonic population, mediated by the Stark-induced mutual tuning of the cavity and excitonic resonances. Strikingly, the MSE manifests even for very dilute exciton concentrations and can yield discontinuous Stark shift jumps of order meV. We expect that the MSE can be realized in readily available transition metal dichalcogenide excitonic systems placed in planar photonic cavities, at modest pump intensities. This phenomenon can provide new means to engineer coupled states of light and matter that can persist even in the single exciton limit.

Strong light-matter interaction can provide a versatile platform for dynamically controlling quantum matter [1]. A striking example is the excitonic optical Stark effect [2][3][4][5]: off-resonant irradiation of an excitonic system, with frequency below the exciton transition energy, continuously blue-shifts the exciton transition to higher frequencies as the light intensity increases [4][5][6][7][8]. In contrast to the fixed Rabi splitting found for polaritons, that is independent of the intensity of light [9][10][11], the optical Stark effect is linear in the irradiation intensity. This dependence grants on-demand tunability of excitonic properties. In transition metal dichalcogenides (TMDs), Stark shifts are furthermore sensitive to light polarization, thereby enabling direct control over the valley excitons necessary for valley opto-electronics [8,12,13].
Here we propose that the optical Stark effect can take on a markedly different character when an excitonic system is placed in a nanophotonic cavity (Fig. 1 inset). In this setting, the Stark shift becomes a dynamical variable, with the cavity field taking on the role of the irradiating field that shifts the excitonic levels. In particular, when the excitonic and cavity modes are simultaneously pumped, the optical Stark effect can become multistable, exhibiting a hysteretic Stark shift that depends on the history of the optical drive. As we explain below and indicate in Fig. 1, this multistable Stark effect (MSE) arises due to a Stark-induced mutual tuning: the excitonic transition frequency (right panel) is sensitive to the cavity mode occupation, while the cavity resonance (left panel) is sensitive to the exciton population. When applied exciton and cavity driving fields are detuned from their respective bare transition frequencies, the mutual tuning sets up a feedback that shifts the exciton and cavity transition frequencies into resonance with their respective drives (dashed to solid lines, Fig. 1). This feedback leads to a highly non-linear, and multistable Stark effect.
The MSE features discontinuous transitions between wherein the exciton transition frequency is sensitive to the cavity mode occupation (and vice versa). When the excitons and the cavity mode are simultaneously pumped (downward arrows indicate the corresponding pump frequencies), the exciton and cavity transitions can shift into resonance with their drives (from dashed to solid curves). These populationinduced shifts generate the feedback loop that gives rise to the multistable Stark effect (MSE). (inset) A two-dimensional excitonic material such as a transition metal dichalcogenide can be readily layered on top of a planar nanophotonic cavity formed by a photonic crystal defect to achieve the conditions for realizing the MSE. multiple distinct steady states of the combined cavityexciton system, and can exhibit large discontinuous Stark shift jumps of order meV. Indeed, we find that the exciton population can take on multiple steady state values (Fig. 2a) with a hysteretic behavior that is controlled by a weak cavity drive far-detuned from the original exciton resonance. Further, the magnitude of the Stark shift jump from one stable state to another can be directly tuned by the drive that pumps the excitonic population. These mechanisms provide in-situ means of tailoring the switching behavior in the exciton/cavity system. We expect that the MSE can be realized in TMDs (e.g., WS 2 ) on currently available high quality factor planar photonic cavities [14][15][16][17] (Fig. 1 inset), even at low optical drive strengths of several to tens of kW/cm 2 . This platform provides new means of constructing hysteretic, nonlinearly coupled states of light and matter that can, in principle, persist even at the single exciton limit.
Stark-induced mutual tuning and nonlinearity.-The key to achieving the MSE is the nonlinearity mediated by strong coupling between cavity photon modes and excitons. As we now explain, this nonlinearity in the cavityexciton system can arise directly through the Stark effect. As a simple and clear illustration of the MSE, we first focus on a single localized excitonic mode interacting with a single cavity mode (of a single polarization). We will discuss the MSE for delocalized excitons in an extended 2D excitonic layer later in the text.
We model the localized exciton mode as a simple twolevel system with bare resonance angular frequency ν (0) ; we denote the ground state (no exciton) by |P = 0 , and the excited state (exciton present) by |P = 1 . The cavity photon mode has angular frequency ω (0) . In the dispersive limit, the dynamics of the system are described by the Hamiltonian H = H X + H 0 + H int with (setting = 1 here and throughout, unless otherwise stated): where a † is the creation operator for the cavity photon mode, andP = s z + 1/2 counts the exciton population viaP |P = P |P , where 1 is the 2 × 2 identity matrix and s z = σ z /2, where σ z is the third Pauli matrix. The last term in Eq. (1) encodes a dispersive coupling, V , between the excitons and cavity photons, that is valid for V |ν (0) − ω (0) |, ω (0) , ν (0) . In this limit, the magnitude of V can be controlled directly through engineering of the microcavity mode profile and its detuning from the exciton resonance, ∆ = ν (0) −ω (0) . Throughout this work we will consider ∆ > 0, which ensures that V > 0; for a derivation of the dispersive coupling V and parameter estimates for a TMD/cavity system, see Supplementary Information (SI).
Crucially, through the dispersive coupling, both the exciton and cavity photon resonances are mutually dependent on the other's occupation. For a state with m cavity photons present, and excitonic state P = {0, 1}, the (cavity-dressed) exciton and (exciton-dressed) cavity photon resonance angular frequencies,ν(m) andω(P ), respectively, are given by: The excitonic resonanceν(m) experiences a blue shift away from its bare resonance frequency that is proportional to the photon number in the cavity -the optical Stark effect [3,4,6,8] Illustrative dimensionless parameters used: the cavity photon resonance frequencyω(P ) also depends on the occupation of the excitonic state, shifting as P changes. The mutual tuning of exciton and cavity photon transitions exhibited in Eq.
(2) provides a natural means of feedback, and as we now discuss, gives rise to nonlinear dynamical phenomena in the system. Multistable Stark effect and cavity-exciton steady states.-To demonstrate the MSE, we consider an exciton-photon microcavity system with laser drives at angular frequencies ν d and ω d . These fields pump the excitonic and cavity photon modes, respectively. This selectivity can be achieved by choosing ν d and ω d to be slightly detuned fromν(0) andω(0), respectively, with their individual detunings much smaller than ∆.
In the presence of these laser driving fields, the Hamiltonian becomes H(t) = H + H (d) where F X and F 0 are the drive amplitudes, and σ + = (σ x + iσ y )/2, where σ x,y are the x, y Pauli matrices. In anticipation of making a rotating wave approximation below, we have discarded counter-rotating terms in Eq. (3).
To explicitly demonstrate the MSE, we track the exciton and cavity photon populations in the driven system in the presence of Markovian dissipation that accounts for exciton relaxation (recombination) and cavity photon loss. As a first step, we transform into a frame that co-rotates with the drives, using U (t) = exp (−iω d t a † a − iν d tP ). In the rotating frame, the system evolves according to the (static) Hamiltoniañ H =H X +H 0 +H int , withH X = (ν (0) −ν d )P +F X σ x /2, andH 0 = (ω (0) − ω d )a † a + F 0 (a † + a)/2. The interactioñ H int = V a † aP does not change under the transformation.
While Eq. (4) can generically encode a variety of complex dynamical regimes of the composite system, as we now discuss, a large separation in the cavity and exciton decay timescales enables direct evaluation of the MSE steady states (cf. Ref. [18] for general discussion). Indeed, the regime wherein the excitonic system relaxes far faster than the cavity photon system can be readily achieved in many exciton-cavity setups, see estimate below. Physically, this separation of timescales means that the reduced density matrix of the excitonic system ρ X (t) ≡ Tr 0ρ (t) rapidly reaches a quasistationary state over a time that is short compared with the characteristic evolution timescale of the cavity photon; here Tr 0 [Tr X ] denotes the partial trace over photonic [excitonic] degrees of freedom. On the timescale of excitonic relaxation, the cavity stateρ 0 (t) ≡ Tr Xρ (t) can be treated as quasistatic, allowing the formation of an excitonic steady state that depends parametrically onρ 0 . On the timescale that the cavity stateρ 0 (t) evolves,ρ X (t) maintains a quasistationary state that adiabatically follows the slow evolution ofρ 0 (t).
Using this separation of timescales, in describing the time evolution of the exciton and cavity photons we adopt a mean-field decoupling [18] of Eq. (4) by replacing the cavity-exciton coupling by its mean-field averages Tr 0 (ρ(t)H int ) → V m(t) P and Tr X (ρ(t)H int ) → V a † a P (t) , where P (t) ≡ Tr[Pρ X (t)] and m(t) ≡ Tr[a † aρ 0 (t)]. This mean-field decoupling is justified in the semiclassical regime where the photon number in the cavity is large and fluctuations are small [18]. With this mean-field decoupling, the (rotating frame) exciton and cavity density matricesρ X (t) andρ 0 (t), respectively, evolve according to: The exciton population dynamics can be obtained by directly evaluating the elements ofρ X (t) in Eq. (5) to obtain effective Bloch equations. Writing s i (t) ≡ Tr [s iρ X (t)] where s i = σ i /2 for i = x, y, z, and noting Tr [ρ X (t)] = 1, we obtain where δν(t) = ν d −ν[ m(t) ]. We solve for the excitonic (quasi)-steady state by setting the three equations above equal to zero, and assuming that the cavity mode occupation m(t) = m is fixed. We thus obtain the (quasi)-steady-state population of the excitonic mode as a function of the cavity occupation, m: where s z is the time independent steady state solution of Eq. (7). As evident from Eq. (8), the steady state excitonic population depends both on the excitonic drive strength, F X , and parametrically on the cavity population through the stark-shifted exciton resonance,ν(m).
The steady state cavity population can be obtained heuristically by first considering the familiar expression for the average population of a driven cavity mode with a fixed resonance frequency, ω:m Due to the Stark-induced mutual tuning described above, the cavity resonance frequency changes with the exciton population, see Eq. (2). As a result, we replace ω →ω[P (m =m)] to yield a self-consistency relation for the cavity mode population: We note that this heuristically-obtained self-consistency relation agrees with results obtained through careful analysis of the evolution of the full density matrix of the joint system [42], in the regime κ/γ 1, and V 2 /γ κ [18]. The steady state cavity occupation thus depends on the steady state exciton population through its mutually-tuned cavity transition in Eq. (2).
We now explicitly exhibit the multistability described by Eqs. (8) and (9). Choosing drive frequencies slightly blue detuned from the bare exciton and cavity resonances (see Fig. 2 and caption for parameter values), Eqs. (8) and (9) yield multiple solutions for P as a function of F 0 (for all other parameters held fixed in this regime). These multiple steady states arise from the MSE, as evidenced by the jumps of the Stark shift δE (on the order of the exciton decay rate γ) displayed in Fig. 2b.
Within the bistable regime, two distinct stable steadystate solutions for P and δE exist for the same drive parameters (solid lines). This enables a hysteretic behavior of the excitonic system that depends on the history of the optical drive. Indeed, as F 0 increases from zero (forward sweep), P (as well asm) jump to the upper branch of solutions (upward arrow) at a forward threshold amplitude. However, when F 0 is then decreased (reverse sweep), both P (and δE) jump to the lower branch of solutions (downward arrow) at a distinct reverse threshold amplitude. This hysteresis enables the system to operate as an optically-controlled "exciton switch" with "off" and "on" states as the lower and upper branches. Strikingly, this excitonic hysteresis occurs even for a single excitonic mode, in sharp contrast to other nonlinearities induced at high exciton density [33].
MSE in transition metal dichalcogenides.-Having exhibited the MSE mechanism for a single excitonic emitter, we now discuss the MSE in readily available twodimensional excitonic systems. A natural class of candidate materials are the atomically thin TMDs, which possess room temperature stable excitons and large Stark effects [8,12,13], and can be easily integrated with planar photonic crystal cavities, as in the inset of Fig. 1.
Here we will focus on zero center of mass momentum (COMM) excitons in a single valley, where excitons obey circular polarization selection rules [19][20][21][22][23][24]: by driving the TMD with circularly polarized light of fixed handedness with frequency close to the exciton resonance, only excitons in the corresponding valley will be excited.
To describe the MSE in TMDs, we consider an extended TMD layer placed on top of a photonic cavity, see Fig. 1 inset. We first note that the TMD excitonic mode at ν (0) can have a large effective degeneracy N . This degeneracy accounts for excitons at distinct exciton center of mass spatial coordinates; these degenerate exciton emitters can form plane wave superpositions that lead to delocalized excitonic modes [25][26][27]. Importantly, the modes with zero COMM interact coherently (in phase) with the same cavity photonic mode [25][26][27] (with a wavelength of a few hundred nanometers); similarly, for exciton pumping fields that have large wavelengths of order several hundred nanometers, multiple excitonic emitters can be driven in phase with each other. As such, in describing the TMD layer excitonic-cavity system, we re-placeP →P tot = jP j in the Hamiltonian Eq. (1), as well as σ +,− → s +,− tot = j σ +,− j in Eq. (3), where the sum over j runs over each of the j = 1, . . . , N degenerate excitonic emitters. Similarly,ω(P tot ) →ω = ω (0) +V P tot in Eq. (2) where P tot = 0, 1, 2, · · · are eigenvalues ofP tot . Since all the emitters interact with the same cavity photon mode,ν(m) in Eq. (2) remains unchanged.
We follow a similar procedure and use the separation of timescales as discussed above for tracking the exciton and cavity photon populations (see SI for full details). In so doing, we take a spin-coherent-state ansatz so that the dynamics of the multiple emitter system can be analyzed in terms of the dynamics of a giant spin s tot = s x totx + s y totŷ + s z totẑ , where s x,y,z tot = i s x,y,z i is summed over the excitonic emitters. For fixed cavity occupation m we obtain the steady-state exciton population in the extended system (see SI) as where Γ is the exciton recombination rate (for the zero COMM excitons) in the extended TMD system; we note, parenthetically, that this rate can be estimated from the recombination rate γ of a single localized exciton emitter as Γ ∼ N γ [25][26][27]. In obtaining Eq. (10) we have taken a large degeneracy N 1 as well as focused on the lowexcitation regime.
Before we exhibit the MSE in TMD systems, we first discuss the parameters for the cavity-exciton system. We note that the excitonic mode degeneracy N can be large and can range from N ∼ 10 2 − 10 4 [25]; this arises from the large number of excitonic modes that can interact coherently within a single wavelength of either the cavity photon mode or the exciton drive [25][26][27]. An estimate of N can be obtained from the ratio of the mode area of the photonic mode (the square of its wavelength) and the effective size of an exciton (the square of its Bohr radius) [27]. Further, recombination times for zero COMM excitons in typical monolayer TMDs can range from Γ −1 ∼ 0.5 to a few picoseconds [28][29][30][31][32], whereas cavity relaxation times can be as long as tens to a hundred picoseconds [14][15][16][17]. As a result, κ Γ, justifying the separation of time scales and meanfield decoupling approach we have used to describe the MSE. Lastly, strong light-matter interaction in monolayer TMDs [8,12] can lead to sizeable values of dispersive coupling V ≈ 0.1 − 0.5 meV, see SI for a detailed estimate. In the plots we have chosen V = 0.2 meV for illustration.
Solving Eq. (10) together with Eq. (9) yields an excitonic multistability and MSE, as shown in Fig. 3. With realistic parameters for monolayer WS 2 and photonic crystal cavities, discontinuous jumps in the excitonic Stark shift can be readily achieved by moderate cavity and exciton drive intensities of order kW/cm 2 .
Interestingly, distinct regimes of multistability can be accessed; at low exciton drive strength a bistable MSE manifests (as cavity drive is swept) whereas larger exciton drives display tristabilities (see Fig. 3a,b). Indeed, the MSE displays hysteretic behavior as either exciton or cavity drives are swept, with Fig. 3b,c displaying sizeable discontinuous δE of order meV. We note that together with multistable δE shown in Fig. 3, the exciton population P tot similarly exhibits multistability and hysteresis (see also Fig. 2). While we have focused on the MSE and its concomitant excitonic multistability, multiple stable states of the cavity mode (so-called "optical multistability," characterized by distinct steady state values ofm) can also arise via the MSE. (Note that in Fig. 2b, δE is directly proportional to the cavity photon occupation.) This effect is similar to dispersive optical multistability in highly nonlinear optical media [35][36][37][38][39][40], and may provide new means for controlling optical states.
From a fundamental perspective, the MSE arises from the fact that the excitonic Stark effect is an inescapable consequence of the partial fermionic nature of excitons [6] -a property that is present even in dilute exciton systems. Indeed, we find in Fig. 3 that a MSE manifests for a steady state excitonic population on the order of P tot ∼ 1, see SI. This indicates that the MSE occurs even as approximately one exciton is excited in the entire photonic cavity (corresponding to a low exciton density of order 10 10 cm −2 ). This distinguishes MSE and its associated nonlinear phenomena from other types of multistable behavior, e.g., optical bistabilities that originate from exciton-exciton interactions that typically require a large density of excitons to enable bistable behavior [33,34]. Perhaps most exciting is how MSE-induced hysteresis in the exciton population as a function of optical drive yields jumps in P tot of order unity; these may provide controllable means of selectively exciting/deexciting a single exciton as well as controlling its emission.
J.C.W.S. acknowledges support from the National Research Foundation (NRF), Singapore under its NRF fellowship programme award number NRF-NRFF2016-05, the Ministry of Education, Singapore under its MOE AcRF Tier 3 Award MOE2018-T3-1-002, and a Nanyang For the convenience of the reader, in this section we review the optical Stark effect and how it arises from the lightmatter interaction between an exciton and light. In so doing, we concentrate on a single localized exciton emitter and its interaction with a single cavity mode. As discussed in the main text, generalization to multiple emitters is straightforward. For illustration, we begin with a Jaynes-Cummings model where σ x,y,z are Pauli matrices, σ ± = (σ x ± iσ y )/2, and g = ME 0 is the coupling constant that captures the dipole interaction between the exciton mode (described by its dipole moment M) and the cavity mode (described by the amplitude of the electric field E 0 ). Here ν bare and ω bare are the bare frequencies of the exciton and cavity mode resonances in the absence of any interaction, i.e., g = 0. The interaction -i.e., the third term of Eq. (S1) -is responsible for the optical Stark effect. To see this we perform a canonical transform T † H JC T with T = exp S; here S = −g(aσ + − a † σ − )/(2∆). The transformed Hamiltonian is given by We note that [ν bare σ z /2 + ω bare a † a, S] = −g(σ − a † + σ + a)/2. By taking the dispersive limit g ∆, the transformed Hamiltonian up to the first order in g/∆ is where V = g 2 /(2∆) is the dispersive coupling constant andP = σ z /2 + 1/2. We note that the last term, −ν bare /2, is a constant offset to the Hamiltonian that does not play any role in the physics we discuss. As a result, in the main text as well as what follows, we drop mention of the constant offset. Using this transformation, Eq. (1) can be readily read off Eq. (S3).

Stability analysis of the steady state solutions
In the plots in the main text, we have shown the multiple branches of steady state solutions that appear at a given set of drive parameters. In these plots, solid lines denote the stable solutions while dashed lines denote unstable solutions. The stability of each steady state solution can be ascertained by examining the time-evolution of the exciton occupation factor ∂ t P (t) = ∂ t s z (t) in Eq. (7). In particular, taking a small deviation δP from the steady state solution P and computing ∂ t P (t) | P +δP , we find a solution is stable if sgn ∂ t P (t) | P+δP = −sgn(δP); i.e., the solution is stable if, after the system is slightly displaced from the steady state, it time evolves back to the original steady state solution P (t) = P .

Steady state exciton population with multiple emitters
In this section, we describe the steady state solution for the MSE with multiple emitters. As in the main text, we replaceP →P tot = iP i = i (s z i + 1/2) in the Hamiltonian, Eq. (1), as well as σ +,− → s +,− tot = i σ +,− i in Eq. (3), where the sum runs over the degenerate excitonic emitters. Similarly,ω(P ) →ω = ω (0) + V P tot in Eq. (2) where P tot = 0, 1, 2, · · · are eigenvalues ofP tot . Since all the emitters interact with the same cavity photon mode, ν(m) in Eq. (2) remains unchanged. For simplicity of notation, we have defined the spin operators s = σ/2, where σ = (σ x , σ y , σ z ) is the vector of Pauli matrices. Similar to our approach described for a single excitonic emitter in the main text, we exploit a separation of time scales between cavity and exciton systems that is easily achieved in readily available systems (see discussion in main text). In this fashion, we can mean-field decouple the excitonic and cavity degrees of freedom. As a result the excitonic system evolves according to where δν(t) = ν d −ν[ m(t) ], and the dissipator for the multiple emitters reads As a sanity check, noting that [ρ X,tot , 1] = 0, we find that probability is conserved: ∂ t Tr(ρ X,tot ) = 0.
We can use Eq. (S4) to obtain the equation of motions governing the "Bloch" (giant spin) dynamics of the multiple emitters. Here we treat the collection of multiple emitters as a giant spin with magnitude S = N /2. When the spin is pointing down such that s z tot = −S, no excitons are excited. When s z tot = S, all the excitonic modes are excited. In this work we focus on the regime far from saturation, where the exciton density remains low.
Writing s i tot (t) = Tr[s i totρX,tot (t)] for i = x, y, z, we find where ijk is the Levi-Civita symbol.

Spin coherent state ansatz
As we now discuss, the dynamics of the multiple emitters can be understood as the dynamics of a giant spin s tot = s x totx + s y totŷ + s z totẑ with a spin magnitude S. Here S = N /2 can be estimated from the degeneracy of the multiple emitters in the system. Take for example the limiting case discussed in detail in the main text: a single two level system (i.e., spin S = 1/2 system). For spin-1/2, we recall that {s z , s y } = {s x , s y } = {s z , s x } = 0; similarly, s 2 − (s z ) 2 = 3/4 − 1/4 = 1/2. Using these identities, we find Eq. (S6) reduces to Eq.
with the giant spin pointing in the directionn: In using this coherent state, it is useful to note that expectation values of various operators (e.g., s z tot s x tot ) in this coherent state can be obtained by applying a suitable rotation to the operator and evaluating the expectation value in a well-known reference state (e.g., |ψẑ ). For example, where U is unitary operator that rotates the spin (determined byn). Using this coherent state ansatz, we obtain where s ⊥ = ( s x tot 2 + s y tot 2 ) 1/2 . In obtaining Eq. (S10) we have used the identity ψẑ|s x tot s z tot |ψẑ = ψẑ|s y tot s z tot |ψẑ = 0. Further we have noted that ψẑ|s x tot s y tot |ψẑ = − ψẑ|s y tot s x tot |ψẑ . This latter identity can be discerned from noting that in s x tot s y tot = ij s x i s y j only the terms where i = j yield non-zero values, and that when i = j, we have {s x , s y } = 0. Substituting Eq. (S10) into Eq. (S6) we get In the absence of a drive, F X = 0, and at steady state, the excitonic system yields its equilibrium value with s z tot eq = −S (similarly, s x tot eq = s y tot eq = 0). Interestingly, Eq. (S11) is a non-linear equation in s z tot . Such nonlinearities become important only when the deviation of s z tot away from its equilibrium (no excitation) value is large (i.e., comparable to s z tot eq ), due to saturation of the exciton resonance. We now analyze the steady state population of the excitons by first takingν[ m(t) ] →ν(m) for a fixed m in the same way as discussed in the main text. In the low-excitation regime where s z tot = −S + P tot for P tot S, we can safely linearize Eq. (S11) in P tot . Solving for P tot at steady state in Eq. (S11), we obtain where Γ ∼ 2Sγ is the recombination rate for excitons in the entire multiple emitter system. In obtaining Eq. (S12), we have used that s 2 ⊥ /S P tot . Writing S = N /2 produces Eq. (10) in the main text.
Alternative derivation of exciton population in the low-density limit: Holstein-Primakoff transformation In this section, we provide an alternative derivation of the steady state exciton population (in the case of multiple degenerate emitters) in the low-density limit. In so doing, we note that the dynamics of the exciton population can be understood from analyzing the evolution of a large spin with magnitude S 1. In the large S limit, the excitations of this system can be analyzed in a bosonic framework, using a Holstein-Primakoff transformation. To begin, we write where b † is a bosonic creation operator so that [b, b † ] = 1. Using Eq. (S13) we can readily verify Note that Eq. (S14) is exact; the factors of [1 − b † b 2S ] 1/2 must be included in order to obtain the b † b dependence on the right-hand side (which will prove to be essential below).
Noting that b † b counts the number of excitons created [see Eq. (S14)] we express the state where P tot excitons have been excited as where the vacuum state (no excitons present) corresponds to the giant spin pointing downwards ( s z tot = −S). In the co-rotating frame of the drive, we find the Hamiltonian that describes the excitons is where we have recalled thatP tot = s z tot + i 1/2 ≈ s z tot +S and s x tot = (s + tot +s − tot )/2. Here we have estimated N = 2S; the arrow labeled HP denotes the Holstein-Primakoff transformation. In obtaining the last term of Eq. (S16) we have Similarly, we find the exciton-cavity interaction reads as (S17) This expression captures the nonlinear interaction between the exciton degrees of freedom (characterized by b's) and the cavity photons (characterized by a's). Note that this interaction is distinct from the bilinear exciton-cavity interaction typically discussed for the polariton effect. The nonlinearity arises due to the underlying fermionic nature of the electrons and holes that make up the excitons [6]. The role of the fermionic Pauli exclusion is reflected in the square root factors in the relation between the collective exciton creation operator s + tot and the bosonic creation operator b † (and similarly for s − tot and b), as well as the b † b dependence of [s + tot , s − tot ] as shown in Eq. (S14). In a similar fashion to that described in the main text and above, we exploit a separation of time scales between the exciton relaxation and cavity relaxation time scales and write the equation of motion of the excitonic excitations as where I HP , with Γ capturing a phenomenological decay rate of the excitonic mode. The steady state solution to Eq. (S18) can be readily obtained using a coherent state ansatzρ tot,HP X = |β β| where b|β = β|β and β ∈ C. By direct calculation we verify that the coherent state ansatz indeed yields a steady state solution of Eq. (6) for fixed oscillator excitation m, with β = (−iF X S/2)/(Γ − i(ν d − iν[m])). The corresponding steady state exciton population is given by As a consistency check, we note that in the limit F X Γ, ν d ,ν, Eq. (S12) reduces to Eq. (S19).

Estimate of parameter values for MSE in a TMD photonic cavity
In this section, we estimate the values of the parameters used in the main text to achieve the MSE in a transition metal dichacolgenide (TMD) sample placed in a photonic crystal cavity. For the purposes of the estimates below, we will restore so that ω, ν as well as κ, Γ are energies.
We consider a cavity with a high quality factor Q = ω (0) /κ and compressed effective mode volume V mode = α(λ/n) 3 , where κ is the cavity mode linewidth, ω (0) is the bare cavity resonance frequency, λ is the free space wavelength of the photons, n is the refractive index of the cavity, and α is a proportionality constant that depends on the geometry of the cavity. In photonic crystal cavities, α values in the range 0.02 ∼ 1 have been achieved with quality factors of order 10 3 ∼ 10 6 [14][15][16][17]. Here for the purposes of a simple demonstration, we choose α = 0.05 and Q = 20000. Similarly, we will consider an exciton resonance at ν (0) = 2 eV and set the cavity mode resonance frequency, ω (0) , to be red-detuned away from it (as discussed below, we choose a detuning of 43 meV). These parameters correspond to cavity linewidth of κ ≈ 0.1 meV; taking the refractive index in the cavity as n ≈ 3, we obtain an effective mode volume of V mode ≈ 4.7 × 10 5 nm 3 .

TMD Stark shift and exciton-cavity interaction
Here we estimate the value for the dispersive coupling V for a TMD material placed on the photonic crystal cavity described in the main text. First, we note that excitons in monolayer TMDs obey circularly polarized optical selection rules wherein the excitons around K (K ) valleys primarily interact with photons of left (right) circular polarization. This valley-circularly polarized light selectivity justifies our use of single flavor of exciton modes (e.g., excitons in the K valley) interacting with a single cavity mode (e.g., left hand circularly polarized cavity mode) in the main text. Indeed, this selectivity is well evidenced in experiment. It not only manifests in selective excitation of excitons in the valleys (by using circularly polarized light), but also yields a valley Stark effect wherein light of left (right) circular polarization only shifts the exciton resonances of K (K ) excitons [8,12,13]. For a left-hand circularly polarized field E(t) = E 0 (x cos ωt +ŷ sin ωt), excitons in the K valley experience a Stark shift [8,12,13] where ∆ = (ν (0) − ω) and M is a material dependent dipole matrix element for the circularly polarized field. In Ref. [8], an optical Stark shift in a single valley of monolayer WS 2 was measured. Ref. [8] reported that a peak Stark shift value δE expt = 18 meV was induced by a circularly polarized laser pulse with fluence of 120 µJ/cm 2 and pulse width 160 fs. The detuning between the exciton and laser pulse was ∆ expt = 180 meV. Recalling that the intensity of a circularly polarized field E(t) (see above) is I = 0 c|E expt 0 | 2 , we estimate peak E expt 0 = 5.3 × 10 7 V/m in Ref. [8], where 0 is the vacuum permittivity, and c is the speed of light. Here we have estimated peak intensity as fluence/pulse width.
The value of V (Stark shift per circularly polarized photon in the cavity) for our TMD in a photonic cavity setup can be estimated by comparing with the above experiments on the Stark effect in WS 2 [8]. To do so, we first note that the electric field amplitude of the circularly polarized cavity mode can be estimated as where eff c is the effective relative permittivity of the cavity. Taking V mode from above, and eff cav = 10, we estimate an electric field strength of the circularly polarized cavity mode as E (0) cav = 2.7 × 10 6 V/m. Comparing the cavity mode electric field amplitude with that used in the experiment above [8] and specifying an exciton-cavity detuning ∆ cav = 43 meV, we obtain V (Stark shift per circularly polarized photon in the cavity) as where we have used the proportionality relation in latter part of Eq. (S20). While we have chosen a specific value of ∆ cav which yields V ≈ 0.2 meV, a range of values of V ∼ 0.1 − 0.5 meV can be readily achieved by tuning ∆ cav . We note that V in the cavity system can also be estimated directly from the dipole matrix element for the circularly polarized field acting on a TMD (in our case, e.g., WS 2 ). Using Eq. (S20) and the experimental parameters from Ref. [8] discussed above, we find M ≈ 73.6 Debye. We note that this experimentally inferred M is within a factor of unity from that obtained from a simple theoretical estimate [8] M theoretical ≈ 56 Debye. For the exciton-cavity photon system, we recall that the exciton-cavity photon interaction in Eq. (S1) is g = ME (0) cav ≈ 4.15 meV. Using V = g 2 /2∆ cav from Eq. (S3) we obtain V ≈ 0.2 meV in agreement with Eq. (S21).

Exciton drive strength, FX
In this section, we connect the exciton driving field intensity (in kW/cm 2 ) to the F X driving strength (in meV) used in the main text. The strength of the exciton drive F X in Eq. (3) of the main text can be obtained in the same fashion as that described above via F X = ME d X , where E d X is the amplitude of a circularly polarized driving electric field E drive . The intensity of the (exciton) driving field is given by where eff is an effective relative permittivity that depends on the geometry of the incident irradiation. For illustration, we have taken eff = 1. For a driving intensity of I X = 0.1 − 0.6 kW/cm 2 , we estimate a circularly polarized driving electric field strength as E d X = 0.1−4.8 × 10 4 V/m. Using the value of M above we obtain F X ≈ 0.002−0.072 meV. We note parenthetically that these values of drive electric field are far smaller (two orders of magnitude smaller) than the electric field amplitude of the cavity mode discussed above.

Cavity drive strength, F0
In this subsection, we relate the incident cavity drive power P to the cavity drive parameter F 0 used in the main text. To do so we first briefly review classical coupled mode theory [41]. In a cavity with cavity eigenmodes e i (r), the electric field profile in the cavity can be described by E cavity (r, t) = Re i C i (t)e i (r) ; here i indexes the cavity eigenmodes. C i (t) is a time-varying amplitude that describes the degree to which each of the cavity eigenmodes are occupied over time. According to coupled mode theory, the classical cavity mode fields coupled to an external drive evolve according to [41] where f + i (t) is the amplitude of the incident drive field that couples to the cavity mode i. Here ω i is the resonance frequency of cavity mode i, κ is the decay rate of the cavity mode, and ξ = √ 2κ represents the strength of the cavity-incident driving field coupling [41].
We now focus on a planar photonic cavity where there are two degenerate modes in the plane, which we label e x and e y , with common frequency ω x = ω y = ω. Here e x and e y are cavity modes that are polarized in the x and y directions. We will also assume that the system is isotropic in-plane so that κ x = κ y = κ. Further, we focus on cavities that primarily decay radiatively, consistent with a high quality factor cavity. We remark that, following the standard convention [41], the cavity modes profiles e x,y (r) can be normalized in such a way that where U (t) is the energy in the cavity and P is the incident driving power. To describe the driving of the cavity by a circularly polarized beam we will use f + x (t) = f 0 cos ω d t and f + y (t) = f 0 sin ω d t. We now turn to the quantum mechanical description of the driven cavity: x a x + a † y a y ) + where a x,y are annihilation operators for the x, y polarization cavity photon modes (in modes e x,y ), and the last two terms describe a circularly polarized driving with amplitude F 0 . The latter terms describe the full Rabi-type drive-cavity coupling, see below for discussion. Here we have dropped the vacuum energy of the cavity mode as it acts as a constant energy off-set; including it does not alter our conclusions/estimate for F 0 below. In the same fashion as described in the main text, the evolution of the cavity can be described by the master equation where I(ρ) = κ a † x a x ρ − 2a x ρa † x + ρa † x a x + (x → y) tracks the decay of the mode. In order to track the evolution of the amplitudes above, we define the mode amplitude operators C i = C (0) i a i , for i = x, y). Following Eq. (S23) above, we have U = i C † i C i = |C (0) x | 2 a † x a x + |C (0) y | 2 a † y a y . This energy must correspond to the energy of the cavity obtained directly from the first term of Eq. (S24): ω a † x a x + a † y a y . As a result, we obtain C (0) = |C (0) ω. The temporal evolution of the cavity amplitude C i (t) = Tr C i ρ cav (t) can be described by multiplying C i into Eq. (S25) and taking the trace. Recalling the commutator identities [a † a, a] = −a, [a, a † ] = 1, and cyclically permuting the operators, we obtain Relating C i (t) with C i (t) in Eq. (S22), we see that Eq. (S26) directly corresponds to Eq. (S22). This yields the association F 0 C (0) / √ 2 = f 0 ξ. As a result, we find that the value of F 0 can be directly estimated from the input driving power where in the last line we have estimated the incident driving power from the laser intensity P = IL 2 where L 2 is the incident area of the cavity; here I is the laser intensity. The incident area can be estimated from the mode volume V mode = L 3 . For the parameters for mode volume chosen, we arrive at L ≈ 78 nm. We remark that Eq. (S24) is the full Rabi Hamiltonian for a cavity mode driven by an external circularly polarized field. We now re-write Eq. (S24) into the basis of circularly polarized operators a L = (a x − ia y )/ √ 2 and a R = (a x + ia y )/ √ 2: This formulation allow us to focus on the component with a single circular polarization by discarding the counterrotating component in the rotating wave approximation. We observe that in the Heisenberg picture, ∂ t a L,R = i/ [H cav , a L,R ], the operator a L,R evolves as e −iωt . Similarly, the operator a † L,R evolves as e iωt . Thus, the right-handed circularly polarized component a † R e iω d t + a R e −iω d t rotates about twice as fast as ω when moving to the rotating frame (as employed in the main text). Therefore, within the rotating wave approximation, we can discard the right-handed component and the resulting Hamiltonian is consistent with Eq. (3) in the main text.
Multistable P tot in a TMD photonic cavity In this section, we plot the multi-stable steady state excitonic population as a function of exciton and cavity drives in Fig. S1. This displays how the average steady state P tot takes on small values on the order of P tot ∼ 1. This indicates that the MSE occurs even as approximately one exciton is excited in the entire photonic cavity (corresponding to a low exciton density of order 10 10 cm −2 ). Further, we find that jumps in the exciton population (see Fig. S1b) can be tuned to be of order unity.