Phase diagram and self-organising dynamics in a strongly-interacting thermal Rydberg ensemble

Dong-Sheng Ding, 2, ∗ Hannes Busche, 4 Bao-Sen Shi, 2, † Guang-Can Guo, 2 and Charles S. Adams ‡ Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei, Anhui 230026, China. Synergetic Innovation Center of Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei, Anhui 230026, China. Department of Physics, Joint Quantum Centre (JQC) Durham-Newcastle, Durham University, South Road, Durham DH1 3LE, United Kingdom Department of Physics, Chemistry and Pharmacy, Physics@SDU, University of Southern Denmark, 5230 Odense M, Denmark (Dated: February 11, 2020)


I. INTRODUCTION
Self-organization and non-equilibrium 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, or economics [1]. Many non-equilibrium 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, produces scale-invariance and 'pink' 1/f -frequency noise, and makes systems insensitive to parameter fluctuations. Quantitative modelling of simple experimental systems displaying non-equilibrium dynamics and dynamical phase transitions is a key to enhancing our understanding of non-equilibrium 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 non-equilibrium 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 non-equilibrium 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, bi-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 roomtemperature Rydberg vapors. In this paper we investigate the non-equilibrium dynamics of a driven-dissipative ensemble of strongly-interacting Rydberg atoms in a room-temperature atomic vapour. 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 non-equilibrium system, which features two thermodynamical phases: one with low Rydberg atom density and thus no interaction (NIphase), 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 two orders of magnitude compared to previous work [24][25][26]. The allows us to probe the dynamical phase diagram with unprecedented precision, and we are able to observe previously unobserved spectral features in the vicinity of the critical point (Sec. IV A and IV B). we show that a combination of interaction-induced lineshifts 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. IV A). 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-organisation 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 fire type model.

