Characterization of the quantum phase transition in a two-mode Dicke model for different cooperation numbers

TL;DR

This paper investigates how the cooperation number influences quantum phase transitions in a multi-mode Dicke model, using variational states to analyze ground state properties and entanglement.

Contribution

It introduces the cooperation number as a key parameter and studies its effect on phase transitions and system observables in the Dicke model.

Findings

Cooperation number significantly affects quantum phase transition behavior.

Variational states provide accurate approximations of the ground state.

Entanglement entropy varies with energy difference and coupling strength.

Abstract

We show how the use of variational states to approximate the ground state of a system can be employed to study a multi-mode Dicke model. One of the main contributions of this work is the introduction of a not very commonly used quantity, the cooperation number, and the study of its influence on the behavior of the system, paying particular attention to the quantum phase transitions and the accuracy of the used approximations. We also show how these phase transitions affect the dependence of the expectation values of some of the observables relevant to the system and the entropy of entanglement with respect to the energy difference between atomic states and the coupling strength between matter and radiation, thus characterizing the transitions in different ways.

Figures

Equations

$$H_D={\omega }_A{J}_z+\Omega a^{†}a-\frac{\gamma }{\sqrt{N}}\left(J_{-}+J_{+}\right)\left(a+a^{†}\right).$$

$$H={\omega }_A{J}_z+\sum _{ı=1}^k{\Omega }_{ı}{a}_{ı}^{†}{a}_{ı}-\frac{1}{\sqrt{N}}\sum _{ı=1}^k{\gamma }_{ı}\left(J_{-}+J_{+}\right)\left(a_{ı}+a_{ı}^{†}\right).$$

$$\left[a_{ı},a_{ȷ}^{†}\right]={\delta }_{ıȷ},\left[a_{ı},a_{ȷ}\right]=\left[a_{ı}^{†},a_{ȷ}^{†}\right]=0.$$

$$\left[S_{+},S_{-}\right]=2S_z,\left[S_z,S_{\pm }\right]=\pm S_{\pm }.$$

$$\begin{gathered} J_{\diamond }=S_{\diamond }\otimes I_2^{\otimes \left(N-1\right)}+I_2\otimes S_{\diamond }\otimes I_2^{\otimes \left(N-2\right)} \\ +\cdots +I_2^{\otimes \left(N-2\right)}\otimes S_{\diamond }\otimes I_2+I_2^{\otimes \left(N-1\right)}\otimes S_{\diamond } \end{gathered}$$

$$\left[J_{+},J_{-}\right]=2J_z,\left[J_z,J_{\pm }\right]=\pm J_{\pm }$$

$$\dim \left\{\mathrm{span}\left\{{\left\{|j,m⟩\right\}}_{|m|\le j}^{j=r,\dots ,\frac{N}{2}}\right\}\right\}=\sum _{j=r}^{\frac{N}{2}}\left(2j+1\right)\le 2^N.$$

$$|\bar{\alpha }⟩:=|{\alpha }_1⟩\otimes \cdots \otimes |{\alpha }_k⟩,$$

$$|{\alpha }_{ı}⟩:=e^{{\alpha }_{ı}{a}_{ı}^{†}-{\alpha }_{ı}^{\ast }{a}_{ı}}|0_{ı}⟩=e^{-\frac{{|{\alpha }_{ı}|}^2}{2}}\sum _{{\nu }_{ı}=0}^{\infty }\frac{{\alpha }_{ı}^{{\nu }_{ı}}}{\sqrt{{\nu }_{ı}!}}|{\nu }_{ı}⟩.$$

$$\begin{gathered} {|\xi ⟩}_j:={|\frac{\upsilon \tan \left|\upsilon \right|}{\left|\upsilon \right|}⟩}_j:=e^{\upsilon J_{+}-{\upsilon }^{\ast }{J}_{-}}|j,0⟩ \\ =\frac{1}{{\left(1+{\left|\xi \right|}^2\right)}^j}\sum _{m=0}^{2j}{\left(\genfrac{}{}{0ex}{}{2j}{m}\right)}^{\frac{1}{2}}{\xi }^m|j,m-j⟩. \end{gathered}$$

