Dynamics of nonequilibrium Dicke models

Motivated by experiments observing self-organization of cold atoms in optical cavities, we investigate the collective dynamics of the associated nonequilibrium Dicke model. The model displays a rich semiclassical phase diagram of long-time attractors including distinct superradiant fixed points, bistable and multistable coexistence phases, and regimes of persistent oscillations. We explore the intrinsic time scales for reaching these asymptotic states and discuss the implications for finite-duration experiments. On the basis of a semiclassical analysis of the effective Dicke model, we find that sweep measurements over 200 ms may be required in order to access the asymptotic regime. We briefly comment on the corrections that may arise due to quantum fluctuations and states outside of the effective two-level Dicke model description.

Figure 1

Experimental setup showing cold atoms (red) in an optical cavity with transverse pumping [20, 21, 22, 23]. (a) Below the threshold pump power only the pump mode is present. (b) Above the threshold the atoms self-organize into a checkerboard lattice and are trapped in the interference pattern of the pump and cavity beams. The self-organization transition for a BEC is described by the onset of superradiance in an effective nonequilibrium Dicke model.

Figure 2

Level scheme corresponding to the experimental setup shown in Fig. 1. The pump beam has frequency ${\omega }_p$, wave vector $k$, and Rabi frequency ${\Omega }_p$. The strength of the cavity coupling is given by $g_0$. The frequencies of the cavity and the atomic transition are denoted by ${\omega }_c$ and ${\omega }_a$, respectively. Here $|k_x,k_z\rangle$ are momentum states of the atoms in the BEC, and excited electronic states are denoted by a prime. The atomic ground state with $|k_x,k_z\rangle =|0,0\rangle$ and the symmetric superposition labeled as $|\pm k,\pm k\rangle$ constitute an effective two-level system governed by an effective nonequilibrium Dicke model. The corresponding level splitting is given by ${\omega }_0=2{\omega }_r$, where ${\omega }_r={\hslash }^2{k}^2/2m$ is the atomic recoil energy resulting from the absorption or emission of a single photon. In general, multiphoton processes are required in order to couple the ground state to higher-momentum states.

Figure 3

Schematic illustration of the different types of behavior displayed by the semiclassical equations of motion in Eq. (3), for trajectories on the Bloch sphere with $\left|\mathbf{S}\right|=N/2$. (a) Evolution from an unstable fixed point (open circle) to a stable attractor (closed circle) of the long-time dynamics for all initial conditions. Both the stable and unstable points have $\dot{\mathbf{S}}=\mathbf{0}$ and $\dot{\psi }=0$. However, for the stable attractor small perturbations decay, while for the unstable fixed point fluctuations grow. (b) As for (a) but with a hyperbolic fixed point (cross) having one stable and one unstable eigenmode. Paths first approach and then leave the vicinity of this hyperbolic point, before eventually reaching a stable attractor with $\dot{\mathbf{S}}=\mathbf{0}$ and $\dot{\psi }=0$. (c) Dynamics exhibiting a stable limit cycle with $\dot{\mathbf{S}}\ne \mathbf{0}$ and $\dot{\psi }\ne 0$ for all initial conditions.

Figure 4

Dynamical phase diagram of the stable attractors as a function of the cavity frequency $\omega$ and the coupling $g=g^{\prime }$ with parameters $\kappa =8.1$ MHz and ${\omega }_0=0.05$ MHz taken from Ref. [22]. The panels represent different values of the feedback term $U$, going from $U=0$ (top) to $UN=-40$ MHz (bottom). The second panel corresponds to the experimental parameters used in Ref. [22]. Points (a)–(f) marked in the bottom panel correspond to the fixed-point illustrations shown in Fig. 6. The characteristic time evolution at points (b), (f), and (g) is given in Figs. 7 and 8.

Figure 5

Vertical slice through the second panel of Fig. 4 with $UN=-10$ MHz, corresponding to the experimental parameters used in Ref. [22]. Variation of (a) ${\left|\psi \right|}^2$, (c) $\mathrm{Arg}\left(\psi \right)$, and (e) $S_z$ on passing from the SRA to the SRB phase. In the vicinity of both of these transitions there is a narrow region of bistability denoted as 2SRA, where both of the SRA solutions given in Eq. (10) coexist; see (b), (d), and (f) for magnified images of the highlighted regions.

Figure 6

