# A Floquet analysis perspective of driven light-matter interaction models

Jonas Larson<EMAIL_ADDRESS>Department of Physics, Stockholm
University, SE-106 91 Stockholm, Sweden

###### Abstract

In this paper, we analyze the harmonically driven Jaynes-Cummings and Lipkin-
Meshkov-Glick models using both numerical integration of time-dependent
Hamiltonians and Floquet theory. For a separation of time-scales between the
drive and intrinsic Rabi oscillations in the former model, the driving results
in an effective periodic reversal of time. The corresponding Floquet
Hamiltonian is a Wannier-Stark model, which can be analytically solved.
Despite the chaotic nature of the driven Lipkin-Meshkov-Glick model, moderate
system sizes can display qualitatively different behaviors under varying
system parameters. Ergodicity arises in systems that are neither adiabatic nor
diabatic, owing to repeated multi-level Landau-Zener transitions. Chaotic
behavior, observed in slow driving, manifests as random jumps in the
magnetization, suggesting potential utility as a random number generator.
Furthermore, we discuss both models in terms of what we call Floquet Fock
state lattices.

##  I Introduction

The Jaynes-Cummings (JC) model stands as a cornerstone in our explorations of
light-matter interaction within the realm of quantum optics Jaynes and
Cummings (1963); Shore and Knight (1993). Despite its apparent simplicity, the
model’s profound impact was evidenced nearly four decades ago when
experimental validation showcased its efficacy, see for example the review
Larson and Mavrogordatos (2022). This model, depicting a two-level system
(representing matter, often an atom) coupled to a harmonic oscillator
(typically a single mode of the electromagnetic field), has since served as a
catalyst for deeper explorations into quantum dynamics. One might expect that
only under extreme conditions would the JC model accurately encapsulate the
complexities of tailored light-matter interaction. However, its validation in
experimental settings and its generalizations to more intricate models have
solidified its significance.

Central to the JC model and similar models in quantum optics is the
limitations few degrees of freedom alongside certain symmetries. The adoption
of the rotating wave approximation (RWA) further enhances the model’s utility,
ensuring the preservation of the number of excitations within the system. This
makes such models ideal for probing fundamental questions within quantum
mechanics Haroche and Raimond (2006). Furthermore, the introduction of
classical driving forces can result in novel phenomena. Again, either a
rotating frame is considered, or an RWA is applied, leading to an effective
time-independent Hamiltonian. For certain drives, however, an RWA might not be
justified, and one must consider the full time-dependent problem.

Early investigations into periodically driven JC models, such as that by
Schlicher Schlicher (1989), gave insight into more complex scenarios. In this
vein, we extend our analysis to encompass sinusoidal coupling within the JC
model and explore its implications. Additionally, we delve into the
sinusoidally driven Lipkin-Meshkov-Glick (LMG) model Lipkin et al. (1965);
Dusuel and Vidal (2004); Ribeiro et al. (2007); Morrison and Parkins (2008);
Larson (2010); Latorre et al. (2005), broadening our scope to encompass
diverse manifestations of light-matter interaction.

To facilitate our study, we employ Floquet theory Holthaus (2015); Eckardt
(2017); Rudner and Lindner (2020), a tool for tackling time-dependent quantum
problems. Floquet theory enables the transformation of such problems into
equivalent time-independent formulations, albeit within an extended state
space termed the Floquet space. While numerically challenging due to the
infinite nature of this space, Floquet analysis offers invaluable insights
into the system’s dynamics.

In this work, we introduce the Floquet Fock state lattice, an extension of the
Fock state lattice (FSL) to the Floquet space. The Floquet FSL serves as a
representation of the Floquet Hamiltonian, shedding light on system evolution
and revealing phenomena such as Bloch oscillations and level repulsion,
indicative of chaos.

The paper is organized as follows: In the next section, we provide an overview
of Floquet theory, particularly as applied to harmonic driving, where the
Floquet Hamiltonian can be easily determined. Insec. III, we examine the
driven JC model, focusing on relatively slow drives where the driving acts as
an effective time reversal. The subsequent section IV addresses the LMG model,
where we conduct an analysis of the time-dependent problem, discussing various
manifestations of chaos. Additionally, we describe the time evolution in the
Floquet FSL. Finally, we conclude in sec. V with a summary. Appendix A
presents the basics of FSLs and how the Floquet FSL emerges.

##  II Floquet theory for a harmonical drive

The analogue of the Bloch theorem for time-periodic Hamiltonians is the
Floquet theorem Shirley (1965). However, in quantum mechanics, time and space
play distinct roles, as time is not an ordinary observable that can be
assigned an operator. Nonetheless, as we will illustrate, it is possible to
establish a correspondence of plane waves for the time ’coordinate’ and
incorporate the time-dependence of the Hamiltonian into a new (Floquet)
Hamiltonian, residing in an expanded space comprising both the original
Hilbert space and an additional Floquet space denoted as
$\mathcal{H}\otimes\mathcal{T}$. This approach reveals that time evolution can
be separated into slow and fast-scale components. The Floquet Hamiltonian
serves as the time-evolution operator for the slow dynamics occurring over a
period $T$, while the short time-scale evolution within one period is
encapsulated by a stroboscopic kick operator. Although this topic is not
commonly addressed in standard textbooks, several comprehensive introductions
to the field of Floquet theory are available, as exemplified by references
Holthaus (2015); Eckardt (2017); Rudner and Lindner (2020). We present
fundamental concepts here and direct interested readers to the aforementioned
references for a more in-depth exploration.

The starting point is a periodic Hamiltonian $\hat{H}(t+T)=\hat{H}(t)$, with
$T$ the period. The theorem states that there exists an orthonormal basis
$\\{|\psi_{n}(t)\rangle\\}$, which under time-evolution of the Hamiltonian
$\hat{H}(t)$, obeying

$|\psi_{n}(t+T)\rangle=e^{-i\varepsilon_{n}T}|\psi_{n}(t)\rangle.$ (1)

The states $|\psi_{n}(t)\rangle$ are called the Floquet states,
$\varepsilon_{n}$ are the quasi-energies, and we use $\hbar=1$ throughout.
Given the Floquet states, the strength of the theorem is that a general
initial state $|\psi(0)\rangle$ will follow the simple evolution

$|\psi(t)\rangle=\sum_{n}a_{n}|\psi_{n}(t)\rangle,$ (2)

where, importantly, the coefficients $a_{n}$ are time independent. To find the
Floquet states we note that they can be decomposed into a time-periodic state
multiplied be a phase factor as111We could shift the phase factor,
$\exp[-i\varepsilon_{n}t]\rightarrow\exp[-i\varepsilon_{n}(t-\tau_{0})]$ for
some $\tau_{0}$, and it would modify the definition of the Floquet
Hamiltonian. This is a gauge freedom we have.

$|\psi_{n}(t)\rangle=e^{-i\varepsilon_{n}t}|\phi_{n}(t)\rangle,\hskip
34.14322pt|\phi_{n}(t+T)\rangle=|\phi_{n}(t)\rangle.$ (3)

By substituting (3) into the Schrödinger equation one derives

$\begin{array}[]{c}\left[\varepsilon_{n}+i\frac{d}{dt}\right]|\phi_{n}(t)\rangle=\hat{H}(t)|\phi_{n}(t)\rangle\\\
\Leftrightarrow\\\
\hat{H}^{\mathrm{F}}|\phi_{n}(t)\rangle\equiv\left[\hat{H}(t)-i\frac{d}{dt}\right]|\phi_{n}(t)\rangle=\varepsilon_{n}|\phi_{n}(t)\rangle.\end{array}$
(4)

Here we have introduced the (time-independent) Floquet Hamiltonian
$\hat{H}_{\mathrm{F}}$. Furthermore, owing to the periodicity of the
Hamiltonian, we can Fourier expand the periodic states as

$|\phi_{n}(t)\rangle=\sum_{m=-\infty}^{\infty}|\varphi_{n}(t),m\rangle\equiv\sum_{m=-\infty}^{\infty}e^{im\omega
t}|\varphi_{n},m\rangle,$ (5)

where $\omega=2\pi/T$, and $|\varphi_{n},m\rangle$ is the $m$’th Fourier
coefficient of $|\phi_{n}(t)\rangle$. We note that
$|\varphi_{n}(t),m\rangle=e^{im\omega t}|\varphi_{n}(t),0\rangle$, and
moreover

$\cos(\omega
t)|\varphi_{n}(t),m\rangle=\frac{1}{2}\left[|\varphi_{n}(t),m+1\rangle+|\varphi_{n}(t),m-1\rangle\right].$
(6)

Thus, $\exp(\pm i\omega t)$ acts as an effective raising/lowering operator for
the states $|\varphi_{n},m\rangle$. It is clear that the Floquet state (3) is
invariant under the shifts
$\varepsilon_{n}\rightarrow\varepsilon_{n}+m^{\prime}\omega$ and
$|\varphi_{n}(t),m\rangle\rightarrow|\varphi_{n}(t),m-m^{\prime}\rangle$. As a
result, the quasi-energies $\varepsilon_{n}$ can be restricted to any interval
of length $\omega$, say $-\omega/2\leq\varepsilon_{n}<\omega/2$, which serves
as the analogue of a Brillouin zone in frequency rather than in momentum
space.

To solve the full problem we still need to find $\varepsilon_{n}$ and
$|\phi_{n}(t)\rangle$. Making use of the Fourier expansion (5) in eq. (4)
gives

