Driven-Dissipative Quantum Dynamics in Ultra Long-lived Dipoles in an Optical Cavity

TL;DR

This paper investigates the complex quantum dynamics of many-body atomic arrays in an optical cavity, revealing how interactions influence transient behaviors and oscillations, especially in systems with long-lived atomic states.

Contribution

It demonstrates that interactions alter transient dynamics without changing steady states, highlighting new oscillatory behaviors in driven-dissipative quantum systems with long-lived dipoles.

Findings

Interactions do not change steady states but affect dynamics.

Transient oscillations become faster with increasing atom number.

Distinct transient behaviors depend on initial states.

Abstract

We study the quantum dynamics of many-body arrays of two-level atoms in a driven cavity subject to collective decay and interactions mediated by the cavity field. We work in the bad cavity limit accessible, for example, using long-lived electronic clock states of alkaline earth atoms, for which the bare atomic linewidth is much less than the cavity linewidth. In the absence of interactions, our system reduces to previously studied models of collective fluorescence. We show that while interactions do not qualitatively change the steady state properties, they lead to a drastic change in the dynamical properties. We find that, for some interval of driving strengths, the system shows two very distinct types of transient behaviors that depend on the initial state of the system. In particular, there is a parameter regime where the system features oscillatory dynamics with a period of oscillation that becomes much shorter than the duration of the overall transient dynamics as the atom number increases. We use both mean field and exact numerical calculations of the quantum system to investigate the dynamics.

Figures

Equations

$$\frac{d\hat{\rho}}{dt}=-i\bigg{[}g\Big{(}\hat{a}\hat{J}^{+}+\hat{a}^{\dagger}\hat{J}^{-}\Big{)}-\Delta\hat{a}^{\dagger}\hat{a},\hat{\rho}\bigg{]}+\kappa\Big{(}\hat{a}\hat{\rho}\hat{a}^{\dagger}-\frac{1}{2}\{\hat{a}^{\dagger}\hat{a},\hat{\rho}\}\Big{)},$$

$$\frac{d\hat{\rho}}{dt}=-i\Big{[}\chi\hat{J}^{+}\hat{J}^{-},\hat{\rho}\Big{]}+\Gamma\bigg{(}\hat{J}^{-}\hat{\rho}\hat{J}^{+}-\frac{1}{2}\big{\{}\hat{J}^{+}\hat{J}^{-},\hat{\rho}\big{\}}\bigg{)},$$

$$\begin{array}{rl}\chi & =\frac{4g^2\Delta }{4{\Delta }^2+{\kappa }^2},\\ \Gamma & =\frac{4g^2\kappa }{4{\Delta }^2+{\kappa }^2}.\end{array}$$

$$\frac{d\hat{\rho}}{dt}=-i\Big{[}\chi\hat{J}^{+}\hat{J}^{-}+\Omega\hat{J}_{x},\hat{\rho}\Big{]}+\Gamma\bigg{(}\hat{J}^{-}\hat{\rho}\hat{J}^{+}-\frac{1}{2}\big{\{}\hat{J}^{+}\hat{J}^{-},\hat{\rho}\big{\}}\bigg{)}=\mathcal{L}\hat{\rho}.$$

$$\hat{H}=\chi\hat{J}^{+}\hat{J}^{-}+\Omega\hat{J}_{x}=\chi\hat{J}^{2}-\chi\hat{J}_{z}^{2}+\chi\hat{J}_{z}+\Omega\hat{J}_{x}\approx\chi\frac{N}{2}\Big{(}\frac{N}{2}+1\Big{)}-\chi\hat{J}_{z}^{2}+\Omega\hat{J}_{x},$$

$$\alpha =\frac{i\Omega }{\Gamma -2i\chi },$$

$$-i\chi\bigg{(}\tilde{J}^{+}\tilde{J}^{-}\hat{\rho}_{ss}-\hat{\rho}_{ss}\tilde{J}^{+}\tilde{J}^{-}\bigg{)}+\Gamma\bigg{(}\tilde{J}^{-}\hat{\rho}_{ss}\tilde{J}^{+}-\frac{1}{2}\tilde{J}^{+}\tilde{J}^{-}\hat{\rho}_{ss}-\frac{1}{2}\hat{\rho}_{ss}\tilde{J}^{+}\tilde{J}^{-}\bigg{)}=0.$$

$${\hat{\rho }}_{ss}\propto \frac{1}{{\tilde{J}}^{-}}\frac{1}{{\tilde{J}}^{+}}$$

$$\hat{\rho}_{ss}=\mathcal{C}|\alpha|^{2}\bigg{(}\frac{1}{\hat{J}^{-}-\alpha^{*}\hat{I}}\bigg{)}\bigg{(}\frac{1}{\hat{J}^{+}-\alpha\hat{I}}\bigg{)},$$

$${⟨{\hat{J}}_z⟩}_{ss}=\{\begin{array}{ll}-\frac{N}{2}\sqrt{1-\frac{{\Omega }^2}{{\Omega }_c^2}}& \Omega <{\Omega }_c\\ \text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}0& \Omega >{\Omega }_c,\end{array}$$

$${\Omega }_c=\frac{N}{2}\sqrt{{\Gamma }^2+4{\chi }^2}.$$

$$\xi^{2}=\underset{\mathbf{n_{\perp}}}{\mathrm{min}}\bigg{(}\frac{N\braket{\Delta\hat{J}_{\mathbf{n_{\perp}}}^{2}}}{|\braket{\mathbf{\hat{J}}}|^{2}}\bigg{)},$$

