# Non-equilibrium dynamics of the Jaynes-Cummings dimer

G. Vivek, Debabrata Mondal, and S. Sinha  Indian Institute of Science
Education and Research-Kolkata, Mohanpur, Nadia-741246, India

###### Abstract

We investigate the non-equilibrium dynamics of a Josephson coupled Jaynes-
Cummings dimer in the presence of Kerr nonlinearity, which can be realized in
the cavity and circuit quantum electrodynamics systems. The semiclassical
dynamics is analyzed systematically to chart out a variety of photonic
Josephson oscillations and their regime of stability. Different types of
transitions between the dynamical states lead to the self trapping phenomenon,
which results in photon population imbalance between the two cavities. We also
study the dynamics quantum mechanically to identify characteristic features of
different steady states and to explore fascinating quantum effects, such as
spin dephasing, phase fluctuation, and revival phenomena of the photon field,
as well as the entanglement of spin qubits. For a particular ‘self trapped’
state, the mutual information between the atomic qubits exhibits a direct
correlation with the photon population imbalance, which is promising for
generating photon mediated entanglement between two non interacting qubits in
a controlled manner. Under a sudden quench from stable to unstable regime, the
photon distribution exhibits phase space mixing with a rapid loss of
coherence, resembling a thermal state. Finally, we discuss the relevance of
the new results in experiments, which can have applications in quantum
information processing and quantum technologies.

##  I Introduction

Recent advancements in cavity and circuit quantum electrodynamics (QED) have
paved the way to study the non-equilibrium dynamics in quantum systems [1, 2,
3], apart from their potential application to quantum information processing
[2, 3]. Moreover, such atom-photon interacting systems exhibit various
fascinating phenomena, some of which include quantum phase transition [4, 5,
6], the onset of chaos [5], thermalization [7, 8], and the formation of
quantum scars [9, 10], which has attracted significant interest in recent
years. In addition, a range of quantum effects associated with photons can
also be explored, for example, the collapse and revival phenomenon [11, 12,
13], formation of Schrödinger’s cat state [14, 15, 16, 17, 18] and non-
classical state of light [19, 20, 21, 22, 23, 24]. Current experiments have
demonstrated that coupling atomic condensates to the cavity mode can lead to
fascinating phenomena like the formation of super solid phase [25, 26] and
non-equilibrium transition [27]. It is also important to note that photon loss
and other natural processes are inherent in the above mentioned systems [28,
29], which give rise to various dissipative effects for sufficiently strong
dissipation [30, 31, 32, 33, 34]. In a single cavity, the atom-photon
interacting systems within a specific regime can be effectively described by
the Jaynes-Cummings [35] or Tavis-Cummings model [36], depending on the number
of atoms in it. Moreover, coupling the cavities in an array opens up the
possibility to explore many body physics with light matter interacting systems
[1, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50], similar to the
Hubbard model. A variety of these models can exhibit quantum phase
transitions, which have been explored theoretically [51, 52, 48, 53, 54, 55,
56]. The simplest configuration of such a many body system is the dimer of two
coupled cavities forming a Jaynes-Cummings Josephson junction (JCJJ), which
has been realized in circuit QED setup [57]. This system can serve as a test
bed to study various non-equilibrium phenomena [57, 58, 59].

In the present work, we investigate the non-equilibrium dynamics and the
various associated quantum phenomena in an atom-photon interacting system
described by JCJJ in the presence of Kerr nonlinearity [60, 61, 62, 63, 64,
65, 66, 67, 68]. An insight into the overall dynamical behavior can be gained
from the semiclassical analysis, which is also useful for identifying a
variety of photonic Josephson oscillations in the JCJJ and transitions between
them. Interestingly, this system exhibits a self-trapping phenomenon, for
which photons are dynamically localized in one of the cavities [57, 58, 69,
59]. Apart from this, other self trapped states also appear as a consequence
of Kerr nonlinearity, which we analyze in detail, focusing on their dynamical
origin and regime of stability. On the other hand, in quantum dynamics, atoms
and photons become entangled, which gives rise to interesting quantum effects,
leading to the deviation from classical behavior. Additionally, the photon
field can lose its coherence as a result of phase fluctuation during the time
evolution. It is a pertinent issue to study the entanglement dynamics and
change in the state of photons due to the combined effect of interaction and
entanglement for different dynamical states, as well as for a rapid quench to
a dynamically unstable regime. We also demonstrate how the self-trapping
phenomena can be employed to control the photon mediated correlation between
the atomic qubits, which are otherwise non interacting. Such dynamical
manipulation of entanglement between the qubits in the Jaynes-Cummings dimer
model can have potential applications in quantum information processing.

The paper is organized as follows. In Sec.II, we describe the JCJJ model and
analyze it semiclassically in Sec.III to obtain different branches of
Josephson dynamics, their stability as well as transitions between them.
Quantum dynamics and its comparison with semiclassical steady states are
presented in Sec.IV. Sec.V contains a detailed discussion on the quantum
nature of the photon field, particularly phase diffusion and revival
phenomena. In this section, we also investigate the entanglement properties of
the spin $1/2$ atomic qubits corresponding to the different steady states, as
well as the signature of phase space mixing of photons in the quench dynamics.
Finally, we summarize the results and conclude in Sec.VI.

##  II THE MODEL

The Jaynes-Cummings Josephson junction formed by coupling two cavities [57],
can be described by the Hamiltonian,

$\hat{\mathcal{H}}=\sum_{i}\left[\hat{\mathcal{H}}_{\rm
JC}^{(i)}+\frac{U}{2}\hat{n}_{i}(\hat{n}_{i}-1)\right]-J\big{(}\hat{a}^{\dagger}_{L}\hat{a}_{R}+h.c.\big{)}-\mu\hat{M}$
(1)

where, the site index $i=L,R$ indicates the left and right cavity, which are
coupled by the Josephson coupling $J$. Each cavity can be modeled by the
Jaynes-Cummings Hamiltonian,

$\displaystyle\hat{\mathcal{H}}_{\rm
JC}^{(i)}=\omega\hat{n}_{i}+\omega_{0}\hat{\sigma}^{+}_{i}\hat{\sigma}_{i}^{-}+g\big{(}\hat{a}_{i}\hat{\sigma}_{i}^{+}+\hat{a}^{\dagger}_{i}\hat{\sigma}_{i}^{-}\big{)},$
(2)

describing the interaction between an atom and single mode cavity field with
frequency $\omega$, represented by the annihilation (creation) operators
$\hat{a}_{i}$($\hat{a}^{\dagger}_{i}$). The two level atom with energy gap
$\omega_{0}$ at each cavity is described by the Pauli spin operators
$\vec{\hat{\sigma}}_{i}$. The last term of $\hat{\mathcal{H}}_{\rm JC}^{(i)}$
describes the atom-photon interaction with strength $g$. In addition, we
consider the effect of Kerr nonlinearity [60, 61, 62, 63, 64, 65, 66, 67, 68]
in each cavity, represented by the second term in Eq.(1), giving rise to the
repulsive interaction of the photon field with strength $U$. The JCJJ
described by the Hamiltonian in Eq.(1) preserves the U(1) symmetry similar to
the Jaynes-Cummings model [60], leading to the conserved total excitation
number,

$\displaystyle\hat{M}=\sum_{i=L,R}(\hat{n}_{i}+\hat{\sigma}^{+}_{i}\hat{\sigma}_{i}^{-}).$
(3)

In the grand canonical ensemble, the $\mu$ in Eq.(1) represents the chemical
potential corresponding to the number of excitations. Such a JCJJ has been
realized in circuit QED setup [57], where the strength of interactions and the
photon hopping amplitude can be tuned. Moving forward, we will discuss the
different Josephson oscillations of the JCJJ, described by the Hamiltonian in
Eq.(1) within the semiclassical method and compare them with the quantum
mechanical dynamics. Throughout the paper, we use units such that
$\hbar,k_{B}=1$ and have scaled the energy (time) by $J$ $(1/J)$.

##  III SEMICLASSICAL ANALYSIS

In this section, we study the dynamics of the JCJJ governed by the Hamiltonian
given in Eq.(1), using the time dependent variational method [70]. The photons
and two level atoms in the cavities can be described semiclassically by their
respective coherent states [71], using which we construct the following time
dependent variational wavefunction,

$\displaystyle\ket{\psi_{c}(t)}=\prod_{i=L,R}\ket{\alpha_{i}(t)}\otimes\ket{\theta_{i}(t),\phi_{i}(t)}.$
(4)

The coherent state of the cavity mode is given by,

$\displaystyle\ket{\alpha_{i}}=\exp(\alpha_{i}\hat{a}_{i}^{\dagger}-\alpha_{i}^{*}\hat{a}_{i})\ket{0}$
(5)

where $\alpha_{i}$ is the eigenvalue of $\hat{a}_{i}$, representing the photon
field classically. The wavefunction for the two level atoms can be expressed
as follows,

$\displaystyle\ket{\theta_{i},\phi_{i}}=\cos(\theta_{i}/2)\ket{\uparrow}+\sin(\theta_{i}/2)e^{i\phi_{i}}\ket{\downarrow},$
(6)

where $\ket{\downarrow}$($\ket{\uparrow}$) represents the ground (excited)
state and the canonically conjugate variables $\phi_{i},z_{i}=\cos\theta_{i}$
describe the orientation of such a spin 1/2 system on the Bloch sphere, for
which
$\langle\vec{\hat{S}}_{i}\rangle=S(\sin\theta_{i}\cos\phi_{i},\sin\theta_{i}\sin\phi_{i},\cos\theta_{i})$
with $S=1/2$. The coherent state representation of the photon field is
appropriate for a large number of photons in each cavity, giving rise to the
substantial number of conserved total excitations, that can be written
semiclassically as $M=\sum_{i}|\alpha_{i}|^{2}+(1+z_{i})/2$. It is evident
from the conservation equation that the amplitude of the classical field
$\alpha_{i}$ scales with $\sqrt{M}$. Therefore, for a large number of
conserved excitations, we define
$\alpha_{i}/\sqrt{M}=\sqrt{n_{i}}\exp(\iota\psi_{i})=(x_{i}+\iota
p_{i})/\sqrt{2}$, where, $n_{i}\in[0,1]$ is the scaled photon number,
$\psi_{i}$ represents its phase, and $x_{i},p_{i}$ are the corresponding
conjugate variables. In terms of the dynamical variables
$\mathbf{x}=\\{n_{i},\psi_{i},z_{i},\phi_{i}\\}$, the Lagrangian scaled by the
total excitation number $M$ can be written as,

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\frac{1}{M}\bra{\psi_{c}}i\frac{\partial}{\partial
t}-\hat{\mathcal{H}}\ket{\psi_{c}}$ (7) $\displaystyle=$
$\displaystyle\sum_{i=L,R}\left[-\dot{\psi}_{i}n_{i}+\frac{\eta}{2}\dot{\phi}_{i}z_{i}-(\omega-\mu)n_{i}-\frac{\eta}{2}(\omega_{0}-\mu)z_{i}\right.$
$\displaystyle-\frac{\tilde{U}}{2}n_{i}^{2}-\left.\tilde{g}\sqrt{n_{i}}\sqrt{1-z_{i}^{2}}\cos(\phi_{i}+\psi_{i})\right]$
$\displaystyle+2\sqrt{n_{L}n_{R}}\cos(\psi_{L}-\psi_{R}),$

