User:Kchinni/sandbox/Dicke Model

Introduction
The Dicke Hamiltonian describes $N$ two-level atoms, with energy splitting $\omega_{0}$, interacting with a single mode $\omega$  of the electromagnetic field. The Hamiltonian is given by$$H=\hbar\omega_{0} J_{z}+\hbar\omega a^{\dagger}a+\frac{\lambda}{\sqrt{N}}\left(a+a^{\dagger}\right)\left(J_{+}+J_{-}\right)$$where are the creation and the annihilation operators for the field satisfying $ \left[a^{\dagger},a\right]=1$ and $ J_{z}=\sum_{i} J_{z}^{(i)}$  ,$J_{\pm}=\sum_{i} J_{\pm}^{(i)}$  are the collective spin operators of the atoms satisfying the usual angular momentum commutation relations $\left[J_{z}, J_{\pm}\right]=\pm J_{\pm} $  and $\left[J_{+}, J_{-}\right]=2J_{z}$.

Since the original dipole coupling strength is proportional to $\frac{1}{\sqrt{V}}$, where $V$ is the volume of the cavity, one can factor out $\frac{1}{\sqrt{N}}$  from the original coupling constant by substituting $V=\frac{N}{\rho}$  where $\rho $ is the density of the atoms in the cavity. This is the reason for the presence of the factor $\frac{1}{\sqrt{N}}$ in the above Hamiltonian. This Dicke Hamiltonian commutes with $J^{2}$, and also commutes with the parity operator $e^{i\pi N} $ , where $N=a^{\dagger}a+J_{z}+j$ is the number of excitations in the system. It is straight-forward to check that $\left[H,\Pi\right]=0$ using the relations  $e^{i\pi a^{\dagger}a}\left(a+a^{\dagger}\right) e^{-i\pi a^{\dagger}a}	=-\left(a+a^{\dagger}\right)$  and $e^{i\pi J_{z}}J_{x}e^{-i\pi J_{z}}	=-J_{x}$. Notice that the $N $ atoms in this model are not interacting with each other. The atoms are also considered to be identical, so each atom in the ensemble has the same energy splitting $$\omega_{0}$$ and the same coupling constant $$\lambda$$ with the field. However, they are considered to be distinguishable, which means that each atom can be given a label.

Some of the interesting phenomena shown by this model include the divergence of atom-field entanglement at $\lambda_{c}=\frac{\sqrt{\omega\omega_{0}}}{2}$ and a quantum phase transition at $\lambda_{c}$  between the “normal phase” and the “super-radiant phase” in the thermodynamic limit. In addition, it has been shown this Hamiltonian in the rotating-wave approximation shows a thermal phase transition between the normal phase and the super-radiant phase. For $\lambda>2\lambda_{c}$, the system undergoes the thermal phase transition at the temperature, $T_{c}=\frac{\omega_{0}}{2\omega k_{B}T \ tanh^{-1}\left(\frac{\omega\omega_{0}}{\lambda^{2}}\right)}$. Below this temperature, the system is in super-radiant phase where the atomic ensemble emits light with the intensity proportional to $N^{2}$, which is stronger emission compared to the case where all atoms emit independently. Above this temperature $T_{c}$, the system is in normal phase. For $$\lambda<2\lambda_{c}$$, the system doesn't exhibit a thermal phase transition. Finally, this Hamiltonian also shows a change in the nearest-neighbor level spacing distribution from the Poissonian distribution to the Wigner-Dyson distribution as $\lambda$ goes through $\lambda_{c}$, indicative of the presence of chaos in the system for $\lambda>\lambda_{c}$.