$$\begin{gathered} {|\bar{\alpha },{\xi }_j⟩}_{\pm }:={\mathcal{N}}_{\pm }{P}_{\pm }|\bar{\alpha }⟩\otimes {|\xi ⟩}_j \\ ={\mathcal{N}}_{\pm }\left(|\bar{\alpha }⟩\otimes {|\xi ⟩}_j\pm |-\bar{\alpha }⟩\otimes {|-\xi ⟩}_j\right), \end{gathered}$$

$$S_{\epsilon }:=-Tr\left\{{\rho }_1\log {\rho }_1\right\}=-Tr\left\{{\rho }_2\log {\rho }_2\right\},$$

$$F(\varphi ,\phi ):={\left|⟨\varphi |\phi ⟩\right|}^2.$$

$$\begin{gathered} {\mathcal{H}}_{j,CS}\left(q_{ı},p_{ı},\theta ,\varphi \right):=⟨\bar{\alpha }|\otimes {⟨\xi |}_{j}H|\bar{\alpha }⟩\otimes {|\xi ⟩}_j \\ =-j{\omega }_A\cos \theta +\sum _{ı=1}^k{\Omega }_{ı}\left(q_{ı}^2+p_{ı}^2\right)-\frac{4j}{\sqrt{N}}\sin \theta \cos \varphi \sum _{ı=1}^k{\gamma }_{ı}{q}_{ı}. \end{gathered}$$

$${\theta }_c=q_{{ı}_c}=p_{{ı}_c}=0,\text{mbox}for{\omega }_A\ge \frac{8j}{N}\sum _{ı=1}^k\frac{{\gamma }_{ı}^2}{{\Omega }_{ı}},$$

$$\left.\begin{array}[]{c}\cos\theta_{c}=\frac{N\omega_{A}}{8j}\left({\displaystyle\sum_{\imath=1}^{k}}\frac{\gamma_{\imath}^{2}}{\Omega_{\imath}}\right)^{-1},\\ \\ \phi_{c}=0,\pi,\\ \\ q_{\imath_{c}}=\frac{2j\gamma_{\imath}}{\Omega_{\imath}\sqrt{N}}\cos\phi_{c}\sin\theta_{c},\\ \\ p_{\imath_{c}}=0\end{array}\right\}\ \ \ \mbox{for}{\displaystyle\;\omega_{A}<\frac{8j}{N}{\displaystyle\sum_{\imath=1}^{k}}\frac{\gamma_{\imath}^{2}}{\Omega_{\imath}}}.$$