where $\eta=2S/M$ and the interaction strengths are scaled as
$\tilde{g}=g/\sqrt{M}$, $\tilde{U}=UM$. Note that, in general $\eta=2S/M$ for
a large spin system with magnitude $S$, which is considered to be small in the
present case of the Jaynes-Cummings model with $S=1/2$ and $M\gg 1$. From the
Euler-Lagrange equation
$\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{\textbf{x}}}\right)-\frac{\partial\mathcal{L}}{\partial\textbf{x}}=0$
of the dynamical variables $\textbf{x}=\\{n_{i},\psi_{i},z_{i},\phi_{i}\\}$,
we obtain the following equations of motion (EOM),

$\displaystyle\dot{n_{i}}$ $\displaystyle=$
$\displaystyle-\tilde{g}\sqrt{n_{i}}\sqrt{1-z_{i}^{2}}\sin{(\phi_{i}+\psi_{i})}$
(8a)
$\displaystyle+2\sqrt{n_{i}n_{\overline{i}}}\sin{(\psi_{i}-\psi_{\overline{i}})}$
$\displaystyle\dot{\psi_{i}}$ $\displaystyle=$
$\displaystyle-(\omega-\mu)-\frac{\tilde{g}}{2\sqrt{n_{i}}}\sqrt{1-z_{i}^{2}}\cos{(\phi_{i}+\psi_{i})}$
(8b)
$\displaystyle+\sqrt{\frac{n_{\overline{i}}}{n_{i}}}\cos{(\psi_{i}-\psi_{\overline{i}})}-\tilde{U}n_{i}$
$\displaystyle\eta\dot{\phi_{i}}$ $\displaystyle=$
$\displaystyle\eta(\omega_{0}-\mu)-\frac{2\tilde{g}z_{i}}{\sqrt{1-z_{i}^{2}}}\sqrt{n_{i}}\cos{(\phi_{i}+\psi_{i})}$
(8c) $\displaystyle\eta\dot{z_{i}}$ $\displaystyle=$ $\displaystyle
2\tilde{g}\sqrt{n_{i}}\sqrt{1-z_{i}^{2}}\sin{(\phi_{i}+\psi_{i})}$ (8d)

where $\bar{i}\neq i$. Conservation of the total excitation number yields the
constraint,

$\displaystyle n_{L}+n_{R}+\frac{\eta}{2}\big{(}z_{L}+z_{R}+2\big{)}=1.$ (9)

We solve Eq.(8) within the grand canonical ensemble, where $\mu$ is fixed by
the Eq.(9). However, in the limit $\tilde{g}\rightarrow 0$, both the photon
number and atomic inversion become conserved individually and therefore, our
formalism can not be continued to this limit. Hence, we exclude the regime of
small $\tilde{g}$ from our discussion. First, we investigate the steady states
corresponding to the fixed point (FP)
$\mathbf{x}^{*}=(n_{i}^{*},\psi_{i}^{*},z_{i}^{*},\phi_{i}^{*})$ of the EOM
given in Eq.(8), for which $\dot{\mathbf{x}}=0$. Next, we perform the linear
stability analysis around the steady states, describing the evolution of small
initial fluctuation $\delta\mathbf{x}(0)$ in the form
$\delta\mathbf{x}(t)=\delta\mathbf{x}(0)e^{i\tilde{\omega}t}$ and determine
the frequency $\tilde{\omega}$. The stability of FPs is ensured if the $\rm
Im(\tilde{\omega})=0$ and the $\tilde{\omega}$ yields the small amplitude
oscillation frequency around the corresponding steady states. In the JCJJ, the
stable steady states describe the different types of photonic Josephson
oscillations with the frequency that can be obtained from the linear stability
analysis mentioned above. Next, we find the different possible steady states
from Eq.(8) and analyze their stability.

###  III.1 STEADY STATE ANALYSIS

In this subsection, we systematically investigate various steady states
obtained from the EOM in Eq.(8) and analyze their stability as outlined above.
As evident from Eq.(8)(a,d), the steady states satisfy the conditions,

$\displaystyle\sin(\phi^{*}_{i}+\psi^{*}_{i})$ $\displaystyle=$ $\displaystyle
0$ (10a) $\displaystyle\sin(\psi^{*}_{L}-\psi^{*}_{R})$ $\displaystyle=$
$\displaystyle 0,$ (10b)

which correspond to the phase relations $\phi^{*}_{i}+\psi^{*}_{i}=0,\pi$ and
$\psi^{*}_{L}-\psi^{*}_{R}=0,\pi$, that is used to classify the steady states.
The relative phase of bosons $\psi^{*}_{L}-\psi^{*}_{R}=0(\pi)$ equivalently
describes the (anti)ferromagnetic spin configuration of the cavities in the
$S_{x}$-$S_{y}$ plane, corresponding to $\phi_{L}^{*}-\phi_{R}^{*}=0(\pi)$. We
categorize the steady states in these two classes, which are represented
schematically in Fig.1(a,b). Note that the transformations
$\phi^{*}_{i}\rightarrow\phi^{*}_{i}+\delta$ and
$\psi^{*}_{i}\rightarrow\psi^{*}_{i}-\delta$ leave the steady state equations
Eq.(8) invariant as a consequence of the U(1) symmetry. This results in a
continuous set of FPs lying on circles in the $x_{i}$-$p_{i}$ and
$S_{ix}$-$S_{iy}$ planes with corresponding radius $\sqrt{2n^{*}_{i}}$ and
$\sqrt{1-z_{i}^{*2}}/2$, respectively (see Fig.1(a,b)).

Figure 1: Schematic diagram illustrating the spin orientations of the various
steady states in the $S_{x}$-$S_{y}$ plane: (a) depicts the steady states
corresponding to the ferromagnetic class, and (b) represents the steady states
corresponding to the antiferromagnetic class.

For a particular class of spin configuration and a given value of $\eta$, the
steady states can be obtained in terms of $\\{n_{i}^{*},z_{i}^{*}\\}$, by
solving Eq.(8)(b,c), subjected to the constraint in Eq.(9), which conserves
the total excitation. The steady states thus obtained, can be categorized in
terms of the relative photon population $f=n^{*}_{R}/n^{*}_{L}$, which we
denote as $symmetric$ ($f=1$) and $self\,\,trapped$ ($f\neq 1$), corresponding
to equal and unequal photon population in the cavities. Note that, once the
photon population $n_{i}^{*}$ is obtained, it also determines the atomic
inversion $z_{i}^{*}$,

$\displaystyle
z^{*}_{i}=\frac{\xi_{2}\eta(\omega_{0}-\mu)}{\sqrt{\eta^{2}(\omega_{0}-\mu)^{2}+4\tilde{g}^{2}n^{*}_{i}}}.$
(11)

It is important to mention that the Hamiltonian given in Eq.(1) remains
invariant under the exchange of the degrees of freedom between the two
cavities ($\hat{a}_{L}\leftrightarrow\hat{a}_{R}$ and
$\hat{\vec{\sigma}}_{L}\leftrightarrow\hat{\vec{\sigma}}_{R}$), which
indicates the discrete left-right symmetry between the cavities. The
spontaneous breaking of this symmetry can give rise to the self trapped state.
Next, we analyze the steady state equations graphically, which provides a
physical picture and qualitative behavior of the steady states as well as the
transitions [72, 73] between them. For small values of $\eta$, from Eq.(9),
the total photon number can be approximately written as
$n_{L}^{*}+n_{R}^{*}=1-\eta$, which yields,

$\displaystyle n_{L}^{*}=\frac{(1-\eta)}{1+f},\quad
n_{R}^{*}=\frac{f(1-\eta)}{1+f}.$ (12)

Using these relations, the steady state equations Eq.(8)(b,c) can be reduced
to a single effective equation in terms of the relative photon population $f$,

$\displaystyle\mathcal{Y}(f)=\,\,$
$\displaystyle\xi_{1}(f-1)-\tilde{U}(1-\eta)\left(\frac{1-f}{f+1}\right)\sqrt{f}$
$\displaystyle-\,\xi_{2}\tilde{g}^{2}\sqrt{f}\left(\frac{1}{\sqrt{\mathcal{F}_{L}(f)}}-\frac{1}{\sqrt{\mathcal{F}_{R}(f)}}\right)=0,$
(13)

where,
$\mathcal{F}_{i}(f)=\eta^{2}(\omega_{0}-\mu)^{2}+4\tilde{g}^{2}n_{i}^{*}(f)$
and the discrete variable $\xi_{1}=\cos(\psi_{L}^{*}-\psi_{R}^{*})=\pm 1$
describes the spin orientation of two qubits while
$\xi_{2}=\cos(\phi^{*}_{i}+\psi_{i}^{*})=\pm 1$. For small $\eta$, the
chemical potential $\mu$ is given by,

$\displaystyle\mu=\omega-\frac{\xi_{1}}{2}\left(\sqrt{f}+\frac{1}{\sqrt{f}}\right)+\frac{\tilde{U}}{2}+\frac{\xi_{2}\tilde{g}}{4}\sqrt{1+f}\left(1+\frac{1}{\sqrt{f}}\right).$
(14)

Note that, as a consequence of the left-right symmetry between the cavities,
Eq.(13) remains invariant under the transformation $f\rightarrow 1/f$, hence
we only consider the steady state solutions for $f\in[0,1]$. The roots of
Eq.(13) yield the possible steady states for a given combination of
$\xi_{1},\xi_{2}$, which we discuss below.

####  III.1.1 Ferromagnetic class $(\boldsymbol{\psi_{L}^{*}-\psi_{R}^{*}=0})$

For the ferromagnetic orientation of the qubits, $\xi_{1}=+1$ while the other
variable can take two values $\xi_{2}=\pm 1$. When $\xi_{2}=-1$, the equation
$\mathcal{Y}(f)$ has only one root for $f=1$, describing a symmetric steady
state corresponding to the ground state configuration (Gs). On the other hand,
$\xi_{2}=1$ is a more interesting scenario since it gives rise to various non
trivial steady states, as shown in Fig.2(a,b). Similar to the previous case,
$f=1$ is always a solution of equation $\mathcal{Y}(f)$ describing a symmetric
state with higher energy density (scaled by the total number of excitation),
which is denoted by FP-F. Interestingly, two new solutions appear above a
critical coupling strength,

