Dissipation-induced instabilities of a spinor Bose-Einstein condensate inside an optical cavity

TL;DR

This paper explores how dissipation causes instabilities in a spinor Bose-Einstein condensate within an optical cavity, revealing non-reciprocal dynamics and phase transitions driven by photon losses.

Contribution

It introduces a detailed analysis of dissipative instabilities in a two-component Dicke model with complex couplings, highlighting the role of non-reciprocal interactions and exceptional points.

Findings

Identification of dynamical limit cycles and instabilities

Dissipation-induced non-reciprocal couplings between spins

Potential for experimental observation of non-equilibrium phenomena

Abstract

We investigate the dynamics of a spinor Bose-Einstein condensate inside an optical cavity, driven transversely by a laser with a controllable polarization angle. We focus on a two-component Dicke model with complex light-matter couplings, in the presence of photon losses. We calculate the steady-state phase diagram and find dynamical instabilities in the form of limit cycles, heralded by the presence of exceptional points and level attraction. We show that the instabilities are induced by dissipative processes which generate non-reciprocal couplings between the two collective spins. Our predictions can be readily tested in state-of-the-art experiments and open up the study of non-reciprocal many-body dynamics out of equilibrium.

Equations

$$\hat{H}=-\Delta {\hat{a}}^{†}\hat{a}+\sum _{m_F}{\omega }_0{\hat{J}}_{z,m_F}+\frac{{\hat{J}}_{x,m_F}}{\sqrt{N_{m_F}}}({\lambda }_{m_F}^{\ast }\hat{a}+{\lambda }_{m_F}{\hat{a}}^{†}),$$

$$\left(\begin{array}{c}\delta {\dot{\hat{J}}}_{x,1}\\ \delta {\dot{\hat{J}}}_{y,1}\\ \delta {\dot{\hat{J}}}_{x,-1}\\ \delta {\dot{\hat{J}}}_{y,-1}\end{array}\right)=\left(\begin{array}{cccc}0& -{\omega }_0& 0& 0\\ {\omega }_0+{\xi }_{+}& 0& {\chi }_{+}& 0\\ 0& 0& 0& -{\omega }_0\\ {\chi }_{-}& 0& {\omega }_0+{\xi }_{-}& 0\end{array}\right)\left(\begin{array}{c}\delta {\hat{J}}_{x,1}\\ \delta {\hat{J}}_{y,1}\\ \delta {\hat{J}}_{x,-1}\\ \delta {\hat{J}}_{y,-1}\end{array}\right).$$

$${\eta }_{\pm ,\pm }=\pm i\sqrt{{\omega }_0({\omega }_0+\xi )\pm {\omega }_0\sqrt{{\chi }_{+}{\chi }_{-}}},$$

$${\chi }_R^{-1}(\omega )\delta \hat{a}(\omega )=\frac{-i\lambda }{2\sqrt{N}}\sum _{\sigma }\left[\delta {\hat{J}}_{\sigma }^{+}(\omega )+\delta {\hat{J}}_{\sigma }^{-}(\omega )\right]e^{i\sigma \varphi }-\sqrt{\kappa }{\hat{a}}_{\text{in}}$$

$$\left(\begin{array}{cc}i({\omega }_0-\omega )+\frac{i}{2}\Sigma ({\omega }_0)& \frac{i}{2}{\Lambda }_{+}({\omega }_0)\\ \frac{i}{2}{\Lambda }_{-}({\omega }_0)& i({\omega }_0-\omega )+\frac{i}{2}\Sigma ({\omega }_0)\end{array}\right)\left(\begin{array}{c}\delta {\hat{J}}_1^{-}\\ \delta {\hat{J}}_{-1}^{-}\end{array}\right)={\hat{\vec{\Gamma }}}_{\text{in}}$$

$${\epsilon }_{\pm }^{-}=-i{\omega }_0-\frac{i}{2}\Sigma ({\omega }_0)\pm \frac{i}{2}\sqrt{{\Lambda }_{+}({\omega }_0){\Lambda }_{-}({\omega }_0)}.$$

$$\hat{H}=$$

$$\displaystyle\times\left(\hat{b}_{\sigma}^{\dagger}\sqrt{N-\hat{b}_{\sigma}^{\dagger}\hat{b}_{\sigma}}+\sqrt{N-\hat{b}_{\sigma}^{\dagger}\hat{b}_{\sigma}}\>\hat{b}_{\sigma}\right)\Big{]}.$$

$$E_{\text{MF}}=$$

$$\displaystyle+\lambda(\alpha e^{-i\sigma\phi}+\alpha^{*}e^{i\sigma\phi})\textrm{Re}[\beta_{\sigma}]\sqrt{1-\frac{|\beta_{\sigma}^{2}|}{N}}\Big{]},$$