$$Q(\theta ,\varphi )=\frac{(N+1)}{4\pi }⟨\theta ,\varphi |\hat{\rho }|\theta ,\varphi ⟩,$$

$$\begin{array}{rl}\frac{d⟨{\hat{J}}_x⟩}{dt}& =\chi ⟨\{{\hat{J}}_y,{\hat{J}}_z\}⟩+\frac{\Gamma }{2}⟨\{{\hat{J}}_x,{\hat{J}}_z\}⟩-\frac{\Gamma }{2}⟨{\hat{J}}_x⟩,\\ \frac{d⟨{\hat{J}}_y⟩}{dt}& =-\Omega ⟨{\hat{J}}_z⟩-\chi ⟨\{{\hat{J}}_x,{\hat{J}}_z\}⟩+\frac{\Gamma }{2}⟨\{{\hat{J}}_y,{\hat{J}}_z\}⟩-\frac{\Gamma }{2}⟨{\hat{J}}_y⟩,\\ \frac{d⟨{\hat{J}}_z⟩}{dt}& =\Omega ⟨{\hat{J}}_y⟩-\Gamma ⟨{\hat{J}}_x^2+{\hat{J}}_y^2⟩-\Gamma ⟨{\hat{J}}_z⟩.\end{array}$$

$$\begin{array}{rl}{\dot{s}}_x& =N\chi s_y{s}_z+\frac{N\Gamma }{2}{s}_x{s}_z,\\ {\dot{s}}_y& =-\Omega s_z-N\chi s_x{s}_z+\frac{N\Gamma }{2}{s}_y{s}_z,\\ {\dot{s}}_z& =\Omega s_y-\frac{N\Gamma }{2}(s_x^2+s_y^2),\end{array}$$

$$\displaystyle\begin{split}s_{z}(t)&=-\dot{b}(t),\\ s_{+}(t)&=s_{x}(t)+is_{y}(t)=s_{+}(0)e^{-\big{(}\frac{N\Gamma}{2}-iN\chi\big{)}b(t)}+\frac{2i\Omega}{N\Gamma-2iN\chi}\bigg{(}1-e^{-\big{(}\frac{N\Gamma}{2}-iN\chi\big{)}b(t)}\bigg{)}.\end{split}$$

$$\begin{array}{rl}b(0)& =0,\\ \dot{b}(t{)}^2& +V(b(t))=1.\end{array}$$

$$\displaystyle\begin{split}V(b)=\frac{4\frac{\Omega^{2}}{N^{2}}}{\Gamma^{2}+4\chi^{2}}\bigg{(}1&-2e^{-\frac{N\Gamma}{2}b}\cos(N\chi b)+e^{-N\Gamma b}\bigg{)}+|s_{+}(0)|^{2}e^{-N\Gamma b(t)}\\[8.0pt] &+2|s_{+}(0)|\sqrt{\frac{4\frac{\Omega^{2}}{N^{2}}}{\Gamma^{2}+4\chi^{2}}}e^{-\frac{N\Gamma b}{2}}\bigg{(}\cos(\phi-N\chi b)+e^{-\frac{N\Gamma b}{2}}\cos(\phi)\bigg{)},\end{split}$$

$$\phi=\arg\bigg{(}\frac{2i\Omega}{N\Gamma-2iN\chi}s_{+}^{*}(0)\bigg{)}.$$

$$\displaystyle\begin{split}b(0)&=0,\\ \dot{b}(t)^{2}&+\frac{4\frac{\Omega^{2}}{N^{2}}}{\Gamma^{2}+4\chi^{2}}\bigg{(}1-2e^{-\frac{N\Gamma}{2}b(t)}\cos(N\chi b(t))+e^{-N\Gamma b(t)}\bigg{)}=1.\end{split}$$

$$\begin{array}{rl}{s}_{+}^{ss}& =\frac{2i\Omega }{N\Gamma -2iN\chi },\\ s_z^{ss}& =\pm \sqrt{1-\frac{4{\Omega }^2}{N^2{\Gamma }^2+4N^2{\chi }^2}}=\pm \sqrt{1-\frac{{\Omega }^2}{{\Omega }_c^2}},\end{array}$$

$$\lambda=N\bigg{(}\frac{\Gamma}{2}\pm i\chi\bigg{)}s_{z}^{ss}.$$

$$\begin{array}{rl}{s}_z^{ss}& =0,\\ s_y^{ss}& =\frac{N\Gamma }{2\Omega },\\ s_x^{ss}& =\pm \sqrt{1-\frac{N^2{\Gamma }^2}{4{\Omega }^2}},\end{array}$$

$$\lambda^{2}=-\Omega^{2}\Bigg{(}\sqrt{1-\frac{N^{2}\Gamma^{2}}{4\Omega^{2}}}+\operatorname{sgn}(s_{x})\frac{N\chi}{\Omega}\Bigg{)}\sqrt{1-\frac{N^{2}\Gamma^{2}}{4\Omega^{2}}}.$$

$${\bar{s}}_z=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T{s}_z(t)dt,$$

$$\frac{d\hat{\rho}^{E}}{dt}=-i\frac{\chi}{2}\big{[}\hat{J}_{x}^{2},\hat{\rho}^{E}\big{]}+\frac{\Gamma}{2}\big{[}\hat{J}_{x},\big{[}\hat{\rho}^{E},\hat{J}_{x}\big{]}\big{]}+\frac{\Gamma}{4}\big{[}\hat{J}_{z},\big{[}\hat{\rho}^{E},\hat{J}_{z}\big{]}\big{]}+\frac{\Gamma}{4}\big{[}\hat{J}_{y},\big{[}\hat{\rho}^{E},\hat{J}_{y}\big{]}\big{]},$$