II. BACKGROUND AND MODEL
We work with the following non-equilibrium system (Fig 1(a)): An probe beam propagates through a thermal vapor of three-level atoms and drives the transition between a ground state |g and a low-lying, short-lived excited state |e . A coupling field couples |e to a highly excited and long-lived Rydberg state |r . 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 |e (|r ). 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 (∆ p = 0) while scanning ∆ c , Ω p , or Ω c .
The interaction between Rydberg atoms introduces a strong dynamical non-linearity as required for selforganisation. Moreover, SOC relies on an interactioninduced avalanches in the non-equilibrium 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 non-interacting (NI-) phase (with Rydberg density N R < N R,(c) ) and a strongly interacting (I-) phase (N R > N R,(c) ). In the latter, interactions cause an effective, N R -dependent detuning of the coupling field, ∆ c → ∆ c + ∆ (N R ) (Fig 1b). The |nD J -Rydberg states used in this work are also subject to broadening and enhanced dephasing Γ r → Γ r + Γ (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] (appendices 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 excitation. This 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 forest fire 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 corresponds to the spreading of fire to adjacent trees, see Fig. 2(a1-a3). Following the excitation avalanche, the Rydberg population in a cell is depleted to zero. Non-facilitated 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 neighbouring cells [13]. The closed boundary of our array reflects finite excitation volume in experiment as defined by the laser beams. Fig. 2(b) shows the distribution of phases in a single plane following 30 iterations. Fig. 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 scaleinvariant, 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 |b| from 1.7 ∼ 4 for p = 0.01, 0.05 and 0.1. For larger p, the occurence of large cluster sizes increases. We also vary f and record |b| against p, see Fig. 2(d), and find |b| close to 2 in the limit f /p → 0. We have also simulated the fire pixel fraction against forest density with assuming the much longer lifetime of the burning trees, which corresponds to the case of fast scanning ∆ c or Ω p in the experiment. We observed an obvious threshold effect in theory, see Fig. 2(e), the nonlinear response in the jump interval 0.3 ∼ 0.35 predicts a critical threshold Rydberg population.

B. Rydberg EIT in a thermal vapour
To model the effect of the dynamical phase transition on the probe transmission, we compute the evolution of the ensembles' denstiy matrix ρ by adopting a thresholdmodified mean-field approach where the additional detuning ∆ = η 1 N R and the enhanced decay Γ = η 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 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 where ∆ D = (ω p − ω c )v/c denotes Doppler shift experienced by an atom moving with velocity v and δ = ∆ c +∆ p the two-photon detuning. Ω ef f 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 wavevector of the probe field. To model the bistability in the EIT spectra, 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 inter-atomic spacing of ≈ 1 µm. The probe beam (1/e 2waist radius ≈ 500 µm) couples |g = |5S 1/2 , F = 3 to |e = |5P 1/2 , F = 2 . The coupling beam (1/e 2 -waist radius ≈ 200 µm) is resonant with the transition from |e to |r = |nD 3/2 . The estimated peak values of Ω c /2π are ≈ 20 MHz for n = 47 and ≈ 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 non-equilibrium 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 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 |g and |e 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 lineshape as ∆ c is scanned through resonance from negative to positive detuning and vice versa. Fig. 3 shows the results for a range of different Ω p as well as results of the threshold-model (b, dashed lines indicate where the system is in the Iphase). At Ω p /2π = 6.3 MHz (Fig. 3 a1) identical EIT spectra are observed for either scan direction without bistability. At 6.7 MHz -and consequently only slightly higher Rydberg density N R -we observe almost symmetric optical bistability above and below resonance (Fig.  3a2). 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 ∆ induced on |r is small compared to both ∆ c and the EIT linewidth. This behavior at low Ω p is distinctively different compared to previous experiments on Rydberginduced 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-(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.5 MHz, two orders of magnitude narrower compared to previous experiments [24,26].
As Ω p is increased (Fig. 3a3 to 3a5), the spectra become asymmetric as an increasingly strong spectral shift ∆ * (see definition in Fig. 3a5) of the peak transmission relative to ∆ c /2π = 0 MHz is observed. The shift increases with Ω p as N R ∝ ρ rr ∝ Ω 2 p . At Ω p /2π = 6.8 MHz (Fig. 3a3), bistability is still observable within the EIT resonance feature, but the phase transition occurs at lower/higher ∆ c value for positive/negative scan directions respectively. The on-resonance phase transition disappears as the system is in the I-phase for either scan direction.
For ∆ 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-is now exceeding that in the NIphase 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 ∆ implies an overall red-shift of the spectra due to the choice of |r . 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. Fig. 3(b), shows the calculated EIT spectra with Doppler-averaging as given by equation 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 calcu-lated 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. 3b. The features occurring as ρ rr ∝ Ω 2 p increases are qualitatively consistent with the experimental data. When the system is in the Iphase, dissipation is increased due to many-body dephasing as discussed before. This results in a lower threshold for the transition from the I-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 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 as reversible on the timescales of the sweeps as atomic motion perpendicular to the excitation volume and recombination effectively repopulate ρ gg .

B. Non-equilibrium phase diagram
The enhanced sensitivity of EIT in the detection of Rydberg induced optical bistability and the underlying phases allows to map out the phase diagram of the drivendisspative system and investigate the character of the non-equilibrium 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 ± 0.5 MHz and obtain the normalized difference in probe transmission between scan directions shown in Fig. 4(a). The phase diagram is emerged by the phase boundaries that are separated by the red/black dashed line within +/-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-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 Ω (c) p /2π = 5.9 ± 0.2 MHz with ∆ (c) c /2π = 0 MHz. The spectra shown in Fig. 3(a1) and (a2) correspond to regimes slightly below and above Ω  Fig. 3a 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 ∆ . In the bistable branch for negative ∆ c , the systems 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 does not correspond to the occupancy of different phases, but the different regimes for the shift ∆ compared the EIT linewidth as discussed above (Fig. 3).
The onset of the phase transition is characterized by a broken symmetry accompanied by a non-zero order parameter [41]. Fig. 4 (b) shows the shift in peak transmission ∆ * (as defined in Fig. 3a5) 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 ultimately related to N R . To determine the critical point for, we fit for Ω 2 p < Ω 2 p,(c) and find Ω 2 p,(c) = 37(2π × MHz) 2 with β = 1. This equation refer 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 appeared is not dependent on the sign of Rydberg polarisability [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 an 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 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 ∆ ∝ ρ rr , is insufficient to describe the non-equilibrium 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 thresholddependent mean-field shift and broadening, the forest  values) when scanning ∆c in either direction, effectively representing the phase diagram of the nonequilibrium system. In the white regions, the system is found in the same phase independent of the scan direction. In the shaded regions, the system is bistable (red if transmission is higher for scans from negative to positive ∆p, blue in the opposite case) and occupies different phases depending on the scan direction. (b) Shift of peak transmission ∆ * vs. Ω 2 p ∝ ρrr. Inset: Transmission difference at the transition point for different Ωp (see Fig. 3a5 for definitions). (c) Measurement of transmission spectra vs. ∆c for Ωc/2π between 2.7 and 6.54 MHz, here Ωp/2π = 6.59 MHz. (d) Susceptibility of the phase transition to variations in Ωc measured in terms of the change in transmission dT /d∆c. The shaded region indicates where the system is in the I-phase and dT /d∆c is no longer continuous at the transition point. The fit is given by y = 1.97e (Ωc/ξ) + 1.82, ξ = 0.38 is the fitted critical exponent fire model described in section 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. 3a5) at the transition point (for negative detunings) against Ω p as ∆ c is scanned. While 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 ∆ to exceed the narrow EIT linewidth.
To further characterize the phase transition, we investigate the dynamics near the critical point. This requires to 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. 4c). 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. Fig. 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 an 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 phase transition cannot be determined from the fitted critical exponent ξ given in Fig. 4(d)) as the non-equilibrium dynamics is such complex that the critical exponents in nature would be more various. The more detailed analysis on this criticality can be found in the Rydberg avalanche dynamics in the Sec. IV E.

