# Quantum chaos in the Dicke model and its variants

Devvrat Tiwari<EMAIL_ADDRESS>Subhashish Banerjee<EMAIL_ADDRESS>Indian Institute of Technology Jodhpur-342030, India

###### Abstract

Recently, the out-of-time-ordered correlator (OTOC) has gained much attention
as an indicator of quantum chaos. In the semi-classical limit, its exponential
growth rate resembles the classical Lyapunov exponent. The quantum-classical
correspondence has been supported for the one-body chaotic systems as well as
realistic systems with interactions, as in the Dicke model, a model of multi-
two-level atoms and cavity field interactions. To this end, we calculate the
OTOC for different variations of the Dicke model in an open quantum system
setting. The connection between the superradiant phase transition of the Dicke
model and the OTOC is studied. Further, we establish a relation between the
OTOC and the second-order coherence function. This becomes important for the
experimental studies of the OTOC and quantum chaos in the models of quantum
optics.

Quantum chaos, out-of-time-ordered-correlator, Dicke model, open quantum
systems

##  I Introduction

The study of quantum chaos seeks to find a relationship between the classical
and quantum dynamics of chaotic systems [1, 2, 3, 4, 5, 6, 7, 8]. In the limit
of large quantum numbers, classical mechanics is shown by the correspondence
principle as the classical limit of quantum mechanics. Therefore, classical
chaos should also have an underlying quantum mechanism. Quantum chaos seeks
answers to this question of the relationship between quantum mechanics and
classical chaos [5, 3]. Random matrix theory ideas [9, 10, 11], formulated by
Wigner [12] and Dyson [13] to understand the spectra of complex atomic nuclei,
were used traditionally to probe quantum systems with chaotic classical
counterparts [14, 3, 15, 16]. Various strategies have been created to bridge
the gap between classical and quantum depictions of chaos by using features
that are similar in both realms. This includes techniques, such as the level
statistics of the quantum states [17, 18], the structure of the eigenstates
[19, 20], the exponential increase of complexity [21, 22], and the exponential
decay of the Loschmidt echo (a measure of the revival that occurs when an
imperfect time-reversal procedure is applied to a complex quantum system) [23,
24, 25].

An idea of quantum information scrambling (a process of spreading and
effective loss of quantum information in a quantum many-body system) by black
holes was brought out [26, 27], which was popularized and given a great fillip
[28]. This scrambling can be measured using out-of-time-ordered correlators
(OTOC), which was conjectured to have a bound on its growth rate [29, 30, 31].
The OTOC, first established in the context of superconductivity [32], is now
presented as a measure of quantum chaos, with its growth rate attributed to
the classical Lyapunov exponent [33, 34, 35, 36, 37, 38]. Recently, OTOC has
started receiving considerable attention because of its appearance as a
valuable tool for generalizing quantum chaos studies beyond the time-
independent, one-body case. Not only that, it has been extensively studied in
a wide variety of systems, including quantum thermalization [39, 40], Sachdev-
Ye-Kitaev model [41, 42, 30, 43], quantum field theories [44, 45, 46, 47],
random unitary models [48, 49, 50, 51], and spin chains [52, 53, 54, 55].
Moreover, the OTOC has also been measured experimentally using techniques of
nuclear magnetic resonance [56, 57, 58, 59].

Atom-cavity field interactions are critical to science and technology and
provide a foundation for the field of quantum optics [60, 61, 62, 63, 64]. The
first quantum mechanical model describing the interaction between a single
atom and a cavity of electromagnetic field was given by Jaynes and Cummings
[65, 66]. Here, we discuss a class of quantum optical models, which are multi-
atom generalizations of the Jaynes-Cummings model. This includes variations of
the model proposed by Dicke in 1954 [67] to describe super-radiance, a field-
mediated collective emission of an ensemble of excited atoms or ions [68, 69,
70, 71, 72, 73]. The superradiant phase transition of the Dicke model has been
studied in recent experiments [74, 75, 76, 77, 78]. The presence of quantum
chaos and quantum phase transition in the Dicke model has long been an area of
interest [79, 80, 81, 82, 83, 84, 85, 86]. A demonstration of quantum-
classical correspondence, using OTOC for the Dicke model, was made in the
unitary case [87]. Further, the Dicke model under dissipation was studied to
explore the connections between the dissipative phase transition and non-
Hermitian random matrix theory in [88]. Here, we endeavor to explore the
chaotic behavior of the variations of the Dicke model in an open quantum
system setting with dissipation.

A realistic quantum system interacts with its surroundings. These interactions
significantly change the dynamics of the quantum system. There is a loss of
information from the system to the environment involved. A framework to study
such systems is provided by the theory of open quantum systems (OQS) [89, 90,
91]. The ideas of open quantum systems find use in a wide variety of
applications [92, 93, 94, 95, 96, 97, 98, 99, 97, 98, 99, 100, 101, 102, 103,
104, 105, 106, 107, 108]. Exploring quantum chaos and OTOC in open quantum
systems is an ongoing area of research. Several techniques have been developed
to address this issue [109, 110, 111, 112, 113, 114]. To this end, we use the
adjoint form of the GKSL master equation [89, 90] to study the operator
evolution in time, which is then used to calculate the OTOC. To calculate the
OTOC in Markovian open quantum systems, an extension of the quantum regression
theorem was derived in [115], where an indirect connection between the second-
order coherence $g^{(2)}(t)$ function and OTOC was given. To this effect, we
try to establish a direct relationship between the characteristics of quantum
optics ($g^{(2)}(t)$ function) and the chaotic behavior of the quantum optical
models.

The paper is structured as follows: in Sec. II, we discuss the preliminaries
of the quantum optical models, such as their quantum phase transitions and the
out-of-time-ordered correlator (OTOC) and its various forms. Section III
discusses the behavior of the OTOC in the models discussed here with the
calculation of the Lyapunov exponents and the effect of quantum phase
transition on the chaotic behavior of these models. The behavior of OTOC in
short, intermediate, and long-time is also discussed. In Sec. V, we highlight
the connection between the $g^{(2)}(t)$ function and the OTOC and discuss the
characteristics of light in conjunction with quantum chaos. This is followed
by the Conclusions.

##  II Preliminaries

In this section, we briefly discuss the many-body quantum optical models,
particularly the models where $N$-qubits interact with QED-cavity. In order to
understand the chaotic properties of these models, we will make use of the
out-of-time-ordered-correlators, which are briefly described here.

###  II.1 Quantum optical models

Here, we discuss the various forms of the multi-atom-field interaction model,
i.e., the Dicke model. In particular, we study the generalized Dicke model,
the $N$-qubit Dicke model, the Tavis-Cummings model, and the Floquet Dicke
model.

####  II.1.1 The generalized Dicke (GD) model

The generalized Dicke (GD) model consists of $N$ two-level atoms with
transition frequencies $\omega_{a}$, coupled to a single mode of a quantized
radiation field of frequency $\omega_{c}$. The Hamiltonian (for $\hbar=1$) of
this model is

$\displaystyle H_{GD}$
$\displaystyle=\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\omega_{a}}{2}\sum_{j=1}^{N}\sigma_{j}^{z}+\frac{\lambda}{\sqrt{N}}\sum_{j=1}^{N}\left(\hat{a}\sigma^{+}_{j}+\hat{a}^{\dagger}\sigma^{-}_{j}\right)$
$\displaystyle+\frac{\lambda^{\prime}}{\sqrt{N}}\sum_{j=1}^{N}\left(\hat{a}\sigma^{-}_{j}+\hat{a}^{\dagger}\sigma^{+}_{j}\right),$
(1)

which, upon using the collective angular momentum operators
$J_{i}=\frac{1}{2}\sum_{j=1}^{N}\sigma_{j}^{i}$ ($\sigma^{i}$ are the Pauli
matrices with $i={x,y,z}$) and $J_{\pm}=\sum_{j=1}^{N}\sigma^{\pm}$, is given
as

$\displaystyle H_{GD}$
$\displaystyle=\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{a}J_{z}+\frac{\lambda}{\sqrt{N}}\left(J_{+}\hat{a}+J_{-}\hat{a}^{\dagger}\right)+\frac{\lambda^{\prime}}{\sqrt{N}}\left(J_{-}\hat{a}+J_{+}\hat{a}^{\dagger}\right).$
(2)

The collective angular momentum operators $J_{i}$ ($i=z,+,-$), in terms of a
pseudospin of length $j=N/2$, describe the ensemble of two-level atoms. The
single mode quantized radiation field is depicted by the number operator
$\hat{a}^{\dagger}\hat{a}$, where $\hat{a}$, and $\hat{a}^{\dagger}$ are the
bosonic creation and annihilation operators. The third and fourth terms depict
the coupling between the quantized radiation field and the atoms for co- and
counter-rotating terms, respectively. Here $\lambda$ denotes the interaction
strength for the co-rotating term, and $\lambda^{\prime}$ denotes the
interaction strength for the counter-rotating term. $\sqrt{N}$ is the scaling
factor.

Note that if we keep the interaction strengths $\lambda^{\prime}$ equal to
$\lambda$, the above model boils down to the Dicke model [67], and for
$\lambda^{\prime}=0$, it becomes the Tavis-Cummings (TC) model [116, 117].
These models are discussed below. The phase transition of the generalized
Dicke model is well studied in the literature [71, 118, 76, 119]. The
superradiant phase transition is one of the prominent features of the Dicke
model. For the generalized Dicke model with dissipation, the phase transition
manifests itself as a dynamical instability of the cavity and, hence, can be
detected by considering its response function, in particular the retarded
Green’s function [118]. At the superradiant phase transition, the response
becomes an exponentially growing function of time, indicating a dynamical
instability. Mathematically, this corresponds to the poles of the Green’s
function, which in turn correspond to the zeros of the inverse Green’s
function. This leads to the critical point given by

$\displaystyle\lambda_{c}^{2}=$
$\displaystyle\frac{(R-1)(\omega_{a}\omega_{c}+R\kappa\gamma)}{4R^{2}s_{z}}\left(1-\sqrt{1-\frac{R^{2}(\omega_{c}^{2}+\kappa^{2})(\gamma^{2}+\omega_{a}^{2})}{(\omega_{c}\omega_{a}+R\kappa\gamma)^{2}}}\right),$
(3)

for a fixed $\lambda/\lambda^{\prime}$ ratio, where
$R=\frac{1-(\lambda^{\prime}/\lambda)^{2}}{1+(\lambda^{\prime}/\lambda)^{2}}$,
$\gamma$ is the decay rate of the atomic system, and $\kappa$ is the decay
rate of the cavity. The value of $s_{z}=\langle\sigma_{z}\rangle$ corresponds
to the initial state, which is $-1/2$ for a fully polarized state. In Fig. 1,
we can see the phase transition of the generalized Dicke model between the
normal and the superradiant phases. It can be observed that at
$\lambda=\lambda^{\prime}$, the minimum critical coupling is achieved. This
corresponds to the case of the $N$-qubit Dicke model discussed below. The
critical coupling is larger for the case when either of the counter- or co-
rotating terms dominate. In the case, when $\lambda\gg\lambda^{\prime}$, the
model becomes equivalent to the TC model, and the super-radiance can not be
achieved (under dissipation) as discussed below in Sec. II.1.3.

Figure 1: Phase transition in the generalized Dicke (GD) model. The region
marked with green color shows the superradiant phase, and the one marked with
black color shows the normal phase. The parameters are $N=7$,
$\omega_{a}=\omega_{c}=1$, $\kappa=1,~{}\gamma=0,~{}s_{z}=-1/2$. This has been
derived from [118].

####  II.1.2 The $N$-qubit Dicke model

It consists of $N$ two-level atoms with transition frequencies $\omega_{a}$,
coupled to a single mode of a quantized radiation field of frequency
$\omega_{c}$. This can be derived from the generalized Dicke model in Sec.
II.1.1, for equal interaction strengths of both co- and counter-rotating
terms. The Hamiltonian (for $\hbar=1$) of the system is given by

$H_{D}=\omega_{a}J_{z}+\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\lambda}{\sqrt{N}}(J_{+}+J_{-})(\hat{a}+\hat{a}^{\dagger}).$
(4)