$$E_{\text{MF}}\simeq$$

$$+$$

$${\partial }_t\alpha$$

$${\partial }_t{S}_{x,\pm 1}$$

$${\partial }_t{S}_{y,\pm 1}$$

$${\partial }_t{S}_{z,\pm 1}$$

$$S_{x,\pm 1}+\frac{2{\lambda }^2}{N({\Delta }^2+\frac{{\kappa }^2}{4})}\{\Delta S_{x,\pm 1}+[\Delta \cos (2\varphi )$$

$$\mp \frac{\kappa }{2}\sin (2\varphi )]S_{x,\mp 1}\}\sqrt{\frac{N^2}{4}-S_{x,\pm 1}^2}=0,$$

$${\xi }_{\pm }$$

$$+\left[\Delta \cos (2\varphi ))\mp \frac{\kappa }{2}\sin (2\varphi )\right]\frac{S_{x,1}{S}_{x,-1}}{S_{z,\pm 1}^2}\}$$

$${\chi }_{\pm }$$

$$\left(\begin{array}{c}\delta \dot{\hat{a}}\\ \delta {\dot{\hat{a}}}^{†}\\ \delta {\dot{\hat{J}}}_{x,1}\\ \delta {\dot{\hat{J}}}_{y,1}\\ \delta {\dot{\hat{J}}}_{x,-1}\\ \delta {\dot{\hat{J}}}_{y,-1}\end{array}\right)=\left(\begin{array}{cccccc}i\Delta -\frac{\kappa }{2}& 0& -i\frac{\lambda }{\sqrt{N}}{e}^{i\varphi }& 0& -i\frac{\lambda }{\sqrt{N}}{e}^{-i\varphi }& 0\\ -i\Delta -\frac{\kappa }{2}& 0& i\frac{\lambda }{\sqrt{N}}{e}^{-i\varphi }& 0& i\frac{\lambda }{\sqrt{N}}{e}^{i\varphi }& 0\\ 0& 0& 0& -{\omega }_0& 0& 0\\ -\frac{\lambda }{\sqrt{N}}{S}_{z,1}{e}^{-i\varphi }& -\frac{\lambda }{\sqrt{N}}{S}_{z,1}{e}^{i\varphi }& {\omega }_0+{\Xi }_{+}& 0& 0& 0\\ 0& 0& 0& 0& 0& -{\omega }_0\\ -\frac{\lambda }{\sqrt{N}}{S}_{z,-1}{e}^{i\varphi }& -\frac{\lambda }{\sqrt{N}}{S}_{z,-1}{e}^{-i\varphi }& 0& 0& {\omega }_0+{\Xi }_{-}& 0\end{array}\right)\left(\begin{array}{c}\delta \hat{a}\\ \delta {\hat{a}}^{†}\\ \delta {\hat{J}}_{x,1}\\ \delta {\hat{J}}_{y,1}\\ \delta {\hat{J}}_{x,-1}\\ \delta {\hat{J}}_{y,-1}\end{array}\right),$$

Full text

Introduction.— Ultracold atomic gases loaded into optical cavities form an ideal set-up for the study of quantum many-body systems far from equilibrium Ritsch et al. (2013). Their large cooperativity allows reaching the strong light-matter coupling regime Brennecke et al. (2007); Colombe et al. (2007) and cavity photon losses enable in-situ monitoring of the many-body dynamics in real time Mekhov et al. (2007); Chen et al. (2007). A representative example of this idea is the experimental realization of the Dicke superradiant phase transition Dicke (1954); Hepp and Lieb (1973); Wang and Hioe (1973), using the motional degrees of freedom of a Bose-Einstein condensate (BEC) Baumann et al. (2010); Klinder et al. (2015a), which provided access to the observation of critical phenomena Mottl et al. (2012) and driven-dissipative dynamics Brennecke et al. (2013). This was also discussed in several theoretical works Nagy et al. (2010); Keeling et al. (2010); Nagy et al. (2011); Öztop et al. (2012); Bhaseen et al. (2012). Further advances have led to the study of competition between short- and long-ranged interactions using optical lattices Klinder et al. (2015b); Landig et al. (2016); Hruby et al. (2018), the simulation of continuous symmetry breaking with multiple cavities Léonard et al. (2017a, b), and the observation of complex many-body phenomena in multi-mode cavities Kollár et al. (2017); Ballantine et al. (2017); Vaidya et al. (2018).

