Dynamics of a Generalized Dicke Model for Spin-1 Atoms

In certain parts of this work we consider a regime where the cavity-assisted Raman transitions used to implement the model are themselves off-resonant, i.e., when $\omega\gg\{\lambda_{\pm},\omega_{0}\}$. In this “dispersive” limit, the cavity mode is sparsely populated and effectively only mediates atom-atom interactions. We may then adiabatically eliminate the cavity mode by tracing over the cavity degrees of freedom and making standard approximations to arrive at the modified master equation,

where $\hat{\rho}_{A}$ is the atomic density operator, $\hat{H}_{S}$ is the adiabatically eliminated Hamiltonian,

and we have define

The full derivation of this adiabatically eliminated Hamiltonian is detailed in the supplemental material of Ref. [10].

We may also represent the atomic operators in terms of bosonic mode annihilation and creation operators (in the Jordan-Schwinger representation), denoted by $\hat{b}_{\pm}$ for the $m=\pm 1$ sublevels, and $\hat{b}_{0}$ for the $m=0$ sublevel. These obey the regular annihilation and creation operator algebra, namely

where $i,j=0,\pm$, and $\delta_{ij}$ is the Kronecker delta. In this representation, summation over an arbitrary atomic operator $\hat{X}$ can be mapped according to [20]

where $i$ and $j$ have the same meaning as above, and $\ket{i}\equiv\ket{m=i}$. For example, $\hat{H}_{QZ}$ is mapped to

One similarly finds that

When treating the system quantum-mechanically, we additionally specialize to the case where all atoms are initially prepared in the $m=0$ sublevel. If we further make the approximation that the ensemble contains a large number of atoms, i.e., essentially that $N\to\infty$, then we can assume that the $m=0$ sublevel remains macroscopically occupied throughout the evolution of the system. This allows us to replace the bosonic mode operator for the $m=0$ sublevel with a constant $c$-number, $\hat{b}_{0}\rightarrow\sqrt{N}$, and consequently to write $\hat{S}_{\pm}$ as

By further defining the operators

which satisfy $[\hat{A},\hat{A}^{\dagger}]=[\hat{B},\hat{B}^{\dagger}]=[\hat{A},\hat{B}]=0$, $[\hat{A},\hat{B}^{\dagger}]=[\hat{B},\hat{A}^{\dagger}]=-2$, we may rewrite the Hamiltonian as

where we have defined the constants

Note that, implicit in the undepleted mode approximation is the assumption that the populations of the $m=\pm 1$ states satisfy $\langle\hat{N}_{\pm}\rangle\ll N$ at all times. However, in certain regimes of parameter space, this condition is violated, which, as we shall see, manifests itself as divergence in the solutions of this simplified model.

Unlike the single-mode spinor BEC systems, our engineered Dicke model has additional freedom with regards to tuning parameters and interactions. For example, if we set $\lambda_{+}=\lambda_{-}=\lambda$ then atom-atom interactions are governed by a term solely proportional to $\hat{S}_{x}^{2}$ (rather than $\hat{S}^{2}$), which, as we shall see, leads to significant changes in behavior. In this subsection we still consider a closed system, where operators evolve subject to the Hamiltonian (II.2), to maintain the single-mode spinor BEC analogy. We find the equations of motion for $\hat{A}$, $\hat{B}$, and their adjoints in the fully general case are given by

The eigenvalues of the above matrix are found to be

where each choice of $\pm$ denotes a different eigenvalue. Two eigenvalues are negatives of the other two, and, consequently, if any eigenvalue has a nonzero real part there must be an eigenvalue with a positive real part, again leading to divergence. Hence, we need only consider two linearly independent eigenvalues to characterize the behavior of the system. We choose these to have the positive outer sign and opposite inner signs in (27), and we denote these by $\alpha_{\pm}$ according to said inner sign.

