Dissipation Induced Nonstationarity in a Quantum Gas

TL;DR

This paper demonstrates that a driven two-component Bose-Einstein condensate coupled to an optical cavity exhibits unavoidable dynamical instability and non-stationary behavior due to dissipation, leading to phenomena like squeezing, entanglement, and dissipative time crystals.

Contribution

The authors solve the model in the thermodynamic limit, revealing persistent instability and non-stationarity not predicted by mean-field theory, and connect these to dissipative time crystal formation.

Findings

System is always dynamically unstable in the thermodynamic limit.

Higher-order correlations show instability, leading to squeezing and entanglement.

Finite-size perturbation theory confirms non-stationary behavior.

Abstract

Non-stationary long-time dynamics was recently observed in a driven two-component Bose-Einstein condensate coupled to an optical cavity [N. Dogra, et al. arXiv:1901.05974] and analyzed in mean-field theory. We solve the underlying model in the thermodynamic limit and show that this system is always dynamically unstable -- even when mean-field theory predicts stability. Instabilities always occur in higher-order correlation functions leading to squeezing and entanglement induced by cavity dissipation. The dynamics may be understood as the formation of a dissipative time crystal. We use perturbation theory for finite system sizes to confirm the non-stationary behaviour.

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)=\hat{\mathcal{L}}\rho (t):=$$

$$-\mathrm{i}[H,\rho (t)]+\kappa \left(2L\rho (t)L^{†}-\{L^{†}L,\rho (t)\}\right),$$

$$H=\hslash \omega a^{†}a+\hslash {\omega }_0(J_{\mathrm{z},+}+J_{\mathrm{z},-})+$$

$$\frac{\hslash }{\sqrt{N}}[{\lambda }_D(a^{†}+a)(J_{\mathrm{x},+}+J_{\mathrm{x},-})+$$

$$\mathrm{i}{\lambda }_S(a^{†}-a)(J_{\mathrm{x},+}-J_{\mathrm{x},-})],$$

$$J_{+,\pm }$$

$$J_{z,\pm }$$

$$a\to a+\alpha \sqrt{N}{b}_{1,2}\to b_{1,2}-\sqrt{{\beta }_{1,2}}\sqrt{N}.$$

$$H=\omega a^{†}a+{\omega }_0(b_1^{†}{b}_1+b_2^{†}{b}_2)+{\lambda }_D(a+a^{†})(b_1+b_1^{†}+b_2+b_2^{†})+\mathrm{i}{\lambda }_S(a-a^{†})(b_1+b_1^{†}-b_2-b_2^{†}).$$

$${\lambda }_{n_1,n_2}=\mathrm{i}{n}_1{\omega }_0+\frac{2n_2}{\kappa }{\Gamma }^2+\mathcal{O}(\frac{1}{{\kappa }^2}),$$

$$⟨\dot{a}(t)⟩=-\frac{i{\lambda }_D\left(q_1^{\ast }(t)+q_1(t)\right)}{\sqrt{N}}+\frac{{\lambda }_S\left(q_2(t{)}^{\ast }+q_2(t)\right)}{\sqrt{N}}-\kappa ⟨a(t)⟩-i\omega ⟨a(t)⟩,$$

$$⟨\dot{b_j}(t)⟩=-\frac{i{\lambda }_D{q}_0(t)}{\sqrt{N}}+(-1{)}^{j-1}\frac{{\lambda }_S{q}_0^{\prime }(t)}{\sqrt{N}}-i{\omega }_0⟨b_j(t)⟩,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\rho (t)=\hat{\mathcal{L}}\rho (t)=-\mathrm{i}[H,\rho (t)]+\kappa \left(2L\rho (t)L^{†}-\{L^{†}L,\rho (t)\}\right),$$

$$H=\hslash \omega a^{†}a+\hslash {\omega }_0(J_{\mathrm{z},+}+J_{\mathrm{z},-})+\frac{\hslash }{\sqrt{N}}[{\lambda }_D(a^{†}+a)(J_{\mathrm{x},+}+J_{\mathrm{x},-})+\mathrm{i}{\lambda }_S(a^{†}-a)(J_{\mathrm{x},+}-J_{\mathrm{x},-})].$$

$$J_{+,\pm }$$

$$J_{z,\pm }$$

$$H=\omega a^{†}a+{\omega }_0(b_1^{†}{b}_1+b_2^{†}{b}_2)+{\lambda }_D(a+a^{†})(b_1+b_1^{†}+b_2+b_2^{†})+\mathrm{i}{\lambda }_S(a-a^{†})(b_1+b_1^{†}-b_2-b_2^{†}).$$

$$a\to a+\alpha \sqrt{N}{b}_{1,2}\to b_{1,2}-\sqrt{{\beta }_{1,2}}\sqrt{N}.$$