Recently, considerable progress has been made, both theoretically Mivehvar et al. (2017, 2019) and experimentally Zhiqiang et al. (2017); Landini et al. (2018); Kroeze et al. (2018); Norcia et al. (2018); Davis et al. (2019), on the coupling of multiple internal atomic states to the cavity modes, given its potential for quantum simulation of magnetism Mivehvar et al. (2019) and for quantum-enhanced metrology Norcia et al. (2018); Davis et al. (2019). The focus of these studies has been however on the coherent effects of the coupling, leaving the impact of dissipative processes largely unexplored. Dissipation can have noticeable effects on the properties of many-body systems, such as modifying the nature of phase transitions Diehl et al. (2010); Sieberer et al. (2013), the form of the phase diagram Soriente et al. (2018); Kirton and Keeling (2017), their dynamical evolution Lang and Piazza (2016); Zheng and Cooper (2018); Rodríguez Chiacchio and Nunnenkamp (2018); Buca et al. (2018), or giving rise to topological effects Diehl et al. (2011). Therefore, it is exciting to explore the impact of dissipation on these complex systems.

In this paper we investigate the driven-dissipative dynamics of a spinor BEC composed of two hyperfine states coupled to a single mode of an optical cavity, as experimentally realized in Ref. Landini et al. (2018), see Fig. 1. This can be captured by an open two-component Dicke model with complex light-matter couplings. We uncover the emergence of a novel unstable region that, as we show, is induced by the photon losses. Our results are to be seen in contrast to previous studies where this type of dissipation leads to only minor quantitative changes in the phase diagram Dimer et al. (2007); Nagy et al. (2010); Keeling et al. (2010); Bhaseen et al. (2012); Nagy et al. (2011); Öztop et al. (2012). By adiabatically eliminating the cavity field, we find that the interplay between dissipation and complex coupling results in level attraction between eigenfrequencies and the appearance of anti-damping, with the emergence of instabilities being heralded by the presence of exceptional points in the spectrum. In the unstable region, the anti-damping prevents the system from approaching a stable steady-state fixed point and leads to limit-cycle oscillations in the long-time limit. We trace this complex phenomena back to dissipative processes of the cavity field mediating non-reciprocal interactions between the spins. Going beyond adiabatic elimination, we find cavity fluctuations to generate an additional anti-damping contribution that renders the normal phase unstable. Nevertheless, we show that this contribution remains negligible for typical parameters in the current generation of experiments Landini et al. (2018); Kroeze et al. (2018), allowing for observation of the aforementioned phenomena.

Model.— We consider a gas of ultracold spin-1 atoms forming a BEC inside an optical cavity, see Fig. 1. The atoms are coupled to a single cavity mode via a linearly-polarized laser that pumps the system transversely. The atoms mediate two-photon scattering processes between the cavity and the pump which lead to transitions between the BEC state $|k_{0}\rangle$ and the excited states $|\vec{k}_{\pm,\pm}\rangle=|\pm(\vec{k}_{c}\pm\vec{k}_{p})\rangle$, where $\vec{k}_{c,p}$ are cavity and pump momenta, respectively. We fix $|\vec{k}_{c}|=|\vec{k}_{p}|=k$ and, in this case, all the states $|\vec{k}_{\pm,\pm}\rangle$ are degenerate, thus for each atom $i$ the transitions take place between $|k_{0}\rangle^{i}$ and the symmetric state $|k\rangle^{i}=\frac{1}{2}\sum_{\mu,\nu=\pm}|\vec{k}_{\mu,\nu}\rangle^{i}$ Baumann et al. (2010). This allows for a description of the system in terms of collective spin operators, which, in the rotating frame of the pump, reads ($\hbar=1$) Landini et al. (2018)

[TABLE]

where $\hat{a}$ is the bosonic annihilation operator for the cavity field, $\Delta=\omega_{p}-\omega_{c}$ is the detuning between the cavity $\omega_{c}$ and the pump $\omega_{p}$ frequency. The operator $\hat{J}_{\alpha,m_{F}}=\sum_{i}\hat{\sigma}^{i}_{\alpha,m_{F}}$ is a collective spin operator, where $\hat{\sigma}_{z,m_{F}}=|k\rangle^{i}{}^{i}\langle k|-|k_{0}\rangle^{i}{}^{i}\langle k_{0}|$ and $\hat{\sigma}^{i}_{x,m_{F}}=\frac{1}{2}(|k_{0}\rangle^{i}{}^{i}\langle k|+\textrm{H.c.})$. The level splitting $\omega_{0}$ equals twice the recoil frequency $\omega_{r}=k^{2}/2m$, and $N_{m_{F}}$ is the number of atoms in spin state $m_{F}$. The third term in (1) describes the scattering of a pump photon into the cavity mode which is accompanied by an atomic transition. Misalignment between pump and cavity polarizations induces a non-vanishing vectorial component in the atomic polarizability, so the spin states couple differently to the cavity Landini et al. (2018). The complex light-matter coupling $\lambda_{m_{F}}=|\lambda_{m_{F}}|e^{i\phi_{m_{F}}}$ have modulus $|\lambda_{m_{F}}|=\sqrt{\lambda_{s}^{2}\cos^{2}\varphi+\lambda_{v}^{2}m_{F}^{2}\sin^{2}\varphi}$, where $\lambda_{s,v}$ is proportional to the scalar and vectorial atomic polarizabilities and $\varphi$ the angle between the pump and cavity polarization vectors $\vec{\epsilon}_{p}$ and $\vec{\epsilon}_{c}$, and $\tan\phi_{m_{F}}=\frac{\lambda_{v}m_{F}}{\lambda_{s}}\tan\varphi$ Landini et al. (2018).