The $N$-qubit Dicke model (Eq. (4)) undergoes a superradiant phase transition
[69, 71] when the atom-field interaction strength reaches a critical value of
$\lambda_{c}=\frac{\sqrt{\omega_{a}\omega_{c}}}{2}$. This model is symmetric
under the parity operator
$\Pi=\exp\left[-i\pi\left(\hat{a}^{\dagger}\hat{a}+J_{z}\right)\right]$, and
supports a $\mathbb{Z}_{2}$ symmetry. It is this symmetry that breaks in the
superradiant phase at the critical value of $\lambda_{c}$. The phase
transition in the ground state ($E_{0}$) of the Dicke model is shown in Fig.
(2). Here, we can observe a discontinuity in the second-order derivative of
the ground state of the Hamiltonian, indicating a second-order phase
transition around the critical point $\lambda_{c}$. This critical limit of the
phase transition in the $N$-qubit Dicke model can also be derived from Eq.
(3). In the limit $\lambda\to\lambda^{\prime}$, the value of $R$ goes to 0 and
the critical coupling reduces to
$\lambda_{c}^{2}=-\frac{(\omega_{c}^{2}+\kappa^{2})(\omega_{a}^{2}+\gamma^{2})}{8\omega_{a}\omega_{c}s_{z}}$.
Further, under the condition of zero dissipation ($\gamma=\kappa=0$) and zero
external temperature, for a fully polarized initial state ($s_{z}=-1/2$), this
critical value becomes $\lambda_{c}=\frac{\sqrt{\omega_{a}\omega_{c}}}{2}$.

Figure 2: Variation of the scaled ground state, $\frac{E_{0}}{\omega_{a}j}$
and its derivatives with the atom-field interaction strength $\lambda$ in case
of the $N$-qubit Dicke model (Eq. (4)). Here we have chosen $N=20,\text{and
}~{}\omega_{a}=\omega_{c}=2$.

####  II.1.3 Tavis-Cummings (TC) model

The Tavis-Cummings (TC) model can be arrived at after removing the counter-
rotating terms from the generalized Dicke model (Eq. (2)). The Hamiltonian of
the TC model is given as

$H_{\text{TC}}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\omega_{a}}{2}\sum_{j=1}^{N}\sigma_{j}^{z}+\frac{\lambda}{\sqrt{N}}\sum_{j=1}^{N}\left(\sigma^{+}_{j}\hat{a}+\sigma^{-}_{j}\hat{a}^{\dagger}\right).$
(5)

The total number of excitations
$\mathcal{N}=\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\sum_{j}\sigma^{z}_{j}$ in
this model remain conserved. Using collective angular momentum operators, we
can rewrite the TC Hamiltonian as

$H_{\text{TC}}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{a}J_{z}+\frac{\lambda}{\sqrt{N}}\left(J_{+}\hat{a}+J_{-}\hat{a}^{\dagger}\right).$
(6)

The TC model was shown to undergo a phase transition at a critical value of
atom-cavity field interaction strength
$\lambda_{c}=\sqrt{\omega_{a}\omega_{c}}$ [120, 70] at zero temperature. This
is 2 times the critical coupling obtained for the $N$-qubit Dicke model
because of the presence of rotating wave approximation in this model. This
phase transition can be observed in the ground state of the TC model as shown
in Fig. 3. Here, we can see that for the given parameters
$\omega_{a}=\omega_{c}=2$, the ground state of the TC model becomes non-
analytical after $\lambda=2$. Interestingly, the TC model shows both first-
and second-order phase transitions. The case considered here, with seven atoms
in the TC model, shows a first-order phase transition. Further, the first-
order derivative of the ground state energy has a step-like discontinuity,
which is consistent with the observations in [121, 122]. This step-like
discontinuity arises due to sudden changes in the slope of the ground state
energy function, as depicted in the inset plot of Fig. 3. The derivative of
the ground state energy becomes a continuous function when we consider the
thermodynamic limit ($N\to\infty$) of the number of atoms. In the
thermodynamic limit, the second-order derivative turns out to be
discontinuous, depicting a second-order phase transition. However, under
dissipation, when we consider the losses due to atomic system and cavity
decay, the criticality of the TC model vanishes [121]. This can also be seen
from Eq. (3), where for $\lambda^{\prime}=0$ (i.e., $R=1$), the critical
coupling $\lambda_{c}$ does not have real solutions.

Figure 3: Variation of the scaled ground state energy
$\frac{E_{0}}{\omega_{a}j}$ and its derivatives with the atom-field
interaction strength $\lambda$ for the TC model. The inset plot shows the
variation of the ground state energy $\frac{E_{0}}{\omega_{a}j}$ zoomed in
between $\lambda=2$ and 2.5.

The parameters, in this case, are $N=7$, $\omega_{a}=\omega_{c}=2$.

####  II.1.4 Floquet Dicke model

In the Floquet Dicke (FD) model, we replace the coupling between the atom and
the cavity field with a time-dependent one; that is, we replace the value of
$\lambda$ in Eq. (4) with, for example,
$\lambda(t)=\lambda_{0}+\Delta\lambda\cos(\Omega t)$. Thus, the Hamiltonian
for this model can be explicitly given as

$H_{FD}=\omega_{a}J_{z}+\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\lambda(t)}{\sqrt{N}}(J_{+}+J_{-})(\hat{a}+\hat{a}^{\dagger}),$
(7)

where $J_{i}\\{i=x,y,z\\}$ are the collective angular momentum operators, and
$\hat{a}~{}(\hat{a}^{\dagger})$ are the bosonic annihilation (creation)
operators. It is opportune here to discuss the phase transition tendencies of
this model. In Fig. 4, we can observe the pattern of the scaled ground state
energy of the Floquet Dicke model, varying in time. With values of
$\lambda_{0}=0.65$ and $\Delta\lambda=0.75$, the superradiant phase transition
in the FD model is evident when the time-varying interaction strength crosses
the critical value of $\sqrt{\omega_{a}\omega_{c}}/2=1$. The value of the
ground state energy is constant at -1.0 in the normal phase, and it decreases
quickly as the phase becomes superradiant. Therefore, the system oscillates
between the normal and the superradiant phases.

Figure 4: Variation of the scaled ground state, $\frac{E_{0}}{\omega_{a}j}$
with time $t$ for the Floquet Dicke model. Here we have chosen $N=7$,
$\omega_{a}=\omega_{c}=2$, $\lambda_{0}=0.65$, $\Delta\lambda=0.75$, and
$\Omega=\pi$.

The numerical simulation of the various Dicke models is important as they have
an infinite bosonic Fock space. This entails a truncation of the bosonic Fock
space. The scheme of truncation depends on the observable we want to
calculate. For example, if we want to calculate the ground state energy
$E_{0}$ of the $N$-qubit Dicke model, we arbitrarily choose the dimension $n$
of the bosonic Fock space. We then calculate the difference
$E_{0}(n+1)-E_{0}(n)$. If this difference is significantly small, we fix the
value of $n$. However, in the case of the strong coupling regime, the
dimension of the Hilbert space increases very rapidly, and for large system
sizes the complexity of the calculation increases manifold. To tackle this
issue and also to find an exact numerical solution, particularly in the case
of the $N$-qubit Dicke model, the formalism of the extended bosonic coherent
basis (or efficient basis) was developed in [123], and further elucidated in
[124, 125]. The idea of the conserved parity operator that commutes with the
$N$-qubit Dicke model was used to arrive at this basis. Using this, the
calculation of the observables converges in the case of the $N$-qubit Dicke
model for a higher number of atoms and in a strong coupling regime.

###  II.2 Out-of-time-ordered-correlator (OTOC)

The out-of-time-ordered correlator (OTOC) is defined as

$\mathcal{C}_{AB}(t)=\braket{[A_{t},B]^{\dagger}[A_{t},B]},$ (8)

where $A_{t}$ and $B$ are arbitrary Heisenberg operators in the Hilbert space.
$\braket{\cdot}={\rm Tr}(\rho\cdot)$ denotes the average with respect to the
thermal density matrix of the system $\rho=e^{-\beta H}/Z$, where $Z={\rm
Tr}(e^{-\beta H})$ and $\beta=1/k_{B}T$. Here $H$ is the Hamiltonian of the
system. The intuition for $C(t)$ to be connected to chaos comes from the study
of the semi-classical limit of a one-particle quantum chaotic system as in
[32]. In the semi-classical limit, for $A_{t}=x(t)$ and $B=p$, we can rewrite
the commutator in Eq. (8) as the Poisson bracket
$i\hbar\\{x(t),p\\}=i\hbar\frac{\partial x(t)}{\partial x(0)}$. This Poisson
bracket is sensitive to the initial conditions and dependent on
$e^{\lambda_{L}t}$ for a classically chaotic system. In the context of quantum
chaos, the quasi-classical limit on the growth of OTOC is given by
$\mathcal{C}_{AB}(t)\sim\hbar^{2}e^{2\lambda_{L}t}$. The exponent
$\lambda_{L}$ is referred to as the Lyapunov exponent. A maximum bound on the
Lyapunov exponent $\lambda_{L}$ was proposed in [29], satisfying

$\lambda_{L}\leq\frac{2\pi}{\beta\hbar}=\frac{2\pi k_{B}T}{\hbar}.$ (9)

For Hermitian operators $A_{t}$ and $B$, the form of OTOC defined in Eq. (8)
becomes,

$\mathcal{C}_{AB}^{\text{herm}}(t)=-\braket{[A_{t},B]^{2}}.$ (10)

In general, we have time evolution of the operator
$A_{t}=e^{\frac{i}{\hbar}Ht}Ae^{-\frac{i}{\hbar}Ht}$. In the context of open
quantum systems, here we make use of the adjoint master equation [89, 126] for
the evolution of the Heisenberg operator $A$ given by

$\displaystyle\frac{dA_{t}}{dt}$ $\displaystyle=\mathcal{L}^{\dagger}A_{t},$
(11)

where $\mathcal{L}^{\dagger}$ is the superoperator corresponding to the
adjoint master equation and can be given by
$\mathcal{L}^{\dagger}(\cdot)=i\left[H,(\cdot)\right]+\sum_{k}\left[\frac{\nu_{k}}{2}\left(N_{th,k}+1\right)\left(2L_{k}^{\dagger}(\cdot)L_{k}-\left\\{L_{k}^{\dagger}L_{k},(\cdot)\right\\}\right)\right]+\sum_{k}\left[\frac{\nu_{k}N_{th,k}}{2}\left(2L_{k}(\cdot)L_{k}^{\dagger}-\left\\{L_{k}L_{k}^{\dagger},(\cdot)\right\\}\right)\right]$,
where $H$ is the Hamiltonian of the respective model and
$N_{th,k}=\frac{1}{e^{\beta\hbar\omega_{k}}-1}$. In the case of the Dicke
class of models, we consider the following decay rates: $i).$ decay due to the
spontaneous emission of the atom, and $ii).$ leakage of photons from the
cavity; this is detailed below in Eq. (23). $L_{k}(L_{k}^{\dagger})$ are the
Lindblad operators. The form of OTOC given in Eq. (8), in a generic case, can
be expressed as

$\mathcal{C}_{AB}(t)=\mathcal{D}_{AB}(t)+\mathcal{I}_{AB}(t)-2\Re\left\\{\mathcal{F}_{AB}(t)\right\\},$
(12)

where $\mathcal{D}_{AB}(t)=\braket{B^{\dagger}(A^{\dagger}A)_{t}B}$,
$\mathcal{I}_{AB}=\braket{A_{t}^{\dagger}(B^{\dagger}B)A_{t}}$, and
$\mathcal{F}_{AB}(t)=\braket{A_{t}^{\dagger}B^{\dagger}A_{t}B}$.
$\mathcal{I_{AB}}$ and $\mathcal{D_{AB}}$ are the time-ordered correlation
functions, and $\mathcal{F}_{AB}$ is the genuine out-of-time-ordered four-
point correlation function. This is the reason why $\mathcal{F_{AB}}$ is
commonly referred to as OTOC. The $\mathcal{D}_{AB}(t),\mathcal{C}_{AB}(t)$,
and $\mathcal{I_{AB}}(t)$ are connected to each other via the following
relation [37],

$\frac{\mathcal{C}_{AB}}{\mathcal{I_{AB}}}=\left(1+\alpha_{t}\sqrt{\frac{\mathcal{D}_{AB}}{\mathcal{I_{AB}}}}\right)^{2},$
(13)

where $\alpha_{t}$ is a time-dependent numerical factor satisfying
$|\alpha_{t}|\leq 1$. It can also be seen that in the case where $A$ and $B$
are unitary operators, Eq. (12) can be rewritten as

$\mathcal{C}_{AB}^{\text{uni}}(t)=2\left(1-\Re\left\\{\mathcal{F}_{AB}(t)\right\\}\right).$
(14)