Bloch spheres with $\left|\mathbf{S}\right|=N/2$ corresponding to the points (a)–(f) in Fig. 4 showing the fixed points where $\dot{\mathbf{S}}=\mathbf{0}$ and $\dot{\psi }=0$. We distinguish between stable fixed points (filled circles), unstable fixed points (open circles), and hyperbolic points with one stable and one unstable eigenmode (crosses). The presence of more than one stable attractor on the Bloch sphere corresponds to a coexistence phase in Fig. 4. The trajectories show the typical time evolution, where the initial conditions are chosen in order to illustrate the attractors. In particular, the time scale for approaching the SRA fixed point in (b) is significantly longer than the time scale for approaching the SRB fixed point in (f). In the language of dynamical systems the transition from (a) to (b) is a pitchfork bifurcation where a single mode goes unstable. The transition from (c) to (a) is a subcritical Hopf bifurcation where two modes go unstable simultaneously. The transition from (c) to (d) involves the appearance of eight additional fixed points; these correspond to two stable and two unstable SRB fixed points and four hyperbolic SRA fixed points. The transitions from (d) to (e) and (e) to (f) are inverse pitchfork bifurcations in which a pair of hyperbolic SRA fixed points coalesce at a previously stable fixed point.

Figure 7

Time evolution in the superradiant regime with initial conditions that are close to the normal state ($\Downarrow$). The top two panels (and insets) correspond to point (b) in Figs. 4 and 6. The bottom two panels correspond to point (f). In each case, the evolution of the semiclassical equations for a single initial condition with $S_x=S_y=\sqrt{N}$ is shown in gray (marked SC). The average evolution for a Wigner-distributed ensemble of initial conditions is shown in black. The insets in the top two panels show magnified images of the highlighted regions.

Figure 8

Time evolution starting close to the normal state ($\Downarrow$) for the parameters at point (g) in Fig. 4. As in Fig. 7, the semiclassical dynamics for both a single starting point with $S_x=S_y=\sqrt{N}$ (marked SC) and a Wigner-distributed set of initial conditions are shown. The latter is indicated in black, and the semiclassical trajectory is shaded to match the Bloch sphere shown as an inset; for these initial conditions the trajectory almost covers the Bloch sphere.

Figure 9

Characteristic time scales of the semiclassical dynamics. We use the same parameters as in Fig. 4 with $UN=-10$ MHz corresponding to the experiments of Ref. [22]. (a) shows the time required for the instability of the initial ground state $(\Downarrow$) to develop; in regions where this state is stable this time scale is taken to $\infty$. (b) shows the asymptotic time scale for approaching the final stable state. In the regions of multistability, the rate of attraction to the final state for the given initial conditions ($\Downarrow$) is shown. In both panels the time scales show a strong variation throughout the phase diagram, with notable implications for finite-duration experiments. A vertical slice along the red dashed line is given in Fig. 10.

Figure 10

Vertical slice through Fig. 9 with $g\sqrt{N}=0.6$ MHz. (a) Initial growth time obtained by linear stability analysis around the normal ($\Downarrow$) state. The open circles show the exact semiclassical results obtained numerically from the quartic Eq. (16). The solid line corresponds to the approximation given in Eq. (18) where the plus sign corresponds to the unstable mode. (b) Asymptotic decay time obtained by linear stability analysis around the appropriate final asymptotic state. The open circles correspond to the eigenvalue of $|\eta \tilde{\mathbb{I}}-\tilde{\mathbb{M}}|=0$ with the largest imaginary part, where $\tilde{\mathbb{M}}$ is given by Eq. (C5). The solid gold (mid gray) line corresponds to the imaginary part of Eq. (20) and the solid yellow (light gray) line corresponds to the exact result in Eq. (22). In both panels, the shaded regions indicate where the normal state ($\Downarrow$) is stable. Both time scales show a strong dependence on the parameters.

Figure 11

Photon intensity maps showing ${\left|\psi \right|}^2$ after two different time intervals, with initial conditions that are close to the normal state ($\Downarrow$) with $S_x=S_y=\sqrt{N}$. We use the same parameters as in the second panel of Fig. 4 with $UN=-10$ MHz, corresponding to the experiments of Ref. [22]. (a) Intensity map obtained in the final asymptotic state with $t\to \infty$, showing the distinct regions of SRA, SRB, and $\Uparrow$ states. (b) Intensity map obtained after 10 ms showing good qualitative agreement with the asymptotic attractors in (a), but with a slower approach in the SRA regions. A cross section of (a) with $g\sqrt{N}=1.0$ MHz is provided in Fig. 5.