$$⟨{\hat{J}}_z⟩=-\frac{N}{2}{e}^{-\frac{3\Gamma t}{4}},$$

$$\braket{\hat{J}_{z}}=-\frac{N}{2}\bigg{(}\cos\frac{\chi t}{2}\bigg{)}^{N-1}\approx-\frac{N}{2}e^{-\frac{N\chi^{2}t^{2}}{8}},$$

$$-\frac{N}{2}{e}^{-\frac{\Gamma {\chi }^2}{48}{N}^2{t}^3}.$$

$$\frac{\partial f(m)}{\partial t}=-\frac{\Gamma}{4}f(m)+\frac{i\chi(2m+1)}{2}f(m)+\frac{\Gamma}{4}\bigg{[}\frac{N}{2}\bigg{(}\frac{N}{2}+1\bigg{)}-m(m+1)\bigg{]}\bigg{[}f(m+1)+f(m-1)-2f(m)\bigg{]}.$$

Full text

I Introduction

Driven-dissipative systems are ubiquitous in nature and are fundamental to modern quantum science. Even in the classical realm non-linear driven dissipative systems can display fascinating behavior. For example, a macroscopic number of intrinsically noisy driven coupled oscillators can exhibit the phenomenon of synchronization Pikovsky et al. (2001), which describes the spontaneous locking of the oscillators to a common phase. At the quantum level their behavior can become extremely complex and intriguing Galve et al. (2017). While dissipation generally tends to degrade quantum correlations, it is now widely appreciated that it can also give rise to many-body physics not possible with strictly coherent dynamics, and can be used explicitly for the creation of entanglement Kraus et al. (2008); Diehl et al. (2008). Harnessing dissipation, and determining its effect on interacting many-body systems is a fundamental step for advancing quantum information and the technological promise of quantum computers and ultra-precise sensors.

Collective spin models are among the most amenable systems to investigate the interplay between dissipation and interactions, since their reduced Hilbert space allows for exact solutions with current computational resources and at the same time their collective nature increases the parameter regime where simple mean field treatments are valid. Recently, collective spin models and generalizations with independent decay and dephasing have regained theoretical interest Dalla Torre et al. (2013); Schneider and Milburn (2002); Morrison and Parkins (2008); Kessler et al. (2012); Iemini et al. (2018); Hannukainen and Larson (2018); Ferreira and Ribeiro (2018); Lang et al. (2018); Maghrebi and Gorshkov (2016); Walls et al. (1978); Lee et al. (2014); Shammah et al. (2017); Foss-Feig et al. (2017) given impressive experimental advances in quantum laboratories. These studies have shown that these systems have steady states with very rich phase diagrams Ferreira and Ribeiro (2018); Walls et al. (1978), can undergo dynamical phase transitions Morrison and Parkins (2008); Žunkovič et al. (2016) including self-organized criticality Helmrich et al. (2018) and bistability Gothe et al. (2018); Ritter et al. (2009); Lugiato (1984), as well as emergent quantum synchronization Bohnet et al. (2012); Zhu et al. (2015); Meiser et al. (2009), and for some parameter regimes in the thermodynamic limit, steady states characterized by oscillating observables Carmichael (1980); Lee et al. (2014); Chan et al. (2015); Buca et al. (2018) in close connection with recently proposed time-crystals Iemini et al. (2018); Gong et al. (2018); Tucker et al. (2018). Spin squeezing has also been a prominent quantity of study in collective systems with strong dissipation Liu et al. (2011); Shen and Duan (2013); Xu et al. (2017); Wolfe and Yelin (2014); Lewis-Swan et al. (2018); Hu et al. (2017).

In this paper we consider a dissipative version of the Lipkin-Meshkov-Glick Lipkin et al. (1965) model that arises in the study of ultra long-lived atoms inside a rapidly decaying QED cavity Norcia et al. (2016), where the bare atomic linewidth is much less than the cavity linewidth. In this limit the rapidly decaying cavity field can be adiabatically eliminated giving rise both to effective elastic interactions between the atoms and to a collective decay process responsible, among other things, for superradiance Norcia et al. (2018). Cavity mediated interactions alone lead to a steady state where all atoms are de-excited, and as such we incorporate an additional coherent drive to pump energy into the system and thus realize a nontrivial steady state. Such a driven-dissipative system has been studied with exact and mean field calculations in the case of zero interactions Drummond and Carmichael (1978); Drummond and Hassan (1980); Drummond (1980); Carmichael (1980), where at a critical value of the driving strength a dissipative phase transition was identified in the steady state. Below the critical driving strength there exists a superradiant phase, which features a steady state with non-zero inversion, whilst above this critical value there is a normal phase characterized by zero net inversion and slowly decaying transient oscillations.

Here, we first show (sec. III) that interactions don’t qualitatively change the properties of the steady state and just induce a renormalization of the critical drive and an overall rotation of the spin observables. In line with other reports in the literature Schneider and Milburn (2002); Dalla Torre et al. (2013); Lee et al. (2014), we also find parameter regimes where the steady state is highly squeezed, scaling with particle number in a similar way to previous studies Lee et al. (2014). However, interactions significantly modify the dynamics (sec. IV). Using a mean field formulation we find that there is a window of driving strengths, which disappears without interactions, in which the transient dynamics is strongly sensitive to the initial conditions. This results in the existence of three distinct regions, instead of two, as compared to the exact steady state or to the dynamics of the non-interacting system (sec. IV.1).