For the remainder of this paper, we will focus on the case $N_{\pm 1}=N$, $|\lambda_{\pm 1}|=\lambda$ and $\phi_{1}=-\phi_{-1}=\phi$. We obtain a two-component variant of the Dicke model Dicke (1954); Hepp and Lieb (1973); Wang and Hioe (1973). We stress that, for the Hamiltonian (1), the two effective atomic spins cannot be encapsulated in a single collective spin due to the phase difference between the couplings and that the phase difference $\phi$ cannot be removed from the Hamiltonian by any gauge transformation. Indeed, we find that it is one of the key ingredients for the effects we discuss below.

The Hamiltonian (1) possesses a $\mathbb{Z}_{2}$ symmetry, associated with invariance under the transformation $\hat{\mathcal{U}}=e^{i\pi\hat{\mathcal{N}}}$, with $\hat{\mathcal{N}}=\hat{a}^{\dagger}\hat{a}+\sum_{\sigma=\pm 1}\hat{J}_{z,\sigma}$, which can be understood as parity conservation of the total number of excitations in the system. For $\phi=0$, spontaneous breaking of this symmetry results in the well-known superradiant phase transition of the Dicke model Dicke (1954); Hepp and Lieb (1973); Wang and Hioe (1973), where the global spins acquire a finite and equal $x$-component. For the atom-cavity system, this corresponds to a transition from the BEC state with no photons inside the cavity, corresponding to the normal phase (NP), into a self-organized, density-wave state (DW), accompanied by the emergence of a macroscopic cavity field Nagy et al. (2010); Keeling et al. (2010); Nagy et al. (2011); Öztop et al. (2012); Bhaseen et al. (2012). It has been shown in Ref. Landini et al. (2018) that for $\phi\neq 0$, the spontaneous breaking of the $\mathbb{Z}_{2}$ symmetry can lead to a different kind of superradiant order where the $x$-component of the collective spins anti-align, corresponding to each cloud of atoms self-organizing in opposite checkerboard patterns, i.e. formation of a spin wave (SW). In Fig. 2(a) we show the phase diagram obtained within mean-field theory SM , in agreement with the observations reported in Landini et al. (2018), where the NP-DW as well as the NP-SW boundaries are given by $\lambda_{\textrm{cr}}^{d/s}=\sqrt{\frac{(-\Delta)\omega_{0}}{1\pm\cos{(2\phi)}}}$ and the DW-SW boundary is located at $\phi=\pi/4$. In the following we will study how this phase diagram is modified when taking into account the dissipative nature of the cavity.

Steady-state phase diagram.— We start by including dissipation in our model via a Lindblad master equation of the form $\partial_{t}\hat{\rho}=(-i)[\hat{H},\hat{\rho}]+\kappa(\hat{a}\hat{\rho}\hat{a}^{\dagger}-\frac{1}{2}\hat{a}^{\dagger}\hat{a}\hat{\rho}-\frac{1}{2}\hat{\rho}\hat{a}^{\dagger}\hat{a})$, where $\kappa$ is the photon loss rate. We first focus on the bad-cavity limit, $(|\Delta|,\kappa)\gg(\omega_{0},\lambda)$, as studied experimentally in Refs. Landini et al. (2018); Kroeze et al. (2018). In this limit, the cavity evolves much faster than the atoms, allowing us to adiabatically eliminate the cavity by considering the cavity field amplitude $\alpha=\langle\hat{a}\rangle$ to be in the steady state $\alpha\approx\frac{\lambda}{\sqrt{N}}\frac{S_{x,1}e^{i\phi}+S_{x,-1}e^{-i\phi}}{\Delta+i\frac{\kappa}{2}}$, with $S_{\sigma,\pm 1}=\langle\hat{J}_{\sigma,\pm 1}\rangle$. By setting $\partial_{t}S_{\sigma,\pm 1}=0$ and factorizing higher-order correlations, we obtain a set of algebraic equations for the steady-state solutions SM . To construct the phase diagram, we determine the stability of these solutions by linearizing the equations of motion around the steady state $\hat{J}_{\sigma,\pm 1}(t)\simeq S_{\sigma,\pm 1}+\delta\hat{J}_{\sigma,\pm 1}(t)$

