Emergence of synchronisation in a driven-dissipative hot Rydberg vapor

We observe synchronisation in a thermal (35-60 {\deg}C) atomic (Rb) ensemble driven to a highly-excited Rydberg state (principle quantum number n ranging from 43 to 79). Synchronisation in this system is unexpected due to the atomic motion, however, we show theoretically that sufficiently strong interactions via a global Rydberg density mean field causes frequency and phase entrainment. The emergent oscillations in the vapor's bulk quantities are detected in the transmission of the probe laser for a two-photon excitation scheme.

Nonlinear systems are abundant in nature, where the nonlinearities introduce a range of rich and varied phenomena.Well known is the ability of nonlinear systems to generate multiple steady states, so that the system's state is determined by its past trajectory and hysteresis loops may form.Such multistable states have been observed numerously in biological [1][2][3][4], mechanical [5][6][7], and atomic systems [8][9][10][11].Nonlinear dynamics and bifurcation theory provide a modelling framework of these phenomena, enabling a fundamental understanding of the underlying processes from within a generalised mathematical framework.
When adding dissipation to a conservative nonlinear system, the resulting dynamics get even richer and the system can support rather unexpected types of stable solutions.Under certain conditions, dissipative systems with nonlinearities can support chaotic behavior [12,13] or limit cycles and time-periodic solutions [14,15].A Hopf bifurcation may cause the appearance of attractive limit cycles, which leads to self-sustained oscillations of the system.This oscillatory behavior is not imprinted by an external drive but arises fundamentally from the system's dynamics.Such self-oscillating systems have been found to model biological processes [16][17][18][19][20] and physical systems [21][22][23][24].
A very curious question regards the behavior of an ensemble of self-sustained oscillators experiencing a form of coupling to another, or to an external force.First studied by Kuramoto for an ensemble of globally coupled oscillators with different natural frequencies [25], it has been found that -under certain conditions -all or a subset of the oscillators begin to lock in frequency and phase [26][27][28].As a result, a transition towards a sychronised state occurs in the ensemble.This synchronisation transition has been used to explain e.g. the strong lateral vibrations of the Millennium bridge, London, on its opening day [29], though this is contested [30], or the Belousov-Zhabotinsky and other chemical reactions [31,32].In nature, synchronisation occurs in ensembles of fireflies flashing in unison [33], the chirps of snowy tree crickets [34], and occasionally in the applause of audiences [35].
To further study the emergence of synchronisation and the resulting non-equilibrium dynamics, a simple and easily controllable system with a macroscopic number of coupled oscillators and tunable properties is highly desirable.In the following, we demonstrate that the occurrence of a synchronised phase is expected in a continuously driven, dissipative three-level system with a power law coupling to a mean field, and report on the observation of synchronisation in a hot Rydberg vapor.A surprising, but expected, feature of this system is that oscillations of the bulk quantities remain observable even though the individual constituents are undergoing random motion.
Rydberg atoms are known to interact strongly with a power law scaling in distance.This translates into a mean-field approach [36] with power law scaling β of the Rydberg level shift in Rydberg density ρ rr .A similar power law scaling also can be used to model the level shift induced by ionization [37] or other mean-field inducing mechanisms.Adopting this mean field approach, the resulting equations of motion (EOMs) are formulated for a three-level basis set with coherent driving by Ω x and dissipation Γ yz , see figure 1 (a).For β ̸ = 0, the EOMs are nonlinear and their steady states are defined by the roots of a polynomial of order max(4β + 1, 1) in ρ i er .The resulting steady-state solutions of the nonlinear EOMs reveal regions of multistability where an odd number of equilibria exist for one set of parameters Ω x , ∆ x , Γ yz , V, β.To extract the stability of the solutions, the spectrum of eigenvalues λ j of the linearisation (Jacobi) is evaluated at the steady state [38].Stability is guaranteed if Re(λ j ) < 0 for the eight non-constant eigenvalues.Consequently, the repulsive branch marked in red in figure 1 (b) is detected by spectral analysis.However, the steady states indicated in green are also unstable.Here, a Hopf bifurcation occurs where a complexconjugate pair of eigenvalues λ j crosses the imaginary axis and renders the steady state unstable.As a result, the system is attracted towards a limit cycle which leads to robust self-sustained oscillations of the system parameters in time.Figure 1 (c) and (d) show that the system is attracted to the same limit cycle for different initial states, but each initial state leads to a different phase in the limit cycle at any fixed time t.This freedom of phase in the limit cycle is indicative of a self-oscillating system and fundamentally distinguishes it from a periodically driven system where the phase in the limit cycle is locked to that of the drive.The freedom of phase in the resulting limit cycle has also been described using the language of continuous time crystals [39,40].The timecrystal interpretation in the context of our experiment is discussed in appendix E of the supplementary material.
Although optical bistability has been found experimentally in driven-dissipative hot Rydberg vapors [10], one would intuitively expect any oscillations in this system to average out due to atomic motion.The motioninduced dephasing for different atomic velocities results in a spread of the natural frequencies of the limit cycles and the phases therein.Although about half of the velocity classes are attracted towards a limit cycle, no macroscopic oscillations can be seen, as shown by the black line in figure 2 (a).
However, above argumentation does not account for the spatial dimension of the situation.The Rydberg level shift of any atom in the vapor depends on the spatial Rydberg density of its local environment so that the different velocity classes do not evolve independent of another.Rydberg atoms of one velocity class experience a level shift depending on the Rydberg population of the other velocity classes in the vapor and, in turn, influence the dynamics of these other velocity classes.When tak- ing this global coupling between the velocity classes into account, the resulting dynamics of the vapor is very different as shown in 2 (b) [see also App.C in Supp.Mat].After an initial transient phase, synchronisation sets in where the velocity classes begin to oscillate in lockstep with a single frequency and fixed phase relation.This is possible because the phase of a velocity class within its limit cycle is free and therefore easily adjusted by the mean field.With a growing number of velocity classes oscillating in phase lock, the mean field strength increases which forces even more velocity classes to align their oscillations until eventually a partially or completely synchronised state is reached.
This transition towards a synchronised state of globally coupled oscillators is known since Christiaan Huygens' time [41] and has since been studied extensively from a mathematical perspective.After the initial work by Winfree [42] and Kuramoto [25], the study of synchronisation has been extended to more general forms of the global coupling force [27,28] and other situations [43].Famous examples where synchronisation is experimentally demonstrated for few oscillators is the synchronisation of pendulum clocks [41] or metronomes [44] fixed to a common support which provides the coupling.However, large numbers of globally coupled oscillators with widely tunable properties are not so easily available.Therefore, a hot Rydberg vapor with ∼ O(10 9 ) atoms in the beam volume, and a somewhat lower number of oscillators, provides an ideal testbed for an experimental study of the synchronisation transition for large numbers of constituent oscillators.
In our experiment, we use 87 Rb number densities of ρ 87Rb ∈ [0.1, 6.1] • 10 11 cm −3 , which corresponds to tem- for Rydberg states with principal quantum numbers n ranging from 43 to 79.Different beam waists of up to w ≤ 1 mm and beam waist ratios of w p /w c ≈ 2, 0.9, 0.5 have been tried, but no direct dependence on the beam waists has been observed.The data presented here was obtained for w p = 390 µm and w c = 440 µm.Setup and relevant level scheme are shown in figure 3 (a) and (b).
Figure 3 (c) shows a typical series of scans for fixed probe and increasing coupling Rabi frequency.After an onset of bistability in the optical response, a window featuring oscillations in the vapor transmission opens.This synchronisation window widens for a further increase in coupling Rabi frequency.When instead setting the coupling Rabi frequency to a fixed value, the width of the oscillation region decreases with increasing probe Rabi frequency (see also App.D in Supp.Mat.).In the various parameter regimes that were explored experimentally, the synchronisation regime is often preceded by bistability but not necessarily so.We find a strong dependence of the onset of oscillations on the Rydberg state and vapor density.Higher atom number densities require lower Rabi frequencies for the oscillations to set in.This behavior is expected from a synchronisation perspective since larger global coupling strengths require lower meanfield strengths to initiate entrainment.
We observe an onset of synchronisation for coupling to both nS and nD Rydberg states, though it is easier to explore the behavior and scaling when coupling to D states due to the stronger dipole coupling at similar n.The oscillations were also observed when coupling a fourth P or F state with an additional rf field in both the weak and strong driving limit, respectively.In the fully Autler-Townes split regime, oscillations occurred as long as the Rydberg population was high enough.The presence of synchronisation is therefore neither a purely three-level phenomenon, nor does it depend on the orbital angular momentum of the Rydberg state.
With all system parameters held constant and fixed laser detunings, the synchronised state persist on timescales on the order of minutes and the oscillations maintain their shape.Analysis of a time trace reveals a narrow frequency peak with a spectrum of weaker, higher harmonics (also shown in App.D of Supp.Mat.).The oscillation frequency ν osc of the first peak was usually observed to lie between 10 kHz and 25 kHz, though persistent oscillations of up to 43 kHz were measured.In Fig. 3 (c) one can see that the oscillation frequency varies along the coupling laser scan.As a general trend, an increase in oscillation frequency ν osc with increasing Rabi frequencies was observed.Additionally, the formation of several separate synchronisation regions, typically with a different range ν osc but similar shapes of the oscillations along the region, has been found.This is also Figure 4 (a) shows the change in oscillation shape and frequency with increasing ∆ c .Each highlighted segment samples the time dependence at a particular detuning as the laser frequency is scanned in time slowly relative to ν osc .The rightmost zoom-in (red) belongs to the next synchronisation region beginning at ∆ c /2π ≈ −45 MHz.It shows again the sawtooth-like shape at its lower frequency end that can also be seen in the two leftmost insets.Figure 4 (b) -(d) show results obtained with the thermal vapor simulation.The imaginary part of the coherence ρ i ge shown in (c) and (d, dashed) is linearly proportional to the probe laser transmission via the probe electric susceptibility [45].Two limit cycle regions appear in the spectrum (b), though a cross-section of phase space shows that the case ∆ c /Γ ge = −3 is not a limit cycle but resembles a system near a strange attractor.Generally, the thermal vapor model shows regions of multistability which implies that the pointwise integration technique in (b) cannot accurately model a laser scan.This is because the thermal vapor system's trajectory depends on its past state and the attractor it is drawn to, which pointwise integration does not account for.
The thermal vapor model reproduces the observed experimental behavior phenomenologically.This includes changes in the width of the synchronisation region with changes in Ω p or Ω c and the earlier onset of oscillations at lower Ω c for increasing interaction strengths V as shown in the data of Appendix D in the Supp.Mat., as well as the expected shape of the oscillations.Therefore, we attribute the emergence of macroscopic oscillations in the bulk response of a hot Rydberg vapor to a Kuramotolike synchronisation transition for sufficiently large global coupling strengths.Possible mechanisms causing the power law scaling of the Rydberg density mean field are Rydberg interactions [36] or charge-induced Stark shifts due to ionisation [37], though other effects could possibly lead to similar power-law scaling behaviors.
In summary, we observe the transition towards synchronisation in a strongly driven, dissipative, hot Rydberg vapor.The observed changes of the synchronised region with variation of the Rabi frequency, vapor density, and interaction strength is reproduced by a theoretical model extended to a thermal vapor simulation.The model's nonlinearity leads to the emergence of attractive limit cycles for individual velocity classes through a Hopf bifurcation.Under the influence of global coupling through the shared Rydberg density, the constituent oscillating velocity classes synchronise in a thermal vapor, which leads to periodic oscillations of the vapor's bulk quantities.The resulting synchronised phase is robust and stable, and therefore ideally suited for an experimental investigation of the emergent non-equilibrium phase of matter.It provides a simple platform for the study of synchronisation in a nonlinear system with a truly macroscopic number of oscillators.
Author's note: During completion of this work, two other reports of oscillations in a continuously driven hot Rydberg vapor were reported.In [46], the oscillations are of a transient nature and the probe Rabi frequency is significantly lower than in this work.The authors attribute the origin of the limit cycles to spatial inhomogeneities and clustering of Rydberg atoms.In [47], the experimental parameter regime is similar to this work.The limit cycles are attributed to a competition for Rydberg population between energetically closely spaced Rydberg states.
The Routh-Hurwitz criterion [49] for a cubic P(x) = x 3 + a 2 x 2 + a 1 x + a 0 states that all eigenvalues of P have a negative real part except for one purely imaginary pair iff a 2 > 0, a 0 > 0 and a 2 a 1 − a 0 = 0 are true.If these conditions are satisfied, a Hopf bifurcation occurs.
From B3 one sees immediately that a 2 = 2Γ > 0 is satisfied for any Γ > 0, as we assume in a dissipative system with non-vanishing population loss from state |r⟩.Looking at the last of the three conditions, the equality, then one finds that a 2 a 1 − a 0 = a 0 + Γ(Ω 2 + 2Γ 2 ) = 0 cannot be satisfied if the condition a 0 > 0 is true.Therefore, no complex conjugate pairs of eigenvalues crosses the imaginary axis and no Hopf bifurcation occurs in the effective two-level model.