Fig. 3 shows the real and imaginary parts of $\alpha_{\pm}$ as functions of $q$ and $\omega_{0}$, for balanced coupling (i.e., $\lambda_{+}=\lambda_{-}$). We see distinct regions where $\alpha_{-}$ has a nonzero real part, which also include the regions where $\alpha_{+}$ has a nonzero real part. The system’s divergence is then conditional purely on the real part of $\alpha_{-}$, a condition which holds for all $\Delta_{-}$.

We can additionally find the boundaries separating the regions of divergence and oscillation, which occur when $\real e(\alpha_{-})=0$. This condition is satisfied in the specific case where $\alpha_{-}=0$, when the determinant of the matrix in Eq. (26) is necessarily zero. Although the converse is not generally true, i.e., a zero determinant does not imply $\alpha_{-}=0$, we may nevertheless find the determinant’s roots and choose those that correspond to the desired boundaries. The determinant is

and its four roots are

The remaining boundaries are given by the roots of the nested square root in (27), when $\alpha_{\pm}$ both become purely imaginary, and are found to be

These are all plotted over the eigenvalue landscapes in Fig. 3, where we see they indeed define all of the boundaries.

The eigenvalue landscape, or “phase diagram” for this generalized model clearly has much more structure than that of the effective BEC model of the previous section. There is now a strong dependence on both $\omega_{0}$ and $q$, with multiple regions of stability and instability as a function of either parameter. The behavior of the populations of the $m=\pm 1$ sublevels also display quite distinct behaviors in the different regions, as demonstrated in Fig. 4, where the populations are plotted as a function of time for each of the seven distinct regions of stability or instability. These are associated with the “slice” of the eigenvalue landscape from Fig. 3 through the line of constant $\omega_{0}=5\Delta_{+}$, which crosses through each of the seven aforementioned regions. In each case from Fig. 4, the population time series were obtained by numerically integrating the Heisenberg equations of motion for the relevant second-order operator expectations, using a standard 4th-order Runge-Kutta method with a sufficiently small stepsize. These equations of motion are given in Appendix A and are discussed further in the case of the general, open system in Sec. IV.2. In principle, however, this is not required; given the system in Eq. (26) is linear, one could solve it exactly. The solutions would be superpositions of complex exponentials, whose oscillation frequencies are given by the imaginary parts of the eigenvalues in Eq. (27). When none of those eigenvalues have a positive real part, as in Figs. 4(a), (c), (e), and (g), the solutions oscillate ad infinitum. Otherwise, as in Figs. 4(b), (d), and (f), the solutions diverge. Hence, although the behaviours in Fig. 4 may appear complicated, they are simple combinations of finitely many oscillation frequencies or exponential increases.

Figure  5 shows the real parts of $\alpha_{-}$ and the divergence boundaries for various choices of $\Delta_{-}$. We see a gradual transition from the landscape in Fig. 3 to a single, vertical band as $\Delta_{-}$ tends to $\Delta_{+}$ (i.e., as $\lambda_{+}\rightarrow 0$), as seen in the previous section. Most importantly, as evident in panel (b), the qualitative structure of the eigenvalue landscape survives through this transition.

Given our system is considered in a quantum optical setting, in order to realistically model it we must take dissipation into account. Doing so will give us a more complete description of the dynamics and thereby a possible explanation for the behavior of the system in the regions of divergence from the previous section.

We begin by considering the master equation in Eq. (7) to find the equations of motion for various operator moments. For an arbitrary operator $\hat{X}$, we have

Given $\hat{H}_{S}$ can be written in terms of just $\hat{S}_{\theta}$ and $\hat{N}_{\pm}$, in order to find the equation of motion for $\langle\hat{X}\rangle$ one simply needs to find the commutators $[\hat{X},\hat{S}_{\theta}]$, $[\hat{X},\hat{S}_{\theta}^{\dagger}]$, and $[\hat{X},\hat{N}_{\pm}]$. Doing so for $\langle\hat{A}\rangle$, $\langle\hat{B}\rangle$, and their conjugates, one finds

where we define the constants