$\left(\sum_{m^{\prime}=-\infty}^{\infty}H^{(m-m^{\prime})}+m\omega\right)|\varphi_{n},m^{\prime}\rangle=\varepsilon_{n}|\varphi_{n},m\rangle,$
(7)

with $H^{(m)}$ the Fourier coefficient of the Hamiltonian, i.e.
$\hat{H}(t)=\sum_{m=-\infty}^{\infty}e^{im\omega t}H^{(m)}$. Equation (7) can
be put on a matrix form
$\hat{H}_{\mathrm{F}}|\varphi_{n}\rangle=\varepsilon_{n}|\varphi_{n}\rangle$,
with the Floquet Hamiltonian expressed in the $m$-basis as

$\hat{H}^{\mathrm{F}}=\left[\begin{array}[]{cccc}\ddots&H^{(-1)}&H^{(-2)}&\\\
H^{(1)}&H^{(0)}-m\omega&H^{(-1)}&H^{(-2)}\\\
H^{(2)}&H^{(1)}&H^{(0)}-(m+1)\omega&H^{(-1)}\\\
&H^{(2)}&H^{(1)}&\ddots\end{array}\right].$ (8)

On the diagonal, $H^{(0)}$ is the time-averaged Hamiltonian. Given that the
dimension of the Hilbert space is $d$, each time-independent block-Hamiltonian
$H^{(m)}$ is also of dimension $d$. The dimension of the Floquet Hamiltonian
is infinite since $m$ runs over all integers, and for numerics an
$m$-truncation must be considered. What we have achieved is to transform the
time-dependent problem into a time-independent eigenvalue problem, which
solved provides the quasi-energies $\varepsilon_{n}$ and the Fourier
coefficients for $|\phi_{n}(t)\rangle$. With these, the Floque states are
obtained from eq. (3).

The aforementioned separation of time-scales is not transparent from the
discussion so far. First we note that for one period, the Floquet state (3) is
multiplied by a phase factor determined by the quasi-energy. Thus, in this
case we can express the time-evolution operator over one period, and starting
at $t_{0}$, in terms of the Floquet Hamiltonian as

$\hat{U}(t_{0}+T,t_{0})=e^{-i\hat{H}^{\mathrm{F}}T}.$ (9)

In principle, the Floquet Hamiltonian can depend on the (gauge) choice
$t_{0}$. It can be shown222This actually defines the operator
$\hat{K}^{\mathrm{F}}$, given the Floquet Hamiltonian $\hat{H}^{\mathrm{F}}$
and the gauge choice. Holthaus (2015), that the full time-evolution operator
could be split into two parts

$\hat{U}(t,t_{0})=e^{-i\hat{K}^{\mathrm{F}}t}e^{-i\hat{H}^{\mathrm{F}}t},$
(10)

where $\hat{K}^{\mathrm{F}}$ is sometimes referred to as stroboscopic kick
operator. It obeys
$e^{-i\hat{K}^{\mathrm{F}}(t+T)}=e^{-i\hat{K}^{\mathrm{F}}t}$, and given this
partitioning of the time-evolution operator we identify $\hat{H}^{\mathrm{F}}$
as generating the slow time evolution over a full period, while
$\hat{K}^{\mathrm{F}}$ provides the short time evolution occurring within one
period. The above result is the Floquet theorem expressed in terms of the
time-evolution operator.

Hence, we have re-written the time-dependent problem, characterized by a time-
evolution operator
$\hat{U}(t,0)=T\exp\left[-i\int_{0}^{t}\hat{H}(t^{\prime})dt^{\prime}\right]$,
with $T$ the time-ordering operator, into an effectively time-independent
problem (10). The new problem is to be solved in the extended Floquet space
that is comprised by the original Hilbert space $\mathcal{H}$ and the space
$\mathcal{T}$ of the time-periodic states $|\phi_{n}(t)\rangle$, i.e.
$\mathcal{F}=\mathcal{H}\otimes\mathcal{T}$. While time-independent problems
are typically far more easy to solve than time-dependent ones, we pay the
price that we are now dealing with a larger state space. Numerically, even if
the dimension of $\mathcal{H}$ is finite, we face an infinite state space
since the Fourier expansion (5) runs from $-\infty$ to $+\infty$, and a
truncation must be introduced Rudner and Lindner (2020). Actually, it is not
uncommon to find that numerically solving the original time-dependent problem
is more time-efficient than solving the time-independent Floquet problem.
Nevertheless, this might give some valuable insight into the dynamics.

In the examples to follow we consider harmonic driving characterized by some
frequency $\omega$, i.e. $\hat{H}(t)=\hat{H}\left(\cos(\omega t)\right)$.
There is a recipe how to construct the Floquet Hamiltonian in this case Klimov
and Chumakov (2009). To this end we note that the cosine function couples the
$m$-state to the $(m\pm 1)$-states in eq. (6). Thus, in the extended space
spanned by the $\\{|\varphi,m\rangle\\}$ basis, $\cos(\omega t)$ is a tri-
diagonal matrix with $1/2$ on its first off-diagonals. For brevity, if we
write $|n,m\rangle\equiv|\varphi_{n},m\rangle$ and introduce the ladder
operators $\hat{E}^{\pm}|n,m\rangle=|n,m\pm 1\rangle$, eq. (6) can be
compactly written as

$\left(\hat{E}^{+}+\hat{E}^{-}\right)|n,m\rangle=|n,m+1\rangle+|n,m-1\rangle.$
(11)

The ladder operators, complemented with the diagonal operator
$\hat{E}_{0}|n,m\rangle=m|n,m\rangle$, form what is called a Euclidean algebra
that obey the commutation relations
$\left[\hat{E}_{0},\hat{E}^{\pm}\right]=\pm\hat{E}^{\pm}$ and
$\left[\hat{E}^{-},\hat{E}^{+}\right]=0$ Klimov and Chumakov (2009).
Furthermore, we may note

$\hat{U}_{0}(t)\hat{E}^{\pm}\hat{U}_{0}^{\dagger}(t)=\hat{E}^{\pm}e^{\mp
i\omega t},$ (12)

where $\hat{U}_{0}(t)=\exp\left(-i\hat{E}_{0}\omega t\right)$. Let us
introduce the phase states
$\hat{E}^{-}|\theta\rangle=e^{i\theta}|\theta\rangle$, and the time-dependent
state $|\theta(t)\rangle=\hat{U}_{0}(t)|\theta\rangle$. We can then construct
a Floquet Hamiltonian $\hat{H}^{\mathrm{F}}$ (in the extended space
$\mathcal{F}$) such that the time-dependent Hamiltonian (in the Hilbert space
$\mathcal{H}$) is
$\hat{H}(t)=\langle\theta(t)|\hat{H}^{\mathrm{F}}|\theta(t)\rangle$. Note that
the phase states $|\theta(t)\rangle\propto\sum_{m}e^{-i\theta m}|m\rangle$ are
also eigenstates of $\hat{E}^{+}+\hat{E}^{-}$ with eigenvalues
$2\cos(\theta)$, and we have a freedom to choose the phase, for example
$\theta=0$.

Give the above, the recipe how to construct this Floquet Hamiltonian follows
from the substitution

$\cos(\omega t)\rightarrow\frac{\hat{E}^{+}+\hat{E}^{-}}{2},\hskip
39.83368pt-i\frac{d}{dt}\rightarrow\omega\hat{E}_{0},$ (13)

and hence

$\hat{H}^{\mathrm{F}}=\hat{H}\left(\left(\hat{E}^{+}+\hat{E}^{-}\right)/2\right)+\omega\hat{E}_{0}.$
(14)

Before proceeding, let us define number, or Fock, states $|n\rangle$ in the
general way, namely as the eigenstates of the number operator. Apart from the
standard ones bosonic and fermionic ones, we complement the these with spin
and Floquet states;

$\begin{array}[]{lll}\mathrm{Bosonic}:&&\hat{n}=\hat{a}^{\dagger}\hat{a},\,\,\,\,\,\hat{n}|n\rangle=n|n\rangle,\,\,\,\,n\in\mathbb{Z}^{+},\\\
\mathrm{Fermionic}:&&\hat{n}=\hat{c}^{\dagger}\hat{c},\,\,\,\,\,\hat{n}|n\rangle=n|n\rangle,\,\,\,\,n=0,\,1,\\\
\mathrm{Spin}:&&\hat{S}_{z}=\frac{\hat{S}^{+}\hat{S}^{-}-\hat{S}^{-}\hat{S}^{+}}{2},\,\,\,\,\,\hat{S}_{z}|m_{s}\rangle=m_{s}|m_{s}\rangle,\\\
\mathrm{Qubit}:&&\hat{\sigma}_{z}=2\hat{\sigma}^{+}\hat{\sigma}^{-}-\mathbb{I},\,\,\,\,\,\hat{\sigma}_{z}|s\rangle=s|s\rangle,\,\,\,\,s=\pm
1,\\\
\mathrm{Floquet}:&&\hat{E}_{0},\,\,\,\,\,\hat{E}_{0}|m\rangle=m|m\rangle,\,\,\,\,m\in\mathbb{Z}.\\\
\end{array}$ (15)

Note that given the properties of the Euclidean algebra, its number operator
$\hat{E}_{0}$ cannot simply be expressed in terms of its ladder operators
$\hat{E}^{\pm}$.

##  III Driven Jaynes-Cummings model