C. Manipulation of the threshold
In the previous experiments the EIT field parameters, 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 allows 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 |g → |e 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. Fig. 5 (a) and (b) compare the phase diagram without (compare Fig. 4a) and with application of the switching field. The general structure of the phase diagram with the switching field applied remains unchanged, however the threshold and thus the critical point of the phase diagram appear at Ω (c) p /2π ≈ 4 MHz instead of ≈ 5 MHz as ρ rr is enhanced by the additional driving. We also measure the spectral shift ∆ * in the I-phase (Fig. 5c,  d) and again find the critical point at a lower value of (Ω (c) p /2π) 2 = 16.2 MHz 2 (Fig. 5d) compared to 24.3 MHz 2 without switching field (Fig. 5c) as the relation ρ rr ∝ Ω 2 p no longer holds. Without effective many-body interaction, all atoms in the vapor could be described individually and the switching field should not affect the threshold. The threshold switch process cannot be demonstrated by varying the coupling beam because we cannot obtain the phase diagram with a threshold by only sweeping Ω c .
The ability to shift the critical value Ω (c) p via the switching field can be used to implement an optically controlled switch between two transmission states for the probe light (Fig. 5d-f) similar to experiments on latching detection of THz-fields using Rydberg optical bistability [27]. To implement the switch, we change the Rydberg state to |r = |70D 3/2 thus increasing the sensitivity while all other parameters remain unchanged. Fig  . 5(d) shows the transmission vs. ∆ c with and without switching field. When present, the transition to the I-phase is reduced by approx. 2 MHz such that a switch can be implemented by fixing ∆ c /2π = −20 MHz in between the two transition points. We record the transmission over time at fixed detuning both with and without the switching field (Fig. 5e). 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 whether the phase transition originates from dipolar interactions or ions, similar to experiments which 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. 5e). 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 selforganization dynamics.