Rotating-wave approximation
In the quantum-optical systems, the frequencies $$\omega_{0}$$ and $$\omega $$ are, in general, very large compared to the coupling strength $\lambda $ . Therefore, the terms that oscillate at the frequency $\omega_{0}+\omega$ average to zero on the time scale set by $$\lambda$$. For example, the system realized in has $\omega_{0}+\omega\approx2\pi\times10^{14}Hz$  compared to $\lambda\approx2\pi\times10^{7}Hz$. Therefore, one could neglect the rapidly oscillating terms(rotating-wave approximation) in the Dicke Hamiltonian to obtain the following Hamiltonian

$$H=\hbar\omega_{0}J_{z}+\hbar\omega a^{\dagger}a+\frac{\lambda}{\sqrt{N}}\left(a^{\dagger}J_{-}+aJ_{+}\right)$$

Also this Hamiltonian has an additional symmetry. The number of excitations $N=a^{\dagger}a+J_{z}+j$ are conserved in this case.

Quantum Phase Transition
The Dicke Hamiltonian shows a $$2^{nd}$$ order QPT at $\lambda_{c}=\frac{\sqrt{\omega_{0}\omega}}{2}$ in the thermodynamic limit, $N\rightarrow\infty$. The presence of this QPT can be seen by solving the Dicke Hamiltonian in the thermodynamic limit, and looking at the second derivative of the ground state energy. Since $\left[H\left(\lambda\right),J^{2}\right]=0$, one knows that the ground state of the system is always in the $j=\frac{N}{2}$ subspace. This can understood based on the following thought experiment. If $\lambda=0$, then the ground state of the system is $|n=0\rangle\otimes|\frac{N}{2},-\frac{N}{2}\rangle$. Suppose the value of $\lambda$ is now adiabatically increased with |$\psi(0)\rangle=|0\rangle\otimes|\frac{N}{2},-\frac{N}{2}\rangle$, then the state of the system at the later time $|\psi(t)\rangle$  will always be the instantaneous ground state of the Hamiltonian. We also know that $\left[H\left(\lambda\right),J^{2}\right]=0$. These two conditions imply that the ground state of the Hamiltonian will always remain in the $j=\frac{N}{2} $ subspace. Therefore, $j=\frac{N}{2}$ will be assumed for the rest of the discussion associated with the QPT. In this case, the system can be treated as a single large pseudo-spin system consisting of $N+1$ levels.

The spin operators in the Hamiltonian can be expressed in terms of the bosonic operators using the Holstein-Primakoff mapping:$$\begin{align} J_{+}&=b^{\dagger}\sqrt{2j-b^{\dagger}b} \\ J_{-}&={\sqrt {2j-b^{\dagger }b}}\ b \\ J_{z}&=b^{\dagger}b-j \\ \end{align}$$

Substituting these expressions into the Dicke Hamiltonian, the Hamiltonian becomes

This Hamiltonian can be now diagonalized in the thermodynamic limit. Two effective Hamiltonians, $$H^{(1)}$$ and $$H^{(2)}$$, will be derived: one for $$\lambda<\lambda_{c}$$ and another for$$\lambda>\lambda_{c}$$ with each one valid in their respective regimes.

Normal Phase
Taylor expanding the square root function in the Hamiltonian mentioned in 1, and then neglecting all the terms with any powers of $j$ in the denominator will produce a Hamiltonian that is bilinear in the bosonic operators. That bilinear Hamiltonian can then be transformed using the Bogoliubov transformation to the following Hamiltonian :$$H^{(1)}=\epsilon_{-}^{(1)}c_{1}^{\dagger}c_{1}+\epsilon_{-}^{(2)}c_{2}^{\dagger}c_{2}+\frac{1}{2}\left(\epsilon_{+}^{(1)}+\epsilon_{-}^{(1)}-\omega-\omega_{0}\right)-j\omega_{0}$$where