The Jaynes-Cummings model, comprising a two-level system, referred to as an
atom, coherently coupled to a single boson mode, playing the role of the
quantized electromagnetic field, is given by Jaynes and Cummings (1963); Shore
and Knight (1993); Larson and Mavrogordatos (2022)

$\hat{H}_{\mathrm{JC}}=\frac{\Delta}{2}\hat{\sigma}_{z}+g(t)\left(\hat{a}^{\dagger}\hat{\sigma}^{-}+\hat{a}\hat{\sigma}^{+}\right).$
(16)

Here, $\Delta$ is the atom-field detuning, $g(t)$ is the light-matter
coupling, $\hat{a}^{\dagger}/\hat{a}$ are the photon creation/annihilation
operators obeying $\left[\hat{a},\hat{a}^{\dagger}\right]=1$, and the
$\hat{\sigma}$-operators are the Pauli $z$\- and raising/lowering matrices,
i.e. $\left[\hat{\sigma}^{+},\hat{\sigma}^{-}\right]=\hat{\sigma}_{z}$ and
$\left[\hat{\sigma}^{\pm},\hat{\sigma}_{z}\right]=\mp 2\hat{\sigma}^{\pm}$. We
denote the eigenstates of the $z$ Pauli matrix as $|g\rangle$ (ground) and
$|e\rangle$ (excited).

The Hamiltonian can be derived from a minimal-coupling Hamiltonian under a
series of approximations, see Larson and Mavrogordatos (2022). Most important
for us is the RWA that accounts for neglecting the counter-rotating terms
$\hat{a}^{\dagger}\hat{\sigma}^{+}$ and $\hat{a}\hat{\sigma}^{-}$. In doing
so, the number of excitations $\hat{N}=\hat{n}+\frac{1}{2}\hat{\sigma}_{z}$,
with $\hat{n}=\hat{a}^{\dagger}\hat{a}$, is preserved such that the model
possess a continuous $U(1)$ symmetry and the Hamiltonian is $2\times 2$ block-
diagonal in the Fock basis $\left\\{|n,g\rangle,|n,e\rangle\right\\}$333Each
block $\hat{h}_{n}$ constitutes a $2\times 2$ matrix apart from the vacuum
$|0,g\rangle$ that is decoupled from any other states in the RWA.. Thus, for
each excitation $N$ sector we find the block Hamiltonian, spanned by the
states $|n,e\rangle$ and $|n+1,g\rangle$, as

$\hat{h}_{n}(t)=\left[\begin{array}[]{cc}\Delta/2&g_{0}\cos(\omega
t)\sqrt{n}\\\ g_{0}\cos(\omega t)\sqrt{n}&-\Delta/2\end{array}\right],$ (17)

where we assumed $g(t)=g_{0}\cos(\omega t)$. Let us for simplicity consider
the case $\Delta=0$ where it is straightforward to decouple the two equations
via a Hadamard transformation with the action
$\hat{\sigma}_{z}\leftrightarrow\hat{\sigma}_{x}$. With
$|n\rangle_{x\pm}=\left(|n-1,e\rangle\pm|n,g\rangle\right)/\sqrt{2}$ (and
$|0\rangle_{x}=|0,g\rangle$) the eigenstates of the above Hamiltonian we can
write down the time-evolved state for a general initial state as

$\begin{array}[]{lll}|\psi(t)\rangle&=&\displaystyle{d_{0}|n\rangle_{x}+\sum_{n=1}^{\infty}\sum_{s=\pm}d_{n\pm}e^{\pm
i\Omega_{n}(t)}|n\rangle_{x\pm}}\\\
&=&\displaystyle{d_{0}|n\rangle_{x}+\sum_{n=1}^{\infty}\sum_{s=g,e}c_{ns}(t)|n,s\rangle,}\end{array}$
(18)

with the time-dependent “Rabi frequency”
$\Omega_{n}(t)=\sqrt{n}g_{0}\sin(\omega t)/\omega$, and

$\begin{array}[]{l}c_{ne}(t)=\frac{1}{\sqrt{2}}\left(d_{n+}e^{i\sqrt{n}\sin(\omega
t)/\omega}+d_{n-}e^{-i\sqrt{n}\sin(\omega t)/\omega}\right),\\\
c_{ng}(t)=\frac{1}{\sqrt{2}}\left(d_{n+}e^{i\sqrt{n}\sin(\omega
t)/\omega}-d_{n-}e^{-i\sqrt{n}\sin(\omega t)/\omega}\right).\end{array}$ (19)

If, for example, the atom is initially in its excited state, $c_{ng}(0)=0$, we
have $d_{n+}=d_{n-}\equiv d_{n}$. If the photon field is initially in a
coherent state with an amplitude $\alpha$, we further have
$\sqrt{2}d_{n}=e^{-|\alpha|^{2}}\frac{\alpha^{n}}{\sqrt{n!}}$. The atomic
inversion reads Prants and Yacoupova (1992); Joshi and Lawande (1993)

$W(t)\equiv\langle\hat{\sigma}_{z}\rangle=\sum_{n}2|d_{n}|^{2}\cos\left[2\Omega_{n+1}(t)\right].$
(20)

As for the regular JC model Meystre et al. (1973); Narozhny et al. (1981), for
very slow driving ($\omega\ll 1$) the above inversion will display the
familiar pattern of collapses and revivals as shown in fig. 1. What is
different from the regular JC model is the sinusoidal dependence of the Rabi
frequency $\Omega_{n}(t)$. Initially the argument of the cosine function
increases from zero, exactly like the time-independent JC model (where the
argument of the cosine $g\sqrt{n}t$), but contrary to that model where the
argument grows indefinitely, in our model it starts to decrease after
$t=t_{\mathrm{tr}}=\pi/2\omega$. After a time $2t_{\mathrm{tr}}$ the
sinusoidal is again zero and we must have that every term in the sum (20)
return back in phase leading to a perfect revival. The sinusoidal time
dependence effectively reverse the evolution and clearly, this behaviour
repeats as time progresses. The first revival time of the JC model is
$t_{\mathrm{r}}=2\pi|\alpha|/g$ Larson and Mavrogordatos (2022), and hence, if
$t_{\mathrm{tr}}<t_{\mathrm{r}}$ there is no time for the JC evolution to
establish a revival, and the only revival will be the one stemming from the
effective time reversal at $2t_{\mathrm{tr}}$. On the other hand, if
$t_{\mathrm{tr}}>t_{\mathrm{r}}$ the system will display at least one revival
before time is effective reversed. We should point out that the effective
time-reversal evolution appears from considering the resonant case $\Delta=0$,
while for $\Delta\neq 0$ this effect is lost.

There are two additional time scales in the problem; the inverse Rabi
frequency $t_{\mathrm{R}}=\pi/g|\alpha|$ and the collapse time
$t_{\mathrm{c}}=1/2g$ Larson and Mavrogordatos (2022), and in principle
$t_{\mathrm{tr}}$ should also be compared to these. For example, for very
rapid driving, $\omega\gg g$ the inversion will be hindered. To lowest order
such dynamical freezing can be understood from time-averaging the Hamiltonian
which would result in a decoupling of the atom from the field. However, more
insight is given by performing a second rotating wave approximation of the
Hamiltonian (17) to find

$\hat{h}_{n}\approx\left[\begin{array}[]{cc}\frac{\Delta-\omega}{2}&\frac{g_{0}}{2}\sqrt{n}\\\
\frac{g_{0}}{2}\sqrt{n}&-\frac{\Delta-\omega}{2}\end{array}\right].$ (21)

Figure 1: The atomic inversion as given in eq. (20) shown as the black curve,
and the corresponding numerical results calculated from the truncated Floquet
Hamiltonian (22) for a lattice of 1401 sites (light blue). The other
parameters are $\omega=0.025$ and $g=1$, and the initial photon mode is in a
coherent state with $\alpha=5$, and the atom in $|g\rangle$. At
$t_{\mathrm{tr}}=\pi/2\omega$ the evolution is effectively reversed to result
in a perfect revival at $t=\pi/\omega$. The discrepancy between the two
results derives from the Floquet truncation; the larger lattices, the more
accurate results at later times.

Let us now take the alternative approach and analyze the problem in terms of
Floquet theory by utilizing the mapping (13). For simplicity we consider the
Fock states $|n\rangle_{x\pm}$, such that the Hamiltonian is diagonal within
the two-level structure, and the Floquet Hamiltonian becomes

$\hat{H}_{\mathrm{JC}}^{\mathrm{F}}(n)=\omega\hat{E}_{0}+\frac{g\sqrt{n}}{2}\left(\hat{E}^{+}+\hat{E}^{-}\right)\hat{\tau}_{z},$
(22)

where $\hat{\tau}_{z}|n\rangle_{x\pm}=\pm|n\rangle_{x\pm}$. We recall that $n$
in the above equation represents the photon number. In the appendix A we
discuss the idea of FSLs Saugmann and Larson (2023). As done in the appendix,
we can generalize FSLs to a Floquet FSL; the composition of the Floquet space
$\mathcal{F}=\mathcal{H\otimes T}$ results in an extended FSL where harmonic
drive provides an additional dimension to the lattice. In the present model,
we see that each copy of the Hamiltonian (22) is a Wannier-Stark Hamiltonian
Hartmann et al. (2004). Thus, the Floquet FSL consists of an infinite set of
(decoupled) 1D chains – each chain represents a specific photon number $n$. We
may note that by adding a photon drive
$\eta\left(\hat{a}+\hat{a}^{\dagger}\right)$, the photon number would no
longer be preserved and the chains would be coupled, resulting in a two square
lattices.