$$H={\lambda }_D\sqrt{\frac{k_1}{2N}}\left(\sqrt{N}{\alpha }^{\ast }+a^{†}+a+\alpha \sqrt{N}\right)\left(-\sqrt{{\xi }_1}\left(\sqrt{N}\sqrt{{\beta }_1{}^{\ast }}+\sqrt{{\beta }_1}\sqrt{N}\right)+b_1\sqrt{{\xi }_1}+\sqrt{{\xi }_1}{b}_1{}^{†}\right)+{\lambda }_D\sqrt{\frac{k_2}{2N}}\left(\sqrt{N}{\alpha }^{\ast }+a^{†}+a+\alpha \sqrt{N}\right)\left(-\sqrt{{\xi }_2}\left(\sqrt{N}\sqrt{{\beta }_2{}^{\ast }}+\sqrt{{\beta }_2}\sqrt{N}\right)+b_2\sqrt{{\xi }_2}+\sqrt{{\xi }_2}{b}_2{}^{†}\right)+i\sqrt{\frac{k_1}{2N}}{\lambda }_S\left(\sqrt{N}{\alpha }^{\ast }+a^{†}-a-\alpha \sqrt{N}\right)\left(-\sqrt{{\xi }_1}\left(\sqrt{N}\sqrt{{\beta }_1{}^{\ast }}+\sqrt{{\beta }_1}\sqrt{N}\right)+b_1\sqrt{{\xi }_1}+\sqrt{{\xi }_1}{b}_1{}^{†}\right)-i\sqrt{\frac{k_2}{sN}}{\lambda }_S\left(\sqrt{N}{\alpha }^{\ast }+a^{†}-a-\alpha \sqrt{N}\right)\left(-\sqrt{{\xi }_2}\left(\sqrt{N}\sqrt{{\beta }_2{}^{\ast }}+\sqrt{{\beta }_2}\sqrt{N}\right)+b_2\sqrt{{\xi }_2}+\sqrt{{\xi }_2}{b}_2{}^{†}\right)+\omega \left(a\sqrt{N}{\alpha }^{\ast }+\alpha N{\alpha }^{\ast }+\alpha \sqrt{N}{a}^{†}+a{a}^{†}\right)+{\omega }_0\left(-b_1\sqrt{N}\sqrt{{\beta }_1{}^{\ast }}+\sqrt{{\beta }_1}N\sqrt{{\beta }_1{}^{\ast }}-\sqrt{{\beta }_1}\sqrt{N}{b}_1{}^{†}+b_1{b}_1{}^{†}-N\right)+{\omega }_0\left(-b_2\sqrt{N}\sqrt{{\beta }_2{}^{\ast }}+\sqrt{{\beta }_2}N\sqrt{{\beta }_2{}^{\ast }}-\sqrt{{\beta }_2}\sqrt{N}{b}_2{}^{†}+b_2{b}_2{}^{†}-N\right),$$

$$\sqrt{\alpha }=\frac{-{\lambda }_S\left(\sqrt{2-{\beta }_1{\beta }_1^{\ast }}\left({\beta }_1^{\ast }+{\beta }_1\right)+\left({\beta }_2^{\ast }+{\beta }_2\right)\sqrt{2-{\beta }_2{\beta }_2^{\ast }}\right)+\mathrm{i}{\lambda }_D\left(\sqrt{2-{\beta }_1{\beta }_1^{\ast }}\left({\beta }_1^{\ast }+{\beta }_1\right)+\left({\beta }_2^{\ast }+{\beta }_2\right)\sqrt{2-{\beta }_2{\beta }_2^{\ast }}\right)}{\sqrt{2}(\kappa +\mathrm{i}\omega )}$$

$$\sqrt{{\alpha }^{\ast }}=\frac{-{\lambda }_S\left(\sqrt{2-{\beta }_1{\beta }_1^{\ast }}\left({\beta }_1^{\ast }+{\beta }_1\right)+\left({\beta }_2^{\ast }+{\beta }_2\right)\sqrt{2-{\beta }_2{\beta }_2^{\ast }}\right)+\mathrm{i}{\lambda }_D\left(\sqrt{2-{\beta }_1{\beta }_1^{\ast }}\left({\beta }_1^{\ast }+{\beta }_1\right)+\left({\beta }_2^{\ast }+{\beta }_2\right)\sqrt{2-{\beta }_2{\beta }_2^{\ast }}\right)}{\sqrt{2}(\kappa -\mathrm{i}\omega )}$$

$$\frac{\sqrt{2}\alpha \left(3{\beta }_1{\beta }_1^{\ast }+{\beta }_1^{\ast }-4\right)\left({\lambda }_S+\mathrm{i}{\lambda }_D\right)+i\left(4{\omega }_0\sqrt{2-{\beta }_1{\beta }_1^{\ast }}{\beta }_1^{\ast }+\sqrt{2}{\alpha }^{\ast }\left(3{\beta }_1{\beta }_1^{\ast }+{\beta }_1^{\ast }-4\right)\left({\lambda }_D+\mathrm{i}{\lambda }_S\right)\right)}{\sqrt{\frac{2-{\beta }_1{\beta }_1^{\ast }}{N}}}=0$$

$$\frac{\sqrt{2}\alpha \left(3{\beta }_1{\beta }_1^{\ast }+{\beta }_1^{\ast }-4\right)\left({\lambda }_S+i{\lambda }_D\right)+i\left(4{\omega }_0\sqrt{2-{\beta }_1{\beta }_1^{\ast }}{\beta }_1^{\ast }+\sqrt{2}\left(3{\beta }_1{\beta }_1^{\ast }+{\beta }_1^{\ast }-4\right){\alpha }^{\ast }\left({\lambda }_D+\mathrm{i}{\lambda }_S\right)\right)}{\sqrt{\frac{2-{\beta }_1{\beta }_1^{\ast }}{N}}}=0$$

$$\frac{\mathrm{i}\sqrt{2}\alpha \left(3{\beta }_2{\beta }_2^{\ast }+{\beta }_2^{\ast }-4\right)\left({\lambda }_D+i{\lambda }_S\right)+\sqrt{2}\left(3{\beta }_2{\beta }_2^{\ast }+{\beta }_2^{\ast }-4\right){\alpha }^{\ast }\left({\lambda }_S+\mathrm{i}{\lambda }_D\right)+4i{\omega }_0\sqrt{2-{\beta }_2{\beta }_2^{\ast }}{\beta }_2^{\ast }}{\sqrt{\frac{2-{\beta }_2{\beta }_2^{\ast }}{N}}}=0$$

