# Phase transitions of the anisotropic Dicke model

Pragna Das  Indian Institute of Science Education and Research Bhopal 462066
India  Devendra Singh Bhakuni  Department of Physics, Ben-Gurion University of
the Negev, Beer-Sheva 84105, Israel  Auditya Sharma  Indian Institute of
Science Education and Research Bhopal 462066 India

###### Abstract

We systematically analyze the various phase transitions of the anisotropic
Dicke model that is endowed with both rotating and counter-rotating light-
matter couplings. In addition to the ground state quantum phase transition
(QPT) from the normal to the super-radiant phase, the anisotropic Dicke model
also exhibits other transitions namely the excited state quantum phase
transition (ESQPT), ergodic to non-ergodic transition (ENET) and the
temperature dependent phase transition. We show that these phase transitions
are profitably studied not only with the standard consecutive level spacing
ratio, but also with the aid of various eigenvector quantities such as von
Neumann entanglement entropy, the participation ratio, multifractal dimension
and mutual information. For ENET, both the statics and dynamics of the
participation ratio offer a consistent and useful picture. An exciting finding
from our work is that the ESQPT and the ENET are closely related to each
other. We show this with the aid of two characteristic energies in the
spectrum corresponding to jumps in von Neumann entropy.

##  I Introduction

The Dicke model [1, 2, 3, 4], which is paradigmatic within the field of cavity
quantum electrodynamics, describes the interaction between $N$ atoms and a
single-mode bosonic field via a dipole coupling strength. In the thermodynamic
limit ($N\to\infty$), the model shows a quantum phase transition from the
normal phase (NP) to the super-radiant phase (SP) [1, 5, 6, 7, 8, 9, 10] at
some critical coupling strength. Along with this quantum phase transition
(QPT), the Dicke model also exhibits two other distinct phase transitions
[11], namely the excited state quantum phase transition (ESQPT) [12, 13, 6,
14, 15, 16] and the thermal phase transition (TPT) [17, 18, 19, 20]. While the
first occurs at finite energy when the coupling strength is sufficiently
large, the second, on the other hand, occurs at a finite temperature [18].
Some of these quantum phase transitions were observed experimentally in Bose-
Einstein condensates [21] and quantum cavity systems [22]. The physics of
systems with light-matter interactions has enjoyed a great deal of interest in
recent times, triggered by a number of experimental works [23, 24, 25].

A generalized version of the Dicke model, namely the anisotropic Dicke model
[26, 27, 28, 29, 30, 31, 32, 33, 34, 16, 35] (ADM), where the coupling
strengths corresponding to the rotating and counter-rotating terms are
different, has gained traction in recent times. While a huge body of
literature has been built around the Dicke model [36, 37, 38, 39, 40, 8], the
anisotropic model has received relatively less attention. The asymmetry in the
coupling brings some novel features in addition to the existing properties of
the Dicke model. One such novel feature that has generated considerable
excitement is that the ADM not only exhibits the normal to super-radiant
quantum phase transition, but also an ergodic-to-nonergodic transition (ENET)
[31, 35]. The model is integrable in the limit where either one of the
couplings is zero. Moreover, while the ground state properties show the
normal-to-super-radiant phase transition, the excited states show signatures
of non-ergodicity [35]. It was also argued that the transition from the normal
to the super-radiant phase is quite different in comparison to the ergodic to
non-ergodic transition [31]. In the present work, focusing on eigenvector
properties, we show that the normal-to-super-radiant phase transition
corresponds to the ground state undergoing a localized-to-multifractal
transition. On the other hand the ergodic-to-nonergodic transition corresponds
to the middle excited state undergoing a delocalized-to-multifractal
transition.

Phase transitions are often characterized by quantum information tools such as
entanglement entropy [2, 41, 42, 43, 44], mutual information [45, 46, 47, 48]
and so on. These quantities have proven useful not only to mark a variety of
phase transitions [41, 49, 42, 43, 44, 50, 51], but in diverse other contexts
[52, 53, 54] where quantum correlations have an important role. Moreover, some
of these quantities are directly related to experimentally measurable
observables and have proven to be useful markers of the phase transitions of
the Dicke model [40]. Thus, a study of these quantities in the context of
phase transitions is of theoretical interest with potential to connect with
experimental work.

In this work, we explore the various phase transitions of the anisotropic
Dicke model and their dependence on the asymmetric coupling strengths. First,
we study the behavior of the well known quantum phase transition of the ADM
using quantum information measures. With the aid of the ground state energy,
average photon number, inverse participation ratio and its scaling with the
Hilbert-space dimension, we highlight that the normal to super-radiant phase
transition is reminiscent of a localization to multifractal phase transition.
Next, we highlight the emergence of the excited state quantum phase transition
and the temperature dependent phase transition and their dependence on the
coupling parameters. While these transitions have been studied extensively for
the Dicke model [10, 38, 55, 18, 12, 40], they are also prominently present in
the ADM [31, 35]. We study the temperature-dependent phase transition as a
function of the rotating and counter-rotating coupling strengths with the aid
of an old analytical result [26] for the transition temperature. Similar to
the Dicke model [10], we find that mutual information between two spins offers
a clear signature of the thermal phase transition, which is benchmarked
against the analytical expression for the critical temperature which has
already been worked out in the literature [26, 30].

Our main result is to show that the ESQPT is profitably studied with the help
of von Neumann entanglement entropy between the bosons and the spins, the
average level spacing ratio and two characteristic energies that define a
central band in the super-radiant phase. With the help of participation ratio
and multifractal dimension, we show that the middle excited state exhibits
multifractal behavior in the nonergodic phase. Thus the middle excited state
behaves in stark contrast to the ground state which shows a change from
localized to multifractal behavior as one goes from the normal to the super-
radiant phase. The correspondence between the multifractal nature of the
middle excited state and the nonergodic phase is also captured dynamically
when we study the participation ratio in a quench dynamical protocol. Another
exciting finding of our work is that the excited state quantum phase
transition and the non-ergodic to ergodic transition are closely related.
Specifically, we find that the phase diagram obtained by keeping track of the
size of the jumps in von Neumann entropy for the different eigenstates of the
system (which carry signatures of the ESQPT), closely resembles the ENET phase
diagram. While the ESQPT in the anisotropic Dicke model was explored in Ref.
[16], it is mainly focused on the properties of the eigenvalues. The authors
of Ref. [32] studied the atom-field and atom-atom entanglement in the
anisotropic Dicke model in which the couplings are restricted to $g_{1}\geq
g_{2}$. Our work considers a more general parameter regime with arbitrary non-
negative $g_{1}$ and $g_{2}$ like in Buijsman $et$ $al$ [31].

The organization of the paper is as follows. In the Sec. II, we introduce the
model Hamiltonian. Next we discuss the various phase transitions (QPT, ESQPT,
ENET, TPT) exhibited by the anisotropic Dicke model and their characterization
via tools from quantum information theory. While the QPT, ESQPT, and ENET are
covered in Sec. III, Sec. IV is dedicated to the TPT. Finally in Sec. V we
summarize the main results.

##  II Model Hamiltonian

The Hamiltonian consists of a single-mode bosonic field coupled to $N$ atoms
with anisotropic couplings of the rotating and counter-rotating terms:

$\displaystyle\mathcal{H}=\omega
a^{\dagger}a+\omega_{0}J_{z}+\frac{g_{1}}{\sqrt{2j}}(a^{\dagger}J_{-}+aJ_{+})$
$\displaystyle+\frac{g_{2}}{\sqrt{2j}}(a^{\dagger}J_{+}+aJ_{-}).$ (1)

Here the operators $a$ and $a^{\dagger}$ are bosonic annihilation and creation
operators respectively, and
$J_{\pm,z}=\sum_{i=1}^{2j}\frac{1}{2}\sigma_{\pm,z}^{(i)}$ are angular
momentum operators of a pseudospin with length $j$, composed of $N=2j$ spin-
$\frac{1}{2}$ atoms described by Pauli matrices $\sigma_{\pm,z}^{(i)}$ acting
on site $i$. The commutation relations (in units where $\hbar=1$) between the
various operators are as follows:

$\displaystyle\left[a,a^{\dagger}\right]=1,\left[J_{z},J_{\pm}\right]=\pm
J_{\pm},\left[J_{+},J_{-}\right]=2J_{z}.$ (2)