The energies $\varepsilon_{j}$ and eigenstates $|\phi_{j}(n)\rangle$ of the
Hamiltonian can be found analytically Fukuyama et al. (1973)

$\varepsilon_{j}=j\omega,\hskip
28.45274pt|\phi_{j}(n)\rangle=\sum_{m}J_{l-m}\left(\pm\frac{g\sqrt{n}}{\omega}\right)|m\rangle.$
(23)

Here, $J_{m-j}(\pm g\sqrt{n}/\omega)$ is the Bessel function, which is
predominantly localized in the interval $|m-j|<g\sqrt{n}/\omega$, and the
$\pm$-sign stems from the value of $\hat{\tau}_{z}$ (for the Bessel functions
$J_{j}(x)$ one has that they are even/odd for even/odd $j$). Note how the
Wannier-Stark ladder spectrum $\varepsilon_{j}$, i.e. the Floquet quasi
energies, is equidistant and does not depend on the coupling $g\sqrt{n}$. This
is a result from the lattice size extending to infinity, while in numerics,
when a trucnation is imposed, one finds a slight dependence.

Given the spectrum and eigenstates of the Floquet Hamiltonian, it is possible
to derive analytically the time-dependence of various expectations Hartmann et
al. (2004). It is known that the time-evolution generated by the Wannier-Stark
Hamiltonian $\hat{H}_{\mathrm{JC}}^{\mathrm{F}}$ results in periodic Bloch
oscillations if the state is initially localized in momentum space, while if
is is initially localized in real space you find instead periodic breathing
modes where the state repeatedly broadens and contracts. However, we must
recall that our initial state in the $\mathcal{T}$ subspace should be a phase
state $|\theta\rangle$, which is extended over the whole Floquet lattice.
Thereby, such oscillations are not clearly manifested in real space (see,
however, next section where they show up in the numerics when a truncation has
been considered). One the other hand, in reciprocal space a phase state has a
definite ‘momentum’ $k=\theta$. For the Wannier-Stark Hamiltonian, the quasi-
momentum obeys $k(t)=k_{0}-\omega t$ Hartmann et al. (2004), and hence we
expect a linear growth as time evolves. Of course, the quasi-momentum is
limited to the first Brillouin zone $-\pi\leq k<\pi$, and over longer times it
repeatedly sweeps the Brillouin zone. We demonstrate this is fig. 2 by
plotting the ‘momentum’ distribution $P(k,t)$ which is obtained by Fourier
transforming the time-evolved state of the $\mathcal{T}$ space, i.e. we first
need to integrate out the photon and spin/atom degrees of freedom. When
simulating the systems numerically we have to truncate the infinite state
space. The initial state of the photon mode will be a coherent state,
$|\alpha=5\rangle$, and for the atom it will be $|g\rangle$. For the
simulation of the Floquet Hamiltonian we will take the initial state as the
ground state of the “bare” Floquet Hamiltonian
$\hat{H}_{0}^{\mathrm{F}}=\frac{\Delta}{2}\hat{\sigma}_{z}-\left(\hat{E}^{+}+\hat{E}^{-}\right)$,
where we impose periodic boundary conditions of the Floquet space, i.e. the
phase state $|\theta=0\rangle$. Hence, in the reciprocal space this is the
$k=0$ quasi-momentum state. The figure makes clear how the system shows the
characteristic Bloch oscillation behaviour. The effect is somewhat less clear
since we plot the logarithm of the distribution. We do this in order to
envision the effects from the Floquet space truncation, which here show up as
populated regions around the $k=0$ quasi-momentum. Without the logarithm these
are hardly visible by bare eye.

Figure 2: The logarithm of the Floquet ‘momentum’ distribution,
$\log\left[P(k,t)\right]$, corresponding to the simulation of fig. 1.
Initially, the Floquet state is the phase state $|\theta=0\rangle$ or
equivalently the $k=0$ quasi-momentum state, and with time it evolves as
$k=-\omega t$ in accordance to Bloch oscillation. The population of the center
light regions are actually greatly exaggerated by plotting the logarithm; this
is actually a vanishingly small part of the total population. The Floqquet
lattice for this numerics comprises 1401 sites, and by increasing this the
center-most population is suppressed.

While fig. 2 expose the characteristic time evolution of Wannier-Stark
problem, the signatures of the space truncation is only clear from the fact
that the logarithm, rather than the true, of the distribution is presented. To
explore this numerical approximation we show the result of the atomic
inversion, again with a finite Floquet lattice size of size 1401, in fig. 1 to
be compared with the exact result. At short times the agreement is good, and
the initial collapse and first revival is evident. However, for longer times
the shortings of the state space truncation sets in and the finite size
Floquet analysis qualitatively fails. This time can be pushed by increasing
the lattice size, but it rapidly becomes computationally costly.

##  IV Driven Lipkin-Meshkov-Glick model

The Lipkin-Meshkov-Glick model (see also Appendix A), first introduced in the
nuclear physics community as a solvable many-body model Lipkin et al. (1965),
contains a single large spin $S$. The physics involves an interplay between a
term attempting to align the spin in the $z$ direction and a coupling of the
$x$ spin component. Specifically, for a negative ferromagnetic coupling, it
becomes critical in the thermodynamic limit of $S\rightarrow\infty$ Latorre et
al. (2005); Dusuel and Vidal (2004); Ribeiro et al. (2007); Larson (2010). It
can be viewed as a transverse field Ising model with an “all-to-all” coupling,
where each spin couples identically to every other spin. Furthermore, since
the total spin is conserved, every spin symmetry sector realizes a spin-$S$
LMG model. As expected from such a long-range model, the critical behavior
belongs to the mean-field universality class [30].

This model has garnered significant attention in the quantum optics community
Larson and Irish (2017), partly because it results from eliminating the boson
degrees of freedom of either the well-known Tavis-Cummings or the Dicke models
Morrison and Parkins (2008); Larson (2010). It has also been considered when
describing the bosonic Josephson effect Milburn et al. (1997).

The periodically driven LMG model, albeit in a different setting, has
undergone thorough analysis dating back to the 1980s Frahm and Mikeska (1985).
More precisely, the so-called “kicked top” is a specific realization of the
driven LMG model. In the case of periodic kicking, either of the two terms in
the Hamiltonian becomes time-dependent as $\delta(t-n\tau)$ (for an integer
$n$, and some period $\tau$). The delta time dependence allows for analytical
explorations, and the kicked top is particularly interesting as it
demonstrates a system transitioning from regular to chaotic behavior as the
kicking strength is increased Ghose et al. (2008); Chaudhury et al. (2009);
Mumford (2023). It also serves as a platform for studying dynamical
localization, where ergodicity can be lost due to a dynamical analogue of
Anderson localization Haake (1991); Wang and Robnik (2023). These observations
have led to the consideration of the model as a candidate for implementing
time crystals Russomanno et al. (2017).

The harmonically driven LMG model, as considered in this paper, is not as
easily analyzed analytically without invoking any approximations. In the past,
it has been examined in terms of dynamical criticality Engelhardt et al.
(2013). We note that criticality can manifest in various ways in non-
equilibrium situations, as briefly discussed in Bento et al. (2023) in terms
of the LMG model.

While there are various versions of the LMG model, we will focus on the driven
LMG model given by the following Hamiltonian:

$\hat{H}_{\mathrm{LMG}}=\Delta\cos(\omega
t)\hat{S}_{z}-\frac{\lambda}{S}\hat{S}_{x}^{2}.$ (24)

Here, $\Delta$ represents the effective drive amplitude, $\omega$ is the drive
frequency as before, $\lambda$ is the coupling strength, and $S$ is the total
spin. The spin operators obey the regular commutators
$\left[\hat{S}_{i},\hat{S}_{j}\right]=i\varepsilon_{ijk}\hat{S}_{k}$, with
$\varepsilon_{ijk}$ being the Levi-Chivita antisymmetric tensor.

In contrast to the driven JC model discussed in the previous section, here the
drive is applied to the “field” and not to the coupling. The simple reason for
this choice is that it leads to more interesting dynamical structures. In ref.
Engelhardt et al. (2013) (similarly in Bastidas et al. (2012); Dasgupta et al.
(2015); Das et al. (2023) for the Dicke model), the periodically driven LMG
model was considered in a regime where a RWA was assumed. The RWA neglects
rapidly rotating terms, analogous to a time averaging over the short time
scale, and the effective Hamiltonian becomes time-independent. Within this
regime, different phases were identified by introducing a mean-field energy
and analyzing its local minima (which correspond to meta-stable states).

Figure 3: Time-averaged magnitization $\langle\hat{S}_{z}\rangle$ (a) and
phase space localization $PR$ (b) of the periodically driven LMG model (24) as
a function of the field strength $\Delta$ and drive frequency $\omega$.
Several features stand out, such as high localization in the state for rapid
and strong driving, with exceptions highlighted for example by the point A. A
slow driving frequency tends to thermalize the system, similarly for a small
$\Delta.$. The two dots mark the parameter values used in the following figure
5. The total spin $S=10$, and the system is initialized in the $|S,-S\rangle$
state.

We will not impose any RWA, and thereby not derive an effective time-
independent Hamiltonian, but analyze the full problem numerically. The
quantities of interest will typically be time-averaged, for example the
magnetization takes the form