$$\frac{\sqrt{2}\alpha \left(3{\beta }_2{\beta }_2^{\ast }+{\beta }_2-4\right)\left({\lambda }_S-\mathrm{i}{\lambda }_D\right)-\mathrm{i}\left(4{\beta }_2{\omega }_0\sqrt{2-{\beta }_2{\beta }_2^{\ast }}+\sqrt{2}\left(3{\beta }_2{\beta }_2^{\ast }+{\beta }_2-4\right){\alpha }^{\ast }\left({\lambda }_D-\mathrm{i}{\lambda }_S\right)\right)}{\sqrt{\frac{2-{\beta }_2{\beta }_2^{\ast }}{N}}}=0$$

$$\frac{1}{4}{f}^2\left(-12{\lambda }_D^2-12{\lambda }_S^2+{\omega }^2+2{\omega }_0^2\right)+\frac{1}{16}f\left[2{\omega }_0^2\left({\omega }^2-6\left({\lambda }_D^2+{\lambda }_S^2\right)\right)-20\omega {\omega }_0\left({\lambda }_D^2+{\lambda }_S^2\right)+36\left({\lambda }_D^2+{\lambda }_S^2\right){}^2+{\omega }_0^4\right]+\frac{1}{64}{\omega }_0^2\left(\omega {\omega }_0-2\left(9{\lambda }_D^2+{\lambda }_S^2\right)\right)\left(\omega {\omega }_0-2\left({\lambda }_D^2+9{\lambda }_S^2\right)\right)+f^3=0,$$

$${\hat{a}}_{0,j}={\hat{a}}_j^L,{\hat{a}}_{0,j}^{\prime }={\widehat{a^{†}}}_j^L-{\widehat{a^{†}}}_j^R,$$

$${\hat{a}}_{1,j}={\widehat{a^{†}}}_j^R,{\hat{a}}_{1,j}^{\prime }={\hat{a}}_j^R-{\hat{a}}_j^L.$$

$$[{\hat{a}}_{\nu ,j}{\hat{a}}_{\mu ,k}^{\prime }]={\delta }_{\nu ,\mu }{\delta }_{j,k}[{\hat{a}}_{\nu ,j}{\hat{a}}_{\mu ,k}]=[{\hat{a}}_{\nu ,j}^{\prime }{\hat{a}}_{\mu ,k}^{\prime }]=0.$$

$$H={\underline{a}}^{†}\cdot \mathbf{H}\cdot \underline{a}+\underline{a}\cdot \mathbf{K}\cdot \underline{a}+{\underline{a}}^{†}\cdot {\mathbf{K}}^{\ast }\cdot {\underline{a}}^{†}+\underline{f}\cdot \underline{a}+{\underline{f}}^{\ast }\cdot {\underline{a}}^{†},$$

$$L=\underline{l}\cdot \underline{a}+\lambda ,$$

$$\hat{\mathcal{L}}=-\mathrm{i}{\hat{H}}^L+\mathrm{i}{\hat{H}}^R+\kappa \left(2{\hat{L}}^L{\widehat{L^{†}}}^R-{\widehat{L^{†}}}^L{\hat{L}}^L-{\hat{L}}^R{\widehat{L^{†}}}^R\right)$$

$$=-\mathrm{i}{\underline{\hat{a}}}_0^{\prime }\cdot \mathbf{H}\cdot {\underline{\hat{a}}}_0+\mathrm{i}{\underline{\hat{a}}}_1^{\prime }\cdot {\mathbf{H}}^{\ast }\cdot {\underline{\hat{a}}}_1+\mathrm{i}{\underline{\hat{a}}}_1^{\prime }\cdot \mathbf{K}\cdot (2{\underline{\hat{a}}}_0+{\underline{\hat{a}}}_1^{\prime })-\mathrm{i}{\underline{\hat{a}}}_0^{\prime }\cdot {\mathbf{K}}^{\ast }\cdot (2{\underline{\hat{a}}}_1+{\underline{\hat{a}}}_0^{\prime })+{\underline{\hat{a}}}_0^{\prime }\cdot \mathbf{M}\cdot {\underline{\hat{a}}}_0-{\underline{\hat{a}}}_1^{\prime }\cdot \mathbf{M}\cdot {\underline{\hat{a}}}_1+\underline{g}\cdot {\underline{a}}^{\prime },$$

$$\hat{\mathcal{L}}=\underline{\hat{b}}\cdot \mathbf{S}\cdot \underline{\hat{b}}-S_0\text{mathbbm}1,$$

$$\mathbf{S}=\left(\begin{array}{cc}0& -\mathbf{X}\\ -{\mathbf{X}}^T& \mathbf{Y}\end{array}\right),$$

$$\mathbf{X}:=\frac{1}{2}\left(\begin{array}{cc}\mathrm{i}{\mathbf{H}}^{\ast }+\mathbf{M}& -2\mathrm{i}\mathbf{K}\\ 2\mathrm{i}{\mathbf{K}}^{\ast }& -\mathrm{i}\mathbf{H}+{\mathbf{M}}^{\ast }\end{array}\right),$$

Full text