While in general mean field theory describes poorly the long time dynamics, specifically for systems with large enough drive strengths and finite $N$, we find that the mean field phases remain present in the exact solution (sec. IV.2) by integrating numerically the master equation of the system. Accounting for quantum correlations we observe that, in contrast to the mean field predictions, Rabi oscillations do decay (for large enough drives). The decay rate scales much faster with the number of atoms compared to the cases when there is no dissipation (Schliemann (2015); Lerose et al. (2018)) or no interactions. In the latter case it reduces to simple exponential decay. Finally, in sec. V we relate the phase boundaries of the mean field phase diagram with the opening or closing of gaps in the real and imaginary parts of the spectrum of the Liouville operator. In this way, we tie features of the mean field dynamics to properties of the exact quantum solution.

While studies of fluorescence of driven atoms in the presence of dissipation, both collective and independent, have a very long history Drummond and Carmichael (1978); Drummond and Hassan (1980); Drummond (1980); Carmichael (1980) and a myriad results about them have been obtained throughout the years, the investigations presented here are timely since only very recently, atoms featuring ultra long-lived optical transitions, such as the clock transition in strontium, have been cooled down and loaded in optical cavities Norcia et al. (2018, 2016). These capabilities are opening a path for the first experimental observation of dynamical phases in driven dipoles. This is thanks to the extremely long lifetime (larger than a hundred seconds) of the clock transition which allows the system to attain the predicted non-trivial steady states before spontaneous decay destroys correlations between atoms. At the same time, the weak dipole moment slows down the collective dynamics making it to occur at times accessible with standard probes in current cavity QED experiments. This is strikingly different to standard dipole allowed transitions which need much faster probes for detection.

II Model

We consider first a system composed of a single cavity mode interacting collectively with $N$ atoms whose dynamics are determined by a master equation:

[TABLE]

where $\{\cdot,\cdot\}$ denotes the anticommutator, $\hat{a}$ is the destruction operator that acts on the cavity mode, $\hat{J}^{\pm}=\hat{J}_{x}\pm i\hat{J}_{y}$ and $\hat{J}_{x,y,z}=\frac{1}{2}\sum_{i}^{N}\hat{\sigma}_{i}^{x,y,z}$ are collective pseudospin operators, $\hat{\sigma}_{i}^{x,y,z}$ are Pauli spin operators that act on the internal two level system of atom $i$, $\Delta$ is the cavity detuning with respect to the atomic transition frequency, $\kappa$ is the cavity decay rate and $g$ is the atom-light coupling constant.

In the bad cavity limit, accessible with long lived dipoles, the spontaneous decay rate of the atoms, $\gamma$, is much smaller than the cavity decay, $\kappa$ and can be neglected, to an excellent approximation. For example in the case of the clock transition in 87Sr, $\gamma$ is of the order of millihertz while $\kappa$ is typically of the order of kilohertz and thus nine orders of magnitude larger. In the bad cavity limit, $\kappa\gg g$, we can adiabatically eliminate the cavity field from Eq. (1) leading to a collective spin model Bonifacio et al. (1971),

[TABLE]

where

[TABLE]

Here, $\chi$ describes effective elastic exchange interactions and $\Gamma$ describes collective coherent decay of the atomic system, responsible for superradiance. As mentioned before, under Eq. (2) the system reaches a steady state with all atoms de-excited so we add a single particle driving term on resonance with the atomic transition and with a Rabi frequency $\Omega$. The full evolution of the system is thus given by

[TABLE]

The Hamiltonian part, when restricted to the collective manifold, can be reduced to a particular realization of the Lipkin-Meshkov-Glick model:

[TABLE]

where we have neglected $\chi\hat{J}_{z}$ on account of it being a factor of $N$ smaller than the other terms in $\hat{H}$. Here we focus on initial states that are eigenstates of $\hat{J^{2}}$, which is a conserved quantity of the dynamics given by Eq. (4). Consequently, this term only leads to a uniform energy shift in the Hamiltonian. It is important to point out, however, that the $\hat{J}^{2}$ term gives some degree of protection against single particle inhomogeneities or slow non-collective noise sources Norcia et al. (2018). Thus, Eq. (5) can also be understood as a one-axis twisting Hamiltonian in the presence of driving. As a final comment, we remark that theoretical studies of the Lipkin-Meshkov-Glick model typically consider $\chi$ and $\Gamma$ to scale as $\frac{1}{N}$ to keep the energy extensive, in contrast to our case in which $\chi$ and $\Gamma$ are independent of $N$. As a consequence, when the number of atoms increases, we observe dynamics that become faster in real time.

III Steady State Properties

The steady state of the system, $\hat{\rho}_{ss}$ is characterized by $\mathcal{L}\hat{\rho}_{ss}=0$. To solve it, we use the methods of Drummond (1980); Puri and Lawande (1979) and introduce displaced spin operators: $\tilde{J}^{+}=\hat{J}^{+}-\alpha\hat{I}$, where $\hat{I}$ is the identity matrix. If we choose:

[TABLE]

then the master equation simplifies to:

[TABLE]

From it, it can be seen that the ansatz:

[TABLE]

solves independently the Hamiltonian and dissipative parts. Thus, a closed form expression can be written for the system’s density matrix:

[TABLE]

where $\mathcal{C}$ is a normalization factor that ensures $\operatorname{Tr}(\hat{\rho})=1$. This allows us to calculate, for example, the expectation value of $\hat{J}_{z}$ in the large $N$ limit (keeping $\frac{\Omega}{N}$ fixed). The result is:

[TABLE]

where

[TABLE]