$S_{z}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\langle\psi(t)|\hat{S}_{z}|\psi(t)\rangle
dt.$ (25)

Another quantity of interest will be the phase space “partition ratio”

$P\\!R=\int Q^{2}(\theta,\phi)\,d\theta d\phi.$ (26)

Here, $Q(\theta,\phi)=\langle\theta,\phi|\hat{\rho}|\theta,\phi\rangle$ is the
Husimi $Q$ function, and
$|\theta,\phi\rangle=\exp\left[i\theta\left(\hat{S}_{x}\cos\phi-\hat{S}_{y}\sin\phi\right)\right]|S,S\rangle$
is a spin-coherent state with $|S,S\rangle$ a Dicke state Zhang et al. (1990).
Often the inverse partition, $1/P$ ratio is instead considered, and more
common is to define the partition ratio from a state $|\psi\rangle$444For some
basis $\\{|\varphi_{n}\rangle\\}$, the partition ratio would be
$P=\sum_{n}|\langle\varphi_{n}|\psi\rangle|^{4}$. rather than the phase space
distribution. $P$, as defined here, serves as a measure of phase space
localization Evers and Mirlin (2008). The closer its value is to unity, the
more localized the state, while a smaller value indicates greater
delocalization in phase space. For completeness, it is worth noting that the
partition ratio is expressed in terms of the square of the Husimi
distribution. However, other powers could also be considered, as is done when
analyzing multi-fractal properties. Additionally, $P$ is related to the Rényi-
Wehrl entropy, which has undergone thorough analysis in terms of the Dicke
model in Romera et al. (2012).

Figure 4: Fig. 4. Upper plot (a): two examples of the Husimi distributions
(more precisely $\log\left[Q(\theta,\phi)\right]$ to demonstrate localization)
at a given time $T=200$, corresponding to points A and B in fig. 3. The final
time $T$ influences the Husimi distributions, but the level of localization
remains roughly constant. In the left plot (point A), the partition ratio is
low, and the distribution is more delocalized (sharp red) than in the right
plot (point B). The total spin is as before $S=10$. Lower plot (b): three
snapshots ($T=1,\,2,\,10$) of $\log\left[Q(\theta,\phi)\right]$ in the
$\theta\phi$-plane, for parameters $\Delta=\omega=1/2$ and $S=25$. Initially,
the Husimi distribution has no zeros (the $|S,-S\rangle$ state), but as the
state evolves, a string of $Z=2S$ zeros develops. For short times, the zeros
(blue or yellow dots) follow nice curves, but already at $t=10$, the zeros are
more or less randomly distributed, indicating chaotic behavior.

In fig. 3, we present the numerical results for the time-averaged
magnetization $S_{z}$ (a) and the time-averaged partition ratio $P\\!R$ (b) as
functions of both the drive frequency $\omega$ and field amplitude $\Delta$.
We numerically integrated the Schrödinger equation, as generated by
$\hat{H}_{\mathrm{LMG}}$, over extended times (sufficient for convergence) and
averaged the results. The plots illustrate a direct correlation between
magnetization and phase space localization. For small $\omega$, the system is
ergodic, and its state populates large fractions of phase space. Large
delocalization is also observed for large $\omega$’s and moderate $\Delta$’s,
but this does not appear to be a chaotic regime, as indicated by the regular
pattern of $P\\!R$ for these parameters. For large $\omega$ and $\Delta$, the
system is highly magnetized and localized to its phase space south pole.

It is observed that small parameter fluctuations can significantly alter the
behavior of the evolution, as exemplified by the two points A and B in fig. 3.
To take a closer look at these points, we plot the logarithm of the Husimi
distributions $Q(\theta,\phi)$ in fig. 4 (a) for these two sets of parameters.
The sharper red regions are more extended for case A, indicating a more
delocalized state as anticipated. As time progresses, and given a sufficiently
long duration, the distributions change, while the number of zeros $Z$ remains
constant, and the $P\\!R$ stays more or less constant as well.

A striking feature of the shown Husimi distributions is its zeros. Analysis of
the zeros of the Husimi distribution has a rather long history Leboeuf and
Voros (1990); Cibils et al. (1992); Korsch et al. (1997); Romera et al.
(2012). It is related to the stellar representation of quantum states, i.e.
the zeros uniquely determines the state. For spin states it is related to the
Majorana representation Serrano-Ensástiga and Braun (2020). It is known that
distribution of these zeros and the number $Z$ of them can be considered as an
indicator of chaotic evolution. For regular motion the zeros tend to align in
regular structures, while if the system displays chaos the zeros appear
randomly spread over phase space. For spin $S$ states, the Majorana
representation of a state is a $2S$ polynomial, from which it follows that the
Husimi distribution can at most have $Z=2S$ zeros Amiet and Weigert (1999);
Serrano-Ensástiga and Braun (2020); Chabaud et al. (2020).

In fig. 4 (b), we illustrate how these zeros evolve for one particular Husimi
distribution. Initially, there are no zeros, and as the state evolves, two
lines of zeros develops from the south pole. These curves of zeros twist, and
rapidly the zeros depart from the line, seemingly randomly spreading over the
full phase space. This clearly demonstrates chaotic behavior. We may mention
(although not shown) that if we let $\omega=0$, we regain the time-independent
LMG model, which is known to be integrable. In that case, one also observes
the buildup of $Z=2S$ zeros, but this time they appear in regular patterns.

Figure 5: Time evolution of the magnetization
$S_{z}(t)=\langle\hat{S}_{z}\rangle$ for the same initial state as in the
previous figures, and this time for the dimensionless parameters $\Delta=20$,
$\omega=0.05$, and $S=10$. At periodic times, determined by $\cos(\omega
t)=0$, the magnetization suddenly make a jump and get “stuck” at a random new
magnetization where it rapidly oscillates around, as explicitly shown in the
inset.

In sec. III, we discussed how the mixing of different time scales alters the
evolution. Naturally, the same applies here. For very slow driving, $\omega\ll
1$, we may encounter adiabatic time evolution as long as $S$ is not too large.
The criticality of the LMG model implies that adiabaticity is doomed in the
thermodynamic limit since the system is gapless at the critical points. For
rapid driving, the evolution is approximately diabatic. This is the regime in
the upper right corner of fig. 3. When $|\Delta|\gg\omega$, the intrinsic LMG
time scale is short in comparison to that of the drive. Here, the chaotic
evolution manifests in a generic feature, as we demonstrate in fig. 5 by
plotting the time evolution of the magnetization $S_{z}$ (hence, not time-
averaged). For $\cos(\omega t)$ away from zero, as the inset shows, the
magnetization oscillates with a frequency set by the instantaneous “field
strength” $\Delta\cos(\omega t)S$. Whenever $\cos(\omega t)\rightarrow 0$, the
magnetization makes a jump to some random new magnetization, where it starts
to oscillate around. For a very large spin $S$, a simple numerical exploration
suggests that these jumps are indeed random 555We integrated the semi-
classical equations of motion over long periods and collected the
magnetizations where the system get stuck. Up to some boundary deviations, the
distribution of these was found to be evenly distributed within the interval
$[-S,S]$.. Thus, the system can serve as a chaos-driven random number
generator Yu et al. (2019). The jumps originate from multi-level Landau-Zener
transitions Landau (1932); Zener (1932); Larson (2014). It is known that
Landau-Zener transitions are highly sensitive to dynamical phases, and the
more involved levels this will generate seemingly random transitions. It is
interesting to note that the same random structure was found in an
optomechanical system in Larson and Horsdal (2011). In that survey, the two
time scales were both intrinsic, while in the present model, one is set by the
classical drive.

Figure 6: The level statistics of the Floquet Hamiltonian (27), for a Floquet
lattice of 201 sites, $\Delta=2$, $\omega=0.25$, and $S=10$. Although the
distribution is neither Poissonian nor of Wigner-Dyson shape, it clearly
exhibits level repulsion.

Up till now, all numerical simulations have been performed for the time-
dependent Hamiltonian (24). We have not analyzed the problem from the
perspective of its Floquet Hamiltonian, which in this case becomes

$\hat{H}_{\mathrm{LMG}}^{\mathrm{F}}=\omega\hat{E}_{0}+\frac{\Delta}{2}\left(\hat{E}^{+}+\hat{E}^{-}\right)\hat{S}_{z}-\frac{\lambda}{S}\hat{S}_{x}^{2}.$
(27)

In terms of a Floquet FSL, see the appendix, this Hamiltonian represents a
tilted square lattice. One dimension corresponds to the Floquet space, and the
other, which is finite, to the spin space. The clear advantage of the Floquet
Hamiltonian is that it is time-independent. However, it increases the state
space, leading to truncation approximation in numerical computations.

It actually turns out to be far more computationally costly to solve the
Floquet problem than the time-dependent one. The main reason for this is that
one must consider very large systems in order to correctly capture the exact
dynamics, which is understood from the fact that our initial state is a phase
state that is delocalized in the Floquet FSL. Nevertheless, insight into the
physics can be gained by analyzing the Floquet Hamiltonian.