$\displaystyle\tilde{g}_{c1}(\tilde{U})=2+3\left(\frac{1+\tilde{U}}{2}\right)^{4/3}\eta^{2/3}-\frac{\eta}{2},$
(15)

giving rise to two self trapped states, one of which is unstable, as seen from
Fig.2(b). The FP with vanishingly small relative photon population $f\approx
0$ corresponds to a stable perfect self trapped (PST) state [57, 58],
describing a situation where almost all the photons are localized in one of
the cavities. As illustrated in Fig.2(a,b), such self trapped states arise as
a result of a saddle-node bifurcation occurring at
$\tilde{g}_{c1}(\tilde{U})$, for which non vanishing small parameter $\eta$
plays a crucial role. The unstable self trapped state STu (with larger value
of $f$) undergoes a subcritical pitchfork bifurcation with FP-F at the
critical point $\tilde{g}_{c2}(\tilde{U})$, which can be approximately written
as,

$\displaystyle\tilde{g}_{c2}(\tilde{U})=\left(\sqrt{8}+\sqrt{2}\tilde{U}\right)-\left(\frac{3}{\sqrt{2}}\tilde{U}+\sqrt{2}\right)\eta,$
(16)

after which the symmetric state FP-F becomes unstable, as depicted in
Fig.2(b). Now we focus on the steady states corresponding to the anti-
ferromagnetic spin configuration.

Figure 2: Graphical analysis of the function $\mathcal{Y}(f)$ to identify
different steady states and their bifurcations: (a) Identification of the
steady states from the roots of the function $\mathcal{Y}(f)$ and (b) the
steady state solutions in terms of $f$ (see the text for details) are shown as
a function of $\tilde{g}$ for $\tilde{U}=3$, corresponding to the
ferromagnetic class ($\xi_{1}=1,\xi_{2}=1$). (c,e) Graphical representation of
$\mathcal{Y}(f)$ and (d,f) bifurcation diagram of the FPs in terms of photon
population imbalance $Z_{p}$ for anti-ferromagnetic class ($\xi_{1}=-1$)
corresponding to $\xi_{2}=\pm 1$, respectively (see the main text), at
$\tilde{g}=1$ for different Kerr nonlinearity $\tilde{U}$. The different
markers in (a,c,e) indicate the various steady states corresponding to the
respective classes. Solid (dashed) lines in (b,d,f) represent stable
(unstable) steady states. Inset of (d,f) shows the behavior of $Z_{p}$ as a
function of $\tilde{g}$ for $\tilde{U}=4,1.5$ respectively. Here and in rest
of the figures, we consider excitation number $M=30$, equivalently,
$\eta=1/M$, unless otherwise mentioned. The energy $E$ and interaction
strengths $\tilde{g},\tilde{U}$ are measured in the units of hopping amplitude
$J$. We set $\hbar,k_{B}=1$ and $\omega=\omega_{0}=2$ for all figures.

####  III.1.2 Anti-ferromagnetic class
$(\boldsymbol{\psi_{L}^{*}-\psi_{R}^{*}=\pi})$

For the anti-ferromagnetic class with $\xi_{1}=-1$, the steady states at
higher energies compared to the ground state and the transitions between them
are very intriguing where the Kerr nonlinearity $\tilde{U}$ plays a crucial
role. The other variable can take two values $\xi_{2}=\pm 1$, and we discuss
the corresponding steady states one by one.

$\boldsymbol{\xi_{2}=1:}$ In this case, a symmetric steady state denoted by
FP-$\pi$ exists, which undergoes a pitchfork bifurcation at a critical Kerr
nonlinearity, as evident from Fig.2(c,d). After the bifurcation, FP-$\pi$
becomes unstable, giving rise to a stable self trapped state ST1. This
bifurcation is indicated by the solid black line in Fig.3(b). This phenomenon
also occurs in the Bose-Josephson junction, in absence of coupling to the spin
($\tilde{g}=0$) [74, 75, 76, 77], which has been detected experimentally [79,
78]. However, in the present case, the critical Kerr nonlinearity also depends
on the coupling strength $\tilde{g}$, which is given by,

$\displaystyle\tilde{U}_{c1}(\tilde{g})=2+\frac{\tilde{g}}{\sqrt{2}}+\left(2+\frac{3\tilde{g}}{2\sqrt{2}}\right)\eta,$
(17)

for small $\eta$. Unlike the perfect self-trapping, the relative photon
population imbalance between the two cavities,

$\displaystyle Z_{p}=\frac{n_{L}-n_{R}}{n_{L}+n_{R}}$ (18)

of the self trapped state ST1 increases continuously after the bifurcation and
approaches to unity with increasing Kerr nonlinearity $\tilde{U}$, which is
shown in Fig.2(d). In contrast, the relative population imbalance $Z_{p}$
decreases with increasing atom-photon coupling strength $\tilde{g}$ (see the
inset of Fig.2(d)), which serves as a characteristic feature of this ST1 state
for its identification. It is evident from Eq.(11), for the self trapped
state, the photon population imbalance $Z_{p}$ leads to the atomic population
imbalance,

$\displaystyle Z_{a}=\frac{|z_{R}-z_{L}|}{2},$ (19)

exhibiting similar behavior with coupling strength.

$\boldsymbol{\xi_{2}=-1:}$ A similar type of phenomenon can also be observed
for $\xi_{2}=-1$. In this case, another symmetric state FP-AF exists, which is
energetically different from FP-$\pi$ but corresponds to the same anti-
ferromagnetic spin orientation. The symmetric state FP-AF undergoes a
pitchfork bifurcation at a critical strength of Kerr interaction,

$\displaystyle\tilde{U}_{c2}(\tilde{g})=2-\frac{\tilde{g}}{\sqrt{2}}+\left(2-\frac{3\tilde{g}}{2\sqrt{2}}\right)\eta,$
(20)

which occurs only for $\tilde{g}\lesssim 2$. Above this critical coupling, the
FP-AF state becomes unstable, giving rise to a new self trapped state denoted
by ST2 (see Fig.2(e,f)). Unlike ST1, this self trapped state loses its
stability above a critical Kerr nonlinearity $\tilde{U}_{\rm I}(\tilde{g})$
(denoted by the black dashed line in Fig.3(b)) without forming a new steady
state, as shown in the bifurcation diagram given in Fig.2(f) for a fixed
coupling $\tilde{g}$. Thus ST2 can exist as a stable state only in the range
$\tilde{U}_{c2}(\tilde{g})\leq\tilde{U}<\tilde{U}_{\rm I}(\tilde{g})$ and for
$\tilde{g}\lesssim 2$. The relative photon population imbalance $Z_{p}$ for
ST2 state increases and approaches unity as both the interaction strengths
$\tilde{U},\tilde{g}$ increases (as depicted in Fig.2(f)). This behavior is
strikingly different from that of the ST1, where $Z_{p}$ diminishes with
$\tilde{g}$. Such qualitatively different features can be employed to
distinguish between the two self trapped states ST1 and ST2 during quantum
dynamics, which we will discuss in the next section. Note that, in addition to
ST1 and ST2, other self trapped states can also appear, exhibiting complicated
scenarios, which we prefer to leave out from the present discussion as they
are less relevant due to their existence within a small range of parameters.
Moreover, the signatures of these states have not been found in quantum
dynamics.

In the limit $\tilde{g}\rightarrow 0$, the steady states corresponding to
$\xi_{2}=\pm 1$ become almost identical (see Eq.(13)), with a small difference
of the order of $\eta$ in the physical quantities. In this regime, both the
self trapped states ST1 and ST2 become practically identical. However, we
exclude the small $\tilde{g}$ regime from our discussion, as the formalism can
not be extrapolated to $\tilde{g}=0$, for which both the atomic excitation and
photon number are conserved separately.

Figure 3: Steady state phase diagram as a function of interaction strengths
$\tilde{U}$ and $\tilde{g}$: (a) Ferromagnetic and (b) anti-ferromagnetic
class. Here, all the phase boundaries are obtained numerically. The boundary
of PST (dashed line) and instability line of FP-F (dashed-dotted line) in (a)
are described approximately by $\tilde{g}_{c1}$ and $\tilde{g}_{c2}$,
respectively, see Eq.(15,16) of the main text. The bifurcation lines between
FP-AF, ST2 (dashed-dotted line) and FP-$\pi$, ST1 (solid line) in (b) can be
written approximately using $\tilde{U}_{c1}$ and $\tilde{U}_{c2}$, as given by
Eq.(17,20) in the main text. The ST2 state becomes unstable outside the dashed
line ($\tilde{U}_{\rm I}(\tilde{g})$) without giving rise to any new steady
state. Note that, the small $\tilde{g}$ regime is kept blank, since our
formalism can not be continued to $\tilde{g}=0$ (see the text for details).

The plethora of steady states obtained from the above analysis are summarized
in the phase diagrams, depicted in Fig.3(a,b), separately for ferromagnetic
and anti-ferromagnetic classes, indicating their region of stability. Here,
the phase diagrams are obtained by solving the steady state equations exactly
for a fixed value of $\eta$ (equivalently, a fixed number of excitations $M$).
The numerically obtained phase boundaries of the steady states PST, FP-F,
(shown in Fig.3(a)) and the transition lines between FP-$\pi$ to ST1 as well
as FP-AF to ST2, (depicted in Fig.3(b)) are in good agreement with the
analytical results given in Eq.(15,16,17,20) for small values of $\eta$. The
appropriate parameter regimes can be identified from the phase diagrams for
observation of different dynamical behavior and transitions.

###  III.2 CLASSICAL DYNAMICS

To this end, we investigate the classical dynamics corresponding to the
different steady states illustrated in the phase diagram of Fig.3, which
provides useful information about various photonic Josephson oscillations and
transitions between them. The time evolution is performed by solving the EOM
given in Eq.(8) numerically for an appropriately chosen initial condition. In
general, if the initial condition is chosen close to a stable fixed point, the
photon number and other physical quantities oscillate around the steady state,
with oscillation frequencies obtained from the linear stability analysis. We
illustrate the oscillation around the symmetric state FP-$\pi$ by computing
the deviation of photon number $\delta n_{i}(t)=n_{i}(t)-n_{i}^{*}$ and atomic
inversion $\delta z_{i}(t)=z_{i}(t)-z_{i}^{*}$ from the corresponding steady
state values, which exhibits small amplitude oscillation around zero, shown in
Fig.4(a,b).