Introduction— Quantum systems composed of many degrees of freedom are expected to relax to stationarity in the long-time limit. This basic principle has been the subject of intense theoretical and experimental research in recent years, in both equilibrium and non-equilibrium settings and for both open and closed systems, e.g. Rossini2009 ; noneq . Moreover, relaxation to stationarity was shown to happen on relatively short timescales ETHReview . In closed systems it is mathematically understood by arguing that every observable evolves as ${\langle O(t)\rangle}=\sum_{n,m}e^{{\rm i}\omega_{nm}t}c_{nm}$, where $\omega_{nm}=E_{n}-E_{m}$. The coefficients $c_{nm}$ are determined by the initial state and $E_{n}$ is the eigenvalue of the Hamiltonian $H$ and eigenstate $\mathinner{|{E_{n}}\rangle}$. For generic observables and initial states the eigenfrequencies $\omega_{nm}$ entering into the time-evolution will be dense and incommensurate. This will lead mutual dephasing (destructive interference) and relaxation to a time-independent (stationary) value of ${\langle O(t)\rangle}$ in the long-time limit Barthel ; Esslerreview .

In contrast, non-stationary dynamics of macroscopic systems is ubiquitous in nature. The emergence of such behavior from the underlying laws of quantum mechanics is an important open question and has far-reaching implications. A possible way for achieving this is for dissipation to dampen all but selected equidistantly spaced frequencies $n\omega_{0}$ with $n=1,2,3,\cdots$ thus preventing eigenstate thermalization as discussed in complex This is to be contrasted with the emergent stationarity due to a dense incommensurate spectrum described above. This mechanism underlies quantum synchronization qsynch , dissipative time crystals timecrystal , and may induced more complex long-time dynamics complex and long-range off-diagonal order JoeyPaper . This type of time crystalline behavior is qualitatively different from the standard discrete (or Floquet) time crystal, which is a phase of matter with external time-dependent driving that shows persistent oscillations with a period different from the one of external driving field and may exist both without discrete1 ; discrete2 and with discretdiss1 ; discretdiss2 ; discretdiss3 ; discretdiss4 dissipation. The phenomenon we will be interested in is emergent non-stationarity in model systems without explicit time-dependent external driving (i.e. in a co-rotating frame where the dynamical equations are not explicitly time-dependent), similar to non-stationarity from certain initial states due to many-body scars scars1 , or the formation of time-crystals in the isolated Heiseberg spin chain Marko .

Long-term non-stationary dynamics was observed in a recent experiment studying a two-component Bose-Einstein condensate (BEC) coherently coupled to two different spatial atomic configurations Esslinger , referred to as density mode (DM) and spin mode (SM). The coherent couplings are mediated by photons scattered by the atomic system from a transverse pump field into an optical cavity note . The experiment showed a rich phase diagram with a dissipation induced region of instability characterized by persistent oscillations. The system was analyzed in mean-field theory and excellent agreement with the experiment was obtained note .

In this Letter, we provide a solution to the long-time dynamics of the underlying model in the thermodynamic limit by employing the approach pioneered by Emary and Brandes for studying the quantum fluctuations around the mean-field solution Dicke2 . Our main result is that this system displays non-stationary dynamics for all choices of external parameters. In regions where mean-field theory predicts stability non-stationarity is confined to higher-order correlation functions and counter-intuitively leads to dissipation induced squeezing and entanglement which is not present in the corresponding closed system. We confirm our findings using perturbation theory for a finite system.

Two-component BEC coupled to a cavity— We model the experimental setup studied in note in the Lindblad master equation framework lindblad ; openbook . By moving to a co-rotating frame we eliminate the (simple) time-dependence coming from the external high-frequency driving and our starting point is a time-independent master equation,

[TABLE]

taking,

[TABLE]

where $a$ ($a^{\dagger}$) is the annihilation (creation) operator of the cavity mode, $\omega$ is the detuning between the cavity resonance and the transverse pump field, $J_{\alpha,+}$($J_{\alpha,-}$) are the collective spin operators of the $+$($-$) Zeeman state separated by angular frequency $\omega_{0}$, and $\lambda_{D,S}$ are the coupling strengths of the atomic spins to the cavity mode. The cavity loss is modelled by a single Lindblad operator $L=a$ with rate $\kappa$.

We extend the approach of Dicke1 ; Cirac1 to study this system in the thermodynamic limit with the number of particles in the BECs $N\rightarrow\infty$. More specifically, the approach will allow us to study the quantum fluctuations in the leading order of large $N$ around the mean-field solutions. We begin by performing a Holstein-Primakoff transformation,

[TABLE]

where $b_{1}$ ($b_{2}$) is the bosonic annihilation operator for the $+$ ($-$) BEC. Anticipating instabilities of the ground state with ${\langle J_{z,\pm}\rangle}=-N/2$ and ${\langle J_{+,\pm}\rangle}=0$ as already obtained in a mean-field treatment note we also perform a shift

[TABLE]

We expand the Liouvillian in orders of $N$ keeping only powers higher than [math]. The values of $\alpha,\beta_{1,2}$ are then determined by demanding that the resulting Liouvillian is quadratic. We call $\alpha=\beta_{1,2}=0$ the normal case, and superradiant otherwise. Finite values of $\alpha,\beta_{1,2}$ physically correspond to ${\langle J_{+,\pm}\rangle}$ acquiring a non-zero macroscopic mean field value. Our study does not capture the transient build up of the mean-field values $\alpha,\beta_{1,2}$ and thus describes the long-time dynamics. Note that unphysical solutions to $\alpha,\beta_{1,2}$ must be discarded by hand.

In the normal case we obtain for the Hamiltonian (up to an irrelevant constant shift)

[TABLE]

In the superradiant case the Hamiltonian also contains squeezing terms like $b_{1,2}^{2}$ and $(b^{\dagger}_{1,2})^{2}$ and is given in supp . The Lindblad operators remain the same following a shift that removes the linear terms in ${\hat{\cal L}}$.