Also, for Hermitian operators $A_{t}$ and $B$, the function
$\mathcal{F}_{AB}(t)$ has a uniform asymptotic expansion given by

$\mathcal{F}_{AB}(t)=c_{0}-\epsilon c_{1}e^{\lambda_{L}t}+O(\epsilon^{2}).$
(15)

In many cases, instead of Eq. (8), a regularized form is adopted for the
definition of OTOC. This is achieved by splitting the thermal density matrices
$\rho$ into two $\rho^{1/2}$’s, one inserted between the commutators and the
other placed in front of them. The definition is given as

$\mathcal{C}_{AB}^{\text{reg}}(t)={\rm
Tr}\left\\{\rho^{1/2}[A_{t},B]^{\dagger}\rho^{1/2}[A_{t},B]\right\\}.$ (16)

In [35], the OTOC was defined in terms of both commutator and anti-commutator
for Hermitian operators $A$ and $B$. The regularized OTOC is redefined in
terms of both (anti-)commutators as

$\mathcal{C}_{AB}^{\text{reg}}(x,y,t,t^{\prime})=\mathcal{C}_{[A,B]_{x}[A,B]_{y}}^{\text{reg}}(t,t^{\prime})={\rm
Tr}\left(\rho^{1/2}[A_{t},B_{t^{\prime}}]_{x}\rho^{1/2}[A_{t},B_{t^{\prime}}]_{y}\right),$
(17)

where $x,y$ can have values $\pm$, and $[.]_{(-)+}$ represents the
(anti-)commutator. The negative of the ordinary OTOC (Eq. (10)) for Hermitian
operators $A$ and $B$ in the form of expectation value ${\rm Tr}(\rho\cdots)$
of products of (anti-)commutator for a given state $\rho$ is given as

$\mathcal{C}_{AB}^{\text{phys}}(x,y,t,t^{\prime})=\mathcal{C}_{[A,B]_{x}[A,B]_{y}}^{\text{phys}}(t,t^{\prime})={\rm
Tr}\left(\rho[A_{t},B_{t^{\prime}}]_{x}[A_{t},B_{t^{\prime}}]_{y}\right).$
(18)

We shall refer to the form of OTOC given above as a physical OTOC, and
$\mathcal{C}_{[A,B]_{x}[A,B]_{y}}^{\text{phys}}(t)=-\mathcal{C}_{AB}^{\text{herm}}(t)$
for x= y = +. The difference between the regularized OTOC and the physical
OTOC is given by the Wigner Yanase (WY) Skew information,

$I_{1/2}(\rho,O)=-\frac{1}{2}{\rm Tr}\left([\rho^{1/2},O]^{2}\right)={\rm
Tr}\left(\rho O^{2}\right)-{\rm Tr}\left(\rho^{1/2}O\rho^{1/2}O\right),$ (19)

for a Hermitian operator $O$. Using the WY skew information, the relation
between the regularized and physical OTOCs is written as

$\displaystyle\mathcal{C}_{[A,B]^{2}}^{\text{reg}}(t,t^{\prime})$
$\displaystyle=\mathcal{C}_{[A,B]^{2}}^{\text{phys}}(t,t^{\prime})+I_{1/2}(\rho,i[A_{t},B_{t^{\prime}}]),$
(20) $\displaystyle\mathcal{C}_{\\{A,B\\}^{2}}^{\text{reg}}(t,t^{\prime})$
$\displaystyle=\mathcal{C}_{\\{A,B\\}^{2}}^{\text{phys}}(t,t^{\prime})-I_{1/2}(\rho,\\{A_{t},B_{t^{\prime}}\\}),$
(21) $\displaystyle\mathcal{C}_{\\{A,B\\}[A,B]}^{\text{reg}}(t,t^{\prime})$
$\displaystyle=\frac{1}{2}\left[\mathcal{C}_{\\{A,B\\}[A,B]}^{\text{phys}}(t,t^{\prime})+\mathcal{C}_{[A,B]\\{A,B\\}}^{\text{phys}}(t,t^{\prime})\right]$
$\displaystyle+\frac{i}{4}I_{1/2}\left(\rho,\\{A_{t},B_{t^{\prime}}\\}+i[A_{t},B_{t^{\prime}}]\right)$
$\displaystyle-\frac{i}{4}I_{1/2}\left(\rho,\\{A_{t},B_{t^{\prime}}\\}-i[A_{t},B_{t^{\prime}}]\right).$
(22)

The physical and regularized forms of OTOC discussed above establish a
relation between the WY Skew information and OTOC. The WY Skew information
metric serves as a measure of information contained in quantum fluctuations of
observables in the context of quantum information theory and provides a
balance between the physical and regularized forms of OTOC. Further, the
regularized form of OTOC has been used in a number of works, such as [29], and
the physical form is defined as the standard expectation value of a quantity.
Moreover, the commutator and anti-commutator forms of OTOC are related via an
out-of-time-ordered Fluctuation-Dissipation relation [35] that relates the
chaotic behavior of the system to the non-linear response function.

##  III Behavior of Out-of-time-ordered-correlator in quantum optical models

Here, we discuss the evolution of out-of-time-ordered correlator (OTOC) for
the quantum optical models discussed in Sec. II.1, which include the
generalized Dicke, $N$-qubit Dicke, Tavis-Cummings, and Floquet Dicke models.
Further, we identify, using OTOC, the regime of these quantum optical models
showing chaotic behavior.

###  III.1 OTOC in the generalized Dicke model

Here, we study the behavior of the OTOC in the generalized Dicke (GD) model.
To this end, we use the hermitian operators
$\hat{q}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ and
$\hat{p}=i(\hat{a}^{\dagger}-\hat{a})/\sqrt{2}$ ($\hat{a}^{\dagger}$ and
$\hat{a}$ are the bosonic creation and annihilation operators) in Eq. (10),
where the operator $\hat{q}$’s evolution in time is governed by the master
equation

$\frac{dq_{t}}{dt}=\mathcal{L}^{\dagger}(q_{t}),$ (23)

where the form of $\mathcal{L}^{\dagger}$ is given in Eq. (11). Here, we
choose the collective angular momentum operators $J_{\pm}=J_{x}\pm iJ_{y}$
(with $J_{x(y)}=\frac{1}{2}\sum_{i}\sigma_{i}^{x(y)}$; $\sigma^{x(y)}$ are the
Pauli spin matrices), and the creation and annihilation operators,
$\hat{a}^{\dagger},~{}\hat{a}$, respectively, as the jump operators. Thus, we
have $L_{1}=J_{-}$, and $L_{2}=\hat{a}$ in Eq. (11), and $\nu_{1}=\gamma$ and
$\nu_{2}=\kappa$ are the decay rates due to spontaneous emission and leakage
of photons from the cavity, respectively. We will now make use of the solution
of this equation $q_{t}$ in the subsequent sections to study the OTOC and its
variants.

Figure 5: Variation of the OTOC
$\mathcal{C}_{\hat{q}\hat{p}}^{\text{herm}}(t)$ (Eq. (10)) or the negative of
physical OTOC in commutator form
$-\mathcal{C}^{\text{phys}}_{[\hat{q}_{t},\hat{p}]^{2}}(t)$, with time for the
generalized Dicke model. Here we have chosen $N=7$, $\omega_{a}=\omega_{c}=2$,
$\lambda=1.5$, $\lambda^{\prime}=1.8$, $\gamma=\kappa=10^{-2}$, and $T=1$.

In Fig. 5, we can observe the variation of the physical OTOC with time. We
have chosen both $\lambda$ and $\lambda^{\prime}$ here to be in the
superradiant phase. An exponential rise in the OTOC can be observed from time
$t=1.5$ to 3.5. We can plot this region on a logarithmic scale to get the
Lyapunov exponent, see above Eq. (9). This is depicted in Fig. 6.

Figure 6: Variation of the OTOC
$\log\left(\mathcal{C}_{\hat{q}\hat{p}}^{\text{herm}}\right)(t)$ (Eq. (10))
and its best fit linear curve from $t=1.5$ to 3.5 for the generalized Dicke
model. Here we have chosen $N=7$, $\omega_{a}=\omega_{c}=2$, and
$\lambda=1.5$, $\lambda^{\prime}=1.8$, $\gamma=\kappa=10^{-2}$, and $T=1$.

The Lyapunov exponent comes out to be 0.30 in the region from $t=1.5$ to 3.5.
This value of the Lyapunov exponent also satisfies the inequality given in Eq.
(9) for $T=1$ in the natural units where $\hbar=k_{B}=1$. We now further
observe the variation of the other forms of OTOC given in Eqs. (17) and (18),
together with the Wigner-Yanase skew information function given in Eq. (19) to
verify the Eqs. (20), (21), and (II.2). This is depicted in Fig. 7.

Figure 7: Variation of the various forms of the OTOCs given in Eqs. (20),
(21), and (II.2), with time for the generalized Dicke model. Here we have
chosen $N=7$, $\omega_{a}=\omega_{c}=2$, and $\lambda=1.5$,
$\lambda^{\prime}=1.8$, $\gamma=\kappa=10^{-2}$, and $T=1$.

In this figure, we can observe that WY Skew information acts as a balance
between the regularized and the physical OTOCs for the commutator, anti-
commutator, and a mix of (anti-)commutator forms.

We further analyze the relation between out-of-time-ordered
$\mathcal{C}_{\hat{q}\hat{p}}(t)$ and time-ordered functions
$\mathcal{D}_{\hat{q}\hat{p}}(t)$, and $\mathcal{I}_{\hat{q}\hat{p}}(t)$ given
in Eq. (13). To this end, we can rewrite the time-varying parameter
$\alpha_{t}$ in terms of $\mathcal{C}_{\hat{q}\hat{p}}(t)$,
$\mathcal{D}_{\hat{q}\hat{p}}(t)$, and $\mathcal{I}_{\hat{q}\hat{p}}(t)$ as

$\alpha_{t}=\sqrt{\frac{\mathcal{I}_{\hat{q}\hat{p}}(t)}{\mathcal{D}_{\hat{q}\hat{p}}(t)}}\left(\sqrt{\frac{\mathcal{C}_{\hat{q}\hat{p}}(t)}{\mathcal{I}_{\hat{q}\hat{p}}(t)}}-1\right).$
(24)

The above relation holds for the value of $|\alpha_{t}|\leq 1$. The variation
of the parameter $|\alpha_{t}|$ with time is plotted in Fig. 8.

Figure 8: Variation of the parameter $|\alpha_{t}|$, and the functions
$\sqrt{\frac{\mathcal{I}_{\hat{q}\hat{p}}(t)}{\mathcal{D}_{\hat{q}\hat{p}}(t)}}$
and
$\sqrt{\frac{\mathcal{C}_{\hat{q}\hat{p}}(t)}{\mathcal{I}_{\hat{q}\hat{p}}(t)}}$
given in Eq. (24), with time in case of generalized Dicke model. Here we have
chosen $N=7$, $\omega_{a}=\omega_{c}=2$, and $\lambda=1.5$,
$\lambda^{\prime}=1.8$, $\gamma=\kappa=10^{-2}$, and $T=1$.

It is observed that the value of $|\alpha_{t}|$ lies between zero and one. The
pattern of the variation $|\alpha_{t}|$ is in contrast to the pattern observed
for the OTOC in Fig. 5. Thus, this behavior of $|\alpha_{t}|$ decreasing
exponentially can point out an exponential increase in the value of OTOC.
Further, we can observe the variation of the various ratios of the out-of-
time-ordered and time-ordered functions with time in Fig. 8. The function
$\mathcal{I}_{\hat{q}\hat{p}}(t)$ is always greater than the function
$\mathcal{D}_{\hat{q}\hat{p}}(t)$. Also, initially the function
$\mathcal{C}_{\hat{q}\hat{p}}(t)$ decreases in comparison to the function
$\mathcal{I}_{\hat{q}\hat{p}}(t)$. Later it increases to approximately match
the value of $\mathcal{I}_{\hat{q}\hat{p}}(t)$ as the graph of
$\sqrt{\frac{\mathcal{C}_{\hat{q}\hat{p}}(t)}{\mathcal{I}_{\hat{q}\hat{p}}(t)}}$
can be seen to approach one. Moreover, the peaks and valleys of the function
$|\alpha_{t}|$ match exactly with the valleys and peaks of the function
$\sqrt{\frac{\mathcal{C}_{\hat{q}\hat{p}}(t)}{\mathcal{I}_{\hat{q}\hat{p}}(t)}}$.