[TABLE]

We observe that the effects of the eliminated cavity field are that of introducing a frequency splitting $\xi_{\pm}$ between the spin components, and inducing effective interactions between the spins $\chi_{\pm}$, with both $\xi_{\pm}$ and $\chi_{\pm}$ functions of the external and order parameters. Note that these are in general different $\chi_{+}\neq\chi_{-}$, resulting in a non-reciprocal coupling. This means that the spins respond differently to the motion of the other one, which turns out to have a strong impact on the driven-dissipative dynamics of the system. The resulting phase diagram is shown in Fig. 2(b), where we identify five different phases, classified by the order parameters $S_{x,\pm 1}$ and the number of stable solutions, which we will now describe in turn.

The most striking difference with the ground-state phase diagram is the emergence of an unstable region inside the NP. To understand this, we look at the spectrum of the dynamical matrix in Eq. (2), which in the NP ($S_{x,\pm 1}=0$) reads

[TABLE]

where $\xi=\xi_{\pm}=\Omega\Delta$ and $\chi_{\pm}=\Omega[\Delta\cos(2\phi)\mp\frac{\kappa}{2}\sin(2\phi)]$, with $\Omega=\lambda^{2}/(\Delta^{2}+\frac{\kappa^{2}}{4})$. From this expression, we see that when the couplings $\chi_{\pm}$ acquire opposite sign, $\tan^{2}{(2\phi)}>4\Delta^{2}/\kappa^{2}$, the imaginary parts (frequencies) coalesce, while the real parts (decay rates) become finite, resulting in the emergence of decay and anti-damping, see Fig. 3(a). This phenomenon is known as level attraction and can only arise in non-Hermitian matrices, such as the dynamical matrix in Eq. (2) Moiseyev (2011). This is signaled by the presence of exceptional points, where the eigenvalues are degenerate and the eigenvectors coalesce. In our case, this corresponds to $\tan^{2}{(2\phi)}=4\Delta^{2}/\kappa^{2}$, where we have $\chi_{\pm}=0$. The emergent anti-damping is what makes the NP unstable. This can be observed in Fig. 4, where we show the time-evolution of the cavity field, from the semi-classical equations of motion SM . Initializing the system in the NP with small fluctuations in the cavity field, we see how the system does not remain in this phase, but becomes unstable and at long times features a limit cycle 111 We have used parameters different from those in the recent experiment Landini et al. (2018) for the phase diagram in Fig. 2(b), so that clean self-sustained oscillations are established within the integration times available through numerical integration. .

We trace the origin of the unstable behavior to the cavity field mediating non-reciprocal interactions. The form of the couplings follows from the equations of motion of the cavity field in the NP, yielding $\chi_{\pm}=-\frac{i}{2}\lambda_{\pm 1}^{*}\lambda_{\mp 1}\chi_{R}+\frac{i}{2}\lambda_{\pm 1}\lambda_{\mp 1}^{*}\chi_{R}^{*}$, with $\chi_{R}=(-i\Delta+\kappa/2)^{-1}$ the cavity response function in the bad-cavity limit. Each term represents an amplitude for a photon scattering process from one spin to another, with $\chi_{\pm}$ being the total scattering amplitude for each pathway. For $\kappa=0$, the total scattering amplitudes are symmetric under the exchange of spins ($\chi_{+}=\chi_{-}$). Conversely, for finite $\kappa$, the phase shift induced by the cavity response results in the interference between the scattering amplitudes being different for each pathway, leading to a non-reciprocal coupling. Thus, we can conclude that emergence of the dynamical instability is a consequence of the dissipative nature of the cavity field. Note that, nevertheless, the phase difference between $\lambda_{\pm 1}$ is a crucial ingredient, as for $\phi=0$ both pathways are equivalent, independently of the value of $\kappa$. This constitutes one of the major findings of this work, which should be contrasted with the impact of photon loss in the standard Dicke model. There, it leads to a shift of the critical point Keeling et al. (2010); Nagy et al. (2011); Öztop et al. (2012); Bhaseen et al. (2012); Dimer et al. (2007) and a change in the critical exponent Nagy et al. (2011), but the ground-state and steady-state phase diagrams are qualitatively similar.