Since the resulting Liouvillian is quadratic we solve it exactly using the method of ’third quantization’ 3rdquant ; 3rdquantbos . The details are in supp . Expanding the eigenvalues for large $\kappa$ we get in the first two leading orders,

[TABLE]

where now $n_{1,2}=0,\pm 1,\pm 2,\dots$. In the normal case $\Gamma^{2}=V^{2}\equiv\left(\lambda_{D}^{2}+\lambda_{S}^{2}\right)$ and can be found numerically in the superradiant case.

Eigenvalues with a positive real part are an unphysical consequence of the unbounded bosonic Liouvillian superoperator. The eigenmodes corresponding to these eigenvalues ’blow up’ and signal an instability. This instability arises from the coupling of the lossy cavity to the BECs. We emphasize that the leading imaginary part of the eigenvalues are equally spaced and thus dephasing of the dynamics is prevented Barthel ; Esslerreview . This is in contrast to the related closed Dicke model, which possesses a dense spectrum and is known to exhibit chaos and thermalization Dicke2 (see supp for more details). It demonstrates how dissipative engineering of the spectrum by coupling to the lossy cavity prevents thermalization.

Our results are consistent with those obtained in mean-field theory in note . Fluctuations around the initial BEC state will be amplified on a time scale given by $\kappa/\Gamma^{2}$. The system will dynamically evolve away from the initial state and may show persistent oscillations with frequency $\approx\omega_{0}$. In order to analyse how this instability manifests itself in the dynamics of observables we now move to the Heisenberg picture.

Equations of motions— The quadratic Liouvillian admits a finite closed set of Heisenberg equations of motion for $\vec{a}=(a,b_{1},b_{2},a^{\dagger},b^{\dagger}_{1},b^{\dagger}_{2})$ of the form $\dot{\vec{a}}={\hat{\cal L}}^{\dagger}\vec{a}$. It preserves the Gaussian nature of quantum states and hence the system dynamics is fully determined by the one and two point observables. For the one-point functions in the normal case the expectation values evolve according to

[TABLE]

[TABLE]

where $q_{0}(t)=\left\langle a^{\dagger}(t)\right\rangle+\langle a(t)\rangle$, $q^{\prime}_{0}(t)=\left\langle a^{\dagger}(t)\right\rangle-\langle a(t)\rangle$, $q_{1}(t)=\left\langle b_{1}(t)\right\rangle+\left\langle b_{2}(t)\right\rangle$, $q_{2}(t)=\left\langle b_{1}(t)\right\rangle-\left\langle b_{2}(t)\right\rangle$. After adiabatic elimination of the cavity mode we obtain agreement with the mean-field treatment in note . The superradiant case is treated analogously supp .

We perform a stability analysis of one and two-point correlators in Fig. 1 showing the maximum rate at which a fluctuation around the stationary solutions determined by $\alpha$ and $\beta_{1,2}$ can be exponentially amplified. The normal case shown in Fig. 1a) exhibits only a small region around $\phi=\arctan(\lambda_{S}/\lambda_{D})$ where the one-point correlators are stable. This region increases with decreasing strength of the cavity coupling $V$. The instability of the one-point function is accompanied by oscillations with a frequency around $\omega_{0}$ shown in Fig. 1b) almost everywhere. In contrast, the superradiant case shown in Fig. 1c) is stable for most values of $\phi$ and $\omega$ for the chosen parameters. In Fig. 1d) we see that there is a small part of the phase diagram where neither the normal nor any of the superradiant solutions are stable in the one-point functions. We find that the region of instability around $\phi=90^{\circ}$ and $\phi=0$ increases with decreasing $V$.

The stability analysis for the two-point correlators gives a different picture. They are also unstable for all of the superradiant solutions except in the points where the model reduces to the Dicke model ($\lambda_{D}=0$ or $\lambda_{S}=0$). This is consistent with some of the eigenvalues of the Liouvillian obtained from the rapidities having a positive real part which implies that some of the observables must be unstable.

The instabilities are observable in connected two-point correlation functions of the form ${\langle X_{1}X_{2}\rangle}_{c}:={\langle X_{1}X_{2}\rangle}-{\langle X_{1}\rangle}{\langle X_{2}\rangle}$, where $X_{1,2}$ are BEC observables and we show examples in Fig. S1. These are related to spin squeezing spinsq1 ; spinsq2 (see supp ). We find that the spin squeezing parameter in the $y$-direction is oscillatory and can be made arbitrarily small with suitable initial choice of ${\langle a^{2}(0)\rangle}$ supp , indicating entanglement. The non-trivial behavior of the connected correlation function is a clear indication of beyond mean-field behaviour which is of purely quantum origin. The system thus exhibits entanglement induced by dissipation, which could also have ramifications on quantum information processing applications.

We see in Fig. 3 that even in the phase where the one-point functions relax to stationarity, the cross-correlation ${\langle J_{x,+}J_{x,-}\rangle}_{c}$ and ${\langle J_{x,+}J_{y,-}\rangle}_{c}$ functions are non-stationary.

Thus the apparent contradiction between mean-field theory showing a phase transition from stable to unstable normal solutions and a full quantum treatment always giving non-decaying eigenvalues with finite imaginary part is resolved: the phase transition takes place for one-point correlators only while higher-order correlations are always unstable. Since there exist eigenmodes of the quantum Liouvillian with eigenvalues that have ${\,{\rm Re}\,}(\lambda)\geq 0,{\,{\rm Im}\,}(\lambda)\neq 0$, there will always be some observables for some initial condition that will persistently oscillate.

