Phase Diagram and Self-Organizing Dynamics in a Thermal Ensemble of Strongly Interacting Rydberg Atoms

Far-from-equilibrium dynamics that lead to self-organization are highly relevant to complex dynamical systems not only in physics but also in life, earth, and social sciences. However, it is challenging to find systems with sufficiently controllable parameters that allow quantitatively modeling of emergent properties. Here, we study a nonequilibrium phase transition and observe signatures of self-organized criticality in a dilute thermal vapor of atoms optically excited to strongly interacting Rydberg states. Electromagnetically induced transparency provides excellent control over the population dynamics and enables high-resolution probing of the driven-dissipative dynamics, which also exhibits phase bistability. Increased sensitivity compared to previous work allows us to reconstruct the complete phase diagram, including in the vicinity of the critical point. We observe that interaction-induced energy shifts and enhanced decay only occur in one of the phases above a critical Rydberg population. This case limits the application of generic mean-field models; however, a modified, threshold-dependent approach is in qualitative agreement with experimental data. Near threshold, we observe self-organized dynamics in the form of population jumps that return the density to a critical value.


I. INTRODUCTION
Self-organization and nonequilibrium dynamics in complex dynamical systems are responsible for a diverse range of phenomena not only in physics but also other fields such as earth sciences, biology, and economics [1]. Many nonequilibrium systems exhibit self-organized criticality (SOC) [2] as they evolve to an attractor that coincides with a critical point in their phase diagram. SOC is considered a source of complexity in nature; it produces scale invariance and "pink" 1=f-frequency noise, and it makes systems insensitive to parameter fluctuations. Quantitative modeling of simple experimental systems displaying nonequilibrium dynamics and dynamical phase transitions is a key to enhancing our understanding of nonequilibrium phenomena [3]. However, finding systems with both strong and controllable dynamical nonlinearities is challenging. Although laser excitation in dilute atomic gases allows precision measurements of state populations and phase diagrams, typical interatomic interactions are too weak to affect nonlinear excitation dynamics. For example, the observation of optical bistability in nonequilibrium light-matter systems [4][5][6][7] initially required cavity feedback or decoupling from the environment in cryogenic setups [7]. Much stronger interactions can be achieved in ensembles of Rydberg atoms [8,9] with strong dipolar interactions and high sensitivity to charges originating from ionized Rydberg atoms. Thanks to these mechanisms, Rydberg gases exhibit rich nonequilibrium dynamics and can form exotic phases of matter [10][11][12][13]. Ultracold atom experiments have observed aggregate formation [14,15] around seed excitations [16][17][18][19][20] or ions [21], which facilitate Rydberg excitation, bistable or metastable thermodynamical phases [14,22], or SOC [23]. Nonequilibrium phase transitions [24][25][26][27][28] and aggregate formation [29] are also observed in experiments using room-temperature Rydberg vapors. In this paper, we investigate the nonequilibrium dynamics of a driven-dissipative ensemble of strongly interacting Rydberg atoms in a room-temperature atomic vapor. We show that the resulting excitation dynamics can be simulated using an analogue of the "forest-fire" model [30], a prominent example of SOC (Sec. II).
Below, we introduce the nonequilibrium system, which features two thermodynamical phases: one with low Rydberg atom density and thus no interaction (NI phase), and another with high Rydberg density, in which strong interactions facilitate excitation (I phase). In our experiments (setup in Sec. II), we exploit the narrow linewidths achievable using electromagnetically induced transparency (EIT) [31,32] to probe the system with MHz-scale frequency resolution-an enhancement of 2 orders of magnitude compared to previous work [24][25][26]. This method allows us to probe the dynamical phase diagram with unprecedented precision and observe previously unobserved spectral features in the vicinity of the critical point (Secs. IVA and IV B). We show that a combination of interaction-induced line shifts and broadening are responsible for this rich structure. The results are suggestive of the importance of ionizing collisions in the bistable dynamics as identified in other recent work [26,28]. We also show that the transition threshold can be manipulated by adding a second probe beam that enhances the Rydberg excitation rate (Sec. IVA).
Finally, we observe behavior indicative of SOC near the phase transition. The system is highly susceptible to small variations, and even short-term fluctuations in the Rydberg density can induce a phase transition, giving rise to self-organization and pink noise (Sec. IV D). We also observe power-law scaling for the size of aggregate clusters (Sec. IV E), which is consistent with our forest-firetype model.

