# Lasing and counter-lasing phase transitions in a cavity QED system

Kevin C. Stitely<EMAIL_ADDRESS>Dodd-Walls Centre for Photonic
and Quantum Technologies, New Zealand  Department of Physics, University of
Auckland, New Zealand  Department of Mathematics, University of Auckland, New
Zealand  Andrus Giraldo  Dodd-Walls Centre for Photonic and Quantum
Technologies, New Zealand  Department of Mathematics, University of Auckland,
New Zealand  Bernd Krauskopf  Dodd-Walls Centre for Photonic and Quantum
Technologies, New Zealand  Department of Mathematics, University of Auckland,
New Zealand  Scott Parkins  Dodd-Walls Centre for Photonic and Quantum
Technologies, New Zealand  Department of Physics, University of Auckland, New
Zealand

###### Abstract

We study the effect of spontaneous emission and incoherent atomic pumping on
the nonlinear semiclassical dynamics of the unbalanced Dicke model – a
generalization of the Dicke model that features independent coupling strengths
for the co- and counter-rotating interaction terms. As well as the ubiquitous
superradiant behavior the Dicke model is well-known for, the addition of
spontaneous emission combined with the presence of strong counter-rotating
terms creates laser-like behavior termed _counter-lasing_. These states appear
in the semiclassical model as stable periodic orbits. We perform a
comprehensive dynamical analysis of the appearance of counter-lasing in the
unbalanced Dicke model subject to strong cavity dissipation, such that the
cavity field can be adiabatically eliminated to yield an effective Lipkin-
Meshkov-Glick (LMG) model. If the coupling strength of the co-rotating
interactions is small, then the counter-lasing phase appears via a Hopf
bifurcation of the de-excited state. We find that if the rate of spontaneous
emission is small, this can lead to resurgent superradiant pulses. However, if
the co-rotating coupling is larger, then the counter-lasing phase must emerge
via the steady-state superradiant phase. Such a transition is the result of
the competition of the coherent and incoherent processes that drive
superradiance and counter-lasing, respectively. We observe a surprisingly
complex transition between the two, associated with the formation of a chaotic
attractor over a thin transitional parameter region.

##  I Introduction

Studies of the dynamics associated with models of trapped ultracold atomic
gases and ion arrays have recently made tremendous progress towards both
experimentally realizing and theoretically modelling systems that present a
broad variety of complex dynamics. Increased technological and computational
capabilities have made experimental and numerical investigation of many-body
quantum systems now tractable [1, 2]. Apart from being of importance to
fundamental theoretical physics, such systems are of interest for their
applications in the fields of metrology, developing state-of-the-art atomic
clocks and other measurement devices [3, 4, 5], and quantum computation and
information processing [6, 7].

One such model that has been garnering much interest of late is the Dicke
model [8], which describes the interaction of an ensemble of atoms confined to
an optical cavity with a single electromagnetic field mode of that cavity. As
later shown by Hepp and Lieb, this model is now renowned for a quantum phase
transition to superradiance, where the emission of light by the atomic
ensemble is enhanced by maintaining atomic coherence [9, 10, 11]. Although
difficult to realize in practice with two-level atoms, it is possible to
engineer _effective_ Dicke models, which present superradiant dynamics in a
wide variety of arrangements. Examples are Bose–Einstein condensates [12, 13,
14, 15], trapped-ion arrays [16], and cavity quantum electrodynamics (QED)
with multilevel atoms driven by Raman transitions [17, 18, 19], which have
provoked much interest recently in studying their dynamics [20, 21, 22, 23,
24, 25, 26].

Figure 1:  Schematic of the atomic energy level configuration with the two
ground states $\ket{\downarrow}$ and $\ket{\uparrow}$, and the two excited
states $\ket{-}$ and $\ket{+}$. Two lasers are driving two Raman transitions
between the ground and excited states, which give rise to two coupling
strengths $\lambda_{-}$ and $\lambda_{+}$, corresponding to the co-rotating
and counter-rotating interactions in the Hamiltonian (II.1), respectively. In
addition, spontaneous emission and incoherent pumping between the two ground
states may occur at rates $\Gamma_{\downarrow}$ and $\Gamma_{\uparrow}$.

We focus here particularly on the cavity QED realizations of the Dicke model.
An effective Dicke model can be produced in this case by applying lasers
stimulating two Raman transitions between two atomic ground states of an
ensemble of four-level atoms. The scheme is illustrated in Fig. 1. There are
two ground states, $\ket{\downarrow}$ and $\ket{\uparrow}$, and two highly
detuned excited states, $\ket{-}$ and $\ket{+}$. Two lasers drive the
transitions $\ket{\downarrow}\leftrightarrow\ket{+}$ and
$\ket{\uparrow}\leftrightarrow\ket{-}$, while the cavity mode couples the
excited and ground states through the transitions
$\ket{\downarrow}\leftrightarrow\ket{-}$ and
$\ket{\uparrow}\leftrightarrow\ket{+}$. If the excited state detunings are
large, then the two atomic excited states can be adiabatically eliminated such
that only the dynamics between the two ground states need to be considered.
This model must be considered as an open quantum system, as photon loss
through cavity decay is unavoidable. The resulting superradiant quantum phase
transition is mediated by the atom-light coupling strength. If the coupling is
sufficiently high, the system transitions from the so-called normal phase,
where the atoms are in their (nominal) ground states $\ket{\downarrow}$ and
the light field is in the vacuum state, to the superradiant phase, where the
atoms and the cavity feature macroscopic excitation in the steady state.

This effective realization presents an interesting generalization of the Dicke
model if the lasers driving the two Raman transitions are not of the same
strength or frequency. In the theoretical model, the possibility of having two
driving lasers with different frequencies and electric field amplitudes leads
to different coupling strengths for the co-rotating and counter-rotating terms
of the interaction Hamiltonian, leading to the _unbalanced Dicke model_ – also
referred to in the literature as the anisotropic Dicke model [27], the
generalized Jaynes-Cummings-Rabi model [28], and the interpolating Dicke-
Tavis-Cummings model [25]. This generalization provides a surprisingly diverse
range of dynamical behavior, including normal-superradiant phase coexistence
[23, 25], hysteresis [29], superradiant oscillations [18], and chaos [26].

Among this behavior is the _pole-flip transition_ [26], which describes a
transition between the two types of atomic configurations in the normal phase:
either all spins down or all spins up. This transition was observed
experimentally in Ref. [18]. However, results indicated that the transition
threshold was not well resolved compared to the expectations provided by
theory. A subsequent theoretical analysis using a fermionic Keldysh path
integral formulation [30] found that the discrepancy can be explained by
considering atomic decay channels that break the conservation of the
collective spin length. The inclusion of atomic dissipation destabilizes the
inverted normal phase (where all spins are up) and can produce a type of
lasing state termed _counter-lasing_ [30, 11], which is “fed” by the counter-
rotating interactions; this is in contrast to usual lasing which is fed by co-
rotating interactions. Other effects brought about by the inclusion of both
individual and collective atomic dissipation have been investigated in Refs.
[31, 32, 33].

With the addition of an incoherent pump exciting the spins, regular lasing
states can be achieved as well [34]. It was found that, in a semiclassical
treatment, the establishment of the lasing phase occurs via a Hopf bifurcation
leading to the appearance of a stable periodic orbit.

In this work, we take a bifurcation theory approach to studying the emergence
of lasing and counter-lasing within a semiclassical treatment. We first study
the emergence of lasing and counter-lasing as a modification of the pole-flip
transition, which occurs in the absence of atomic dissipation or pumping. We
find that the transition is most naturally described in the bad-cavity or
dissipative limit, where the cavity field mode can be adiabatically eliminated
to yield a spin-only description of the relevant dynamics. For this model, we
perform a detailed analysis of the periodic orbits describing lasing and
counter-lasing states, resulting from Hopf bifurcations. We then investigate
the transition from superradiance to the counter-lasing phase, which we find
occurs via a small transitional parameter region across which infinitely many
(global) bifurcations occur, leading to the emergence of chaotic dynamics.

In the counter-lasing phase the spontaneous emission rate acts as a _pump_ for
the lasing process, in conjunction with the coherent counter-rotating
interactions between the atomic ensemble and the light field. The superradiant
phase is driven by both of the coherent co- and counter-rotating interactions,
while the counter-lasing phase is seeded by incoherent atomic dissipation and
the counter-rotating interactions. As such, the transition from superradiance
to the counter-lasing phase allows us to explore the competition between
coherent and incoherent processes of the model. This interplay results in a
remarkably complex transition, which we probe in detail by analyzing the
symbolic dynamics of the model via a topological quantity called the kneading
invariant.

The paper is organized as follows. In Sect. II we introduce the model and its
semiclassical approximation which we will study throughout the rest of the
paper. In Sect. III we study the transition from the normal phase to lasing
and counter-lasing states in the “strong dissipative limit” where the cavity
decay rate is large in comparision to the effective cavity resonance
frequency. In Sect. IV we consider the transition from superradiance to the
counter-lasing phase, which reveals an abundance of complex chaotic dynamics
within a small transitional region of parameter space. Here, we present a
parameter sweeping of a kneading invariant to showcase the intricate
organization of global bifurcations. Our conclusions are presented in Sect. V.

##  II Model

We now introduce the unbalanced Dicke model with photon loss, atomic decay and
pumping. We consider the model in the dissipative limit for large photon decay
rates, where the field mode can be adiabatically eliminated to yield another
effective Hamiltonian. We then introduce the semiclassical limit, given by a
system of nonlinear ordinary differential equations.

###  II.1 Quantum Mechanical Model

The unbalanced Dicke model considers an ensemble of $N$ two-level atoms
interacting with a single mode of the electromagnetic field of an optical
cavity. The Hamiltonian (with $\hbar=1$) is

$\displaystyle\hat{H}=\ $
$\displaystyle\omega\hat{a}^{\dagger}\hat{a}+\omega_{0}\hat{J}_{z}+\frac{\lambda_{-}}{\sqrt{N}}\left(\hat{a}\hat{J}_{+}+\hat{a}^{\dagger}\hat{J}_{-}\right)$
$\displaystyle+\frac{\lambda_{+}}{\sqrt{N}}\left(\hat{a}\hat{J}_{-}+\hat{a}^{\dagger}\hat{J}_{+}\right),$
(1)

where $\hat{a}$ is the annihilation operator of the cavity field mode, and

$\hat{J}_{\pm,x,y,z}=\sum_{j=1}^{N}\hat{\sigma}^{(j)}_{\pm,x,y,z}$ (2)

are collective angular momentum operators; here
$\hat{\sigma}^{(j)}_{\pm,x,y,z}$ are the spin-$\frac{1}{2}$ Pauli operators
for the $j$th atom. These operators obey the commutation relations
$[\hat{\sigma}^{(j)}_{m},\hat{\sigma}^{(j)}_{n}]=i\varepsilon_{mn\ell}\hat{\sigma}^{(j)}_{\ell}$
(following the Einstein summation convention), where the indices $m$, $n$,
$\ell$ take on the values $x$, $y$, $z$ and $\varepsilon_{mn\ell}$ is the
Levi-Civita symbol. Here, $\hat{\sigma}^{(j)}_{\pm}=\hat{\sigma}^{(j)}_{x}\pm
i\hat{\sigma}^{(j)}_{y}$ are the angular momentum raising and lowering
operators. The effective frequency of the cavity mode is $\omega$, and the
effective frequency splitting of the atomic energy levels is $\omega_{0}$. The
co-rotating and counter-rotating coupling strengths are $\lambda_{-}$ and
$\lambda_{+}$, respectively.

With the addition of cavity decay, and incoherent (individual) atomic
dissipation and pumping, the quantum dynamics are governed by the master
equation

$\frac{d\hat{\rho}}{dt}=-i[\hat{H},\hat{\rho}]+\kappa\mathcal{D}\\{\hat{a}\\}\hat{\rho}+\sum_{j=1}^{N}\Gamma_{\downarrow}\mathcal{D}\\{\hat{\sigma}^{(j)}_{-}\\}\hat{\rho}+\Gamma_{\uparrow}\mathcal{D}\\{\hat{\sigma}^{(j)}_{+}\\}\hat{\rho},$
(3)

where $\hat{\rho}$ is the density operator, $\kappa$ is the cavity field decay
rate and $\Gamma_{\downarrow}$ and $\Gamma_{\uparrow}$ are the atomic
dissipation and pumping rates, respectively (see Fig. 1). The operator

$\mathcal{D}\\{\hat{X}\\}\hat{\rho}=2\hat{X}\hat{\rho}\hat{X}^{\dagger}-\hat{X}^{\dagger}\hat{X}\hat{\rho}-\hat{\rho}\hat{X}^{\dagger}\hat{X}$
(4)

is the standard Lindblad dissipator for a jump operator $\hat{X}$. Equation
(3) is invariant under the exchange

$\hat{a}\leftrightarrow-\hat{a},\
\hat{a}^{\dagger}\leftrightarrow-\hat{a}^{\dagger},\
\hat{\sigma}_{-}^{(j)}\leftrightarrow-\hat{\sigma}_{-}^{(j)},\
\hat{\sigma}_{+}^{(j)}\leftrightarrow-\hat{\sigma}_{+}^{(j)},$ (5)

which gives rise to a $\mathbb{Z}_{2}$ (parity) symmetry.

In either the dissipative limit $\kappa\gg\omega,\omega_{0},\lambda_{\pm}$
[35] or the dispersive limit $\omega\gg\kappa,\omega_{0},\lambda_{\pm}$ [36,
37] excitations in the cavity mode are short-lived. In both limits, the
timescales for the atomic and photonic dynamics can be separated and the field
mode adiabatically eliminated. That is, the cavity field is treated as
mediating (long-range) atomic interactions. The resulting Hamiltonian

$\hat{H}_{\mathrm{LMG}}=\omega_{0}^{\prime}\hat{J}_{z}-\frac{\xi}{N}\left[(\lambda_{-}+\lambda_{+})^{2}\hat{J}_{x}^{2}+(\lambda_{-}-\lambda_{+})^{2}\hat{J}_{y}^{2}\right]$
(6)

is a generalized form of a Lipkin-Meshkov-Glick (LMG) model [38] describing
the spin interactions. Including dissipative and incoherent pumping terms
introduces the master equation

$\displaystyle\frac{d\hat{\rho}}{dt}=\ $
$\displaystyle-i[\hat{H}_{\mathrm{LMG}},\hat{\rho}]+\frac{\eta}{N}\mathcal{D}\\{\lambda_{-}\hat{J}_{-}+\lambda_{+}\hat{J}_{+}\\}\hat{\rho}$
$\displaystyle+\sum_{j=1}^{N}\Gamma_{\downarrow}\mathcal{D}\\{\hat{\sigma}_{-}^{(j)}\\}\hat{\rho}+\Gamma_{\uparrow}\mathcal{D}\\{\hat{\sigma}_{+}^{(j)}\\}\hat{\rho},$
(7)

where

$\omega_{0}^{\prime}=\omega_{0}-\frac{\xi}{N}\left(\lambda_{-}^{2}-\lambda_{+}^{2}\right),\
\xi=\frac{\omega}{\kappa^{2}+\omega^{2}},\
\eta=\frac{\kappa}{\kappa^{2}+\omega^{2}}.$ (8)

Equation (II.1) constitutes the quantum mechanical generalized LMG model that
we consider in this paper.

###  II.2 Semiclassical Model

We take a semiclassical approach valid for large numbers of atoms
$N\rightarrow\infty$, where the expectations of operator products factorize,
that is, quantum fluctuations are negligible. Moreover, the atoms couple
uniformly to the cavity mode and feature identical atomic decay and pumping
rates, so that the density operator site-decouples as
$\hat{\rho}=\bigotimes_{j=1}^{N}\hat{\rho}_{j}$, where $\hat{\rho}_{j}$ is the
density operator for the $j$th atom [31, 33]. Hence, the dynamics described by
the master equation (II.1) can be represented as a closed system of nonlinear
differential equations. Defining the variables

$\beta=\frac{\braket{\hat{J}_{-}}}{N}=b_{x}-ib_{y}\in\mathbb{C},\
\gamma=\frac{\braket{\hat{J}_{z}}}{N}\in\mathbb{R},$ (9)

the semiclassical equations of motion are then

$\displaystyle\frac{d\beta}{dt}={}$
$\displaystyle-i\omega_{0}\beta-2\eta(\lambda_{+}^{2}-\lambda_{-}^{2})\beta\gamma-2i\xi(\lambda_{+}^{2}+\lambda_{-}^{2})\beta\gamma$
$\displaystyle-4i\xi\lambda_{-}\lambda_{+}\beta^{*}\gamma-(\Gamma_{\uparrow}+\Gamma_{\downarrow})\beta,$
(10) $\displaystyle\frac{d\gamma}{dt}={}$ $\displaystyle
2\eta(\lambda_{+}^{2}-\lambda_{-}^{2})\beta^{*}\beta+2i\xi\lambda_{-}\lambda_{+}[(\beta^{*})^{2}-\beta^{2}]$
$\displaystyle-2\Gamma_{\downarrow}\left(\frac{1}{2}+\gamma\right)+2\Gamma_{\uparrow}\left(\frac{1}{2}-\gamma\right),$
(11)