The basis of the full Hilbert space of the system is
$\\{|n\rangle\otimes|j,m\rangle\\}$ where $|n\rangle$ are the number states of
the field satisfying $a^{\dagger}a|n\rangle=n|n\rangle$ and $|j,m\rangle$ are
the Dicke states satisfying $J_{\pm}|j,m\rangle=\sqrt{j(j+1)-m(m\pm 1)}|j,m\pm
1\rangle$. $\omega$ is the single-mode frequency of the bosonic field while
$\omega_{0}$ is the level splitting of the atoms. $g_{1}$ and $g_{2}$ are the
time independent coupling strengths corresponding to the rotating and counter-
rotating light-matter interaction terms. In the thermodynamic limit, the
system shows a second order quantum phase transition from normal to super-
radiant phase at $g_{1}+g_{2}=1$. For $g_{1}+g_{2}<1$, the system is in the
normal phase with $\langle a^{\dagger}a\rangle/j\approx 0$ (the bosonic mode
is microscopically excited) and for $g_{1}+g_{2}>1$, it is in the super-
radiant phase with a positive value of $\langle a^{\dagger}a\rangle/j$
(bosonic mode of the system is macroscopically excited). Here the expectation
value is calculated with respect to the ground state of the system
Hamiltonian. The ADM possesses a parity symmetry $[H,\Pi]=0$ with
$\Pi=\exp(i\pi[a^{\dagger}a+J_{z}+j])$ having eigenvalues $\pm 1$ [31]. Here,
we restrict ourselves to the $+1$ eigenvalue and an even atom number $N$ and
the symmetric subspace where $j=\frac{N}{2}$. Hence the $(N+1)$ values that
$m$ can take are: ($-\frac{N}{2},\ .\ .\ .\ ,\ 0,\ .\ .\ .\,\ \frac{N}{2}$).
For our numerics we truncate the boson number to take the values $n=0,\ 1,\ .\
.\ .\,\ n_{\text{max}}$. Thus the total Hilbert space dimension for the
truncated model is: $N_{D}=(n_{\text{max}}+1)(N+1)$. We have checked that all
our numerical results (for the specified $n_{\text{max}}$) are robust against
further increase in $n_{\text{max}}$, thus the truncation is carried out in
such a way that our numerical results are reliable. Furthermore, we fix
$\omega_{0}=1$ as the dimension of energy and the other observables are
calculated in units of $\omega_{0}$.

In the thermodynamic limit (when the number of atoms $N\to\infty$) the model
is analytically solvable using the Holstein-Primakoff representation [56, 57]
of the angular momentum operators $J_{z}=(b^{\dagger}b-j)$,
$J_{+}=b^{\dagger}\sqrt{2j-b^{\dagger}b}$, $J_{-}=J_{+}^{\dagger}$. Here $b$
and $b^{\dagger}$ are bosonic operators that convert the system Hamiltonian
into a two-mode bosonic problem. This allows us to obtain effective
Hamiltonians that are exact in the thermodynamic limit, by neglecting terms
from expansions of the Holstein-Primakoff square roots [3]. In the normal
phase $g_{1}+g_{2}<\sqrt{\omega\omega_{0}}$, the square roots can be expanded
directly and the effective Hamiltonian is

$\displaystyle\mathcal{H}^{(1)}=\omega
a^{\dagger}a+\omega_{0}(b^{\dagger}b-j)+g_{1}(a^{\dagger}b+ab^{\dagger})$
$\displaystyle+g_{2}(a^{\dagger}b^{\dagger}+ab),$ (3)

which is bilinear in the bosonic operators. In this representation, the parity
operator $\Pi$ becomes
$\Pi=\exp\Big{(}i\pi\Big{[}a^{\dagger}a+b^{\dagger}b\Big{]}\Big{)}$. In the
super-radiant phase $g_{1}+g_{2}>\sqrt{\omega\omega_{0}}$, both the field and
the atomic ensemble acquire macroscopic occupations and for that one needs to
displace the bosonic modes:

$a^{\dagger}\rightarrow c^{\dagger}+\sqrt{\alpha},\hskip
8.53581ptb^{\dagger}\rightarrow d^{\dagger}-\sqrt{\beta},$

where the undetermined parameters $\alpha$ and $\beta$ are of order $O(j)$.
Now considering the thermodynamic limit, the Hamiltonian can be written as:

$\displaystyle\mathcal{H}^{(2)}$ $\displaystyle=$ $\displaystyle\omega
c^{\dagger}c+\left[\omega_{0}+\frac{(g_{1}+g_{2})}{k}\sqrt{\frac{\alpha\beta
k}{2j}}\right]d^{\dagger}d-\left[(g_{1}+g_{2})\frac{\beta
k}{2j}-\omega\sqrt{\alpha}\right]\left(c^{\dagger}+c\right)+\left[\frac{2(g_{1}+g_{2})}{k}\sqrt{\frac{\alpha
k}{2j}}(j-\beta)-\omega_{0}\sqrt{\beta}\right]$ (4)
$\displaystyle\times\left(d^{\dagger}+d\right)+\frac{(g_{1}+g_{2})}{4k^{2}}\sqrt{\frac{\alpha\beta
k}{2j}}(2k+\beta)\left(d^{\dagger}+d\right)^{2}+\frac{(g_{1}+g_{2})\beta}{2k}\sqrt{\frac{k}{2j}}\left(c^{\dagger}+c\right)\left(d^{\dagger}+d\right)+\sqrt{\frac{k}{2j}}\Big{[}g_{1}(c^{\dagger}d+cd^{\dagger})$
$\displaystyle+g_{2}(c^{\dagger}d^{\dagger}+cd)\Big{]}+\left[\omega_{0}(\beta-j)+\omega\alpha-2\frac{(g_{1}+g_{2})}{k}\sqrt{\frac{\alpha\beta
k}{2j}}\right]$

Figure 1:  (a) Ground state energy of ADM shows NP to SP QPT. (b) Average
number of boson (which is scaled by the pseudospin length $j$) is close to
zero for NP ($g_{1}+g_{2}<1$) and nonzero for SP ($g_{1}+g_{2}>1$). (c) The
ground state IPR is one in NP (localized) and close to zero in the SP
(delocalized) (d) Scaling of ground state $PR$ with the full Hilbert space
dimension $N_{D}$ at some point $g_{1}=1.2$, $g_{2}=0.8$, in the SP for ADM
for $n_{\text{max}}=600$ and changing the atom number: $N\in[20,40]$.
$PR_{\text{gs}}$ scales as $\sqrt{N_{D}}$. Here
$N_{D}=(n_{\text{max}}+1)(N+1)$ is increasing due to increase of atom number
$N$. In the inset we consider four $n_{\text{max}}$ values:
$n_{\text{max}}=200$ (blue solid circles), $300$ (purple stars), $400$ (green
triangles), $600$ (red squares) and changing the atom number: $N\in[20,40]$.
$PR_{\text{gs}}$ scales as $\sqrt{N}$. For each fixed $n_{\text{max}}$ dashed
lines indicate the fitting with $\sqrt{N}$. The four sets of $P_{\text{gs}}$
values are overlap. For (a)-(c) the parameters are: $\omega=\omega_{0}=1$,
$N=40$. The bosonic cut-off is set to be $n_{\text{max}}=200$, and this is
shown to be large enough. Energy is calculated in units of $\omega_{0}$ and we
fix $\omega_{0}=1$ throughout this paper.

where $k\equiv 2j-\beta$. This effective Hamiltonian in SP is also bilinear in
the bosonic operators. The global symmetry $\Pi$ is broken at the phase
transition and two new local symmetries appear corresponding to the operator
$\Pi^{(2)}=\exp\Big{(}i\pi\Big{[}c^{\dagger}c+d^{\dagger}d\Big{]}\Big{)}$ [3].

##  III QPT, ESQPT and ENET

In this section we discuss three types of transitions in the anisotropic Dicke
model: QPT, ESQPT and ENET, in separate subsections. To do the numerics we
perform exact diagonalization of the system Hamiltonian. Due to the bosonic
mode, the Hilbert space dimension of the ADM is infinite dimensional, however
for numerics one has to truncate the Hilbert space by cutting off the bosonic
mode at some finite $n_{\text{max}}$. We checked that our results remain
robust on increasing $n_{\text{max}}$ (see the Appendix).

###  III.1 Quantum phase transition (QPT)

In Fig. 1 we show the ground state energy density (energy scaled by the atom
number, $E_{\text{gs}}/N$) of the system Hamiltonian as a function of $g_{1}$
and $g_{2}$. While in the normal phase, the energy is almost constant close to
$=-0.5$ (which is very similar to the value seen for the isotropic Dicke
model), the super-radiant phase has a broad energy spectrum with the density
ranging from $-4.5\leq\frac{E_{gs}}{N}\leq-0.5$. This clearly distinguishes
the normal and super-radiant phases. Furthermore, the mean photon number given
by the operator $\langle a^{\dagger}a\rangle/j$ [3] is almost zero in the
normal phase, while in the super-radiant phase it has a non-zero value with a
continuous change across the transition line $g_{1}+g_{2}=1$ (Fig. 1), and
indicates a second order phase transition.