II. BACKGROUND AND MODEL
We work with the following nonequilibrium system [ Fig. 1(a)]: A probe beam propagates through a thermal vapor of three-level atoms and drives the transition between a ground state jgi and a low-lying, short-lived excited state jei. A coupling field couples jei to a highly excited and long-lived Rydberg state jri. The probe (coupling) Rabi frequency, angular frequency, and detuning are denoted by Ω pðcÞ , ω pðcÞ , and Δ pðcÞ , respectively, while Γ eðrÞ are the decay rates of jei (jri). Since Γ r ≪ Γ e , we observe ladder EIT [32], where for Ω p ≪ Ω c , the ensemble is rendered transparent to the probe in a narrow frequency window around two-photon resonance Δ p ¼ −Δ c . In the experiment, we observe the probe transmission on resonance The interaction between Rydberg atoms introduces a strong dynamical nonlinearity as required for selforganization. Moreover, SOC relies on interaction-induced avalanches in the nonequilibrium dynamics [33] to trigger a dynamical phase transition. In a thermal Rydberg vapor, such an avalanche can be induced by either dipolar interactions [29] or ionizing collisions with electrons, ions, or other atoms [34][35][36][37]. In either case, the avalanche occurs above a critical Rydberg density N R;ðcÞ , and we observe two distinct thermodynamic phases: a noninteracting phase (with Rydberg density N R < N R;ðcÞ ) and a strongly interacting phase (N R > N R;ðcÞ ). In the latter, interactions cause an effective, N R -dependent detuning of the coupling field, . The probe beam is overlapped with a counterpropagating EIT coupling beam. Their transmission signals are detected on a differencing photodetector (DD). The probe and reference fields drive the jgi ¼ j5S 1=2 ; F ¼ 3i → jei ¼ j5P 1=2 ; F 0 ¼ 2i, and the coupling light couples the jei → jri ¼ jnD 3=2 i transition in 85 Rb. A switch beam can be added to enhance the Rydberg population without changing the probe intensity. PBS is for polarizing beam splitter, DM is for dichroic mirror, and λ=2 stands for half-wave plate. (b) Excitation and interaction dynamics in the different phases. In the noninteracting (NI) phase, the probe transmission under EIT conditions remains unaffected. Above a critical Rydberg atom density N R;ðcÞ , the system is found in an interacting (I) phase, where strong interactions lead to energy shifts and enhanced decay of jri. Fig. 1(b)]. The jnD J i-Rydberg states used in this work are also subject to broadening and enhanced dephasing, Γ r → Γ r þ Γ 0 ðN R Þ [38,39], resulting from the motional averaging over dipolar interaction potentials or position-dependent Stark shifts and atom loss due to ionization [37,40] (Appendixes A and B). These result in facilitated Rydberg excitation in the I phase and changes in macroscopic properties of the system. Here, we observe the changes of the optical response, i.e., discrete changes in the probe transmission, as the optical susceptibilities of the EIT systems are strongly affected by the level shifts and broadening in the I phase.
The Rydberg density N R depends not only on the current parameters of the driving fields but also on the phase of the system determined by previous excitations. This case effectively introduces an element of system memory. Hence, the effective values of Δ c and Γ r depend on the scan direction and give rise to hysteresis effects that result in a bistable optical response. Near N R;ðcÞ , even small fluctuations can disturb the equilibrium and trigger a dynamical phase transition due to the avalanche effect [1].