Equation (10) shows that there is a phase transition at $\Omega=\Omega_{c}$ for which $\braket{\hat{J}_{z}}_{ss}$ is a good order parameter. We refer to these phases as superradiant $({\mathcal{S}})$ and normal $({\mathcal{N}})$ for $\Omega<\Omega_{c}$ and $\Omega>\Omega_{c}$ respectively (see Fig. 1), consistent with previous studies in the absence of elastic interactions Carmichael (1980).

These works also found phenomena such as spin squeezing Kitagawa and Ueda (1993), which we find persists in the presence of elastic interaction as shown in Fig. 2. There, we use the definition of the squeezing parameter $\xi$ introduced by Wineland et al. Wineland et al. (1992):

[TABLE]

\braket{\mathbf{\hat{J}}}=\big{(}\braket{\hat{J}_{x}},\braket{\hat{J}_{y}},\braket{\hat{J}_{z}}\big{)} is the expectation value of the spin vector, $\braket{\Delta\hat{J}^{2}_{n}}=\braket{(\mathbf{\hat{J}}\cdot\mathbf{n})^{2}}-(\braket{\mathbf{\hat{J}}\cdot\mathbf{n}})^{2}$ is the variance of the projection of the spin in direction $\mathbf{n}$, $\mathbf{n_{\perp}}$ is a direction orthogonal to $\braket{\mathbf{\hat{J}}}$ (that is, $\mathbf{n_{\perp}}\cdot\braket{\mathbf{\hat{J}}}=0$) and we are minimizing the quantity in parentheses in Eq. (12) among all possible $\mathbf{n_{\perp}}$ directions. If $\xi^{2}<1$, then the system is entangled Sørensen et al. (2001). In Fig. 2 we plot $-\log_{10}\xi^{2}$ and it clearly shows that the phases defined by Eq. (10) have very different squeezing properties: it is positive in the ${\mathcal{S}}$ phase (thus entangled) and negative in the ${\mathcal{N}}$ phase. The best squeezing is obtained in a very small region around the transition line and is, quite surprisingly, independent of the interaction strength $\chi$. This should become clear by looking at Eq. (9), from which we can deduce that the squeezing can only be a function of $|\alpha|$ so that any change in $\chi$ can be compensated by a change in $\Omega$ that keeps $|\alpha|$ constant.

To get more insight into the squeezing we calculate the Husimi distribution, defined by

[TABLE]