D. Sensitivity to fluctuations
In the following, we further investigate the sensitivity of the phase transition to fluctuations and the resulting self-organizing 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 sub-levels 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 + Γ (N R ) of |r in the I-phase is further increased as previously uncoupled m J -states become more populated. Fig. 6(a) shows transmission scans vs. ∆ c for various angles θ of the half-wave plate controlling the polarization (θ = 0 • corresponds to parallel, 45 • to orthogonal polarization) at |r = |70D 3/2 with m J ∈ {±3/2, ±1/2}. As the parallel polarization component is reduced, the bistability window becomes narrower and is almost closed for orthogonal polarization (Fig. 6b). Interestingly, the polarization dependence of the hysteresis width cannot be observed at lower Rydberg states (|r = |47D 3/2 , see appendix D), 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 is on stark contrast to the fact that Γ r is reduced 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. 6c). 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 non-continuous 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. 6d). 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 quickly reverts to the I-phase again (see also Fig. 5e). We investigate the distribution of response times required to restore the phase. Fig. 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 C). The transition from the I-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 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 restores to the I-phase and the histogram is biased to higher transmission as predicted in [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 characterised by ∆E = 0 in Fig. 7, 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 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. 7b and f as a result of balanced transfer between the non-equilibrium steady states.
Fitting the experimental data, we find power law exponents |b| 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 system's macroscopic phases of the system are characterised by a low and high number of NI-clusters, respectively. If we tune the system to ∆ c /2π = −19 MHz near the transition point (Fig. 7c,  d), the evolution exhibits two peaks as shown in Fig. 7d. Approaching the transition point, the system dynamics become more complex, and the exponents |b| increase to values between 2 and 11 as non-facilitated 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 two orders of magnitude compared to previous experiments which allows to map out the phase diagram in the vicinity of the system's critical point where a continuous phase transitions can be observed. In particular, the observed threshold behavior for the interaction induced shift ∆ and broadening Γ 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 it has been confirmed at 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 electro-magnetic 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 helps 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] 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. Rydberg atomic dipoles are randomly aligned, resulting in the I-phase. Dipoles align if normalized rate of energy change is less than broadening. Broadening is dipoledipole interaction at distance of closest approach C 3 /r 3 , where r = 5/9N −1/3 . Scanning the system from the I-to NI-phase resulting in bifurcation, and vice versa. We also calculate the stark shift of Rydberg D-state against ions density, and find that there are also a broadening structure in energy diagram. In the ionization picture of the bifurcation of the phase diagram, the shifts and broadening result from ionized Rydberg atoms vanish from the excitation volume at a different rate compared to the Rydberg atoms and are not explicitly included in ρ. If the phase transition is indeed caused by ions, the model in Sec. II induced in the beginning could be extended by introducing a charge population in ρ that is incoherently populated by including ionization rates in L.
Although the pair dipole-dipole interaction potentials for S state has no obvious splitting structure. This shows that the I-phase of S state is narrow than the D state, which is consistent with the experimental observations. The broadened spectra of S state in I-phase can be caused by the Rydberg atoms with distinct interaction distance and different density and so on. minimum amount high-ρ rr population to generate phase transition. The duration of interaction is slower than the natural response time, thus it appears optical instability like in Fig. 5e and Fig. 6d. Fig. 10(b) shows an example of the change in probe transmission as Ω p is scanned either way away from the from transition point for while ∆ c /2π ∼ −15 MHz remains fixed. Unlike at the critical point, the difference in transmission is not continuous. It displays a sudden fall in the transmission. A decrease in transmission with increasing intensity-a kind of self-induced opacity that is contrary to the behavior of most media, which normally exhibits self-induced transparency; this behavior is indicative of a certain type of many-body blockade effect from many-body interaction, not analogy to the character of single-body blockade effect. An interesting feature noted in Fig. 10(b) is that the noise in the transmission signal is stronger for high Ω p , and thus higher Rydberg fractions, indicating a higher degree of instability beyond the transition to the high-ρ rr phase.
Finally, we measure the width of hysteresis associated with the bistability as a function of scan frequency, the results are shown in Fig. 10(c). The reason why there is a hysteresis cycle is that the system undergoes two dissipate processes with distinct decay rates, causing energy loss in the I-phase. The dependence on sweep speed implies the transition has a long relaxation time on the millisecond scale, which characterizes an amount of high-ρ rr population, supporting the structured phase transition. This is attributed to a self-organization dynamics, a nonequilibrium transition from the NI-to I-phase [1]. Alter- natively, there is a saturation point occurred around 50 Hz, the corresponding width of hysteresis is ∆T ∼ 0.48 ms. We thus estimate the maximum amount of attractor for our system ∆T /τ 0 ≈ 4.4, where τ 0 = 0.11 ms represents the measured natural response time of nonequilibrium phase transition from minimum interacted atoms, giving four times larger than the threshold of attractor in our case. The natural response time of nonequilibrium phase transition characterizes the amount of minimum high-ρ rr population, the more the high-ρ rr population, the corresponding hysteresis width is larger, thus the system is more robust to the external interference or perturbation. This timescale of τ 0 is comparable to that of avalanche ionization processes in cold Rydberg ensembles [37] and could indicate that the phase transition is induced by ionization processes. 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 |r and thus a much weaker modification Γ (N R ) of the decay rate. This is reflected in a much weaker dependence of the width of the bistability window on the probe and