A. Self-organization
To model self-organization in our Rydberg system, we adapt a model first employed to describe the spreading of forest fires based on cellular automata that exhibits self-organized criticality [30]. We consider a 3D array of 100 × 100 × 100 cells. First, the cells are randomly filled with N R Rydberg and N G ground-state atoms. A cell is in the NI phase if N R < N R;c [green in Fig. 2(a)], in the I phase for N R > N R;c (red), or depleted of Rydberg excitations (black). These situations correspond to a healthy tree, a burning tree, or an empty site in the original forestfire model. The thermodynamic phase of each cell is then iterated according to the rules illustrated in Fig. 2(a): Cells in the I phase facilitate Rydberg excitation and trigger an excitation avalanche in adjacent cells in the NI phase such that they undergo a transition to the I phase. This case corresponds to the spreading of fire to adjacent trees; see Figs. 2(a1)-2(a3). Following the excitation avalanche, the Rydberg population in a cell is depleted to zero. Nonfacilitated Rydberg excitation induces transitions of depleted cells to the NI phase (0 < N R < N R;c ) with probability p, corresponding to the growth of a new tree on an empty site and of cells in the NI to the I phase with probability f even if not adjacent to cells in the I phase as N R , corresponding to trees catching fire due to lightning. In practice, the probabilities p and f that a cell undergoes a phase transition depend on probe and control detunings and Rabi frequencies and the phase in neighboring cells [13]. The closed boundary of our array reflects the finite excitation volume in the experiment as defined by the laser beam geometry. Figure 2 phases in a single plane following 30 iterations. Figure 2(c) depicts the result of our model. We record the number of cells n t in each isolated cluster of cells in the I phase following each iteration. The cumulative distribution of cluster sizes following all iterations has a power-law behavior with exponent b. The shape of the distribution is scale invariant, that is, where C is a constant that characterizes the scaling transformation. The scaling invariance implies that the cells organize themselves into a critical state where events of all sizes occur. We obtain jbj from 1.7-4 for p ¼ 0.01, 0.05, and 0.1. For larger p, the occurrence of large cluster sizes increases. We also vary f and record jbj against p [see

B. Rydberg EIT in a thermal vapor
To model the effect of the dynamical phase transition on the probe transmission, we compute the evolution of the ensembles' density matrix ρ by adopting a thresholdmodified mean-field approach where the additional detuning Δ 0 ¼ η 1 N R and the enhanced decay Γ 0 ¼ η 2 N R are proportional to the Rydberg density above a critical Rydberg population ρ ðcÞ rr corresponding to N R;ðcÞ . The full master equation and Hamiltonian are given in Appendix C. This approach is applicable independent of whether interactions originate from dipolar interactions or ion-induced Stark shifts. The complex susceptibility of the EIT medium, including the Doppler effect due to atomic motion, is Þv=c denotes Doppler shift experienced by an atom moving with velocity v and δ ¼ Δ c þ Δ p the two-photon detuning. Note that Ω eff is the effective Rabi frequency of coupling laser. The transmission of the probe beam through the EIT medium can be obtained from the susceptibility via where L is the medium length and k the wave vector of the probe field. To model the bistability in the EIT spectra, we compute χðvÞdv, substituting Δ c → Δ c þ Δ 0 ðN R Þ and Γ r → Γ r þ Γ 0 ðN R Þ once the critical threshold population ρ rr > ρ ðcÞ rr is reached.

III. EXPERIMENTAL SETUP
The experimental setup and a scheme of the relevant 85 Rb energy levels are depicted in Fig. 1. The probe and coupling fields counterpropagate through a 5-cm-long Rb vapor cell at a temperature of T ≈ 50°C. The atomic density is 1.5 × 10 11 cm −3 , corresponding to a mean interatomic spacing of approximately 1 μm. The probe beam (1=e 2 -waist radius of approximately 500 μm) couples The estimated peak values of Ω c =2π are approximately 20 MHz for n ¼ 47 and approximately 10 MHz for n ¼ 70 unless stated otherwise. A reference beam identical to the probe, but not overlapping with the coupling light, is sent through the cell in parallel, and its transmission signal is subtracted from the probe signal following detection on a pair of balanced amplified photodiodes.
The nonequilibrium phase transition is investigated as three parameters are varied: the coupling detuning Δ c and the Rabi frequencies Ω p and Ω c . Scanning Δ c allows us to obtain EIT transmission spectra, while the probe light is kept on resonance, Δ p =2π ¼ 0 MHz, such that the Doppler-broadened absorption background can be neglected. Scanning Ω p alters the Rydberg population, and thus mean interaction strength, as N R ∝ ρ rr ∝ Ω 2 p (see above). Here, we make an approximation that all atoms undergo the same Ω p and Ω c along probe propagation by ignoring its attenuation.
In addition to the probe, an additional switching beam can be applied to alter the Rydberg density and control the threshold of the phase transition without changing the properties of the probe itself. It also couples jgi and jei with detuning Δ p but an independent Rabi frequency Ω s . It intersects with probe and coupling beams at the center of the cell at 2°to avoid crosstalk between the switching, probe, and reference beams at the detector.

A. Optical bistability
In order to observe Rydberg-mediated optical bistability in the optical response of the three-level EIT medium, we initially observe the EIT line shape as Δ c is scanned through resonance from negative to positive detuning and vice versa. , identical EIT spectra are observed for either scan direction without bistability. At 6.7 MHzand, consequently, only slightly higher Rydberg density N R -we observe almost symmetric optical bistability above and below resonance [ Fig. 3(a2)]. A sudden drop in the transmission occurs on resonance, where the excitation rate is highest, independent of the scan direction indicating the transition from the NI to the I phase. In this regime, the energy shift Δ 0 induced on jri is small compared to both Δ c and the EIT linewidth. This behavior at low Ω p is distinctively different compared to previous experiments on Rydberg-induced optical bistability that used shelving techniques with higher Ω p and linewidths [24,26]. These observed bistability at large red detuning as the peak transmission was strongly shifted due to higher ρ rr , and transmission in the I phase (scanning from positive to negative Δ c ) exceeded the NI phase. The ability to probe this new regime and detection of the phase transition with sub-MHz frequency resolution is a consequence of the narrow EIT resonance. The bistability windows observed in Fig. 3(a2) have widths <0.5MHz, 2 orders of magnitude narrower compared to previous experiments [24,26].
, the spectra become asymmetric as an increasingly strong spectral shift Δ Ã [see definition in Fig. 3(a5)] of the peak transmission relative to Fig. 3(a3)], bistability is still observable within the EIT resonance feature, but the phase transition occurs at lower or higher Δ c value for positive or negative scan directions, respectively. The onresonance phase transition disappears, as the system is in the I phase for either scan direction.
For Δ 0 comparable to the EIT linewidth, near-resonant bistability features disappear entirely. Instead, new bistability windows appear on both sides of the EIT window, as in shelving experiments [24,26]. Unlike above, the transmission in the I phase is now exceeding that in the NI phase as the former's enhanced decay leads to a broadened transmission window. In the bistability window below Δ c ¼ 0, the system is found in the I phase for scans from positive to negative Δ c , where the scan has already crossed through resonance and sufficient Rydberg population has been built up to maintain the state. Above Δ c ¼ 0, the roles are reversed and system is in the I phase if scanned from negative to positive detuning, again after crossing through resonance. The bistability window is narrower for Δ c > 0 as the sign of Δ 0 implies an overall redshift of the spectra due to the choice of jri. Bistability for blue detuning was also observed by Weller et al. [24,26] but not in the original experiment by Carr et al. [24] where the probe was significantly stronger. Figure 3(b) shows the calculated EIT spectra with Doppler averaging as given by Eq. (4). The critical density is estimated to be N r;ðcÞ ¼ 2.9 × 10 10 cm −3 (assuming that each probe photon creates a Rydberg atom). In the model, we have set η 1 =2π ¼ 1.27 × 10 −2 MHz μm 3 and η 2 =2π ¼ 7.96 × 10 4 MHz μm 3 , to match the calculated EIT spectra to the experimental results. As dissipation is increased in the I phase, we reduce the threshold population κρ ðcÞ rr below which the system reverts to the NI phase. Setting κ ¼ 0.78 gives the best fit with the experimental results shown in Fig. 3(b). The features occurring as ρ rr ∝ Ω 2 p increases are qualitatively consistent with the experimental data. When the system is in the I phase, dissipation is increased due to many-body dephasing, as discussed before. This case results in a lower threshold for the transition from the I phase to the NI phase compared to the reverse process. We also note an increase in the peak transmission with Ω p in the experimental data that is not reproduced by the model. This increase may indicate that ionization processes lead to a depletion of the atomic density and hence the optical depth as previously observed in cold atoms [37]. Unlike decay to the ground state, ionization leads to overall atoms loss without repopulation of ρ gg and is thus not accounted for in the model. In practice, ionization could be considered reversible on the timescales of the sweeps as atomic motion perpendicular to the excitation volume and recombination effectively repopulate ρ gg .