Figure 4: Small amplitude oscillations of the different physical quantities
and their characteristic frequencies corresponding to FP-$\pi$ state: Dynamics
of (a) photon number $\delta n_{L}(t)=n_{L}(t)-n_{L}^{*}$, (b) atomic
inversion $\delta z_{L}(t)=z_{L}(t)-z_{L}^{*}$, subtracted from their steady
state values. The respective Fourier spectra,
$\delta\tilde{n}_{L}(\tilde{\omega})/\delta\tilde{n}_{L}(\tilde{\omega}_{0})$
and
$\delta\tilde{z}_{L}(\tilde{\omega})/\delta\tilde{z}_{L}(\tilde{\omega}_{0})$,
scaled by the amplitude of the lowest frequency are shown in (c) and (d). The
lowest frequency $\tilde{\omega}_{0}$ obtained from the linear stability
analysis is present in both the quantities. Here and in the remaining figures,
time $t$ and frequencies $\tilde{\omega},\tilde{\omega}_{0}$ are scaled by
$1/J$ and $J$ respectively.

Numerically, the Fourier transform of the time evolution of photon population
and atomic inversion yields the relevant frequencies present in the dynamics.
As observed from Fig.4(c,d), the lowest frequency $\tilde{\omega}_{0}$
obtained from the linear stability analysis of the steady state FP-$\pi$
corresponds to the highest amplitude of the Fourier transform, indicating its
dominant role in both the photon and spin (atom) dynamics. However, as evident
from Fig.4(d), the higher frequency modes also contribute to the spin degree
with small amplitude, resulting in fast dynamics (see Fig.4(b)). Such dynamics
around a stable fixed point as depicted in Fig.4(a,b), usually lie on an
invariant KAM torus [80] with off resonant frequencies. Consequently, a
trajectory over a long duration of time densely fills the region around the
ring of FPs.

Figure 5: Classical phase portrait and variation of oscillation frequency
across the transition between FP-$\pi$ and ST1 states: Dynamics of photon
field in $x_{L}$-$p_{L}$ plane of left cavity corresponding to the (a) stable
regime with $(\tilde{g},\tilde{U})=(1.0,1.5)$ and (b) unstable regime with
$(\tilde{g},\tilde{U})=(1.0,6.5)$, of FP-$\pi$ state. The red (blue)
trajectory in (b) is obtained by starting the dynamics near the stable ST1
(unstable FP-$\pi$) state. (c) Time evolution of the relative phase of photon
$\psi_{r}=\psi_{L}-\psi_{R}$ for stable FP-$\pi$ state. (d) Dynamics of the
relative photon population imbalance $Z_{p}$ in the unstable regime of
FP-$\pi$. The colored lines carry the same meaning as in (b). (e) Variation of
the lowest frequency $\tilde{\omega}_{0}$ (of small amplitude oscillation)
with $\tilde{U}$ across the transition, for $\tilde{g}=1$. The black circles
and the solid lines denote the lowest frequency obtained from the Fourier
transform of the dynamics, and the linear stability analysis, respectively.

It is fascinating to study the dynamics across the bifurcation of the steady
states, particularly the emergence of the self trapped states. Here, we focus
on the classical dynamics across the pitchfork bifurcation of the symmetric
state FP-$\pi$ to the self trapped state ST1, which occurs by tuning the Kerr
nonlinearity $\tilde{U}$. Before the bifurcation, since the stable FP-$\pi$ is
a symmetric state, we study the dynamics of the photon field in $x$-$p$ plane
for one of the cavities, depicted in Fig.5(a). As mentioned before, due to the
U(1) symmetry, the continuous FPs lie on a circle in the $x$-$p$ plane of the
photon field (black line in Fig.5(a)). Ideally, the small amplitude dynamics
is expected to be confined around one of the FPs, which only occurs if we
consider the initial condition for $\phi+\psi$ to be the same as the value of
the FPs, without any fluctuation around it. However, for an arbitrary initial
condition around one of the FPs, the trajectory surrounds all the fixed points
on the ring, as depicted in Fig.5(a). As the main characteristic feature of
the FP-$\pi$ mode, the relative phase of photons $\psi_{r}=\psi_{L}-\psi_{R}$
oscillates around the value $\pi$, which is shown in Fig.5(c). Above the
critical coupling $\tilde{U}_{c1}$, FP-$\pi$ becomes unstable, and depending
on the initial condition, the dynamics is attracted towards one of the stable
self trapped states. It is evident from Fig.5(b) that the trajectory is
repelled from the FP-$\pi$ state and attracted towards the ring of FPs
corresponding to the ST1 state. Consequently, the photon imbalance $Z_{p}$
oscillates around a finite value corresponding to the steady state (see
Fig.5(d)). The signature of this transition can be observed from the
oscillation frequencies of FP-$\pi$ and ST1 state, both of which vanish at the
critical coupling strength $\tilde{U}_{c1}$, as evident from Fig.5(e). A
similar phenomenon also occurs for the bifurcation of FP-AF to ST2 state. Such
behavior is similar to the mode softening phenomenon associated with the
quantum phase transition [81].

Next, we focus on the dynamics of the self trapped states ST1 and ST2. As a
distinguishing feature between them, the relative photon population imbalance
$Z_{p}$ decreases with increasing atom-photon coupling $\tilde{g}$ for ST1
(see Fig.6(a)) whereas it increases for ST2, as depicted in Fig.6(b).

Figure 6: Classical dynamics around different self trapped states: (a,b)
Relative photon population imbalance $Z_{p}$ and (c,d) atomic imbalance
$Z_{a}$ for the ST1 and ST2 states, respectively, at different coupling
strengths $\tilde{g}$. The Kerr interaction strength for the dynamics of ST1
and ST2 are $\tilde{U}=5$ and $\tilde{U}=2$, respectively. The horizontal
lines with different line styles indicate the classical steady-state values of
both $Z_{p}$ as well as $Z_{a}$.

Since the atomic inversion is directly related to the photon population in
each cavity, as given in Eq.(11), the relative photon population imbalance
$Z_{p}$ can also induce an atomic inversion imbalance, $Z_{a}$ for the self
trapped states. The variation of $Z_{a}$ with $\tilde{g}$ can also distinguish
between the two self trapped states ST1 and ST2, exhibiting opposite behavior,
which is illustrated in Fig.6(c,d). However, its variation is small for the
ST2 state as compared to that of ST1.

So far, we have analyzed the classical dynamics based on a simplified
description, neglecting the atom-photon correlation. Hence, it is important to
investigate the signature of such dynamical states in quantum dynamics and the
effect of atom-photon entanglement, which we consider in the next sections.

##  IV QUANTUM DYNAMICS

In this section, we study the full quantum dynamics of the JCJJ and compare
them with the classical dynamics to investigate the effect of Kerr
nonlinearity as well as the atom-photon correlation. We evolve the initial
state $\ket{\Psi(0)}$, with a fixed number of excitations $M$, within the
Schrödinger prescription, which is performed numerically by truncating the
basis up to a sufficiently large number $N_{\rm max}$. To compare with
classical dynamics, we choose the initial state as the product of coherent
states of photons and spins, described in Eq.(5,6) respectively which
represents the classical phase space point. To investigate the signature of
different branches of the dynamical states, we time evolve the appropriately
chosen initial state and obtain the dynamics of different physical quantities
such as the population of photons and the atom in different cavities as well
their imbalance, characterizing those states. The parameters are also chosen
from the stability region of the corresponding states from the phase diagram,
given in Fig.3.

First, we study the dynamics of the symmetric states FP-F and FP-AF
corresponding to the ferromagnetic and the anti-ferromagnetic classes
respectively. To characterize these states quantum mechanically, we obtain the
photon population imbalance,
$Z_{p}=(\langle\hat{n}_{L}\rangle-\langle\hat{n}_{R}\rangle)/(\langle\hat{n}_{L}\rangle+\langle\hat{n}_{R}\rangle)$
where $\langle\hat{n}_{i}\rangle$ is computed from the time evolved state
$\ket{\Psi(t)}$ starting from the initial coherent state. For both FP-F and
FP-AF states, we obtain the time evolution of $Z_{p}$ and compare them with
that obtained from the classical dynamics, as shown in Fig.7(a,b). It is clear
from Fig.7(a,b) that the simple classical analysis is able to capture the full
quantum dynamics reasonably well, however, there are certain deviations as $t$
increases.

Figure 7: Comparison between quantum dynamics of the symmetric states FP-F and
FP-AF: The small amplitude oscillation of the photon imbalance $Z_{p}$ around
(a) FP-F, (b) FP-AF, are compared with their respective classical dynamics
(black dashed lines). The spin orientations of these states are described by
$C_{LR}$ in (c,d) respectively. The spin trajectory in the
$\langle\hat{S}_{Lx}\rangle$-$\langle\hat{S}_{Ly}\rangle$ plane of the left
cavity is shown for (e) FP-F and (f) FP-AF state. The black dashed circles in
(e,f) represent the ring of classical FPs. The evolution of spins around the
FP-AF state exhibits dephasing phenomenon. Parameter chosen:
($\tilde{g},\tilde{U}$)=($0.5,0.5$).

To reveal the relative spin orientation in the two cavities, we introduce the
quantity

$\displaystyle
C_{LR}=\frac{\langle\hat{S}_{Lx}\hat{S}_{Rx}+\hat{S}_{Ly}\hat{S}_{Ry}\rangle}{\sqrt{\left(\frac{1}{4}-\langle\hat{S}_{Lz}\rangle^{2}\right)\left(\frac{1}{4}-\langle\hat{S}_{Rz}\rangle^{2}\right)}},$
(21)

which in the classical limit takes the value -1(+1) corresponding to the
(anti)ferromagnetic class of steady states. As shown in Fig.7(c,d), the
quantum dynamics of $C_{LR}$ also approaches these values for FP-F and FP-AF,
which is consistent with their classification based on classical analysis. On
the other hand, in quantum dynamics, the correlation (entanglement) between
spins and photons gives rise to interesting effects leading to the deviation
from classicality. In the spin dynamics of the FP-F state, the average values
of the spin components in the $x$-$y$ plane evolve around a circle
corresponding to the classical FPs. Whereas for the FP-AF state, the spin
trajectory deviates from the ring of classical FPs and spirals to the center
corresponding to $\langle\hat{S}_{x}\rangle=\langle\hat{S}_{y}\rangle=0$,
exhibiting spin dephasing phenomena [82], as seen from Fig.7(f). Typically for
spin $1/2$ qubits, the classical description fails due to the enhanced quantum
fluctuations and entanglement with photons, which we analyze later.

Figure 8: Quantum dynamics corresponding to the perfect self trapped state
PST: (a) Time evolution of the relative photon population imbalance $Z_{p}$
for $\tilde{g}=5$ and different Kerr nonlinearity $\tilde{U}$. The imbalance
decays above a certain large value of Kerr interaction strength. (b) Variation
of the decay rate $\Gamma$ of the imbalance with $\tilde{U}$.