Finally, we study the ground state properties by looking at the degree of
localization using the (inverse) participation ratio. The participation ratio
[58, 59] (PR) of an eigenstate
$|\psi\rangle=\sum_{j}^{N_{D}}\psi_{j}|j\rangle$ (where $N_{D}$ is the Hilbert
space dimension) is defined as:

$P=\frac{1}{\sum_{j=1}^{N_{D}}|\psi_{j}|^{4}}.$ (5)

The inverse of the participation ratio, called the inverse participation ratio
($\text{IPR}=1/\text{PR}$) is often a useful measure in its own right. From
the inset of Fig. 1, it is clear that, in the normal phase, IPR is close to
unity suggesting that the ground state is localized. On the other hand, a
careful study of the scaling of PR in the super-radiant phase with the
Hilbert-space dimension $N_{D}$ reveals interesting features. The PR scales as
PR $\sim\sqrt{N_{D}}$ and suggests that the ground state in the super-radiant
phase exhibits multi-fractal behavior (see Fig. 1).

In Fig. 1 we fix the parameters to be $g_{1}=1.2,g_{2}=0.8$ to show a
representative example in the SP. While the bosonic cut-off is fixed at
$n_{\text{max}}=600$ we vary the atom number $N$ and hence the total Hilbert
space dimension $N_{D}$ of the truncated model also varies. The change in
$N_{D}$ is entirely due to the change in atom number $N$, since
$n_{\text{max}}$ remains fixed. In this figure the red squares are the data
for $P_{\text{gs}}$ whereas the red dashed line denotes the fit to the
functional form $N_{D}^{0.5}$. In the inset of Fig. 1 we choose four different
$n_{\text{max}}$ values: $n_{\text{max}}=200$ (blue solid circles), $300$
(purple stars), $400$ (green triangles), $600$ (red squares) and for each
fixed $n_{\text{max}}$ we show data as the atom number $N$ (and hence the
total Hilbert space dimension $N_{D}$) varies exactly like in the main figure.
The four sets of data corresponding to different $n_{\text{max}}$ exactly
overlap with each other. The dashed lines denote the data fitting of
$P_{\text{gs}}$ with $N^{0.5}$. Hence we can conclude that for a fixed atom
number $N$, $P_{\text{gs}}$ is independent of $n_{\text{max}}$, and only
depends on $N$.

###  III.2 Excited state quantum phase transition (ESQPT)

Figure 2: (a), (e), and (i) Von-Neumann entanglement entropy as a function of
the eigenstate energies of the ADM Hamiltonian at the point $g_{1}=1.0,\
g_{2}=1.1$ [for (a) and (e)] and at the point $g_{1}=4,\ g_{2}=2$ [for (i)].
In the figure the orange dash-dotted line denotes $E_{\text{lower}}$ [given in
Eqn. 7] and the blue dashed line is for $E_{\text{upper}}$ (given in Eqn. 8).
(b), (f), and (j) $\chi_{\text{lower}}$ [given in Eqn. 9] as a function of
$g_{1}$ and $g_{2}$, (c), (g), and (k) $\chi_{\text{upper}}$ [given in Eqn.
10] as a function of $g_{1}$ and $g_{2}$. (b) and (c) shows similar phase
transition from non-ergodic to ergodic phase, which suggest that ESQPT is
related to ENET. (d), (h), and (l) Consecutive level spacing ratio, $\langle
r\rangle$ of the energies lying between $E_{\text{normal}}^{0}$,
$E_{g_{1},g_{2}=0}^{\text{max}}$. For all the above figures atom number,
$N=40$, and we take the bosonic cut-off: $n_{\text{max}}=200$. Top panel is
for resonant case $\omega=\omega_{0}=1$, middle and bottom panels are for two
off-resonant cases: $\omega=\frac{\omega_{0}}{4}$ and $\omega=4\omega_{0}$
respectively.

We now look at the properties of the excited states and find that similar to
the Dicke model, the anisotropic Dicke model also exhibits an excited state
quantum phase transition in the super-radiant phase (for $g_{1}+g_{2}>1$).
While the literature on ESQPT has mainly considered eigenvalue properties [60,
40, 38, 12, 13, 14, 15], we showned in our recent paper [10] the profitability
of studying eigenvector properties such as the von Neumann entanglement
entropy [2], the average bosonic number [3], concurrence [61, 62] and
participation ratio [58, 59]. We also aruged that there is not only a lower
cut-off energy but also an upper cut-off and energies between these two cut-
offs behave differently from the upper and lower bands. In the Dicke model,
the lower cut-off energy [60, 40] is the ground state energy at the critical
coupling strength $g_{c}$ and the upper cut-off energy [10] is the maximum
energy at $g=0$ (for finite $n_{\text{max}}$). On the other hand, in the ADM,
we find that the lower cut-off energy is around the ground state energy of the
system on the critical line of quantum phase transitions ($g_{1}+g_{2}=1$)
while the upper cut-off energy is the maximum energy at $g_{1}=g_{2}=0$ (for
finite $n_{\text{max}}$). Energies below (above) the lower (upper) cut-off
form the lower (upper) energy band and the energies in between the two cut-
offs form the central band.

The entanglement entropy between the spins and the bosons is simply the von
Neumann entropy of the reduced density matrix of the spins:

$S_{\text{spins}}=-\text{Tr}[\rho_{\text{spins}}\ln(\rho_{\text{spins}})]$ (6)

where $\rho_{\text{spins}}$ is the reduced density matrix of the spins
obtained by tracing over the bosonic degrees of freedom. In Fig. 2 we show the
von Neumann entanglement entropy (between spins and bosons) as a function of
the energy eigenvalues at some point $g_{1}=1.0$ and $g_{2}=1.1$. We observe
from this representative plot that there are three characteristic parts: (i)
an increasing part upto some energy, followed by (ii) a plateau region, and
then (iii) a region where VNEE decreases. To quantify the beginning and the
end of the plateau region, we define two characteristic energies
$E_{\text{lower}}$ and $E_{\text{upper}}$ as:

$E_{\text{lower}}=\left[\frac{\sum_{n=0}^{\frac{N_{D}}{2}}E_{n}|\Delta
S_{n}|}{\sum_{n=0}^{\frac{N_{D}}{2}}|\Delta S_{n}|}\right]$ (7)

and

$E_{\text{upper}}=\left[\frac{\sum_{n=\frac{N_{D}}{2}}^{N_{D}}E_{n}|\Delta
S_{n}|}{\sum_{n=\frac{N_{D}}{2}}^{N_{D}}|\Delta S_{n}|}\right]$ (8)

where $\Delta S_{n}=S_{n+1}-S_{n}$ is the VNEE difference between that of the
$(n+1)^{\text{th}}$ eigenstate and the $n^{\text{th}}$ eigenstate, and $E_{n}$
is the $n^{\text{th}}$ energy. The above energies are obtained by using the
jumps in the VNEE as weights. The change in VNEE is taken as weights for the
different energies, and we would thus expect these quantities to signal the
two ends of the plateau region. In Fig. 2, these two energies are marked by
two vertical lines: the orange dash-dotted line denotes $E_{\text{lower}}$
whereas the blue dashed line denotes $E_{\text{upper}}$ and they more or less
match with the two end energies of the plateau region. Now we define two more
new quantities, $\chi_{\text{lower}}$ and $\chi_{\text{upper}}$:

$\chi_{\text{lower}}=\frac{E_{\text{lower}}}{E_{\text{normal}}^{0}}$ (9)
$\chi_{\text{upper}}=\frac{E_{\text{upper}}}{E_{g_{1},g_{2}=0}^{\text{max}}}$
(10)

where $E_{\text{normal}}^{0}$ is the minimum energy in the normal phase
($g_{1}+g_{2}<1$) and $E_{g_{1},g_{2}=0}^{\text{max}}$ is the maximum energy
at $g_{1}=g_{2}=0$ (for finite $n_{\text{max}}$). We plotted these two
quantities as a function of $g_{1}$ and $g_{2}$ in Fig. 2 and 2. Remarkably
this reveals a clear visual correlation between the ESQPT and what is called
the ergodic to non-ergodic phase transition (ENET) [31], which we describe in
greater detail in the following subsection. Moreover we notice that, along the
diagonal line [it is more clear in Fig. 2] there is a relatively dark line
which indicates that the symmetric Dicke model is special. A quantitative way
of identifying the central band of energies corresponding to the plateau
region in Fig. 2 is to consider the energies that lie between
$\chi_{\text{lower}}=1$ and $\chi_{\text{upper}}=1$, i.e
$E_{\text{lower}}=E_{\text{normal}}^{0}$ and
$E_{\text{upper}}=E_{g_{1},g_{2}=0}^{\text{max}}$. We perform a study of the
average level spacing ratio [63] $\langle r\rangle$ for the energies lying
between $E_{\text{lower}}$ and $E_{\text{upper}}$. Let ${s_{n}}$ denote the
level spacing between two consecutive energies $E_{n+1}$ and $E_{n}$, then the
$\langle r\rangle$ is defined as the average over $n$ of the ratio of
consecutive level spacings:

$r_{n}=\frac{\text{min}(s_{n-1},\hskip
2.84526pts_{n})}{\text{max}(s_{n-1},\hskip 2.84526pts_{n})}.$ (11)

For an ergodic system, the value of $\langle r\rangle=0.53$ and the
probability distribution of the consecutive gaps shows Wigner-Dyson statistics
[64] while for a non-ergodic system it is $\langle r\rangle=0.386$ and the
probability distribution becomes Poissonian.

Figure 3: Participation ratio and multifractal dimension, $D_{1}$ of ADM for
different eigenstate of the system Hamiltonian as a function of coupling
parameters $g_{1}$ and $g_{2}$. (a), (b), (c), (d) are the figure for
participation ratio corresponding to ground state (gs), tenth excited state,
$100{\text{th}}$ excited state, and middle excited state or thermal state (TS)
respectively. (e) $D_{1}$ for the ground state shows QPT from localization
(NP) to multifractal (SP) behavior, (f) $D_{1}$ for the tenth excited state,
(g) $D_{1}$ for the $100{\text{th}}$ excited state, (h) $D_{1}$ for the middle
excited state shows ENET from non-ergodic extended (multifractal) to ergodic
(delocalized) behavior. For all the above plots parameters are:
$\omega=\omega_{0}=1$, atom number $N=40$. We take the bosonic cut-off to be
$n_{\text{max}}=200$.

Figure 2 shows the phase diagram based on the level spacing ratio. It can be
seen that a small value of either $g_{1}$ or $g_{2}$ leads to a non-egrodic
phase where the eigenvalue statistics obey the Poisson distribution with
$\langle r\rangle\approx 0.386$. On the other hand when both the couplings are
significantly large, we see a crossover to Wigner-Dyson distribution where the
level-spacing ratio becomes $\langle r\rangle\approx 0.53$. It is worth
mentioning that this non-ergodic to ergodic transition is quite different from
the normal to super-radiant phase transition. We perform a careful analysis of
this phase transition in the next sub-section. We show that the Wigner-Dyson
behavior corresponds to the energy band that lies between the two cut-offs
[10] ($E_{\text{lower}}=E_{\text{normal}}^{0}$,
$E_{\text{upper}}=E_{g_{1},g_{2}=0}^{\text{max}}$).

So far we restricted ourselves to the resonant case $\omega=\omega_{0}$. Now
we also study ESQPT for the off-resonant cases considering
$\omega=\frac{\omega_{0}}{4}$ (middle panel of Fig. 2) and
$\omega=4\omega_{0}$ (bottom panel of Fig. 2). Studying the consecutive level
spacing ratio of the energy band sandwiched between $E_{\text{normal}}^{0}$
and $E_{g_{1},g_{2}=0}^{\text{max}}$ [see Figs. 2(h), and 2(l)], we see a
transition from $\langle r\rangle\approx 0.386$ to $\langle r\rangle\approx
0.53$ that is a non-ergodic to ergodic transition. From a careful observation
of Figs. 2(f), and 2(g) and 2(j), and 2(k) we infer that while a clear
correspondence between ENET and ESQPT is present for $g_{1}\neq g_{2}$ even in
the off-resonant scenario, the diagonal direction ($g_{1}=g_{2}$)
corresponding to the Dicke model exhibits special behavior [6].

###  III.3 Ergodic to non-ergodic transition (ENET)

####  III.3.1 Statics

As discussed in the previous subsection, different eigenstates play a role in
different phase transitions. While the ground state shows the normal to super-
radiant phase transition or equivalently from a localized phase to a
multifractal phase, the middle excited states exhibit a non-ergodic to ergodic
transition. Here, we study the phase diagram on the $g_{1},\ g_{2}$ plane of
the ADM for different eigenstates with the help of participation ratio and
multifractal dimension and explore the possibility of multifractal states in
the excited states.

While the participation ratio quantifies the degree of localization and
delocalization of a quantum state, a study of its scaling with the system size
offers further insights. When the Hilbert space dimension is large ($N_{D}$ is
large) the multifractal dimension [65, 66] can be represented as:

$\displaystyle D_{q}$ $\displaystyle=$
$\displaystyle\frac{1}{1-q}\frac{\ln\Big{(}\sum_{j=1}^{N_{D}}|\psi_{j}|^{2q}\Big{)}}{\ln(N_{D})}$
(12)

where $|\psi\rangle$ is an eigenstate of the Hamiltonian and
$S_{q}=\frac{1}{1-q}\ln\Big{(}\sum_{j=1}^{N_{D}}|\psi_{j}|^{2q}\Big{)}$ is
known as the $q$-dependent participation entropy. In the Shannon limit,
$S_{1}=-\sum_{j}|\psi_{j}|^{2}\ln\Big{(}|\psi_{j}|^{2}\Big{)}$, while the
$q=2$ participation entropy is connected to the usual participation ratio as
$S_{2}=\ln(P)$. For a perfectly delocalized state $S_{q}=\ln(N_{D})$, (when
$N_{D}$ is large) and hence $D_{q}=O(1)$, for all $q$. On the other hand for a
localized state $S_{q}=$ constant and $D_{q}=0$. In an intermediate situation,
wave functions are extended but non-ergodic with $S_{q}=D_{q}\ln(N_{D})$ where
$0<D_{q}<1$ and the state is multifractal in that particular basis.

Figure 4:  (a) Quench dynamics of participation ratio of ADM for
$\omega=\omega_{0}=1$, $N=20$ and $n_{\text{max}}=100$, at different times:
(a) $t=0.01$, (b) $t=0.2$, (c) $t=1$, (d) $t=1000$. Here the initial state is
the middle excited state of the decoupled Hamiltonian of the system
($H_{0}=\omega a^{\dagger}a+\omega_{0}J_{z}$). The time $t$ is in units of
$\omega_{0}^{-1}$ and we fix $\omega_{0}=1$ throughout this paper.

In Fig. 3 we show the phase diagrams of the ADM based on participation ratio
[Figs. 3(a) to 3(d)] and the multifractal dimension $D_{1}$ [Fig. 3(e) to
3(h)] for different eigenstates of the system Hamiltonian. For a fixed system
size, we pick a few states including the ground state and the middle excited
state. While for the ground state, we clearly see the normal to super-radiant
phase transition along the line $g_{1}+g_{2}=1$ as shown by the phase diagram
of the participation ratio [Fig. 1] and the multifractal dimension [Fig. 3],
as the states become more and more excited, it is the non-ergodic to ergodic
transition that is highlighted. Studying the PR of the middle excited state
[Fig. 3], we see a similar phase diagram as that of the level spacing ratio of
the system as a function of $g_{1}$ and $g_{2}$ [see Fig. 2], signifying a
transition from the non-ergodic phase to the ergodic phase. In the non-ergodic
phase the PR value is low whereas in the ergodic phase its value is relatively
higher.

In Figs. 3(e) to 3(h) we study the multifractal dimension $D_{1}$ for the
different eigenstates of the ADM as a function of coupling parameters $g_{1}$
and $g_{2}$ similar to the participation ratio. Here in Figure 3 we show the
nature of $D_{1}$, for the ADM ground state. In the NP, $D_{q}\approx 0$
whereas in the SP, $0<D_{1}<1$. This suggests a transition from a localized to
a multifractal phase. Fig. 3 shows $D_{1}$ for the middle excited state of the
ADM. In this figure we see the regions (depicted by the blue color) where
$0<D_{1}<1$ which behaves like a multifractal phase, whereas the region where
$D_{1}\approx 0.9$ (depicted by the white color), behaves as more like a
delocalized phase. Hence the middle excited state shows a transition from an
extended non-ergodic (multifractal) phase to an ergodic (delocalized) phase.

####  III.3.2 Dynamics