One question that comes to mind is how the chaos is manifested in
$\hat{H}_{\mathrm{LMG}}^{\mathrm{F}}$. The Hamiltonian bears similarities with
the Dicke model, which is known to display a transition from regular to
chaotic evolution across the normal-superradiant phase transition Emary and
Brandes (2003). What differs between $\hat{H}_{\mathrm{LMG}}^{\mathrm{F}}$ and
the Dicke Hamiltonian is that here a spin couples to the ‘Euclidean’ field
(represented by the $\hat{E}$-operators) rather than a boson field, and the
bare spin Hamiltonian is quadratic ($\hat{S}_{X}^{2}$ instead of
$\hat{S}_{x}$). Chaos in the Dicke model has been analyzed in various ways,
for example, by the level statistics Haake (1991).

Given the spectrum $E_{n}$ (ordered $E_{1}\leq E_{2}\leq E_{2}\leq\dots$) we
introduce the level spacings $S_{n}=E_{n+1}-E_{n}$. With $S_{n}$, one
constructs the normalized distribution function $P(S)$666The level spacing $S$
should not be confused with the spin $S$.. In case of an integrable
Hamiltonian we expect a Poisson distribution

$P_{\mathrm{P}}(S)=e^{-S},$ (28)

while for chaotic motion we expect a Wigner surmise

$P_{\mathrm{W}}(S)=a_{\beta}S^{\beta}e^{-b_{\beta}S^{2}},$ (29)

where $a_{\beta}$ and $b_{\beta}$ are constants and the exponent $\beta$ is
set by the symmetry classes Haake (1991). The fact that $\lim_{S\rightarrow
0}P_{\mathrm{W}}(S)=0$ implies that that the energy levels do not cross,
something referred to as level repulsion.

Taking the Wigner surmise as a chaos indicator, one should note that fully
developed chaos is not to be expected for low energies. Therefore, one
typically considers intermediate energies. For a finite-dimensional system,
the central part of the spectrum is used. Furthermore, the spectrum is
typically unfolded before calculating the distribution function. This involves
normalizing the level distances $S$ relative to a mean distance Haake (1991).
The procedure of unfolding follows these steps:

  1. 1.

Compute $N(E)$, which counts the number of states with energies smaller or
equal to $E$.

  2. 2.

Fit a smooth curve $R(E)$ to the step-like function $N(E)$ (we use a seventh
order polynomial).

  3. 3.

Identifies energies $e_{n}=R(E_{n})$, and calculate the corresponding
‘unfolded’ distribution function $P(s)$.

An example of the distribution function P(s)P(s) is shown in fig. 6. It
demonstrates level repulsion, as one would expect for chaotic evolution, but
it is not fully developed. Part of this could be explained by the fact that we
consider the full spectrum without removing degeneracies arising from
symmetries. The Hamiltonian $\hat{H}_{\mathrm{LMG}}^{\mathrm{F}}$ possess two
$\mathcal{Z}_{2}$-parity symmetries,
$(\hat{E}_{0},\hat{E}^{\pm},\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})\rightarrow(\hat{E}_{0},-\hat{E}^{\pm},-\hat{S}_{x},\hat{S}_{y},-\hat{S}_{z})$
and
$(\hat{E}_{0},\hat{E}^{\pm},\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})\rightarrow(\hat{E}_{0},-\hat{E}^{\pm},\hat{S}_{x},-\hat{S}_{y},-\hat{S}_{z})$.
We should also mention that the actual shape of $P(s)$ depends strongly on the
system parameters, and sometimes the distribution looks more Poisson and
sometimes nothing like Poisson or Wigner.

Figure 7: Time evolution of the magnitization $\langle\hat{S}_{z}\rangle$ (a)
and $\langle\hat{E}_{0}\rangle$ (b) for parameters corresponding to the case A
(red curve) and B (blue curve) of fig. 3. In case B, the system remains highly
magnetized, while in case A, the magnetization oscillates over large periods.
In the lower plot, it is evident that in the Floquet space, the two cases
behave qualitatively differently. Both curves show rapid oscillations
corresponding to Bloch oscillations, but the red curve exhibits collapses and
revivals on top of these. The number of lattice sites (states in the Floquet
space) is taken to be 101.

Let us conclude our explorations by revisiting points A and B in fig. 3. Our
initial Floquet state is delocalized over the full lattice, i.e.,
$\langle\hat{E}_{0}\rangle=0$. However, our truncation implies that
$\langle\hat{E}_{0}\rangle$ may actually display non-trivial evolution, as
demonstrated in Fig. 7. In (a), we display the evolution of the magnetization
$\langle\hat{S}_{z}\rangle$, where a clear difference is observed. For case B,
the system rapidly settles to a certain magnetization and thereafter only
shows small oscillations around it. It resembles rapid thermalization, but we
should remember that the state is actually fairly localized in phase space
here, as seen in fig. 4. One could speculate whether this could be a result of
dynamical localization Haake (1991), but this seems highly unlikely since the
kicked top does not exhibit dynamical localization, and this should hold for a
harmonically driven top as well. For case A, we observe large-amplitude
sustained oscillations spanning over long periods. We believe that these are
instances where constructive interference in the repeated multi-level Landau-
Zener transitions occurs. An indicator supporting this is that by increasing
$S$, it seems as if these ’resonances’ become rarer.

Part of the evolution in the Floquet space is captured in fig. 3 (b) by
showing the evolution of the expectation $\langle\hat{E}_{o}\rangle$ for the
two cases. Apart from short oscillations resulting from the tilt, i.e.
manifestations of Bloch oscillations, case A also displays a collapse-revival
pattern. For one oscillation of $\langle\hat{S}_{z}\rangle$,
$\langle\hat{E}_{0}\rangle$ makes two oscillations. We should recall that the
Floquet Hamiltonian is truncated, and at these long times we cannot expect
convergence of the results. But even so, it is clear that the two cases behave
very different.

##  V Conclusion

We have examined two periodically driven models, the JC and LMG models, which
frequently appear in the quantum optics community. In both cases, the dynamics
strongly depend on the involved time scales: those set by the drive and the
intrinsic ones. In the JC model, multiple intrinsic time scales exist, such as
Rabi oscillations, collapse, and revival times. When the time scale of the
drive is comparable to the Rabi time, a rather complex dynamics emerges.
Conversely, if this time scale is comparable to the revival time, the dynamics
can be understood from an adiabatic perspective. In such cases, the drive can
act as an effective reversal of time, akin to spin-echoes.

In our Floquet analysis of the problem, we introduce a Floquet FSL, which are
generalizations of FSLs for time-independent models Wang et al. (2016); Cai
and Wang (2021); Saugmann and Larson (2023); Yuan et al. (2024). The time-
independent Floquet Hamiltonian determines the Floquet lattice in the same way
as for regular FSLs, with the Floquet space adding an additional dimension to
the lattice. For the JC model, counterparts of Bloch oscillations are found in
the Floquet FSL.

The LMG model becomes chaotic when a periodic drive is included, giving rise
to several novel features. While the undriven LMG model is not chaotic, it is
not immediately clear whether the driven one should generally lead to heating
and ergodicity. Indeed, the kicked LMG has been studied in terms of Floquet
time-crystals, which evade heating up Russomanno et al. (2017). Despite the
states being fairly localized in phase space, we do observe several signatures
of chaos, such as the distribution of zeros of the Husimi function and random
spin jumps. Many of these results can be analyzed in terms of repeated multi-
level Landau-Zener transitions, where strong chaos leading to ergodicity
arises from complex multi-level interferences. The random jumps occur for slow
driving, but not during fully adiabatic evolution. Additionally, the
interferences may lead to ‘resonances’ where the system behaves qualitatively
differently for slight parameter changes. These resonances also manifest in
the Floquet FSL in the form of collapses and revivals.

We have limited our analysis to two models, but one could envision exploring
many other models as well. It would be intriguing to investigate different
Floquet FSL geometries. In principle, by employing two harmonic drives with
incommensurate frequencies, it is possible to create Floquet FSLs with two
additional dimensions, one for each drive. However, a drawback of the Floquet
FSL, compared to regular FSLs, is that the initial states are delocalized.
This aspect makes them less practical for ”lattice engineering” and
simulations involving synthetic dimensions.

The LMG model is integrable, unlike the Dicke model, for instance.
Consequently, one might anticipate qualitatively different behaviors for a
driven Dicke model. For example, it could exhibit dynamical localization,
which differs from the localization discussed in Wang and Robnik (2020).

##  Appendix A Fock state lattices

In this appendix we briefly discuss Fock state lattices and consider it in
terms of the LMG model. The FSL for generalized JC models was first introduced
in Wang et al. (2016). Each Fock state represents a quantum state of one
particular lattice site. Tunneling between the synthetic lattice sites occurs
via the action of the ladder operators, but due to the square-root
normalization, e.g. for bosons $\hat{a}|n\rangle=\sqrt{n}|n\rangle$, the FSL
will typically not be translationally invariant even within the bulk.
Nevertheless, as demonstrated in Cai and Wang (2021) for the three-mode JC
model and in Saugmann and Larson (2023) for various other setups, lattice
properties like topology may survive. Most recently, the topological aspects
has been experimentally explored Deng et al. (2022).

Characteristic of many models of light-matter interaction, as found in the
quantum optics setting, is that they are typically restricted to few degrees
of freedom. The JC model is a prime example Jaynes and Cummings (1963); Shore
and Knight (1993); Larson and Mavrogordatos (2022). Furthermore, stemming from
dipole interactions, the interactions are usually linear in the boson and spin
degrees of freedom, meaning that in the Fock basis (15), the Hamiltonian
becomes sparse with a structure similar to tight-binding lattice models. In
Saugmann and Larson (2023), this was discussed for numerous models, and it was
demonstrated that most known light-matter interaction models in the quantum
optics community resulted in FSLs with the same geometry as well-known lattice
models found in solid-state physics.