where $\ket{\theta,\phi}$ is a spin coherent state, defined as the eigenstate of $\mathbf{n}\cdot\mathbf{\hat{J}}$ with positive eigenvalue, where $\mathbf{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$. We show the Husimi distributions of states in different regions of the phase diagram in Fig. 3.

Inside the $\mathcal{S}$ region, the state is very close to being a spin coherent state, with the noise ellipse close to being isotropic. When $\Omega=0$, the steady state is clearly a coherent state pointing in the $-\hat{z}$ direction, so the inclusion of a small driving just lifts this state a little bit but doesn’t change its qualitative properties. The steady state is reached in a situation where the rate at which the state is lifted by $\Omega$ is compensated by the rate at which the state relaxes back to $-\hat{z}$ due to dissipation. If there were no interactions, the state would be pointing in the $yz$ plane (perpendicular to the driving axis), but $\chi$ induces a rotation about the $\hat{z}$ axis. As we increase $\Omega$ and move closer to the transition region, the driving starts winning over the dissipation shearing the noise ellipse in the vertical direction and moving the collective Bloch vector closer to the equator since a larger dissipation is needed to equilibrate the system. Note that dissipation is the greatest at the equator. Finally, when increasing $\Omega$ past the transition point, dissipation is not strong enough to equilibrate the drive at any point in the sphere and the Bloch vector keeps rotating. Dissipation then just decoheres the state, which results in the state spreading over a large area of the Bloch sphere.

IV Dynamics

Given that a simple closed form solution of the exact many-body dynamics is not possible, we instead use a mean field analysis and then compare it to the exact steady state presented in the previous section and the results from exact numerical solutions. The exact equations of motion for the expectation values of the spin observables derived from Eq. (2) are:

[TABLE]

IV.1 Mean field analysis

The mean field approximation assumes the spin variables are uncorrelated and replaces $\braket{\hat{J}_{i}\hat{J}_{k}}\approx\braket{\hat{J}_{i}}\braket{\hat{J}_{k}}$ in Eq. (14). Furthermore, for convenience we define $s_{i}=\frac{\braket{\hat{J}_{i}}}{N/2}$. At the mean field level the equations for $s_{i}$ are then:

[TABLE]

and we have assumed here that $N\Gamma$, $N\chi$ and $\Omega$ are comparable. Regardless of the nonlinear character of the mean field equations, the solution can be parameterized in terms of a variable $b(t)$:

[TABLE]

The first line of Eq. 16 defines $b(t)$, which leaves us free to choose its initial value. The second line of Eq. 16 is obtained through an integration of the equation for $s_{+}(t)$ and setting $b(0)=0$. Thus, $b(t)$ satisfies:

[TABLE]

where

[TABLE]

and

[TABLE]

The time evolution of $b(t)$ can then be understood as the coordinate of a particle of unit mass and energy moving in the effective potential $V(b)$, which clearly depends on the initial conditions. This analogy makes it easier to visualize the solution.

In the case of starting with all the atoms de-excited (that is, the south pole), $s_{+}(0)=0$ and the equations simplify to

[TABLE]

As is shown in Fig. 4, there are two different situations depending on whether $\Omega$ is greater or smaller than some $\Omega_{c}^{1}$ (not equal to the critical value $\Omega_{c}$ derived before for the steady state). For $\Omega<\Omega_{c}^{1}$: $b$ is unbounded, so $\dot{b}$ never changes sign. Thus $s_{z}$ remains negative throughout the evolution and tends to a constant value in the infinite future which agrees with the exact value calculated from the steady state density matrix. For $\Omega>\Omega_{c}^{1}$: $b$ keeps oscillating indefinitely and so does $s_{z}$. There is no equilibration and the system never reaches a steady state.

The importance of $\Omega_{c}^{1}$ is not its precise value but instead two facts: (1) it depends on the initial conditions (since $V(b)$ itself depends on them), and (2) as all different initial conditions are explored, $\Omega_{c}^{1}$ acquires values in the finite interval \big{[}N\Gamma/2,\Omega_{c}\big{]}. As a result, when $\Omega<N\Gamma/2$, the system relaxes to the steady state for any initial condition. When $N\Gamma/2<\Omega<\Omega_{c}$, whether the system relaxes or oscillates indefinitely depends on the initial condition chosen, and when $\Omega_{c}<\Omega$, the system oscillates indefinitely no matter the initial condition (see Fig. 5). To justify this, we computed the possible steady states predicted by Eqs. (15), analyzed the behaviour of small perturbations and found that:

1. For $\Omega<\frac{N\Gamma}{2}$, there is one set of steady states with spin observables given by

[TABLE]

and stability eigenvalues:

[TABLE]

Thus, only the state with $s_{z}^{ss}<0$ is stable and we denote it by $SS_{1}$.

2. For $\frac{N\Gamma}{2}<\Omega<\Omega_{c}$, $SS_{1}$ is still a valid steady state but there is another set of solutions, with spin observables given by

[TABLE]

and stability eigenvalues:

[TABLE]

One of the solutions, which we call $SS_{2}$, is unstable in the region $\frac{N\Gamma}{2}<\Omega<\Omega_{c}$ while the other one ($SS_{3}$) has imaginary eigenvalues. Thus, if the system starts close to $SS_{3}$, then it will keep oscillating indefinitely (as found numerically), and if it starts close to $SS_{1}$ it will decay towards it.

3. For $\Omega_{c}<\Omega$, $SS_{1}$ is no longer a valid steady state, so only $SS_{2}$ and $SS_{3}$ remain. In this case, however, both have imaginary stability eigenvalues, so the system will oscillate, no matter where it starts.

This allows the construction of the phase diagram shown in Fig. 5. For $\Omega<N\Gamma/2$ (blue) we have the superradiant region, where the system always decays to a steady state with a time constant $t_{\mathrm{decay}}\sim N^{-1}$ (see Eq.(22)), hence the name superradiant. For $N\Gamma/2<\Omega<\Omega_{c}$ (red), we have a region of multistability, with more than one possible steady state. Finally for $\Omega>\Omega_{c}$ (yellow), we have a region where the driving term dominates and generates oscillations in the spin inversion that can be understood from a single particle perspective, hence the name normal. Note that the boundary between the red and the yellow regions is the same as the boundary between the superradiant ${\mathcal{S}}$ and the normal ${\mathcal{N}}$ steady state phases, with the yellow phase exactly corresponding to the ${\mathcal{N}}$ phase.

Furthermore, if we do a time average of $s_{z}$,

[TABLE]

as is usually done in the context dynamical phase transitions Žunkovič et al. (2016), we find a jump at $\Omega_{c}^{1}$ (shown in Fig. 6), characteristic of a first order transition. This can be understood as follows: in the decaying region, the average is always going to be dominated by the steady state value of $s_{z}$, so both the average and steady state calculations agree. This does not hold all the way to $\Omega_{c}$, however. Instead, once the oscillating behaviour sets in, $s_{z}$ will average to 0.

The mean field predictions and in particular the oscillatory dynamics in the normal phase seem to contradict the steady state exact solution. In fact, while mean field works well for $\Omega<\frac{N\Gamma}{2}$ and accurately predicts the steady state values of the spin observables, it completely misses the relaxation process for $\Omega>\Omega_{c}$. In between these two values, depending on the initial condition, the mean field treatment may or may not describe the relaxation towards the steady state.

IV.2 Exact solution

To investigate the extent to which the mean field behaviour is correct, we numerically integrate the master equation for $N=400$. We fix $\chi$ and $\Gamma$, vary $\Omega$ and start at the south pole of the Bloch sphere. The resulting plots of inversion, shown in Fig. 7, reveal many interesting features:

In the metastable phase and for drive strengths much less than $\Omega_{c}^{1}$ (A in Fig. 7), mean field offers a good description of the dynamics, including the relaxation towards the steady state, which is shown to occur rather quickly compared to $\Gamma^{-1}\approx 1500$ s. The decay time predicted by mean field is $\sim 10$ s [Eq. (22)], which is faster than observed in Fig. 7. However, the discrepancy is explained by the fact that the system is initialized far from the steady state.

As $\Omega$ is increased, the mean field transition is reached (B and C in Fig. 7), which is signaled in the exact solution by the fact that for $\Omega>\Omega_{c}^{1}$, the inversion $\braket{\hat{J}_{z}}$ crosses 0 at some time during the evolution. Traces of oscillations can also be seen, but the qualitative behaviour of the exact solution is very similar in B and C (in contrast to the mean field predictions). Note, however, that the relaxation towards the steady state is very slow in this region (especially after crossing the boundary towards C).

In the vicinity of $\Omega_{c}$ (D and E in Fig. 7), oscillations have already kicked in, and to a great extent the mean field solution is able to capture the oscillation frequency (as can be seen by looking at the zero crossings). However, there is an envelope that describes relaxation towards $\braket{\hat{J}_{z}}=0$. This is the correct steady state for E since it is already in the $\mathcal{N}$ phase, but for D, which is still in the metastable phase the real steady state is attained after a much longer time (set by $\Gamma$, see Fig. 8, right panel).

Finally, deep in the normal phase (F in Fig. 7), the picture of oscillations modulated by an envelope is accurate. Note that it seems to be independent of $\Omega$ and thus the effect of increasing $\Omega$ is just to increase the oscillation frequency compared to its characteristic decay rate. This envelope is a beyond mean field effect and, as shown in Appendix A, is characterized by $t_{\mathrm{decay}}\sim N^{-2/3}$.

Relaxation in C and D is very slow, so we extend the simulations to long times to be able to observe the decay. The results are shown in Fig. 8. The decay is still exponential, but, in contrast to the superradiant region, there is no $N$-fold enhancement. The decay constants are almost independent of $N$, giving an interesting separation of scales for $\braket{\hat{J}_{z}}$, especially in point $D$: oscillations occur on a timescale of $O(N^{-1})$ (this is the mean field prediction), modulated by an envelope that decays to 0 on a timescale $O(N^{-2/3})$ and then relaxes to its true steady state in a time $\sim O(1)$.

We can also perform a time average. The time window is chosen to be short enough such that $D$ has not yet reached its true quantum steady state but long enough that it has relaxed to its metastable value. We take this characteristic time to be $\sim 200$ s and plot in Fig. 9 the corresponding time average value against the prediction of the mean field and exact steady state (see Fig. 6). In contrast to the first order transition observed at the mean field level, quantum fluctuations soften the sharp jump in the time average, especially at finite $N$. The sharp behaviour is only recovered in the thermodynamic limit $N\rightarrow\infty$.

IV.3 Squeezing dynamics

It is also worthwile to consider the dynamical behaviour of spin squeezing for two reasons: (1) To characterize the amount of spin squeeezing generated during transient dynamics, and (2) to analyze how fast the steady state squeezing described in sec. III can be attained. We show a plot of the spin squeezing as a function of time in Fig. 10, where we have chosen the value of $\Omega$ to be identical to the one that results in the largest steady state squeezing (the fact that it occurs at $\Omega\neq\Omega_{c}$ is a finite size effect). We also show the Husimi distribution of the state at specific points in time to better visualize the nature of the squeezing. At short times, we find a minimum, marked by $\alpha$ in Fig. 10, which rapidly evolves into an oversqueezed state, indicated by $\beta$. This is illustrated by the associated Husimi function for these points in Fig. 10: for $\alpha$, the state remains strongly polarized in one direction with a sheared noise ellipse, whereas for $\beta$ the state wraps around the Bloch sphere. At later times, the distribution relocalizes on the Bloch sphere (point $\gamma$) and slowly relaxes to its steady state at much longer times (notice the timescale in Fig. 10).

V Spectrum of Liouville operator

To better understand the different decay timescales and the onset of oscillatory behaviour, we study the spectrum of the Liouville operator, Eq. (2), and try to find signatures of the mean field transition in its eigenvalues. A gap in the real part of the spectrum translates into exponential decay of observables to their steady state. On the other hand, the imaginary part of the spectrum leads to nontrivial transient oscillatory dynamics.

For a system with $N=200$, in Fig. 11 we plot the real and imaginary parts of the first 150 eigenvalues (ordered ascendingly with respect to their magnitude).

Very clearly, a gap closes in the real part at a value close to $N\Gamma/2$, which implies that the decay times for $\Omega>N\Gamma/2$ can become very long. This is consistent with mean field since this is the region where oscillating nondecaying solutions begin to exist. The imaginary part of the eigenvalues remains gapless, however, indicating that nonoscillating behaviour is also allowed (metastable region). Furthermore, a gap opens in the imaginary part close to $\Omega_{c}$, which implies that observables of the system should generally show oscillatory behaviour for $\Omega>\Omega_{c}$. Once again, this is consistent with the mean field analysis since in this region only oscillatory solutions exist. From the computed spectra, we expect that, as $N$ increases, the gap in the real (imaginary) part should open(close) progressively closer to $N\Gamma/2$ ($\Omega_{c}$) and in a more abrupt way.

In Fig. 12 we confirm this hypothesis. In panel (a) we show the region around which the gap in the real part of the spectrum closes. As $N\rightarrow\infty$, the curve becomes sharper and moves closer to $N\Gamma/2$. In panel (b), we show the gap in the imaginary part of the spectrum for different values of $N$ and observe evidence that, effectively, the point at which it opens gets closer to $\Omega_{c}$ as $N\rightarrow\infty$ [see inset in panel (b)].

VI Conclusions

In this work, we have shown that the addition of interactions to a driven-dissipative collective system doesn’t change qualitatively the steady state properties but does affect significantly its dynamical behaviour in the limit of a large number of atoms. While in the noninteracting case both a slow relaxation and oscillations of the atomic inversion occur beyond a critical driving strength, in the presence of interactions a new metastable region appears. We showed that this emergent behavior could be understood both via a mean field analysis, in which case the metastable region manifests in the dynamics, and via an exact treatment, where for a finite system fingerprints of the metastable region can be observed in the spectrum of the Liouville operator. Finally, we also demonstrated that the slow relaxation present in one of the phases follows a scaling law different from the cases of pure interactions or pure dissipation, i.e. both must be present with comparably equal strengths to generate a non-trivial relaxation dynamics.

Given recent progress in loading arrays of strontium atoms in an optical cavity Norcia et al. (2016), and the capability to probe the associated electronic clock transition ${}^{1}S_{0}\rightarrow{}^{3}P_{0}$ we believe experiments are reaching the level of control to observe the predicted behavior. Nevertheless given the ultra-long relaxation time in some of the phases it is worthwhile to re-examine in the future our predictions including possible extra decoherence processes such as light scattering and atom loss, which are unavoidable in real experiments.

An interesting parallel direction for future work would be to use the ${}^{1}S_{0}\rightarrow{}^{3}P_{1}$ transitions, which is still quite narrow compared to other dipole allowed transition (lifetime of the order of $20$ $\mu$s) but at the same time six orders of magnitude faster that the clock transition. The faster dynamics relaxes the stringent requirements on atom loss and decoherence on one hand, but on the other might invalidate the adiabatic elimination of the cavity photons. This is due to the comparable decay rates of the atomic transition and cavity. The necessity of including the photons during the dynamics might nevertheless make the dynamical phases richer and the phase diagram more complex than the one studied here, thus opening interesting new avenues of research.

VII Acknowledgements

We thank Juan A. Muniz Silva, Zhe-Xuan Gong and Itamar Kimchi for feedback on the manuscript and Julia Cline for interesting discussions. This work is supported by the Air Force Office of Scientific Research grants FA9550-18-1-0319 and its Multidisciplinary University Research Initiative grant(MURI), by the Defense Advanced Research Projects Agency (DARPA) and Army Research Office grant W911NF-16-1-0576, the National Science Foundation grant PHY-1820885, JILA-NSF grant PFC-173400, and the National Institute of Standards and Technology.

Appendix A Envelope of oscillations in N phase

To understand the envelope in Fig. 7(F), we go into the interaction picture generated by the driving term, $\Omega\hat{J}_{x}$, and time average the resulting Liouville superoperator over one oscillation. The ensuing equation for the time averaged density matrix, $\hat{\rho}^{E}$, is:

[TABLE]

and, as shown in Fig. 13(a), describes correctly the envelope of $\hat{J}_{z}$. The situation for $\hat{J}_{x}$ and $\hat{J}_{y}$ is worse, however, because for finite $\Omega$ they have a nonzero steady state value. To describe them accurately, higher order terms in the averaging expansion are needed (which is effectively an expansion in $\frac{1}{\Omega}$). The fact that Eq. (26) is independent of $\Omega$ translates directly into the previously mentioned fact that the envelope is $\Omega$ independent.

It would be useful to understand how does the envelope scale with $N$, but even when working with Eq. (26), there is no closed form solution. The decay of $\hat{J}_{z}$ is given by:

[TABLE]

in the case of $\chi=0$ with $t_{decay}\sim\Gamma^{-1}$ and

[TABLE]

for $\Gamma=0$ with $t_{\mathrm{decay}}\sim\chi^{-1}N^{-\frac{1}{2}}$. To determine the scaling when both terms are present and $\chi\sim\Gamma$ (as is our case), we plotted the envelope for $N=200$, $400$, $800$, extracted the $e^{-\frac{3\Gamma t}{4}}$ time dependence (since it always dominates at short times) and then scaled the graphs to put them on top of each other, as shown in Fig. 13(b). The scaling exponent found was $t_{\mathrm{decay}}\sim N^{-0.66}$, which is corroborated by perturbation theory calculations, which indicate a leading dependence of the form of Eq. (29) to sixth order. The exponent is then predicted to be $\frac{2}{3}$, very much in accordance with the numerical result,

[TABLE]

A more satisfactory description of this behaviour can be obtained in the following way: we multiply Eq. (26) by $\tilde{J}^{+}=\hat{J}_{z}-i\hat{J}_{y}$ on the right and then do a projection to obtain the equation satisfied by the matrix element $\bra{m}\rho\tilde{J}^{+}\ket{m}\equiv f(m)$, where $\ket{m}$ are eigenstates of $\hat{J}_{x}$.

[TABLE]

Since $f(t=0,m)\approx\sqrt{\frac{N}{2\pi}}e^{-\frac{m^{2}}{N/2}}$ we introduce a variable $x=\frac{m}{\sqrt{N/2}}$, which becomes continuous in the limit $N\rightarrow\infty$, and define the small quantity $\Delta=\frac{1}{\sqrt{N/2}}$. In terms of these, Eq. (30) takes the more sugestive form of Eq. (31).

[TABLE]

and we are now considering $f$ as a function of $x$ too. The limit $N\rightarrow\infty$ corresponds to $\Delta\rightarrow 0$ upon which we obtain a partial differential equation (omitting terms of order $\Delta$),

[TABLE]

In principle, we would have to solve this equation with the initial condition $f(0,x)=\sqrt{\frac{N}{2\pi}}e^{-x^{2}}$, corresponding to a spin pointing in the $z$ direction, and then integrate with respect to $x$ to obtain the average value of $\braket{\hat{J}_{z}-i\hat{J}_{y}}$ at later times, since:

[TABLE]

This becomes easier to treat in Fourier space (with respect to $x$), in which case the relevant equations are:

[TABLE]

If we neglect the second order partial derivative (since $\tilde{f}(0,k)$ is slowly varying) a closed form expression can be written for the solution (we also neglect the $k^{4}$ term for simplicity),

[TABLE]

For $\chi t,\Gamma t\ll 1$, this becomes

[TABLE]

which agrees with the results from perturbation theory, up to to an irrelevant minus sign. We can compare both expressions, Eqs. (35) and (36), to the exact solution for $N=400$ and $N=10000$ (Fig. 14). Finally, we remark that this formulation of long time corrections to mean field using a single partial differential equation is always achievable after doing an averaging over a single particle constant drive, as in the case of Eq. (26). If we tried to do this with the original master equation, we would find, instead, an infinite set of coupled partial differential equations.