For $\lambda>\lambda^{\textrm{SS}}_{\textrm{cr}}=\sqrt{\frac{(-\omega_{0})\left(\Delta^{2}+\kappa^{2}/4\right)}{\Delta\pm\sqrt{\Delta^{2}\cos^{2}{(2\phi)}-(\kappa^{2}/4)\sin^{2}{(2\phi)}}}}$, the system becomes unstable favoring two different steady-state superradiant phases, which we denote as SDW I and SDW II. These are different from the DW and SW phases in Fig. 2(a) as the effects of dissipation in the steady-state equations leads to $|S_{x,1}|\neq|S_{x,-1}|$, resulting in simultaneous presence of density- and spin-wave order. The SDW I phase is a reminiscent of the DW phase with two of solutions corresponding to the spins being almost aligned and the SDW II phase is instead reminiscent of the SW phase, with a pair of solutions associated with the spins being almost anti-aligned. Finally, we identify a fifth phase at large coupling $\lambda$, where both SDW I and II are steady states of the system. This is analogous to the top middle part in Fig. 2(a) where both the DW and SW phases are local minima of the mean-field energy.

Beyond adiabatic elimination.— Finally, we go beyond adiabatic elimination and include cavity field fluctuations. The steady-state equations for the cavity field $\alpha$ and spins $S_{x,\pm 1}$ remain unchanged. In the linear stability analysis we now have to include the dynamics of the cavity-field fluctuations $\delta\hat{a}$ and $\delta\hat{a}^{\dagger}$, leading to a dynamical matrix that is a 6$\times$6 matrix which does not allow for an analytical expression for the eigenvalues. By solving these equations numerically, we find that the resulting phase diagram is qualitatively similar to the one presented in Fig. 2(b), with the important exception of the NP being unstable for all $\phi\neq 0,\pm\frac{\pi}{2}$ SM . This is due to a pair of complex conjugate eigenvalues $\varepsilon^{\pm}_{+}$ with finite real positive part. In Fig. 3(b), we show the real part of these eigenvalues. These are finite for all $\phi\neq 0,\pm\frac{\pi}{2}$ and they reduce to expression (3) we obtained above in the limit $(|\Delta|,\kappa)\gg(\omega_{0},\lambda)$. As a consequence the system is driven into limit cycles all throughout the region associated with the solutions $S_{x,\pm 1}=0$. This is shown in the inset of Fig. 4, where the time-evolution is considered at a point where adiabatic elimination predicts the NP to be stable, i.e. $\lambda<\lambda_{\textrm{cr}}^{\textrm{SS}}$ and $\tan^{2}{(2\phi)}<4\Delta^{2}/\kappa^{2}$. We observe how the system initially remains in the NP, but at longer times, the system dynamics features limit-cycle behavior.

We investigate this further by calculating the eigenvalues perturbatively, exploiting that this phenomenon is also present at small $\lambda$. In Fourier space, the linearized equations of motion for the cavity fluctuations are

[TABLE]

with $\left[\delta\hat{a}(\omega)\right]^{\dagger}=\delta\hat{a}(-\omega)^{\dagger}$, $\chi_{R}(\omega)=[-i(\Delta+\omega)+\kappa/2]^{-1}$ the cavity response function at finite frequency, $\hat{J}^{\pm}=\hat{J}_{x}\pm i\hat{J}_{y}$ the spin raising and lowering operators and $\hat{a}_{\textrm{in}}$ the cavity input noise. For $\lambda=0$, $\delta\hat{J}^{+}_{\pm 1}$ and $\delta\hat{J}^{-}_{\pm 1}$ rotate with frequencies $\omega_{0}$ and $-\omega_{0}$, respectively. Using a rotating-wave approximation, we can consider these pairs of modes to be effectively uncoupled and focus only on the dynamics of $\delta\hat{J}^{-}_{\pm 1}$. Substituting Eq. (4) into the equation of motion for $\delta\hat{J}^{-}_{\pm 1}(\omega)$, we obtain

[TABLE]

with $\Sigma(\omega)=\frac{\lambda^{2}}{2}[-i\chi_{R}(\omega)+i\chi_{R}^{*}(-\omega)]$ the self-energy and $\Lambda_{\pm}(\omega)=\frac{\lambda^{2}}{2}[-i\chi_{R}(\omega)e^{\mp 2i\phi}+i\chi_{R}^{*}(-\omega)e^{\pm 2i\phi}]$ the non-reciprocal coupling, where in the sprit of Fermi’s Golden Rule we have evaluated the energy-dependent self-energy and coupling at the unperturbed frequency of the mode $\omega_{0}$ Marquardt et al. (2007). We have incorporated all noise terms in $\hat{\vec{\Gamma}}_{\textrm{in}}$. The spectrum follows from the determinant of the dynamical matrix (5) as

[TABLE]