which constitute the semiclassical generalized LMG model that we consider
throughout the paper. The nonlinear dynamics described by these equations is
the central focus of this paper. They form a dynamical system with a three-
dimensional phase space $(b_{x},b_{y},\gamma)\in\mathbb{R}^{3}$ since $\beta$
is complex. As is the case for Eq. (3), Eqs. (10)–(11) feature
$\mathbb{Z}_{2}$ (parity) symmetry given by the equivariance [39] under the
transformation

$\mathcal{T}:\beta\mapsto-\beta,$ (12)

which corresponds to a rotation by an angle $\pi$ about the $\gamma$-axis in
$(b_{x},b_{y},\gamma)$-space. This symmetry produces trajectories that either
come in pairs, with the symmetric counterpart related by $\mathcal{T}$, or are
themselves symmetric under $\mathcal{T}$.

When the atomic pumping and damping terms vanish (i.e.,
$\Gamma_{\uparrow}=\Gamma_{\downarrow}=0$) angular momentum is conserved and
trajectories are confined onto the Bloch sphere defined by

$b_{x}^{2}+b_{y}^{2}+\gamma^{2}=R^{2}.$ (13)

Here $0<R\leq 1/2$ is the radius of the Bloch sphere, determined by the
initial atomic polarization. The maximal atomic polarization $R=1/2$ is
determined by the total spin length $N/2$ of $N$ spin-$\frac{1}{2}$ atoms.
Hence, the phase space of Eqs. (10)–(11) in the absence of atomic dissipation
($\Gamma_{\downarrow}=0$) and pumping ($\Gamma_{\uparrow}=0$) is the sphere
$\mathbb{S}^{2}(R)$ of radius $R$.

##  III Strong Dissipative Limit

Analysis of the semiclassical dynamics of the unbalanced Dicke model without
atomic dissipation or pumping indicates two distinct types of pole-flip
transition: the normal and superradiant cases [26]. The normal pole-flip
transition describes the transition from the spin-down normal phase to the
inverted (spin-up) normal phase, while the superradiant pole-flip transition
describes the transition from superradiance to the inverted normal phase. In
this section we focus on the simpler normal pole-flip transition. To this end,
we consider the _strong dissipative limit_ $\xi=0$ ($\kappa\gg\omega$).

In this regime, Eqs. (10)–(11) are most conveniently represented in spherical
coordinates

$\displaystyle b_{x}$ $\displaystyle=r\sin(\theta)\cos(\phi),$ (14)
$\displaystyle b_{y}$ $\displaystyle=r\sin(\theta)\sin(\phi),$ (15)
$\displaystyle\gamma$ $\displaystyle=r\cos(\theta),$ (16)

with radius $r\geq 0$, polar angle $\theta\in[0,\pi]$ (measured from the
positive $\gamma$-axis), and azimuthal angle $\phi\in[0,2\pi)$ (measured from
the positive $b_{x}$-axis). This leads to the set of ordinary differential
equations:

$\displaystyle\frac{dr}{dt}$ $\displaystyle=\delta\cos(\theta)-\sigma
r\cos^{2}(\theta)-\sigma r,$ (17) $\displaystyle\frac{d\theta}{dt}$
$\displaystyle=-\left[\mu
r-\sigma\cos(\theta)+\frac{\delta}{r}\right]\sin(\theta),$ (18)
$\displaystyle\frac{d\phi}{dt}$ $\displaystyle=\omega_{0},$ (19)

where we have introduced the parameters

$\mu=2\eta(\lambda_{+}^{2}-\lambda_{-}^{2}),\ \
\sigma=\Gamma_{\uparrow}+\Gamma_{\downarrow},\ \
\delta=\Gamma_{\uparrow}-\Gamma_{\downarrow}.$ (20)

Figure 2: Normal pole-flip transition on the Bloch sphere in the generalized
semiclassical LMG model in the strong dissipative limit, given by Eqs.
(10)–(11). Panel (a1) shows a trajectory on the Bloch sphere in the spin-down
normal phase and panel (a2) the trajectory in the $(b_{x},\gamma)$-half-plane
of the rotating frame. Panel (b1) illustrates a selection of periodic orbits
that appear during the normal pole-flip transition for $\lambda_{+}=0.5$,
where in panel (b2) these periodic orbits appear as equilibria in the rotating
frame. In panel (c1) the system is now in the spin-up normal phase for
$\lambda_{+}=0.8$, and its rotating frame trajectory is shown in panel (c2).
Other parameters for all plots are $\kappa=5$, $\omega=0$, $\omega_{0}=1$,
$\lambda_{-}=0.5$, $\Gamma_{\uparrow}=\Gamma_{\downarrow}=0$.

In these spherical coordinates Eq. (19) for $\phi$ decouples from Eqs. (17)
and (18) for $r$ and $\theta$. Furthermore, $d\phi/dt$ is constant, meaning
that system (10)–(11) features a constant rotation about the $\gamma$-axis in
$(b_{x},b_{y},\gamma)$-space at the frequency $\omega_{0}$. Thus, in the
strong dissipative limit, the semiclassical equations of motion are invariant
under rotations of any angle about the $\gamma$-axis. That is, Eqs. (10)–(11)
have the continuous $U(1)\cong\mathbb{S}^{1}$ symmetry, given by the
equivariance under the map

$\mathcal{S}:\beta\mapsto e^{i\phi}\beta,$ (21)

which contains the discrete $\mathbb{Z}_{2}$ symmetry as the special case
$\phi=\pi$. This has the effect that the overall dynamics of Eqs. (10)–(11) is
given by the dynamics in the $(r,\theta)$-half-plane with $r\geq 0$, subject
to a constant rotation about the $\gamma$-axis at frequency $\omega_{0}$.

###  III.1 Normal Pole-Flip Transition

In the case of zero atomic damping or pumping,
$\Gamma_{\downarrow}=\Gamma_{\uparrow}=\sigma=\delta=0$, the dynamics become
even simpler. Here, angular momentum is now conserved [Eq. (13)], and the
equations of motion in spherical coordinates become

$\displaystyle\frac{dr}{dt}$ $\displaystyle=0,$ (22)
$\displaystyle\frac{d\theta}{dt}$ $\displaystyle=-\mu r\sin(\theta),$ (23)
$\displaystyle\frac{d\phi}{dt}$ $\displaystyle=\omega_{0}.$ (24)

Equation (22) means the dynamics are confined to invariant spherical surfaces
parameterized by $R$ with $0<R\leq 1/2$. Additionally, from Eq. (11) it
follows that the entire $\gamma$-axis is a one-dimensional _manifold of
equilibria_ , since every value of $\gamma$ for $\beta=0$ is an equilibrium
point. This manifold’s intersection with one of the spheres creates the two
(spin-down and spin-up) normal phase equilibria at the South ($\theta=\pi$)
and North ($\theta=0$) poles of the Bloch sphere, respectively. The dynamics
on the sphere are essentially one-dimensional, determined by Eq. (23).

Since $\sin(\theta)>0$ for $0\leq\theta\leq\pi$, the sign of $\mu$ then
determines which equilibrium is approached, and thus determines the phase. If
$\mu<0$ (i.e., $\lambda_{-}>\lambda_{+}$) then $\theta(t)\rightarrow\pi$ as
$t\rightarrow\infty$ and the system is in the spin-down normal phase with
$\gamma\rightarrow-\frac{1}{2}$. Hence, trajectories on the Bloch sphere
approach the South pole in a rotational fashion. This is illustrated on the
Bloch sphere in Fig. 2(a1) for an initial condition on the equator in order to
highlight the direction of the flow. It is illustrative here to consider the
dynamics without the constant frequency rotation of $\phi$, that is, consider
only Eqs. (17)–(18); the corresponding trajectory in the rotating frame in the
$(b_{x},\gamma)$-half-plane is shown in Fig. 2(a2). As the atoms become de-
excited during their traversal of the Bloch sphere, they rapidly release a
pulse of photons in a manner referred to as _transient_ or _Dicke_
superradiance [40].

In contrast, if $\mu>0$ (i.e., $\lambda_{+}>\lambda_{-}$) then
$\theta(t)\rightarrow 0$ as $t\rightarrow\infty$ and the system is in the
inverted spin-up normal phase; see Fig. 2(c1) and the flow in the rotating
frame in Fig. 2(c2). Here, trajectories approach the equilibrium point at the
North pole of the Bloch sphere.

The transition between these two cases for $\mu=0$ is the pole-flip
transition. The dynamics here are highly degenerate, as $d\theta/dt=0$. Then,
subject to the constant rotation of the angle $\phi$, the Bloch sphere (again
for a fixed radius $0<R\leq 1/2$) is foliated by infinitely many circular
periodic orbits, which lie in planes of constant $\gamma$. This is illustrated
in Fig. 2(b1), with the flow frozen in the rotating frame forming a circle of
equilibria in Fig. 2(b2). Physically, this means that these periodic orbits
conserve the semiclassical energy per atom
$E=\braket{\hat{H}_{\mathrm{LMG}}}/N=\omega_{0}\gamma$. That is, in the
(three-dimensional) LMG model description, after adiabatic elimination, the
entire phase space $\mathbb{S}^{2}$ is foliated by periodic orbits, meaning
that the conservation of energy is satisfied _globally_. In contrast, in the
full (four-dimensional) description of the Dicke model [26], the periodic
orbits foliate a two-dimensional (attracting) invariant surface within the
full phase space $\mathbb{R}^{2}\times\mathbb{S}^{2}$. As observed in Ref.
[26], the conservation of energy at the pole-flip transition in the Dicke
model is satisfied within the two-dimensional surface, but not globally. This
implies that the two-dimensional surface observed in the Dicke model at the
moment of the pole-flip transition is a signature of the semiclassical LMG
model behavior discussed above.

###  III.2 Effect of Atomic Dissipation and Pumping

The addition of a pump that excites the atomic states removes the possibility
for a normal phase with all spins down, and produces lasing behavior if the
co-rotating coupling strength $\lambda_{-}$ is sufficiently large. Similarly,
the inclusion of atomic dissipation through spontaneous emission removes the
possibility of an inverted normal phase and can seed counter-lasing behavior
with a strong counter-rotating coupling strength $\lambda_{+}$. We will now
show that these lasing phases appear in the semiclassical model as stable
periodic orbits arising from Hopf bifurcations.

####  III.2.1 Hopf Bifurcations to Counter-Lasing

With nonzero atomic pumping $\Gamma_{\uparrow}$ and dissipation
$\Gamma_{\downarrow}$ in the strongly dissipative model, angular momentum
conservation (13) is broken ($dr/dt\neq 0$); thus, freeing the dynamics from
(the surface of) the Bloch sphere and removing the inverted normal phase
equilibrium. Crucially, this means that system (10)–(11) now has a _single_
equilibrium at

$\beta_{\mathrm{eq}}=0,\
\gamma_{\mathrm{eq}}=\frac{1}{2}\frac{\Gamma_{\uparrow}-\Gamma_{\downarrow}}{\Gamma_{\uparrow}+\Gamma_{\downarrow}}=\frac{\delta}{2\sigma}.$
(25)

Note that if $\Gamma_{\uparrow}>\Gamma_{\downarrow}$, then
$\gamma_{\mathrm{eq}}$ is strictly positive, and it is strictly negative if
$\Gamma_{\downarrow}>\Gamma_{\uparrow}$.

Figure 3:  Hopf bifurcation curves in the $(\lambda_{-},\lambda_{+})$-plane.
Panel (a) shows a selection of the hyperbolae for the counter-lasing
transition when $\Gamma_{\uparrow}=0$; here
$\Gamma_{\downarrow}=0.01,0.05,0.1,0.15$. Panel (b) shows the case of the
lasing transition for $\Gamma_{\downarrow}=0$ for
$\Gamma_{\uparrow}=0.01,0.05,0.1,0.15$. Panels (a1) and (b1) show the
extension of the Hopf bifurcation curves for negative coupling strengths to
highlight the symmetry properties of the curves. Panels (a2) and (b2) are
enlargements and their insets show on which sides of the curves the counter-
lasing/lasing and normal phases reside. Other parameters are $\kappa=5$,
$\omega_{0}=1$.

The Jacobian matrix $J$ of system (10)–(11) in the strongly dissipative case
at the equilibrium $\beta_{\mathrm{eq}}=0$, $\gamma=\gamma_{\mathrm{eq}}$ has
the simple (block diagonal) structure

$J=\begin{bmatrix}-\mu\gamma_{\mathrm{eq}}-\sigma&-\omega_{0}&0\\\
\omega_{0}&-\mu\gamma_{\mathrm{eq}}-\sigma&0\\\ 0&0&-2\sigma\end{bmatrix},$
(26)

with eigenvalues

$\displaystyle v_{x,y}$ $\displaystyle=-\mu\gamma_{\mathrm{eq}}-\sigma\pm
i\omega_{0},$ (27) $\displaystyle v_{\gamma}$ $\displaystyle=-2\sigma.$ (28)

The eigenvectors associated with the eigenvalues $v_{x,y}$ span the
$\beta$-plane at a height $\gamma_{\mathrm{eq}}$. On the other hand, the
eigenvector associated with $v_{\gamma}$ spans the $\gamma$-axis. Equation
(27) for the eigenvalues $v_{x,y}$ indicates the presence of a Hopf
bifurcation when $v_{x,y}$ form a purely imaginary pair $\pm i\omega_{0}$,
that is, when

$\mu=\frac{-2\sigma^{2}}{\delta}.$ (29)

Or in terms of the LMG model parameters [Eqs. (8)],

$\lambda_{+}^{2}-\lambda_{-}^{2}=\frac{-(\Gamma_{\uparrow}+\Gamma_{\downarrow})^{2}}{\eta(\Gamma_{\uparrow}-\Gamma_{\downarrow})},$
(30)

which describes a family of hyperbolae in the
$(\lambda_{-},\lambda_{+})$-plane that asymptote onto the diagonal
$\lambda_{-}=\lambda_{+}$. The Hopf bifurcations along these curves create
circular periodic orbits, since system (17)–(19) is $U(1)$-symmetric, that is,
invariant under $\mathcal{S}$.

The periodic orbits formed from these Hopf bifurcations represent _lasing
states_ (L) for $\Gamma_{\uparrow}>\Gamma_{\downarrow}$, and _counter-lasing
states_ (CL) for $\Gamma_{\uparrow}<\Gamma_{\downarrow}$ [34]. The Hopf
bifurcation then marks the phase transition from the normal phase (N) to the
lasing/counter-lasing phase. A selection of these Hopf bifurcation curves are
shown in Fig. 3(a) for the counter-lasing phase transition for
$\Gamma_{\uparrow}=0$ and $\Gamma_{\downarrow}=0.01,0.05,0.1,0.15$, while Fig.
3(b) shows the bifurcation curves for the lasing transition for
$\Gamma_{\downarrow}=0$ and $\Gamma_{\uparrow}=0.01,0.05,0.1,0.15$. These
hyperbolae divide the $(\lambda_{-},\lambda_{+})$-plane into regions where the
counter-lasing (CL), lasing (L), and normal (N) phases may be found.
Importantly, the addition of atomic dissipation or pumping,
$\Gamma_{\uparrow},\Gamma_{\downarrow}>0$, not only creates the possibility of
lasing or counter-lasing but lifts the phase transition curve from the origin.
This means that, for the case $\Gamma_{\uparrow}=0$, there is now a nonzero
minimum threshold value of $\lambda_{+}$ for the transition to occur. This was
noticed experimentally in Ref. [18] where, above a certain threshold of the
ratio $\lambda_{+}/\lambda_{-}$ corresponding to the line of the pole-flip
transition, no phase transition is expected if atomic damping is neglected.
The minimum (so-called single-beam) threshold for the transition was also
investigated experimentally in Ref. [19]. This minimum threshold occurs when
either $\lambda_{-}$ or $\lambda_{+}$ is zero, depending on the sign of
$\delta=\Gamma_{\uparrow}-\Gamma_{\downarrow}$. In the semiclassical LMG
model, the minimum threshold occurs at

$\lambda_{\pm}=\frac{\Gamma_{\uparrow}+\Gamma_{\downarrow}}{\sqrt{\mp\eta(\Gamma_{\uparrow}-\Gamma_{\downarrow})}}.$
(31)

We remark that in terms of Eqs. (17)–(18) the Hopf bifurcation is manifested
as a _pitchfork bifurcation_ of relative equilibria. This is easily seen since
the normal phase equilibrium

$(r,\theta)=\begin{cases}(|\gamma_{\mathrm{eq}}|,0),\ \text{if }\delta>0,\\\
(|\gamma_{\mathrm{eq}}|,\pi),\ \text{if }\delta<0,\\\ (0,0),\ \text{if
}\delta=0,\ \sigma>0,\end{cases}$ (32)

always exists, and the relative equilibria, which represent periodic orbits,
bifurcate from this equilibrium.