We may now perform a similar eigenvalue analysis as in the previous section to find the regions in parameter space where this system diverges (or not). One should note that this is a system of moment equations, not operator equations, meaning the stability of (32) does not directly determine the stability of higher-order moments. Nevertheless, since the higher-order moments inherently involve the same operators as the first-order moments, we expect them to diverge or oscillate in the same regions of parameter space.

We should also note that although the matrix in (32) is not significantly different from that in (26), its eigenvalues have significantly more complicated forms. We therefore evaluate them numerically at each point in parameter space and extract the maximal real part from them, which is sufficient to determine the stability.

The determinant of the matrix is however comparatively simple, and given by

where we define

The roots of this determinant are

and

which, given the analysis in the previous section, we expect will correspond to boundaries of regions of divergence. In Fig. 6 we plot these roots alongside our numerical eigenvalue landscapes for various values of $\kappa$, where we see this is actually not the case.

The addition of a nonzero $\kappa$ adds a finite background to the landscape, which increases in magnitude with $\kappa$ and eventually dominates over the landscape. The boundaries between oscillation and divergence consequently become blurred, and the roots in (IV.2) just appear to separate regions of fast divergence (i.e., large real part) and slow divergence (i.e., small real part). In fact, the roots do not separate these regions perfectly, especially for small $q$ and $\omega_{0}$, but instead only give their approximate outline.

Lastly, we see that the boundaries disappear when $\kappa$ becomes sufficiently large. This occurs when $K$, and consequently the roots, acquire a nonzero imaginary part, or, more specifically, when $\kappa$ reaches the value

By relaxing the undepleted mode approximation, we inherently need to work with the full form of the collective spin operators, given by Eq. (14). Given products of these operators appear in the Hamiltonian, inclusion of the $m=0$ mode leads to a non-quadratic Hamiltonian, which leads to a coupling between increasingly higher-order moments in the equations of motion, which never form a closed set of equations. For example, first order moment equations couple to third order moments, which themselves couple to fifth order moments, etc.

We therefore take the thermodynamic limit, where $N\to\infty$ and the non-commuting operators can be replaced with commuting $c$-numbers. This allows us to factorize high order moments into products of first order moments, and thus arrive at a closed set of semiclassical nonlinear equations for the first order moments. We define the $c$-numbers

and use the Hamiltonian and master equation in Eqs.  (4) and (7) to arrive at the (closed) set of equations

where the equation of motion for $\beta_{-}$ can be found by applying the transformation $\mathcal{T}$ to Eq. (41).

However, this system of equations does not conserve the number of atoms, given by $N=|\beta_{+}|^{2}+|\beta_{-}|^{2}+|\beta_{0}|^{2}$, when $\kappa$ is nonzero. This is an artefact of the adiabatic elimination of the cavity mode, where cavity dissipation is essentially incorporated into the atomic dynamics, and thereby expressed as a changing total population. Since we already assumed and relied on number conservation, the system cannot be relied upon to give physical predictions. To consistently and properly include cavity dissipation we must include the cavity dynamics itself, and we do this in the following subsection. For the moment, however, we set $\kappa=0$ and the above equations reduce to

This system of equations can be readily integrated, and some example results are shown in Fig. 8 for several parameter choices. We see various types of population oscillations, ranging from roughly sinusoidal, to sharply peaked, to seemingly irregular. They can, however, be separated into two broad categories: small or large amplitude oscillations. The former involves sparse population of the $m=\pm 1$ sublevels and essentially sinusoidal oscillation, while the latter involves much larger population of those sublevels, accompanied by large depletion of the $m=0$ sublevel. These large-amplitude oscillations also tend to have the same general shape – an initial, roughly exponential increase, followed by a sharp peak, and exponential decrease.

Furthermore, the locations of these oscillations agree well with the eigenvalue landscapes of the earlier sections. The small-amplitude oscillations lie in regions of oscillation in the landscapes, while the large-amplitude oscillations lie in regions of divergence. Although here we can only confirm this result through numerical integration of the equations of motion (i.e., not through detailed analysis of its bifurcations), it is compelling and offers an explanation for the behavior of the system in regions of divergence.