Next, we investigate different types of self-trapping phenomena from quantum
dynamics. We search for a perfect self trapped state, where almost all the
photons become localized in one of the cavities. It is evident from the
classical phase diagram that atom-photon interaction is crucial for perfect
self-trapping of photons. For small Kerr nonlinearity, we identify the perfect
self trapped state quantum mechanically, for which the relative imbalance of
photon $Z_{p}$ remains close to unity for a sufficiently long time (see the
red line in Fig.8(a)).

Figure 9: Comparison between the quantum dynamics of the two distinct self
trapped states ST1 and ST2: The plots (a,c,e) and (b,d,f) correspond to
dynamics of various quantities around ST1 and ST2 states for $\tilde{U}=5$ and
$\tilde{U}=2$, respectively. Time evolution of the (a,b) relative photon
imbalance $Z_{p}$ and (e,f) atomic imbalance $Z_{a}$, for different coupling
strength $\tilde{g}$. (c,d) Dynamics of the spin correlation function $C_{LR}$
at $\tilde{g}=1.2$.

However, for sufficiently large Kerr nonlinearity, the imbalance becomes
significantly lower than unity and decays with time, which is shown in
Fig.8(a). The rate of exponential decay $\Gamma$ can be obtained by
numerically fitting the time evolution of the imbalance. The variation of the
decay rate with Kerr nonlinearity exhibits an interesting feature as it grows
rapidly above certain Kerr nonlinearity, which is depicted in Fig.8(b). This
indicates that sufficiently large Kerr nonlinearity induces an instability in
perfect self-trapping. Although, classically, the stable perfect self trapped
state exists for $\tilde{g}>\tilde{g}_{c1}$, the sufficiently large Kerr
nonlinearity gives rise to instability in such a state during quantum
dynamics.

We also analyze the self trapped states ST1 and ST2 of the anti-ferromagnetic
class, shown in Fig.9(a,b). As a characteristic feature of these states, we
study the dynamics of the photon imbalance $Z_{p}$ for increasing values of
atom-photon coupling strength $\tilde{g}$. As depicted in Fig.9(a), for ST1,
the imbalance decreases with $\tilde{g}$, whereas, it increases for ST2 (see
Fig.9(b)), which is consistent with the classical analysis and can be used to
distinguish between these two self trapped states. Since both these states
belong to the anti-ferromagnetic class, during dynamics, the quantity $C_{LR}$
acquires a negative value for both states, however, its deviation from the
classical value is large for ST2 (depicted in Fig.9(c,d)). We also study the
dynamics of the atomic imbalance
$Z_{a}=|\langle\hat{S}_{Lz}\rangle-\langle\hat{S}_{Rz}\rangle|$, which
decreases with increasing coupling strength $\tilde{g}$ for ST1, as shown in
Fig.9(e), that is in agreement with the classical analysis (see Fig.6). For
the ST2 state, the evolution of $Z_{a}$ always saturates to a very small
value, exhibiting a weak variation with $\tilde{g}$ (see Fig.9(f)), which is
in stark contrast with the ST1 state. The above analysis reveals that the
deviation from classicality is significantly large for the ST2 state as
compared to ST1. In spite of such quantitative differences, the main features
of the semiclassical steady states are observed in the quantum dynamics.

It is worth mentioning that new effects in dynamics can arise, since
dissipation is inevitable in such systems, specifically due to photon loss and
spontaneous emission from atoms, which can however be controlled by the
appropriate experimental setup [83, 84, 85, 86]. For weak dissipation, it is
expected that in JCJJ the photonic oscillations corresponding to the symmetric
states will be damped, nevertheless, their signature can be detected from the
oscillation frequency. Similarly, the self-trapping phenomenon can still be
detected within its lifetime due to a decay in the photonic imbalance. Even in
the presence of dissipation, the self trapping phenomenon has already been
observed in different experimental setups [87, 57, 78]. It is possible to
stabilize the steady states through the incoherent pumping processes,
balancing the photon loss [34, 59].

In addition to the characteristic features of the steady states, it is crucial
to investigate other quantum effects manifested in the dynamics. These include
the generation of entanglement between the qubit and the photonic degree of
freedom, as well as the quantum state of the photon field, which we discuss in
the next section.

##  V ENTANGLEMENT AND QUANTUM FLUCTUATIONS

The semiclassical formalism presented in Sec.III is based on the product
coherent state representation, which is appropriate for describing the phase
coherent photonic Josephson dynamics. However, the presence of interactions
and Kerr nonlinearity can destroy such coherent dynamics due to enhanced phase
fluctuations, which in turn gives rise to the deviation from classicality due
to a change in the nature of the quantum state. To this end, we study the
phase fluctuations of the photon field by constructing the phase states [79,
88],

$\displaystyle\ket{\psi_{m}}=\frac{1}{\sqrt{N_{max}+1}}\sum_{n=0}^{N_{max}}\exp(in\psi_{m})\ket{n}$
(22)

with $\psi_{m}=\psi_{0}+2\pi m/(N_{max}+1)$, where $m$ is an integer
$m\in[0,N_{max}]$ and $\psi_{m}\in[-\pi,\pi]$. These phase states are
eigenstates of the phase operator $\hat{\psi}$, as given by,

$\exp(\pm i\hat{\psi})\ket{\psi_{m}}=\exp(\pm i\psi_{m})\ket{\psi_{m}}.$ (23)

The phase distribution corresponding to the photon field in one of the
cavities ($i=L,R$) is given by,

$\displaystyle{\rm P}(\psi_{m}^{i})={\rm
Tr}(\hat{\rho}^{i}_{p}\ket{\psi_{m}}\bra{\psi_{m}})$ (24)

with $\sum_{m}{\rm P}(\psi_{m}^{i})=1$, where $\hat{\rho}_{p}^{i}$ is the
reduced density matrix corresponding to the photon field of the $i$th cavity,
obtained by tracing out the other degrees. The average value and the
fluctuation of the phase of the photon field in each cavity can be computed
from the phase distribution as,

$\displaystyle\langle\hat{\psi}_{i}\rangle$ $\displaystyle=$
$\displaystyle\sum_{m}\psi_{m}{\rm P}(\psi^{i}_{m})$ (25a)
$\displaystyle(\Delta\hat{\psi}_{i})^{2}$ $\displaystyle=$
$\displaystyle\sum_{m}(\psi_{m}-\langle\hat{\psi}_{i}\rangle)^{2}\,{\rm
P}(\psi^{i}_{m}).$ (25b)

Using the above prescription, we compute the mean phase difference between the
cavity modes
$\psi_{r}=\langle\hat{\psi}_{L}\rangle-\langle\hat{\psi}_{R}\rangle$ and its
time evolution.

Figure 10: Classical-quantum correspondence for the symmetric states FP-F and
FP-$\pi$: Quantum dynamics of the small amplitude oscillation of the relative
phase $\psi_{r}=\langle\hat{\psi}_{L}\rangle-\langle\hat{\psi}_{R}\rangle$ of
the photon field corresponding to (a) FP-F and (b) FP-$\pi$ states for the
coupling strengths ($\tilde{g}$,$\tilde{U}$)=(0.5,0.5). The black dashed lines
represent the classical dynamics around these states.

The dynamics of the relative phase of the photon modes for the symmetric FP-F
and FP-$\pi$ states are shown in Fig.10(a) and (b) respectively, which exhibit
coherent oscillations around their steady state values $0$ and $\pi$.

Figure 11: Phase diffusion dynamics: Time evolution of the relative phase
fluctuation $(\Delta\psi)^{2}_{\rm N}$ for (a) symmetric state FP-$\pi$ with
$(\tilde{g},\tilde{U})=(1.0,0.5)$ and (b) self trapped state ST1 for
$(\tilde{g},\tilde{U})=(1.0,6.5)$. The horizontal dashed line indicates the
maximum value $(\Delta\psi)^{2}_{\rm N}=1$. (c,d) Snapshots of the
corresponding phase distributions at different times.

To quantify the degree of coherence, we calculate the normalized phase
fluctuation of photons $(\Delta\psi_{i})^{2}_{\rm
N}=(\Delta\psi_{i})^{2}/(\Delta\psi_{i})^{2}_{\rm max}$ in one of the cavities
($i=L,R$) where the maximum phase fluctuation $(\Delta\psi_{i})^{2}_{\rm
max}=\pi^{2}/3$ corresponds to a uniform phase distribution [79]. For the
FP-$\pi$ state, the phase fluctuation remains small during time evolution, as
illustrated in Fig.11(a), due to which the coherent phase oscillation is
retained. However, the phase fluctuation of the FP-$\pi$ state increases
slowly, and after a sufficiently long time, it approaches the maximum value.
On the other hand, an enhancement in the growth of phase fluctuation can be
observed for the self trapped state ST1, arising for large $\tilde{U}$, as
seen from Fig.11(b). In general, the phase fluctuation increases with Kerr
nonlinearity, which is evident from the above comparison. Such enhanced phase
fluctuation during the time evolution is associated with the broadening of the
phase distribution, indicating the deviation of the photon field from its
classical representation in terms of the coherent state. Spreading of the
phase distribution of FP-$\pi$ and ST1 states during time evolution is
apparent from Fig.11(c,d). Even though the phase fluctuation attains its
maximum value almost immediately for the ST1 state, the appearance of dips in
the time evolution of $(\Delta\psi_{L})^{2}_{\rm N}$, as observed from
Fig.11(b), corresponds to the revival of the phase of the photon field, which
we discuss later.

Figure 12: Entanglement entropy of the spin degree for various dynamical
states: Dynamics of the scaled entanglement entropy $\mathcal{S}_{L}$ of the
spin in the left cavity corresponding to (a) FP-F and (b) FP-AF for
$(\tilde{g},\tilde{U})=(0.5,0.5)$. Entanglement entropy of the spins in both
the cavities corresponding to (c) ST1 state for
$(\tilde{g},\tilde{U})=(0.5,5.0)$ and (d) ST2 state for
$(\tilde{g},\tilde{U})=(1.2,2.0)$. (e,f) Variation of the relative
entanglement entropy $\Delta\mathcal{S}=\mathcal{S}_{L}-\mathcal{S}_{R}$ with
coupling strength $\tilde{g}$ for ST1 and ST2 states, respectively. The Kerr
nonlinearity in (e,f) is chosen the same as that in (c,d).

Apart from the phase fluctuation, the entanglement between the photon field
and spins during the time evolution gives rise to interesting quantum effects
and deviation from classicality. Starting from the total density matrix
$\hat{\rho}=\ket{\Psi(t)}\bra{\Psi(t)}$, computed from the full wavefunction
$\ket{\Psi(t)}$, the reduced density matrix of a subsystem (such as the
spin/photon field of each cavity) can be obtained by integrating out the rest
of the degrees of freedom. Following this prescription, we compute the
entanglement entropy of the subsystem (corresponding to the cavities) as,