Perturbation theory for finite system size — The positive real eigenvalues of the Liouvillian in the thermodynamic limit are unphysical and an artefact of the unboundedness of the Liouvillian. Furthermore, since we kept only the linearized quantum fluctuations in the Holstein-Primakoff expansion it is possible that higher order correlators qualitatively modify the main results above. To address these issues we study the system using perturbation theory. We find that the short-time dynamics of a system of any size are well-described by the results of the Holstein-Primakoff approach supp . More specifically, we obtain eigenfrequencies that are integer multiples of $\omega_{0}$ in leading order for the longest lived oscillations. The corresponding eigenmodes are givenin supp and are perturbatively close to the vacuum state of the cavity.

Using simple large deviation arguments devi , it is trivial to argue that the number of photons leaving the cavity should also be low. This is consistent with experimental results in the dynamical instability region (Fig. 2(a) of note ). Thus we recover the eigenfrequencies of the persistent oscillations, but without the unphysical unbounded increase in the expectation values.

Conclusion— By taking the thermodynamic limit we exactly solved a model of a driven two-component Bose-Einstein condensate coupled to an optical cavity undergoing dissipation note . We identified that the system has both normal and superradiant behavior (in which the modes acquire macroscopic mean values) Dicke1 ; Dicke2 . We find long-time oscillations in the model due to existence of eigenvalues with non-negative real part and imaginary parts which are close to integer multiples of ${\rm i}\omega_{0}$. The imaginary part corresponds to the frequencies of the persistent oscillations of observables. The fact that these frequencies are almost equidistant means that eigenstate dephasing Barthel ; Esslerreview of the dynamics is impossible. In contrast, if the cavity were closed, we would have densely spaced and incommensurate frequencies in the Hamiltonian that could mutually dephase and lead to equilibration like in the related Dicke model Haake1 ; Haake2 . (see also limitcycle ; semiclassics for related semiclassical treatment of equilibration).

We thus conclude that we have an example of dissipation induced dynamics in a system that would otherwise equilibrate, akin to a dissipative time crystal timecrystal . The system in the laboratory frame is time-dependent due to the the (simple) high-frequency driving note . Moving to the co-rotating frame leaves us with a time-independent master equation (S1). The persistent oscillations arise in in the rotating frame from the interplay between the external drive inputing energy into the system and the cooling from the cavity loss. Without the driving the cavity would simply be empty and without the dissipation the system would heat and thermalize. A sketch of such a mechanism as a possibility for realizing discrete time crystals was proposed in PRX , but here we have identified an actual physical system with this property. The system we study is better understood in the co-rotating frame as an example of continuous time symmetry breaking induced by dissipation as the period of the oscillatory response in the laboratory frame frame bears no fixed relation to the driving period. In other words, the system does not realize a discrete-time crystalline structure in the laboratory frame.

We discussed evidence of spin squeezing in supp which should imply entanglement spinsq2 . We validated the thermodynamic results by taking the large cavity limit and performing perturbation theory for finite systems. We find in the leading order of this perturbation the dark Hamiltonian criteria of Ref. complex are trivially satisfied. The persistent oscillations may thus be also understood as an example of quantum Zeno dynamics (e.g.Facchi , metastability , Zanardi ).

In addition to the persistent oscillations at frequencies close to $\omega_{0}$ the dynamics also entangles the condensates leading to strong squeezing through the dissipative coupling. These features are not captured by mean field theory which instead assumes the condensates to be decoupled.

In the future we plan to investigate the implications of dissipation induced non-stationary squeezing and entanglement e.g. for quantum enhanced metrology carlos .

Note Added: While this manuscript was under preparation, the article andreas appeared which discusses extensions of the mean-field results of note . We go beyond both and study the quantum model.

Acknowledgments— We thank N. Dogra and T. Esslinger for both sharing their preliminary results with us and for useful discussion. We also thank for A. Dietrich and C. Sanchez Munoz for useful discussions. The work here has been supported by EPSRC programme grant EP/P009565/1 and the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 319286 (Q-MAC)..

I Holstein-Primakoff transform

Let us recall the model studied in the main text. We study the dynamics of the two-component BEC in a driven lossy cavity with the following Lindblad master equation,

[TABLE]

where,

[TABLE]

and the cavity loss is given by a single Lindblad operator $L=a$. The superoperator ${\hat{\cal L}}$ is the generator of the dynamics. We will be interested in its spectral properties. In this section we will study two possible Holstein-Primakoff transformations. The standard one, corresponding to the normal phase SDicke1 ; SDicke2 ,

[TABLE]

gives, after expanding in $1/N$ and keeping only terms up to $N^{0}$,

[TABLE]

In certain phases the model may acquire macroscopic stationary values, much like in the related Dicke model SDicke1 ; SDicke2 . To treat this we will follow SDicke1 ; SDicke2 and perform a shift prior to doing the expansion in $1/N$,

[TABLE]

Following the Holstein-Primakoff transformation (S4) and the shift (S6) the Hamiltonian becomes,

[TABLE]

where $k_{i}=2j-\sqrt{\beta^{*}_{i}\beta_{i}}j$, and $\xi_{i}=\sqrt{1-\frac{b^{\dagger}_{i}b_{i}-\sqrt{\beta_{i}j}b^{\dagger}_{i}-\sqrt{\beta^{*}_{i}j}b_{i}}{k_{i}}}$. We then expand again in $1/N$ and keep only the zeroth and higher orders. We arrive at a Hamiltonian that contains both quadratic and linear terms. Likewise, the Lindblad operators now contain constant terms. We must remove these non-quadratic terms by suitable choices of $\alpha$ and $\beta_{i}$ so that we can diagonalize the Liouvillian ${\hat{\cal L}}$ later on S3rdquantbos . The corresponding equations for this can be reduced to the following ones,