B. Nonequilibrium phase diagram
The enhanced sensitivity of EIT in the detection of Rydberg-induced optical bistability and the underlying phases allows us to map out the phase diagram of the driven-dissipative system and investigate the character of the nonequilibrium phase transition. We scan Δ c =2π over 96 MHz around resonance at a frequency of 10 Hz for various Ω p with Ω c =2π ¼ 13.8 AE 0.5 MHz and obtain the normalized difference in probe transmission between scan directions shown in Fig. 4(a). The phase diagram emerges from the phase boundaries that are separated by the red or black dashed line within þ or − scan direction. As the bistability is a consequence of the scan-direction-dependent occupation of the two phases, the corresponding regions represent the boundaries between the NI phase and the I phase. Outside these areas, there is no transmission difference, indicating that the system is in the same phase irrespective of scan direction, i.e., in the I phase in between the bistable regions and the NI phase outside. We observe a critical point at a threshold Rabi frequency of Ω  Fig. 3(a) and Fig. 4 results from a slight temperature difference between data sets. As above, the bistable branch is wider for negative Δ c due to the sign of Δ 0 . In the bistable branch for negative Δ c , the system occupies the I phase for scans from blue to red detuning and vice versa in the branch for positive Δ c as discussed above. Note that in the direct vicinity of the critical point, the sign of the transmission difference is switched compared to larger Ω p . This case does not correspond to the occupancy of different phases but the different regimes for the shift Δ 0 compared to the EIT linewidth as discussed above (Fig. 3).
The onset of the phase transition is characterized by a broken symmetry accompanied by a nonzero order parameter [41]. Figure 4(b) shows the shift in peak transmission Δ Ã [as defined in Fig. 3(a5)] vs Ω 2 p ∝ ρ rr ∝ N R . Depending on the origin of the underlying interactions, either the charge density (ion-induced interactions) or the mean spacing between Rydberg atoms (dipolar interactions) would represent the order parameter, both of which are ultimately related to N R . To determine the critical point, we fit and find Ω 2 p;ðcÞ ¼ 37ð2π × MHzÞ 2 with β ¼ 1. This equation refers to the critical regime near Δ c ¼ 0. The continuity at Ω 2 p;ðcÞ indicates that the system undergoes a continuous phase transition.
The observation of a clear threshold for Δ Ã has consequences for possible physical origins of the phase transition and theoretical descriptions of the system. The physical origin of the phase transition is the broadening of the Rydberg excitation above the critical Rydberg population ρ ðcÞ rr . It can be observed both in blue and red detuning; thus, the direction of bistability is not dependent on the sign of Rydberg polarizability [26]. The threshold behavior implies that the shift cannot originate directly from effects on individual Rydberg atoms, e.g., ionization via collisions with ground-state atoms, as one would expect Δ Ã ∝ ρ rr for any ρ rr in this case. Hence, the transition must result from processes that involve multiple Rydberg atoms, i.e., dipolar interactions or plasma formation. Rydberg ground-state atom collisions have been observed in thermal Rydberg ensembles as a charged mean field under EIT conditions [32] and as a dominant ionization mechanism in a beam of thermal Sr atoms, but the absence of significant shifts and broadening below ρ ðcÞ rr rules them out as the direct origin of the phase transition. However, even though they have no immediate effect on the spectra, the resulting ions and electrons are crucial for ionization avalanches and plasma formation [28].
For models describing the system, the threshold behavior also implies that a standard mean-field model [10], where the interaction is incorporated by introducing Δ 0 ∝ ρ rr , is insufficient to describe the nonequilibrium dynamics near the critical point but should remain a valid approximation for Ω p ≫ Ω p;ðcÞ as in previous experiments [24][25][26]. However, by introducing a threshold-dependent mean-field shift and broadening, the forest-fire model described in Sec. II reproduces the characteristic features of transmission spectra qualitatively well (Fig. 3). Data showing analogous behavior for the linewidth of the transmission feature at Ω 2 p;ðcÞ can be found in Appendix B. The inset in Fig. 4(b) shows the transmission difference between the phases at the transition point. We also measure the transmission difference [definition in Fig. 3(a5)] at the transition point (for negative detunings) against Ω p as Δ c is scanned. The transmission difference between the phases increases only slowly for larger Ω p , and the change in sign (see also Fig. 3) near the critical point becomes, once again, evident as Ω p becomes large enough for Δ 0 to exceed the narrow EIT linewidth.
To further characterize the phase transition, we investigate the dynamics near the critical point. This process requires that we systematically vary further parameters such as the number and excitation rate of Rydberg atoms or the interaction strength, which influence the modified decay rates, energy shift, and the transmission level. Here, we change the coupling Rabi frequency Ω c while keeping Ω p =2π ¼ 6.59 MHz constant and record the probe transmission vs Δ c [ Fig. 4(c)]. As Ω c is increased, the change in transmission dT=dΔ c for red detuning becomes steeper up to the point where the phase transition appears. Figure 4(d) shows dT=dΔ c vs Ω c and thus the susceptibility of the system to Ω c . The slope diverges at the transition threshold as expected for a continuous phase transition. In principle, it would be possible to extract a critical exponent from a similar measurement, but we refrain from this as the errors on dT=dΔ c are large and, unlike before, the condition Ω c > Ω p for EIT is no longer met, which may significantly affect the excitation dynamics compared to above. Here, the universality class for this type of phase transition cannot be determined from the fitted critical exponent ξ given in Fig. 4(d) as the nonequilibrium dynamics is so complex that the critical exponents in nature would be more various.
A more detailed analysis on this criticality can be found in the Rydberg avalanche dynamics in Sec. IV E.