Despite the insight we gained through the semiclassical description above, the issue of incorporating dissipation remains. To do so we must reintroduce the cavity mode and account for dissipation through cavity loss. If we once again take the thermodynamic limit, define the additional $c$-number

and use the full (i.e., not adiabatically eliminated) Hamiltonian and master equation, given by Eqs. 4) and (5) respectively, we arrive the semiclassical equations of motion

In addition to conserving the number of atoms, this system allows us to study the cavity dynamics. For example, by setting $\kappa=0$ and taking $\omega$ to be large, we can approach the dispersive limit and emulate our previous semiclassical non-dissipative results. In doing so we not only confirm those results, but gain additional insight into the previous sections’ regions of divergence. Figure 9 shows the cavity and atomic populations given by numerical integration of Eqs. (46)-(48) under a parameter regime which approximates the dispersive limit. The atomic populations once again follow either small or large amplitude oscillation, in agreement with the eigenvalue landscapes. Furthermore, as one might expect, at every significant depletion of the $m=0$ sublevel during large-amplitude oscillations a burst of photons is generated in the cavity.

Moreover, we now have the freedom to explore dissipative scenarios. Figs. 10(a-d) show the system evolution for four different parameter sets that correspond to the known phases of the Dicke model. In the normal phase, atoms may undergo a one-way transition from the $m=0$ sublevel to the other sublevels, accompanied by a burst of photons and subsequent decay towards a steady state with no photons in the cavity (Fig. 10(a)). In the superradiant phase, both the cavity and atomic modes settle into a steady state with constant populations, corresponding to balanced atomic transitions (Fig. 10(b)). In oscillatory superradiance the steady state is akin to the superradiant phase, but involves oscillating atomic and cavity populations (Fig. 10(c)). Lastly, we see signatures of normal-superradiant bistability, where the system appears to pass close to the normal phase equilibrium, but eventually reaches the superradiant equilibrium (Fig. 10(d)).

These are all in agreement with previous analyses of the two-level-atom version of the Dicke model [11], which not only mapped these phases, but also regions of chaotic behaviour. It therefore comes as no surprise that we also observe chaos, which specifically occurs when the counter-rotating terms in the Hamiltonian dominate, i.e., when $\lambda_{+}>\lambda_{-}$. Figs. 10(e-g) show variations of the chaos we observe: two are very similar to the two-level system chaos, where we see sudden spikes and slower decay of the cavity population. However, the other, with $q\neq 0$, is quite different; such spikes and decays are not as easily visible, and the cavity population appears to fluctuate randomly.

These chaotic behaviors are also manifest in the atomic dynamics. Although population time-series give a useful context to some behaviors, they become less useful when the dynamics are chaotic, given they only provide a one-dimensional projection of our eight-dimensional dynamical system. A more useful representation of the atomic dynamics is a plot of the system trajectory in the space spanned by the spin expectations (hence dubbed the “spin-space”) $s_{x}=\langle\hat{S}_{x}\rangle$, $s_{y}=\langle\hat{S}_{y}\rangle$, and $s_{z}=\langle\hat{S}_{z}\rangle$. Since such a representation is a three dimensional projection of the full system, its trajectories are allowed to self-intersect (though they do not self-intersect in the full eight-dimensional phase space), but they nevertheless provide information on attractors and equilibria in the system, and have indeed been employed in investigations of the two-level system [11].

In Fig. 11, we plot the spin-space trajectories of the chaotic dynamics shown in Fig. 10, which characterize three forms of chaos in our system. First, in Fig. 11(a), we see a trajectory typical of the two-level system: a symmetric (about reflections along the $s_{x}=s_{y}$ diagonal) chaotic attractor is present and the system maintains a constant spin-length (i.e., the trajectory lies on the Bloch sphere). In Fig. 11(b), we see an asymmetric chaotic attractor that also follows the Bloch sphere, which most likely arises due to the additional atomic degree of freedom compared with the two-level model. Lastly, in Fig. 11(c) the attractor is not confined to the Bloch sphere.