###  III.2 OTOC in the $N$-qubit Dicke model

We now calculate the OTOC for the $N$-qubit Dicke model, introduced in Sec.
II.1.2, using the efficient coherent basis introduced in [123, 124, 125]. In
Fig. 9, we can observe the variation of OTOC $\mathcal{C}_{\hat{q}\hat{p}}(t)$
for different values of the interaction strength between the spin and the
cavity systems $\lambda$.

Figure 9: Variation of the OTOC
$\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ with time for different values
of the spin-cavity system interaction strength $\lambda$ for the $N$-qubit
Dicke model. Here we have chosen $N=25$, $\omega_{a}=\omega_{c}=2$,
$\gamma=\kappa=10^{-2}$, and $T=1$.

It is evident here that the exponential behavior of the OTOC kicks in when the
Dicke Hamiltonian is in the superradiant region, i.e., when $\lambda$ is
greater than the critical value of phase transition in the Dicke model
$\lambda>\lambda_{c}$ ($\approx 1$ for $\gamma=\kappa=10^{-2}$). The Lyapunov
exponent is now obtained from the variation of OTOC for the value of
$\lambda=1.2$ using the plot of the function
$\log\left(\mathcal{C}_{\hat{q}\hat{p}}^{\text{herm}}\right)(t)$ and fitting
it with a straight line. This can be seen in Fig. 10 from time $t=2.0$ to 4.5.
In this case, the value of the Lyapunov exponent comes out to be 0.36, which
satisfies the inequality given in Eq. (9) for $T=1$.

Figure 10: Variation of the OTOC
$\log\left(\mathcal{C}_{\hat{q}\hat{p}}^{\text{herm}}\right)(t)$ (Eq. (10))
and its best fit linear curve from $t=2.0$ to 4.5 for the $N$-qubit Dicke
model. Here we have chosen $N=25$, $\omega_{a}=\omega_{c}=2$, $\lambda=1.2$,
$\gamma=\kappa=10^{-2}$, and $T=1$.

Further, we bring out the effect of dissipation on the growth of the OTOC in
this model. To this end, for different values of the cavity dissipative factor
$\kappa$, see Eq. (23), we plot the variation of OTOC with time for a given
value of spin-cavity system interaction strength $\lambda=1.2$. This is
depicted in Fig. 11.

Figure 11: Variation of the OTOC
$\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ with time for different values
of the dissipative factor $\kappa$ (Eq. (23)), for the $N$-qubit Dicke model.
Here we have chosen $N=25$, $\omega_{a}=\omega_{c}=2$, $\lambda=1.2$,
$\gamma=10^{-2}$, and $T=1$.

It is observed that the OTOC is impacted by dissipation. An increment in the
value of the cavity dissipative factor $\kappa$ significantly reduces the
exponential behavior of the OTOC. For $\kappa=0.5$, the exponential behavior
is completely diminished, and the value of OTOC tends to zero at longer times
due to higher dissipation.

###  III.3 OTOC in the Tavis-Cummings model

Here, we study the behavior of the out-of-time-ordered correlator in the
Tavis-Cummings model, discussed in Sec. II.1.3. The time evolution of the
Hermitian OTOC $\mathcal{C}_{\hat{q}\hat{p}}^{\text{herm}}(t)$ in this model
is depicted in Fig. 12.

Figure 12: Variation of the OTOC
$\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ with time for different values
of the spin-cavity interaction strength $\lambda$ in the case of the TC model.
Here we have chosen $N=7$, $\omega_{a}=\omega_{c}=2$, $\gamma=\kappa=10^{-2}$,
and $T=1$.

The critical value of the system-cavity interaction strength for the phase
transition of the TC model for the given parameters is $\lambda_{c}=2$ for
negligible values of $\gamma$ and $\kappa$, as shown in Fig. 3. Interestingly,
we observe that in this model, the OTOC grows exponentially when the spin-
cavity interaction strength $\lambda$ is greater than the critical value of
$\lambda_{c}=2$. This brings out the connection of the OTOC with the phase
transition, in general, as we observed a similar pattern of OTOC in the case
of the generalized Dicke and $N$-qubit Dicke models, where OTOC showed
exponential behavior in the superradiant phase characterized by second-order
phase transition.

The Lyapunov exponent can be found using curve-fitting, as done previously for
the generalized and the $N$-qubit Dicke models. It comes out to be 0.12 for
variation of OTOC $\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ at the value
of interaction strength $\lambda=4.5$ in Fig. 12.

###  III.4 OTOC in the Floquet Dicke model

In the Floquet Dicke (FD) model, discussed in Sec. II.1.4, the system-cavity
interaction strength is time-dependent, as given in Eq. (7). We plot the
variation of OTOC with time in Fig. 13.

Figure 13: Variation of the OTOC
$\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ with time for the Floquet
Dicke model. Here we have chosen $N=7$, $\lambda=0.65$, $\Delta\lambda=0.75$,
$\Omega=\pi$, $\omega_{c}=2$, $\gamma=\kappa=10^{-2}$, and $T=1$.

It can be observed that the OTOC keeps on increasing to an exponential value
for the time period considered. Upon fitting the curve with a logarithmic plot
of the OTOC, the Lyapunov exponent comes out to be 0.325, which again follows
the inequality of the Eq. (9). An interesting observation can be made when we
compare this figure with the plot given in Fig. 4. The occasional dips in the
value of the OTOC come at the points where the ground state energy of the FD
model is constant, or the model is in a normal phase. At the points in time
where the model shows superradiant phase transitions, the OTOC picks up and
grows exponentially. This observation is in accordance with the behavior of
OTOC observed in the previous models, where phase transition acted as the
cause of the exponential rise of OTOC in the system. Further, the FD model has
correspondence with the $N$-qubit Dicke model if the system-cavity interaction
strength is assumed to be a constant. Here, due to the choice of the
parameters and time dependence of the system-cavity interaction strength, the
FD model oscillates between the normal and the superradiant phases. However,
the maximum value of the OTOC in the case of $N$-qubit Dicke model’s
superradiant phase is lesser than the same in the FD model’s superradiant
phase at longer times. This value keeps rising after every transition from the
normal to the superradiant phase, which is not the case in the $N$-qubit Dicke
model. The rapid oscillations between the normal and the superradiant phases
could be the root cause of the rapid increment in the value of OTOC.

###  III.5 Short-time, long-time, and intermediate-time behavior of OTOC

We have calculated the OTOC of the various quantum optical models above. We
now observe the variation of the OTOC in various time regimes together with
the effect of dissipation at longer times. In Fig. 14, we have plotted the
long-time behavior of the OTOC.

Figure 14: Variation of the OTOC
$\mathcal{C}^{\text{herm}}_{\hat{q}\hat{p}}(t)$ with time in case of (a):
generalized Dicke (with $\lambda=1.5$, and $\lambda^{\prime}=1.2$), (b):
$N$-qubit Dicke (with $\lambda=1.5$), (c): TC model (with $\lambda=4.5$), and
(d): Floquet Dicke models (with $\lambda=0.65$, $\Delta\lambda=0.75$, and
$\Omega=\pi$) at longer times. Here we have chosen
$N=7,\omega_{a}=\omega_{c}=2$, $\gamma=\kappa=10^{-2}$, and $T=1$.

It can be observed that due to the dissipation in the evolution of the system
with time, the OTOC decays after reaching an exponential value. Furthermore,
the OTOC in the case of the Floquet Dicke model reaches the maximum value in
comparison to the other models. The initial time up to which the OTOCs grow
exponentially depicts the short-time behavior of the OTOCs, and these attain a
global maximum value at the scrambling time $t^{*}$. In the case of the TC
model, the scrambling time is around 9. The exponential growth of the OTOCs
stops at the scrambling time. A scrambled system corresponds to a maximal
entanglement entropy of any subsystem when the system is prepared in a pure
state. The time regime in which the OTOC transitions from exponential growth
to decay, around the scrambling time, represents the intermediate time regime.
However, the intermediate time regime for the models considered here can not
be explicitly pointed out because the transition from exponential growth to
decay is very quick. In the long-time regime, we observe that due to
dissipation in the dynamics of the system, the OTOCs decay. At longer times,
the value of the OTOC for the $N$-qubit and generalized Dicke models is higher
than the value of OTOCs for the Floquet Dicke and TC models, which are
approximately zero.

The long-time behavior of the OTOC is important to identify the chaotic
dynamics of a system. The short-time exponential growth as an identifier of
quantum chaos is not universal [36] as it detects the instabilities in the
system as well, and it may not be used as an absolute measure of the quantum
chaos [127, 128, 129]. This can be the reason that the TC model, which is an
integrable model, shows an initial exponential growth of the OTOC. It was
shown in [127] that the long-time behavior of the OTOC for a regular model is
oscillatory, whereas, in the case of a chaotic system, the OTOC saturates to a
fixed value with very small asymptotic oscillations. In this context, for the
models studied here, we observe that the OTOC for the TC model drops quickly
to zero due to dissipation. However, for the generalized Dicke and the
$N$-qubit Dicke models, the OTOC drops after initial exponential growth but
does not go quickly to zero; instead, it tries to saturate to a non-zero value
with small oscillations.

##  IV OTOC in the $N$-qubit Dicke model for all angular momentum sectors

The Dicke model is shown to exhibit both the excited-state quantum phase
transition and the thermal phase transition [130, 131]. This depends on the
angular momentum sectors of the Dicke model. When the maximum angular momentum
sector is taken into account, that is, $j=N/2$, the excited-state quantum
phase transition is observed, and a thermal phase transition occurs when all
the sectors of the angular momentum are used. Consider the model given in Eq.
(1) for $\lambda=\lambda^{\prime}$, at $T\to 0$. The model exhibits quantum
phase transition at
$\lambda=\lambda_{c}=\frac{\sqrt{\omega_{a}\omega_{c}}}{2}$. However, for
finite temperatures, there exists a thermal phase transition. The critical
temperature at which this transition occurs is accounted for using the
following equation [130]

$\beta_{c}=\frac{2}{\omega_{a}}\tanh^{-1}\left(\frac{\omega_{a}\omega_{c}}{4\lambda^{2}}\right),$
(25)

where $\beta_{c}=1/k_{B}T_{c}$. Incidentally, the value of critical
temperature is one if, for example, we consider the values of
$\omega_{a}=\omega_{c}=2$ and $\lambda=1.15$. This value of $\lambda$ is close
to the value of $\lambda_{c}=1.0$ for the $T\to 0$ case.

An interesting thing to do would be to study the behavior of OTOC for the
$N$-qubit Dicke model without restricting to the $j=N/2$ maximally symmetric
angular momentum sector. This would entail taking into account all the sectors
of the angular momentum, which in turn makes the computation very cumbersome.

Figure 15: Variation of the OTOC with time for the $N$-qubit Dicke model
considering all the angular momentum sectors. In subplot (a), we plot the
variation of OTOC for different values of interaction strength (and for
$\kappa=10^{-2}$), and in (b), we plot the variation of OTOC for different
cavity decay rates $\kappa$, keeping interaction strength $\lambda=1.5$. The
parameters chosen are: $N=4$, $\omega_{a}=\omega_{c}=2$, $\gamma=10^{-2}$, and
$T=1$.

The corresponding analysis is depicted in Fig. 15. We observe that the
behavior of OTOC is similar to that observed in the case of the $N$-qubit
Dicke model in Sec. III B for the collective angular momentum operator. In the
regime where the thermal phase transition has occurred, corresponding to
$\lambda\geq 1.15$, the OTOC shows an initial exponential behavior, whereas,
below this critical atom-cavity interaction strength, the exponential behavior
is suppressed (as in Fig. 15(a)). The initial exponential behavior further
diminishes for lower interaction strengths. We further observe the effect of
dissipation on the variation of OTOC for fixed interaction strength, where the
OTOC decays as we increase the cavity decay rate $\kappa$, as depicted in Fig.
15(b).

##  V Connection of OTOC with the characterizers of quantum optics

Here, we discuss the connection between the OTOC and the quantities
characterizing quantum optics, primarily the second-order coherence function
$g^{(2)}(t)$. Based on various conditions on this function, the Poissonian,
sub-Poissonian, bunching, and anti-bunching effects of the cavity field can be
characterized. The second-order coherence $g^{(2)}(t)$ function is one of the
most important characterizers of a light source. It is defined as