C. Manipulation of the threshold
In the previous experiments, the EIT field parameter, i.e., Ω p , and its current value determine ρ rr and thus N R and the onset of the phase transition. The bistability illustrates that N R , and not any specific value of the EIT fields, is the critical quantity. This quantity allows us to further explore the threshold behavior and demonstrate control of the threshold as we apply an additional weak switching field with Rabi frequency Ω s =2π ¼ 4.1 MHz (see Fig. 1) to enhance the Rydberg population and observe the effect on the phase diagram. The switching field increases the driving strength of the jgi → jei transition leading to a higher Rydberg excitation rate without increasing Ω p itself. In the following, Δ c =2π is scanned over 88 MHz at 20 Hz and T ≈ 55°C.   [27]. To implement the switch, we change the Rydberg state to jri ¼ j70D 3=2 i, thus increasing the sensitivity while all other parameters remain unchanged. Figure 5(d) shows the transmission vs Δ c with and without switching field. When present, the transition to the I phase is reduced by approximately 2 MHz such that a switch can be implemented by fixing Δ c =2π ¼−20MHz in between the two transition points. We record the transmission over time at fixed detuning, both with and without the switching field [ Fig. 5(e)]. Figure 5(f) shows the corresponding histograms of the probe transmission for the two transmission states that highlight the change in the statistical distribution of the transmission level. Bimodal counting statistics depending on interaction strength have been predicted by theoretical investigations of Rydberg-mediated optical bistability [13]. Multiple switching fields that only intersect with parts of the excitation volume could be used to realize switching between more than two transmission states cascading multiple bistability regions with different thresholds. Analogously, a setup addressing different Rydberg states in different regions could give further insight about whether the phase transition originates from dipolar interactions or ions, similar to experiments that used different isotopes excited to different states as field probes [26] and to study energy exchange processes [42,43].
While the transmission levels are clearly distinguishable, short-term fluctuations between them can be observed [ Fig. 5(e)]. We attribute these to short-term fluctuations in the Rydberg density that result in a phase transition. Remarkably, the system quickly reverts to its initial phase, further supporting the presence of self-organization dynamics.