This may be explained by considering the total ensemble spin, which is only conserved when atomic transitions to the $m=+1$ magnetic sublevel result in energy changes equal in magnitude but opposite in sign to transitions to the $m=-1$ sublevel. In other words, the spin-length is only conserved when $q=0$. One can also see this by noting that $\hat{S}^{2}=\hat{S}^{2}_{z}+\hat{S}_{z}+\hat{S}_{-}\hat{S}_{+}$ commutes with all of the Hamiltonian terms except the term pertaining to the quadratic Zeeman shift. More specifically, the operator $\hat{b}_{+}^{\dagger}\hat{b}_{+}+\hat{b}_{-}^{\dagger}\hat{b}_{-}$ commutes with $\hat{S}_{z}=\hat{b}_{+}^{\dagger}\hat{b}_{+}-\hat{b}_{-}^{\dagger}\hat{b}_{-}$, but does not commute with $\hat{S}_{-}\hat{S}_{+}\propto(\hat{b}_{0}^{\dagger}\hat{b}_{+}+\hat{b}_{-}^{\dagger}\hat{b}_{0})(\hat{b}_{+}^{\dagger}\hat{b}_{0}+\hat{b}_{0}^{\dagger}\hat{b}_{-})$. Therefore, the two trajectories in Figs. 11(a,b) follow the Bloch sphere, while the trajectory in Fig. 11(c) does not.

In order to construct a preliminary map of the various phases and chaotic behaviour in parameter space, we use a simple numerical scheme to analyse the mean and any nonzero, steady-state oscillation amplitude. At each point in a grid of $\lambda_{+}$ and $\lambda_{-}$ values, we numerically integrate the system of equations, extract the cavity population over a late portion of the trajectory, and calculate the mean and oscillation amplitude (if any) of the sample. The mean is calculated by simply averaging over the sample, and the amplitude calculated by taking the difference between the maximum and minimum values and dividing it by two. It should be noted that the mean and oscillattion amplitude of the steady-state cavity population are sufficient to characterize the phase of the system: in the normal phase both the mean and amplitude will be zero, in the superradiant phase the mean will be nonzero but the oscillation amplitude is zero, and in the oscillatory phase both will be nonzero. Since no steady-state is reached in chaotic behaviour, we expect it to display seemingly random fluctuations in both the mean and amplitude.

We must note, however, that these characterizations only hold true under idealized conditions, and hinge on the system reaching the steady-state (when one exists). Therefore, if a particular trajectory takes too long to reach steady-state, or if the portion of the trajectory over which we perform our average is too short, the calculated mean and amplitude may not faithfully represent the true phase. In addition, since the calculation evolves the system in discrete timesteps, if the timestep is exactly a half-integer multiple of an oscillation period, the amplitude may register incorrectly as zero. Lastly, the calculation is somewhat reliant on the initial conditions, so the emergence of certain phases in a regime of multistability may not be detected. Nevertheless, these caveats will not influence the entirety of our calculations, so we can rely on them to give us at least a rough map for the location of each phase.

Fig. 12 shows the results of three such calculations, for choices of $q$ and $\omega_{0}$ pertaining to different energy-level structures. Regardless of the specific values of $q$ and $\omega_{0}$, we see regions where each phase is displayed, as per the criteria above, with the most notable separation between the normal phase and superradiance. In Figs. 12(b1,b2) and 12(c1,c2) there is also a seemingly sharp boundary between chaos and the other phases; in plots 12(b1,2) specifically, where $q=0$, the boundary appears to agree well with the equivalent phase diagram for the two-level model [11]. Lastly, we note that plots 12(a1,2) show no chaotic regions, possibly owing to the symmetrical energy level structure (i.e., $\omega_{0}=0$) being considered.