$c_{1}=\frac{1}{2}\left[\frac{cos\gamma^{(1)}}{\sqrt{\omega\epsilon_{-}^{(1)}}}\left[\left(\epsilon_{-}^{(1)}-\omega\right)a^{\dagger}+\left(\epsilon_{-}^{(1)}+\omega\right)a\right]-\frac{sin\gamma^{(1)}}{\sqrt{\omega_{0}\epsilon_{-}^{(1)}}}\left[\left(\epsilon_{-}^{(1)}-\omega_{0}\right)b^{\dagger}+\left(\epsilon_{-}^{(1)}+\omega_{0}\right)b\right]\right]$ $c_{2}=\frac{1}{2}\left[\frac{sin\gamma^{(1)}}{\sqrt{\omega\epsilon_{+}^{(1)}}}\left[\left(\epsilon_{+}^{(1)}-\omega \right)a^{\dagger}+\left(\epsilon_{+}^{(1)}+\omega \right)a\right]+\frac{cos\gamma^{(1)}}{\sqrt{\omega_{0}\epsilon_{+}^{(1)}}}\left[\left(\epsilon_{+}^{(1)}-\omega_{0}\right)b^{\dagger}+\left(\epsilon_{+}^{(1)}+\omega_{0}\right)b\right]\right]$ $\gamma^{(1)}	=\frac{1}{2}tan^{-1}\left(\frac{4\lambda\sqrt{\omega \omega_{0}}}{\omega_{0}^{2}-\omega^{2}}\right)$

and

$\left(\epsilon_{\pm}^{(1)}\right)^{2}=\frac{1}{2}\left(\omega_{0}^{2}+\omega^{2}\pm\sqrt{\left(\omega_{0}^{2}-\omega^{2}\right)^{2}+16\lambda^{2}\omega\omega_{0}}\right)$

From the above equation, one can notice that the excitation energy $$\epsilon_{-}^{(1)}$$is real only when $\lambda<\lambda_{c}$. Hence, this Hamiltonian is valid only when $\lambda<\lambda_{c}$. Also, the ground state energy in this phase is $$-j\omega_{0}$$.

Super-radiant Phase
Since the field and the atomic ensemble acquire macroscopic occupations in the super-radiant phase, the bosonic modes needs to be displaced in one of the following ways:

and where $$\alpha$$ and $$\beta$$ are assumed to be of the $$\mathcal{O}\left(j\right)$$. After substituting these displaced modes(choosing the positive displacement) in equation 1, the Hamiltonian becomes $$\begin{align} H^{(2)}_{+} & =\omega_{0}\left(d^{\dagger}d-\sqrt{\beta}\left(d^{\dagger}+d\right)+\beta-j\right)+\omega\left(c^{\dagger}c+\sqrt{\alpha}\left(c+c^{\dagger}\right)+\alpha\right) \\ & +\lambda\sqrt{\frac{k}{2j}}\left(c^{\dagger}+c+2\sqrt{\alpha}\right)\left(d^{\dagger}\sqrt{\xi}+\sqrt{\xi}d-2\sqrt{\beta}\sqrt{\xi}\right) \\ \end{align}

$$where $\sqrt{\xi}=\sqrt{1-\frac{d^{\dagger}d-\sqrt{\beta}\left(d^{\dagger}+d\right)}{k}}$ and $$k=2j-\beta$$. Now, Taylor expanding the $\sqrt{\xi}$ and neglecting all the terms with any powers of $j$  in the denominator will result in a Hamiltonian that contains terms both linear and bilinear in the bosonic modes. Then setting $\sqrt{\alpha}=\frac{2\lambda}{\omega}\sqrt{\frac{j}{2}\left(1-\left(\frac{\lambda_{c}}{\lambda}\right)^{4}\right)}$ and $$\sqrt{\beta}=\sqrt{j\left(1-\mu\right)}$$ will eliminate the linear terms, and the Hamiltonian can then be diagonalized using the Bogoliubov transformation. The final Hamiltonian is shown below.