$\mathcal{S}_{i}=-\sum_{l}\lambda_{l}^{i}\log(\lambda_{l}^{i})$ (26)

where, $\lambda_{l}^{i}$ represents the eigenvalue with index $l$ of the
reduced density matrix corresponding to the subsystem denoted by $i$ (for
example, $i=L,R$ is the cavity index). In a similar manner, we can also
compute the reduced density matrix and entanglement entropy $\mathcal{S}_{LR}$
for the total photon and spin degree separately. Ideally, $\mathcal{S}_{i}$
vanishes for product state, but due to atom-photon interactions, the
entanglement entropy increases during time evolution. We obtain the
entanglement entropy $\mathcal{S}_{i}$ of the spin in each cavity,
corresponding to the symmetric states FP-F and FP-AF, which are compared in
Fig.12(a,b). Unlike the FP-F state, $\mathcal{S}_{i}$ grows rapidly and
saturates to its maximum value $k_{B}\ln 2$ for the FP-AF state, due to which
the spin dynamics deviates from the classical steady states, exhibiting
dephasing phenomenon (see Fig.7(f)), as discussed in Sec.IV.

Figure 13: Photon mediated correlation between the spins in the different
cavities for the ST1 state: The dynamics of (a) mutual information
$\mathcal{I}$ and (b) relative entanglement entropy $\Delta\mathcal{S}$
between the spins in the cavities for different Kerr strengths $\tilde{U}$.
The saturation value of the photon population imbalance $\overline{Z}_{p}$
obtained after long time, corresponding to these values of $\tilde{U}$ are
mentioned in the figure. Here, we set $\tilde{g}=1$.

We also compare the entanglement entropy of spins in both the cavities for the
self trapped states ST1 and ST2, which reveals contrasting features between
them. For ST2, the $\mathcal{S}_{i}$ is almost the same for both the cavities
and saturates to the maximum value. On the contrary, for ST1, the entanglement
entropy is larger, corresponding to the cavity containing more number of
photons, as seen from Fig.12(c,d). In addition, we also study the difference
between the entanglement entropy of spins in the two cavities
$\Delta\mathcal{S}=\mathcal{S}_{L}-\mathcal{S}_{R}$ and their variation with
coupling strength $\tilde{g}$, as shown in Fig.12(e,f). For ST1 state, similar
to the photon imbalance $Z_{p}$, the saturation value of $\Delta\mathcal{S}$
decreases with increasing $\tilde{g}$, which is in stark contrast to ST2
state, for which $\Delta\mathcal{S}$ vanishes over time, irrespective of the
values of $\tilde{g}$. However, the timescale at which $\Delta\mathcal{S}$
vanishes exhibits a weak variation with $\tilde{g}$ within a small range. Such
contrasting feature of entanglement dynamics of two qubits can also
distinguish the self trapped states ST1 and ST2.

Apart from the interaction induced entanglement between spins and photons in
each cavity, two apparently non interacting spins can also be entangled, which
is mediated by photons. Such photon induced hidden correlation between two
spins can be analyzed from the mutual information [89, 90, 91, 92, 93, 94],

$\mathcal{I}=\mathcal{S}_{L}+\mathcal{S}_{R}-\mathcal{S}_{LR},$ (27)

which reveals very interesting behavior for the self trapped states. For the
ST2 state, both $\Delta\mathcal{S}$ and $\mathcal{I}$ are very small,
exhibiting almost no variation with interaction strengths, which indicates
that the reduced density matrix corresponding to the two spins approaches to
the maximally mixed state [89]. On the other hand, in the case of ST1,
increasing the photon population imbalance leads to an increase in
$\Delta\mathcal{S}$, while the mutual information $\mathcal{I}$ decreases, as
evident from Fig.13.

Figure 14: Evolution of the quantum state and its signature in the phase
fluctuation and mutual information: Dynamics of (a) mutual information
$\mathcal{I}$ of the spins and (b) relative phase fluctuation
$(\Delta\psi_{L})^{2}_{\rm N}$ of the photon field, corresponding to ST1
state. (c) Husimi distribution $Q$ in the $x_{L}$-$p_{L}$ plane of the photon
field at different time instances of the revival cycle, marked in (a,b) by
different symbols. The white dashed lines in (c) represent the ring of FPs
corresponding to the ST1 state. Parameter chosen:
$(\tilde{g},\tilde{U})=(1.0,6.5)$.

Such tunability of quantum correlations between the two non interacting spins
in the cavities can have potential applications in quantum information
processing. Additionally, for the ST1 state, the mutual information
$\mathcal{I}$ exhibits dip and spike like structure during the time evolution,
as seen from Fig.14(a). Such dips in the mutual information correspond to the
phase revival phenomenon [11, 12, 13, 95], resulting in a sudden drop in phase
fluctuation, as seen in Fig.14(b). This revival cycle can be analyzed from the
evolution of the semiclassical phase space density of the photon field,
described by the Husimi distribution,

$Q(\alpha,\alpha^{*})=\frac{1}{\pi}\bra{\alpha}\hat{\rho}_{i}^{p}\ket{\alpha},$
(28)

where $\hat{\rho}_{i}^{p}$ represents the reduced density matrix of the photon
field in the cavities. In order to compare the Husimi distribution with the
corresponding classical dynamics, we introduce the scaled phase space
variables $x_{i}=(\alpha_{i}+\alpha_{i}^{*})/\sqrt{2M}$ and
$p_{i}=(\alpha_{i}-\alpha_{i}^{*})/i\sqrt{2M}$, where $M$ represents the
conserved total excitation number associated with the initial state.
Initially, the density is localized around one of the FPs, exhibiting the
coherent structure of the photon field. As time evolves, the phase space
density spreads over the ring of fixed points, describing the loss of
coherence and finally, it is reconstructed at a point in the phase space when
another dip in $(\Delta\psi)^{2}_{\rm N}$ occurs, exhibiting the revival
phenomenon (see Fig.14(c)). Interestingly, in the middle of the cycle, the
phase space density splits and becomes localized around two diagonally
opposite phase space points. Although such bimodal phase space distribution of
photons resembles that of a cat state [14, 15, 16, 17, 18, 60], they differ in
terms of their coherence property, which can be captured from the Wigner
function, as discussed in Appendix.A. Such structure of phase space density is
associated with the appearance of a spike in the mutual information, as seen
from Fig.14(a).

The above analysis elucidates fascinating quantum effects and entanglement
associated with the evolution of the quantum state, corresponding to different
dynamical branches, which can also be relevant in the context of quantum
information processing. Such quantum effects give rise to the deviation from
classicality, nevertheless, the qualitative behavior of the system can still
be captured from the coherent state description. Apart from the steady state
dynamics, it would also be interesting to investigate the evolution of the
quantum state, particularly that of the photon field, when the system is
driven to the unstable regime.

##### Quench dynamics to unstable regime:

Figure 15: Quench dynamics under sudden change in the Kerr nonlinearity
$\tilde{U}$: The initial state is prepared corresponding to the stable
FP-$\pi$ state for $\tilde{U}=0.5$ before the quench. (a,c,e) and (b,d,f)
correspond to the dynamics after a quench to the stable regime of FP-$\pi$
with $\tilde{U}=1.5$ and unstable regime at $\tilde{U}=6.5$, respectively.
(a,b) the time evolution of the scaled potential and kinetic energies, PE
$=m\omega\langle\hat{x}_{L}^{2}\rangle/2\hbar M$ (red) and KE
$=\langle\hat{p}_{L}^{2}\rangle/2\hbar m\omega M$ (blue), respectively. (c,d)
depicts the dynamics of scaled photon number variance $(\Delta
n_{L})^{2}/\langle n_{L}\rangle$. The pink dashed line in (c,d) represents
$(\Delta n_{L})^{2}/\langle n_{L}\rangle=1$, corresponding to the coherent
state. (e,f) Husimi distribution $Q$ of the photon field in the left cavity at
time $t=4$. The white dashed line in (e) represents the ring of steady states
corresponding to the FP-$\pi$ state. Parameter chosen: $\tilde{g}=0.5$.

Next, we investigate the quench dynamics corresponding to an abrupt change in
the Kerr nonlinearity $\tilde{U}$, starting from the initial coherent state
corresponding to the stable FP-$\pi$ state. First, we consider a small change
in $\tilde{U}$, for which the FP-$\pi$ state remains stable. Under this sudden
change, the system still follows the stable FP-$\pi$ branch, exhibiting
oscillations around it. After quench, the initial coherent state begins to
move around the ring of FPs, as depicted in Fig.15(e). During the time
evolution, the wavefunction initially remains fairly localized and then slowly
spreads along the ring of FPs. As a consequence, the scaled potential energy
(PE) and kinetic energy (KE)
$\left(m\omega\langle\hat{x}_{i}^{2}\rangle/2\hbar
M,\langle\hat{p}_{i}^{2}\rangle/2\hbar m\omega M\right)$ of the photon field
in each cavity oscillate coherently, keeping the average photon number fixed
(see Fig.15(a)). Moreover, the photon phase fluctuation increases slowly and
finally saturates to its maximum value after a sufficiently long time, while
its phase space density remains localized around the ring of FPs.

On the contrary, when the interaction strength $\tilde{U}$ is quenched above
the transition point where the $\pi$-mode becomes unstable, the system
exhibits incoherent dynamics, dominated by large fluctuations, instead of
following any stable branch of the self trapped state. After quenching, the
photon field loses its coherence rapidly, as the phase fluctuation attains the
maximum value (See inset of Fig.16(b)). Moreover, the scaled kinetic and
potential energies approach the same steady value, without exhibiting large
amplitude oscillations (see Fig.15(b)), analogous to the equipartitioning of
energy. In this regime, the reduced density matrix of the photon field in each
cavity exhibits a predominant contribution from the diagonal elements.
Simultaneously, the associated Husimi distribution spreads widely across the
phase space, revealing substantial fluctuations in photon number $(\Delta
n_{i})^{2}=\langle\hat{n}_{i}^{2}\rangle-\langle\hat{n}_{i}\rangle^{2}\gg\langle
n_{i}\rangle$, notably in contrast to the initial coherent state for which
$(\Delta n_{i})^{2}=\langle n_{i}\rangle$. The entanglement entropy of the
photon field grows linearly with time and finally saturates, as seen from
Fig.16(b). Moreover, the reduced density matrix of the spins in two cavities
approaches the maximally mixed state. Such scenarios of quench dynamics to an
unstable regime resemble thermalization, which led us to compare the state of
the photon field with that of thermal gas.

The density matrix of the thermal photon gas is given by [96],

$\hat{\rho}_{\rm
Th}=\sum_{n=0}^{\infty}\frac{\langle\hat{n}\rangle^{n}}{(1+\langle\hat{n}\rangle)^{1+n}}\ket{n}\bra{n},$
(29)