The self-energy $\Sigma(\omega_{0})$ only provides a frequency shift and a finite damping rate, given that for $\Delta<0$, $\textrm{Im}[\Sigma(\omega_{0})]<0$. On the contrary, the couplings $\Lambda_{\pm}(\omega_{0})$ always yield an anti-damping contribution, which cannot be compensated by the self-energy damping, as they emerge in a $\pm$ pair and due to $\textrm{Im}[\Lambda_{+}(\omega_{0})\Lambda_{-}(\omega_{0})]\neq 0$ for all $\phi\neq 0,\pm\frac{\pi}{2}$ and $\kappa\neq 0$. Thus, the finite-frequency response of the cavity fluctuations is responsible for the emergence of anti-damping and for the NP becoming unstable to self-sustained oscillations. As expected, in the limit $\kappa\rightarrow 0$, we obtain $\textrm{Im}[\Lambda_{+}(\omega_{0})\Lambda_{-}(\omega_{0})]\rightarrow 0$ and $\textrm{Re}[\Lambda_{+}(\omega_{0})\Lambda_{-}(\omega_{0})]>0$, thus restoring the stability of the NP and confirming again the dissipative nature of the instability. Interestingly, in this limit, the interaction still remains non-reciprocal $\Lambda_{+}(\omega_{0})\neq\Lambda_{-}(\omega_{0})$, meaning that outside the bad-cavity regime the presence of non-reciprocity does not imply unstable behavior. A second pair of eigenvalues $\varepsilon^{+}_{\pm}$ is obtained if one considers the dynamics of $\delta\hat{J}^{+}_{\pm 1}$ instead, which together with (6) provides an approximate form for the eigenvalues $\varepsilon^{\pm}_{+}$ shown in Fig. 3(b) in the limit of small $\lambda$. The form of (6) also explains why in the bad-cavity limit the instability is confined to a finite region. More specifically, the bad-cavity limit is equivalent to $\omega\rightarrow 0$, corresponding to the zero-frequency response of the cavity fluctuations being the only component playing a role in the dynamics. This leads to $\Sigma(\omega)\rightarrow\xi$ and $\Lambda_{+}(\omega)\Lambda_{-}(\omega)\rightarrow\chi_{+}\chi_{-}$, i.e. the instability occurs if $\chi_{+}\chi_{-}<0$, in agreement with our previous result.

Outlook.— Our work opens exciting avenues for future investigations. First, finding an exact solution similar to Emary and Brandes (2003) or efficient numerics Shammah et al. (2018) would allow one to explore the instability beyond the semi-classical approximation employed here. Second, non-reciprocity has recently been investigated with several platforms Jalas et al. (2013); Lodahl et al. (2017); Bernier et al. (2017); Verhagen and Alù (2017) that have been specifically engineered. Here, it emerges naturally as a consequence of the dissipative nature of the cavity field, offering a testbed for non-reciprocal phenomena in a highly-controlled environment. In particular, interesting directions include the impact of non-reciprocity on higher-order photon correlations Harder et al. (2018) or on synchronization behavior Xu et al. (2014). Third, the effects of interatomic interactions can be investigated with an additional optical lattice and lead to complex many-body behavior Mazzucchi et al. (2016). Finally, following Refs. Landini et al. (2018); Kroeze et al. (2018), we expect the emergence of the unstable regime and the steady-state phase diagram in Fig. 2 to be experimentally observable, by means of photon counting and heterodyne detection.

Note added.— Recently, a preprint Dogra et al. (2019) has reported the observation of the dissipation-induced instability we discuss.

Appendix A Mean-field ground-state phase diagram

In this section we describe the main steps taken to obtain the form of the ground-state phase diagram presented in the main text. This was first discussed in Ref. Landini et al. (2018). We start by introducing the Holstein-Primakoff transformation $\hat{J}_{z,\pm 1}=\hat{b}_{\pm}^{\dagger}\hat{b}_{\pm}-\frac{N}{2}$ and $\hat{J}_{x,\pm 1}=\frac{1}{2}\left(\hat{b}_{\pm}^{\dagger}\sqrt{N-\hat{b}_{\pm}^{\dagger}\hat{b}_{\pm}}+\sqrt{N-\hat{b}_{\pm}^{\dagger}\hat{b}_{\pm}}\>\hat{b}_{\pm}\right)$ and inserting it in the Hamiltonian (1) defined in the main text, which leads to (up to a constant shift)

[TABLE]

We obtain the ground-state phase diagram of the Hamiltonian (A) by studying the mean-field energy

[TABLE]