D. Sensitivity to fluctuations
In the following, we further investigate the sensitivity of the phase transition to fluctuations and the resulting selforganizing response. So far, we have investigated the dynamics for probe and coupling fields with parallel linear polarization. By rotating the polarization of the two beams with respect to each other, the rate at which different magnetic sublevels of the Rydberg state are populated is enhanced. As a result, the system becomes more sensitive to fluctuations in the N R -dependent interaction strength, and Γ r þ Γ 0 ðN R Þ of jri in the I phase is further increased as previously uncoupled m J states become more populated. angles θ of the half-wave plate controlling the polarization (θ ¼ 0°corresponds to parallel, 45°to orthogonal polarization) at jri ¼ j70D 3=2 i with m J ∈ fAE3=2; AE1=2g. As the parallel polarization component is reduced, the bistability window becomes narrower and is almost closed for orthogonal polarization [ Fig. 6(b)]. Interestingly, the polarization dependence of the hysteresis width cannot be observed at lower Rydberg states (jri ¼ j47D 3=2 i (see Appendix E), which further indicates that enhanced dephasing occurs at higher n due to a wider range of polarizabilities (∝ n 7 ) among the m J states that lead to a wider spread of level shifts. This case is in stark contrast to the fact that Γ r is reduced to an isolated Rydberg atom with higher n due to their longer lifetimes (∝ n 3 ).
For orthogonal polarization, the hysteresis loop is almost closed; the phase of the system is particularly sensitive to fluctuations. Transmission fluctuates between the levels associated with the two phases when scanning Δ c (Δ c =2π is scanned over 96 MHz at 10 Hz) over the bistability window [ Fig. 6(c)]. Here, even small fluctuations in ρ rr can induce a phase transition due to the strong dynamic nonlinearities in excitation and decay rates, which are even noncontinuous at the threshold ρ rr ¼ ρ ðcÞ rr and characteristic for self-organizing systems.
To characterize the response time to fluctuations, we fix Δ c =2π ¼ 25.5 MHz and record the probe transmission over time [ Fig. 6(d)]. The system is mostly found in the I phase, but occasionally, fluctuations induce a transition to the NI phase, resulting in a drop in transmission; yet it quickly reverts to the I phase again [see also Fig. 5(e)]. We investigate the distribution of response times required to restore the phase. Figure 6(e) shows the distribution of the time between adjacent phase jumps. We find that most intervals last below 100 μs, corresponding to the brief transitions to the NI phase apparent in Fig. 6(d). Measuring the width of one of the spikes yields a similar result of 110 μs (Appendix D). The transition from the I phase to the NI phase is slower than the NI to I phase transition, which occurs as an avalanche process. The inset in Fig. 6(e) shows a histogram of the transmission difference between adjacent data points. The bimodal distribution confirms that the observed fluctuations are not only artifacts originating from noise induced in data acquisition but also correspond to fluctuations between two distinct transmission states. Whether the I or NI phase is predominantly occupied depends on the value of Δ c , with detuning closer to resonance favoring the I phase. Here, the system is restored to the I phase, and the histogram is biased to higher transmission as predicted in Ref. [13]. Similar behavior can also be observed in spontaneous recovery processes in dynamical networks [44], where the network switches back and forth between two collective modes characterized by high and low network activity in analogy to the weakly and strongly interacting phases.
Finally, we analyze the Fourier spectrum of the probe transmission [ Fig. 6(f)]. For low frequencies f, the noise spectrum follows a 1=f power law corresponding to pink noise, which occurs in systems exhibiting self-organized criticality [2]. If we split the data over time into individual data sets each only including data points corresponding to the same phase, a separate analysis of the Fourier spectra yields white noise instead, confirming that the pink noise is induced by the dynamical phase transition.

E. Power-law scaling
Besides the sensitivity to perturbations and pink noise observed above, scale invariance, and power-law scaling of the avalanche size, here the fraction of the system undergoing the transition from the NI to the I phase is another characteristic feature of self-organized criticality (Sec. II). It can be experimentally observed in the magnitude of fluctuations in probe transmission, as the two phases have distinct transmission levels.
In its simplest form, the ensemble displays two nonequilibrium steady states corresponding to predominant occupation of the NI and I phases, respectively. In the experiment, the switching between these two states is controlled by varying Δ c . In the NI phase (Δ c =2π ¼ −15 MHz), transmission is saturated at the level characterized by ΔE ¼ 0 in Fig. 7, and in the I phase by ΔE ¼ 1 (Δ c =2π ¼ −25 MHz). The level of the transmission relative to the steady states, 0 ≤ ΔE ≤ 1, thus indicates the fraction of the ensemble that has transitioned to the I phase. In contrast to our simulations, we can only observe the macroscopic transmission of the entire ensemble, corresponding to the cumulated size of all microscopic I-phase clusters in the model introduced in Sec. II. The bounded distribution on avalanche sizes differs from the original Oslo model of SOC, where there is no upper cutoff [45,46], but is similar to the forest-fire model [30], which we have identified as analogous to our system. As predicted, we still observe a power-law dependence in the switching dynamics as shown in Figs. 7(b) and 7(f) as a result of balanced transfer between the nonequilibrium steady states.
Fitting the experimental data, we find power-law exponents jbj of 2.2 and 2 for Δ c =2π ¼ −25 MHz and −15 MHz, respectively, close to the regime of f=p → 0 in the model of Sec. II. Both cases are far away from the point of the macroscopic transition, and the predominance of a single phase results from a low probability for nonfacilitated transition to the I phase, f ≪ p. The signs of the exponents are opposite, as the macroscopic phases of the system are characterized by a low and high number of NI clusters, respectively. If we tune the system to Δ c =2π ¼ −19 MHz near the transition point [Figs. 7(c) and 7(d)], the evolution exhibits two peaks as shown in Fig. 7(d).
Approaching the transition point, the system dynamics become more complex, and the exponents jbj increase to values between 2 and 11 as nonfacilitated transition to the I phase becomes more likely, f=p ≥ 1, as predicted in Fig. 2(d) in Sec. II.

V. DISCUSSION AND CONCLUSION
In this work, we have studied an optically driven nonequilibrium phase transition and optical bistability in a thermal vapor of Rydberg atoms. By applying the driving fields in an EIT configuration, the sensitivity in the frequency domain was enhanced by 2 orders of magnitude compared to previous experiments, which allows us to map out the phase diagram in the vicinity of the system's critical point, where a continuous phase transition can be observed. In particular, the observed threshold behavior for the interaction-induced shift Δ 0 and broadening Γ 0 limited to the I phase at low Ω p suggests that the shift originates not directly from ionizing collisions between Rydberg and ground-state atoms but either from avalanche ionization and plasma formation, as has been confirmed for stronger excitation fields [28], or dipolar interactions. Further insight on the role of dipolar interactions or avalanche ionization on the phase transition could be made in experiments that combine the sensitivity of a weak EIT probe demonstrated here with simultaneous charge detection [26,47,48] or spectroscopy of a resulting plasma [28].
Independent of the exact mechanism, the many-body nature of the resulting Rydberg interaction is the key ingredient for the phase transition and self-organization. Besides fundamental studies, the susceptibility to noise also plays a crucial role in phase-transition-enhanced sensors for electromagnetic fields such as microwave or THz fields [27].
Our results highlight that a rich range of nonequilibrium phenomena can be studied even in a relatively simple experiment and provide observational data that help us to refine theoretical models as highlighted by the modifications needed to a standard mean-field approach. The underlying physics of self-organization dynamics observed suggest that the system could provide a controllable environment to study analogous dynamical systems from, e.g., biology [49], economics [44], and many-body physics in condensed matter [50]. Finally, the sensitivity of the phase to small fluctuations could be applied to the sensing of not only optical fields as demonstrated here but also THz fields [27], system noise, microwaves, or, more generally, any weak fluctuations of parameters linked to the EIT laser fields or external perturbations.

APPENDIX A: ENERGY SHIFTS AND BROADENING OF RYDBERG TRANSITIONS
The many-body interaction underlying the nonequilibrium dynamics may either result from pairwise dipolar interactions or ionization of Rydberg atoms, i.e., avalanche ionization [37]. Either mechanism can introduce the broadening of Rydberg transitions that accompanies the nonequilibrium phase transition between the NI and I phases. Figure 8(a) shows dipole-dipole interaction potentials for 47D 3=2 vs distance for relevant pair states. The broadening is characterized by dipole-dipole interaction energy C 3 =r 3 at the distance of closest approach r ¼ 5=9N −1=3 . Similarly, we calculate the Stark shift against ion density based on the resulting mean spacing and find that there is also an increasing spread of energies over different m J states. In the ionization case given by Fig. 8(b), shifts and broadening resulting from ionized Rydberg atoms vanish from the excitation volume at a different rate compared to the Rydberg atoms and are not explicitly included in the density matrix ρ. If the phase transition is indeed caused by ions, the model in Sec. II could be extended by introducing a charge population in ρ that is incoherently populated by including ionization rates in L.

APPENDIX B: THRESHOLD-DEPENDENT MODIFICATION OF THE DECAY RATE
In addition to the interaction-induced shift in peak transmission (Fig. 4), we also measure the width Δω EIT of the EIT transmission window as Ω p , and thus N R is increased (Fig. 9). For Ω 2 p ≥ Ω 2 p;ðcÞ ¼ 37 ð2π × MHzÞ 2 , the nonequilibrium phase transition appears. As for Δ c → Δ c þ Δ 0 ðN R Þ, we observe an increase in the I phase scaling linearly with N R ∝ ρ rr ∝ Ω 2 p , indicating a modified decay rate Γ r → Γ r þ Γ 0 ðN R Þ above threshold.

APPENDIX C: MASTER EQUATION
In the absence of Rydberg-mediated interactions, the interaction of the EIT light fields with an ensemble of threelevel atoms is governed by the master equation where ρ is the atomic ensemble's density matrix and H ¼ P k H½ρ ðkÞ the atom-light interaction Hamiltonian summed over all the single-atom Hamiltonians The Lindblad superoperator L ¼ P k L½ρ ðkÞ comprised of the single-atom superoperators accounts for the finite lifetime of jei and jri.

APPENDIX D: MEASUREMENT OF RESPONSE TIME
The transition between the two nonequilibrium phases has a characteristic response time. Treating these two assuming an electric field amplitude E ∼ ρ ion proportional to the mean ion density based on the corresponding mean spacing between ions. Calculations were made using Alkali Rydberg Calculator [51].
FIG. 9. Width of EIT transmission window vs Ω 2 p ∝ N R for the same data set as in Fig. 4(b) (n ¼ 47). As for the shift in peak transmission Δ Ã , a linear increase can be observed above the phase transition threshold Ω 2 p;ðcÞ ¼ 37 ð2π × MHzÞ 2 as expected for a modified decay rate Γ r → Γ r þ Γ 0 ðN R Þ.
phases as two energy states of the system, we use a twolevel model of the instability oscillation. The enlarged collective jump given in Fig. 10(a) illustrates the competing dynamics of decay and recovery. We solve the master equation of a simplified two-level system with a decay rate of γ ¼ 9.2 kHz corresponding to the lifetime of the 70D 3=2 Rb Rydberg state at 320 K. The corresponding lifetime τ 0 ¼ 1=γ characterizes the minimum ρ rr population required to trigger a phase transition. The interaction duration is slower than the natural response time, leading to the appearance of optical instability as in Figs. 5(e) and 6(d). Figure 10(b) shows an example of the change in probe transmission as Ω p is scanned away from the transition point while Δ c =2π ∼ −15 MHz remains fixed. Unlike at the critical point, the difference in transmission is not continuous. The sudden decrease in transmission with increasing intensity represents a self-induced opacity that is contrary to the behavior of most media, which commonly exhibit self-induced transparency. This behavior indicates a blockade effect resulting from many-body interactions. The transmission signal is noisier for high Ω p and thus higher Rydberg fractions, indicating a higher degree of instability beyond the transition to the I phase.
Finally, we measure the width of hysteresis associated with the bistability as a function of scan frequency [ Fig. 10(c)]. The hysteresis loop results from the distinct decay rates of the two phases, causing higher energy loss in the I phase. This dependence implies that the transition has a long relaxation time on a ms scale, as the higher ρ rr population supports the I phase in a self-organizing manner [1]. The width is highest near 50 Hz, with ΔT ∼ 0.48 ms. We thus estimate that the facilitation support in the I phase is strongest for ΔT=τ 0 ≈ 4.4, where τ 0 ¼ 0.11 ms represents the measured natural response time of the nonequilibrium phase transition. The natural response time of the nonequilibrium phase transition characterizes the amount of the minimum value of ρ rr population. For stronger facilitation, the corresponding hysteresis width is larger as the system is more robust to external perturbations. The timescale of τ 0 is comparable to avalanche ionization processes in cold Rydberg ensembles [37] and further indicates that the phase transition is induced by ionization processes.

APPENDIX E: POLARIZATION DEPENDENCE OF OPTICAL BISTABILITY FOR n = 47
In addition to n ¼ 70, we also investigate the influence of the relative polarization between the probe and coupling fields at a lower principle quantum number n ¼ 47. At lower n, the weaker interaction-induced shifts result in lower relative shifts between the m J sublevels of jri and thus a much weaker modification Γ 0 ðN R Þ of the decay rate. This result is reflected in a much weaker dependence of the width of the bistability window on the probe and coupling polarizations (Fig. 11).  11. (a) Recorded spectra for different polarization of the coupling beam (n ¼ 47). The red or black line corresponds the þ (−) scan direction from red to blue detunings (red) and vice versa (black). Note that θ is defined as the angle of a half-wave plate inserted in the optical path of the coupling beam. (b) Angular dependence of hysteresis width for n ¼ 47.