To study the quench dynamics of a closed quantum system, one prepares the
system in some eigenstate of the initial Hamiltonian $\mathcal{H}_{0}$. The
Hamiltonian is suddenly changed to
$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{1}$, and the system is allowed to
evolve under the corresponding unitary time evolution operator. Here, we take
the middle excited state of the decoupled Hamiltonian ($\mathcal{H}_{0}=\omega
a^{\dagger}a+\omega_{0}J_{z}$), as our initial state which can be written as
$|\psi_{\text{in}}\rangle=\sum_{\alpha}C_{\alpha}|\alpha\rangle$, with
$|\alpha\rangle=|n,j,m\rangle$ being a computational basis state and
$C_{\alpha}$ being the corresponding coefficient. The time-evolution of the
state is given by
$|\psi_{t}\rangle=e^{-i\mathcal{H}t}|\psi_{\text{in}}\rangle=\sum_{\alpha}C_{\alpha}(t)|\alpha\rangle$,
where $\mathcal{H}$ is the ADM Hamiltonian. To study the dynamical properties,
we calculate the dynamics of the participation ratio
$\text{PR}(t)=1/{\sum_{\alpha}|C_{\alpha}(t)|^{4}}$, at different times:
$t=0.01$, $0.2$, $1$, $1000$. From Fig. 4 we see that, for a very small
duration of time say $t=0.01$, the participation ratio has low value for all
$g_{1}$ and $g_{2}$. As we evolve the system a bit, say at $t=0.2$, in Fig. 4,
one can notice that in the red part, the PR value is increasing, for higher
values of $g_{1}$ and $g_{2}$. In fact this portion in the central region
increases with time also. In Fig. 4 we see that in the red part the value of
PR is significantly higher than in the blue part. When the dynamics is carried
out over a very long time [say $t=1000$ as in Fig. 4], we get a phase diagram
which exactly looks like the non-ergodic (blue color) to ergodic (red color)
phase diagram [see Fig. 4 and 3]. Figure 4 suggests that, if the system is
initially in the non-ergodic phase, and we evolve it for a long time, the
system will stay in the non-ergodic phase, i.e., the participation ratio is
relatively low no matter how long the time is. This is in contrast to the
ergodic phase where the PR value is higher for long times even if the inital
PR is low.

##  IV Thermal phase transition (TPT)

Figure 5: Mutual information of two spins $I_{12}$ as a function of
temperature. Inset of the figure shows the numerical differentiation of MI
with respect to temperature $\frac{dI_{12}}{dT}$, at $g_{1}=1.0$, $g_{2}=0.5$.
The vertical line represents the theoretical value of the critical
temperature, $T_{c}\approx 0.46$. Here $\omega=\omega_{0}=1$,
$n_{\text{max}}=40$.

Figure 6: Mutual information between two spins as a function of $g_{1}$ and
$g_{2}$ at different temperature: (a) $T=0$ (QPT), (b) $T=0.1$, (c) $T=1.0$,
(d) $T=1.5$. As temperature increases, the region corresponding to the normal
phase (denoted by black color) also increases. The red dashed line denotes the
value of $g_{1}+g_{2}$ for a fixed temperature following Eq. 19. Zero
temperature case denotes QPT along the line $g_{1}+g_{2}=1$, for $T=0.1$, $1$
and $1.5$, $g_{1}+g_{2}=1.01$, $2.16$, $2.63$ respectively. The parameters
are: $\omega=\omega_{0}=1$, $N=6$. The bosonic cut-off is taken to be
$n_{\text{max}}=40$.

Another type of phase transition exhibited by the ADM is the finite
temperature phase transition [26]. As with the Dicke model, it is known that a
finite critical temperature can destroy the super-radiant phase and a
transition back to the normal phase is obtained going beyond the critical
temperature [26]. We can derive an exact analytical expression for the
transition temperature as a function of the parameters $g_{1}$ and $g_{2}$. We
start by rewriting the Hamiltonian (in units of $\omega$) as:

$\displaystyle\tilde{\mathcal{H}}=\frac{\mathcal{H}}{\omega}=a^{\dagger}a+\sum_{j=1}^{N}\frac{\epsilon}{2}\sigma_{j}^{z}+\frac{\lambda_{1}}{2\sqrt{N}}\sum_{j=1}^{N}(a\sigma_{j}^{+}+a^{\dagger}\sigma_{j}^{-})$
$\displaystyle+\frac{\lambda_{2}}{2\sqrt{N}}\sum_{j=1}^{N}(a^{\dagger}\sigma_{j}^{+}+a\sigma_{j}^{-})$
(13)

where, $\epsilon=\frac{\omega_{0}}{\omega}$,
$\lambda_{1}=\frac{g_{1}}{\omega}$, $\lambda_{2}=\frac{g_{2}}{\omega}$. The
partition function for the full ADM is given by:

$\displaystyle Z(N,T)=\sum_{s_{1},...,s_{N}=\pm
1}\int\frac{d^{2}\alpha}{\pi}\langle
s_{1},...,s_{N}|\langle\alpha|e^{-\beta\tilde{\mathcal{H}}}|\alpha\rangle|s_{1},...,s_{N}\rangle.$
(14)

The expectation value of the Hamiltonian with respect to the bosonic modes is:

$\displaystyle\langle\alpha|\tilde{\mathcal{H}}|\alpha\rangle=\alpha^{*}\alpha+\sum_{j=1}^{N}\Big{[}\frac{\epsilon}{2}\sigma_{j}^{z}+\frac{\lambda_{1}}{2\sqrt{N}}(\alpha\sigma_{j}^{+}+\alpha^{*}\sigma_{j}^{-})$
$\displaystyle+\frac{\lambda_{2}}{2\sqrt{N}}(\alpha^{*}\sigma_{j}^{+}+\alpha\sigma_{j}^{-})\Big{]}.$
(15)

Defining

$\displaystyle
h_{j}=\frac{\epsilon}{2}\sigma_{j}^{z}+\frac{\lambda_{1}}{2\sqrt{N}}(\alpha\sigma_{j}^{+}+\alpha^{*}\sigma_{j}^{-})+\frac{\lambda_{2}}{2\sqrt{N}}(\alpha^{*}\sigma_{j}^{+}+\alpha\sigma_{j}^{-})$
(16)

the expectation value with respect to the spins becomes a product of single-
spin expectation values:

$\displaystyle\langle
s_{1}...s_{N}|\langle\alpha|e^{-\beta\tilde{\mathcal{H}}}|\alpha\rangle|s_{1}...s_{N}\rangle$
$\displaystyle=$ $\displaystyle e^{-\beta|\alpha|^{2}}\Pi_{j=1}^{N}\langle
s_{j}|e^{-\beta h_{j}}|s_{j}\rangle.$ (17)

Thus the computation of the partition function reduces to the evaluation of a
double integral:

$\displaystyle
Z(N,T)=\int\frac{d^{2}\alpha}{\pi}e^{-\beta|\alpha|^{2}}\Big{[}\text{Tr}e^{-\beta
h}\Big{]}^{N}$
$\displaystyle=\int\frac{d^{2}\alpha}{\pi}e^{-\beta|\alpha|^{2}}\Big{(}2\cosh\Big{[}\frac{\beta\epsilon}{2}\Big{[}1+\frac{4(\lambda_{1}+\lambda_{2})^{2}\alpha^{2}}{\epsilon^{2}N}\Big{]}^{1/2}\Big{]}\Big{)}^{N}$
(18)

which in the thermodynamic limit ($N\to\infty$) can be solved using the method
of steepest descent, within the super-radiant phase. Tracking the point at
which the method breaks down (see the appendix), we have an exact expression
for the transition temperature:

$T_{c}=\frac{1}{\beta_{c}}=\Big{(}\frac{\omega_{0}}{2\omega}\Big{)}\frac{1}{\tanh^{-1}\Big{(}\frac{\omega\omega_{0}}{(g_{1}+g_{2})^{2}}\Big{)}}.$
(19)

Thus, in the super-radiant phase ($g_{1}+g_{2}>1$ at $T=0$), raising the
temperature to a value larger than the critical temperature ($T_{c}$), causes
the system to go back to the normal phase. We now provide numerical evidence
of this phase transition with the help of mutual information between two
spins. While the entanglement entropy is a good measure to capture a QPT, it
is unsuitable for a TPT since mixed states are involved. Mutual information
[45, 46, 47, 48], between two spins is a good measure to capture the TPT, as
we showned in an earlier work for the Dicke model [10]. We show here that the
usefulness of mutual information as a marker of the thermal phase transition
extends to the ADM.

When the overall state is mixed, the correlations between two subsystems can
be quantified with the help of the mutual information defined as

$I_{12}=S_{1}+S_{2}-S_{12},$ (20)

where $S_{1,2}=-\text{Tr}[\rho_{1,2}\ln(\rho_{1,2})]$ ,
$S_{12}=-\text{Tr}[\rho_{12}\ln(\rho_{12})]$. Here $\rho_{1}$, $\rho_{2}$ are
the reduced density matrices for the two subsystems, $S_{1}$, $S_{2}$ are the
corresponding von Neumann entropies, $\rho_{12}$ is the density matrix of the
overall system, and $S_{12}$ is the corresponding entropy. When the overall
state is in a pure state, $S_{12}=0$ and the mutual information become twice
the entanglement entropy since $S_{1}=S_{2}$. For our model we study the
mutual information between any two spins - (due to the symmetry of the system
Hamiltonian, it does not matter which spin pair is chosen). Here we use the
spin product space hence we have to diagonalize the system Hamiltonian with
dimension $(n_{\text{max}}+1)2^{N}$. To calculate $I_{\text{12}}$, we have to
take a partial trace of the total density matrix over the bosonic part first
and then over the $N-2$ atoms.