where $\alpha=\langle\hat{a}\rangle$ and $\beta_{\pm}=\langle\hat{b}_{\pm}\rangle$. Minimizing the energy with respect to both fields yields $\beta_{\pm}=\textrm{Re}[\beta_{\pm}]$ and $\alpha=\frac{\lambda}{\Delta}\sum_{\sigma}e^{i\sigma\phi}\beta_{\sigma}\sqrt{1-\frac{\beta_{\sigma}^{2}}{N}}$. We expand the square-root terms to first order, since close to the phase transition $\frac{\beta_{\pm}}{N}\ll 1$, and introduce new order parameters $\beta_{d,s}=\frac{1}{\sqrt{2}}(\beta_{+}\pm\beta_{-})$, which signal the presence of density- and spin-wave order, respecitvely. This leads to

[TABLE]

From the global minima of (A), we distinguish three different phases as function of $\lambda$ and $\phi$, see Fig. 2(a) in the main text. First, for $\lambda<(\lambda_{\textrm{cr}}^{d},\lambda_{\textrm{cr}}^{s})$, with $\lambda_{\textrm{cr}}^{d}=\sqrt{\frac{(-\Delta)\omega_{0}}{1+\cos{(2\phi)}}}$ and $\lambda_{\textrm{cr}}^{s}=\sqrt{\frac{(-\Delta)\omega_{0}}{1-\cos{(2\phi)}}}$, both order parameters vanish $\beta_{d,s}=0$, corresponding to the NP, where both collective spins point down in the $z$ direction, meaning that all atoms remain in the BEC state. Second, for $\lambda>\lambda_{\textrm{cr}}^{d}$ and $0\leq\phi\leq\frac{\pi}{4}$, we obtain $\beta_{d}=\pm\sqrt{N(1-(\lambda_{\textrm{cr}}^{d}/\lambda)^{2})}$, while $\beta_{s}=0$. This is associated with both collective spins aligning in the $x$-$z$ plane, meaning that both atomic species self-organize in the same pattern, thus realizing a density-wave state (DW). Lastly, for $\lambda>\lambda_{\textrm{cr}}^{s}$ and $\frac{\pi}{4}\leq\phi\leq\frac{\pi}{2}$, we obtain $\beta_{s}=\pm\sqrt{N(1-(\lambda_{\textrm{cr}}^{s}/\lambda)^{2})}$ and $\beta_{d}=0$, where now the $x$-components of the spins point in opposite directions. This can be understood as the atomic species self-organizing in the opposite checkerboard patterns, resulting in the emergence of a spin-wave state (SW).

Appendix B Steady-state equations

The semi-classical equations of motion for the system in the presence of cavity loss are

[TABLE]

Setting $\alpha\approx\frac{\lambda}{\sqrt{N}}\frac{S_{x,1}e^{i\phi}+S_{x,-1}e^{-i\phi}}{\Delta+i\frac{\kappa}{2}}$ and $\partial_{t}S_{\sigma,\pm 1}=0$, leads to a set of algebraic equations that define the steady-state solutions for the spins

[TABLE]

where we used $S_{y,\pm 1}=0$ and $S_{z,\pm 1}=-\sqrt{\frac{N^{2}}{4}-S_{x,\pm 1}^{2}}$.

The linearized equations of motion in the main text follow from inserting the expansion $\hat{J}_{\sigma,\pm 1}(t)\simeq S_{\sigma,\pm 1}+\delta\hat{J}_{\sigma,\pm 1}(t)$ in the equations of motion (10) and keeping only the terms which are linear in $\delta\hat{J}_{\sigma,\pm 1}(t)$, where the $z$-component is fixed by spin conservation $S_{z,\pm 1}\delta\hat{J}_{z,\pm 1}=-S_{x,\pm 1}\delta\hat{J}_{x,\pm 1}$. We obtain the following expression for the self-energy and the coupling

[TABLE]

Appendix C Analysis beyond adiabatic elimination

Outside the bad-cavity regime, to obtain the phase diagram, we need to perform a stability analysis of the solutions of Eq. (B) taking into account cavity fluctuations $\delta\hat{a}$ and $\delta\hat{a}^{\dagger}$. The linearized equations of motion read

[TABLE]

where $\Xi_{\pm}$ is related to $\xi_{\pm}$, as defined in Eq. (B), through $\Xi_{\pm}=\xi_{\pm}+\frac{2\Delta\lambda^{2}S_{z,\pm 1}}{N(\Delta^{2}+\frac{\kappa^{2}}{4})}$. The eigenvalues of this dynamical matrix cannot be obtained analytically, independently of the values of the order parameters, and thus need to be computed numerically. The eigenvalues $\varepsilon^{\pm}_{+}$ plotted in the main text are obtained by numerically diagonalizing the dynamical matrix setting $S_{x,\pm 1}=0$ and $S_{z,\pm 1}=-\frac{N}{2}$. As in depicted in Fig. 2(b) in the main text, due to the cavity fluctuation dynamics, the unstable region is enhanced, with the NP remaining stable only for $\phi=0$ and $\phi=\pm\frac{\pi}{2}$.