[TABLE]

The trivial solution to (S8) $\alpha=\beta_{i}=0$ corresponds to the normal phase from earlier in this section. The other solutions with finite $\beta_{i}$ correspond to the superradiant phase. After setting $\beta_{i}$ and $\alpha$ to the non-trivial solutions of (S8) the Liouvillian ${\hat{\cal L}}$ is again quadratic. In the leading order of ${\cal O}(\frac{1}{\kappa^{0}})$ the solutions to (S8) are always $\alpha^{(0)}=\beta^{(0)}_{1}=\beta^{(0)}_{2}=0$. This implies that both the normal and superradiant cases coincide in the leading order. In particular, the perturbative results for the eigenfrequencies being integer multiples of $\omega_{0}$ (purely imaginary leading order contribution to the rapidities) from the main text hold in the leading order of the superradiant case, as well. Note that certain solutions $\alpha,\beta_{1,2}$ may be unphysical, e.g. correspond to expectation values larger than the maximum spin $N$. Such solutions must be discarded by hand.

II Bogolioubov transformation for the closed system

When the system is closed we may perform a simple Bogolioubov transformation and find that the frequencies are given in terms of roots of a third-order polynomial, which read for the normal case,

[TABLE]

We find analogous results for the superradiant case upon solving (S8).

These roots must be real for the Hamiltonian to be physical. If they are not, we may try applying a shift of the form (S6) prior to the Holstein-Primakoff transformation. If even following the shift the frequencies are not real, then we cannot treat that phase with the Holstein-Primakoff. This is how we find the phase diagram of the closed model in Fig. 1a. Using these equations we also find the spectrum of the Liouvillian for the closed system and show a part of it in Fig. 1b. More precisely, we solve the equations for the rapidities in both the normal and superradiant phases.

We also show a plot of the spectrum of the quantum Liouvillian in the normal phase for both the closed an open case in Fig. 1b. The closed case has a dense set of frequencies on the imaginary axis, whereas the open one does not. This demonstrates the discussion in the Introduction of the main text and presents the basic idea of dissipative spectral engineering.

III Third quantization

In this section we briefly recall the method of third quantization and apply it to our problem. We follow the work of Prosen S3rdquant . More specifically, we will now recall the form of third quantization applied to bosons given by Prosen and Seligman S3rdquantbos .

For sake of brevity we use vector notation $\underline{a}:=(a_{1},a_{2},a_{3})=(a,b_{1},b_{2})$. Superoperators will be denoted with a hat $\hat{O}$. We also define left and right multiplication superoperators, $\hat{O}^{L}x:=Ox$ and $\hat{O}^{R}x:=xO$. We now define the following superoperator maps,

[TABLE]

Crucially, these operator satisfy almost-canonical commutation relations,

[TABLE]

We write the Hamiltonian as,

[TABLE]

where ${\bf H}$ and ${\bf K}$ are $6\times 6$ matrices acting on the space of $(a_{1},a_{2},a_{3},a^{\dagger}_{1},a^{\dagger}_{2},a^{\dagger}_{3})$, and $\underline{f}$ is a 6-dimensional vector. These matrices and this vector can be easily read out of (S5) (or (S7) for the shifted case). We write the Lindblad operator as,

[TABLE]

where $\underline{l}=(1,0,0,0,0,0)$ and $\lambda$ is a constant shift coming from (S6). The full Liouvillian is then,

[TABLE]

where $\underline{\hat{a}}_{0,1}$($\underline{\hat{a}}^{\prime}_{0,1}$) are the corresponding vectors whose components are defined in (S10), ${\bf M}=\kappa(l\otimes l)$, and $\underline{g}$ is a linear shift coming from $\underline{f}$ and $\lambda$ (S12), (S13). In order to remove this linear term and make the Liouvillian fully quadratic we will need to eliminate $\underline{g}$. This is done via solving (S8). Following this the Liouvillian can be written as,

[TABLE]

where $\underline{b}$ is a 12-dimensional vector $\underline{b}:=(\hat{a}_{0},\hat{a}_{1},\hat{a}^{\prime}_{0},\hat{a}^{\prime}_{1})$, and,

[TABLE]

where we further have,

[TABLE]

and,

[TABLE]

The scalar $S_{0}={{\rm tr}}{\bf M}$. Assuming that ${\bf X}$ is diagonalizable,

[TABLE]

where $\chi_{j}$ we call the rapidities. It can be shown that diagonalizing ${\bf X}$ and solving the following matrix equation,

[TABLE]

for ${\bf Z}$ provides the full diagonalization of ${\hat{\cal L}}$ S3rdquantbos . Namely, the Liouvillian can be written in terms of normal master modes,

[TABLE]

as,

[TABLE]

As the normal modes satisfy almost-canonical commutation relations,

[TABLE]

we may construct the entire spectrum of ${\hat{\cal L}}$ from the stationary state $\rho_{\infty}$ (${\hat{\cal L}}\rho_{\infty}=0$),

[TABLE]

where we have slightly abused notation by indexing with $k$ all the possible values of $n_{r}=0,1,2,\ldots$. This follows from the fact that $\hat{\zeta}^{\prime}_{r}$ function as a raising superoperators for ${\hat{\cal L}}$, $[{\hat{\cal L}},\hat{\zeta}_{r}^{\prime}]=\chi_{r}\hat{\zeta}_{r}^{\prime}$. The stationary state itself is Gaussian and can be calculated from the 2-point correlation functions only,