Figure 4:  Convergence of trajectories towards counter-lasing periodic orbits
as $\lambda_{+}$ is varied, keeping $\Gamma_{\downarrow}=0.01$ constant. Panel
(a1) shows a trajectory initiated at $(r,\theta,\phi)=(1/2,\pi+0.001,0)$ in
$(b_{x},b_{y},\gamma)$-space for $\omega_{0}=1$. Panel (a2) shows the
trajectory in the rotating frame, where the periodic orbit appears as a
relative equilibrium point. Panel (a3) shows the temporal trace of $b_{x}$ of
the counter-lasing periodic orbit. Panel (b1) shows a trajectory approaching a
counter-lasing periodic orbit further from the Hopf bifurcation with
$\omega_{0}=1$, with the rotating frame trajectory shown in panel (b2). In
panel (a) $\lambda_{+}=0.6$ and in panel (b) $\lambda_{+}=1$. For both,
$\lambda_{-}=0.5$, $\kappa=5$, $\Gamma_{\uparrow}=0$. Figure 5: Surfaces of
bifurcating periodic orbits. Panel (a) shows two oblate spheroids of periodic
orbits formed after the Hopf bifurcation. On the top in red is the spheroid
formed from the lasing orbits for $\lambda_{+}=0.5$, $\Gamma_{\uparrow}=0.01$,
$\Gamma_{\downarrow}=0$, with increasing $\lambda_{-}$ down the surface
towards the origin, shown as a white dot. On the bottom in blue is the
spheroid formed from counter-lasing periodic orbits for $\lambda_{-}=0.5$,
$\Gamma_{\uparrow}=0$, $\Gamma_{\downarrow}=0.01$, with increasing
$\lambda_{+}$ up the surface towards the origin. Panel (b) shows the
projections of these surfaces of periodic orbits on the half-plane
$(b_{x},\gamma)$. Panel (c) shows $\max(b_{x,\mathrm{PO}})$ of the periodic
orbits for the counter-lasing orbits as $\lambda_{+}$ is increased. The dashed
line indicates where the normal phase equilibrium point is unstable; vertical
gray lines indicate the periodic orbits highlighted on the spheroids in panel
(a). The curve of $\max(b_{x,\mathrm{PO}})$ for the lasing orbits is the same
as for the counter-lasing orbits, but for increasing $\lambda_{-}$ rather than
$\lambda_{+}$. Other parameters are $\kappa=5$, $\omega=0$, $\omega_{0}=1$.

####  III.2.2 Transient Dynamics

The system’s convergence to a counter-lasing periodic orbit formed just after
the Hopf bifurcation is shown in Fig. 4(a1) by numerical integration of Eqs.
(17)–(19), with initial conditions perturbed from the (linearly unstable)
equilibrium point $(r,\theta)=(0.5,\pi)$. Due to the constant rotation by the
frequency $\omega_{0}$ about the $\gamma$-axis, the “density” of intersection
points on an (arbitrary) plane defined by $b_{y}=0$ gives information about
the speed of convergence of the trajectory onto the periodic orbit, since the
trajectory will intersect the plane every $T=2\pi/\omega_{0}$ units of time.
Thus, it can be seen that, just after the Hopf bifurcation, the spiralling
trajectory in Fig. 4(a1) approaches the periodic orbit slowly. The
corresponding trajectory in the rotating frame is shown in Fig. 4(a2),
highlighting the tightly-wound spiral structure of the flow, while Fig. 4(a3)
is a time series of the attracting periodic orbit. Further away from the Hopf
bifurcation, as $\lambda_{+}$ is increased, convergence to the periodic orbit
is initially faster in Fig. 4(b1) and the spiralling structure in Fig. 4(b2)
is more pronounced.

####  III.2.3 Geometry of Counter-Lasing Orbits

The geometry of the periodic orbits describing the lasing and counter-lasing
phases can be examined with Eqs. (17)–(19) as $\mu$ is varied. As mentioned
previously, the invariance under the transformation $\mathcal{S}$ ensures the
periodic orbits are circular and have period $T=2\pi/\omega_{0}$. A periodic
orbit with radius $r_{\mathrm{PO}}$ and polar angle $\theta_{\mathrm{PO}}$
must have $dr/dt=d\theta/dt=0$, that is, periodic orbits are described by
relative equilibria of the equations in the rotating frame, Eqs. (17)–(18).
Solving for them as common roots gives

$\cos(\theta_{\mathrm{PO}})=\frac{\sigma}{\sqrt{-(\delta\mu+\sigma^{2})}}.$
(33)

This equation sets the condition for relative equilibria in addition to Eq.
(32) to exist, since $-1\leq\cos(\theta_{\mathrm{PO}})\leq 1$. Either equality
is attained, for $\delta>0$ or $\delta<0$, respectively, when Eq. (29) is
satisfied. This describes a pitchfork bifurcation of Eqs. (17)–(18) as the
relative equilibria $(r_{\mathrm{PO}},\theta_{\mathrm{PO}})$ emerge from the
normal phase equilibrium with an exchange of stability. The radius
$r_{\mathrm{PO}}$ is

$r_{\mathrm{PO}}=\frac{\delta\cos(\theta_{\mathrm{PO}})}{\sigma\left[1+\cos^{2}(\theta_{\mathrm{PO}})\right]}=\frac{\sqrt{-(\delta\mu+\sigma^{2})}}{\mu},$
(34)

which also illustrates the Hopf bifurcation condition Eq. (29) when
$r_{\mathrm{PO}}=0$.

Figure 6:  Bifurcation diagram of periodic orbits formed from Hopf
bifurcations as $\lambda_{+}$ is increased. The blue colored curves on the
right are curves in the counter-lasing phase:
$\Gamma_{\downarrow}=0.001,0.01,0.05,0.1$, $\Gamma_{\uparrow}=0$. The orange
colored curves on the left are for the lasing phase: $\Gamma_{\downarrow}=0$,
$\Gamma_{\uparrow}=0.001,0.01,0.05,0.1$. Black dots indicate Hopf
bifurcations. Other parameters are $\kappa=5$, $\lambda_{-}=1$.

The first equality in Eq. (34) describes an _oblate spheroid_ centered at
$(b_{x},b_{y},\gamma)=(0,0,\gamma_{\mathrm{eq}}/2)$, an ellipsoid with equal
axes in planes of constant $\gamma$, parameterized by $\mu$. Two surfaces of
periodic orbits for the lasing (L) and counter-lasing (CL) cases are shown in
Fig. 5(a), for constant $\lambda_{+}$ and varying $\lambda_{-}$ (L surface)
and for constant $\lambda_{-}$ and varying $\lambda_{+}$ (CL surface), the
projection of the surface onto the half-plane $(b_{x},\gamma)$ is shown in
Fig. 5(b). For the surface L of lasing orbits, the Hopf bifurcation occurs at
the top where $\gamma=\gamma_{\mathrm{eq}}>0$; as $\lambda_{-}$ is increased
the orbits move down the surface towards the origin as
$\lambda_{-}\rightarrow\infty$. For the counter-lasing orbits, the Hopf
bifurcation occurs at the bottom of the surface CL where
$\gamma=\gamma_{\mathrm{eq}}<0$; as $\lambda_{+}$ is increased the orbits move
up the surface, again towards the origin as $\lambda_{+}\rightarrow\infty$.
The semi-minor and semi-major axes of the spheroid are

$a=\frac{\delta}{4\sigma}\ \text{and}\ b=\frac{\delta}{2\sqrt{2}\sigma},$ (35)

respectively. The poles along the $\gamma$-axis are
$\gamma=\gamma_{\mathrm{eq}}$, where the Hopf bifurcation occurs, and
$\gamma=0$, where the periodic orbits accumulate as $\mu\rightarrow\infty$. We
remark that even though the adiabatic approximation in the large coupling
limit is not strictly satisfied, the consideration of such limits give insight
into how the system behaves when the coherent processes strongly dominate.

From Eqs. (33)–(34), the family of periodic orbits emerging from the Hopf
bifurcation can be described analytically as $\mu$ varies. The periodic orbits
are characterized by the maxima (or minima) of either $b_{x}$ or $b_{y}$ (we
only need one due to $U(1)$ symmetry) and $\gamma$. Firstly, we find

$\max(b_{x,\mathrm{PO}})=r_{\mathrm{PO}}\sin(\theta_{\mathrm{PO}})=\frac{\sqrt{2\eta\delta(\lambda_{-}^{2}-\lambda_{+}^{2})-2\sigma^{2}}}{2\eta(\lambda_{+}^{2}-\lambda_{-}^{2})},$
(36)

which is illustrated, as $\lambda_{+}$ is varied, in Fig. 5(c) for the
counter-lasing transition. Here, we see a very fast initial growth with
$\lambda_{+}$ in the size of the periodic orbits (as a function of
$\lambda_{+}$) after the Hopf bifurcation, followed by a slow relaxation as
the periodic orbits shrink towards the origin as
$\lambda_{+}\rightarrow\infty$. This indicates that, once the Hopf bifurcation
is reached, the periodic orbits move along the surface of periodic orbits in
Fig. 5(a) very quickly upon variation of $\lambda_{+}$, before slowing down as
they approach the origin, at the top of the surface where $\gamma=0$. This
applies to the lasing transition as well, where
$\Gamma_{\uparrow}>\Gamma_{\downarrow}$ and the curves are now functions of
$\lambda_{-}$, with $\lambda_{+}$ held constant.

###  III.3 Atomic Dissipation and Pumping near the Pole-Flip Limit

Figure 7: Hopf bifurcation curves in the
$(\Gamma_{\downarrow},\Gamma_{\uparrow})$-plane. In panel (a), the green
colored curves on the right are the Hopf bifurcation curves that lead to the
creation of counter-lasing periodic orbits for $\lambda_{-}=0.5$,
$\lambda_{+}=0.75,1,1.5,2$, while the purple colored curves on the left lead
to the creation of lasing periodic orbits for
$\lambda_{+}=0.5,\lambda_{-}=0.75,1,1.5,2$. Panels (b) and (c) show how these
curves separate the $(\Gamma_{\downarrow},\Gamma_{\uparrow})$-plane into
regions where the normal (N), lasing (L), and counter-lasing (CL) phases
reside. All of the curves are for $\kappa=5$.

Information about how the system approaches the pole-flip transition in the
limit $\delta,\sigma\rightarrow 0$ can be obtained from the
$\gamma$-components of the periodic orbits,

$\gamma_{\mathrm{PO}}=r_{\mathrm{PO}}\cos(\theta_{\mathrm{PO}})=\frac{\sigma}{\mu}=\frac{\Gamma_{\uparrow}+\Gamma_{\downarrow}}{2\eta(\lambda_{+}^{2}-\lambda_{-}^{2})},$
(37)

which are shown in Fig. 6 as $\lambda_{+}$ is varied in the lasing phase for
$\Gamma_{\downarrow}=0$, $\Gamma_{\uparrow}=0.001,0.01,0.05,0.1$, and in the
counter-lasing phase for $\Gamma_{\uparrow}=0$,
$\Gamma_{\downarrow}=0.001,0.01,0.05,0.1$. The curves, which are hyperbolae,
in the upper left hand side correspond to the periodic orbits in the lasing
phase, where $\Gamma_{\uparrow}>\Gamma_{\downarrow}$ and
$\lambda_{-}>\lambda_{+}$. As $\lambda_{+}$ is increased these orbits shrink
and disappear at Hopf bifurcations, causing the disappearance of the lasing
phase and the appearance of the normal phase. On the other hand, the curves on
the bottom right-hand side correspond to periodic orbits in the counter-lasing
phase, where $\Gamma_{\downarrow}>\Gamma_{\uparrow}$ and
$\lambda_{+}>\lambda_{-}$. In this case, upon increase of $\lambda_{+}$ the
normal phase is destabilized and the counter-lasing phase emerges. In the
limit $\Gamma_{\downarrow},\Gamma_{\uparrow}\rightarrow 0$ all of these curves
accumulate onto the vertical dashed line at $\lambda_{-}=\lambda_{+}$, which
corresponds to the pole-flip transition, and indicates coexistent periodic
orbits for all values of $\gamma$ between -0.5 and 0.5 [see Fig. 2(b)].

The two-parameter bifurcation diagram in the
$(\Gamma_{\downarrow},\Gamma_{\uparrow})$-plane can be obtained from Eq. (30),
solved with the coupling strengths $\lambda_{-}$ and $\lambda_{+}$ kept
constant with $\Gamma_{\downarrow}$ and $\Gamma_{\uparrow}$ now allowed to
vary. Here the Hopf bifurcation curves form a family of parabolae in the
$(\Gamma_{\downarrow},\Gamma_{\uparrow})$-plane, which intersect the origin
tangentially to the diagonal $\Gamma_{\downarrow}=\Gamma_{\uparrow}$ and
intersect $\Gamma_{\uparrow}=0$ at $\Gamma_{\downarrow}=\mu$ in the counter-
lasing phase; and similarly intersect $\Gamma_{\downarrow}=0$ at
$\Gamma_{\uparrow}=-\mu$ in the lasing phase. A selection of these bifurcation
curves for various values of the coupling strengths are shown in Fig. 7(a).
The curves in the bottom right portion of the diagram are the Hopf bifurcation
curves that give rise to the periodic orbits in the counter-lasing phase,
likewise the curves in the top left portion create the lasing phase periodic
orbits. Along the (positive) diagonal
$\Gamma_{\downarrow}=\Gamma_{\uparrow}>0$ the normal phase equilibrium is
located at the origin and is always stable. The counter-lasing Hopf
bifurcation curves accumulate onto the diagonal in the limit
$\lambda_{+}\rightarrow\infty$, as do the lasing Hopf bifurcation curves for
$\lambda_{-}\rightarrow\infty$.

Figure 8:  Transient dynamics showing how the counter-lasing periodic orbits
are approached. Panels (a1) and (b1) show a simulation of Eqs. (17)–(19) in
atomic space $(b_{x},b_{y},\gamma)$ for $\omega_{0}=1$. In panels (a)
$\Gamma_{\downarrow}=0.01$, in panels (b) $\Gamma_{\downarrow}=0.001$. Panels
(a2) and (b2) show the trajectories in the rotating frame that describe the
convergence to the periodic orbit, corresponding to (a1) and (b1),
respectively. In panels (a3) and (b3) we show the temporal trace of the photon
number $|\alpha|^{2}$ for $\omega_{0}=0.2$. For both $\lambda_{-}=0.5$,
$\lambda_{+}=0.8,\kappa=5$, $\Gamma_{\uparrow}=0$.

The Hopf bifurcation curves approach the origin
$\Gamma_{\downarrow}=\Gamma_{\uparrow}=0$. At this point, the situation
described in Sect. III.1 is regained when angular momentum is conserved and
the $\gamma$-axis is a one-dimensional manifold of equilibria. The dynamics
here are highly degenerate. According to Eqs. (27) and (28), the eigenvalues
of these equilibria are $v_{x,y}=-\mu\gamma\pm i\omega_{0}$ and
$v_{\gamma}=0$, indicating a zero eigenvalue for all $\mu$. There are two
configurations this line of equilibria can take, depending on the sign of
$\mu$. For $\mu<0$, the equilibria with $\gamma<0$ are stable (in the
$(b_{x},b_{y})$-plane), while for $\gamma>0$ they are unstable. Here the
system is in the spin-down normal phase. The change in stability that takes
place at the origin $b_{x}=b_{y}=\gamma=0$ is a bifurcation parameterized not
by a system parameter, but by the _coordinate_ $\gamma$. This is an example of
a _bifurcation without parameters_ [41]. If $\mu>0$ the opposite configuration
occurs: equilibria with $\gamma<0$ are unstable, and those with $\gamma>0$ are
stable. Here the system is in the spin-up normal phase. The transition between
these two configurations occurs for $\mu=0$, when each equilibrium features
the eigenvalues $v_{x,y}=\pm i\omega_{0}$ and a zero eigenvalue
$v_{\gamma}=0$. Unlike the generic case of the zero-Hopf bifurcation [42],
which satisfies this eigenvalue configuration, the particular situation we are
considering leads to the creation of infinitely many coexisting periodic
orbits that foliate phase space, as we saw before in Fig. 2(b). This means
that the case we are observing in the semiclasscal LMG model corresponds to an
unfolding of a degenerate zero-Hopf bifurcation.

The bifurcation analysis of the emergence of the periodic orbits corresponding
to the lasing and counter-lasing phases above considers only the final states
as $t\rightarrow\infty$. However, in many situations (barring the addition of
laser fields tuned to transitions between the atomic ground states) the
spontaneous emission rates $\Gamma_{\downarrow}$ are much smaller than the
rates that determine the interactions $\lambda_{\pm}$ or the cavity decay rate
$\kappa$. As a result, there are two relevant timescales when the atomic
damping and pumping is small.

To examine how trajectories converge onto the periodic orbits in the counter-
lasing phase, we numerically integrate Eqs. (17)–(19) with an initial
condition perturbed from the (linearly unstable) normal phase equilibrium
$(r,\theta)=(0.5,\pi)$. The results are plotted in Fig. 8. In panel (a1) we
show the transient convergence to the counter-lasing periodic orbit for
$\Gamma_{\downarrow}=0.01$. The trajectory spirals radially inwards to the
periodic orbit, accompanied by angular rotation at frequency $\omega_{0}$.
Again, the spiralling motion is more clearly seen in the rotating frame in
panel (a2). This spiralling behavior of the atomic dynamics also appears in
the cavity field dynamics, as the cavity field amplitude
$\alpha=\braket{\hat{a}}/\sqrt{N}$ is directly related to $\beta$. The photon
number is

$|\alpha|^{2}=(\xi^{2}+\eta^{2})\left[(\lambda_{-}^{2}+\lambda_{+}^{2})\beta^{*}\beta+\lambda_{-}\lambda_{+}(\beta^{2}+{\beta^{*}}^{2})\right],$
(38)

which is shown in Fig. 8(a3) (with $\omega_{0}=0.2$ for clarity of
illustration). The inwards spiralling of the rotating frame trajectories
modulates the amplitude of oscillations of the photon number, which eventually
settles to a periodic oscillation in the counter-lasing phase. The photon
number in the rotating frame provides the amplitude envelope.