Figure 5 shows the mutual information between two spins. For the parameters
$g_{1}=1.0,g_{2}=0.5$, the system is in the super-radiant phase at zero
temperature where the value of mutual information takes significantly large
values (Fig. 5). On increasing the temperature, we see that the mutual
information starts to decrease, signifying a change in the direction of the
normal phase. To check that mutual information does indeed capture the exact
transition between the super-radiant-to-normal phase, we plot the derivative
of the mutual information with respect to the temperature:
$\frac{d{I_{12}}}{dT}$ (inset of the Fig. 5). We observe that the temperature
at which the derivative is minimum corresponds to the transition temperature
of the TPT from SP to NP for $g_{1}+g_{2}>1$. For this particular choice of
$g_{1}$ and $g_{2}$, the critical temperature is $T_{c}=0.46$, which we denote
by the vertical straight line in the inset figure.

In Fig. 6 we showned $I_{\text{12}}$ as a function of $g_{1}$ and $g_{2}$ at
different temperatures: $T=0$, $0.1$, $1.0$, and $1.5$. Figure 6 ($T=0$) shows
a clear QPT from NP (black color) to SP (white color) along the line
$g_{1}+g_{2}=1$. We emphasize that at zero temperature, the ground state is a
pure state, and so the mutual information is really the same as twice the
entanglement entropy. $I_{12}$ is close to zero in the normal phase and close
to unity in the super-radiant phase. Another way of saying this is that the
total correlation between two spins is almost zero in the normal phase whereas
it is maximum in the super-radiant phase. Hence the QPT in the anisotropic
Dicke model is similar to the isotropic Dicke model (for which $g_{1}=g_{2}$)
but here an additional parameter is introduced. On the other hand, one can
notice that as temperature increases [in Fig. 6 ($T=0.1$), 6 ($T=1$) and 6
($T=1.5$)] the region corresponding to the NP (black portion) also expands -
the phase boundary is highlighted by a red dashed line
($g_{1}+g_{2}=\text{constant}$) which depends on the temperature. Using Eq. 19
we see that for $T=0.1$, $1$, and $1.5$, the phase boundary is given by
$g_{1}+g_{2}=1.01$, $g_{1}+g_{2}=2.16$, and $g_{1}+g_{2}=2.64$, respectively;
beyond the phase boundary, the system is in the SP where the correlation
between two spins is close to unity. These figures indicate that the mutual
information between spins is an excellent measure of the thermal phase
transition of the system.

Figure 7: (a), (f), (k), (p), and (u): $\chi_{\text{lower}}$ [given in Eqn.
9], (b), (g), (l), (q), and (v): $\chi_{\text{upper}}$ (given in Eqn. 10),
(c), (h), (m), (r), and (w): average consecutive level spacing ratio $\langle
r\rangle$; (d), (i), (n), (s), and (x): participation ratio ($PR$) and (e),
(j), (o), (t), and (y): the multifractal dimension $D_{1}$ for the middle
excited state of the system for the resonance condition
($\omega=\omega_{0}=1$) considering $20$ spins and truncating the bosonic mode
at gradually increasing values: $n_{\text{max}}=200,\ 300,\ 400,\ 500,\ 600$
(panels $1$ to $5$, respectively).

##  V Summary and conclusions

We first discuss the ground state phase transition from a normal phase to the
super-radiant phase showing that a critical line $g_{1}+g_{2}=1$ separates the
two phases. We see that the ground state energy density is almost constant in
the NP whereas in the SP we find a broad range of energy densities, which are
lower than that in the normal phase. The ground state number operator is
almost zero in the NP while it is non-zero in the SP, indicating macroscopic
excitations in the bosonic mode. By studying the participation ratio, we
conclude that the ground state of the system exhibits multifractal features in
the super-radiant phase with the participation ratio scaling as
$PR_{\text{gs}}\propto\sqrt{N_{D}}$.

Next, we explore the excited state features and find that the ADM also
exhibits the excited state phase transition both for the resonant
($\omega=\omega_{0}$) and the off-resonant ($\omega\neq\omega_{0}$) cases. The
ESQPT is nicely captured by the von Neumann entanglement entropy (between
spins and bosons) as a function of eigenstate energies. We observe that for
the ADM there exist two cut-off energies separating the different phases: a
lower cut-off energy (corresponding to the ground state energy along the line,
$g_{1}+g_{2}=1$) and an upper cut-off energy (corresponding to the maximum
energy at $g_{1}=g_{2}=0$ for finite $n_{\text{max}}$). Between these two
cutoff energies, we find that the level statistics exhibits either Poisson
statistics or Wigner-Dyson statistics depending on the values of the coupling
parameters $g_{1}$ and $g_{2}$ suggesting a non-ergodic to ergodic phase
transition. A study of the consecutive level spacing ratio of the system for
the middle energy band (energy band between the lower and the upper cut-off
energies) supports these findings. It is convenient to introduce two new
quantities (having the dimensions of energy) that correspond to the lower and
upper cut-off energies of the spectrum in the super-radiant phase. The above
energies are obtained by using the jumps in the VNEE as weights. These
characteristic energies, which are a measure of the energies at which the
corresponding von Neumman entropies of the eigenstates begin and end their
plateau-like behavior, signal the ESQPT of this model. Remarkably when a phase
diagram is obtained using these characteristic energies, we get a picture that
looks very similar to the phase diagram obtained using level spacing ratios.
Thus from our study, we conclude that the ESQPT and ENET are intimately
related to each other for the anisotropic Dicke model ($g_{1}\neq g_{2}$). We
checked that this connection is robust both for the resonant
($\omega=\omega_{0}$) and off-resonant ($\omega\neq\omega_{0}$) cases for the
generic ADM; the diagonal direction ($g_{1}=g_{2}$) corresponds to the Dicke
model exhibits special behavior [6]. From a study of eigenstate properties
with the aid of participation ratio and multifractal dimension $D_{1}$, we see
that the normal-to-super-radiant phase transition corresponds to the ground
state undergoing a localized-to-multifractal transition. On the other hand the
ergodic-to-nonergodic transition corresponds to the middle excited state
undergoing a delocalized-to-multifractal transition. The correspondence
between the multifractal nature of the middle excited state and the nonergodic
phase is also captured dynamically when we study participation ratio in a
quench dynamical protocol.

Finally, the ADM also exhibits yet another phase transition namely the
temperature dependent phase transition. For $g_{1}+g_{2}>1$, there exists a
critical temperature, $T_{c}$ above which the super-radiant phase disappears
and the system comes back to the normal phase. Following the old work of Hioe
[26] we write down an analytical expression for $T_{c}$ as a function of
$g_{1}$, $g_{2}$. We show that the mutual information between two spins as a
function of temperature proves handy to obtain an independent characterization
of the thermal phase transition. A study of the mutual information suggests
that for $g_{1}+g_{2}<1$, the system lies in the normal phase for all
temperatures with a relatively lower value of mutual information, but for
$g_{1}+g_{2}>1$ there exists a $T_{c}$ such that when $T<T_{c}$ the system
lies in the SP with a relatively higher value of mutual information and for
$T>T_{c}$ the system goes back to NP showing a temperature dependent phase
transition.

## Acknowledgments

We are grateful to the High Performance Computing(HPC) facility at IISER
Bhopal, where large-scale calculations in this project were run. P.D. is
grateful to IISERB for the PhD fellowship. A.S acknowledges financial support
from SERB via Grant No: CRG/2019/003447, and from DST via the DST-INSPIRE
Faculty Award No. DST/INSPIRE/04/2014/002461.

##  Appendix A Expression for the critical temperature for the TPT

In this appendix we include a brief derivation of the expression for the
critical temperature for TPT similar to the previous work [26]. We begin by
recalling the expression of the partition function in the form of a double
integral [Eqn. 18]:

$\displaystyle
Z(N,T)=\int\frac{d^{2}\alpha}{\pi}e^{-\beta|\alpha|^{2}}\Big{[}\text{Tr}e^{-\beta
h}\Big{]}^{N}$
$\displaystyle=\int\frac{d^{2}\alpha}{\pi}e^{-\beta|\alpha|^{2}}\Big{(}2\cosh\Big{[}\frac{\beta\epsilon}{2}\Big{[}1+\frac{4(\lambda_{1}+\lambda_{2})^{2}\alpha^{2}}{\epsilon^{2}N}\Big{]}^{1/2}\Big{]}\Big{)}^{N}.$
(21)