[TABLE]

In the normal phase the rapidities $\chi_{k}$ are given in terms of roots of the following 6th order polynomial,

[TABLE]

Then expanding for large $\kappa$ we arrive to the simple result of the main text for the eigenvalues $\lambda$. Analogous equations for the rapidities may be obtained in the superradiant phase. We show in Fig. S2 how the mean field found in Snote forms from the quantum fluctuations.

It is crucial to note that since the equations of motion for the observables are linear for all one and two-point correlators choosing initial values of all observables to be 0 will lead to them staying 0 for all time. It is for this reason that as the initial state we must pick a non-zero value for ${\langle a(0)\rangle}$, which is physically reasonable in this driven system. We find, however, by integrating the differential equations numerically that the precise value of ${\langle a(0)\rangle}$ only determines the rate of divergence of the unstable observables. We also find that by varying ${\langle a^{2}(0)\rangle}$ we can make the spin squeezing Sspinsq1 ; Sspinsq2 in the y-direction $\xi_{y,\pm}=\frac{N(\Delta J_{y,\pm})^{2}}{{\langle J_{z,\pm}^{2}\rangle}+{\langle J_{x,\pm}^{2}\rangle}}$ arbitrarily small.

IV Semiclassical limit

Here we give the semiclassical approximation of the master equation from the main text. We start from the Heisenberg picture for ${\langle a(t)\rangle}=\alpha,{\langle J_{-\pm}(t)\rangle}=\beta_{1,2},{\langle J_{z,\pm}(t)\rangle}=w_{1,2}$, which we assume factor ${\langle O_{1}O_{2}\rangle}={\langle O_{1}\rangle}{\langle O_{2}\rangle}$ (noting that ${\langle O^{\dagger}\rangle}={\langle O\rangle}^{*}$),

[TABLE]

The crucial difference to the mean-field approximation used in Snote is that the semiclassical equations allow for dynamics in the ${\langle J_{z,\pm}\rangle}$, too.

V Perturbation theory for finite system size

As in the main text we take large $\kappa$, i.e., $\kappa=\gamma\kappa^{\prime}$ and that $\gamma\gg 1$, $\kappa^{\prime}\approx\omega_{0}$. This allows us to do perturbation theory. We split ${\hat{\cal L}}$ into ${\hat{\cal L}}=\gamma({\hat{\cal L}}^{(0)}+1/\gamma{\hat{\cal L}}^{(1)})$, where,

[TABLE]

We are interested in both the right and left eigenvectors of ${\hat{\cal L}}$,

[TABLE]

We then formally expand $\rho=\rho^{(0)}+1/\gamma\rho^{(1)}+1/\gamma^{2}\rho^{(2)}+{\cal O}(1/\gamma^{3})$ and likewise for $\lambda$ and $\sigma$. We insert this into the eigenvalue equations (S33) and collect the same orders in $\gamma$. The leading order for $\lambda^{(0)}=0$,

[TABLE]

,

[TABLE]

where $A_{\pm}=\mathbbm{1}_{c}\otimes J_{+,\pm}$, and $J_{+,\pm}:=1/2(J_{x,\pm}+{\rm i}J_{y,\pm})$ (with the subscript $c$ denoting the cavity part of the Hilbert space). The NESS subspace is highly degenerate, as the $A_{\pm}$ operators Scomplex trivially commute with both $L$ and $H^{(0)}$ SBucaProsen . The left stationary states are solved similarly with,

[TABLE]

In the next order we find that the degenerate eigenvalue 0 is split in analogy with standard perturbation theory. By acting with $(\sigma^{(0)}_{\infty,n_{\pm},m_{\pm}})^{\dagger}$ from the left and using biorthogonality ${{\rm tr}}(\sigma^{\dagger}_{j}\rho_{k})=\delta_{j,k}$.,

[TABLE]

This may be easily evaluated and gives that the longest lived oscillating observables will oscillate with eigenfrequencies that are integer multiples of $\omega_{0}$ in the leading order. The fact that if the perturbation is purely Hamiltonian, the eigenvalues coming from the splitting of eigenvalue 0 of ${\hat{\cal L}}$ are purely imagainary has been appreciated before, e.g. Smetastability ; SZanardi . The real correction to these imaginary eigenvalues will be of order $1/\gamma$ smaller than the imaginary part. The corresponding eigenmodes are given by (S34) and are perturbatively close to the vacuum state in the cavity.

We are also interested in the statistics of the number of photons leaving the cavity, which is another experimentally relevant parameter Snote . To compute this we apply the method of large deviations Sdevi . Mathematically, this corresponds to introducing a counting field $\chi$ in the Liouvillian that ’counts’ how much photons enter ($e^{+{\rm i}\chi}$), or leave ($e^{-{\rm i}\chi}$) the cavity. The cumulant generating function for the flow of photons will be given by the leading eigenvalue of this ’deformed’ Liouvillian Sdevi .

In the leading order, in ${\hat{\cal L}}^{(0)}$ the only term changing the number of photons is $L$. It is straightforward to show that the counting field does not change the eigenspectrum of ${\hat{\cal L}}^{(0)}$, meaning that there is no flow of photons in the leading order of strong cavity loss.

Both the result for the frequencies of the oscillations and the number of photons leaving the cavity are consistent with the experimental observation in the dynamical instability region (Fig. 2 of note ) where we may assume that $\kappa\gg\omega_{0},\lambda_{d},\lambda_{S}$