$g^{(2)}(t)=\frac{\langle\hat{a}^{\dagger}(t^{\prime})\hat{a}^{\dagger}(t)\hat{a}(t)\hat{a}(t^{\prime})\rangle}{\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle^{2}},$
(26)

where $\hat{a}$ and $\hat{a}^{\dagger}$ are the bosonic annihilation and
creation operators, respectively, and $t=t^{\prime}+\tau$. To this effect, in
this section, we try to elaborate on the behavior of OTOC around these
effects. Up to this point, we computed the OTOC using the operators $\hat{q}$
and $\hat{p}$ as given in Eq. (23). We now make use of the annihilation
operator $\hat{a}$ as the operators $A$ and $B$ given in Eq. (8). The OTOC in
this operator form can now be rewritten as

$\mathcal{C}_{\hat{a}\hat{a}}(t,t^{\prime})=\braket{[\hat{a}_{t},\hat{a}_{t^{\prime}}]^{\dagger}[\hat{a}_{t},\hat{a}_{t^{\prime}}]},$
(27)

where the the annihilation operator $\hat{a}$, is evolved to times
$t^{\prime}$ and $t=t^{\prime}+\tau$, and then their commutator is calculated.
Using Eq. (11), the time evolution of $\hat{a}$ is given by the master
equation

$\frac{d\hat{a}_{t}}{dt}=\mathcal{L}^{\dagger}(\hat{a}_{t}),$ (28)

where $\mathcal{L}^{\dagger}$ is same as given in Eq. (23), where Lindblad
jump operators are $L_{1}=J_{-}$ and $L_{2}=\hat{a}$ with $\nu_{1}=\gamma$ and
$\nu_{2}=\kappa$ as the spontaneous and cavity decay rates, respectively. We
use the above master equation for the time evolution of the operators
$\hat{a}$ and $\hat{a}^{\dagger}$ in the subsequent sections. We now move on
to derive a direct relation between the OTOC defined in Eq. (27) and the
second-order coherence $g^{(2)}(t)$ function.

###  V.1 OTOC and the second-order coherence $g^{(2)}(t)$ function with
superradiant phase transition

To bring out an exact relation between the out-of-time-ordered-correlator
function given in Eq. (8) and the $g^{(2)}(t)$ function, we first rewrite the
OTOC in Eq. (27) in its time-ordered and out-of-time-ordered components. This
is similar to Eq. (12) and can be given as