$\begin{align} H_{+}^{(2)} & =\epsilon_{-}^{(2)}e_{1}^{\dagger}e_{1} +\epsilon_{+}^{(2)}e_{1}^{\dagger}e_{1}-j\left(\frac{2}{\omega}\left(\frac{\lambda_{c}}{\lambda}\right)^{4} + \frac{\omega\omega_{0}^{2}}{8\lambda^{2}}\right) \\ & +\frac{1}{2}\left(\epsilon_{+}^{(2)}+\epsilon_{-}^{(2)}-\frac{\omega_{0}}{2}\left(\frac{\lambda}{\lambda_{c}}\right)^{2}\left(1+\left(\frac{\lambda_{c}}{\lambda}\right)^{2}\right) -\omega-\frac{2\lambda^{2}}{\omega}\left(1-\left(\frac{\lambda_{c}}{\lambda}\right)^{2}\right)\right) \end{align}

$

where

$e_{1} =\frac{1}{2}\left[\frac{cos\gamma^{(2)}}{\sqrt{\omega\epsilon_{-}^{(2)}}}\left[\left(\epsilon_{-}^{(2)}-\omega\right)c^{\dagger}+\left(\epsilon_{-}^{(2)}+\omega\right)c\right]-\frac{sin\gamma^{(2)}}{\sqrt{\tilde{\omega}\epsilon_{-}^{(2)}}}\left[\left(\epsilon_{-}^{(2)}-\tilde{\omega}\right)d^{\dagger}+\left(\epsilon_{-}^{(2)}+\tilde{\omega}\right)d\right]\right]

$

$e_{2}=\frac{1}{2}\left[\frac{sin\gamma^{(2)}}{\sqrt{\omega_{c}\epsilon_{+}^{(2)}}}\left[\left(\epsilon_{+}^{(2)}-\omega_{c}\right)c^{\dagger}+\left(\epsilon_{+}^{(2)}+\omega_{c}\right)c\right]+\frac{cos\gamma^{(2)}}{\sqrt{\tilde{\omega}\epsilon_{+}^{(2)}}}\left[\left(\epsilon_{+}^{(2)}-\tilde{\omega}\right)d^{\dagger}+\left(\epsilon_{+}^{(2)}+\tilde{\omega}\right)d\right]\right]

$

$\gamma^{(2)}=\frac{1}{2}tan^{-1}\left(\frac{2\omega\omega_{0}}{\omega_{0}^{2}-\omega^{2}\left(\frac{\lambda_{c}}{\lambda}\right)^{4}}\left(\frac{\lambda_{c}}{\lambda}\right)^{4}\right)

$, $\tilde{\omega}=\frac{\omega_{0}}{2}\left(1+\left(\frac{\lambda_{c}}{\lambda}\right)^{4}\right)

$

and

$\left(\epsilon_{\pm}^{(2)}\right)^{2}=\frac{1}{2}\left(\omega_{0}^{2}\left(\frac{\lambda}{\lambda_{c}}\right)^{4}+\omega^{2}\pm\sqrt{\left(\omega_{0}^{2}\left(\frac{\lambda}{\lambda_{c}}\right)^{4}-\omega^{2}\right)^{2}+4\omega^{2}\omega_{0}^{2}}\right)

$

Notice that $$\epsilon_{-}^{(2)}

$$ is real only if $\lambda\geq\lambda_{c}

$, and hence the Hamiltonian is valid only for $\lambda\geq\lambda_{c}

$ . The ground state energy is phase is then given by $-j\left(\frac{2}{\omega}\left(\frac{\lambda_{c}}{\lambda}\right)^{4}+\frac{\omega\omega_{0}^{2}}{8\lambda^{2}}\right)

$ . If one had chosen the displacements with the negative signs in equation 2, then the spectrum of $H_{-}^{(2)}$ would have been identical to the one obtained above. This implies that the energy spectrum of the Hamiltonian above $$\lambda_{c}$$ becomes degenerate. Also, note that this Hamiltonian $H^{(2)}$ doesn't commute with $$\Pi $$, $\left[H^{(2)},\Pi\right]\neq0$, but $\left[H^{(2)},\Pi^{(2)}\right]=0$  where $\Pi^{(2)}\equiv e^{i\pi\left(c^{\dagger}c+d^{\dagger}d\right)}$. This means the parity symmetry $\Pi$ becomes spontaneously broken and two new symmetries appeared, $\Pi^{(2)}$ (one associated with the positive displacement of the mode and the other one associated with the negative displacement of the mode).