Figure 8(b) shows the full and rotating frame trajectories as well as the
temporal trace of the photon number when the rate of spontaneous emission is
even smaller at $\Gamma_{\downarrow}=0.001$. In Fig. 8(b1) the trajectory is
initially perturbed from the normal phase equilibrium and makes a rapid
transition closely following along the Bloch sphere of radius $1/2$, to its
North pole. Here the dynamics slow down dramatically and the trajectory slowly
descends along the $\gamma$-axis; this axis acts as a slow manifold, since the
timescale of motion in this region of phase space is given only by
$\Gamma_{\downarrow}$ and $\Gamma_{\uparrow}$. After the slow descent along
the $\gamma$-axis, the trajectory lifts off and makes another rapid traversal
along another Bloch sphere “remnant.” Again, the geometry of how the periodic
orbit is approached is more clear in the rotating frame trajectory illustrated
in Fig. 8(b2). There is a sequence of rapid movement up Bloch sphere remnants
as the periodic orbit is approached, each followed by slow descents along the
$\gamma$-axis. The contrast between the slow movements along the $\gamma$-axis
and the fast movements away is clearly demonstrated in the photon number,
shown in Fig. 8(b3). Each rapid transition along a Bloch sphere remnant is
associated with a pulse of photons, followed by a relaxation period. As a
result, the lift-off of the trajectory from the $\gamma$-axis produces
resurgent superradiant pulses. Each subsequent pulse has a larger width and a
smaller amplitude as the counter-lasing phase is approached.

In the parameter regime we are considering here, the light field in the full
Dicke model description adiabatically follows the motion of the spins, such
that each rapid traversal along each shell produces a pulse of photons as the
spins are excited. The first traversal reaches the highest level of
excitation, and so produces the highest amplitude light pulse. Moreover, each
pulse arrives sooner than the last, until there are no identifiable pulses
left when the periodic orbit is reached. It is the first of these photon
pulses that are observed experimentally in Ref. [18] to determine when the
system is in the inverted phase. Experiments typically occur on a timescale
such that only one pulse can be observed, due to atom loss or dephasing.
However, when the counter-lasing periodic orbit is large (closer to the Hopf
bifurcation), trajectories are not compressed along the slow manifold [see
Fig. 4(a2)], so that the separation of timescales that leads to the creation
of clearly separated photon pulses is absent. This leads to a lengthening in
time of the photon pulses and a sharp peak in the Fourier spectrum near the
phase transition threshold. A similar effect has been noticed experimentally
in Ref. [18].

##  IV Transition between counter-lasing and superradiance

We now consider the full Eqs. (10)–(11) for $\xi\neq 0$, where the system
loses its $U(1)$-symmetry and is only $\mathbb{Z}_{2}$ invariant. Hence, there
is now the possibility to spontaneously break the $\mathbb{Z}_{2}$ symmetry in
a superradiant phase transition upon variation of $\lambda_{+}$. In this
Section, we will show that for fixed coupling $\lambda_{-}$ and increasing
$\lambda_{+}$, the counter-lasing phase must eventually stabilize even in the
presence of superradiance. It is then imperative to understand how the
counter-lasing phase emerges out of superradiance.

###  IV.1 Emergence of Superradiance

The emergence of the superradiant phase corresponds to the appearance of
stable asymmetric equilibria, which we will refer to as _superradiant
equilibria_. Much like its Dicke model counterpart, the transition to
superradiance in the LMG model with unbalanced coupling can occur in two ways
when varying $\lambda_{+}$, depending on the values of the other parameters
(see Ref. [26] for details):

  1. Case (i):

A supercritical pitchfork bifurcation $\mathrm{P}_{\mathrm{spr}}$ leads to the
creation of two stable superradiant equilibria, whilst simultaneously turning
the normal phase equilibrium unstable.

  2. Case (ii):

A saddle-node bifurcation SN leads to the creation of two pairs of
superradiant equilibria, two of which are stable and two of which are
unstable. Subsequently, the two unstable superradiant equilibria collide with
the stable normal phase equilibrium in a subcritical pitchfork bifurcation
$\mathrm{P}_{\mathrm{sub}}$, where the unstable superradiant equilibria
disappear and the normal phase becomes unstable; this leaves the other two
superradiant equilibria as the only remaining stable objects.

Case (ii) leads to a parameter range of coexisting superradiant and normal
phases, in between the saddle-node and pitchfork bifurcations, with an
associated hysteresis loop and switching dynamics between the phases in the
presence of quantum fluctuations [29]. Such a hysteresis loop between the
coexistent normal and superradiant phases has recently been observed
experimentally in Ref. [43].

###  IV.2 Superradiant Pole-Flip Transition

In the absence of atomic dissipation and pumping, there is a similar
transition to that as described in Sec. III.1 where the inverted normal phase
arises from the superradiant phase via a foliation of periodic orbits on the
Bloch sphere [26]. This transition, called the _superradiant pole-flip
transition_ , is illustrated in Fig. 9. In Fig. 9(a), the system begins in the
superradiant phase, where the asymmetric equilibria are stable. At the
superradiant pole-flip transition point, the Bloch sphere is again foliated by
periodic orbits, but with a different foliation: there are now periodic orbits
established around each superradiant equilibrium point, as well as a pair of
_homoclinic orbits_ , trajectories that converges to the (saddle-type) normal
phase equilibrium forwards and backwards in time; see Fig. 9(b). As
$\lambda_{+}$ is increased, the superradiant equilibria become unstable and
the spin-up normal phase equilibrium (North pole of the Bloch sphere) becomes
stable; see Fig. 9(c).

Figure 9: Superradiant pole-flip transition on the Bloch sphere in the
generalized semiclassical LMG model. In panel (a) the system is in the
superradiant phase for $\lambda_{+}=1.45$. During the superradiant pole-flip
transition in panel (b), periodic orbits oscillating around the spin-up normal
phase equilibrium and superradiant equilibria foliate the Bloch sphere,
separated by a homoclinic orbit (shown in red) forming a figure eight through
the spin-down normal phase equilibrium for $\lambda_{+}=1.5$. In panel (c) the
system is in the spin-up normal phase for $\lambda_{+}=1.55$. In panels (a)
and (c) trajectories are initialized along the unstable eigenspace of the
normal phase equilibrium point. Other parameters are $\lambda_{-}=1.5$,
$\kappa=4$, $\omega=0.5$, $\omega_{0}=0.2$,
$\Gamma_{\uparrow}=\Gamma_{\downarrow}=0$.

###  IV.3 Inclusion of Atomic Dissipation

Figure 10:  Two-parameter phase diagram in the
($\lambda_{-}$,$\lambda_{+}$)-plane for $\kappa=4$, $\omega=0.5$,
$\omega_{0}=0.2$, $\Gamma_{\uparrow}=0$, $\Gamma_{\downarrow}=0.02$. Shown are
curves of Hopf bifurcations $\mathrm{H}_{\mathrm{sub}}$ and
$\mathrm{H}_{\mathrm{spr}}$ (gray and black curves), pitchfork bifurcations
$\mathrm{P}_{\mathrm{sub}}$ and $\mathrm{P}_{\mathrm{spr}}$ (light and dark
blue curves), and saddle-node bifurcations SN (red curves). Coloring shows the
stable regions of the normal phase N (red region), the superradiant phase SR
(blue region), coexistent normal and superradiant phases
$\mathbf{N}+\mathbf{SR}$ (magenta region), and the counter-lasing phase CL
(light-green region). Figure 11:  One-parameter bifurcation diagram
illustrating phase transitions from the normal phase, through the superradiant
phase, into the counter-lasing phase for $\kappa=4$, $\omega=0.5$,
$\omega_{0}=0.2$, $\Gamma_{\uparrow}=0$, $\Gamma_{\downarrow}=0.02$,
$\lambda_{-}=1.5$. Branches with a change in stability are unstable on the
lighter-colored portion of the branch. The bifurcation diagram is divided up
into colored regions where the system is in the normal phase (red), the
superradiant phase (blue), a transitional region (white) between superradiance
and counter-lasing, and where the system is entirely counter-lasing (light
green). Each subsequent panel is a magnification of the rectangular region of
the previous panel.

As in Sec. III.2, the addition of nonzero $\Gamma_{\downarrow}$ destroys the
spin-up normal phase equilibrium and breaks the conservation law Eq. (13).
However, when $\xi\neq 0$ the situation is further complicated by the
possibility of breaking the $\mathbb{Z}_{2}$ symmetry, leading to
superradiance. If $\lambda_{-}$ is small, then varying $\lambda_{+}$ can lead
to a transition from the normal phase to the counter-lasing phase as in Sect.
III.1, that is, via a Hopf bifurcation of the normal phase equilibrium point.
However, the counter-lasing periodic orbits are no longer $U(1)$-symmetric.

####  IV.3.1 Phase Diagram in the Coupling Plane

The addition of nonzero $\Gamma_{\downarrow}$ does not qualitatively change
either case of the superradiant phase transition; it only shifts the phase
transition boundaries. The two cases of the superradiant phase transition, as
well as the Hopf bifurcation leading to the counter-lasing phase for small
$\lambda_{-}$ can be summarized in one figure by computing the curves of
pitchfork, saddle-node, and Hopf bifurcations in the
$(\lambda_{-},\lambda_{+})$-plane. The details of the calculations of these
curves are given in Appendix A. The resulting bifurcation or phase diagram is
illustrated in Fig. 10. We note here that for the case $\Gamma_{\downarrow}=0$
the pitchfork and saddle-node bifurcation curves are symmetric about the
diagonal $\lambda_{-}=\lambda_{+}$, but the addition of
$\Gamma_{\downarrow}\neq 0$ breaks this parameter symmetry also. Here we
choose the parameter values $\kappa=4$, $\omega=0.5$, $\omega_{0}=0.2$,
$\Gamma_{\downarrow}=0.02$, and $\Gamma_{\uparrow}=0$ so that the
$\lambda_{+}$ values of the bifurcations we are concerned with are not too
large. Otherwise, the adiabatic elimination procedure used to derive the LMG
Hamiltonian Eq. (6) from the Dicke Hamiltonian Eq. (II.1) is no longer valid.
The phase diagram shows the super- and subcritical Hopf bifurcation curves,
$\mathrm{H}_{\mathrm{spr}}$ and $\mathrm{H}_{\mathrm{sub}}$, the super- and
subcritical pitchfork bifurcation curves, $\mathrm{P}_{\mathrm{spr}}$ and
$\mathrm{P}_{\mathrm{sub}}$, and a pair of saddle-node bifurcation curves SN.
Also shown in Fig. 11 is an additional “transitional region,” which we will
discuss later in the section. These bifurcation curves divide the phase
diagram into regions where the normal N, superradiant SR, and counter-lasing
CL phases are stable.

Figure 10 outlines several mechanisms by which the counter-lasing phase CL is
arrived at via a sequence of phase transitions from the normal phase N when
$\lambda_{-}$ is kept constant and $\lambda_{+}$ is varied. When $\lambda_{-}$
is small, the transition from the normal phase to the counter-lasing phase
occurs via a supercritical Hopf bifurcation $\mathrm{H}_{\mathrm{spr}}$. This
is the same transition as in Sect. III.2, but now with the $U(1)$ symmetry
explicitly broken. If $\lambda_{-}$ is larger, then upon increasing
$\lambda_{+}$ the $\mathbb{Z}_{2}$-symmetry is spontaneously broken in a
supercritical pitchfork bifurcation $\mathrm{P}_{\mathrm{spr}}$ after which
the system is in the superradiant phase SR. This is case (i) of the
superradiant phase transition. For even larger $\lambda_{-}$, we encounter
case (ii). Here the superradiant phase emerges via a pair of saddle-node
bifurcations SN and the system becomes multistable, with simultaneously stable
normal and superradiant phases, in the region $\textbf{N}+\textbf{SR}$ shown
in the inset panel of Fig. 10. The normal phase is then destabilized at the
subcritical pitchfork bifurcation $\mathrm{P}_{\mathrm{sub}}$, after which the
system is entirely in the superradiant phase SR. Our results indicate that in
every case as $\lambda_{+}$ becomes large the counter-lasing phase must take
hold at some point since the Raman transition described by $\lambda_{+}$
becomes much stronger than that of $\lambda_{-}$.

We remark that the ranges of $\lambda_{\pm}$ in Fig. 10 are large for
illustrative purposes to highlight the structure of the bifurcation curves. In
regions where the couplings are large, the behavior of the LMG model may
depart somewhat from the Dicke model when the cavity field no longer
adiabatically follows the atomic dynamics. In Appendix B, we compare the
results of the semiclassical LMG model with the full semiclassical Dicke model
to show that, in the parameter regime we are concerned with in the following
discussion, the adiabatic elimination procedure is valid.

Figure 12:  Phase space configuration of some relevant orbits at homoclinic
bifurcations and inside the chaotic region. The pair of symmetry-related
homoclinic orbits at the homoclinic bifurcation $\mathrm{Hom}^{1}$ is shown in
panel (a1) in projection onto the $(b_{x},b_{y})$-plane and in panel (a2) in
projection onto the $(b_{x},\gamma)$-plane. Panels (b1) and (b2) show the same
projections of the pair of homoclinic orbits at $\mathrm{Hom}^{2}$, which
feature an additional loop around the opposite superradiant equilibrium point;
only one of the homoclinic orbits is shown in panel (b2). Panels (c1) and (c2)
show in the same way the trajectory with initial condition
$(b_{x},b_{y},\gamma)=(-0.035,-0.023,-0.495)$ which accumulates onto the
chaotic attractor. In panels (a): $\lambda_{+}\approx 1.53208$, (b):
$\lambda_{+}\approx 1.53286$, (c): $\lambda_{+}=1.53292$; other parameters are
set as in Fig. 11.

Once in the superradiant phase, upon increase of $\lambda_{+}$ the system
enters the counter-lasing region CL through a thin transitional region where
many bifurcations occur. This transition represents the mechanisms of how the
superradiant pole-flip transition is perturbed with the addition of atomic
dissipation, where the inverted normal phase is no longer viable. In contrast
to the normal pole-flip transition, which goes from a transition involving
infinitely many periodic orbits to a Hopf bifurcation when atomic dissipation
is included, the superradiant pole-flip transition is made substantially more
complex with the inclusion of atomic dissipation.

Figure 13:  Phase space configurations after the inverse period-doubling
cascade (a) and after the subcritical Hopf bifurcation (b). For (a):
$\lambda_{+}=1.537$, for (b):$\lambda_{+}=1.54$. Other parameters are as in
Fig. 11.

####  IV.3.2 One-Parameter Slice Through the Transitional Region

To investigate the series of bifurcations that lead to the emergence of the
counter-lasing phase inside the transitional region, we take a one-parameter
slice through the phase diagram in Fig. 10 along the vertical line for
$\lambda_{-}=1.5$. The bifurcation diagram illustrating the transition as
$\lambda_{+}$ is varied along this slice is shown in Fig. 11(a); it has been
obtained by numerical continuation of equilibria and periodic orbits with the
program auto-07p [44, 45]. Due to the system’s $\mathbb{Z}_{2}$ equivariance,
this bifurcation diagram is symmetric about the horizontal axis, but for
clarity we are only showing the upper half. The system starts for
$\lambda_{+}$ small in the normal phase $\mathrm{N}_{s}$; here the subscript
denotes stability: $s$ for stable, $u$ for unstable. As $\lambda_{+}$ is
increased the system undergoes the supercritical pitchfork bifurcation
$\mathrm{P}_{\mathrm{spr}}$ from which two stable superradiant equilibria
$\mathrm{SR}_{s}$ emerge, and the normal phase becomes unstable, now labelled
$\mathrm{N}_{u}$. As $\lambda_{+}$ is further increased, the equilibria
$\mathrm{SR}_{s}$ undergo a pair of (subcritical) Hopf bifurcations at
$\mathrm{H}_{\mathrm{sub}}$, where they bifurcate with a pair of unstable
superradiant periodic orbits $\mathrm{SRO}_{u}$. Here, the superradiant
equilibria become unstable ($\mathrm{SR}_{u}$); they disappear, as
$\lambda_{+}$ is increased further, in a subcritical pitchfork bifurcation
$\mathrm{P}_{\mathrm{sub}}$ where they meet the normal phase equilibrium
$\mathrm{N}_{u}$.

Also shown in Fig. 11(a) is a branch of stable counter-lasing periodic orbits
$\mathrm{CL}_{s}$ emerging from a small transitional region where they undergo
interactions with the branch of superradiant periodic orbits
$\mathrm{SRO}_{u}$. A magnification of these branches is shown in Fig. 11(b).
Upon further enlargement in Fig. 11(c), we can glean some information on how
this transitional region is organized. The beginning of the transition is
marked by the creation of two branches of (saddle) periodic orbits,
$\mathrm{SRO}_{u}$ and $\mathrm{CL}_{u}$, in a homoclinic bifurcation labelled
$\mathrm{Hom}^{1}$, where there is a pair of homoclinic orbits to the normal
phase equilibrium $\mathrm{N}_{u}$. As $\lambda_{+}$ is increased, the
(unstable) counter-lasing branch of periodic orbits $\mathrm{CL}_{u}$
undergoes a pitchfork bifurcation of periodic orbits (PPO) upon collision with
another branch of (stable, asymmetric) counter-lasing periodic orbits,
$\mathrm{CL}_{s}^{\mathrm{asym}}$. This bifurcation renders the symmetric
counter-lasing branch stable ($\mathrm{CL}_{s}$).