where $\langle n\rangle$ is the average number of photons. The corresponding
Husimi distribution of thermal photons expressed in terms of the dimensionless
classical variables $x,p$ takes the form of a symmetric Gaussian function,
which can be written as,

$\displaystyle Q_{\rm Th}=\frac{1}{\pi(1+\langle
n\rangle)}\exp\left(-\frac{r^{2}}{2(1+\langle n\rangle)}\right),$ (30)

where $r^{2}=M(x^{2}+p^{2})$. For a comparison with the thermal state, we
obtain the angular averaged Husimi function
$\overline{Q}(r)=\int_{0}^{2\pi}Q(r,\theta)d\theta$ of the photon field after
the quench dynamics. Although the spread in the phase space distribution is
comparable to that of the thermal state, their detailed structure differs,
which is illustrated in Fig.16(a). Notably, a peak appears in the phase space
distribution of the photon field near the ring of fixed points. Interestingly,
such a scenario indicates the retention of memory of the underlying steady
states, even when they become unstable after quench, which is similar to
scarring phenomena [9, 10, 97, 98]. Consequently, the entanglement entropy
$\mathcal{S}_{L}^{p}$ saturates to a slightly lower value compared to the
entropy of the thermal state, as depicted in Fig.16(b).

Figure 16: Comparison of the photon field with the thermal state, after the
quench to the unstable regime. (a) Angular averaged Husimi distribution of the
photon field in the left cavity $\overline{Q}(r_{L})$ (solid red line)
compared with the thermal distribution (black dashed line) at time $t=4$. (b)
Dynamics of the entanglement entropy $\mathcal{S}_{L}^{p}$ of the photon field
in the left cavity. The entropy of the thermal state is marked by the black
dashed line. The inset of (b) shows the time evolution of the relative phase
fluctuation $(\Delta\psi)^{2}_{\rm N}$, which saturates to unity.

The approach to a thermal state can be investigated experimentally through
non-equilibrium dynamics of photon field in JCJJ by initially preparing
identical coherent states in both cavities with a total energy corresponding
to the $\pi$ mode. A subsequent quench to the unstable domain of $\pi$ mode
can be performed by suitably choosing the strength of Kerr nonlinearity or
photon hopping. Both the cavity and circuit QED can serve as ideal platforms
to achieve such tunable parameters [37, 67]. Moreover, the optomechanical
systems can also be used to generate Kerr nonlinearity in a controlled manner
[68]. In this context, it is noteworthy that the thermalization of the
photonic gas has already been observed in experimental setups such as optical
microwave resonators [99, 20]. Apart from the coherence properties, such non-
equilibrium dynamics of the photon fields can also reveal interesting
phenomena, which can be probed in experiments [25, 27].

##  VI Conclusion

To summarize, we explore the non-equilibrium dynamics of the Jaynes-Cummings
dimer model in the presence of Kerr nonlinearity, focusing on the quantum
states of photons as well as entanglement properties corresponding to the
different dynamical states. Within the semiclassical approach, we
systematically study the dynamics to chart out a variety of steady states and
their regime of stability for different atom-photon coupling strengths and
Kerr nonlinearity. Moreover, the stability analysis yields the frequency of
photonic Josephson oscillation that can be probed in experiments. The various
branches of dynamical states are categorized based on the relative spin
orientation and photon population imbalance between the cavities. Self trapped
states with unequal photon populations in the two cavities emerge as a
consequence of the different transitions. Apart from a perfect self trapped
state, arising from a saddle-node bifurcation, we also identify two different
self trapped states for which the Kerr interaction plays an important role.
From quantum dynamics, we also observe the characteristic features of
different steady states obtained semiclassically. However, interactions and
atom-photon entanglement give rise to intriguing quantum effects, leading to a
deviation from classicality. In contrast to the classical motion, dephasing is
observed in spin dynamics as a result of relatively large quantum fluctuations
in the spin $1/2$ qubits. During the time evolution, the state of the photon
field deviates from the initial coherent state and gradually loses its
coherence due to phase fluctuation, which is typically enhanced by the Kerr
nonlinearity. Apart from the phase fluctuations, we identify a periodic
revival phenomenon for the self trapped state, exhibiting fascinating phase
space structures of the photon field, particularly the appearance of a bimodal
density distribution, resembling a cat state of the photon. Interestingly,
photon mediated entanglement between the two atomic qubits, which are
otherwise non interacting, makes JCJJ a promising candidate for quantum
information processing. Using mutual information, we demonstrate how the
quantum correlation between the atomic qubits in the two cavities can be
manipulated by changing the photon population imbalance. Finally, we
investigate the quench dynamics starting from a stable steady state to an
unstable regime, which results in the formation of an incoherent gas of
photons spread over the phase space, resembling a thermal state.

The Jaynes-Cummings dimer has already been realized in a circuit QED setup
[57], as well it can be engineered by coupling the optical cavities [37, 39,
42]. The signature of the self-trapping phenomena has also been observed
experimentally in micro-cavities [87], which is promising for the observation
of different types of photonic Josephson oscillations discussed in this work.
The Kerr nonlinearity can be realized in circuit QED [62, 67] as well as in
optical cavities [63, 64, 65, 66], which is the key ingredient for the
observation of various quantum phenomena related to the steady states, such as
the revival cycle in self trapped regimes. A rich variety of collective
phenomena can also be observed in cavities containing many atoms, which have
been implemented in experiments by coupling the condensates of ultracold atoms
with a cavity mode [25, 26, 27].

Since dissipation is inherent in such atom-photon interacting systems, they
serve as an ideal platform to explore out of equilibrium dissipative phenomena
[30, 31, 32, 33, 34], which has gained momentum due to recent experiments [28,
29]. Consequently, the dynamical states of JCJJ can acquire a finite lifetime,
primarily due to weak dissipation arising from photon loss, which requires
further investigation.

In conclusion, the Josephson coupled Jaynes-Cumming dimer can serve as a test
bed to study the fascinating non-equilibrium phenomena, as well as
manipulation of entanglement between the two atomic qubits, which can have
potential applications in quantum information processing.

###### Acknowledgements.

We thank Nirmalya Ghosh, Sudip Sinha and Sayak Ray for comments and fruitful
discussions.

##  Appendix A Comparison between the cat state and dynamically generated
bimodal distribution of photons

Here, we present a comparison between the photonic cat state and the bimodal
phase space distribution of photons, which appears during the revival cycle of
the self trapped (ST1) state, shown in Fig.14(c) of Sec.V. The Husimi
distribution corresponding to this photonic state in a particular cavity
exhibits sharply peaked densities at the diametrically opposite points in
phase space, which resembles the probability distribution of a cat state [60,
14, 15, 16, 17, 18] given by,

$\displaystyle\ket{\psi_{\rm
cat}}=\frac{\ket{\alpha}+\ket{-\alpha}}{\sqrt{2(1+e^{-2|\alpha|^{2}})}},$ (31)

where the density is peaked at the two phase space points $\alpha,-\alpha$.
Although the cat state is a superposition of two coherent states, their
interference effect in the probability distribution almost disappears when
they are far apart in phase space with $|\alpha|^{2}\gg 1$. On the other hand,
the Wigner function of the cat state reveals the interference pattern at the
center, as seen in Fig.17(b).

Figure 17: Distinguishing the dynamical state of the photon field from the cat
state in terms of the Wigner function: The distribution (a) corresponds to the
quantum state representing the photon field in the left cavity in the middle
of the phase revival cycle (at a time marked by a square in Fig.14 of the main
text) and (b) the cat state given in Eq.(31). The inset of (b) represents the
Wigner function of the mixed density matrix $\hat{\rho}_{L}^{m}$ given in
Eq.(32), which does not show any coherent oscillation between the two peaks,
unlike the cat state.

In contrast, the Wigner function of the density matrix representing the
dynamical state of the photons in a particular cavity does not exhibit any
interference pattern. Nevertheless, it is localized at two diagonally opposite
points in the phase space, similar to a cat state (see Fig.17(a)).
Furthermore, the reduced density matrix of the photon field in each cavity
represents a mixed state, as well as the phase fluctuation approaches close to
its maximum limit (See Fig.14(b)). Consequently, due to the lack of coherence,
such a dynamical state of the photons can not be described by a simple linear
combination of two coherent states similar to the cat state. Instead, this
photonic state closely resembles an incoherent mixture of two coherent states
represented by the density matrix,

$\displaystyle\hat{\rho}_{i}^{m}$ $\displaystyle=$
$\displaystyle\frac{\ket{\alpha_{i}}\bra{\alpha_{i}}+\ket{-\alpha_{i}}\bra{-\alpha_{i}}}{2},$
(32)

which exhibits the bimodal phase space distribution without any interference
pattern (see inset of Fig.17(b)).