For clarity, let us explore the time-independent example of sec. IV, namely
the LMG Hamiltonian (24) that takes the form

$\hat{H}_{\mathrm{LMG}}=\Delta\hat{S}_{z}-\frac{\lambda}{S}\hat{S}_{x}^{2}.$
(30)

The model supports a $\mathbb{Z}_{2}$-parity symmetry characterized by
$\left(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z}\right)\rightarrow\left(-\hat{S}_{x},-\hat{S}_{y},\hat{S}_{z}\right)$,
and furthermore the total spin $S$ is preserved. The denominator in the last
term warrants that the two terms scale as $\mathcal{O}(S)$ in the large $S$
limit. As mentioned in the main text, in the thermodynamic limit,
$S\rightarrow\infty$, the model is critical at the coupling strength
$\lambda_{c}=|\Delta|$; below the critical point the ground state is
ferromagnetic with $S_{x}=\langle\hat{S}_{x}\rangle=0$ (and $S_{z}\neq 0$),
while above it is paramagnetic $S_{x}\neq 0$ Dusuel and Vidal (2004); Ribeiro
et al. (2007); Larson (2010); Latorre et al. (2005).

The Fock states, according to (15), are $|m\rangle$ with $-S\leq m\leq S$, and

$\begin{array}[]{l}\hat{S}_{z}|m\rangle=m|m\rangle,\\\
\hat{S}_{x}|m\rangle=\alpha_{m+}|m+1\rangle+\alpha_{m-}|m-1\rangle,\end{array}$
(31)

with the normalization $\alpha_{m\pm}=\sqrt{S(S+1)-m(m\pm 1)}$ (note that
hermiticity implies $\alpha_{m+}=\alpha_{m+1-}$). The parity symmetry implies
that the Hamiltonian couples states $|m\rangle$ with odd (even) $m$ to new
states $|m^{\prime}\rangle$ with odd (even) $m^{\prime}$. Thus, we expect the
FSL to be comprised of two 1D chains, as visualized in fig. 8. In the figure,
we do not emphasize the onsite energies set by the diagonal terms of the
Hamiltonian. In particular, the field $h\hat{S}_{z}$ induces a linear tilt,
$hm$, of the lattice, while the interaction term causes an onsite shift
$\alpha_{m+}^{2}+\alpha_{m-}^{2}=2S(S+1)+2m^{2}$. If the
$\mathbb{Z}_{2}$-symmetry is broken, by adding for example a field
$\lambda\hat{S}_{x}$ to the Hamiltonian (30), the two chains would be coupled
and one instead finds a triangular ladder.

Figure 8: The FSL of the LMG model (30). The parity symmetry divide the
lattice into two 1D chains, one of even parity comprised of all states
$|m\rangle$ with even quantum numbers $m$, and one with odd parity for the odd
$m$’s. The blue lines connecting neighboring sites represent non-vanishing
tunnelings. The corresponding tunneling rates $t_{m}$ are identified with the
square-root prefactors given in eq. (31), i.e.
$t_{m}=\alpha_{m+}\alpha_{m+1+}$.

For harmonically driven systems, the method outlined in this work allows the
introduction of a Floquet FSL. In particular, we identify the Fourier modes
$m$ with lattice sites, i.e. $|m\rangle$ represents the state for a particle
residing on site $m$. Hence, $\cos(\omega t)$ translates into a tight-binding
Hamiltonian which induces hopping between neighbouring lattice sites. As a
most simple example we can take the Rabi model with its corresponding
Hamiltonian

$\hat{H}_{\mathrm{R}}(t)=\Omega\hat{\sigma}_{z}+g\cos(\omega
t)\hat{\sigma}_{x}.$ (32)

Here, $\hat{\sigma}_{\alpha}$ ($\alpha=x,\,y,\,z$) are the Pauli matrices,
$\Omega$ is the transition frequency of our two-level “atom”, and $g$ is the
amplitude of the periodic drive. The Floquet Hamiltonian takes the form Duan
(2023)

$\hat{H}_{\mathrm{R}}^{\mathrm{F}}=\omega\hat{E}_{0}+\Omega\hat{\sigma}_{z}+\frac{g}{2}\left(\hat{E}^{+}+\hat{E}^{-}\right)\hat{\sigma}_{x}.$
(33)

Following the notations of (15), the Fock states are $|s,m\rangle$. It is
possible to give a physical picture of the various terms forming the
Hamiltonian

  * •

The first term $\omega\hat{E}_{0}$ is diagonal in the Fock basis and thereby
only shifts the onsite ‘energies’. This shift is linear and can be seen as a
constant “electric” field for a charged particle, akin to the case of Bloch
oscillations Dahan et al. (1996); Morandotti et al. (1999); Holthaus (2000) or
the Wannier-Stark ladder problem Wilkinson et al. (1996); Glück et al. (1998);
Hartmann et al. (2004).

  * •

The second term $\Omega\hat{\sigma}_{z}$ is again diagonal and provides an
onsite shift depending on the internal spin state.

  * •

The last term
$\frac{g}{2}\left(\hat{E}^{+}+\hat{E}^{-}\right)\hat{\sigma}_{x}$ represents
the nearest neighbor tunneling between the $m$-sites. The tunneling is,
however, accompanied by a spin-flip.