The pair of homoclinic orbits associated with the homoclinic bifurcation
$\mathrm{Hom}^{1}$ responsible for creating the periodic orbits
$\mathrm{SRO}_{u}$ and $\mathrm{CL}_{u}$ are shown projected onto the
$(b_{x},b_{y})$-plane in Fig. 12(a1) and the $(b_{x},\gamma)$-plane in Fig.
12(a2). The two homoclinic orbits are symmetric counterparts under
$\mathcal{T}$, and are homoclinic to the real saddle $\mathrm{N}_{u}$, which
has one unstable and two stable directions. In the projection in Fig. 12(a1),
the pair of orbits forms a figure eight, while their configuration in three-
dimensional space is that of an _expanding butterfly_ [46, 47]. Let $v_{s}$
and $v_{u}$ be the leading stable and unstable eigenvalues of the Jacobian
matrix evaluated at $\mathrm{N}_{u}$, respectively, with corresponding
eigenvectors $\vec{y}_{s}$ and $\vec{y}_{u}$. The butterfly configuration of
homoclinic orbits occurs when both homoclinic orbits approach $\mathrm{N}_{u}$
from the same side along the direction given by $\vec{y}_{s}$ (the positive
side of the $\gamma$-axis), but leave alongside opposite directions,
$\vec{y}_{u}$ and $-\vec{y}_{u}$. This is illustrated in Fig. 12(a2), where it
can be seen that the pair of homoclinic orbits actually approach
$\mathrm{N}_{u}$ from the same direction along the $\gamma$-axis. This
configuration is referred to as expanding when $|v_{u}|>|v_{s}|$ [47].

Figure 11(d) shows an additional magnification of Fig. 11(c) near the branch
$\mathrm{CL}_{s}^{\mathrm{asym}}$. Here there are infinitely many more
homoclinic bifurcations $\mathrm{Hom}^{2^{k}},\ k=2,3,\cdots$, labelled by the
number of loops $2^{k}$, which each create a pair of saddle periodic orbits.
This leads to the formation of a chaotic attractor upon their accumulation,
shown in Fig. 12(c). As $\lambda_{+}$ is increased further, this chaotic
attractor disappears in a period-doubling cascade (with period-doubling
bifurcations $\mathrm{PD}_{k},\ k=1,2,\cdots$). In the process the saddle
periodic orbits coalesce into the branch of periodic orbits
$\mathrm{CL}_{s}^{\mathrm{asym}}$, which then collides with the branch
$\mathrm{CL}_{u}$ in the bifurcation PPO to create the stable counter-lasing
branch $\mathrm{CL}_{s}$.

The homoclinic orbits at the second homoclinic bifurcation $\mathrm{Hom}^{2}$
are shown in Fig. 12(b1) in the $(b_{x},b_{y})$-plane, and only one of them is
shown in the $(b_{x},\gamma)$-plane in Fig. 12(b2). This bifurcation occurs
along the branch $\mathrm{CL}^{\mathrm{asym}}$ in Fig. 11(d), which is formed
at the pitchfork bifurcation of periodic orbits PPO and becomes unstable at
the first period-doubling bifurcation $\mathrm{PD}_{1}$. Again, there are two
coexistent homoclinic orbits in an expanding butterfly configuration, but now
each features an additional loop in Fig. 12(b1) around the symmetric
counterpart of the superradiant equilibrium point $\mathrm{SR}_{s}$. In the
same fashion, the homoclinic orbits at the subsequent homoclinic bifurcation
$\mathrm{Hom}^{2^{k}},\ k=2,3,\cdots$ feature $2^{k}-1$ additional loops
around $\mathrm{SR}_{s}$. Figure 12(c) shows the chaotic attractor that exists
after the accumulation point of homoclinic bifurcations
$\mathrm{Hom}^{2^{k}},\ k=0,1,\cdots$, also in the same two projections. This
attractor is localized in a small neighborhood of the union of the two initial
homoclinic orbtis; compare with Fig. 12(a). As such, it “surrounds” the two
(stable) superradiant equilibria $\mathrm{SR}_{s}$, and the chaotic
trajectories switch irregularly between circling around (yet well away from)
either of these stable equilibria.

Figure 14: Demonstration of phase coexistence in the semiclassical LMG model.
Panel (a1) shows a trajectory with the same initial condition as Fig. 12(d)
which converges to the chaotic attractor (purple), while a trajectory (blue)
with another initial condition $(b_{x},b_{y},\gamma)=(-0.1,0,-0.475)$ becomes
superradiant. Panel (a2) shows the temporal trace of the two trajectories in
$b_{x}$. Panel (b1) shows a trajectory with initial condition
$(b_{x},b_{y},\gamma)=(-0.1,0.065,-0.475)$ which exhibits counter-lasing
(green), and another trajectory which becomes superradiant (blue), with the
same initial condition as in panel (a1). Panel (b2) also shows the temporal
trace of these two trajectories in $b_{x}$. In panels (a1) and (a2),
$\lambda_{+}=1.53292$, and for panels (b1) and (b2) $\lambda_{+}=1.534$. Other
parameters are set to those in Fig. 11.

As $\lambda_{+}$ becomes larger, the chaotic attractor in Fig. 12(c)
disappears in the (inverse) period-doubling cascade $\mathrm{PD}_{k},\
k=1,2,\cdots$, shown in Fig. 11(d). The pair of periodic orbits remaining
after this process are the (stable) asymmetric counter-lasing orbits
$\mathrm{CL}_{s}^{\mathrm{asym}}$, which subsequently disappear in the
pitchfork bifurcation of periodic orbits PPO; here they meet the unstable
periodic orbit $\mathrm{CL}_{u}$, which then transitions to become stable
($\mathrm{CL}_{s}$). After this final transition in Fig. 11(d), the system now
supports solutions that lase _regularly_ , as opposed to chaotically. Again,
the superradiant periodic orbits $\mathrm{SRO}_{u}$ and $\mathrm{SR}_{s}$
remain unaffected. This configuration of equilibria and periodic orbits is
shown in Fig. 13(a). Further along the bifurcaton diagram in Fig. 11(a), the
periodic orbits $\mathrm{SRO}_{u}$ contract around the equilibria
$\mathrm{SR}_{s}$ until the subcritical Hopf bifurcation
$\mathrm{H}_{\mathrm{sub}}$ is reached, where the periodic orbits disappear
and the equilibria become unstable ($\mathrm{SR}_{u}$). After this
bifurcation, the counter-lasing periodic orbit $\mathrm{CL}_{s}$ is the only
attractor, so the system is entirely in the counter-lasing phase. This
situation is shown in Fig. 13(b).

####  IV.3.3 Coexistent Superradiance and Counter-Lasing

The appearance of a chaotic attractor emerging from lasing states in the
semiclassical model arising from the master equation (3) was observed in Ref.
[34]. However, the situation presented here is slightly different, as the
chaotic attractor here is associated with the appearance of _counter-lasing_
states; that is, the coherent driving terms occur with a greater strength than
the incoherent pumping. As a result, in the chaotic regime, the system
simultaneously supports chaotic counter-lasing and steady-state superradiance.
The superradiant equilibria $\mathrm{SR}_{s}$ in Fig. 11(a) remain stable
until the Hopf bifurcation $\mathrm{H}_{\mathrm{sub}}$, so there is a small
region of overlap leading to phase coexistence of chaotic counter-lasing and
superradiance.

Figure 14(a1) illustrates the dynamics in this region in phase space. Here, we
show two trajectories:

  * •

one which accumulates onto the chaotic attractor, which lies beyond, yet
surrounds the basins of attraction of the two superradiant equilibria
$\mathrm{SR}_{s}$,

  * •

and another, which approaches a superradiant state.

The dramatic difference in the behavior generated by the coexistence of these
different types of attractors is highlighted in Fig. 14(a2), which shows the
temporal trace of $b_{x}$ of the two trajectories shown in panel (a1). The
trajectory on the chaotic attractor exhibits an irregular and nonrepeating
dynamics with intermittent transitions between larger oscillations around
either of the superradiant states, i.e., the trajectory chaotically switches
between positive and negative $b_{x}$. This trajectory corresponds to the
atomic ensemble lasing chaotically. Meanwhile, the trajectory in the basin of
the superradiant equilibrium simply decays into the steady state. We have
observed qualitatively that the basin of attraction of the chaotic attractor
is quite small and includes the normal phase equilibrium $\mathrm{N}_{u}$. In
particular, this means that trajectories with initial conditions sufficiently
close to $\mathrm{N}_{u}$ accumulate onto the chaotic attractor and begin to
lase chaotically. On the other hand, initial conditions with higher levels of
excitation, that is, further away from $\mathrm{N}_{u}$, cause the atoms to
begin to emit in a superradiant fashion, and eventually, regularly.

In the parameter regime after the (inverse) period-doubling cascade
$\mathrm{PD}_{k},\ k=1,2,\cdots$ and the pitchfork bifurcation of periodic
orbits PPO, the chaotic attractor disappears and is replaced by a (symmetric
under $\mathcal{T}$) stable counter-lasing periodic orbit $\mathrm{CL}_{s}$.
As mentioned previously, the superradiant equilibria are unaffected by this
transition. Thus, before the Hopf bifurcation $\mathrm{H}_{\mathrm{sub}}$, the
system features phase coexistence of superradiance and regular counter-lasing.
Figure 14(b1) shows two trajectories corresponding to the simultaneously
existing superradiant and counter-lasing phases. Temporal traces of the
trajectories’ evolution are shown in Fig. 14(b2). Similarly to the case
discussed above, initial conditions near $\mathrm{N}_{u}$ lead to periodic
counter-lasing, while for most initial conditions atomic ensembles are
emitting in a superradiant manner.

###  IV.4 Global Bifurcations via Kneading Sequences

In order to investigate the dynamics leading to the establishment of chaos in
more detail, we now turn our attention to characterizing relevant trajectories
of Eqs. (10)–(11) with the use of parameter sweeping techniques based on the
calculation of a _kneading sequence_. Since their inception in Ref. [48],
these techniques have been used successfully to characterize the complicated
arrangement of homoclinic (and other global) bifurcations in parameter space
close to special organizing bifurcation points [48, 49, 50].

The kneading sequence encodes how a particular trajectory explores phase
space. In particular, this facilitates an understanding of how complex
switching patterns change as parameters are varied, thus, elucidating complex
bifurcation orderings. To characterize kneading sequences, we will also make
use of a mapping from the space of kneading sequences to the interval $[0,1]$,
known as the _kneading invariant_ $K$. Thus, by sweeping $\lambda_{+}$ while
holding other parameters constant, we compute (an approximation of) the
kneading invariant $K$ for each $\lambda_{+}$ to show how the transitional
region is organized. In doing so, we discover a complex recursive fractal
structure, arranged around various global bifurcations, that gives insight
into the emergence of complicated dynamics in this model.

To compute the kneading sequence, we consider the _unstable manifold_
$W^{u}(\mathrm{N}_{u})$ of the saddle equilibrium $\mathrm{N}_{u}$, which is
the set of all points in phase space that converge towards $\mathrm{N}_{u}$ in
backwards time, that is,

$W^{u}(\mathrm{N}_{u})=\bigg{\\{}x\in\mathbb{R}^{3}\bigg{|}\lim_{t\rightarrow-\infty}\varphi_{t}x=\mathrm{N}_{u}\bigg{\\}},$
(39)

where $\varphi_{t}$ is the flow or time-evolution operator of Eqs. (10)–(11).
In the parameter regime we are considering here, the unstable manifold
$W^{u}(\mathrm{N}_{u})$ is one-dimensional, and changes in its kneading
sequence allows one to study global bifurcations [48] as $\lambda_{+}$ is
varied. The manifold $W^{u}(\mathrm{N}_{u})$ has two branches,
$W^{u}_{+}(\mathrm{N}_{u})$ and $W^{u}_{-}(\mathrm{N}_{u})$, which leave
$\mathrm{N}_{u}$ tangent to the unstable eigenvector $\vec{y}_{u}$ in the
positive and negative direction, respectively. As the normal phase equilibrium
$\mathrm{N}_{u}$ lies within the parity symmetry subspace, the two branches of
the unstable manifold are each other’s symmetric counterparts; that is,
$\mathcal{T}[W^{u}_{+}(\mathrm{N}_{u})]=W^{u}_{-}(\mathrm{N}_{u})$. As such,
we only need to consider a single branch of the manifold to fully categorize
the dynamics.

We consider here the kneading sequence defined by $W^{u}_{-}(\mathrm{N}_{u})$.
More specifically, we label its successive minima in $b_{x}$ with $b_{x}<0$
with the symbols $\mathtt{0}$, and maxima in $b_{x}$ with $b_{x}>0$ with the
symbol $\mathtt{1}$. To avoid misidentifying minima and maxima due to small
oscillations near $\mathrm{N}_{u}$, we only consider those that fall outside
of an _exclusion zone_ $-0.2<b_{x}<0.2$. This procedure defines the symbolic
sequence of $\mathtt{0}$s and $\mathtt{1}$s generated by
$W^{u}_{-}(\mathrm{N}_{u})$, known as the _kneading sequence_. It is unique
for every point in parameter space and encodes the order in which
$W^{u}_{-}(\mathrm{N}_{u})$ visits the two regions of phase space where
$b_{x}$ is either negative or positive. Each kneading sequence is denoted
$S=(a_{1}a_{2}a_{3}\cdots)$, where $a_{k}\in\\{\mathtt{0},\mathtt{1}\\}$ for
$k\in\mathbb{N}$, and is a member of the set of all binary sequences, i.e.,
$S\in\\{\mathtt{0},\mathtt{1}\\}^{\mathbb{N}}$. Kneading sequences that end in
a repeating infinite string of symbols will be denoted with an overbar, e.g.,
$\bar{\mathtt{0}}\equiv\mathtt{0}\mathtt{0}\mathtt{0}\cdots$ and
$\overline{\mathtt{0}\mathtt{1}}\equiv\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\cdots$.
Due to $\mathbb{Z}_{2}$ equivariance, the kneading sequences generated by
either branch $W^{u}_{-}(\mathrm{N}_{u})$ or $W^{u}_{+}(\mathrm{N}_{u})$ are
mapped to one another in the following way: for the kneading sequence
$S=(a_{1}a_{2}a_{3}\cdots)$ of $W^{u}_{-}(\mathrm{N}_{u})$ the associated
kneading sequence of $W^{u}_{+}(\mathrm{N}_{u})$ takes the form

$\neg S=(\neg a_{1}\neg a_{2}\neg a_{3}\cdots),$ (40)

where $\neg$ is Boolean negation, that is,

$\neg a_{k}=\begin{cases}\mathtt{1},&\mathrm{if}\ a_{k}=\mathtt{0},\\\
\mathtt{0},&\mathrm{if}\ a_{k}=\mathtt{1}.\end{cases}$ (41)

As an illustration, Fig. 15 shows the example where the unstable manifold
generates the kneading sequence
$S=(\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$,
which we also write simply as
$\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}$ for
notational convenience.

Figure 15:  Kneading sequence
$S=(\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$
generated by $W^{u}_{-}(\mathrm{N}_{u})$. Each negative minimum or positive
maximum of $b_{x}$ outside of the exclusion zone $(-0.2,0.2)$ creates a symbol
0 or 1, respectively. Here, $\lambda_{+}=1.532895$ and other parameters are
set to those in Fig. 11.

The key point here is that the kneading sequence of
$W^{u}_{-}(\mathrm{N}_{u})$ is a topological invariant that allows one to
identify global bifurcations [48]. This means that the kneading sequence is
constant in open regions of parameter space and, if it is different at two
parameter values, then a homoclinic or heteroclinic bifurcation must have
occurred. Thus, as $\lambda_{+}$ is varied, the kneading sequence of
$W^{u}_{-}(\mathrm{N}_{u})$ can be used to identify the relevant bifurcations
on the route to chaotic dynamics [48, 49, 50].

Typically (in open regions of parameter space), $W^{u}_{-}(\mathrm{N}_{u})$
has infinitely many extrema and so the kneading sequence indeed contains
infinitely many symbols. However, homoclinic orbits, such as those shown in
Fig. 12, make excursions exploring phase space only for some time before
eventually closing back up at $\mathrm{N}_{u}$. We therefore describe a
homoclinic orbit with a kneading sequence of only finitely many symbols, which
are derived in the same way and encode how $W^{u}_{-}(\mathrm{N}_{u})$ visits
the regions with negative and positive $b_{x}$, respectively, before returning
to $\mathrm{N}_{u}$; in particular, the bit-length of such a finite kneading
sequence is given by the number of loops of the homoclinic orbit. Hence, the
finite kneading sequence of the homoclinic orbit formed by
$W^{u}_{-}(\mathrm{N}_{u})$ in Fig. 12(a) is $\mathtt{0}$ and that in Fig.
12(b) is $\mathtt{0}\mathtt{1}$; their symmetric counterparts formed by
$W^{u}_{+}(\mathrm{N}_{u})$ have kneading sequences $\mathtt{1}$ and
$\mathtt{1}\mathtt{0}$, respectively.

To facilitate an analysis of the kneading sequences, we consider the _kneading
invariant_

$K=\sum_{k=1}^{\infty}\frac{a_{k}}{2^{k}}\in[0,1],$ (42)

which assigns a real number to each kneading sequence. In this way, performing
a parameter sweep across $\lambda_{+}$ and computing the manifold
$W_{-}^{u}(\mathrm{N}_{u})$, we can associate a kneading invariant $K$ for
each value of $\lambda_{+}$ and identify global bifurcations for values where
$K$ changes. Note that $K$ is not a one-to-one mapping, since the two kneading
sequences

$a_{1}a_{2}\cdots a_{k}\mathtt{1}\bar{\mathtt{0}}\ \text{and}\
a_{1}a_{2}\cdots a_{k}\mathtt{0}\bar{\mathtt{1}}$

have the same kneading invariant $K$. Similarly, the finite kneading invariant
of a homoclinic orbit is equal to the kneading invariant of the same sequence
followed by an infinite string of $\mathtt{0}$s.

In practice we can only consider kneading invariants arising from sequences of
finitely many computed symbols, so we introduce the kneading invariant

$K_{n}=\sum_{k=1}^{n}\frac{a_{k}}{2^{k}}$ (43)

of sequences of bit-length $n$, which we use as an approximation of the full
kneading invariant $K$. Note that $K_{n}$ divides the interval $[0,1]$ into
$2^{n}$ subintervals, whose values of $K_{n}$ differ by at least $2^{-n}$; for
our computations we use $n=12$ throughout.

Figure 16:  Parameter sweep through $\lambda_{+}$ of the 12-symbol kneading
invariant $K_{12}$. Panel (a) shows the parameter values for the bifurcations
of Fig. 11(d), panel (b) is a magnification of the relevant
$\lambda_{+}$-range, and panel (c) shows the associated kneading invariant
$K_{12}$, where the two frames indicate parts that are enlarged in Fig. 17.
Other parameters are as in Fig. 11. Figure 17:  Magnifications of the kneading
invariant sweep in the two frames of Fig. 16(c). Panel (a) shows the structure
near the homoclinic $\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$, and panel (b) shows
the staircase formed between the EtoP connections
$\mathrm{EtoP}_{\mathtt{0}\bar{\mathtt{1}}}$ and
$\mathrm{EtoP}_{\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}}$.