Rewriting the double integral using polar coordinates and defining
$y=\frac{r^{2}}{N}$ and

$\displaystyle\phi(y)=-\beta
y+\ln\Big{(}2\cosh\Big{[}\frac{\beta\epsilon}{2}\Big{[}1+\frac{4(\lambda_{1}+\lambda_{2})^{2}y}{\epsilon^{2}}\Big{]}^{1/2}\Big{]}\Big{)}$
(22)

we can write:

$\displaystyle Z(N,T)=N\int_{0}^{\infty}dy\exp\Big{(}N\phi(y)\Big{)}.$ (23)

Since we are interested in the thermodynamic limit where $N\to\infty$, we can
invoke Laplace’s method [67] to evaluate the integral as:

$\displaystyle Z(N,T)=N\frac{C}{\sqrt{N}}\max_{0\leq
y\leq\infty}\exp\Big{(}N\Big{[}\phi(y)\Big{]}\Big{)}$ (24)

where $C$ is some constant. To find the maximum of the function $\phi(y)$, we
compute its derivative:

$\displaystyle\phi^{\prime}=-\beta+\frac{\beta(\lambda_{1}+\lambda_{2})^{2}}{\epsilon}\frac{1}{\eta}\tanh\Big{(}\frac{\beta\epsilon\eta}{2}\Big{)}$
(25)

where

$\displaystyle\eta=\Big{[}1+\frac{4(\lambda_{1}+\lambda_{2})^{2}y}{\epsilon^{2}}\Big{]}^{1/2}.$
(26)

Putting

$\displaystyle\phi^{\prime}=0,$ (27)

we have

$\displaystyle\frac{\epsilon\eta}{(\lambda_{1}+\lambda_{2})^{2}}=\tanh\Big{(}\frac{\beta\epsilon\eta}{2}\Big{)}.$
(28)

The hyperbolic tangent funtion is a monotonically increasing function and is
bounded above by unity. Since $\eta\geq 1$ by definition [Eqn. 26], if
$(\lambda_{1}+\lambda_{2})^{2}<\epsilon$, there is no solution for Eq. 28. On
the other hand, for $(\lambda_{1}+\lambda_{2})^{2}>\epsilon$, the solution
depends on the value of $\beta$. The critical value of the inverse temperature
$\beta_{c}$ can be computed by putting $\eta=1$ and is given by:

$\displaystyle\beta_{c}=\frac{2}{\epsilon}\tanh^{-1}\Big{(}\frac{\epsilon}{(\lambda_{1}+\lambda_{2})^{2}}\Big{)}.$
(29)

Substituting $\epsilon=\frac{\omega_{0}}{\omega}$ and
$\lambda_{1}=\frac{g_{1}}{\omega}$, $\lambda_{2}=\frac{g_{2}}{\omega}$, we
have an exact expression for the transition temperature:

$T_{c}=\Big{(}\frac{\omega_{0}}{2\omega}\Big{)}\frac{1}{\tanh^{-1}\Big{(}\frac{\omega\omega_{0}}{(g_{1}+g_{2})^{2}}\Big{)}}.$
(30)

##  Appendix B Robustness of our results against increasing bosonic cut-off

In Fig. 7 we show the data corresponding to the ESQPT and ENET for $20$ spins
and considering gradually increasing values of the bosonic truncation number:
$n_{\text{max}}=200,\ 300,\ 400,\ 500,\ 600$ (panel $1$ to $5$ respectively)
to check the robustness of our results. In each panel of this figure, we show
$\chi_{\text{lower}}$ (the lower cut-off energy which is scaled by the minimum
energy in the nornal phase), $\chi_{\text{upper}}$ (upper cut-off energy which
is scaled by the maximum energy at $g_{1}=g_{2}=0$), consecutive level spacing
ratio $\langle r\rangle$, participation ratio ($PR$) and the multifractal
dimension ($D_{1}$) of the middle excited state respectively as a function of
the coupling parameters $g_{1}$ and $g_{2}$. We notice that the results are
qualitatively unchanged against increasing values of the bosonic truncation
number $n_{\text{max}}$ and are also converging.