Appendix C: Thermal vapor integration scheme
In a thermal vapor model, the different velocity classes v j do not evolve independent of another.This is because the Rydberg atom density generates an effective mean field that all velocity classes couple to and are influenced by.Simple integration of the individual velocity classes therefore cannot account for the Rydberg density produced by other velocity classes, and the integration scheme has to be adapted accordingly.The mean field Hamiltonian for a given velocity class v j changes from H vj ∝ (ρ vj rr ) β |r vj ⟩ ⟨r vj | to H tot vj ∝ (ρ tot rr ) β |r vj ⟩ ⟨r vj |.To this end, an adapted rk4 integration scheme [50] has been implemented.After preforming a single rk4 time step for all velocity classes, the Rydberg population of the vapor is computed as the weighted sum over the Rydberg population of all velocity classes.Then, the density-dependent level shift V (ρ tot rr ) β is adjusted in the EOMs for all velocity classes and the next integration step is performed.A matrix-based implementation of the integration scheme in Python allows for a simple parallel computation of one time step for all velocity classes at once.An example time trace for all system parameters held constant is shown in (c) and the corresponding spectral density in (d).This Fourier transform analysis was used to determine the oscillation frequency ν osc .
The upper row of figure 6 shows the synchronisation window for three different vapor densities but otherwise identical system parameters with a coupling to the |79D⟩ states.The vapor temperatures were varied from T 1 = (40.(blue) to T 2 = (51.0± 0.5) • C (yellow) and T 3 = (52.0± 0.5) • C (red).By varying the vapor temperature, the critical coupling Rabi frequency Ω c where synchronisation sets in is shifted for fixed probe Rabi frequencies and Rydberg state.Therefore, a direct dependence of the coupling strength V on the vapor density can be inferred.A higher vapor density leads to stronger effective Rydberg interactions and ion densities at the same Rabi frequencies due to the reduced interatomic spacing.For higher vapor densities, the onset of synchronisation is therefore shifted to lower Rabi frequencies where the critical Rydberg atom density is then located.
In the lower row of figure 6, the synchronisation window is shown for different combinations of Rydberg state but similar vapor densities.At the various probe Rabi frequencies Ω p , the minimum coupling Rabi frequency required for an onset of synchronisation is different for all three Rydberg states.Higher principal quantum numbers n need lower Rabi frequencies for meeting the synchronisation threshold.The states were selected by Rydberg interaction strength since the |43D 5/2 ⟩ state is on a Förster resonance, leading to higher Rydberg interactions than for the |50D 5/2 ⟩ state [51], while the ionisation cross-section of the lower-n Rydberg state is lower [37].The |63D 5/2 ⟩ state has the highest interaction strength and collisional ionisation rates.A direct effect of the Förster resonance of the |43D 5/2 ⟩ state cannot be seen in the data.
In the green |63D 5/2 ⟩ data, the occurrence of spectrally separate synchronisation regions is clearly visible.

FIG. 1 .
FIG. 1. Single velocity class model.The basic model with the relevant parameters is shown in (a).An example steadystate solution of the resulting nonlinear OBEs is shown in (b) where the dark-red steady-state branch is repulsive and green indicates the limit cycle region.For a fixed detuning ∆c/Γge = −1, indicated by the dashed line, the time evolution from an initial state |Ψ⟩ t=0 = (1 − x) |g⟩ + x |r⟩ with x ∈ [0, 1] towards a limit cycle is shown in (c).For the same time traces, a phase space projection of the limit cycle in the ρge-plane is shown in (d).The other model parameters were set to: ∆p = 0, Ωp/Γge = 3.8, Ωc/Γge = 2, V /Γge = −12, Γer/Γge = 10 −5 , Γgr/Γge = 10 −2 and β = 3.

FIG. 3 .
FIG. 3. Setup and example onset of oscillations.(a) The counterpropagating probe and coupling lasers are polarisation cleaned with a polarising beamsplitter (PBS) after exiting the fibers.The subsequent acousto-optic modulator (AOM) and aperture are used to remote control the laser powers incident on the heated, 4 cm long rubidium cell.The probe light is detected by a photodetector (PD).(b) shows the relevant level scheme for two-photon spectroscopy in rubidium.The coupling laser addresses either the |n ′ S 1/2 ⟩ or the |nDm J ⟩ state(s).In (c), an example set of traces obtained for fixed probe Rabi frequency Ωp/2π = 191 MHz and increasing coupling Rabi frequencies is shown.The Rydberg laser is coupled to the |43D 5/2 ⟩ state, and the number density is ρ 87Rb = (4.7 ± 0.2) • 10 10 cm −3 .Here, the oscillatory regime is preceded by an onset of bistability.

FIG. 4 .
FIG. 4. Change in oscillation shape and frequency along coupling laser scan.(a) shows the oscillation region for a scan across resonance with |43D 5/2 ⟩ at T = (52.0±0.5)• C with Ωp/2π = 191 MHz, Ωc/2π = 37 MHz, and a scan rate of 2π ×10 MHz/ms.The colored insets show a zoom-in of the trace in the color-shaded regions, each of width 2π × 4.8 MHz.Different shapes of the oscillations can be distinguished.(b) Pointwise integrated spectrum with errorbars denoting the amplitude of the oscillations.The time evolution towards a limit cycle is shown in (c) with the inset showing only the limit cycles approached after an integration time of t = 5000Γ −1ge .In (d), the oscillations in Rydberg population ρ r rr (solid) and in the imaginary part of the coherece ρ i ge (dashed) are shown.The case ∆ = −3Γge did not approach a limit cycle within the maximum integration time but behaves similar to a system near a strange attractor.The simulation assumes a thermal vapor with N vel = 101 velocity classes with equal populations, Ωp = 1.5, Ωc = 1, ∆p = 0, Γer = 10 −6 , Γgr = 10 −3 and V = −300, in units of Γge, and β = 2.
Appendix D: Supplementary experimental data

Figure 5
Figure 5 shows the observed behavior of the synchronisation window as well as the change in oscillation frequency when keeping coupling (a) or probe (b) Rabi frequency constant and changing the respective other.The observed closing and opening of the synchronisation window, as well as the change in oscillation frequency, are reproduced by the thermal vapor model.The scaling behavior in the thermal vapor model is indirectly inherited from the three-level single velocity class model, which initially seeds the oscillations.An example time trace for all system parameters held constant is shown in (c) and the corresponding spectral density in (d).This Fourier transform analysis was used to determine the oscillation frequency ν osc .The upper row of figure6shows the synchronisation window for three different vapor densities but otherwise identical system parameters with a coupling to the |79D⟩ states.The vapor temperatures were varied from T 1 = (40.5 ± 0.5) • C

FIG. 5 .FIG. 6 .
Figure 5 shows the observed behavior of the synchronisation window as well as the change in oscillation frequency when keeping coupling (a) or probe (b) Rabi frequency constant and changing the respective other.The observed closing and opening of the synchronisation window, as well as the change in oscillation frequency, are reproduced by the thermal vapor model.The scaling behavior in the thermal vapor model is indirectly inherited from the three-level single velocity class model, which initially seeds the oscillations.An example time trace for all system parameters held constant is shown in (c) and the corresponding spectral density in (d).This Fourier transform analysis was used to determine the oscillation frequency ν osc .The upper row of figure6shows the synchronisation window for three different vapor densities but otherwise identical system parameters with a coupling to the |79D⟩ states.The vapor temperatures were varied from T 1 = (40.5 ± 0.5) • C

FIG. 7 .
FIG. 7. Phase drift in the limit cycle during time sequence.(a) The frequency spectrum of oscillations on the |63D 5/2 ⟩ state with an oscillation frequency νosc of 43.1 kHz.Two segments of the 0.1 s long data trace are shown in (b) together with a triangular reference waveform of frequency νosc in gray.The bottom time segment shows a phase shift relative to the reference waveform and a change in oscillation frequency along the length of the segment.This leads to a drift of the phase of the signal relative to the reference, also shown in (c), where the extracted phases of the entire sequence are shown in gray.The radius is proportional to the measured oscillation frequency.