To see how the transitional region in Fig. 11 is organized, we now execute a
parameter sweep that determines the kneading invariant $K_{12}$ as a function
of $\lambda_{+}$. To indicate the $\lambda_{+}$-interval of the sweep, we plot
in Fig. 16(a) the values of $\lambda_{+}$ for which there are homoclinic,
period-doubling, and pitchfork bifurcations of periodic orbits already found
by numerical continuation; see Fig. 11. An enlargement of the relevant range
between the bifurcations $\mathrm{Hom}^{2}$ and $\mathrm{PPO}$ is shown in
Fig. 16(b). Over a fine grid of this $\lambda_{+}$-interval, we compute, for
each value of $\lambda_{+}$ in the sweep, the one-dimensional unstable
manifold $W^{u}_{-}(\mathrm{N}_{u})$ accurately up to a length that allows us
to determine the associated 12-symbol initial part of the kneading sequence;
this is achieved with the software package tides [51], which uses a Taylor
series method for high-precision numerical integration of ordinary
differential equations.

Figure 16(c) shows the resulting plot of the kneading invariant $K_{12}$, and
we first discuss its overall features. At the start of the sweep there is an
interval where the kneading sequence is $\mathtt{0}\bar{\mathtt{1}}$ with a
kneading invariant $\frac{1}{2}$, and at the end of the sweep there is an
interval where it is $\overline{\mathtt{0}\mathtt{1}}$ with a kneading
invariant $\frac{1}{3}$. In between, the kneading invariant is constant on
smaller and smaller intervals and undergoes a remarkably complex series of
changes that correspond to the sequence of global bifurcations that give rise
to homoclinic chaos [48]. On the left of Fig. 16(c) (in the larger frame)
there is a descending “staircase” structure; it will be discussed later in the
section. Following this, as $\lambda_{+}$ increases further, we find a
“plateau” where the kneading invariant is constant over an interval, followed
by a sharp, spiked increase and decrease of $K_{12}$ (in the smaller frame).
It is centered around a homoclinic bifurcation $\mathrm{Hom}^{2}$, which we
now refer to by its finite kneading sequence as
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$; similarly, we refer to any homoclinic
bifurcation as $\mathrm{Hom}_{S}$ from now on, where $S$ is its finite
kneading sequence. Notice that there is another spiked increase and decrease
of $K_{12}$ near
$\mathrm{Hom}^{4}=\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}}$.
The nature of this large local change of $K_{12}$ will be discussed below. As
$\lambda_{+}$ is increased further, near the central range in Fig. 16(c) where
we expect a chaotic attractor to appear, $K_{12}$ flattens out; this is due to
changes in the kneading sequences occurring only at later and later symbols,
with correspondingly smaller changes to the kneading invariant. As
$\lambda_{+}$ is further increased, the kneading invariant makes two large
jumps. First down to $\frac{7}{24}$ and an interval where the kneading
sequence is $\mathtt{0}\mathtt{1}\mathtt{0}\overline{\mathtt{0}\mathtt{1}}$,
after which there is a second jump at the pitchfork bifurcation of periodic
orbits PPO, up to $\frac{1}{3}$ and the final interval with kneading sequence
$\overline{\mathtt{0}\mathtt{1}}$, which corresponds to the stabilization of
the symmetric counter-lasing phase.

Figure 17 shows magnifications of the kneading invariant plot of Fig. 16(c) in
the two highlighted frames. We first turn our attention to the plateau and the
spike near the homoclinic bifurcation $\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$,
which is shown over an enlarged $\lambda_{+}$-range in Fig. 17(a). This closer
inspection reveals that the sharp spike in the kneading invariant is a smaller
structure with additional plateaux with spikes at further homoclinic
bifurcations, thus, revealing an interesting type of self-similarity. Notice
that there are two kneading sequences associated with each plateau, one to the
left of the homoclinic bifurcation, and one to the right. Specifically, the
start and end of the main plateau in Fig. 17(a) is generated by the kneading
sequences $S_{\ell}=(\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$ and
$S_{r}=(\mathtt{0}\mathtt{1}\mathtt{0}\bar{\mathtt{1}})$, respectively, which
both have the same kneading invariant $K=\frac{3}{8}$ (but differ by $2^{-12}$
in $K_{12}$). The ‘spike’ is identified here as consisting of two additional
staircase structures, similar to the one seen in Fig. 16(c) — an increasing
one to the left of $\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$ and a decreasing one
to the right of $\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$, which meet at the
central organizing homoclinic bifurcation
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$ with a discontinuity in the kneading
invariant $K$. Each staircase structure in Fig. 17(a) has its own central,
largest plateau with a spike at the homoclinic bifurcations
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}\mathtt{0}}$ and
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}}$,
respectively, and within each of the plateau of the staircase we find further
staircases with infinitely more plateaux and staircases, and so on.

The staircase structure we find highly resembles the Cantor function [52] —
also known as the Cantor staircase or “Devil’s staircase” — which is a
continuous and monotone (increasing or decreasing) function that is locally
constant. However, the graph of $K$ we find here has the additional spikes in
the middle of each plateau, so is neither monotone nor continuous and features
additional self-similarity since each spike consists of two such staircase
structures. We therefore refer to the graph of $K$ shown in Fig. 17(a) as a
_spiked (Cantor) staircase_ , and we now proceed to consider in more detail
how it is organized recursively by sequences of homoclinic orbits.

Firstly, all kneading sequences of the spike in Fig. 17(a) start with the
common first two symbols $\mathtt{0}\mathtt{1}$ of $S_{\ell}$ and $S_{r}$, the
kneading sequences at the plateau we start with. Moreover, for the entire left
spiked staircase for
$\lambda_{+}<\lambda_{+}^{\mathrm{Hom_{\mathtt{0}\mathtt{1}}}}$ this is
followed by the symbol $\mathtt{1}$, which means that the kneading invariant
$K$ is larger than the value of the plateau itself. In contrast, for the right
spiked staircase with
$\lambda_{+}>\lambda_{+}^{\mathrm{Hom}_{\mathtt{0}\mathtt{1}}}$ any kneading
sequence always starts with $\mathtt{0}\mathtt{1}\mathtt{0}$, meaning the
kneading invariant is always less than the value of the plateau. This
combinatorial structure is repeated at each sub-plateau, as we checked for the
respective largest one, where we find spikes (that is, additionally left and
right spiked staircases) near the homoclinic orbits with the finite kneading
sequences $\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}\mathtt{0}$ and
$\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}$. We remark that we have
identified these and other additional homoclinic orbits by computing them with
a Lin’s method approach [53]; see Appendix C for details.

As a general observation, we find the following for each homoclinic
bifurcation $\mathrm{Hom}_{a_{1}a_{2}\cdots a_{k}}$. For each plateau of its
left spiked staircase with kneading sequence

$S_{\ell}=(a_{1}a_{2}\cdots a_{k}\,\mathtt{1}\,a_{k+2}\cdots)$

there is a corresponding plateau on its right spiked staircase with the
kneading sequence

$S_{r}=(a_{1}a_{2}\cdots a_{k}\ \neg[\mathtt{1}\,a_{k+2}\cdots]).$

This relationship of negating/inverting the symbols after the common prefix
determines a map from the left to the right spiked staircase near each
homoclinic bifurcation and, thus, explains the observed symmetry in Fig.
17(a). As a concrete example, note that the kneading sequence of the main
plateau to the left of the homoclinic bifurcation
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$ is
$S_{\ell}=(\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$, which hence maps
to its right counterpart with kneading sequance
$S_{r}=(\mathtt{0}\mathtt{1}\neg[\mathtt{1}\bar{\mathtt{0}}])=(\mathtt{0}\mathtt{1}\mathtt{0}\bar{\mathtt{1}})$.
Note further, that this observation also explains the self-similar structure
of the graph of $K$: apart from the common prefix and the fact that
$a_{k+1}=\mathtt{1}$ for any plateau of the left spiked staircase, there is no
restriction on the symbols $a_{i}\in\\{\mathtt{0},\mathtt{1}\\}$ for $i\geq
k+2$. Hence, at every level, there is a full binary tree of further spiked
plateaux of $K$, each with its own central homoclinic bifurcation.

Figure 17(b) is an enlargement of the first part of the staircase structure in
Fig. 16(c) between the initial plateau with kneading sequence
$\mathtt{0}\bar{\mathtt{1}}$ with $K=\frac{1}{2}$ and the next largest plateau
with kneading sequence $\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}$ with
$K=\frac{3}{2}$. At first glance, it looks like a regular descending Devil’s
staircase, but it does also feature spikes at every sub-plateau, the most
prominent of which is the spike in Fig. 17(a) we just discussed. Hence, it is
also a spiked staircase, and we find that it features plateaux and associated
homoclinic bifurcations with kneading sequences composed of all sequences
arising from the binary tree starting with the symbols
$\mathtt{0}\mathtt{1}\mathtt{1}$. Specifically, the largest plateau in panel
(b) is centered around $\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}}$, with
the next two largest inheriting the prefix and centered around
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}}$ and
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{1}}$, which are
homoclinic orbits that we also identified numerically (see again Appendix C).
While the spikes in panel (b) near these homoclinic bifurcations are smaller
(extend less in $K$ from the value of the respective plateau) compared to
panel (a), each indeed consists of a left and a right spiked sub-staircase.

It is clear from the above discussion that the number of homoclinic
bifurcations $\mathrm{Hom}_{S}$, as given by their finite symbol sequences
$S$, is countably infinite. In particular, this means that there must be
accumulation points of homoclinic bifurcations in the finite
$\lambda_{+}$-range where they occur. As we have checked and show now, the two
boundaries of each plateau, that is, the left and right end points of a spike,
are limits of specific sequences of homoclinic bifurcations, namely those with
more and more symbols $\mathtt{0}$ and $\mathtt{1}$ after the common prefix,
respectively. For example, in Fig. 17(a) the sequence
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}\mathtt{0}}$, $\cdots$
converges to the end point of the part of the plateau with kneading sequence
$S_{\ell}=(\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$, and the sequence
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}\mathtt{1}\mathtt{1}}$, $\cdots$
converges to the end point of the plateau with kneading sequence
$S_{r}=(\mathtt{0}\mathtt{1}\mathtt{0}\bar{\mathtt{1}})$. Similarly, in Fig.
17(b) the sequence $\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{1}\mathtt{1}}$, $\cdots$
limits on the plateau with kneading sequence $\mathtt{0}\bar{\mathtt{1}}$ at
$K=\frac{1}{2}$, and the sequence
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{0}\mathtt{0}}$, $\cdots$
on that with kneading sequence
$\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}$ at $K=\frac{3}{8}$.

It turns out that the limits of the respective sequences of homoclinic orbits
exist as objects in phase space: they are heteroclinic connections from the
saddle equilibrium $\mathrm{N}_{u}$ to either of the symmetric pair of saddle
periodic orbits labelled $\mathrm{SRO}_{u}$ in Fig. 13(a). We refer to such
connecting orbits as _EtoP connections_ , and they occur at isolated
$\lambda_{+}$-values where the one-dimensional unstable manifold
$W^{u}(\mathrm{N}_{u})$ lies in the two-dimensional stable manifold
$W^{s}(\mathrm{SRO}_{u})$. We distinguish each limiting EtoP connection
$\mathrm{EtoP}_{S}$ by its associated kneading sequence $S$, which either ends
in $\bar{\mathtt{0}}$ or $\bar{\mathtt{1}}$ when $W^{u}_{-}(\mathrm{N}_{u})$
accumulates on the periodic orbit with $b_{x}<0$ or with $b_{x}>0$,
respectively. We remark that EtoP connections can also be identified reliably
with a Lin’s method approach [53]; see Appendix C. In this way, we found the
limiting EtoP connections bounding the plateau with kneading sequences
$S_{\ell}=(\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}})$ and
$S_{r}=(\mathtt{0}\mathtt{1}\mathtt{0}\bar{\mathtt{1}})$ and with kneading
sequences $\mathtt{0}\bar{\mathtt{1}}$ and
$\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}$, the $\lambda_{+}$-values of
which are indicated in panels (a) and (b) of Fig. 17.

Figure 18:  Convergence of homoclinic orbits to an EtoP connection. Panels
(a)–(c) show the temporal traces of the homoclinic orbits $\mathrm{Hom}_{01}$,
$\mathrm{Hom}_{011}$ and $\mathrm{Hom}_{0111}$ as they accumulate on
$\mathrm{EtoP}_{0\bar{1}}$, shown as a temporal trace in panel (d) and in
$(b_{x},b_{y},\gamma)$-space in panel (e).

The convergence of homoclinic orbits is illustrated in Fig. 18 for the example
of $\mathrm{EtoP}_{\mathtt{0}\bar{\mathtt{1}}}$. As the temporal traces of
$W^{u}_{-}(\mathrm{N}_{u})$ in panels (a)–(c) show, the homoclinic orbits
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}}$, and
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{1}\mathtt{1}}$ each have one
further excursion around the periodic orbit $\mathrm{SRO}_{u}$ with $b_{x}>0$
before $W^{u}_{-}(\mathrm{N}_{u})$ returns to the saddle equilibrium
$\mathrm{N}_{u}$. In the limit, $W^{u}_{-}(\mathrm{N}_{u})$ accumulates on
$\mathrm{SRO}_{u}$, generating the EtoP connection
$\mathrm{EtoP}_{\mathtt{0}\bar{\mathtt{1}}}$, which is shown as a temporal
profile in Fig. 18(d) and as a connecting orbit in
$(b_{x},b_{y},\gamma)$-space in panel (e).

The limiting EtoP connections exist as part of the overall self-similar
structure of the spiked staircase, where each plateau is bounded by the
respective two EtoP connections. The tree of kneading sequences of a plateau
and its spike is determined by its central homoclinic orbit $\mathrm{Hom}_{S}$
with given finite sequence $S$, which provides the prefix for all further
homoclinic bifurcations and EtoP connections, as well as the associated
plateaux. Note that plateaux with prefixes of large bit-lengths (before
periodic repetition) correspond to trajectories in phase space that feature
long chaotic transients before they settle down to either of the pair of
superradiant equilibria $\mathrm{SR}_{s}$; see Fig. 13(a). These chaotic
transients may become arbitrarily long as part of the self-similar structure
of the spiked Cantor staircase in Figs. 16 and 17. We therefore suggest that
in this range of $\lambda_{+}$-values prior to when the chaotic attractor
emerges, Eqs. (10)–(11) displays _preturbulence_ — in recognition of the
analogy with the preturbulent regime of the Lorenz system prior to its
transition to a chaotic attractor due to a single EtoP connection [54].

##  V Conclusions

We analyzed the nonlinear semiclassical physics arising from phase transitions
leading to lasing and counter-lasing phases in a generalized Lipkin-Meskov-
Glick model derived from a generalized Dicke model after adiabatic elimination
in the dissipative limit. We first performed a study of the model in what we
term the strong dissipative limit, where the dissipative terms of the quantum
master equation greatly outweigh the unitary dynamics. In doing so, we have
described the effect of adding spontaneous emission to the dynamics of a
collective spin-flip transition. We call this the _normal pole-flip
transition_ , and it manifests itself as a Hopf bifurcation to a counter-
lasing periodic orbit. This can lead to resurgent superradiant pulses driven
by spontaneous emission. We then presented a detailed analysis of the dynamics
associated with the appearance of these counter-lasing periodic orbits, as
well as how they can be reconciled with the inverted (spin-up) normal phase
that appears in the limit of disappearing spontaneous emission.