## References

  * [1] Changsuk Noh and Dimitris G. Angelakis, Rep. Prog. Phys. 80, 016401 (2016).
  * [2] J. M. Raimond, M. Brune, and S. Haroche, Rev. Mod. Phys. 73, 565 (2001).
  * [3] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * [4] K. Hepp and E. H. Lieb, Ann. Phys. (N.Y.) 76, 360 (1973).
  * [5] C. Emary and T. Brandes, Phys. Rev. Lett. 90, 044101 (2003); Phys. Rev. E 67, 066203 (2003).
  * [6] A. Baksic and C. Ciuti, Phys. Rev. Lett. 112, 173601 (2014).
  * [7] A. Altland and F. Haake, Phys. Rev. Lett. 108, 073601 (2012); New J. Phys. 14, 073011 (2012).
  * [8] 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).
  * [9] S. Pilatowsky-Cameo, D. Villaseñor, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, L. F. Santos and J. G. Hirsch, Nat Commun 12, 852 (2021).
  * [10] S. Sinha, S. Ray, and S. Sinha, Journal of Physics: Condensed Matter 33, 174005 (2021).
  * [11] J. Braumüller, M. Marthaler, A. Schneider, A. Stehli, H. Rotzinger, M. Weides, and A. V. Ustinov, Nat Commun 8, 779 (2017).
  * [12] J. I. Cirac, R. Blatt, A. S. Parkins, and P. Zoller, Phys. Rev. A 49, 1202 (1994).
  * [13] A. De and R. Joynt, Phys. Rev. A 87, 042336 (2013).
  * [14] L. Sun, A. Petrenko, Z. Leghtas, B. Vlastakis, G. Kirchmair, K. M. Sliwa, A. Narla, M. Hatridge, S. Shankar, J. Blumoff, L. Frunzio, M. Mirrahimi, M. H. Devoret and R. J. Schoelkopf, Nature 511, 444–448 (2014).
  * [15] X. Pan, J. Schwinger, Ni-Ni Huang, P. Song, W. Chua, F. Hanamura, A. Joshi, F. Valadares, R. Filip, and Y. Y. Gao, Phys. Rev. X 13, 021004 (2023).
  * [16] J. Lebreuilly, C. Aron, and C. Mora, Phys. Rev. Lett. 122, 120402 (2019).
  * [17] S. M. Girvin, arXiv:1710.03179.
  * [18] V. Bužek, H. Moya-Cessa, P. L. Knight, and S. J. D. Phoenix, Phys. Rev. A 45, 8190 (1992).
  * [19] C. Navarrete-Benlloch, J. J. García-Ripoll, and D. Porras, Phys. Rev. Lett. 113, 193601 (2014).
  * [20] M. J. Gullans, J. Stehlik, Y.-Y. Liu, C. Eichler, J. R. Petta, and J. M. Taylor, Phys. Rev. Lett. 117, 056801 (2016).
  * [21] G. Gasse, C. Lupien, and B. Reulet, Phys. Rev. Lett. 111, 136601 (2013).
  * [22] S. J. Masson and S. Parkins, Phys. Rev. A 99, 023822 (2019).
  * [23] A. Majumdar, M. Bajcsy, and J. Vučković, Phys. Rev. A 85, 041801(R) (2012).
  * [24] K. Weiher, E. Agudelo, and M. Bohmann, Phys. Rev. A 100, 043812 (2019).
  * [25] J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Nature (London) 543, 87 (2017).
  * [26] J. Léonard, A. Morales, P. Zupancic, T. Donner, and T. Esslinger, Science 358, 1415 (2017).
  * [27] J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, Proc. Natl. Acad. Sci. U.S.A. 112, 3290 (2015).
  * [28] F. Brennecke, R. Mottl, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Proc. Natl. Acad. Sci. USA 110, 11763 (2013).
  * [29] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch and A. A. Houck, Phys. Rev. X 7, 011016 (2017).
  * [30] H. J. Carmichael, Phys. Rev. X 5, 031028 (2015).
  * [31] F. Reiter, T. L. Nguyen, J. P. Home, and S. F. Yelin, Phys. Rev. Lett. 125, 233602 (2020).
  * [32] K. C. Stitely, A. Giraldo, B. Krauskopf, and S. Parkins, Phys. Rev. Research 2, 033131 (2020).
  * [33] J. Li, R. Fazio, and S. Chesi, New J. Phys. 24, 083039 (2022).
  * [34] W. Kopylov, M. Radonjić, T. Brandes, A. Balaž, and A. Pelster, Phys. Rev. A 92, 063832 (2015).
  * [35] E.T. Jaynes, F.W. Cummings, Proc. IEEE. 51 (1): 89–109 (1963).
  * [36] M. Tavis and F. W. Cummings, Phys. Rev. 170, 379 (1968).
  * [37] M. Hartmann, F. G. S. L. Brandão and M. B. Plenio, Nature Phys 2, 849–855 (2006).
  * [38] A. Greentree, C. Tahan, J. Cole et al., Nature Phys 2, 856–861 (2006).
  * [39] G. Lepert, M. Trupke, M. J. Hartman, M. B. Plenio, and E. A. Hinds, New J. Phys. 13, 113002 (2011).
  * [40] Michael J. Hartmann, J. Opt. 18, 104005 (2016).
  * [41] I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * [42] D. L. Underwood, W. E. Shanks, Jens Koch and A. A. Houck, Phys. Rev. A 86, 023837 (2012).
  * [43] P. Leboeuf and S. Moulieras, Phys. Rev. Lett. 105, 163904 (2010).
  * [44] L. Dominici, R. Carretero-González, A. Gianfrate et al., Nat Commun 9, 1467 (2018).
  * [45] K. Lagoudakis, M. Wouters, M. Richard et al., Nature Phys 4, 706–710 (2008).
  * [46] M. J. Hartmann, F. G. S. L. Brandão, and M. B. Plenio, Phys. Rev. Lett. 99, 160501 (2007).
  * [47] T. Byrnes, N. Kim and Y. Yamamoto, Nature Phys 10, 803–813 (2014).
  * [48] D. Rossini and R. Fazio, Phys. Rev. Lett. 99, 186401 (2007).
  * [49] J. Cho, D. G. Angelakis, and S. Bose, Phys. Rev. Lett. 101, 246809 (2008).
  * [50] J. Koch, A. A. Houck, K. Le Hur, and S. M. Girvin, Phys. Rev. A 82, 043811 (2010).
  * [51] S. Schmidt and G. Blatter, Phys. Rev. Lett. 103, 086403 (2009).
  * [52] J. Koch and K. Le Hur, Phys. Rev. A 80, 023811 (2009).
  * [53] M. Knap, E. Arrigoni, and W. von der Linden, Phys. Rev. B 82, 045126 (2010).
  * [54] L. Guo, S. Greschner, S. Zhu, and W. Zhang, Phys. Rev. A 100, 033614 (2019).
  * [55] D. G. Angelakis, M. F. Santos, and S. Bose, Phys. Rev. A 76, 031805(R) (2007).
  * [56] N. Na, S. Utsunomiya, L. Tian, and Y. Yamamoto, Phys. Rev. A 77, 031803 (2008).
  * [57] J. Raftery, D. Sadri, S. Schmidt, H. E. Türeci, and A. A. Houck, Phys. Rev. X 4, 031043 (2014).
  * [58] S. Schmidt, D. Gerace, A. A. Houck, G. Blatter, and H. E. Türeci, Phys. Rev. B 82, 100507(R) (2010).
  * [59] A. Dey and M. Kulkarni, Phys. Rev. A 101, 043801 (2020).
  * [60] D. F. Walls and G. J. Milburn (1995). Quantum optics. Berlin; New York: Springer-Verlag.
  * [61] S. D. Du and C. D. Gong, Phys. Rev. A 50, 779 (1994).
  * [62] S. Rebić, J.Twamley, and G. J. Milburn, Phys. Rev. Lett. 103, 150503 (2009).
  * [63] A. Imamoğlu, H. Schmidt, G. Woods, and M. Deutsch, Phys. Rev. Lett. 79, 1467 (1997).
  * [64] H. Schmidt and A. Imamoğlu, Opt. Lett. 21, 1936 (1996).
  * [65] H. Rokhsari and K. J. Vahala, Opt. Lett. 30, 427 (2005).
  * [66] S. Rebic, S. M. Tan, A. S. Parkins, and D. F. Walls, J. Opt. B 1, 490 (1999).
  * [67] M. Kounalakis, C. Dickel, A. Bruno, N. Langford, and G. Steele, npj Quantum Inf. 4, 38 (2018).
  * [68] Z. R. Gong, H. Ian, Yu-xi Liu, C. P. Sun, and F. Nori, Phys. Rev. A 80, 065801 (2009).
  * [69] H. Shapourian and D. Sadri, Phys. Rev. A 93, 013845 (2016).
  * [70] P.A.M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930); J. Frenkel, Wave Mechanics, Claredon Press, Oxford, 1934.
  * [71] J. M. Radcliffe, J. Phys. A: Gen.Phys., 4, 313 (1971).
  * [72] J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J. Young, J. R. K. Cline, A. M. Rey, and J. K. Thompson, Nature 580, 602 (2020).
  * [73] S. A. Weidinger, M. Heyl, A. Silva, and M. Knap, Phys. Rev. B 96, 134313 (2017).
  * [74] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 (1997).
  * [75] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, Phys. Rev. A 59, 620 (1999).
  * [76] E. Boukobza, M. Chuchem, D. Cohen, and A. Vardi, Phys. Rev. Lett. 102, 180403 (2009).
  * [77] E. Boukobza, M. G. Moore, D. Cohen, and A. Vardi, Phys. Rev. Lett. 104, 240402 (2010).
  * [78] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
  * [79] R. Gati and M. K. Oberthaler, J. Phys. B 40, R61 (2007).
  * [80] M. Tabor, “Chaos and Integrability in nonlinear Dynamics: An Introduction,” Wiley-Interscience, USA, 1989.
  * [81] R. Mottl, F. Brennecke, K. Baumann, R. Landig, T. Donner, and T. Esslinger, Science 336, 1570 (2012).
  * [82] S. Pramanik, S. Bandyopadhyay, and M. Cahay, Phys. Rev. B 68, 075313 (2003).
  * [83] V. Sevriuk, K. Y. Tan, E. Hyyppä, M. Silveri, M. Partanen, M. Jenei, S. Masuda, J. Goetz, V. Vesterinen, L. Grönberg, and M. Möttönen, Appl. Phys. Lett. 115, 082601 (2019).
  * [84] A. J. Fleisher, D. A. Long, Q. Liu, and J. T. Hodges, Phys. Rev. A 93, 013833 (2016).
  * [85] T. Lauprêtre, C. Proux, R. Ghosh, S. Schwartz, F. Goldfarb, and F. Bretenaker, Opt. Lett. 36, 1551 (2011).
  * [86] Y. Okawachi, B. Y. Kim, Y. Zhao, X. Ji, M. Lipson, and A. L. Gaeta, Optica 8, 1458–1461 (2021).
  * [87] M. Abbarchi, A. Amo, V. G. Sala, D. D. Solnyshkov, H.Flayac, L. Ferrier, I. Sagnes, E. Galopin, A. Lemaître, G. Malpuech, and J. Bloch, Nat. Phys. 9, 275 (2013).
  * [88] D. T. Pegg and S. M. Barnett, Phys. Rev. A 39, 1665 (1989).
  * [89] M. Nielsen and I. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, England, 2000).
  * [90] D. P. DiVincenzo, M. Horodecki, D. W. Leung, J. A. Smolin, and B. M. Terhal, Phys. Rev. Lett. 92, 067902 (2004).
  * [91] G. Adesso and A. Datta, Phys. Rev. Lett. 105, 030501 (2010).
  * [92] X.-M. Lu, J. Ma, Z. Xi, and X. Wang, Phys. Rev. A 83, 012327 (2011).
  * [93] L. Henderson and V. Vedral, J. Phys. A: Math. Gen. 34, 6899 (2001).
  * [94] P. Das, D. S. Bhakuni, and A. Sharma, Phys. Rev. A 107, 043706 (2023).
  * [95] M. Greiner, O. Mandel, T. W. Hänsch, I. Bloch, Nature 419, 51–54 (2002).
  * [96] C. Gerry and P. Knight, Introductory Quantum Optics (Cambridge University Press, Cambridge, 2004).
  * [97] S. Sinha and S. Sinha, Phys. Rev. Lett. 125, 134101 (2020).
  * [98] D. Mondal, S. Sinha, S. Ray, J. Kroha, and S. Sinha, Phys. Rev. A 106, 043321 (2022).
  * [99] J. Klaers, F. Vewinger and M. Weitz, Nature Phys 6, 512–515 (2010).