Properties of the QPT
From the above analysis, the ground state energy is given by $$E_{G} = \left\{ \begin{array}{ll} -j\omega_{0} & \quad \lambda < \lambda_{c} \\ -j\left(\frac{2}{\omega}\left(\frac{\lambda_{c}}{\lambda}\right)^{4}+\frac{\omega\omega_{0}^{2}}{8\lambda^{2}}\right) & \quad \lambda > \lambda_{c} \end{array} \right.$$and the second derivative of the ground state energy is given by $$\frac{1}{j} \frac{d^{2}}{d\lambda^{2}}E_{G} = \left\{ \begin{array}{ll} 0 & \quad \lambda < \lambda_{c} \\ -\left(\frac{40}{\omega}\left(\frac{\lambda_{c}^{4}}{\lambda^{6}}\right)+\frac{3 \omega\omega_{0}^{2}}{4\lambda^{4}}\right) & \quad \lambda > \lambda \end{array} \right.$$ It can be seen from the above expression and the plot that there is a discontinuity in the second-derivative of the ground state energy at $\lambda_{c}$. The plot of the excitation energies indicate that $$\epsilon_{-}$$ is the energy gap between the ground state and the first excited state, and this gap goes to zero at $$\lambda=\lambda_{c}$$ as expected for a second-order QPT. Specifically, this gap closes as $\Delta\sim\sqrt{\frac{32\lambda_{c}^{3}\omega_{0}^{2}}{\omega_{0}^{4}+16\lambda_{c}^{4}}}|\lambda-\lambda_{c}|^{\frac{1}{2}} $ , and the characteristic length scale diverges as $l_{-}=\frac{1}{\sqrt{\epsilon_{-}}}=\left(\frac{\omega_{0}^{4}+16\lambda_{c}^{4}}{32\lambda_{c}^{3}\omega_{0}^{2}}\right)^{\frac{1}{4}}|\lambda-\lambda_{c}|^{-\frac{1}{4}}   $ [3]. These expressions imply that the critical exponents are given by $\nu=\frac{1}{4} $  and $z=2$, as the characteristic length diverges as $|\lambda-\lambda_{c}|^{-\nu}  $  and the energy gap closes as $|\lambda-\lambda_{c}|^{z\nu}  $.

Chaos
We know that the Dicke Hamiltonian is integrable in the thermodynamic limit, $$j\rightarrow\infty$$, as it was explicitly diagonalized above in the quantum phase transition section. The Dicke Hamiltonian is also integrable in the $$\lambda\rightarrow\infty$$ limit because in that case, the Hamiltonian can be written as follows:$$H=\omega a^{\dagger}a+2\frac{\lambda}{\sqrt{N}}\left(a+a^{\dagger}\right)J_{x}$$and, this Hamiltonian can be diagonalized by simply displacing the bosonic modes. The eigenvalues in this case are given by $E_{km}=\frac{\omega}{j}k-\frac{2\lambda^{2}}{\omega j^{2}}m^{2}$  where $$k =0,1,...$$ and $$m$$ can be $$\{j,j-1,...-j\}$$. For the general case, the Dicke Hamiltonian can be numerically diagonalized in the basis $\left\{ |n\rangle\otimes|j,m\rangle\right\}$. To determine the chaoticity in this general case, the authors in have looked at the nearest-neighbor level spacing distribution, $P\left(S\right)$, where $S_{n}=E_{n+1}-E_{n}$. Then, the chaoticity of the system was quantified using the quantity $$\eta$$, which is defined as $\eta\equiv\left|\frac{\int_{0}^{S_{0}}dS\left[P(s)-P_{W}(S)\right]}{\int_{0}^{S_{0}}dS\left[P_{P}(s)-P_{W}(S)\right]}\right|$ where $P_{P}(S)$  is the Poissonian distribution, $$P_{w}(S)$$ is the Wigner-Dyson distribution, and $S_{0}=0.472913$, which is the value at which $P_{P}(S) $  and $P_{W}(S)$  first intersect. Notice that $\eta=0$ if $P(S)=P_{W}(S)$  and $\eta=1$  if $P(S)=P_{P}(S)$, so the system becomes more chaotic as $$\eta$$ goes towards 0. If one plots the quantity $$\eta$$ as a function of $$\lambda$$ for the case of the Dicke Hamiltonian, it can be seen from the plot that the value of $$\eta$$ will be decreasing as $$\lambda$$ is increased $$0$$ to $$\lambda_{c}$$, and above $$\lambda_{c}$$, it is almost zero. In summary, the Hamiltonian seems to become chaotic as $$\lambda$$ passes through $$\lambda_{c}$$. This behavior of the chaoticity seems to agree well with the chaoticity of the semi-classical Hamiltonian.

Semi-classical Analysis
As there is no direct analog for the spin operators in the classical mechanics, there are many different approaches to perform semiclassical analysis. Here two different approaches are presented.

Method 1: In this method, the differential equations associated with the expectation values of the operators are obtained, and these expectation values are replaced by the classical variables to obtain a set of nonlinear equations of motion for these observables. More specifically, the time evolution of the expectation values of the observables is given by the equation $\frac{d}{dt}\langle A\rangle=i\langle[H,A]\rangle+\langle\frac{\partial A}{\partial t}\rangle$. Using this equation, one can obtain$$\begin{align} \frac{d\langle a\rangle}{dt}	&=-i\omega\langle a\rangle-i\frac{\lambda}{\sqrt{2j}}\langle J_{x}\rangle \\

\frac{d\langle J_{-}\rangle}{dt}	&=-i\omega_{0}\langle J_{-}\rangle+\frac{2i\lambda}{\sqrt{2j}}\langle\left(a+a^{\dagger}\right)J_{z}\rangle \\

\frac{d\langle J_{z}\rangle}{dt}	&=i\frac{\lambda}{\sqrt{2j}}\langle\left(a+a^{\dagger}\right)\left(J_{-}-J_{+}\right)\rangle \\

\end{align}$$

Assuming that the expectation values of the product of the observables is the product of expectation values of the observables, we have $$\begin{align} \frac{d\alpha}{dt}&=-i\omega\alpha-i\frac{\lambda}{\sqrt{2j}}\left(\beta+\beta^{*}\right)\\ \frac{d\beta}{dt}	&=-i\omega_{0}\beta+\frac{2i\lambda}{\sqrt{2j}}\left(\alpha+\alpha^{*}\right)w \\

\frac{dw}{dt}	&=i\frac{\lambda}{\sqrt{2j}}\left(\alpha+\alpha^{*}\right)\left(\beta-\beta^{*}\right)\\

\end{align}$$where $\alpha\equiv\langle a\rangle$, $\beta\equiv\langle J_{-}\rangle$ and $w\equiv\langle J_{z}\rangle$. The above set of equations also conserve the pseudo angular momentum,$w^{2}+|\beta|^{2}=\langle J_{x}\rangle^{2}+\langle J_{y}\rangle^{2}+\langle J_{z}\rangle^{2}=\frac{N^{2}}{4}$. For $\lambda<\frac{\sqrt{\omega_{0}\omega}}{2}$, there are two fixed points given by$\alpha=\beta=0$ and $w=\pm\frac{N}{2}$  with the positive population inverse being unstable and the negative population inverse being stable, which could have been expected by looking at the ground state of the Dicke Hamiltonian at  $$\lambda=0$$. For $\lambda>\frac{\sqrt{\omega_{0}\omega}}{2}$, the fixed points $\alpha=\beta=0$ and $w=\pm\frac{N}{2}$  become unstable and two new stable fixed points appear $\alpha	=\pm\sqrt{N}\frac{\lambda}{\omega}\sqrt{1-\left(\frac{\lambda_{c}}{\lambda}\right)^{4}}$ , $\beta	=\mp\frac{N}{2}\sqrt{1-\left(\frac{\lambda_{c}}{\lambda}\right)^{4}}$  and $w	=-\frac{N}{2}\left(\frac{\lambda_{c}}{\lambda}\right)^{2}$. These stable points could have also been expected as well because the Hamiltonian in the $$\lambda\longrightarrow\infty $$ has a degenerate spectrum with the eigenvectors that are superpositions of $$|0\rangle\otimes|m_{x}\rangle $$ and $$ |0\rangle\otimes|-m_{x}\rangle $$.

Method 2: Another method to perform semiclassical analysis is to obtain the semiclassical Hamiltonian by writing the bosonic modes in equation 1 in terms of the corresponding position and the momentum operators, and then setting the commutators between the position and the corresponding momentum operators to zero. This will result in the Hamiltonian given by
 * $$H=\frac{1}{2}\left(p_{x}^{2}+p_{y}^{2}\right)+\frac{1}{2}\left(\omega^{2}x^{2}+\omega_{0}^{2}y^{2}\right)+2 \lambda \sqrt{\omega \omega_{0}} x y \sqrt{1-\frac{\omega_{0}^{2}y^{2}+p_{y}^{2}-\omega_{0}}{4j \omega_{0}}} -\frac{1}{2} \left(\omega + \omega_{0}+2j \omega_{0} \right) $$$$

where $a\equiv\sqrt{\frac{\omega}{2}}\left(x+\frac{i}{\omega_{0}}p_{x}\right) $ and $b\equiv\sqrt{\frac{\omega_{0}}{2}}\left(y+\frac{i}{\omega_{0}}p_{y}\right)$. The Hamilton's equations of motion associated with the above Hamiltonian are then given by
 * $$\begin{align}

\dot{x}	&=p_{x} \\ \dot{y}	&=p_{y} \left(1-\frac{\lambda}{2j} \sqrt{\frac{\omega_{c}}{\omega_{0}}} \frac{xy}{\sqrt{1-\eta}}\right)\\ \dot{p_{x}}	&=-\omega^{2}x-2\lambda\sqrt{\omega_{0}\omega}y\sqrt{1-\eta}\\

\dot{p_{y}}&	=-\omega_{0}^{2}y-2\lambda\sqrt{\omega_{0}\omega}x\sqrt{1-\eta}\left(1-\frac{\omega_{0}y^{2}}{4j\left(1-\eta\right)}\right)\\ \end{align}$$$$ where $\eta=\frac{1}{4j\omega_{0}}\left(\omega_{0}^{2}y^{2}+p_{y}^{2}-\omega_{0}\right) $ . One of the fixed point is given by $$x=y=p_{x}=p_{y}=0 $$, which is equivalent to the fixed point $\alpha=0$ and $w=-\frac{N}{2}$  obtained by the first method. However, this fixed point is stable only for $\lambda<\frac{\lambda_{c}}{\sqrt{1+\frac{1}{4j}}} $ unlike the condition $$\lambda<\lambda_{c}$$ in method 1. The other two fixed points $p_{x}=p_{y}=0 $ and $x_{0}=\pm\frac{2\lambda}{\omega}\sqrt{\frac{j}{\omega}\left[\left(1+\frac{1}{4j}\right)^{2}-\left(\frac{\lambda_{c}}{\lambda}\right)^{4}\right]} $

and $y_{0}=\mp\sqrt{\frac{2j}{\omega}\left[\left(1+\frac{1}{4j}\right)-\left(\frac{\lambda_{c}}{\lambda}\right)^{2}\right]} $ . These two fixed point exist and are stable only for

$\lambda>\frac{\lambda_{c}}{\sqrt{1+\frac{1}{4j}}} $.

From the Hamiltonian in equation 3, it can be noticed that this semiclassical system can be treated as a particle moving in the following potential well: $U=\frac{1}{2}\left(\omega^{2}x^{2}+\omega_{0}^{2}y^{2}\right)+2\lambda\sqrt{\omega\omega_{0}}xy\sqrt{1-\frac{\omega_{0}^{2}y^{2}+p_{y}^{2}-\omega_{0}}{4j\omega_{0}}} $ .The plots of this potential well are shown for $$p_{y}=0 $$. As mentioned in, these contour plots don't show bifurcation for all the values of $$p_{y} $$. Using the solutions obtained from the Hamilton's equations in equation 4, one can plot the Poincare' sections by setting $$p_{x}$$ to a certain value. Note that the value of $$p_{y}$$ is fixed by the energy. The Poincare' sections have been plotted in. In these figures, one can notice that as the value of $$\lambda$$ is increased from 0 to 1, the trajectories in the phase space appear to become chaotic. Note that this increase in chaoticity as $$\lambda$$ crosses through $$\lambda_{c}$$ is also expected from the results of the random matrix theory for the corresponding quantum Hamiltonian.

Experimental Implementation
The Dicke model could be implemented using multiple methods. To realize the Dicke model in the absence of a cavity, one would need a large sample of atoms, and the distance between the atoms should be smaller than the wavelength of the field for the collective spontaneous emission. Since all the atoms in the Dicke model should experience the same amount of electric field, large number of atoms has to be arranged in a regular pattern within a small distance for the realization of the model. On the other hand, the atoms should not be too close else there will be dipole-dipole interactions between the atoms. However, in the presence of the cavity, high cavity finesse allows one to use small sample of atoms to achieve superradiance.

Ring Cavity and Neutral Atoms
This experimental scheme is discussed in. In this implementation scheme, an ensemble of atoms are confined in the ring cavity. Two laser fields with Rabi frequencies $$\Omega_{r}$$ and $$\Omega_{s}$$ enter through one cavity mirror, and the quantized field enter and exit through one other cavity mirror as shown in the figure. The two level atom is realized through a pair of Raman channels, whose atomic excitation scheme is shown below.

In, it is shown that this system evolves under the effective Dicke Hamiltonian:$$H=\omega a^{\dagger}a+\omega_{0}J_{z}+\frac{\lambda}{\sqrt{N}}\left(a+a^{\dagger}\right)\left(J_{+}+J_{-}\right)$$

where $$\omega=\frac{Ng_{r}^{2}}{\Delta_{r}}+\delta_{cav}$$, $$\omega_{0}	=\omega_{1}-\omega_{1}^{'}$$, $$\lambda	=\frac{\sqrt{N}g_{r}\Omega_{r}}{2\Delta_{r}}$$, $$\omega_{1}^{'}=\frac{\omega_{ls}-\omega_{lr}}{2}$$ and $$\delta_{cav}=\omega-(\omega_{ls}-\omega_{1}^{'})$$. Here $$\omega_{lr}$$ and $$\omega_{ls}$$ are the laser frequencies, and $$\omega_{r}$$,$$\omega_{s}$$ and $$\omega_{1}$$ are the atomic frequencies. Notice that this model also takes cavity loss into account, but the parameters $$\omega$$, $$\omega_{0}$$ and $$\lambda$$ can be chosen appropriately for the Hamiltonian dynamics to dominate.

Circuit QED
Circuit QED has been used to realize the Dicke model in the RWA for $N=3$. In this system, superconducting qubits act as atoms and the cavity mode will be the resonant waveguide mode. With CQED, one can precisely place the effective atoms and cavity frequency is also very controllable. For large $$N$$, the physical size of the device could be a problem.