We then focussed on the transition to the counter-lasing phase once the
system’s parity symmetry is broken in a superradiant phase transition, as well
as its relation to another collective spin-flip transition involving
superradiant states, which we call the _superradiant pole-flip transition_. We
discovered a surprisingly complex sequence of bifurcations that organize this
transition, which involves the establishment of a chaotic attractor.
Specifically, we identified sequences of homoclinic bifurcations to the now
unstable symmetric equilibrium, as well as sequences of connecting EtoP orbits
from it to saddle periodic orbits, and showed that these are organized in an
intricate structure we call a spiked Cantor staircase when one considers their
associated topological quantity known as the kneading invariant. As we showed,
this phenomenon entails long chaotic transients, which we characterized as a
preturbulent regime that exists before the chaotic attractor emerges.

Overall, this work demonstrates the appearance of highly complicated nonlinear
dynamics from a dissipative many-body quantum optical model in the
semiclassical limit. In this description, the dissipative phase transitions
were shown to result from the interplay between the coherent processes that
underpin superradiance and the incoherent processes that seed counter-lasing,
and can be very complex in the transitional region where the two are in
competition.

A possible future direction of research that this work suggests is an analysis
of the semiclassical dynamics of the full unbalanced Dicke model with atomic
dissipation, outside of the dissipative limit where the use of the LMG model
description is justified. Here, rather than three-dimensional, the
semiclassical dynamics are actully five-dimensional by virtue of the dynamics
of the light field no longer adiabatically following the atomic states. It
would be interesting to investigate how, as the validity of the adiabatic
elmination procedure is departed, the dynamics described in this paper are
modified by the states of the light field lagging further and further behind
the atomic states. Moreover, once entirely outside of the dissipative regime,
there is the possibility of the appearance of additional complex dynamics and
phase transition mechanisms between superradiance and counter-lasing in
addition to those analyzed here.

##  Appendix A Calculation of Local Bifurcation Conditions in the LMG Model

The Jacobian matrix corresponding to Eqs. (10)–(11) evaluated at the normal
phase equilibium $b_{x}=b_{y}=0$, $\gamma=\gamma_{\mathrm{eq}}$ is

$J=\begin{bmatrix}2\eta\gamma_{\mathrm{eq}}(\lambda_{-}^{2}-\lambda_{+}^{2})-(\Gamma_{\downarrow}+\Gamma_{\uparrow})&\omega_{0}+2\xi\gamma_{\mathrm{eq}}(\lambda_{-}-\lambda_{+})^{2}&0\\\
-\omega_{0}-2\xi\gamma_{\mathrm{eq}}(\lambda_{-}+\lambda_{+})^{2}&2\eta\gamma_{\mathrm{eq}}(\lambda_{-}^{2}-\lambda_{+}^{2})-(\Gamma_{\downarrow}+\Gamma_{\uparrow})&0\\\
0&0&-2(\Gamma_{\downarrow}+\Gamma_{\uparrow})\end{bmatrix}.$ (44)

###  A.1 Pitchfork Bifurcation Condition

With the nonzero $\xi$-terms included, the eigenvalues of this matrix are
complicated and are not of much help in determining the locations of the
pitchfork or Hopf bifurcations in the $(\lambda_{-},\lambda_{+})$-plane.
Instead, a necessary condition for a pitchfork bifurcation curve describing
the breaking of the LMG model’s $\mathbb{Z}_{2}$ symmetry is that

$\det(J)=-2(\Gamma_{\uparrow}+\Gamma_{\downarrow})\bigg{[}(\Gamma_{\downarrow}+\Gamma_{\uparrow})^{2}+4\gamma_{\mathrm{eq}}\eta(\Gamma_{\downarrow}+\Gamma_{\uparrow})(\lambda_{+}^{2}-\lambda_{-}^{2})+4\gamma_{\mathrm{eq}}^{2}(\eta^{2}+\xi^{2})(\lambda_{+}^{2}-\lambda_{-}^{2})^{2}+4\gamma_{\mathrm{eq}}\xi\omega_{0}(\lambda_{+}^{2}+\lambda_{-}^{2})+\omega_{0}^{2}\bigg{]}=0.$
(45)

This equation is quadratic in $\lambda_{+}^{2}$, so an analytic solution is
readily found and subsequently confirmed to be a pitchfork bifurcation.
Asserting $\Gamma_{\downarrow}+\Gamma_{\uparrow}>0$, this gives the locations
of pitchfork bifurcation shown in Fig. 10.

###  A.2 Saddle-node Bifurcation Condition

The Jacobian matrix $J$ also yields a necessary condition for saddle-node
bifurcations when it is evaluated at the superradiant equilibria (obtained
with a computer algebra package), giving from $\det(J)=0$ the equation

$2\eta\xi\omega_{0}(\Gamma_{\downarrow}+\Gamma_{\uparrow})(\lambda_{-}^{4}-\lambda_{+}^{4})+\bigg{[}\xi^{2}(\Gamma_{\downarrow}+\Gamma_{\uparrow})^{2}+\eta^{2}\omega_{0}^{2}\bigg{]}(\lambda_{-}^{2}-\lambda_{+}^{2})^{2}-4\xi^{2}\omega_{0}^{2}\lambda_{-}^{2}\lambda_{+}^{2}=0,$
(46)

which turns out to be linear in $\lambda_{-}$, yielding

$\lambda_{+}=\mp\left\\{\frac{\sqrt{\xi^{2}(\Gamma_{\uparrow}+\Gamma_{\downarrow})^{2}+\omega_{0}^{2}(\eta^{2}+2\xi^{2})\pm
2\xi\omega_{0}\sqrt{(\eta^{2}+\xi^{2})\big{[}(\Gamma_{\uparrow}+\Gamma_{\downarrow})^{2}+\omega_{0}^{2}\big{]}}}}{\xi(\Gamma_{\uparrow}+\Gamma_{\downarrow})-\eta\omega_{0}}\right\\}\
\lambda_{-}.$ (47)

We validated with numerical continuation that Eq. (47) indeed gives the curves
of saddle-node bifurcations of the superradiant equilibria.

###  A.3 Hopf Bifurcation Condition

A necessary condition for the Hopf bifurcation is formulated in terms of the
_bialternate product_ $\odot$ of the Jacobian matrix $J$ with twice the
identity $I$ [42] (evaluated at the normal phase equilibrium), which is

$J\odot
2I=\begin{bmatrix}4\eta\gamma_{\mathrm{eq}}(\lambda_{-}^{2}-\lambda_{+}^{2})-2(\Gamma_{\uparrow}+\Gamma_{\downarrow})&0&0\\\
0&2\eta\gamma_{\mathrm{eq}}(\lambda_{-}^{2}-\lambda_{+}^{2})-3(\Gamma_{\uparrow}+\Gamma_{\downarrow})&\omega_{0}+2\xi\gamma_{\mathrm{eq}}(\lambda_{-}-\lambda_{+})^{2}\\\
0&-\omega_{0}-2\xi\gamma_{\mathrm{eq}}(\lambda_{-}+\lambda_{+})^{2}&2\eta\gamma_{\mathrm{eq}}(\lambda_{-}^{2}-\lambda_{+}^{2})-3(\Gamma_{\downarrow}+\Gamma_{\uparrow})\end{bmatrix}.$
(48)

Its determinant is

$\det(J\odot
2I)=-2\bigg{[}2\eta\gamma_{\mathrm{eq}}(\lambda_{+}^{2}-\lambda_{-}^{2})+(\Gamma_{\uparrow}+\Gamma_{\downarrow})\bigg{]}F(\lambda_{-},\lambda_{+}),$
(49)

where $F(\lambda_{-},\lambda_{+})$ is a quartic in $\lambda_{-}$ and
$\lambda_{+}$. The eigenvalues of the matrix $J\odot 2I$ are the pairwise sums
of the eigenvalues of $J$ [42]. Thus, there are two ways that $\det(J\odot
2I)=0$, which is the condition of interest: the matrix $J$ either has a
complex conjugate pair of eigenvalues with zero real part (that is, a Hopf
bifurcation), or it has a pair of real eigenvalues of $J$ that sum to zero. We
identify the latter condition to correspond to $F(\lambda_{-},\lambda_{+})=0$,
so the defining condition for the Hopf bifurcation is

$2\eta\gamma_{\mathrm{eq}}(\lambda_{+}^{2}-\lambda_{-}^{2})+(\Gamma_{\uparrow}+\Gamma_{\downarrow})=0.$
(50)

This is exactly the condition (30) in the strong dissipative limit $\xi=0$,
which was again confirmed by numerical continuation. We remark that Eq. (50)
is a necessary condition, as the zeros of Eq. (49) can also be used to compute
Hopf-degeneracies such as the pole-flip transition [26].

##  Appendix B Validity of Adiabatic Elimination

Here we will briefly discuss the validity of the adiabatic elimination process
that leads to the semiclassical Lipkin-Meshkov-Glick (LMG) model [Eqs.
(10)–(11)] in the parameter regime of Fig. 10. For comparison, the
semiclassical equations of motion arising from the Dicke model master equation
(3) are

$\displaystyle\frac{d\alpha}{dt}={}$
$\displaystyle-\kappa\alpha-i\omega\alpha-i\lambda_{-}\beta-i\lambda_{+}\beta^{*},$
(51) $\displaystyle\frac{d\beta}{dt}={}$
$\displaystyle-i\omega_{0}\beta+2i\lambda_{-}\alpha\gamma+2i\lambda_{+}\alpha^{*}\gamma-(\Gamma_{\uparrow}+\Gamma_{\downarrow})\beta,$
(52) $\displaystyle\frac{d\gamma}{dt}={}$ $\displaystyle
i\lambda_{-}(\alpha^{*}\beta-\alpha\beta^{*})+i\lambda_{+}(\alpha\beta-\alpha^{*}\beta^{*})$
$\displaystyle-2\Gamma_{\downarrow}\left(\frac{1}{2}+\gamma\right)+2\Gamma_{\uparrow}\left(\frac{1}{2}-\gamma\right),$
(53)

where $\alpha=\braket{\hat{a}}/\sqrt{N}$. The adiabatic elimination procedure
asserts that the cavity field amplitude $\alpha$ is effectively in the steady
state ($d\alpha/dt=0$) when the system is dissipative
$\kappa\gg\omega,\omega_{0},\lambda_{-},\lambda_{+}$. The cavity field is then
determined solely by the spin dynamics.

Figure 19:  Comparison of semiclassical dynamics of the Dicke model (blue) and
the LMG model (red) in the superradiant phase (a) and the counter-lasing phase
(b). Parameters are as in Fig. 11, with $\lambda_{+}=1.5$ for (a) and
$\lambda_{+}=2$ for (b).

In Fig. 19 we show the results of numerical integration of Eqs. (51)–(53) and
Eqs. (10)–(11) for the same parameter values. Figure 19(a) shows the
semiclassical Dicke and LMG model trajectories as they approach a superradiant
equilibium point. The LMG model trajectory slightly underestimates the
effective decay rate into the superradiant state and, as a result, converges
to the state more slowly. However and crucially, the two trajectories reach
the same steady-state. Figure 19(b) shows the trajectories converging to the
counter-lasing phase. Again, the initial transients of the two trajectories
differ slightly: the Dicke model trajectory accumulates a larger phase which
causes the LMG model trajectory to slightly lag behind the Dicke model
trajectory once the counter-lasing periodic orbit is reached. However, the
periodic orbits of the two are essentially the same. Overall, there is
qualitative agreement between the two models in this parameter regime, and
good quantitative agreement of the long-term dynamics.

##  Appendix C Lin’s Method

We now briefly outline the numerical implementation of _Lin’s method_ [53, 55]
we used to find homoclinic orbits and EtoP connections at specific values of
$\lambda_{+}$. It is based on the continuation of a suitable two-point
boundary value problem defining two orbit segments, which can be implemented
and solved within the software package auto-07p [44, 45].

An orbit segment with an associated integration time $T$ is a solution of the
time-rescaled differential equation

$\frac{d\mathbf{x}}{dt}=T\vec{f}(\mathbf{x},\lambda_{+}).$ (54)

Here the function $f$ is given by (the real and imaginary parts of) the right-
hand side of Eqs. (10)–(11), $\mathbf{x}(t)=[b_{x}(t),b_{y}(t),\gamma(t)]$,
and $\lambda_{+}$ is the continuation parameter. In this formulation, any
orbit segment is represented over the unit time interval as

$\mathbf{x}:[0,1]\rightarrow\mathbb{R}^{3},$ (55)

and any boundary conditions imposed during a contination concern the begin
point $\mathbf{x}(0)$ and/or the end $\mathbf{x}(1)$. The integration time $T$
is solved for during a computation and, in some formulations, it is allowed to
vary.

###  C.1 Finding Homoclinic Orbits

For the computation of homoclinic orbits to $\mathrm{N}_{u}$, we first compute
two orbit segments: $\mathbf{x}_{1}$ with integration time $T_{1}$ in the one-
dimensional unstable manifolds of $\mathrm{N}_{u}$ that ends in a suitable
plane $\Sigma\subset\mathbb{R}^{3}$, and $\mathbf{x}_{2}$ with integration
time $T_{2}$ in the two-dimensional stable manifolds of $\mathrm{N}_{u}$ that
starts in $\Sigma$; for the purposes of the computations here, we choose
$\Sigma=\\{(b_{x},b_{y},\gamma)\in\mathbb{R}^{3}|\gamma=-0.4\\}$. We
subsequently continue both orbit segments in the system parameter
$\lambda_{+}$ until their points $\mathbf{x}_{1}(1)$ and $\mathbf{x}_{2}(0)$
meet in $\Sigma$, meaning that $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ jointly
form (a numerical approximation of) a homoclinic orbit; this is achieved by
requiring that $\mathbf{x}_{1}(1)$ and $\mathbf{x}_{2}(1)$ lie (for any
$\lambda_{+}$) along a fixed direction, given by the _Lin vector_ , which
defines the _Lin gap_ ; see Ref. [55] for further technical details. The
formulation of the respective boundary conditions below assumes knowledge of
the eigenvalues and eigenvectors of $\mathrm{N}_{u}$ as a function of
$\lambda_{+}$.

Figure 20:  Demonstration of Lin’s method to compute a homoclinic orbit. Panel
(a) shows the orbit segments $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ for
$\lambda_{+}=1.5$ after Steps 1A and 1B with a nonzero Lin gap $\varepsilon$
between them. Panel (b) shows them forming the homoclinic orbit
$\mathrm{Hom}^{1}=\mathrm{Hom}_{\mathtt{0}}$ for $\lambda_{+}=1.53286$ from
Fig. 12(a), found in Step 2 when $\varepsilon=0$. Other parameters are as in
Fig. 11.

  1. Step 1A:

Compute, via continuation in $T$, the unstable manifold orbit segment
$\mathbf{x}_{1}\subset W^{u}_{-}(\mathrm{N}_{u})$. The boundary conditions are

$\displaystyle\mathbf{x}_{1}(0)$
$\displaystyle=\mathrm{N}_{u}+\delta_{1}\vec{y}_{u},$ (56)
$\displaystyle\mathbf{x}_{1}(1)$ $\displaystyle\in\Sigma.$ (57)

Condition (56) represents a first-order approximation of the one-dimensional
unstable manifold of $\mathrm{N}_{u}$, where $\vec{y}_{u}$ is its unstable
eigenvector and $\delta_{1}$ is sufficiently small. Initially,
$\mathbf{x}_{1}(1)$ does not lie in $\Sigma$ and the continuation in $T$ stops
when it does and (57) is satisfied.

  2. Step 1B:

Compute, also via continuation in $T$, a stable manifold orbit segment
$\mathbf{x}_{2}\subset W^{s}(\mathrm{N}_{u})$. The boundary conditions are

$\displaystyle\hskip 34.1433pt\mathbf{x}_{2}(0)$ $\displaystyle\in\Sigma,$
(58) $\displaystyle\mathbf{x}_{2}(1)$
$\displaystyle=\mathrm{N}_{u}+\delta_{2}\left[\vec{y}_{s}\cos(\Phi)+\vec{y}_{ss}\sin(\Phi)\right].$
(59)

Condition (59) represents the first-order approximation of the two-dimensional
unstable manifold of $\mathrm{N}_{u}$, spanned by the eigenvectors
$\vec{y}_{s}$ and $\vec{y}_{ss}$ of its two real stable eigenvalues; here the
chosen angle $\Phi$ determines the direction of approach of the orbit segment
to $\mathrm{N}_{u}$ [56]. Again, the continuation in $T$ stops when condition
(58) is satisfied.

  3. Step 2:

Close the Lin gap via continuation in $\lambda_{+}$. We define the Lin vector
as

$\vec{v}\equiv\mathbf{x}_{2}(0)-\mathbf{x}_{1}(1)$

at the beginning of this step and, importantly, then keep $\vec{v}$ fixed
throughout (even though $\mathbf{x}_{1}(0)$ and $\mathbf{x}_{2}(1)$
subsequently change with $\lambda_{+}$). We then continue both orbit segments
$\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ subject to (56), (57), (59), and the
additional boundary condition

$\mathbf{x}_{2}(0)=\mathbf{x}_{1}(1)+\varepsilon\vec{v},$ (60)