$\mathcal{C}_{\hat{a}\hat{a}}(t,t^{\prime})=\mathcal{D}_{\hat{a}\hat{a}}(t,t^{\prime})+\mathcal{I}_{\hat{a}\hat{a}}(t,t^{\prime})-2\Re(\mathcal{F}_{\hat{a}\hat{a}}(t,t^{\prime}),$
(29)

where for $t=t^{\prime}+\tau$, we have

$\displaystyle\mathcal{D}_{\hat{a}\hat{a}}(t,t^{\prime})$
$\displaystyle=\langle\hat{a}^{\dagger}(t^{\prime})\hat{a}^{\dagger}(t)\hat{a}(t)\hat{a}(t^{\prime})\rangle,$
$\displaystyle\mathcal{I}_{\hat{a}\hat{a}}(t,t^{\prime})$
$\displaystyle=\langle\hat{a}^{\dagger}(t)\hat{a}^{\dagger}(t^{\prime})\hat{a}(t^{\prime})\hat{a}(t)\rangle,~{}\text{and}$
$\displaystyle\mathcal{F}_{\hat{a}\hat{a}}(t,t^{\prime})$
$\displaystyle=\langle\hat{a}^{\dagger}(t)\hat{a}^{\dagger}(t^{\prime})\hat{a}(t)\hat{a}(t^{\prime})\rangle.$
(30)

Upon a close inspection of the above equation and Eq. (26), it can be seen
that the time-ordered function $\mathcal{D}_{\hat{a}\hat{a}}(t,t^{\prime})$
has a direct correspondence with the $g^{(2)}(t)$ function, and in fact,
$\mathcal{D}_{\hat{a}\hat{a}}(t,t^{\prime})=\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle^{2}g^{(2)}(t)$.
Using this, a direct relation between the out-of-time-ordered correlator
$\mathcal{C}_{\hat{a}\hat{a}}(t,t^{\prime})$ and the second-order coherence
$g^{(2)}(t)$ function can be written as

$\mathcal{C}_{\hat{a}\hat{a}}(t,t^{\prime})=\langle\hat{a}^{\dagger}(t)\hat{a}(t)\rangle^{2}g^{(2)}(t)+\mathcal{I}_{\hat{a}\hat{a}}(t,t^{\prime})-2\Re(\mathcal{F}_{\hat{a}\hat{a}}(t,t^{\prime}).$
(31)

The $g^{(2)}(t)$ function is a four-point correlation function of the field
operators $\hat{a}$ and $\hat{a}^{\dagger}$, and in the model considered here,
this function tells us the behavior of the light interacting with the atomic
system. The OTOC characterizes the dynamics of the system; in particular, it
measures how the evolution of an operator affects the dynamics of the system.
Here, we are able to make a direct comparison between the $g^{(2)}(t)$ and the
OTOC. Practically, the $g^{(2)}(t)$ can be measured numerically, and the
correspondence between this function and the OTOC can pave the way to observe
the chaotic properties of a quantum optical system experimentally.

Figure 16: Variation of the OTOC $\mathcal{C}_{\hat{a}\hat{a}}(t)$ given in
Eq. (27), and the second-order coherence $g^{(2)}(t)$ function (Eq. (26)) with
time $t$ for $t^{\prime}=0$, for the $N$-qubit Dicke model. The parameters are
chosen to be $N=7$, $\omega_{a}=\omega_{c}=2$, $\gamma=\kappa=10^{-2}$, and
$T=1$.

Moreover, we try to re-establish the connection of the OTOC with the
superradiant phase transition of the $N$-qubit Dicke model as given in Sec.
III.2. Here, the OTOC is calculated using the time evolution of the creation
and annihilation operators $\hat{a}$ and $\hat{a}^{\dagger}$, respectively,
instead of $\hat{q}$ and $\hat{p}$ operators used in Sec. III. For the choice
of parameters here, the critical value of atom-field interaction strength
$\lambda_{c}$ is one for superradiant phase transition. In Fig. 16, it can
again be observed that the OTOC has an exponential growth in the superradiant
regime, and it is absent in the normal phase. An interesting comparison
between the $g^{(2)}(t)$ function, OTOC, and the superradiant phase transition
can be made in this case. The frequency of the oscillations in the variation
of the $g^{(2)}(t)$ function and the initial value of this function decreases
as the $N$-qubit Dicke model enters the superradiant phase. In the normal
phase, the dips and the blips of the OTOC and $g^{(2)}(t)$ function do not
match. However, in the superradiant phase, the fluctuations in the OTOC and
$g^{(2)}(t)$ function are in tandem. Further, we explore the optical
properties of light attached to the $g^{(2)}(t)$ function and analyze the
behavior of OTOC in conjunction with these properties next.

###  V.2 Characteristics of light and the OTOC

The second-order $g^{(2)}(t)$ coherence function can be used to characterize
various properties of light. This includes bunching and anti-bunching of
light, as well as super-Poissonian (or classical), Poissonian (or coherent),
and sub-Poissonian (or non-classical) distributions of light. The conditions
on the equal time second-order coherence $g^{(2)}(0)$ function, that is,
whether $g^{(2)}(0)>1$, $g^{(2)}(0)=1$, or $g^{(2)}(0)<1$ identifies whether
the photon distribution of the light is super-Poissonian, Poissonian, or sub-
Poissonian, respectively. Furthermore, the conditions on the $g^{(2)}(t)$
function, that is, whether $g^{(2)}(t)<g^{(2)}(0)$, or $g^{(2)}(t)>g^{(2)}(0)$
determines the photon-counting statistics of light. In the time scale where
$g^{(2)}(t)<g^{(2)}(0)$, the photons distribute themselves in bunches, known
as the bunching effect of the light, whereas when $g^{(2)}(t)>g^{(2)}(0)$,
there are fewer photons that can be detected in bunches; this is the anti-
bunching effect of the light.

Figure 17: Variation of the OTOC $\mathcal{C}_{\hat{a}\hat{a}}(t)$, and the
second-order coherence $g^{(2)}(t)$ function with time in case of the
$N$-qubit Dicke model. The red line shows the value of equal time second-order
coherence $g^{(2)}(0)$ function at $t^{\prime}=0$. Here we have chosen $N=7$,
$\omega_{a}=\omega_{c}=2$, $\lambda=1.5$, $\gamma=\kappa=10^{-2}$, and $T=1$.

In Fig. 17, the variation of OTOC together with the second-order coherence
$g^{(2)}(t)$ function is shown with time for the $N$-qubit Dicke model. It can
be observed that the initial time behavior of OTOC (from $t=1$ to 3.5) in this
case is exponential, and we get a value of $0.51$ for the Lyapunov exponent
upon fitting the curve in this region. The light in the case considered has a
super-Poissonian distribution of photons as $g^{(2)}(0)$ is greater than 1.
Further, here at any time $t$, we have $g^{(2)}(t)<g^{(2)}(0)$. Therefore, the
light is bunched in nature. This is expected as the state of the system is
considered to be in the thermal state. Moreover, a connection between the
$g^{(2)}(t)$ and OTOC $\mathcal{C}_{\hat{a}\hat{a}}(t)$ is observed. The
changes of local maxima and minima in the plot of the $g^{(2)}(t)$ match with
those of the OTOC. Further, in the long-time regime, both $g^{(2)}(t)$
function and OTOC saturates to a fixed value.

##  VI Conclusion

We have studied the quantum chaotic properties of a class of quantum optical
models, viz., the Dicke model, and its variants. To this end, we studied the
generalized Dicke model, which is a multi-atom atom-cavity field interaction
model. From this model, one can derive the two famous atom-cavity field
interaction models, one using equal contributions from the co- and counter-
rotating interaction terms, that is, the $N$-qubit Dicke model. The other one
can be obtained using the rotating wave approximation, that is, the Tavis-
Cummings model. Further, we studied the case of the time-dependent interaction
strength between the atom and cavity field using the Floquet Dicke model. In
all of the above models, we studied the boundaries of the phase transitions.
In the case of the generalized, $N$-qubit, Tavis-Cummings, and Floquet Dicke
models, we discussed the phase transition between the normal and the
superradiant phases. We discussed the various forms of the out-of-time-
ordered-correlators (OTOCs) present in the literature and the connection
between them. This included the commutator and anti-commutator forms of OTOCs,
regularized and physical OTOCs, and the relationship between them using
Wigner-Yanase skew information. Further, we explored the connection between
the out-of-time-ordered and the time-ordered correlators.

We considered a non-unitary evolution of the operators in time to calculate
OTOC using the adjoint form of the GKSL quantum master equation. The
motivation for this was to bring out a more realistic picture of evolution,
which is affected by factors like spontaneous emission and cavity decay. The
OTOC for the class of quantum optical models discussed in this paper was
calculated. We benchmarked the relationship between the various forms of OTOCs
using quantum optical models. The Lyapunov exponent for the exponential
behavior of OTOC was calculated, and the bound on the Lyapunov exponents was
verified. It was observed that the exponential behavior of OTOC kicks in after
the quantum phase transition of the models. The effect of dissipation could
also be seen in the growth of OTOC; as we increase the dissipative factor, the
OTOC decays. The long-time behavior of the OTOC was also explored. For the TC
and the Floquet Dicke models, the OTOC decayed approximately to zero after the
initial exponential rise because of dissipation. However, for the generalized
and $N$-qubit Dicke models, the OTOC exhibited a tendency to saturate after an
initial exponential rise. Taken in conjunction with the short-time exponential
behavior of OTOC, this would suggest quantum chaotic behavior in the
generalized and $N$-qubit Dicke models. We also checked the variation of OTOC
for the disordered Dicke model and its variant, the Richardson-Gaudin model.
In these models, an exponential rise in the growth of OTOC was not observed.

We also established a relationship between the OTOC and the second-order
coherence ($g^{(2)}(t)$) function using the time evolution of bosonic creation
and annihilation operators through the adjoint GKSL master equation. The
relationship between the $g^{(2)}(t)$ function and OTOC is important from an
experimental point of view. In this case, we also studied the behavior of the
OTOC with quantum phase transition of the $N$-qubit Dicke model together with
the changes in the $g^{(2)}(t)$ function. It was observed that the frequency
of the oscillations and the initial value of $g^{(2)}(t)$ function decreases
as we increase the atom-cavity field interaction strength. Below the critical
point of phase transition, where the Dicke model is in the normal phase, the
dynamics of OTOC exhibit non-exponential behavior and becomes exponential in
the superradiant phase. Further, the peaks and dips of the OTOC match with
those of the $g^{(2)}(t)$ function in the superradiant region, which is not
the case in the normal phase. Moreover, for the case considered in this paper,
we observed that the behavior of light is super-Poissonian, and its photons
are bunched in nature.

## Acknowledgements

S.B. acknowledges support from the Interdisciplinary Cyber-Physical Systems
(ICPS) programme of the Department of Science and Technology (DST), India,
Grant No.: DST/ICPS/QuST/Theme-1/2019/6 and DST/ICPS/QuST/Theme-1/2019/13.

## References

  * Casati  _et al._ [1979] G. Casati, B. Chirikov, F. Izraelev, and J. Ford,  _Stochastic Behavior in Classical and Quantum Hamiltonian Systems: Volta Memorial Conference, Como 1977_, edited by G. Casati and J. Ford, Lecture Notes in Physics (Springer Berlin Heidelberg, 1979).
  * Shepelyansky [1983] D. Shepelyansky, Physica D: Nonlinear Phenomena 8, 208 (1983).
  * Haake [2001] F. Haake,  _Quantum Signatures of Chaos_, Physics and astronomy online library (Springer, 2001).
  * Adachi  _et al._ [1988] S. Adachi, M. Toda, and K. Ikeda, Phys. Rev. Lett. 61, 659 (1988).
  * Gutzwiller [1991] M. Gutzwiller,  _Chaos in Classical and Quantum Mechanics_, Interdisciplinary Applied Mathematics (Springer New York, 1991).
  * Zurek and Paz [1994] W. H. Zurek and J. P. Paz, Phys. Rev. Lett. 72, 2508 (1994).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016), https://doi.org/10.1080/00018732.2016.1198134 .
  * Ammann  _et al._ [1998] H. Ammann, R. Gray, I. Shvarchuck, and N. Christensen, Phys. Rev. Lett. 80, 4111 (1998).
  * Mehta [2004] M. Mehta,  _Random Matrices_, Pure and applied mathematics No. v. 142 (Elsevier/Academic Press, 2004).
  * Brody  _et al._ [1981] T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981).
  * Porter [1965] C. Porter,  _Statistical Theories of Spectra: Fluctuations_, Perspectives in Physics Series (Elsevier Science & Technology Books, 1965).
  * Wigner [1967] E. P. Wigner, SIAM Review 9, 1 (1967), https://doi.org/10.1137/1009001 .
  * Dyson [1962] F. J. Dyson, Journal of Mathematical Physics 3, 140 (1962), https://pubs.aip.org/aip/jmp/article-pdf/3/1/140/8157927/140_1_online.pdf .
  * Bohigas and Giannoni [1984] O. Bohigas and M.-J. Giannoni, in  _Mathematical and Computational Methods in Nuclear Physics_ , edited by J. S. Dehesa, J. M. G. Gomez, and A. Polls (Springer Berlin Heidelberg, Berlin, Heidelberg, 1984) pp. 1–99.
  * Kriecherbauer  _et al._ [2001] T. Kriecherbauer, J. Marklof, and A. Soshnikov, Proceedings of the National Academy of Sciences 98, 10531 (2001), https://www.pnas.org/doi/pdf/10.1073/pnas.191366198 .
  * Pandey  _et al._ [2019] A. Pandey, A. Kumar, and S. Puri, “Quantum chaotic systems and random matrix theory,” (2019), arXiv:1905.10596 [cond-mat.stat-mech] .
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52, 1 (1984).
  * Ullmo [2016] D. Ullmo, Scholarpedia 11, 31721 (2016).
  * Chirikov [1985] B. V. Chirikov, Physics Letters A 108, 68 (1985).
  * Flambaum  _et al._ [1999] V. Flambaum, A. Gribakina, G. Gribakin, and I. Ponomarev, Physica D: Nonlinear Phenomena 131, 205 (1999), classical Chaos and its Quantum Manifestations.
  * Peres [1996] A. Peres, Phys. Rev. E 53, 4524 (1996).
  * Borgonovi  _et al._ [2019] F. Borgonovi, F. M. Izrailev, and L. F. Santos, Phys. Rev. E 99, 010101 (2019).
  * Jalabert and Pastawski [2001] R. A. Jalabert and H. M. Pastawski, Phys. Rev. Lett. 86, 2490 (2001).
  * Gorin  _et al._ [2006] T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Physics Reports 435, 33 (2006).
  * Cucchietti  _et al._ [2002] F. M. Cucchietti, C. H. Lewenkopf, E. R. Mucciolo, H. M. Pastawski, and R. O. Vallejos, Phys. Rev. E 65, 046209 (2002).
  * Hayden and Preskill [2007] P. Hayden and J. Preskill, Journal of High Energy Physics 2007, 120 (2007).
  * Sekino and Susskind [2008] Y. Sekino and L. Susskind, Journal of High Energy Physics 2008, 065 (2008).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Journal of High Energy Physics 2014, 67 (2014).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, Journal of High Energy Physics 2016, 106 (2016).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Phys. Rev. D 94, 106002 (2016).
  * Xu and Swingle [2020] S. Xu and B. Swingle, Nature Physics 16, 199 (2020).
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Soviet Journal of Experimental and Theoretical Physics 28, 1200 (1969).
  * Hashimoto  _et al._ [2017] K. Hashimoto, K. Murata, and R. Yoshii, Journal of High Energy Physics 2017, 138 (2017).
  * Rozenbaum  _et al._ [2017] E. B. Rozenbaum, S. Ganeshan, and V. Galitski, Phys. Rev. Lett. 118, 086801 (2017).
  * Tsuji  _et al._ [2018] N. Tsuji, T. Shitara, and M. Ueda, Phys. Rev. E 97, 012101 (2018).
  * García-Mata  _et al._ [2022] I. García-Mata, R. A. Jalabert, and D. A. Wisniacki,  (2022), 10.48550/ARXIV.2209.07965.
  * Hamazaki  _et al._ [2018] R. Hamazaki, K. Fujimoto, and M. Ueda,  (2018), 10.48550/ARXIV.1807.02360.
  * Haehl  _et al._ [2019] F. M. Haehl, R. Loganayagam, P. Narayan, and M. Rangamani, SciPost Phys. 6, 001 (2019).
  * Srednicki [1994] M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Srednicki [1999] M. Srednicki, Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * Kitaev [2015] A. Kitaev,  _A simple model of quantum holography_ , Talk given at the KITP program: entanglement in strongly-correlated quantum matter (2015).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Phys. Rev. Lett. 70, 3339 (1993).
  * Zhang  _et al._ [2021] P. Zhang, S.-K. Jian, C. Liu, and X. Chen, Quantum 5, 579 (2021).
  * Roberts and Stanford [2015] D. A. Roberts and D. Stanford, Phys. Rev. Lett. 115, 131603 (2015).
  * Stanford [2016] D. Stanford, Journal of High Energy Physics 2016, 9 (2016).
  * Chowdhury and Swingle [2017] D. Chowdhury and B. Swingle, Phys. Rev. D 96, 065005 (2017).
  * Patel  _et al._ [2017] A. A. Patel, D. Chowdhury, S. Sachdev, and B. Swingle, Phys. Rev. X 7, 031047 (2017).
  * Nahum  _et al._ [2018] A. Nahum, J. Ruhman, and D. A. Huse, Phys. Rev. B 98, 035118 (2018).
  * Khemani  _et al._ [2018] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev. X 8, 031057 (2018).
  * Rakovszky  _et al._ [2018] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk, Phys. Rev. X 8, 031058 (2018).
  * Weinstein  _et al._ [2022] Z. Weinstein, S. P. Kelly, J. Marino, and E. Altman, arXiv preprint arXiv:2210.14242 (2022), 10.48550/arXiv.2210.14242.
  * Luitz and Bar Lev [2017] D. J. Luitz and Y. Bar Lev, Phys. Rev. B 96, 020406 (2017).
  * Fortes  _et al._ [2020] E. M. Fortes, I. García-Mata, R. A. Jalabert, and D. A. Wisniacki, Europhysics Letters 130, 60001 (2020).
  * Heyl  _et al._ [2018] M. Heyl, F. Pollmann, and B. Dóra, Phys. Rev. Lett. 121, 016801 (2018).
  * Lin and Motrunich [2018] C.-J. Lin and O. I. Motrunich, Phys. Rev. B 97, 144304 (2018).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Nature Physics 13, 781 (2017).
  * Li  _et al._ [2017] J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, Phys. Rev. X 7, 031011 (2017).
  * Wei  _et al._ [2018] K. X. Wei, C. Ramanathan, and P. Cappellaro, Phys. Rev. Lett. 120, 070501 (2018).
  * Niknam  _et al._ [2020] M. Niknam, L. F. Santos, and D. G. Cory, Phys. Rev. Res. 2, 013200 (2020).
  * Loudon [2000] R. Loudon,  _The Quantum Theory of Light_ (Oxford University Press, 2000).
  * Mandel and Wolf [1995] L. Mandel and E. Wolf,  _Optical Coherence and Quantum Optics_, EBL-Schweitzer (Cambridge University Press, 1995).
  * Scully and Zubairy [1997] M. Scully and M. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Agarwal [2013] G. S. Agarwal,  _Quantum Optics_ (Cambridge University Press, 2013).
  * Larson and Mavrogordatos [2021] J. Larson and T. Mavrogordatos,  _The Jaynes–Cummings Model and Its Descendants_, 2053-2563 (IOP Publishing, 2021).
  * Jaynes and Cummings [1963] E. Jaynes and F. Cummings, Proceedings of the IEEE 51, 89 (1963).
  * Shore and Knight [1993] B. W. Shore and P. L. Knight, Journal of Modern Optics 40, 1195 (1993), https://doi.org/10.1080/09500349314551321 .
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Hepp and Lieb [1973a] K. Hepp and E. H. Lieb, Annals of Physics 76, 360 (1973a).
  * Hepp and Lieb [1973b] K. Hepp and E. H. Lieb, Phys. Rev. A 8, 2517 (1973b).
  * Wang and Hioe [1973] Y. K. Wang and F. T. Hioe, Phys. Rev. A 7, 831 (1973).
  * Hioe [1973] F. T. Hioe, Phys. Rev. A 8, 1440 (1973).
  * Carmichael  _et al._ [1973] H. Carmichael, C. Gardiner, and D. Walls, Physics Letters A 46, 47 (1973).
  * Garraway [2011] B. M. Garraway, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369, 1137 (2011), https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2010.0333 .
  * Baumann  _et al._ [2010] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Nature 464, 1301 (2010).
  * Dimer  _et al._ [2007] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael, Phys. Rev. A 75, 013804 (2007).
  * Zhiqiang  _et al._ [2017] Z. Zhiqiang, C. H. Lee, R. Kumar, K. J. Arnold, S. J. Masson, A. S. Parkins, and M. D. Barrett, Optica 4, 424 (2017).
  * Domokos and Ritsch [2002] P. Domokos and H. Ritsch, Phys. Rev. Lett. 89, 253003 (2002).
  * Baden  _et al._ [2014] M. P. Baden, K. J. Arnold, A. L. Grimsmo, S. Parkins, and M. D. Barrett, Phys. Rev. Lett. 113, 020408 (2014).
  * Emary and Brandes [2003a] C. Emary and T. Brandes, Phys. Rev. Lett. 90, 044101 (2003a).
  * Emary and Brandes [2003b] C. Emary and T. Brandes, Phys. Rev. E 67, 066203 (2003b).
  * Altland and Haake [2012] A. Altland and F. Haake, Phys. Rev. Lett. 108, 073601 (2012).
  * Alavirad and Lavasani [2019] Y. Alavirad and A. Lavasani, Phys. Rev. A 99, 043602 (2019).
  * Bastarrachea-Magnani  _et al._ [2015] M. A. Bastarrachea-Magnani, B. López-del Carpio, S. Lerma-Hernández, and J. G. Hirsch, Physica Scripta 90, 068015 (2015).
  * Chávez-Carlos  _et al._ [2016] J. Chávez-Carlos, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, Phys. Rev. E 94, 022209 (2016).
  * Villaseñor  _et al._ [2022] D. Villaseñor, S. Pilatowsky-Cameo, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Entropy 25, 8 (2022).
  * Wang [2022] Q. Wang, Entropy 24, 1415 (2022).
  * Chávez-Carlos  _et al._ [2019] J. Chávez-Carlos, B. López-del Carpio, M. A. Bastarrachea-Magnani, P. Stránský, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Phys. Rev. Lett. 122, 024101 (2019).
  * Prasad  _et al._ [2022] M. Prasad, H. K. Yadalam, C. Aron, and M. Kulkarni, Phys. Rev. A 105, L050201 (2022).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Banerjee [2018] S. Banerjee,  _Open Quantum Systems_ (Springer Singapore, 152 Beach Rd, 21-01 Gateway East, Singapore 189721, 2018).
  * Weiss [1999] U. Weiss,  _Quantum Dissipative Systems_ , 2nd ed., Series in Modern Condensed Matter Physics, Vol. 10 (World Scientific, 1999).
  * Louisell [1990] W. H. Louisell,  _Quantum Statistical Properties of Radiation_ (1990).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Annals of Physics 149, 374 (1983).
  * Grabert  _et al._ [1988] H. Grabert, P. Schramm, and G.-L. Ingold, Physics Reports 168, 115 (1988).
  * Banerjee and Ghosh [2003] S. Banerjee and R. Ghosh, Phys. Rev. E 67, 056120 (2003).
  * Banerjee and Ghosh [2000] S. Banerjee and R. Ghosh, Phys. Rev. A 62, 042105 (2000).
  * Banerjee  _et al._ [2016] S. Banerjee, A. K. Alok, and R. MacKenzie, The European Physical Journal Plus 131, 1 (2016).
  * Naikoo  _et al._ [2018] J. Naikoo, A. K. Alok, and S. Banerjee, Phys. Rev. D 97, 053008 (2018).
  * Dixit  _et al._ [2019] K. Dixit, J. Naikoo, S. Banerjee, and A. K. Alok, The European Physical Journal C 79, 1 (2019).
  * Bhattacharya  _et al._ [2018] S. Bhattacharya, S. Banerjee, and A. K. Pati, Quantum Information Processing 17, 236 (2018).
  * Utagi  _et al._ [2020] S. Utagi, R. Srikanth, and S. Banerjee, Scientific Reports 10, 15049 (2020).
  * Kumar  _et al._ [2018] N. P. Kumar, S. Banerjee, R. Srikanth, V. Jagadish, and F. Petruccione, Open Systems & Information Dynamics 25, 1850014 (2018), https://doi.org/10.1142/S1230161218500142 .
  * Tanimura [2020] Y. Tanimura, The Journal of Chemical Physics 153, 020901 (2020), https://doi.org/10.1063/5.0011599 .
  * Hughes  _et al._ [2009] K. H. Hughes, C. D. Christ, and I. Burghardt, The Journal of Chemical Physics 131, 124108 (2009), https://doi.org/10.1063/1.3226343 .
  * Iles-Smith  _et al._ [2014] J. Iles-Smith, N. Lambert, and A. Nazir, Phys. Rev. A 90, 032114 (2014).
  * Hu and Matacz [1994] B. L. Hu and A. Matacz, Phys. Rev. D 49, 6612 (1994).
  * Srikanth and Banerjee [2008] R. Srikanth and S. Banerjee, Phys. Rev. A 77, 012318 (2008).
  * Huelga and Plenio [2013] S. Huelga and M. Plenio, Contemporary Physics 54, 181 (2013), https://doi.org/10.1080/00405000.2013.829687 .
  * Zhang  _et al._ [2019] Y.-L. Zhang, Y. Huang, and X. Chen, Phys. Rev. B 99, 014303 (2019).
  * Tuziemski [2019] J. Tuziemski, Phys. Rev. A 100, 062106 (2019).
  * Touil and Deffner [2021] A. Touil and S. Deffner, PRX Quantum 2, 010306 (2021).
  * Zanardi and Anand [2021] P. Zanardi and N. Anand, Phys. Rev. A 103, 062214 (2021).
  * Han  _et al._ [2022] L.-P. Han, J. Zou, H. Li, and B. Shao, Entropy 24, 1532 (2022).
  * Andreadakis  _et al._ [2023] F. Andreadakis, N. Anand, and P. Zanardi, Phys. Rev. A 107, 042217 (2023).
  * Blocher and Mølmer [2019] P. D. Blocher and K. Mølmer, Phys. Rev. A 99, 033816 (2019).
  * Tavis and Cummings [1968] M. Tavis and F. W. Cummings, Phys. Rev. 170, 379 (1968).
  * Tavis and Cummings [1969] M. Tavis and F. W. Cummings, Phys. Rev. 188, 692 (1969).
  * Shchadilova  _et al._ [2020] Y. Shchadilova, M. M. Roses, E. G. Dalla Torre, M. D. Lukin, and E. Demler, Phys. Rev. A 101, 013817 (2020).
  * Kirton  _et al._ [2019] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Advanced Quantum Technologies 2, 1800043 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qute.201800043 .
  * Hepp and Lieb [1973c] K. Hepp and E. H. Lieb, Annals of Physics 76, 360 (1973c).
  * Larson and Irish [2017] J. Larson and E. K. Irish, Journal of Physics A: Mathematical and Theoretical 50, 174002 (2017).
  * Castaños  _et al._ [2009] O. Castaños, R. López-Peña, E. Nahmad-Achar, J. G. Hirsch, E. López-Moreno, and J. E. Vitela, Physica Scripta 79, 065405 (2009).
  * Chen  _et al._ [2008] Q.-H. Chen, Y.-Y. Zhang, T. Liu, and K.-L. Wang, Phys. Rev. A 78, 051801 (2008).
  * Bastarrachea-Magnani and Hirsch [2014] M. A. Bastarrachea-Magnani and J. G. Hirsch, Physica Scripta 2014, 014005 (2014).
  * Hirsch and Bastarrachea-Magnani [2014] J. G. Hirsch and M. A. Bastarrachea-Magnani, Physica Scripta 2014, 014018 (2014).
  * Sǎndulescu and Scutaru [1987] A. Sǎndulescu and H. Scutaru, Annals of Physics 173, 277 (1987).
  * Novotný and Stránský [2023] J. Novotný and P. Stránský, Phys. Rev. E 107, 054220 (2023).
  * Hummel  _et al._ [2019] Q. Hummel, B. Geiger, J. D. Urbina, and K. Richter, Phys. Rev. Lett. 123, 160401 (2019).
  * Pilatowsky-Cameo  _et al._ [2020] S. Pilatowsky-Cameo, J. Chávez-Carlos, M. A. Bastarrachea-Magnani, P. Stránský, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Phys. Rev. E 101, 010202 (2020).
  * Pérez-Fernández and Relaño [2017] P. Pérez-Fernández and A. Relaño, Phys. Rev. E 96, 012121 (2017).
  * Bastarrachea-Magnani  _et al._ [2016] M. A. Bastarrachea-Magnani, S. Lerma-Hernández, and J. G. Hirsch, Journal of Statistical Mechanics: Theory and Experiment 2016, 093105 (2016).