Figure 12

${\left|\psi \right|}^2$ found by increasing $g^2\propto t$ and recording the value achieved as a function of time. We use the same parameters as in Fig. 4 with $UN=-10$ MHz, corresponding to the experiments of Ref. [22]. The sweep is chosen so that $g^{2}N=(t/t_0)\times 2.5{\text{MHz}}^2$, where $t_0=10$ ms in (a) and (b) and $t_0=200$ ms in (c). The top two panels correspond to the experimental sweep duration used in Ref. [22]. (a) Comparison of the steady-state value of ${\left|\psi \right|}^2$ with that obtained by a semiclassical evolution (marked SC) and Wigner-distributed initial conditions, for the value of $g\sqrt{N}$ reached at a given time with $\omega =40$ MHz. (b) Photon intensity map obtained after a 10 ms sweep. (c) Photon intensity map obtained after a 200 ms sweep. In comparing these figures to Fig. 5 of Ref. [22], one should note that the vertical scale on our intensity plots is the photon frequency $\omega$ appearing in Eq. (2). In comparison, the authors of Ref. [22] use the detuning of the pump from the bare-cavity frequency as the vertical scale, hence the inverted and shifted axis.

Figure 13

Dynamical phase diagram as a function of $\omega$ and $g=g^{\prime }$ for $U⩾0$. The phase boundaries at ${\omega }_{\pm }=0$ separate with increasing $U$, and a regime of persistent oscillations emerges. The dynamics in this regime is shown in Fig. 14. In contrast to Fig. 4, the SRB phase is absent here since $U⩾0$.

Figure 14

Persistent oscillations at $\omega =10$ MHz, $UN=+40$ MHz, and $g\sqrt{N}=1$ MHz, starting close to the normal state ($\Downarrow$) with $S_x=S_y=\sqrt{N}$. The upper panel shows the attraction toward persistent oscillations, illustrating the transient behavior at short times and the persistent oscillations at later times. The inset shows the same data on the Bloch sphere using the same shading scheme. For all times after $\sim$15 ms, the trajectory on the Bloch sphere is restricted to a circle at constant polar angle. The lower panel shows the time dependence of ${\left|\psi \right|}^2$ and $\mathbf{S}$ in the persistent-oscillation regime.

Figure 15

Dynamical phase diagram at $UN=-40$ MHz as a function of $\delta g\equiv g^{\prime }-g$ and $\omega$ for a number of values of $\overline{g}=\frac{1}{2}(g+g^{\prime })$. Vertical dashed lines indicate cuts shown in Fig. 18.

Figure 16

Persistent oscillations for $UN=-40$ MHz, $\omega =40$ MHz, $g\sqrt{N}=0.998$ MHz, and $g^{\prime }\sqrt{N}=1.002$ MHz, corresponding to $\delta g/\overline{g}=0.004$ or $g^{\prime }/g=1.004$. The panels mirror those in Fig. 14. The shading in the top panel match those in the inset, highlighting the initial and intermediate trajectories, and the final persistent oscillations.

Figure 17

Magnified portion of the bottom panel of Fig. 4 in the vicinity of the tricritical point where three phase boundaries cross. In addition to the phases visible in Fig. 4, there is a narrow region denoted as 2SRA, where the two distinct SRA solutions given by Eq. (10) coexist; see inset.

Figure 18

Evolution of the phase diagram shown in the bottom panel of Fig. 4 as one varies the ratio $g^{\prime }/g$, increasing from 1.000 (top) to 1.004 (bottom). We use the same phase labeling conventions as in Fig. 4. Vertical dashed lines indicate cuts shown in Fig. 15.

Figure 19

Analog of the phase diagram in Fig. 15 for $UN=+40$ MHz. We use the same phase labeling conventions as in Fig. 4. The panels show the dependence on $\delta g\equiv g^{\prime }-g$ for fixed $\overline{g}=\frac{1}{2}(g+g^{\prime })$. The white region corresponds to a regime of persistent oscillations. Vertical dashed lines indicate cuts shown in Fig. 20.

Figure 20

Analog of the phase diagram shown in Fig. 18 for $UN=+40$ MHz. We use the same phase labeling conventions as in Fig. 4. Vertical dashed lines indicate cuts shown in Fig. 19.