which ensures that the Lin direction remains constant and defines the Lin gap
$\varepsilon$. Here $\lambda_{+}$, $T_{1}$, $T_{2}$, $\Phi$, and $\varepsilon$
are free to vary, and the continuation stops when $\varepsilon=0$ is detected,
which means that $\mathbf{x}_{2}(0)=\mathbf{x}_{1}(1)$ and a homoclinic orbit
has been found.

Figure 20 illustrates Lin’s method with the example of finding the homoclinic
orbit $\mathrm{Hom}^{1}=\mathrm{Hom}_{\mathtt{0}}$ from Fig. 12(a). Figure
20(a) shows the situation for $\lambda_{+}=1.5$ after Steps 1A and 1B, when
the two orbit segments $\mathbf{x}_{1}$ from $\mathrm{N}_{u}$ to the section
$\Sigma$ and $\mathbf{x}_{2}$ from $\Sigma$ to $\mathrm{N}_{u}$ have been
found. The Lin gap $\varepsilon$ between the two endpoints
$\mathbf{x}_{1}(1),\mathbf{x}_{2}(0)\in\Sigma$ is initially nonzero, and it is
then closed during Step 2. This zero of $\varepsilon$ is found for
$\lambda_{+}=1.53286$, and Fig. 20(b) shows that then
$\mathbf{x}_{1}(1)=\mathbf{x}_{2}(0)$ and the two orbit segments
$\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ jointly form the homoclinic orbit
$\mathrm{Hom}_{\mathtt{0}}$.

Further intersection points $\mathbf{x}_{1}(1)\in\Sigma$ of
$W^{u}_{-}(\mathrm{N}_{u})$ with the section $\Sigma$ can be found with Step
1A by further continuation in $T$. Each such point defines the associated Lin
gap with $\mathbf{x}_{2}(0)\in\Sigma$, which is then closed in Step 2. In this
way, we computed the other homoclinic orbits
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}}$,
$\mathrm{Hom}_{\mathtt{0}\mathtt{1}\mathtt{0}}$, etc. Note that the symmetric
counterparts of any of these homoclinic orbits, formed by
$W^{u}_{+}(\mathrm{N}_{u})$, are given by the images of $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$ under rotation by $\pi$ about the $\gamma$-axis.

###  C.2 Numerical Continuation of EtoP Connections

The numerical implementation of Lin’s method can also be used to find an EtoP
connection; in fact, it was developed for this purpose in Ref. [55]. We use it
here to find EtoP connections from the equilibrium $\mathrm{N}_{u}$ to either
of the periodic orbits $\mathrm{SRO}_{u}$ when they are of saddle type. To
this end, we consider an orbit segment $\mathbf{x}_{2}$ that now lies in the
two-dimensional stable manifold $W^{s}(\mathrm{SRO}_{u})$. Such an orbit
segment is defined by the boundary condition

$\mathbf{x}_{2}(1)=p_{0}+\delta_{2}\vec{w}_{s},$ (61)

where $p_{0}$ is a point on the periodic orbit ($\mathrm{SRO}_{u}$ with either
$b_{x}<0$ or $b_{x}>0$, depending on the connection sought) and $\vec{w}_{s}$
is the vector at $p_{0}$ of the stable linear Floquet bundle of the periodic
orbit. By replacing boundary condition (59) with (61) in Step 1B, an orbit
$\mathbf{x}_{2}\in W^{s}(\mathrm{SRO}_{u})$ with $\mathbf{x}_{2}(0)\in\Sigma$
is found by continuation in $T$; here $\delta_{2}$ is a chosen small distance
of $\mathbf{x}_{2}(1)$ from the periodic orbit.

The thus found point $\mathbf{x}_{2}(0)$ defines the Lin gap $\varepsilon$
with $\mathbf{x}_{1}(1)$ as before and Step 2 is then used to close it, where
the distance $\delta_{2}$ is free to vary (instead of the angle $\Phi$ as for
homoclinic orbits). We remark that the periodic orbit as well as its stable
bundle must be computed, and continued in Step 2 when the parameter
$\lambda_{+}$ changes. The technical details of the required extended boundary
value problem setup can be found in Ref. [55]. As for homoclinic orbits, we
considered further intersection points $\mathbf{x}_{1}(1)\in\Sigma$ of
$W^{u}_{-}(\mathrm{N}_{u})$ found in Step 1A to find the more complicated EtoP
connections, such as
$\mathrm{EtoP}_{\mathtt{0}\mathtt{1}\mathtt{1}\bar{\mathtt{0}}}$; their
symmetric counterparts are again given by rotation of $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$ by $\pi$ about the $\gamma$-axis.

## References

  * Muniz  _et al._ [2020] J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J. Young, J. R. K. Cline, A. M. Rey, and J. K. Thompson, “Exploring dynamical phase transitions with cold atoms in an optical cavity,”  _Nature_ ,  580, 602 (2020).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, “Ultrastrong coupling between light and matter,”  _Nat. Rev. Physics_ , 1, 19 (2019).
  * Ludlow  _et al._ [2015] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, “Optical atomic clocks,”  _Rev. Mod. Phys._ ,  87, 637 (2015).
  * Qu and Rey [2019] C. Qu and A. M. Rey, “Spin squeezing and many-body dipolar dynamics in optical lattice clocks,”  _Phys. Rev. A_ ,  100, 041602 (2019).
  * Jiang  _et al._ [2011] Y. Y. Jiang, A. D. Ludlow, N. D. Lemke, R. W. Fox, J. A. Sherman, L.-S. Ma, and C. W. Oates, “Making optical atomic clocks more stable with $10^{-16}$-level laser stabilization,”  _Nat. Photonics_ ,  5, 158 (2011).
  * Häffner  _et al._ [2008] H. Häffner, C. Roos, and R. Blatt, “Quantum computing with trapped ions,”  _Physics Reports_ ,  469, 155 (2008).
  * Daley [2011] A. J. Daley, “Quantum computing and quantum simulation with group-II atoms,”  _Quantum Inf. Process._ ,  10, 865 (2011).
  * Dicke [1954] R. H. Dicke, “Coherence in Spontaneous Radiation Processes,”  _Phys. Rev._ , 93, 99 (1954).
  * Hepp and Lieb [1973a] K. Hepp and E. H. Lieb, “Equilibrium Statistical Mechanics of Matter Interacting with the Quantized Radiation Field,”  _Phys. Rev. A_ ,  8, 2517 (1973a).
  * Hepp and Lieb [1973b] K. Hepp and E. H. Lieb, “On the Superradiant Phase Transition for Molecules in a Quantized Radiation Field: the Dicke Maser Model,”  _Ann. Physics_ ,  76, 360 (1973b).
  * Kirton  _et al._ [2019] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, “Introduction to the Dicke Model: From Equilibrium to Nonequilibrium, and Vice Versa,”  _Adv. Quantum Technol._ ,  2, 1970013 (2019).
  * Baumann  _et al._ [2010] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, “Dicke quantum phase transition with a superfluid gas in an optical cavity,”  _Nature_ ,  464, 1301 (2010).
  * Baumann  _et al._ [2011] K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger, “Exploring Symmetry Breaking at the Dicke Quantum Phase Transition,”  _Phys. Rev. Lett._ ,  107, 140402 (2011).
  * Hamner  _et al._ [2014] C. Hamner, C. Qu, Y. Zhang, J. Chang, M. Gong, C. Zhang, and P. Engels, “Dicke-type phase transition in a spin-orbit-coupled Bose-Einstein condensate,”  _Nat. Commun._ ,  5, 1 (2014).
  * Klinder  _et al._ [2015] J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, “Dynamical phase transition in the open Dicke model,”  _Proc. Natl. Acad. Sci. USA_ ,  112, 3290 (2015).
  * Safavi-Naini  _et al._ [2018] A. Safavi-Naini, R. J. Lewis-Swan, J. G. Bohnet, M. Gärttner, K. A. Gilmore, J. E. Jordan, J. Cohn, J. K. Freericks, A. M. Rey, and J. J. Bollinger, “Verification of a Many-Ion Simulator of the Dicke Model Through Slow Quenches across a Phase Transition,”  _Phys. Rev. Lett._ ,  121, 040503 (2018).
  * Dimer  _et al._ [2007] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael, “Proposed realization of the Dicke-model quantum phase transition in an optical cavity QED system,”  _Phys. Rev. A_ ,  75, 013804 (2007).
  * Zhiqiang  _et al._ [2017] Z. Zhiqiang, C. H. Lee, R. Kumar, K. J. Arnold, S. J. Masson, A. S. Parkins, and M. D. Barrett, “Nonequilibrium phase transition in a spin-1 Dicke model,”  _Optica_ , 4, 424 (2017).
  * Zhang  _et al._ [2018] Z. Zhang, C. H. Lee, R. Kumar, K. J. Arnold, S. J. Masson, A. L. Grimsmo, A. S. Parkins, and M. D. Barrett, “Dicke-model simulation via cavity-assisted Raman transitions,”  _Phys. Rev. A_ ,  97, 043858 (2018).
  * Emary and Brandes [2003a] C. Emary and T. Brandes, “Chaos and the quantum phase transition in the Dicke model,”  _Phys. Rev. E_ , 67 (2003a).
  * Emary and Brandes [2003b] C. Emary and T. Brandes, “Quantum Chaos Triggered by Precursors of a Quantum Phase Transition: The Dicke Model,”  _Phys. Rev. Lett._ ,  90, 044101 (2003b).
  * Bhaseen  _et al._ [2012] M. J. Bhaseen, J. Mayoh, B. D. Simons, and J. Keeling, “Dynamics of nonequilibrium Dicke models,”  _Phys. Rev. A_ , 85 (2012).
  * Keeling  _et al._ [2010] J. Keeling, M. J. Bhaseen, and B. D. Simons, “Collective Dynamics of Bose-Einstein Condensates in Optical Cavities,”  _Phys. Rev. Lett._ ,  105, 043001 (2010).
  * Grimsmo and Parkins [2013] A. L. Grimsmo and A. S. Parkins, “Dissipative Dicke model with nonlinear atom-photon interaction,”  _J. Phys. B: At. Mol. Opt. Phys._ , 46 (2013).
  * Soriente  _et al._ [2018] M. Soriente, T. Donner, R. Chitra, and O. Zilberberg, “Dissipation-Induced Anomalous Multicritical Phenomena,”  _Phys. Rev. Lett._ , 120, 183603 (2018).
  * Stitely  _et al._ [2020a] K. C. Stitely, A. Giraldo, B. Krauskopf, and S. Parkins, “Nonlinear semiclassical dynamics of the unbalanced, open Dicke model,”  _Phys. Rev. Research_ ,  2, 033131 (2020a).
  * Buijsman  _et al._ [2017] W. Buijsman, V. Gritsev, and R. Sprik, “Nonergodicity in the Anisotropic Dicke Model,”  _Phys. Rev. Lett._ ,  118, 080601 (2017).
  * Gutiérrez-Jáuregui and Carmichael [2018] R. Gutiérrez-Jáuregui and H. J. Carmichael, “Dissipative quantum phase transitions of light in a generalized Jaynes-Cummings-Rabi model,”  _Phys. Rev. A_ ,  98, 023804 (2018).
  * Stitely  _et al._ [2020b] K. C. Stitely, S. J. Masson, A. Giraldo, B. Krauskopf, and S. Parkins, “Superradiant switching, quantum hysteresis, and oscillations in a generalized Dicke model,”  _Phys. Rev. A_ ,  102, 063702 (2020b).
  * Shchadilova  _et al._ [2020] Y. Shchadilova, M. M. Roses, E. G. Dalla Torre, M. D. Lukin, and E. Demler, “Fermionic formalism for driven-dissipative multilevel systems,”  _Phys. Rev. A_ ,  101, 013817 (2020).
  * Gelhausen  _et al._ [2017] J. Gelhausen, M. Buchhold, and P. Strack, “Many-body quantum optics with decaying atomic spin states: $(\kappa,\gamma)$ Dicke model,”  _Phys. Rev. A_ ,  95, 063824 (2017).
  * Gelhausen and Buchhold [2018] J. Gelhausen and M. Buchhold, “Dissipative Dicke model with collective atomic decay: Bistability, noise-driven activation, and the nonthermal first-order superradiance transition,”  _Phys. Rev. A_ ,  97, 023807 (2018).
  * Kirton and Keeling [2017] P. Kirton and J. Keeling, “Suppressing and Restoring the Dicke Superradiance Transition by Dephasing and Decay,”  _Phys. Rev. Lett._ ,  118, 123602 (2017).
  * Kirton and Keeling [2018] P. Kirton and J. Keeling, “Superradiant and lasing states in driven-dissipative Dicke models,”  _New J. Phys._ ,  20, 015009 (2018).
  * Morrison and Parkins [2008a] S. Morrison and A. S. Parkins, “Dissipation-driven quantum phase transitions in collective spin systems,”  _J. Phys. B: At. Mol. Opt. Phys._ ,  41, 195502 (2008a).
  * Masson  _et al._ [2017] S. J. Masson, M. D. Barrett, and S. Parkins, “Cavity QED Engineering of Spin Dynamics and Squeezing in a Spinor Gas,”  _Phys. Rev. Lett._ ,  119, 213601 (2017).
  * Morrison and Parkins [2008b] S. Morrison and A. S. Parkins, “Collective spin systems in dispersive optical cavity QED: Quantum phase transitions and entanglement,”  _Phys. Rev. A_ ,  77, 043810 (2008b).
  * Lipkin  _et al._ [1965] H. Lipkin, N. Meshkov, and A. Glick, “Validity of many-body approximation methods for a solvable model: (I). Exact solutions and perturbation theory,”  _Nuclear Physics_ ,  62, 188 (1965).
  * Golubitsky and Schaeffer [1985] M. Golubitsky and D. G. Schaeffer,  _Siingularities and Groups in Bifurcation Theory: Volume I_ (Spring-Verlag, New York, 1985).
  * Gross and Haroche [1982] M. Gross and S. Haroche, “Superradiance: An essay on the theory of collective spontaneous emission,”  _Physics Reports_ ,  93, 301 (1982).
  * Liebscher [2015] S. Liebscher,  _Bifurcation without Parameters_ , Lecture Notes in Mathematics No. 2117 (Springer International Publishing, 2015).
  * Kuznetsov [2004] Y. A. Kuznetsov,  _Elements of Applied Bifurcation Theory_ , 3rd ed., Applied Mathematical Sciences (Springer, 2004).
  * Ferri  _et al._ [2021] F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, “Emerging Dissipative Phases in a Superradiant Quantum Gas with Tunable Decay,”  _Phys. Rev. X_ ,  11, 041046 (2021).
  * Doedel [1981] E. J. Doedel, “AUTO: A program for the automatic bifurcation analysis of autonomous systems,”  _Congr. Numer._ , 30 (1981).
  * Doedel and Oldeman [2010] E. J. Doedel and B. E. Oldeman,  _AUTO-07p: Continuation and bifurcation software for ordinary differential equations_ , Tech. Rep. (Department of Computer Science, Concordia University, Montreal, Canada, 2010) with major contributions from A. R. Champneys, F. Dercole, T. F. Fairgrieve, Y. Kuznetsov, R. C. Paffenroth, B. Sandstede, X. J. Wang, and C. H. Zhang, available at http://cmvl.cs.concordia.ca/auto/.
  * Homburg and Sandstede [2010] A. J. Homburg and B. Sandstede, “Homoclinic and Heteroclinic Bifurcations in Vector Fields,” in  _Handbook of Dynamical Systems_ , Vol. 3, edited by H. Broer, B. Hasselblatt, and F. Takens (Elsevier Science, 2010) pp. 379–524.
  * Shilnikov [2002] L. P. Shilnikov, “Bifurcations and Strange Attractors,” in  _Proceedings of the International Congress of Mathematicians Volume III, Beijing 2002_ (Higher Education Press of China, 2002) pp. 349–372.
  * Barrio  _et al._ [2012] R. Barrio, A. Shilnikov, and L. Shilnikov, “Kneadings, symbolic dynamics and painting Lorenz chaos,”  _Int. J. Bifurc. Chaos_ , 22, 1230016 (2012).
  * Xing  _et al._ [2014] T. Xing, R. Barrio, and A. Shilnikov, “Symbolic Quest into Homoclinic Chaos,”  _Int. J. Bifurc. Chaos_ , 24, 1440004 (2014).
  * Giraldo  _et al._ [2021] A. Giraldo, N. G. R. Broderick, and B. Krauskopf, “Chaotic switching in driven-dissipative Bose-Hubbard dimers: When a flip bifurcation meets a T-point in $\mathbb{R}^{4}$,”  _Discrete Continuous Dyn. Syst. Ser. B_ , Online First (2021).
  * Abad  _et al._ [2012] A. Abad, R. Barrio, F. Blesa, and M. Rodríguez, “Algorithm 924: TIDES, a Taylor Series Integrator for Differential EquationS,”  _ACM Trans. Math. Softw._ , 39 (2012).
  * Dovgoshey  _et al._ [2006] O. Dovgoshey, O. Martio, V. Ryazanov, and M. Vuorinen, “The Cantor function,”  _Expo. Math._ ,  24, 1 (2006).
  * Lin [1990] X.-B. Lin, “Using Melnikov’s method to solve Silnikov’s problems,”  _Proc. Roy. Soc. Edinburgh Sect. A_ ,  116, 295 (1990).