$${\mathcal{E}}_{CS}\left({\omega }_A,{\gamma }_{ı}\right)=\{\begin{array}{ccc}-j{\omega }_A& ,& \text{mbox}for\delta \ge 1\\ -\frac{j{\omega }_A}{2}\left(\frac{1}{\delta }+\delta \right)& ,& \text{mbox}for\delta <1,\end{array}$$

$${⟨J_z⟩}_{CS}\left({\omega }_A,{\gamma }_{ı}\right)=\{\begin{array}{ccc}-j& ,& \text{mbox}for\delta \ge 1\\ -j\delta & ,& \text{mbox}for\delta <1,\end{array}$$

$${⟨{\nu }_{ı}⟩}_{CS}\left({\omega }_A,{\gamma }_{ı},{\gamma }_{ȷ}\right)=\{\begin{array}{ccc}0& ,& \text{mbox}for\delta \ge 1\\ \frac{{\gamma }_{ı}^2}{{\Omega }_{ı}^2}\frac{j{\omega }_A}{2\varsigma }\left(\frac{1}{\delta }-\delta \right)& ,& \text{mbox}for\delta <1.\end{array}$$

$$\begin{gathered} {\mathcal{H}}_{j,SAS}\left(q_{ı},p_{ı},\theta ,\varphi \right):={⟨\bar{\alpha },{\xi }_j|}_{+}H{|\bar{\alpha },{\xi }_j⟩}_{+} \\ =\left(\frac{1+E{\left(-\cos \theta \right)}^{2j-2}}{1+E{\left(-\cos \theta \right)}^{2j}}\right)\left(-j{\omega }_A\cos \theta \right) \\ +\left(\frac{1-E{\left(-\cos \theta \right)}^{2j}}{1+E{\left(-\cos \theta \right)}^{2j}}\right)\sum _{\ell =1}^k{\Omega }_{\ell }\left(q_{\ell }^2+p_{\ell }^2\right) \\ -\frac{4j}{\sqrt{N}}\sin \theta \sum _{\ell =1}^k\left\{\frac{\cos \varphi {\gamma }_{\ell }{q}_{\ell }+E{\left(-\cos \theta \right)}^{2j-1}\sin \varphi {\gamma }_{\ell }{p}_{\ell }}{1+E{\left(-\cos \theta \right)}^{2j}}\right\}. \end{gathered}$$

$$\mathcal{E}_{SAS}\left(\omega_{A},\gamma_{\imath}\right)\\ =\left\{\begin{array}[]{cc}{\displaystyle-j\omega_{A}}&,\,\mbox{for}\;{\displaystyle\delta\geq 1}\\ \\ {\displaystyle-j\omega_{A}\left[\delta\left(\frac{1+\varepsilon\left(-\delta\right)^{2j-2}}{1+\varepsilon\left(-\delta\right)^{2j}}\right)\right.}\\ {\displaystyle\left.+\frac{1}{2}\left(\frac{1}{\delta}-\delta\right)\right]}&,\,\mbox{for}\;{\displaystyle\delta<1,}\end{array}\right.$$

$$\begin{gathered} {⟨J_z⟩}_{SAS}\left({\omega }_A,{\gamma }_{ı}\right) \\ =\{\begin{array}{ccc}-j& ,& \text{mbox}for\delta \ge 1\\ -j\delta \left(\frac{1+\epsilon {\left(-\delta \right)}^{2j-2}}{1+\epsilon {\left(-\delta \right)}^{2j}}\right)& ,& \text{mbox}for\delta <1,\end{array} \end{gathered}$$

$$\begin{gathered} {⟨{\nu }_{ı}⟩}_{SAS}\left({\omega }_A,{\gamma }_{ı},{\gamma }_{ȷ}\right) \\ =\{\begin{array}{ccc}0& ,& \text{mbox}for\delta \ge 1\\ \frac{{\gamma }_{ı}^2}{{\Omega }_{ı}^2}\frac{j{\omega }_A}{2\varsigma }\left(\frac{1}{\delta }-\delta \right)\left(\frac{1-\epsilon {\left(-\delta \right)}^{2j}}{1+\epsilon {\left(-\delta \right)}^{2j}}\right)& ,& \text{mbox}for\delta <1,\end{array} \end{gathered}$$

Full text

Introduction

Quantum phase transitions (QPTs) are informally seen as sudden, drastic changes in the physical properties of the ground state of a system at zero temperature due to the variation of some parameter involved in the modeling Hamiltonian. One model of particular interest for the study of such phenomena is the Dicke model key-1 , as it describes, in a simplified way (electric dipole approximation), the interaction between matter and electromagnetic radiation. In 1973, Hepp and Lieb key-2 ; key-3 , and Wang and Hioe key-4 first theoretically proved the existence of a second-order QPT in the Dicke model. Wang and Hioe also treated the multi-mode radiation case, where they reduce it to a single-mode case by using an effective coupling constant. To date, this QPT has been experimentally observed in a Bose-Einstein Condensate coupled to an optical cavity key-5 ; key-6 and it has been shown to be relevant to quantum information and quantum computing key-7 ; key-8 ; key-9 ; key-10 .

Even though the formal definition of a QPT requires us to compute the ground state’s energy as a function of any desired parameter in order to find its transition values, one of the main contributions of this work is to show how the QPT in the Dicke model influences the behavior of other quantities relevant to the system, thus characterizing the transition in different, simpler ways.

Quantum Phase Transitions

The formal definition of the concept of “quantum phase” that we will be using throughout this paper is that of an open region $\mathcal{R}\subseteq\mathbb{R}^{\ell}$ where the ground state’s energy $\mathcal{E}_{0}$, as a function of $\ell$ parameters involved in the modeling Hamiltonian, is analytic. Thus a QPT is identified by the boundary $\partial\mathcal{R}$ of the region at which $\frac{\partial^{n}\mathcal{E}_{0}}{\partial x^{n}}$ is discontinuous for some $n$ (known as the order of the transition).

Notice that in the previous definition, for the sake of generality, we did not consider the thermodynamic limit, as it has been shown that interesting phenomena regarding QPTs occur even for a finite number of particles key-11 ; key-12 .

Modeling Hamiltonian

The Hamiltonian (Dicke’s Hamiltonian) describing the interaction, in a dipolar approximation, between N two-level identical atoms (same energy difference between the two levels) and one-mode of an electromagnetic field in an ideal cavity, has the expression ($\hbar=1$)

[TABLE]

Here, $\omega_{A}$ is the energy difference between the atomic levels, $\Omega$ is the frequency of the field’s mode, $\gamma$ is the dipolar coupling constant, $J_{z}$, $J_{-}$, $J_{+}$ are the collective spin operators and $a$, $a^{\dagger}$ are the annihilation and creation operators of the harmonic oscillator. The multi-mode Hamiltonian is obtained summing over the number $k$ of modes key-13 , and has the expression

[TABLE]

The $k$ modes of the electromagnetic field are described in terms of annihilation and creation operators for each mode $a_{\imath}$, $a_{\imath}^{\dagger}$, acting on the tensor product of $k$ copies of the Fock space $\bigotimes\limits_{\imath=1}^{k}\mathcal{F}_{\imath}$ and satisfying the commutation relations

[TABLE]

A two-level atom is described using the $\frac{1}{2}$-spin matrices $S_{z}=\frac{1}{2}\sigma_{z}$, $S_{\pm}=\frac{1}{2}\left(\sigma_{x}\pm i\sigma_{y}\right)$ ($\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$ being the Pauli matrices), which act on a two-dimensional complex Hilbert space $\mathbb{C}^{2}$ and satisfy the commutation relations

[TABLE]

When considering a system of N two-level atoms, we use the collective spin operators $J_{z}$, $J_{-}$, $J_{+}$ defined as

[TABLE]

where $I_{2}$ is the identity operator on $\mathbb{C}^{2}$ and $\diamond\in\left\{z,-,+\right\}$. These collective spin operators satisfy the commutation relations

[TABLE]

and act, in principle, on the complex Hilbert space $\left(\mathbb{C}^{2}\right)^{\otimes N}$; however, working with this space is physically equivalent to studying a system of N fully distinguishable atoms, which we don’t usually have in the experimental setups used in the study of the QPT in the Dicke model. To overcome this issue, we must use the common set of eigenvectors $\left\{\left|j,m\right\rangle\right\}$ of the two commuting observables $J_{z}$ and $J^{2}=\frac{1}{2}\left(J_{+}J_{-}+J_{-}J_{+}\right)+J_{z}^{2}$, where the label $j$ is limited to the values $j\in\left\{r,r+1,\ldots,\frac{N}{2}\right\}$ ($r=0$ for even N and $r=\frac{1}{2}$ for odd N) and the label $m\in\mathbb{Z}$ is constricted by $\left|m\right|\leq j$. These vectors do not form a basis of $\left(\mathbb{C}^{2}\right)^{\otimes N}$ for $N>2$, as the dimension of their linear span is

[TABLE]

We will denote by $\mathcal{H}_{A}$ the subspace of $\left(\mathbb{C}^{2}\right)^{\otimes N}$ generated by the states $\left\{\left|j,m\right\rangle\right\}_{\left|m\right|\leq j}^{j=r,\ldots,\frac{N}{2}}$.

There are two main results concerning the states $\left\{\left|j,m\right\rangle\right\}_{\left|m\right|\leq j}^{j=r,\ldots,\frac{N}{2}}$ and the space $\mathcal{H}_{A}$: the first comes from noticing that $\left[H,J^{2}\right]=0$, which means that the label $j$ of the eigenvalues of $J^{2}$ remains constant during the system’s evolution; the second is the decomposition $\mathcal{H}_{A}=\bigoplus\limits_{j=r}^{\frac{N}{2}}\mathcal{H}_{j}$, where each $\mathcal{H}_{j}$ is the subspace of dimension $\mathrm{dim}\left\{\mathcal{H}_{j}\right\}=2j+1$ generated by the states $\left\{\left|j,m\right\rangle\right\}_{\left|m\right|\leq j}$ with a fixed $j$. In this treatment, in order to study indistinguishable atoms, we are ignoring the multiplicities $g(j)$ of the irreducible representations of $SU\left(2\right)$, i.e. the number of times that each $\mathcal{H}_{j}$ appears in the full decomposition $\left(\mathbb{C}^{2}\right)^{\otimes N}=\bigoplus\limits_{j=r}^{\frac{N}{2}}g\left(j\right)\mathcal{H}_{j}$.

To make it clear that the space $\mathcal{H}_{A}$ is the one we must work with when indistinguishable atoms are considered we should inquire into the physical interpretation of the labels $j$ and $m$. In order to give a physical interpretation to the label $j$ we must notice that the energy of the atomic system is bounded by $\pm j\omega$ independently of the number of atoms $N$ (but with the restriction $j\leq\frac{N}{2}$), this leads us to interpret the quantity $2j$ as the effective number of atoms in the system and define it as the cooperation number. To make the notion of the cooperation number more intuitive, Dicke, in his original paper key-1 , compares a state with $j=0$, which exists only for an even number of atoms, with a classical system of an even number of oscillators swinging in pairs oppositely phased. The interpretation of the label $m$ is clear from the definition of $J_{z}$: $m=\frac{1}{2}\left(n_{e}-n_{g}\right)$, where $n_{e}$ and $n_{g}$ are the number of atoms in the excited and ground states, respectively.

In this paper we restrict our analysis to the space $\mathcal{H}_{A}$, as it allows us to choose $j$ as an initial condition (which will remain constant) and work in $\mathcal{H}_{j}$, where the atoms are indistinguishable.

Methodology

There have been various contributions to the study of the phase transition in the Dicke model (and other two-level models) key-14 ; key-15 ; key-16 ; key-17 ; key-18 and different approaches such as Husimi function analysis key-19 , entropic uncertainty relations key-20 and energy surface minimization key-21 ; key-22 ; key-23 ; key-24 ; key-25 ; key-26 , have been used for its investigation.

In this work we use the energy surface minimization method, which consists on minimizing the surface that is obtained by taking the expectation value of the modeling Hamiltonian with respect to some trial variational state. The strength of this method lies on the choice of the trial state, as it is the latter, after minimization, the one that will be modeling the ground state of the system.

Here we take a variational approach for both matter and radiation fields, and show how to calculate the QPT of the system modeled by the Hamiltonian $H$ given in eq. (2) via four means:

Using a tensor product of Heisenberg-Weyl HW(1) coherent states for each mode of the electromagnetic field and SU(2) coherent states for the atomic field as trial states, and analytically minimizing the obtained energy surface with respect to its parameters. 2. 2.

Using a projection operator on HW(1) coherent states and SU(2) coherent states to obtain trial states that preserve the parity symmetry of the Hamiltonian with respect to the total excitation number of the system (symmetry adapted states), and numerically minimize the obtained energy surface with respect to its parameters. 3. 3.

Using symmetry adapted states, as in (2) above, to obtain the energy surface and “minimize” it with the minimizing parameters obtained in (1) above, thus allowing us to have analytic expressions for the ground state. 4. 4.

Numerically diagonalizing the Hamiltonian, which gives us the exact quantum solution.

Coherent states (CS)

For each mode of the electromagnetic field the annihilation and creation operators $a_{\imath}$ and $a_{\imath}^{\dagger}$, appearing in the modeling Hamiltonian $H$, satisfy the commutation relations (3) of the Lie algebra generators of the Heisenberg-Weyl group HW(1); hence, a natural choice of a trial state for the radiation field is a tensor product of $k$ (number of modes) coherent states of HW(1)

[TABLE]

where each $\left|\alpha_{\imath}\right\rangle$ is defined as

[TABLE]

Furthermore, the commutation relations of the collective spin operators $J_{-}$, $J_{+}$ and $J_{z}$ (6) are the same as the ones of the Lie algebra generators of the special unitary group SU(2). Thus, analogously as for the radiation field, we use the coherent states of SU(2)

[TABLE]

as trial states for the matter field.

Symmetry adapted states (SAS)

The modeling Hamiltonian we are considering has a parity symmetry given by $\left[e^{i\pi\Lambda},H\right]=0$, where $\Lambda=\sqrt{J^{2}+\frac{1}{4}}-\frac{1}{2}+J_{z}+{\displaystyle\sum_{\imath=1}^{k}}a_{\imath}^{\dagger}a_{\imath}$ is the excitation number operator with eigenvalues $\lambda=j+m+{\displaystyle\sum_{\imath=1}^{k}}\nu_{\imath}$. This symmetry allows us to classify the eigenstates of $H$ in terms of the parity of the eigenvalues $\lambda$; however, as states with opposite symmetry are strongly mixed by the CS defined in the previous section, we should then adapt this symmetry to the CS by projecting them with the operator $P_{\pm}=\frac{1}{2}\left(I\pm e^{i\pi\Lambda}\right)$, i.e.

[TABLE]

with $\mathcal{N}_{\pm}=\left(2\pm 2E\left(-\cos\theta\right)^{2j}\right)^{-\frac{1}{2}}$ the normalization factors for the even (+) and odd (-) states (where ${\displaystyle E=\exp\left\{-2{\displaystyle\sum_{\imath=1}^{k}}\left|\alpha_{\imath}^{2}\right|\right\}}$).

As we are interested in the ground state of the system, which has an even parity, we only focus on the state $\left|\bar{\alpha},\xi_{j}\right\rangle_{+}$.

Entropy of entanglement ($S_{\varepsilon}$)

Entropy of entanglement is defined for a bipartite system as the Von Neumann entropy of either of its reduced states, that is, if $\rho$ is the density matrix of a system in a Hilbert space $\mathcal{H}=\mathcal{H}_{1}\otimes\mathcal{H}_{2}$, its entropy of entanglement is defined as

[TABLE]

where $\rho_{1}=Tr_{2}\left\{\rho\right\}$ and $\rho_{2}=Tr_{1}\left\{\rho\right\}$.

Our Hamiltonian $H$ models a bipartite system formed by matter and radiation subsystems, which means that their entropy of entanglement can be used to see the influence of the QPT on its behavior; this we do below.

Fidelity between neighboring states ($F$)

Fidelity is a measure of the "distance" between two quantum states; given $\left|\phi\right\rangle$ and $\left|\varphi\right\rangle$ it is defined as

[TABLE]

Across a QPT the ground state of a system suffers a sudden, drastic change, thus it is natural to expect a drop in the fidelity between neighboring states near the transition. This drop has been, in fact, already shown to happen key-15 ; key-25 for the case $2j=N$. We study it here also, and its behavior with the cooperation number.

Results

Writing the complex labels $\alpha_{\imath}$ and $\xi$ as $\alpha_{\imath}=q_{\imath}+ip_{\iota}$ and $\xi=\tan\left(\frac{\theta}{2}\right)e^{i\phi}$, with $q_{\imath},p_{\imath}\in\mathbb{R}$, $\theta\in\left[0,\pi\right)$, $\phi\in\left[0,2\pi\right)$, the CS’s energy surface is obtained by taking the expectation value of the modeling Hamiltonian $H$ with respect to the state $\left|\bar{\alpha}\right\rangle\otimes\left|\xi\right\rangle_{j}$, and has the form

[TABLE]

The critical points which minimize it are then found to be

[TABLE]

[TABLE]

Substituting these values into (13) we obtain the energy of the coherent ground state as a function of the Hamiltonian parameters,

[TABLE]

where we have defined ${\displaystyle{\displaystyle\delta=\frac{N\omega_{A}}{8j\varsigma}}}$ with ${\displaystyle\varsigma={\displaystyle\sum_{\imath=1}^{k}}\frac{\gamma_{\imath}^{2}}{\Omega_{\imath}}}$. Using the information of this coherent ground state we also obtain the expectation values of the atomic relative population operator $J_{z}$ and of the number of photons of mode $\imath$ operator $\nu_{\imath}:=a_{\imath}^{\dagger}a_{\imath}$:

[TABLE]

[TABLE]

Analogously as for the CS’s energy surface, the SAS’s energy surface is obtained by taking the expectation value of the modeling Hamiltonian $H$ with respect to the state $\left|\bar{\alpha},\xi_{j}\right\rangle_{+}$, and has the more complicated form

[TABLE]

As a first approximation, we may substitute the critical values obtained for the CS’s energy surface into (17), we obtain the trial state which approximates the lowest symmetry-adapted energy state, and with respect to which we evaluate the expectation values of the observables $H$, $J_{z}$ and $\nu_{\imath}$:

[TABLE]

[TABLE]

[TABLE]

where ${\displaystyle\varepsilon=\exp\left\{{\displaystyle\frac{-j\omega_{A}\sigma}{\varsigma}}\right\}}$ with ${\displaystyle\sigma={\displaystyle\sum_{\imath=1}^{k}}\frac{\gamma_{\imath}^{2}}{\Omega_{\imath}^{2}}}$.

Of course, we can minimize eq. (17) numerically for the SAS and obtain the expectation value of the relevant matter and field observables.

In our numerical analysis we study the case with two modes of the radiation field, as it is the maximum number of orthogonal modes that can be present in a 3D cavity with the restrictions that the modes interact with the electric dipole moment of the atoms and to be in resonance with the frequency associated with the energy difference between the two levels of the atoms. This latter restriction is just considered to have the maximum transition probability between states.

For the exact quantum solution we must resort to numerical diagonalization of the Hamiltonian and use the lowest eigenstate to compute the expectation values of the relevant observables.

The results, properties of the ground state related to the CS, those related to the SAS using the critical points of the CS (which have the advantage of also providing analytical solutions), those of the SAS minimized numerically and those of the quantum solution through numerical diagonalization, are shown in figures 1 - 15 and are discussed below.

One advantage of having analytical solutions is, of course, that the order of the transition may be easily found. Equations (14) and (18) show a second-order QPT at $\delta=\frac{N\omega_{A}}{8j\varsigma}=1$ with the CS and SAS using CS’s minima (SASc) approximations. In figure 1 it can be seen that the data of the SAS using numerical minimization (SASn) has a small discontinuity (the QPT) at $\gamma_{2}\approx 1.485$ for $j=5$ and $\gamma_{2}\approx 1.015$ for $j=9$, while in figure 2 this discontinuity is at $\omega_{A}\approx 0.975$ for $j=5$ and $\omega_{A}\approx 1.965$ for $j=9$. Note that the SASn solution always approximate better the exact quantum result, as the cooperation number increases this approximation gets better, in fact, for $2j=18=N$ the loci of the separatrix between the normal and collective regions for the quantum and SASn solutions are indistinguishable (except in the zoomed inset). The true loci of the QPT may be found through the fidelity: figures 10 and 11 show the fidelity between neighboring states of the quantum solution, where the exact QPT is characterized by the minimum, which is localized at $\gamma_{2}\approx 1.550$ for $j=5$, $\gamma_{2}\approx 1.031$ for $j=9$ in figure 10; and $\omega_{A}\approx 0.817$ for $j=5$, $\omega_{A}\approx 1.870$ for $j=9$ in figure 11.

The discrepancies between the transition values of the SASc approximation and the exact quantum solution become obvious when looking at figures 12 and 13, where the fidelity between SASc and the quantum solution drops (and oscillates) in a vicinity of the separatrix. Therefore, we conclude that SASc offer a good approximation (with an analytic expression) to the exact quantum solution far from the QPT for low cooperation numbers, but as $j\rightarrow\infty$, the interval where the SASc fail to reproduce the correct behavior, becomes smaller.

Figures 14 and 15 show the fidelity drop at the separatrix for the SASn. The resemblance to figures 10 and 11 is uncanny, showing the benefits of restoring the Hamiltonian symmetry into the trial variational states. This improvement comes with the disadvantage of losing the analytic expression, but still has an advantage over the quantum solution: the computational time. SASn are obtained by numerically minimizing a real function, which is far easier to do (computationally speaking) than numerically diagonalizing the Hamiltonian matrix.

Figures 1 - 7 show the comparison between the different approximations to the ground state: CS, SASc, SASn and quantum solution. We show the behavior of $\mathcal{E}:=\left\langle H\right\rangle$, $\left\langle J_{z}\right\rangle$ and $\left\langle\nu_{\imath}\right\rangle$ as functions of the atomic frequency $\omega_{A}$ and one of the coupling constants $\gamma_{2}$, for different cooperation numbers. It can be noticed that the discontinuity in the second derivative of the energy (as modeled with CS and SASc) translates into a discontinuity in the first derivative of $\left\langle J_{z}\right\rangle$ and $\left\langle\nu_{\imath}\right\rangle$, thus characterizing the QPT by means of an abrupt change in the expectation values of the observables. In general, it can be observed that the four methods (CS, SASc, SASn and quantum solution) converge in the limit $\delta\rightarrow 0$, where the case $j\rightarrow\infty$ is particularly interesting as the interval around the QPT, where all the approximations fail to reproduce the correct behavior, becomes smaller.

It is worth mentioning the significance and importance of figures 7 and 16 as they show aspects of the multi-mode Dicke model which are not present in the single-mode case. In figure 7 it is shown how the different modes of radiation (orthogonal in principle) interact through the matter field, analogously as it occurs with different atoms interacting through the radiation field. On the other hand, figure 16 shows (pictorically) the phase diagrams of the two-mode system, in which it can be observed that any two points in the super-radiant region can be joined by a trajectory that does not cross the normal region, a characteristic that the single-mode system does not have.

Figures 8 and 9 show the comparison between SASc, SASn and the quantum solution for the entropy of entanglement $S_{\varepsilon}$ as a function of the atomic frequency $\omega_{A}$ and one of the coupling constants $\gamma_{2}$, using different cooperation numbers. A characterization of the QPT can be made by observing that the entropy of entanglement obtained using the quantum solution shows a maximum at the transition, an attribute that SASc and SASn approximations fail to reproduce.

Discussion and Conclusions

From figures presented we conclude that SASc offer a good approximation (with an analytic expression) to the exact quantum solution far from the QPT for low cooperation numbers, but as $j\rightarrow\infty$, the interval where the SASc fail to reproduce the correct behavior, becomes smaller.

On the other hand, the SASn provide a better approximation to the quantum solution. This improvement comes with the disadvantage of losing the analytic expression, but still has the advantage over the quantum solution of the computational time. SASn are obtained by numerically minimizing a real function, which is far easier to do (computationally speaking) than numerically diagonalizing the Hamiltonian matrix.

A characterization of the QPT can be made by looking at the entropy of entanglement; that obtained using the quantum solution shows a maximum at the transition, an attribute that SASc and SASn approximations fail to reproduce.

The behavior of the expectation values of the relevant observables of the system $\left\langle H\right\rangle$, $\left\langle J_{z}\right\rangle$ and $\left\langle\nu_{\imath}\right\rangle$, is also affected by the QPT (figures 1 - 7), thus allowing us to characterize the QPT by means of its influence over the observables. In general, it can be observed that the four methods (CS, SASc, SASn and quantum solution) converge in the limit $\delta\rightarrow 0$, where the case $j\rightarrow\infty$ is particularly interesting as the interval around the QPT, where all the approximations are weaker, becomes smaller.

In conclusion, we have shown how the use of variational states to approximate the ground state of a system can be useful to characterize the QPT in a multi-mode Dicke model using the expectation value of the observables relevant to the system and the entropy of entanglement between matter and radiation. We have also introduced a not very commonly used dependence: the cooperation number, showing its influence over the behavior of the system, paying particular attention to the QPT and the accuracy of the used approximations. Some aspects of the multi-mode Dicke model which are not present in the single-mode case were also briefly discussed.

Acknowledgements.

We thank R. López-Peña, O. Castaños and S. Cordero for their comments and discussion. This work was partially supported by DGAPA-UNAM under project IN101217. L. F. Q. thanks CONACyT-México for financial support (Grant #379975).