# Debye entropic force and modified Newtonian dynamics

Xin Li1,3<EMAIL_ADDRESS>Zhe Chang2,3<EMAIL_ADDRESS>1Institute of
Theoretical Physics, Chinese Academy of Sciences, 100190 Beijing, China
2Institute of High Energy Physics, Chinese Academy of Sciences, 100049
Beijing, China
3Theoretical Physics Center for Science Facilities, Chinese Academy of
Sciences

###### Abstract

Verlinde has suggested that the gravity has an entropic origin, and a
gravitational system could be regarded as a thermodynamical system. It is
well-known that the equipartition law of energy is invalid at very low
temperature. Therefore, entropic force should be modified while the
temperature of the holographic screen is very low. It is shown that the
modified entropic force is proportional to the square of the acceleration,
while the temperature of the holographic screen is much lower than the Debye
temperature $T_{D}$. The modified entropic force returns to the Newton’s law
of gravitation while the temperature of the holographic screen is much higher
than the Debye temperature. The modified entropic force is connected with
modified Newtonian dynamics (MOND). The constant $a_{0}$ involved in MOND is
linear in the Debye frequency $\omega_{D}$, which can be regarded as the
largest frequency of the bits in screen. We find that there do have a strong
connection between MOND and cosmology in the framework of Verlinde’s entropic
force, if the holographic screen is taken to be bound of the Universe. The
Debye frequency is linear in the Hubble constant $H_{0}$.

###### pacs:

98.80.Cq,95.35.+d

##  I Introduction

Gravity, the most universal force in nature, is still unclear in quantum
level. The profound relation between black hole entropy and thermodynamics,
which proposed by Hawking et al.Hawking , throws new light to the genius of
gravity. This fact leads some physicists investigating gravity from a point of
view of thermodynamics. JacobsonJacobson  derived Einstein’s gravitational
field equation from the first law of thermodynamics. The further investigation
of PadmanabhanPadmanabhan1  reinterpreted entropy as the equipartition rule of
energy. It also provided a thermodynamic interpretation of gravity. These
results support that the gravity has an entropic originPadmanabhan2 .
Recently, VerlindeVerlinde  proposed an interesting scenario that interpret
gravity as entropic force caused by the changes in the information associated
with the positions of material bodies. Verlinde’s scenario has two important
consequences. One is that, by making use of the the entropic force together
with the Unruh temperatureUnruh , one could derive the second law of Newton.
The other is that, by making use of the the entropic force together with the
holographic principleSusskind  and the equipartition law of energy, one could
derive Newton’s law of gravitation and Komar massKomar  in relativistic case.

In the framework of entropic force, various aspects of gravitational physics
have been discussed. A succinct list includes: a Friedmann equation derived
from entropic forceShu ; Cai1 , the implication for loop quantum gravitySmolin
, the construction of holographic actions from black hole entropyMakela ;
Caravelli , the holographic dark energy derived from the entropic forceLiM ,
the correct entropy/area relation in the relativistic caseTian , and the
modified entropic force derived from the modified the equipartition law of
energyGao . More related works can be seen in Cai ; Easson ; Banerjee ;
Padmanabhan3 ; Konoplya ; Sheykhi ; Wei ; Zhao

It should be noticed that Verlinde used the free particle’s the equipartition
law of energy to derive the Newton’s law of gravitation. It is well known that
the equipartition law of energy does not hold at very low temperatures, and
the equipartition law of energy derived from Debye model is in good agreement
with experimental results for most of the solid objects in very low
temperatures. Verlinde got the Newton’s law of gravitation with the assumption
that each bit on holographic screen is free of interaction. It should be more
general that the bits on holographic screen interact each others. In such
case, one could anticipate that the Newton’s law of gravitation must be
modified. For example, GaoGao  used the three dimensional Debye model, which
modified the equipartition law of energy, to modify the entropic force. Such
modification can interpret the current acceleration of the Universe while
without invoking any kind of dark energy.

Since the equipartition law of energy in Verlinde’s scenarioVerlinde  is one
dimensional free particle’s equipartition law of energy, we use the one
dimensional Debye model to modified the entropic force. We find that this
modified entropic force is able to derive the famous modified Newtonian
dynamics(MOND). MOND was constructed to explaining the flat rotational curves
of spiral galaxies.

There are a great variety of observations which show that the rotational
velocity curves of all spiral galaxies tend to some constant valuesTrimble .
These include the Oort discrepancy in the disk of the Milky WayBahcall , the
velocity dispersions of dwarf Spheroidal galaxiesVogt , and the flat rotation
curves of spiral galaxiesRubin . These facts violate sharply the prediction of
Newtonian dynamics or Newton’s gravity.

The most widely adopted way to resolve these difficulties is the dark matter
hypothesis. It is assumed that all visible stars are surrounded by massive
nonluminous matters. Another approach is the famous MONDMilgrom . It assumed
that the Newtonian dynamics does not hold in galactic scale. The particular
form of MOND is given as

$\begin{array}[]{l}m\mu\left(\displaystyle\frac{a}{a_{0}}\right)\mathbf{a}=\mathbf{F},\\\\[11.38092pt]
\displaystyle\lim_{x\gg 1}\mu(x)=1,~{}~{}~{}\lim_{x\ll 1}\mu(x)=x,\end{array}$
(1)

where $a_{0}$ is at the order of $10^{-8}$ cm/s2. At beginning, as a
phenomenological model, MOND explains well the flat rotation curves with a
simple formula and a new parameter. In particular, it deduces naturally a
well-known global scaling relation for spiral galaxies, the Tully-Fisher
relationTF . By introducing several scalar, vector and tensor fields,
BekensteinBekenstein  rewrote the MOND into a covariant formulation. He showed
that the MOND satisfies all four classical tests on Einstein’s general
relativity in Solar system.