## References

  * Jaynes and Cummings (1963) Edwin T Jaynes and Frederick W Cummings.  Comparison of quantum and semiclassical radiation theories with application to the beam maser.  _Proceedings of the IEEE_ , 51(1):89–109, 1963.
  * Shore and Knight (1993) Bruce W Shore and Peter L Knight.  The jaynes-cummings model.  _Journal of Modern Optics_ , 40(7):1195–1238, 1993.
  * Larson and Mavrogordatos (2022) Jonas Larson and Th K Mavrogordatos.  The jaynes-cummings model and its descendants.  _IOP ebooks; ISBN: 978-0-7503-3447-1 arXiv:2202.00330_ , 2022.
  * Haroche and Raimond (2006) Serge Haroche and J-M Raimond.  _Exploring the quantum: atoms, cavities, and photons_.  Oxford university press, 2006.
  * Schlicher (1989) Rainer R Schlicher.  Jaynes-cummings model with atomic motion.  _Optics communications_ , 70(2):97–102, 1989.
  * Lipkin et al. (1965) Harry J Lipkin, N Meshkov, and AJ Glick.  Validity of many-body approximation methods for a solvable model:(i). exact solutions and perturbation theory.  _Nuclear Physics_ , 62(2):188–198, 1965.
  * Dusuel and Vidal (2004) Sébastien Dusuel and Julien Vidal.  Finite-size scaling exponents of the lipkin-meshkov-glick model.  _Physical review letters_ , 93(23):237204, 2004.
  * Ribeiro et al. (2007) Pedro Ribeiro, Julien Vidal, and Rémy Mosseri.  Thermodynamical limit of the lipkin-meshkov-glick model.  _Physical review letters_ , 99(5):050402, 2007.
  * Morrison and Parkins (2008) S Morrison and AS Parkins.  Dynamical quantum phase transitions in the dissipative lipkin-meshkov-glick model with proposed realization in optical cavity qed.  _Physical review letters_ , 100(4):040403, 2008.
  * Larson (2010) Jonas Larson.  Circuit qed scheme for the realization of the lipkin-meshkov-glick model.  _Europhysics letters_ , 90(5):54001, 2010.
  * Latorre et al. (2005) José I Latorre, Román Orús, Enrique Rico, and Julien Vidal.  Entanglement entropy in the lipkin-meshkov-glick model.  _Physical Review A_ , 71(6):064101, 2005.
  * Holthaus (2015) Martin Holthaus.  Floquet engineering with quasienergy bands of periodically driven optical lattices.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ , 49(1):013001, 2015.
  * Eckardt (2017) André Eckardt.  Colloquium: Atomic quantum gases in periodically driven optical lattices.  _Reviews of Modern Physics_ , 89(1):011004, 2017.
  * Rudner and Lindner (2020) Mark S Rudner and Netanel H Lindner.  The floquet engineer’s handbook.  _arXiv preprint arXiv:2003.08252_ , 2020.
  * Shirley (1965) Jon H Shirley.  Solution of the schrödinger equation with a hamiltonian periodic in time.  _Physical Review_ , 138(4B):B979, 1965.
  * Klimov and Chumakov (2009) Andrei B Klimov and Sergei M Chumakov.  _A group-theoretical approach to quantum optics: models of atom-field interactions_.  John Wiley & Sons, 2009.
  * Prants and Yacoupova (1992) SV Prants and LS Yacoupova.  The jaynes-cummings model with modulated field-atom coupling in resonator quantum electrodynamics.  _Journal of Modern Optics_ , 39(5):961–971, 1992.
  * Joshi and Lawande (1993) Amitabh Joshi and SV Lawande.  Generalized jaynes-cummings models with a time-dependent atom-field coupling.  _Physical Review A_ , 48(3):2276, 1993.
  * Meystre et al. (1973) P Meystre, E Geneux, A Faist, and A Quattropani.  Destruction of coherence by scattering of radiation on atoms.  _Lettere al Nuovo Cimento (1971-1985)_ , 6(8):287–291, 1973.
  * Narozhny et al. (1981) NB Narozhny, JJ Sanchez-Mondragon, and JH Eberly.  Coherence versus incoherence: Collapse and revival in a simple quantum model.  _Physical Review A_ , 23(1):236, 1981.
  * Saugmann and Larson (2023) Pil Saugmann and Jonas Larson.  Fock-state-lattice approach to quantum optics.  _Physical Review A_ , 108(3):033721, 2023.
  * Hartmann et al. (2004) Timo Hartmann, F Keck, HJ Korsch, and S Mossmann.  Dynamics of bloch oscillations.  _New Journal of Physics_ , 6(1):2, 2004.
  * Fukuyama et al. (1973) Hidetoshi Fukuyama, Robert A Bari, and Hans C Fogedby.  Tightly bound electrons in a uniform electric field.  _Physical Review B_ , 8(12):5579, 1973.
  * Larson and Irish (2017) Jonas Larson and Elinor K Irish.  Some remarks on ‘superradiant’phase transitions in light-matter systems.  _Journal of Physics A: Mathematical and Theoretical_ , 50(17):174002, 2017.
  * Milburn et al. (1997) GJ Milburn, J Corney, Ewan M Wright, and DF Walls.  Quantum dynamics of an atomic bose-einstein condensate in a double-well potential.  _Physical Review A_ , 55(6):4318, 1997.
  * Frahm and Mikeska (1985) H Frahm and HJ Mikeska.  On the dynamics of a quantum system which is classically chaotic.  _Zeitschrift für Physik B Condensed Matter_ , 60(1):117–126, 1985.
  * Ghose et al. (2008) Shohini Ghose, Rene Stock, Poul Jessen, Roshan Lal, and Andrew Silberfarb.  Chaos, entanglement, and decoherence in the quantum kicked top.  _Physical Review A_ , 78(4):042318, 2008.
  * Chaudhury et al. (2009) S Chaudhury, A Smith, BE Anderson, S Ghose, and Poul S Jessen.  Quantum signatures of chaos in a kicked top.  _Nature_ , 461(7265):768–771, 2009.
  * Mumford (2023) J Mumford.  Many topological regions on the bloch sphere of the spin-1/2 double-kicked top.  _Physical Review A_ , 107(5):053316, 2023.
  * Haake (1991) Fritz Haake.  _Quantum signatures of chaos_.  Springer, 1991.
  * Wang and Robnik (2023) Qian Wang and Marko Robnik.  Statistics of phase space localization measures and quantum chaos in the kicked top model.  _Physical Review E_ , 107(5):054213, 2023.
  * Russomanno et al. (2017) Angelo Russomanno, Fernando Iemini, Marcello Dalmonte, and Rosario Fazio.  Floquet time crystal in the lipkin-meshkov-glick model.  _Physical Review B_ , 95(21):214307, 2017.
  * Engelhardt et al. (2013) G Engelhardt, VM Bastidas, C Emary, and T Brandes.  ac-driven quantum phase transition in the lipkin-meshkov-glick model.  _Physical Review E_ , 87(5):052110, 2013.
  * Bento et al. (2023) Pedro HS Bento, Adolfo del Campo, and Lucas C Céleri.  Krylov complexity and dynamical phase transition in the quenched lmg model.  _arXiv preprint arXiv:2312.05321_ , 2023.
  * Bastidas et al. (2012) VM Bastidas, C Emary, B Regler, and T Brandes.  Nonequilibrium quantum phase transitions in the dicke model.  _Physical review letters_ , 108(4):043003, 2012.
  * Dasgupta et al. (2015) Sayak Dasgupta, Utso Bhattacharya, and Amit Dutta.  Phase transition in the periodically pulsed dicke model.  _Physical Review E_ , 91(5):052129, 2015.
  * Das et al. (2023) Pragna Das, Devendra Singh Bhakuni, Lea F Santos, and Auditya Sharma.  Periodically and quasiperiodically driven anisotropic dicke model.  _Physical Review A_ , 108(6):063716, 2023.
  * Zhang et al. (1990) Wei-Min Zhang, Robert Gilmore, et al.  Coherent states: Theory and some applications.  _Reviews of Modern Physics_ , 62(4):867, 1990.
  * Evers and Mirlin (2008) Ferdinand Evers and Alexander D Mirlin.  Anderson transitions.  _Reviews of Modern Physics_ , 80(4):1355, 2008.
  * Romera et al. (2012) E Romera, R Del Real, and M Calixto.  Husimi distribution and phase-space analysis of a dicke-model quantum phase transition.  _Physical Review A_ , 85(5):053831, 2012.
  * Leboeuf and Voros (1990) P Leboeuf and André Voros.  Chaos-revealing multiplicative representation of quantum eigenstates.  _Journal of Physics A: Mathematical and General_ , 23(10):1765, 1990.
  * Cibils et al. (1992) MB Cibils, Y Cuche, P Leboeuf, and WF Wreszinski.  Zeros of the husimi functions of the spin-boson model.  _Physical Review A_ , 46(8):4560, 1992.
  * Korsch et al. (1997) HJ Korsch, C Müller, and H Wiescher.  On the zeros of the husimi distribution.  _Journal of Physics A: Mathematical and General_ , 30(20):L677, 1997.
  * Serrano-Ensástiga and Braun (2020) E Serrano-Ensástiga and D Braun.  Majorana representation for mixed states.  _Physical Review A_ , 101(2):022332, 2020.
  * Amiet and Weigert (1999) Jean-Pierre Amiet and Stefan Weigert.  Coherent states and the reconstruction of pure spin states.  _Journal of Optics B: Quantum and Semiclassical Optics_ , 1(5):L5, 1999.
  * Chabaud et al. (2020) Ulysse Chabaud, Damian Markham, and Frédéric Grosshans.  Stellar representation of non-gaussian quantum states.  _Physical Review Letters_ , 124(6):063605, 2020.
  * Yu et al. (2019) Fei Yu, Lixiang Li, Qiang Tang, Shuo Cai, Yun Song, and Quan Xu.  A survey on true random number generators based on chaos.  _Discrete Dynamics in Nature and Society_ , 2019:1–10, 2019.
  * Landau (1932) Lev Landau.  Zur theorie der energieubertragung. ii.  _Physikalische Zeitschrift der Sowjetunion_ , 2:46, 1932.
  * Zener (1932) Clarence Zener.  Non-adiabatic crossing of energy levels.  _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ , 137(833):696–702, 1932.
  * Larson (2014) Jonas Larson.  Interaction-induced landau-zener transitions.  _Europhysics Letters_ , 107(3):30007, 2014.
  * Larson and Horsdal (2011) Jonas Larson and Mats Horsdal.  Photonic josephson effect, phase transitions, and chaos in optomechanical systems.  _Physical Review A_ , 84(2):021804, 2011.
  * Emary and Brandes (2003) Clive Emary and Tobias Brandes.  Chaos and the quantum phase transition in the dicke model.  _Physical Review E_ , 67(6):066203, 2003.
  * Wang et al. (2016) Da-Wei Wang, Han Cai, Ren-Bao Liu, and Marlan O Scully.  Mesoscopic superposition states generated by synthetic spin-orbit interaction in fock-state lattices.  _Physical review letters_ , 116(22):220502, 2016.
  * Cai and Wang (2021) Han Cai and Da-Wei Wang.  Topological phases of quantized light.  _National Science Review_ , 8(1):nwaa196, 2021.
  * Yuan et al. (2024) Jiale Yuan, Han Cai, and Da-Wei Wang.  Quantum simulation in fock-state lattices.  _Advances in Physics: X_ , 9(1):2325611, 2024.
  * Wang and Robnik (2020) Qian Wang and Marko Robnik.  Statistical properties of the localization measure of chaotic eigenstates in the dicke model.  _Physical Review E_ , 102(3):032212, 2020.
  * Deng et al. (2022) Jinfeng Deng, Hang Dong, Chuanyu Zhang, Yaozu Wu, Jiale Yuan, Xuhao Zhu, Feitong Jin, Hekang Li, Zhen Wang, Han Cai, et al.  Observing the quantum topology of light.  _Science_ , 378(6623):966–971, 2022.
  * Duan (2023) Liwei Duan.  Periodic jumps in binary lattices with a static force.  _Physical Review B_ , 108(17):174306, 2023.
  * Dahan et al. (1996) Maxime Ben Dahan, Ekkehard Peik, Jakob Reichel, Yvan Castin, and Christophe Salomon.  Bloch oscillations of atoms in an optical potential.  _Physical Review Letters_ , 76(24):4508, 1996.
  * Morandotti et al. (1999) R Morandotti, U Peschel, JS Aitchison, HS Eisenberg, and Y Silberberg.  Experimental observation of linear and nonlinear optical bloch oscillations.  _Physical Review Letters_ , 83(23):4756, 1999.
  * Holthaus (2000) Martin Holthaus.  Bloch oscillations and zener breakdown in an optical lattice.  _Journal of Optics B: Quantum and Semiclassical Optics_ , 2(5):589, 2000.
  * Wilkinson et al. (1996) SR Wilkinson, CF Bharucha, KW Madison, Qian Niu, and MG Raizen.  Observation of atomic wannier-stark ladders in an accelerating optical potential.  _Physical review letters_ , 76(24):4512, 1996.
  * Glück et al. (1998) Markus Glück, AR Kolovsky, Hans Jürgen Korsch, and N Moiseyev.  Calculation of wannier-bloch and wannier-stark states.  _The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics_ , 4:239–246, 1998.