## References

  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. 93, 99 (1954).
  * Lambert  _et al._ [2004] N. Lambert, C. Emary, and T. Brandes, Entanglement and the phase transition in single-mode superradiance, Phys. Rev. Lett. 92, 073602 (2004).
  * Emary and Brandes [2003a] C. Emary and T. Brandes, Chaos and the quantum phase transition in the dicke model, Phys. Rev. E 67, 066203 (2003a).
  * Emary and Brandes [2003b] C. Emary and T. Brandes, Quantum chaos triggered by precursors of a quantum phase transition: The dicke model, Phys. Rev. Lett. 90, 044101 (2003b).
  * Kadantseva  _et al._ [1990] E. Kadantseva, W. Chmielowski, and A. Shumovsky, Superradiance in the dicke model, in  _Nonlinear Optics in Solids_ (Springer, 1990) pp. 37–41.
  * Chávez-Carlos  _et al._ [2016] J. Chávez-Carlos, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, Classical chaos in atom-field systems, Phys. Rev. E 94, 022209 (2016).
  * Kirton and Keeling [2018] P. Kirton and J. Keeling, Superradiant and lasing states in driven-dissipative dicke models, New Journal of Physics 20, 015009 (2018).
  * Kirton  _et al._ [2019] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Introduction to the dicke model: from equilibrium to nonequilibrium, and vice versa, Advanced Quantum Technologies 2, 1800043 (2019).
  * Prasad  _et al._ [2022] M. Prasad, H. K. Yadalam, C. Aron, and M. Kulkarni, Dissipative quantum dynamics, phase transitions, and non-hermitian random matrices, Phys. Rev. A 105, L050201 (2022).
  * Das and Sharma [2022] P. Das and A. Sharma, Revisiting the phase transitions of the dicke model, Phys. Rev. A 105, 033716 (2022).
  * Yamamoto  _et al._ [2020] D. Yamamoto, C. Suzuki, G. Marmorini, S. Okazaki, and N. Furukawa, Quantum and thermal phase transitions of the triangular su (3) heisenberg model under magnetic fields, Phys. Rev. Lett. 125, 057204 (2020).
  * Brandes [2013] T. Brandes, Excited-state quantum phase transitions in dicke superradiance models, Phys. Rev. E 88, 032133 (2013).
  * Bastarrachea-Magnani  _et al._ [2014] M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, Comparative quantum and semiclassical analysis of atom-field systems. i. density of states and excited-state quantum phase transitions, Phys. Rev. A 89, 032101 (2014).
  * Stránskỳ and Cejnar [2016] P. Stránskỳ and P. Cejnar, Classification of excited-state quantum phase transitions for arbitrary number of degrees of freedom, Physics Letters A 380, 2637 (2016).
  * Stránskỳ  _et al._ [2017] P. Stránskỳ, M. Kloc, and P. Cejnar, Excited-state quantum phase transitions and their manifestations in an extended dicke model, in  _AIP Conference Proceedings_, Vol. 1912 (AIP Publishing LLC, 2017) p. 020018.
  * Cejnar  _et al._ [2021] P. Cejnar, P. Stránskỳ, M. Macek, and M. Kloc, Excited-state quantum phase transitions, Journal of Physics A: Mathematical and Theoretical 54, 133001 (2021).
  * Carmichael  _et al._ [1973] H. Carmichael, C. Gardiner, and D. Walls, Higher order corrections to the dicke superradiant phase transition, Physics Letters A 46, 47 (1973).
  * Wang and Hioe [1973] Y. K. Wang and F. T. Hioe, Phase transition in the dicke model of superradiance, Phys. Rev. A 7, 831 (1973).
  * Duncan [1974] G. C. Duncan, Effect of antiresonant atom-field interactions on phase transitions in the dicke model, Phys. Rev. A 9, 418 (1974).
  * Bastarrachea-Magnani  _et al._ [2017] M. A. Bastarrachea-Magnani, A. Relaño, S. Lerma-Hernández, B. L. del Carpio, J. Chávez-Carlos, and J. G. Hirsch, Adiabatic invariants for the regular region of the dicke model, Journal of Physics A: Mathematical and Theoretical 50, 144002 (2017).
  * Rubeni  _et al._ [2012] D. Rubeni, A. Foerster, E. Mattei, and I. Roditi, Quantum phase transitions in bose–einstein condensates from a bethe ansatz perspective, Nuclear Physics B 856, 698 (2012).
  * Greentree  _et al._ [2006] A. D. Greentree, C. Tahan, J. H. Cole, and L. C. Hollenberg, Quantum phase transitions of light, Nature Physics 2, 856 (2006).
  * Stranius  _et al._ [2018] K. Stranius, M. Hertzog, and K. Börjesson, Selective manipulation of electronically excited states through strong light–matter interactions, Nature communications 9, 1 (2018).
  * Zhou  _et al._ [2019] J. Zhou, B. Huang, Z. Yan, and J.-C. G. Bünzli, Emerging role of machine learning in light-matter interaction, Light: Science & Applications 8, 1 (2019).
  * Mueller  _et al._ [2020] N. S. Mueller, Y. Okamura, B. G. Vieira, S. Juergensen, H. Lange, E. B. Barros, F. Schulz, and S. Reich, Deep strong light–matter coupling in plasmonic nanoparticle crystals, Nature 583, 780 (2020).
  * Hioe [1973] F. T. Hioe, Phase transitions in some generalized dicke models of superradiance, Phys. Rev. A 8, 1440 (1973).
  * de Aguiar  _et al._ [1991] M. de Aguiar, K. Furuya, and M. C. Nemes, The classical analogue of the super-radiant phase transition in the dicke model, Quantum Optics: Journal of the European Optical Society Part B 3, 305 (1991).
  * De Aguiar  _et al._ [1992] M. De Aguiar, K. Furuya, C. Lewenkopf, and M. Nemes, Chaos in a spin-boson system: classical analysis, Annals of Physics 216, 291 (1992).
  * Furuya  _et al._ [1992] K. Furuya, M. de Aguiar, C. Lewenkopf, and M. Nemes, Husimi distributions of a spin-boson system and the signatures of its classical dynamics, Annals of Physics 216, 313 (1992).
  * Bastarrachea-Magnani  _et al._ [2016] M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, Thermal and quantum phase transitions in atom-field systems: a microcanonical analysis, Journal of Statistical Mechanics: Theory and Experiment 2016, 093105 (2016).
  * Buijsman  _et al._ [2017] W. Buijsman, V. Gritsev, and R. Sprik, Nonergodicity in the anisotropic dicke model, Phys. Rev. Lett. 118, 080601 (2017).
  * Kloc  _et al._ [2017] M. Kloc, P. Stránskỳ, and P. Cejnar, Quantum phases and entanglement properties of an extended dicke model, Annals of Physics 382, 85 (2017).
  * Aedo and Lamata [2018] I. Aedo and L. Lamata, Analog quantum simulation of generalized dicke models in trapped ions, Phys. Rev. A 97, 042317 (2018).
  * Shapiro  _et al._ [2020] D. S. Shapiro, W. V. Pogosov, and Y. E. Lozovik, Universal fluctuations and squeezing in a generalized dicke model near the superradiant phase transition, Phys. Rev. A 102, 023703 (2020).
  * Hu and Wan [2021] J. Hu and S. Wan, Out-of-time-ordered correlation in anisotropic dicke model, Communications in Theoretical Physics 73, 125703 (2021).
  * Baumann  _et al._ [2011] K. Baumann, R. Mottl, F. Brennecke, and T. Esslinger, Exploring symmetry breaking at the dicke quantum phase transition, Phys. Rev. Lett. 107, 140402 (2011).
  * Bhattacharya  _et al._ [2014] U. Bhattacharya, S. Dasgupta, and A. Dutta, Exploring chaos in the dicke model using ground-state fidelity and loschmidt echo, Phys. Rev. E 90, 022920 (2014).
  * Pérez-Fernández and Relaño [2017] P. Pérez-Fernández and A. Relaño, From thermal to excited-state quantum phase transition: The dicke model, Phys. Rev. E 96, 012121 (2017).
  * Ashhab  _et al._ [2019] S. Ashhab, Y. Matsuzaki, K. Kakuyanagi, S. Saito, F. Yoshihara, T. Fuse, and K. Semba, Spectrum of the dicke model in a superconducting qubit-oscillator system, Phys. Rev. A 99, 063822 (2019).
  * Lewis-Swan  _et al._ [2019] R. Lewis-Swan, A. Safavi-Naini, J. J. Bollinger, and A. M. Rey, Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the dicke model, Nature communications 10, 1 (2019).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * Roy and Sharma [2018] N. Roy and A. Sharma, Entanglement contour perspective for “strong area-law violation” in a disordered long-range hopping model, Phys. Rev. B 97, 125116 (2018).
  * Nehra  _et al._ [2018] R. Nehra, D. S. Bhakuni, S. Gangadharaiah, and A. Sharma, Many-body entanglement in a topological chiral ladder, Phys. Rev. B 98, 045120 (2018).
  * DiVincenzo  _et al._ [2004] D. P. DiVincenzo, M. Horodecki, D. W. Leung, J. A. Smolin, and B. M. Terhal, Locking classical correlations in quantum states, Phys. Rev. Lett. 92, 067902 (2004).
  * Lu  _et al._ [2011] X.-M. Lu, J. Ma, Z. Xi, and X. Wang, Optimal measurements to access classical correlations of two-qubit states, Phys. Rev. A 83, 012327 (2011).
  * Henderson and Vedral [2001] L. Henderson and V. Vedral, Classical, quantum and total correlations, Journal of physics A: mathematical and general 34, 6899 (2001).
  * Adesso and Datta [2010] G. Adesso and A. Datta, Quantum versus classical correlations in gaussian states, Phys. Rev. Lett. 105, 030501 (2010).
  * Roósz  _et al._ [2014] G. m. H. Roósz, U. Divakaran, H. Rieger, and F. Iglói, Nonequilibrium quantum relaxation across a localization-delocalization transition, Phys. Rev. B 90, 184202 (2014).
  * Sirker  _et al._ [2014] J. Sirker, M. Maiti, N. Konstantinidis, and N. Sedlmayr, Boundary fidelity and entanglement in the symmetry protected topological phase of the ssh model, Journal of Statistical Mechanics: Theory and Experiment 2014, P10032 (2014).
  * Cho and Kim [2017] J. Cho and K. W. Kim, Quantum phase transition and entanglement in topological quantum wires, Scientific reports 7, 1 (2017).
  * Ma  _et al._ [2013] Z. Ma, Z. Chen, and F. F. Fanchini, Multipartite quantum correlations in open quantum systems, New Journal of Physics 15, 043023 (2013).
  * Sharma and Rabani [2015] A. Sharma and E. Rabani, Landauer current and mutual information, Phys. Rev. B 91, 085121 (2015).
  * Sable  _et al._ [2018] H. S. Sable, D. S. Bhakuni, and A. Sharma, Landauer current and mutual information in a bosonic quantum dot, Journal of Physics: Conference Series 964, 012007 (2018).
  * Hepp and Lieb [1973] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model, Annals of Physics 76, 360 (1973).
  * Chang [1995] Z. Chang, Generalized holstein-primakoff transformation, spin-charge separation and quantum nonlinear $\sigma$ model in a 2d lattice, Physics Letters A 204, 405 (1995).
  * Ressayre and Tallet [1975] E. Ressayre and A. Tallet, Holstein-primakoff transformation for the study of cooperative emission of radiation, Phys. Rev. A 11, 981 (1975).
  * Edwards and Thouless [1972] J. T. Edwards and D. J. Thouless, Numerical studies of localization in disordered systems, Journal of Physics C: Solid State Physics 5, 807 (1972).
  * Roy  _et al._ [2020] N. Roy, A. Ramachandran, and A. Sharma, Interplay of disorder and interactions in a flat-band supporting diamond chain, Phys. Rev. Research 2, 043395 (2020).
  * Pérez-Fernández  _et al._ [2011] P. Pérez-Fernández, A. Relaño, J. M. Arias, P. Cejnar, J. Dukelsky, and J. E. García-Ramos, Excited-state phase transition and onset of chaos in quantum optical models, Phys. Rev. E 83, 046208 (2011).
  * Hill and Wootters [1997] S. Hill and W. K. Wootters, Entanglement of a pair of quantum bits, Phys. Rev. Lett. 78, 5022 (1997).
  * Wootters [1998] W. K. Wootters, Entanglement of formation of an arbitrary state of two qubits, Phys. Rev. Lett. 80, 2245 (1998).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Kamenev and Mézard [1999] A. Kamenev and M. Mézard, Wigner-dyson statistics from the replica method, Journal of Physics A: Mathematical and General 32, 4373 (1999).
  * Macé  _et al._ [2019] N. Macé, F. Alet, and N. Laflorencie, Multifractal scalings across the many-body localization transition, Phys. Rev. Lett. 123, 180601 (2019).
  * Lindinger  _et al._ [2019] J. Lindinger, A. Buchleitner, and A. Rodríguez, Many-body multifractality throughout bosonic superfluid and mott insulator phases, Phys. Rev. Lett. 122, 106603 (2019).
  * Jeffreys  _et al._ [1999] H. Jeffreys, B. Jeffreys, and B. Swirles,  _Methods of mathematical physics_ (Cambridge university press, 1999).