In this paper, we will derive the MOND in the framework of Verlinde’s entropic
force.

##  II MOND from modified entropic force

In verlinde’s scenario, the change of a particle’s position ($\triangle x$)
corresponds the change of entropy ($\triangle S$) associated with the
information on the boundary. Motivated by Bekenstein’s entropic bound,
Verlinde postulated that the change of entropy is proportional to the change
of a particle’s position,

$\triangle S=2\pi k_{B}\frac{mc}{\hbar}\triangle x.$ (2)

Such particle will experience an effective force when the holographic screen
carries a temperature $T$, together with the first law of thermodynamics, the
force equals to

$F\triangle x=T\triangle S.$ (3)

According to the Unruh formula, the temperature in (3) corresponds the
acceleration

$k_{B}T=\frac{1}{2\pi}\frac{\hbar a}{c}$ (4)

where $a$ denotes the acceleration. Since the holographic screen carries
temperature, the equipartition law of energy gives the relation between the
total energy of the screen and the temperature

$E=\frac{1}{2}Nk_{B}T,$ (5)

where $N$ denotes the number of degrees of freedom (bits) on the screen. The
number of bits $N$ on the screen is assumed to proportional to the area $A$ of
the screen

$N=\frac{Ac^{3}}{G\hbar}.$ (6)

Verlinde assumed that the energy of the screen is proportional to the mass
that would emerge in the part of space enclosed by the screen

$E=Mc^{2}.$ (7)

By making use of the Eq. (2), (3) and (4), one could obtain the second law of
Newton $F=ma$. By making use of the Eq. (2–7) together with the relation
$A=4\pi R^{2}$, one could obtain Newton’s law of gravitation
$F=G\frac{Mm}{R^{2}}$.

In a statistic system, there are interactions among molecules. Thus the
equipartition law of energy for free molecule only valid for the situation
that the kinetic energy of molecule is much larger than the effective
potential of the interaction between each molecule. Therefore, the
equipartition law of energy is invalid at very low temperatures. It is found
that Debye model, which modified the equipartition law of energy, is in good
agreement with experimental results for most solid objects. Following
Verlinde’s scenario, we know that the gravity can be explained as an entropic
force, it means that the gravity may have a statistical thermodynamics
explanation. If so, the gravity should modified while the corresponded
statistics changes.

Here, we modified the equipartition law of energy as one dimensional Debye
model

$E=\frac{1}{2}Nk_{B}T\mathfrak{D}(y),$ (8)

where the one dimensional Debye function is defined as

$\mathfrak{D}(y)\equiv\frac{1}{y}\int^{y}_{0}\frac{z}{e^{z}-1}dz.$ (9)

$y$ is related to the Debye frequency $\omega_{D}$, and defined as

$y\equiv\frac{\hbar\omega_{D}}{k_{B}T}=\frac{2\pi c\omega_{D}}{a}.$ (10)

Following Verlinde’s method, we obtain the modification of Newton’s law of
gravitation

$\frac{GM}{R^{2}}=a\mathfrak{D}(\frac{2\pi c\omega_{D}}{a}).$ (11)

There lies two limit case for the modification of Newton’s law of gravitation
(11). One is the high temperature limit $y\ll 1$, the Debye function
$\mathfrak{D}(y)$ reads

$\mathfrak{D}(y)\approx\frac{1}{y}\int^{y}_{0}dy=1.$ (12)

Then, the modified equipartition law of energy (8) returns to (5). Therefore,
at very high temperature, the Newtonian gravity is recovered. The other is the
low temperature limit $y\gg 1$, the Debye function $\mathfrak{D}(y)$ reads

$\mathfrak{D}(y)\approx\frac{1}{y}\int^{\infty}_{0}\frac{z}{e^{z}-1}dz=\frac{\pi^{2}}{6y}.$
(13)

Then, the modification of Newton’s law of gravitation (11) reads

$\frac{GM}{R^{2}}=\frac{\pi}{12c\omega_{D}}a^{2}=\frac{a^{2}}{a_{0}},$ (14)

where the constant $a_{0}$ is defined as

$a_{0}=\frac{12c\omega_{D}}{\pi}.$ (15)

Combining the Eq. (11) and (15), we find that

$\displaystyle\mathfrak{D}(y)=\mathfrak{D}(\frac{\pi^{2}a_{0}}{6a})$
$\displaystyle=$
$\displaystyle\frac{6}{\pi^{2}}\frac{a}{a_{0}}\int_{0}^{\frac{\pi^{2}a_{0}}{6a}}\frac{z}{e^{z}-1}dz$
(16)

Defining the function

$\mu(x)=\frac{6}{\pi^{2}}x\int_{0}^{\frac{\pi^{2}}{6x}}\frac{z}{e^{z}-1}dz,$
(17)

where $x=\frac{a}{a_{0}}$, one could find that this function has properties

$\lim_{x\gg 1}\mu(x)=1,~{}~{}~{}\lim_{x\ll 1}\mu(x)=x.$ (18)

Therefore, the modification of Newton’s law of gravitation (11) reads

$\frac{GM}{R^{2}}=a\mathfrak{D}(\frac{2\pi c\omega_{D}}{a})=a\mu(x).$ (19)

Equation (19) is just the combination of the MOND relation (1) and Newton’s
law of gravitation. Thus, the modified entropic force (11) could account for
the MOND. To specific this point, by making use of the formula (19) we give
the modified Poisson equation for the gravitational potential $\phi$

$\nabla\cdot(\mu(|\nabla\phi|/a_{0})\nabla\phi)=4\pi G\rho,$ (20)

where $\rho$ is the energy density of matter sources. This equation is just
modified Poisson equation in MONDMilgrom1 .

MilgromMilgrom ; Milgrom1  observed the relation between the MOND constant
$a_{0}$ and the Hubble constant $H_{0}$ that

$2\pi a_{0}\approx cH_{0}.$ (21)

In MONDMilgrom1 , this relation (21) is just a coincidences, and it may point
to a strong connection between MOND and cosmology. In Verlinde’s scenario, we
show that this coincidence is just a consequence of modified entropic force.
By making use of the formula (15) and (21), we find that the Debye frequency
$\omega_{D}$ is proportional to the Hubble constant $H_{0}$

$H_{0}\approx 24\omega_{D}.$ (22)

Since the Debye frequency $\omega_{D}$ is the largest oscillate frequency of
the particle on holographic screen, its wave length $\lambda<c/\omega_{D}$
should be less than scale of the screen. On the other hand, the cosmological
horizon $L_{h}\sim c/H_{0}$ has the same scale as the holographic screen, for
the screen could be regarded as the bound of the Universe. The relation (22)
tells us that

$\lambda<c/\omega_{D}=c/24H_{0}<c/H_{0}.$ (23)

Thus, the relation (22) is reasonable.

##  III conclusion

In this paper, we find that a modified entropic force is connected with
modified Newtonian dynamics. In Verlinde’s scenario, gravity has an entropic
origin and gravitational system could be regarded as a thermodynamical system.
However, it is well-known that the equipartition law of energy is invalid at
very low temperature. Therefore, entropic force should be modified while the
temperature of the holographic screen is very low. It is shown that the
modified entropic force is proportional to the square of the acceleration,
while the temperature of the holographic screen is much lower that the Debye
temperature $T_{D}=\frac{\hbar\omega_{D}}{k_{B}}$. And the modified entropic
force returns to the Newton’s law of gravitation while the temperature of the
holographic screen is much larger that the Debye temperature. This
thermodynamical approach gives a possible origin of the famous MOND. Also, the
constant $a_{0}$ involved in MOND is linear in the Debye frequency
$\omega_{D}$, which can be regarded as the largest frequency of the bits in
screen. We find that there do have a strong connection between MOND and
cosmology in the framework of Verlinde’s entropic force, if the holographic
screen is the bound of the Universe. It is shown that the Debye frequency is
linear in the Hubble constant $H_{0}$. A specific microscopic statistical
thermodynamical model of spacetime may give the origin of our modified
entropic force model and the microscopic origin for the MOND and the
cosmology. We will study it in the future.

###### Acknowledgements.

We would like to thank Prof. C. J. Zhu, Y. Tian and X. N. Wu for useful
discussions. The work was supported by the NSF of China under Grant No.
10525522 and 10875129\.

## References

  * (1) J.M. Bardeen, B. Carter and S.W. Hawking, Commun. Math. Phys. 31, 161 (1973). J.D. Bekenstein, Phys. Rev. D 7, 949 (1973). J.D. Bekenstein, Phys. Rev. D 7, 2333 (1973). S. W. Hawking, Commun. Math. Phys. 43, 199 (1975) [Erratum-ibid. 46, 206 (1976)].
  * (2) T. Jacobson, Phys. Rev. Lett. 75, 1260 (1995).
  * (3) T. Padmanabhan, Class. Quantum Grav. 21, 4485 (2004).
  * (4) T. Padmanabhan, Thermodynamical Aspects of Gravity: New Insights, arXiv:gr-qc/0911.5004.
  * (5) E. P. Verlinde, On the Origin of Gravity and the Laws of Newton, arXiv:hep-th/1001.0785.
  * (6) W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * (7) L. Susskind, J. Math. Phys. 36, 6377 (1995). G. t Hooft, arXiv:gr-qc/9310026. E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998).
  * (8) A. Komar, Phys. Rev. 113, 934 (1959).
  * (9) F. W. Shu and Y. Gong, “Equipartition of energy and the first law of thermodynamics at the apparent horizon, ” arXiv:gr-qc/1001.3237.
  * (10) R. G. Cai, L. M. Cao and N. Ohta, “Friedmann Equations from Entropic Force,” arXiv:hep-th/1001.3470.
  * (11) L. Smolin, “Newtonian gravity in loop quantum gravity,” arXiv:gr-qc/1001.3668.
  * (12) J, Mäkelä,“Notes Concerning ‘On the Origin of Gravity and the Laws of Newton’ by E. Verlinde,” arXiv:gr-qc/1001.3808.
  * (13) F. Caravelli and L. Modesto, “Holographic actions from black hole entropy,” arXiv:gr-qc/1001.4364.
  * (14) M, Li and Y. Wang, “Quantum UV/IR Relations and Holographic Dark Energy from Entropic Force,” arXiv:hep-th/1001.4466.
  * (15) Y. Tian and X. N. Wu, “Thermodynamics of Black Holes from Equipartition of Energy and Holography,” arXiv:hep-th/1002.1275.
  * (16) C. J. Gao, “Modified Entropic Force,” arXiv:hep-th/1001.4585.
  * (17) Y. F. Cai, J. Liu and H. Li, arXiv:astro-ph/1003.4526.
  * (18) D. A. Easson, P. H. Frampton and G. F. Smoot, arXiv:hep-th/1002.4278; arXiv:hep-th/1003.1528.
  * (19) R. Banerjee and B. R. Majhi, arXiv:gr-qc/1003.2312.
  * (20) T. Padmanabhan, arXiv:gr-qc/0912.3165; arXiv:gr-qc/1003.5665.
  * (21) R. A. Konoplya, arXiv:hep-th/1002.2818.
  * (22) A. Sheykhi, Phys. Rev. D 81, 104011 (2010).
  * (23) S. W. Wei, Y. X. Liu and Y.Q. Wang, arXiv:hep-th/1001.5238; arXiv:hep-th/1002.1062.
  * (24) Y. Zhao, arXiv:hep-th/1002.4039.
  * (25) V. T. Trimble, Ann. Rev. Astron. Astrophys. 25, 425 (1987).
  * (26) J. N. Bahcall, C. Flynn, and A. Gould, Astrophys. J. 389, 234 (1992).
  * (27) S. S. Vogt, M. Mateo, E. W. Olszewski, and M. J. Keane, Astron. J. 109, 151 (1995).
  * (28) V. C. Rubin, W. K. Ford, and N. Thonnard, Astrophys. J. 238, 471 (1980).
  * (29) M. Milgrom, Astrophys. J. 270, 365 (1983).
  * (30) R. B. Tully and J. R. Fisher, Astr. Ap. 54, 661 (1977).
  * (31) J. D. Bekenstein, Phys. Rev. D 70, 083509 (2004).
  * (32) M. Milgrom, “The MOND paradigm,” arXiv:astro-ph/0801.3133.