# A Fock state lattice approach to quantum optics

Pil Saugmann and Jonas Larson  Department of Physics, Stockholm University,
AlbaNova University Center, 106 91 Stockholm, Sweden

(September 3, 2024)

###### Abstract

We analyze a set of models frequently appearing in quantum optical settings by
expressing their Hamiltonians in terms of Fock-state lattices (FSLs). The few
degrees-of-freedom of such models, together with the system symmetries, make
the emerging FSLs rather simple such that they can be linked to known lattice
models from the condensed matter community. This sheds new light on known
quantum optical systems. While we provide a rather long list of models and
their corresponding FSLs, we pick a few ones in order to demonstrate the
strength of the method. The three-mode boson model, for example, is shown to
display a fractal spectrum, and chiral evolution in the FSL characterized by
localized distributions traversing along symmetric trajectories. In a second
example we consider the central spin model which generates a FSL reminiscent
of the Su-Schrieffer–Heeger model hosting topological edge states. We further
demonstrate how the phenomena of flat bands in lattice models can manifest in
related FSLs, which can be linked to so called dark states.

###### pacs:

42.50.Ct, 42.50.Pq, 63.90.+t

##  I Introduction

In physics, lattice models prove a very powerful tool for the investigation of
condensed matter systems altland . They arise naturally in many situations,
most notably in crystal structures, and at the same time they provide a direct
link to computational physics as all continuous models turn into lattice
models when discretized. In the quantum realm, lattice models are in one way
or another a key ingredient in many suggestions for realization of quantum
simulators mb . While some quantum single particle lattice models can be
solved analytically, others can effectively be simulated on a classical
computer. However, they quickly grow in complexity when more particles are
added to the system even when they are non interacting, no longer allowing for
an analytical solution to be found nor even being numerically simulated.
Indeed, it was recently shown how quadratic bosonic lattice models with only a
few degrees-of-freedom can become computationally hard to simulate, and
furthermore how such models can be mapped to interacting many-body systems
FSsimulation . Thus, alternative physical realizations of the lattice models
are in place. Here, the power of lattice models becomes particularly visible,
as many properties are not linked to details of the specific physical
realization of the lattice model, but instead linked to the underlying
symmetries of the lattice model – the idea of universality. As long as these
stay the same, so do the universal features. Keeping this in mind, it is of
course interesting to look for novel ways to realize both well known and new
lattice models.

Traditionally, the geometry of the lattice is rooted in the system’s spatial
dimensions, i.e. we consider lattices in three or lower dimensions. In
mathematics there are no such limitations and we may even have a non-integer
number of dimensions. Synthetic lattice dimensions occur in various proposals
for quantum simulators as a route towards new exotic models synt ; synt1 ;
synt2 ; synt3 . Examples on this are utilizing atomic internal degrees-of-
freedom synt , or vibrational states of the harmonic trapping potential synt2
, which can be employed to create synthetic dimensions for ultracold atoms
loaded in optical lattices. More generally, if we give up the idea of a
lattice living in real space, any Hamiltonian $\hat{H}$, single-particle or
many-body, can be described as a single-particle system if mapped into state
space. In such a view, a lattice-like structure emerges by identifying the
different basis states as the lattice sites. Then, $\hat{H}$ is a matrix in
which the diagonal terms serve as on-site energies in the lattice, and the
off-diagonal elements give the tunneling rates between these sites. Such a
lattice will rapidly grow in complexity as the number of degrees-of-freedom is
increased, and it might well be that then not much new insight is gained by
such a viewpoint. Furthermore, for practical reasons, we require that the
chosen basis should be physically relevant. Paradigmatic light-matter models
found in quantum optics seem perfect for this; it is possible to isolate a few
relevant degrees-of-freedom, and the natural basis comprises the bare Fock
states. Hence, their Hamiltonians give rise to Fock state lattices (FSLs). The
FSL emerging from the Jaynes-Cummings (JC) model was already introduced in
Ref. FSL1 , and studied further in Ref. jctop  to explore a new type of
topological matter formed by quantized light. We note that, historically,
similar lattice ideas were discussed in terms of so-called Glauber-Fock
lattices in wave-guide systems of classical light synt3 .

In this work we bring forward a plethora of light-matter type models, even
though we shall focus on a few of them. It turns out that these capture the
general ideas, at the same time as the physics is extremely rich despite their
simplicity. Taking the multi-mode JC model as an example, in the high detuned
limit the spin degree-of-freedom freezes out and one is left with an effective
quadratic Hamiltonian consisting only of bosonic degrees-of-freedom. The Fock
states provide a natural basis to study such systems, and the sparseness of
the Hamiltonian makes the FSLs convenient for getting insight into the models.
For three bosonic modes, a triangular FSL emerges, and by further imposing a
synthetic magnetic flux it is shown how the system’s spectrum becomes fractal,
akin to the Hofstadter butterfly spectrum of two dimensional (2D) tight-
binding lattice models exposed to homogeneous perpendicular magnetic fluxes
hof . It should be noted, though, that the FSLs studied in this paper lack
translational invariance, and hence the fractal structure is not deriving from
the same mechanism as for the Hofstadter models. The properties of our fractal
spectrum imply that initial Fock states will follow very symmetric curves
within the FSL. As another example we consider the central spin model in which
a single spin-1/2 particle interacts identically to $N$ other spin-1/2
particles. By tuning the interaction it is possible to create a FSL bearing
similarities to the Su-Schrieffer–Heeger (SSH) model, and like for the SSH
model one finds exponentially localized edge states. Related to this we also
discuss a few models which support a large number of degenerate $E=0$
eiegenstates, which serve as the counterpart of flat bands in for example Lieb
lattices.

The structure of the paper is the following. In the next section we
systematically introduce the FSLs. First for the resonant JC models, and from
there we elaborate on how one may design interesting FSLs by looking at
related models. In Sec. III we discuss what role the underlying symmetries has
for the lattice dimensions. Multi-mode JC models in the large detuned limit,
and in the presence of a synthetic flux, are considered in Sec. IV. For these
systems, we see how a fractal like structure of energy spectrum emerges for
certain fluxes, which motivates us to in Sec. V study how initially localized
states evolve within the FSL. Section VI discusses a connection between the
central spin model and the SSH model, namely the appearance of exponentially
localized edge sates in the FSL. Finally, in Sec. VIII we summarize our
findings.

##  II Resonant Jaynes-Cummings models and their Fock state lattices

###  II.1 Multi-mode Jaynes-Cummings models

In this section we introduce the JC model, which, as the simplest example of a
spin-boson model, will serve as a starting point for building other models
with interesting underlying FSLs. The JC model, first introduced by Edwin
Jaynes and Fred Cummings in 1963 JC , provides a fully quantum mechanical
treatment of an atom interacting with an electromagnetic field Themis . The
original model consists of a two-level atom, pseudo spin-1/2 particle,
interacting with a single bosonic mode representing the electromagnetic field.
The interaction can be understood as an exchange of excitations between the
atom and the field; when the ground state is excited to the excited state, the
number of bosons is lowered by one, and vice versa for de-excitation to the
ground state. This gives rise to a Hamiltonian on the form

$\hat{H}_{\mathrm{JC}}=\hat{H}_{0}+\hat{H}_{\mathrm{int}},$ (1)

with the non-interacting, or bare, part of the Hamiltonian $\hat{H}_{0}$
comprises the free energies of the field and atom, i.e.

$\hat{H}_{0}=\omega\hat{n}+\frac{\Omega}{2}\hat{\sigma^{z}},$ (2)

where we have set $\hbar=1$, and $\omega$ is the photon frequency and $\Omega$
the transition frequency for the two-level atom. Further, for the boson number
operator $\hat{n}|n\rangle=n|n\rangle$ for some Fock state $|n\rangle$, and
$\hat{\sigma}^{z}$ is the Pauli $z$-matrix which upon acting on the bare
atomic states; $\hat{\sigma}^{z}|g\rangle=-|g\rangle$ and
$\hat{\sigma}^{z}|e\rangle=|e\rangle$, with $|g\rangle$ and $|e\rangle$ the
lower and upper atomic states respectively. Throughout, following the logic of
bosonic states, we will also call the atomic basis states, like $\ket{g}$ and
$\ket{e}$, as Fock states. Returning to the interaction part
$\hat{H}_{\mathrm{int}}$, both the atom and the bosonic mode emit and absorb
excitations, and we may write the interaction Hamiltonian as Themis

$\hat{H}_{\mathrm{int}}=g\left(\hat{\sigma}^{+}+\hat{\sigma}^{-}\right)\left(\hat{a}+\hat{a}^{\dagger}\right),$
(3)

with $g$ the light-matter interaction strength. It includes four different
contributions; two of which preserve the number of excitations in the system,
and two corresponding to simultaneous excitation/de-excitation of the two
degrees-of-freedom. The latter do not conserve the bare energy of the system,
and for now we shall drop them both. This is what is called the rotating wave
approximation (RWA) and it leads to an interaction on the form Themis

$\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}\hat{\sigma}^{-}+\hat{\sigma}^{+}\hat{a}\right).$
(4)

Moving to the interaction picture with respect to $\omega\hat{N}$, with
$\hat{N}=\hat{n}+\frac{1}{2}\hat{\sigma}^{z}$ the excitation number operator,
the number of free parameters can be reduced from three to two. Defining the
atom-field detuning as $\Delta=\Omega-\omega$, the Hamiltonian (in this
interaction picture) can be written as

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

which is the celebrated JC Hamiltonian.

The JC model can be expanded in numerous ways, of which the most common ones
are Themis : multi-modes where more than one boson mode is considered, multi-
level atom in which more than two electronic levels couple to the
electromagnetic field, multi-atom where $N$ identical two-level atoms couple
to the boson mode (the so called Dicke model), or pumped/driven models.
Focusing on the first extension for now (and return to the other below), and
labelling each boson mode by an index $i$ and assuming that the atom couples
to each bosonic mode with a corresponding coupling strength $g_{i}$, the
multi-mode JC model takes the from

$\hat{H}=\frac{\Omega}{2}\hat{\sigma}^{z}+\sum_{i}\omega_{i}\hat{n}_{i}+\sum_{i}g_{i}\left(\hat{\sigma}^{+}\hat{a}_{i}+\hat{\sigma}^{-}\hat{a}_{i}^{\dagger}\right).$
(6)

Like for the regular JC model we could turn to an interaction picture with
respect to the total excitation number, which is especially convenient in the
degenerate case when $\omega_{i}\equiv\omega,\,\,\,\forall\,\,i$.

###  II.2 Fock state lattices for the resonant multi-mode Jaynes-Cummings
models

From the single mode JC model we first see how the lattice like structure
emerges in state space. As long as the coupling $g\neq 0$, the Fock states are
in general not the eigenstates of the JC model, instead they provide a natural
basis for analyzing the dynamics of the system. Consider first the Fock states
for the single mode JC model, $\ket{n,e}$ and $\ket{n,g}$. Letting the
Hamiltonian in (5) act on these states one finds that they, except for the
vacuum state, will couple pairwise as

$\ket{n,g}\longleftrightarrow\ket{n-1,e}.$ (7)

Thus, the JC Hamiltonian is $2\times 2$ block diagonal in the Fock basis.

It is important to appreciate that the Fock states are coupled by the
interaction Hamiltonian (4) alone, while the bare part $\hat{H}_{0}$ of the
Hamiltonian acts only as onsite energy shifts in the FSL. Throughout this
paper we break up the Hamiltonian of the problem as
$\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{int}}$, where the first part only
contributes to onsite energy shifts but not any kinematics, which instead is
contained within the interacting part $\hat{H}_{\mathrm{int}}$. Indeed, the
lattice geometry and its tunneling rates are defined via
$\hat{H}_{\mathrm{int}}$, and we will thereby focus on the lattices deriving
from the interaction Hamiltonians without paying any attention on the energy
shifts. In some occasions the bare part will only contribute to an overall
energy shift, while in others its role is more delicate and in such cases we
imagine working in an interaction picture.

Returning to the JC model, the Hamiltonian decouples into blocks, and the
related lattice structure becomes rather trivial as shown in Fig. 1 (a). We
find an infinite set of two-site lattices, alternatively a ladder lattice with
the tunneling rates along the legs being zero. One should note that as
$\hat{a}^{\dagger}\ket{n}=\sqrt{n}\ket{n+1}$, the tunneling rate scales with
$\sqrt{n}$ for the two-site lattice $n$. This is general for all our models,
i.e. the tunnelings are not constant throughout the lattices, and hence the
lattices do not possess translational invariance. A consequence is that
Bloch’s theorem does not apply and we typically find a point spectrum, rather
than a band spectrum.

For the two-mode JC model a similar, but still conceptually different, picture
arises. The interaction Hamiltonian in Eq. (6) reduces to

$\hat{H}_{\mathrm{int}}=\left(g_{A}\hat{\sigma}^{+}\hat{a}+g_{B}\hat{\sigma}^{+}\hat{b}+h.c\right).$
(8)

As there are now two boson modes the Fock states will take the form
$\ket{n_{A},n_{B},g}$ and $\ket{n_{A},n_{B},e}$. According to Eq. (8) the
states couple as

$\begin{split}..\ket{n_{A},n_{B}-1,e}\leftrightarrow\ket{n_{A},n_{B},g}\leftrightarrow\ket{n_{A}-1,n_{B},e}..\,.\end{split}$
(9)

Hence, for the two-mode JC model, a lattice like structure of two identical
copies of 1D chains results, see Fig.1 (b). As for the single mode JC model,
the couplings scale with $\sqrt{n}$. In the two-mode JC model the atom can
mediate photon transfer between the two modes Themis , which for the FSL
implies moving left or right along the 1D lattice.

$(\textbf{a})$$g\sqrt{n-1}$$g\sqrt{n}$$|n-1,g\rangle$$|n,g\rangle$$|n-2,e\rangle$$|n-1,e\rangle$

$(\textbf{b})$$|0,n-1,e\rangle$$|n-1,0,e\rangle$$|0,n,g\rangle$$|n,0,g\rangle$$|n-1,1,g\rangle$$|1,n-1,g\rangle$$\cdots$

Figure 1: The FSL for the single-mode JC (a) and two-mode JC (b) models, with
red/blue dots marking the atomic $|e/g\rangle$ states. In the former one finds
an infinite series of two-site lattices, and in the latter one has a 1D chain.
The black lines indicate the non-zero tunnelings of the lattice in the
resonant limit, whereas the orange ones (in the two-mode model) are the same
but in the large detuned limit when the excited atomic level $\ket{e}$ has
been adiabatically eliminated. In both cases the tunneling rates scale such
that the lattice is not translationally invariant ($\sqrt{n_{i}}$ in the
resonant case – black lines, and $\sqrt{n_{a}n_{b}}$ in the large detuning
limit – orange lines).

Finally we consider the FSL for the three-mode JC model, which was first
introduced in Ref. FSL1 . For a single mode there is only a single
neighbouring state to tunnel to (see Eq. (7)), while for two modes there are
two neighbouring states, and we then build up a 1D chain. For three modes,
$a$, $b$, and $c$, we instead have three neighbouring sites, and as such the
resulting lattice will be hexagonal as shown in Fig. 2. The hexagonal lattice
has a sublattice triangular structure seen as the red and blue dots in the
figure, which further represent the different atomic states; $|g\rangle$
(blue) and $|e\rangle$ (red). These triangular sublattices will be those
surviving in the large detuning limit as shown in the next section. The
translationally invariant hexoganl lattice has several interesting features
jonasbook , e.g. the presence of Dirac cones at the corners of the Brillouin
zone. Furthermore, if tunneling beyond nearest neighbour is considered
(Haldane model), the model may host topological features which was analyzed
for the FSL in jctop . You could generalize this construction of the FSLs to
the $n$-mode JC model which will result in a lattice in an $(n-1)$-dimensional
space.

$|0,3,0,g\rangle$$|3,0,0,g\rangle$$|0,0,3,g\rangle$ Figure 2: Same as for Fig.
1 but for the three-mode JC model. Here we show it for the $N=3$ excitation
sector, which results in a finite lattice comprising 16 sites (increasing $N$
makes the lattice larger). In the resonant case, the FSL forms a hexagonal
structure in which each site couples to three neighbours (envisioned by the
solid black lines), while in the detuned limit one of the atomic states (red
or blue dotes) are adiabatically eliminated and a triangular FSL derives,
where every site couples to six neighbouring sites instead (shown by the solid
orange lines).

###  II.3 Beyond the rotating wave equation – the quantum Rabi model

$\sigma^{+}$$\hat{a}\sigma^{+}$ (I)$\hat{a}^{\dagger}\sigma^{+}$ (II)$\hat{a}$
(III)$\hat{a}^{\dagger}$
(III)$|n-1,g\rangle$$|n,g\rangle$$|n+1,g\rangle$$|n-1,e\rangle$$|n,e\rangle$$|n+1,e\rangle$
Figure 3: Schematic representation of the zeroth and first order couplings of
Fock states. The order of the coupling refers to the number of participating
boson creation or annihilation operators in each tunneling process.

Returning to the Hamiltonian of Eq. (3), we are immediately presented with
another lattice model if we keep the counter rotating terms rather than
dismiss them. This model is called the quantum Rabi model Themis ; parity .
More compactly, the interaction Hamiltonian can be expressed as
$\hat{H}_{int}=g(\hat{a}^{\dagger}+\hat{a})\hat{\sigma}_{x}$, where
$\hat{\sigma}_{x}$ is the Pauli-$x$ matrix;
$\hat{\sigma}_{x}|g(e)\rangle=|g(e)\rangle$. This Hamiltonian does not
preserve the number of excitations as the JC Hamiltonian of Eq. (5). The
counter rotating terms break down the continuous symmetry of the JC model to a
discrete $\mathbb{Z}_{2}$-parity symmetry parity . We will discuss the role of
symmetries in more general terms below. For the FSL, these terms change the
number of excitations by one, i.e. they connect the $2\times 2$ blocks of the
JC Hamiltonian. In addition, each Fock state now couples to two neighboring
Fock states;

$\dots\ket{e,-2}\longleftrightarrow\ket{g,n-1}\longleftrightarrow\ket{e,n}\longleftrightarrow\ket{g,n+1}\dots\,,$
(10)

With two neighbours the FSL becomes a 1D chain (see Fig. 4). The
aforementioned $\mathbb{Z}_{2}$-symmetry implies two decoupled parity sectors,
meaning that we actually get two copies of the chain. If we consider the
anisotropic quantum Rabi model Anrabi2 , we have that the counter rotating
terms couple with a different coupling amplitude than the excitation
preserving terms, i.e. the interaction Hamiltonian takes the form

$\hat{H}_{\mathrm{int}}=g_{\mathrm{jc}}\left(\hat{\sigma}_{+}\hat{a}+\hat{a}^{\dagger}\hat{\sigma}_{-}\right)+g_{\mathrm{ajc}}\left(\hat{\sigma}_{-}\hat{a}+\hat{a}^{\dagger}\hat{\sigma}_{+}\right).$
(11)

This does not change the lattice geometry, but an alternation between the
tunneling rates of every second site occurs. For a translational invariant
system this would mean that the lattice had a bi-partite structure where every
unit cell would be comprised of two sites. This is the well known Su-
Schrieffer–Heeger (SSH) model ssh ; jonasbook  which has exponetially
localized topologically protected zero energy edge states. Similar edge states
are also found in the anisotropic quantum Rabi model, but an obstacle here is
that the FSL of the anisotropic quantum Rabi model has only one and not two
edges (the boson Fock space has no upper bound). Furthermore, even though the
edge state is localized to the edge it is not exponentially localized.
However, one can introduce another model supporting a finite FSL, which we do
in Sec. VI, in which the edge states come out more natural and they are
exponentially localized.

The above serves also as a demonstration how the different coupling terms
manifest as tunneling rates in the FSL. This can be done systematically in
order to construct other FLSs. In Fig. 3 the three different couplings of a
Fock state consisting of a single boson mode and a spin half excitation are
shown. We label the couplings $(I)-(III)$, and refer to them as first order
couplings since they all consist of a single boson creation or annihilation
operator, and therefore all scale with $\sqrt{n}$. The last coupling shown in
Fig. 3 only comprises a coupling of the spin mode - and no boson operator -
and is therefore a zeroth order. We do not consider the purely spin-coupling
any further; higher order couplings typically lead to couplings beyond nearest
neighbours, and they may lead to new interesting lattice like structure. The
$(I)$ term is the JC coupling and the $(II)$ the counter rotating coupling,
which we both have already been discussed. The last couplings, marked by
$(III)$, stem from driving the boson mode. The inclusion of the drive term is
discussed for the JC model in the section below and for the Quantum Rabi model
in Sec. II.

###  II.4 Driven Jaynes-Cummings like systems

As mentioned above, the counter rotating terms lowers the symmetry of the JC
model, and the $2\times 2$ block structure of the Hamiltonian in the Fock
basis is lost. Looking at the JC FSL in Fig. 1 (a), this implies additional
coupling terms such that we find two 1D chains; one for each parity.
Naturally, we can break the continuous $U(1)$ number conservation symmetry in
different ways. A natural way to break this symmetry is to consider the driven
JC model, where we include a field drive

$\hat{H}_{\mathrm{fd}}=\eta\left(\hat{a}+\hat{a}^{\dagger}\right),$ (12)

or an atom drive

$\hat{H}_{\mathrm{ad}}=\eta\hat{\sigma}_{x}.$ (13)

For either of these, the excitation number is no longer preserved, while both
still support a $\mathbb{Z}_{2}$-parity symmetry. More precisely, the field
drive causes non-zero horizontal tunneling in the FSL of Fig. 1 (a), and the
atom drive results in diagonal tunneling terms. As for the anisotropic quantum
Rabi model, which mimicked the SSH model with broken translational invariance,
the driven quantum Rabi model realizes a Creutz ladder cl  with broken
translational invariance, see Fig. 4 (d). The Creutz ladder has attracted
attention as it may host dispersionless, i.e. flat, bands. Naturally, such
flat bands do not occur in the driven quantum Rabi model since it is no longer
translationally invariant. However, remnants of these bands are still present.
For field driving, the interaction Hamiltonian can be written as

$\hat{H}_{\mathrm{int}}^{(\pm)}=\left(g\pm\eta\right)\left(\hat{a}+\hat{a}^{\dagger}\right),$
(14)

which results in two energy branches for which one is flat provided
$g=\pm\eta$. This, however, is true in the interaction picture; the bare
Hamiltonian $\hat{H}_{0}$ would lift this degeneracy as well as causing the
Hamiltonian to be bounded from below.

##  III Note on symmetries and lattice dimensions

We have seen a couple of examples in which the structure, and especially the
dimensionallity, of the lattice is changed if the symmetry of the system is
altered. Here we give some general rules for how this comes about.

Returning to Fig. 3 we may also use it to prove a recipe for how one may
construct even more elaborate FSLs. First we notice that in this figure we
consider a system of a single boson mode and a two-level atom, but we could
expand it to multi-level atoms, several boson modes or a combination of these.
If we replace the two-level atom by an $S$-level atom, it will actually
provide a pseudo dimension that will be $2S+1$ sites deep, while including
another boson mode adds an extra full dimension (infinitely deep). Generally,
the bosonic modes, constituting one continuous degree-of-freedom, each gives
rise to one lattice dimension, such that with each additional bosonic mode the
dimension of the FSL grows with one. The atom, which provides a discrete
degree-of-freedom, causes an extra lattice pseudo dimension, i.e. the lattice
is finite in the new dimension.

A continuous $U(1)$ symmetry reduces the lattice dimension by unity, while a
discrete $\mathbb{Z}_{n}$ symmetry decouple the lattice into sublattices. For
example, as the JC model consists of one continuous boson degree-of-freedom
and a discrete atomic degree-of-freedom, one might expect a ladder with two
legs. However, number conservation reduces the lattice dimension by one, which
prohibits tunneling along the legs of the ladder. The quantum Rabi model, on
the other hand, does not support particle conservation and we allow tunneling
along the legs, while the parity symmetry causes the rung tunnelings to be
zero. Note, further, that a continuous symmetry can make the size of the
lattice finite, apart from reducing the dimension. Examples of these are given
in Sec. IV.

Given some Hamiltonian, its FSL typically supports some discrete point
symmetries. Take for example the triangular FSL of Fig. 2, which clearly has a
$2\pi/3$ rotational symmetry around its center site, and three mirror
reflections. These lattice symmetries can be traced back to symmetries of the
Hamiltonian, e.g. the rotational symmetry translates to the unitary
transformation

$\hat{a}\rightarrow\hat{b},\,\hat{b}\rightarrow\hat{c},\quad\text{and}\quad\hat{c}\rightarrow\hat{a}$
(15)

The lattice geometry, e.g. number of neighbours, is determined by the specific
form of the interaction Hamiltonian. When expressed in terms of
raising/lowering operators (e.g. $\hat{a}$ and $\hat{\sigma}^{-}$), the number
of terms in $\hat{H}_{\mathrm{int}}$ gives the maximum number of possible
neighbours. In the quantum Rabi model there are four terms, but only two cause
coupling to other states since
$\hat{\sigma}^{+}\ket{e}=\hat{\sigma}^{-}\ket{g}=0$. The interaction terms may
represent tunneling either between nearest neighbours in the FSL or beyond,
like for the Creutz ladder as shown in Fig. 4 (d).

$(a)$$|n-1,g\rangle$$|n,g\rangle$$|n+1,g\rangle$$|n-2,g\rangle$$|n-2,e\rangle$$|n-1,e\rangle$$|n,e\rangle$$|n-1,e\rangle$

$(b)$$|n-1,g\rangle$$|n,g\rangle$$|n+1,g\rangle$$|n+2,g\rangle$$|n-1,e\rangle$$|n,e\rangle$$|n+1,e\rangle$$|n+2,e\rangle$

$(c)$$|n-1,g\rangle$$|n,g\rangle$$|n+1,g\rangle$$|n+2,g\rangle$$|n-2,e\rangle$$|n-1,e\rangle$$|n,e\rangle$$|n+1,e\rangle$

$(d)$$|n-1,g\rangle$$|n,g\rangle$$|n+1,g\rangle$$|n+2,g\rangle$$|n-1,e\rangle$$|n,e\rangle$$|n+1,e\rangle$$|n+2,e\rangle$

Figure 4: Summary of of FSLs for models consisting of only a single boson and
two-level atomic degree-of-freedom. In (a) the lattice of the JC model (4) is
shown. Here, the conservation of $\hat{N}$ reduces the dimension of the system
such that an infinite set of two-site lattices emerges. The second plot (b)
gives the FSL for the quantum Rabi model (11), with red lines marking the
tunneling due to the counter rotating terms and the black lines the JC
interaction terms. Note that this is actually two copies of 1D chains, even
though they are drawn in this zig-zag way. The FSL of the driven JC model is
shown in (c). The tunneling stemming from the driving are marked by the blue
lines, i.e. by removing them we return to JC lattice (a). Finally in (d) we
display the FSL of the driven quantum Rabi model (again, the lattice in (b) is
regained by omitting the blue lines). This is a type of Creutz lattice with
known interesting properties (see main text). A drive of the spin instead of
the boson would result in vertical blue lines rather than horizontal ones.

##  IV The large detuning limit – quadratic bosnic models

###  IV.1 Effective bosonic models

Quadratic spin-boson models are intrinsically non-linear due to the finite
spin Hilbert spaces. However, there is a natural regime in which the models
become approximately linear, namely in the large detuning limit. Here we
assume $\left|\Delta\right|\gg g\sqrt{n}$. The fast timescale is then set by
$\Delta$, and we face the situation where the atomic degree-of-freedom follows
that of the boson. Hence, it is enough to only consider the evolution of the
bosons. The effective Hamiltonian that dictates the bosonic dynamics can be
derived by adiabatic elimination adel . Here one first writes down the
Heisenberg equations for both the boson and spin operators, which for the JC
model read

$\begin{array}[]{lll}\partial_{t}\hat{a}&=&-i\left[\hat{a},\hat{H}_{\mathrm{JC}}\right]=-ig\hat{\sigma}^{-},\\\
\\\
\partial_{t}\hat{\sigma}^{-}&=&-i\Delta\hat{\sigma}^{-}+ig\hat{a}\hat{\sigma}_{z},\\\
\\\
\partial_{t}\hat{\sigma}_{z}&=&2ig\left(\hat{a}^{\dagger}\hat{\sigma}^{-}\right).\end{array}$
(16)

As the fast timescale is set by $\Delta^{-1}$, the spin operators can be
replaced by their steady state solutions, i.e.
$\hat{\sigma}^{-}=\frac{g}{\Delta}\hat{a}\hat{\sigma}_{z}$, and after
inserting this into the first equation above, the effective dynamical equation
for $\hat{a}$ becomes

$\partial_{t}\hat{a}=-\frac{g^{2}}{\Delta}\hat{a}\hat{\sigma}^{z}\equiv-i\left[\hat{a},\hat{H}_{\mathrm{eff}}\right],$
(17)

where in the second step we define the effective bosonic Hamiltonian
$\hat{H}_{\mathrm{eff}}$. Solving this operator equation gives

$\hat{H}_{\mathrm{eff}}=\frac{2g^{2}}{\Delta}\hat{n}\hat{\sigma}^{z}.$ (18)

In the case of more than one boson modes, the adiabatic elimination is
straightforward to generalize. With the summation over an index $i$ on the
bosonic operator, the effective Hamiltonian stays the same such that

$\begin{array}[]{lll}\partial_{t}\hat{a}_{i}&=&-ig_{i}\hat{\sigma}^{-},\\\ \\\
\partial_{t}\hat{\sigma}^{-}&=&\displaystyle{-i\Delta\hat{\sigma}^{-}+i\sum_{i}g_{i}\hat{a}_{i}\hat{\sigma}_{z}},\\\
\\\
\partial_{t}\hat{\sigma}_{z}&=&\displaystyle{2i\sum_{i}g_{i}\hat{a}_{i}^{\dagger}\hat{\sigma}^{-}}.\end{array}$
(19)

Hence, with $\partial_{t}\hat{\sigma}^{-}=0$ we obtain
$\hat{\sigma}^{-}=\sum_{i}\frac{g_{i}\hat{a}_{i}\hat{\sigma}_{z}}{\Delta}$
such that

$\hat{H}_{\mathrm{eff}}=\sum_{i\neq
j}\frac{2g_{i}g_{j}}{\Delta}\hat{a}^{\dagger}_{i}\hat{a}_{j}\hat{\sigma}^{z}.$
(20)

From a perturbative perspective, the effective Hamiltonians describe the
virtual processes of absorbing and emitting a photon – the Stark shift. In the
multi-mode case, such two-photon processes act as beam-splitters among the
involved modes such that bosons will tunnel from one mode to another. This
manifests as non-vanishing tunneling rates in the corresponding FSLs.

###  IV.2 Fock state lattices for multi-mode JC models in the large detuning
limit

With the effective Hamiltonian in Eq. (20) we set out to study their
properties. Noting that in the detuned limit no population is transferred
between the atomic states, the spin part of the Fock state becomes trivial and
we can leave it out for now. Generally, the interaction Hamiltonian for the
detuned multi-mode case can be written on the form
$\hat{H}_{\mathrm{int}}=\sum_{i\neq
j}\tau_{ij}\hat{a}_{i}^{\dagger}\hat{a}_{j}$. This describes non-interacting
bosons living on a lattice defined by the subscripts $i$ and $j$, e.g. for two
and three modes we get a bosonic “dimer” or “trimer” respectively. We,
however, are interested in their respective FSLs. A priori, the tunneling
strengths $\tau_{ij}$ depend on the system details. For degenerate modes and
appropriate polarizations it is natural to assume that they all have the same
amplitudes, i.e. $\tau_{ij}=\tau e^{i\phi_{ij}}$. The phases $\phi_{ij}$ can,
in principle, be experimentally controlled jctop . We include these phases
since we will see that their presence leads to a plethora of interesting
phenomena. Taking this into account, we reach the final Hamiltonian

$\hat{H}_{\mathrm{int}}=\tau\sum_{i\neq
j}e^{i\phi_{ij}}\hat{a}^{\dagger}_{i}\hat{a}_{j}.$ (21)

Note that hermiticity implies that $\phi_{ij}=-\phi_{ji}$, and since the
tunneling rate only sets the time scale we pick $\tau=1$ throughout. For two
modes (label them with $a$ and $b$) this Hamiltonian reduces to
$\hat{H}_{\mathrm{int}}=\left(e^{i\phi}\hat{a}^{\dagger}\hat{b}+e^{-i\phi}\hat{b}^{\dagger}\hat{a}\right)$
in which the phase can be trivially omitted by a gauge transformation.
However, for three or more modes this will not be true in general gauge . The
phases will result in complex tunneling rates in the FSLs, and according to
the idea of the Peierls substitution these phases mimic a synthetic magnetic
flux penetrating the lattice perpendicularly jonasbook . If a loop in the
lattice is formed, the (gauge invariant) magnetic flux $\phi$ through the loop
equals the sum of the phases of the corresponding tunneling rates. Typically
we are interested in the smallest such loop, a lattice plaquette, and its
corresponding flux, see Fig. 5. Consider, for example, three boson modes $a$,
$b$ and $c$, from Eq. (21) we have the three phases $\phi_{ab}$, $\phi_{bc}$,
and $\phi_{ca}$. Each of these can be changed via a gauge transforation, but
the sum $\phi=\phi_{ab}+\phi_{bc}+\phi_{ca}$ remains constant under such
transformations. The corresponding FSL is a triangular lattice with a
staggered magnetic flux, i.e. alternating $\pm\phi$ fluxes through each
triangular plaquette, see Fig. 5. If we make the choice that for any $i<j$
then $\phi_{ij}=\phi$, then for all higher dimensional multi-mode FSLs the
same pattern will emerge, i.e. any 2D plane/cut of the FSL will display a
triangular lattice with a staggered $\pm\phi$ field. We may note in passing,
that it is possible to also have a homogeneous field in the FSL, but it turns
out that this is only possible for non-quadratic boson models, i.e.
interacting many-body systems.

As the Hamiltonian in Eq. (21) preserves the number of particles, there is a
corresponding $U(1)$-symmetry. This lowers the lattice dimension by one, e.g.
the three-mode model generates a 2D FSL and so on. Generally for any number of
modes, the unit cell (in a strict sense we do not have unit cells since the
lattice is not translationally invariant, but here we think only of the
lattice geometry) has the form of a $d$-simplex, where $d$ indicates the
number of dimensions, e.g in 2D it is a triangle, in 3D it is a tetrahedron
and so on. Furthermore, for a given particle number $N$, the lattice is
finite, with the corners represented by the Fock states with all $N$ bosons
occupying a single mode and the remaining ones in vacuum. For example, for
three modes with $N$ bosons, the number of sites of the FSL becomes
$\mathcal{S}=(N+1)(N+2)/2$. Finally, we note that the introduction of a
synthetic field breaks time-reversal symmetry. However, for $\phi=\pi/2$ the
spectrum becomes symmetric around $E=0$ which is a manifestation of a chiral
symmetry given by the unitary
$\hat{U}_{\mathrm{C}}=\hat{K}\exp\left(i\pi\hat{n}_{c}\right)$, where
$\hat{K}$ stands for complex conjugation, which anti-commutes with the
Hamiltonian.

$\phi$$\phi$$\phi$$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$-$\phi$$\phi$$|N,0,0\rangle$$|0,N,0\rangle$$|0,0,N\rangle$
Figure 5: The triangular FSL deriving from the three-mode boson model (21). A
synthetic magnetic flux $\pm\phi$ penetrates each plaquette in a staggered
manner. The corner states are represented by the Fock states with two modes in
vacuum. Along the straight lines, two diagonal and one horizontal, the boson
number is fixed in one of the modes, e.g. for the horizontal lines the third
$c$-mode keeps a fixed particle number. These observations holds for higher
dimensional lattices as well – in the four-mode model one finds a tetrahedral
lattice.

###  IV.3 Closer look at the three-mode model – a fractal spectrum

In this section we limit ourselves to the three-mode model with a varying
flux. Later we comment also on higher dimensional models. In particular, we
are interested in the energy spectrum, for the Hamiltonian given in Eq. (21),
as a function of the flux $\phi$. The translational invariant case of a 2D
lattice exposed to a constant perpendicular magnetic field results in a
fractal spectrum known as the Hofstadter butterfly hof . It is relevant to
analyze similar situations for the FSLs. For this objective we use that the
Hamitlonian is quadratic and can be written as
$\hat{H}_{\mathrm{int}}=\hat{\bf a}^{\dagger}h\hat{\bf a}$, where

$h=\left[\begin{array}[]{ccc}0&e^{i\phi}&e^{i\phi}\\\
e^{-i\phi}&0&e^{i\phi}\\\ e^{-i\phi}&e^{-i\phi}&0\end{array}\right]$ (22)

and $\hat{\bf
a}^{\dagger}=\left(\hat{a}^{\dagger},\hat{b}^{\dagger},\hat{c}^{\dagger}\right)$.
If we diagonalize $h$, i.e. $d=UhU^{-1}$, the full Hamiltonian is diagonalized
as
$\hat{H}_{\mathrm{int}}=\gamma_{0}\hat{\alpha}_{0}^{\dagger}\hat{\alpha}_{0}+\gamma_{1}\hat{\alpha}_{1}^{\dagger}\hat{\alpha}_{1}+\gamma_{2}\hat{\alpha}_{2}^{\dagger}\hat{\alpha}_{2}$,
with $\hat{\alpha}_{k}$ ($k=0,\,1,\,2$) the transformed boson operators, and
$\gamma_{k}$ the eigenvalues of $h$. Thus, the spectrum for a single boson is
simply the eigenvalues of $h$, while the full spectrum, given $N$ bosons, can
be found combinatorically from these single boson energies $\gamma_{k}$. We
may note that the structure of $h$ in Eq. (22) remains for higher dimensions.
Numerically its eigenvalues can be found for large mode numbers, but
analytical solutions are only found for the three- and four-mode cases. For
the three-mode model, finding the single boson energies reduces to solving the
characteristic polynomial $\gamma^{3}-3\gamma-2\cos\phi=0$, which leads to the
three energies

$\gamma_{k}=2\cos\left(\frac{\phi-2\cdot\pi\cdot k}{3}\right),\qquad
k\in\left\\{0,1,2\right\\}$ (23)

Hence, the full spectrum for a general boson number is written as

$E_{\bf n}=n_{0}\gamma_{0}+n_{1}\gamma_{1}+n_{2}\gamma_{2},$ (24)

where ${\bf n}=(n_{0},n_{1},n_{2})$ subject to the constraint
$n_{0}+n_{1}+n_{2}=N$. To gain further insight of the spectral structure we
consider the energy differences between any two energy levels

$\delta E_{\bf n}=E_{\bf n}-E_{\tilde{\bf
n}}=\left(\gamma_{1}-\gamma_{0}\right)\\!\left(n_{1}-\tilde{n}_{1}\right)+\left(\gamma_{2}-\gamma_{0}\right)\\!\left(n_{2}-\tilde{n}_{2}\right)$
(25)

where we have used that we may rewrite $E_{\bf
n}=N\gamma_{0}+(\gamma_{1}-\gamma_{0})n_{1}+(\gamma_{2}-\gamma_{0})n_{2},$. In
three dimensions there exist certain fluxes $\phi$ for which the spectrum
becomes equidistant or quasi equidistant, i.e. there then exists a few
energies characterizing the involved time-scales. This is opposed to some
general $\phi$ which reproduces a whole range of different time-scales.

Figure 6: Upper two plots, (a) and (b), show the spectra and the distances
between nearby energies as a function of the flux $\phi$, for the three-mode
model with $N=25$ bosons, while the lower two plots, (c) and (d), give the
same for the four-mode case with $N=16$. The fractal structure of the three-
mode case is evident, while the four-mode spectrum does not display the same
fractal properties.

From the expression (25) we see that the energy spectrum becomes quasi
equidistant whenever

$\frac{\gamma_{0}-\gamma_{1}}{\gamma_{1}-\gamma_{2}}=q,\quad q\in\mathbb{Q}.$
(26)

Writing $q=\frac{j}{j^{\prime}}$ with $j,j^{\prime}\in\mathbb{N}$, we find
that given the solutions for $\gamma_{k}$ in (23), that Eq. (26) can be
written as

$j^{\prime}\\!\left[\\!\sqrt{3}\cos\\!\left(\frac{\phi}{3}\right)\\!-\\!\sin\\!\left(\frac{\phi}{3}\right)\\!\right]=j\\!\left[\\!\sqrt{3}\cos\\!\left(\frac{\phi}{3}\right)\\!+\\!\sin\\!\left(\frac{\phi}{3}\right)\\!\right]\\!.$
(27)

As all three single boson energies are continuous functions of $\phi$, so are
$(\gamma_{0}-\gamma_{1})$ and $(\gamma_{1}-\gamma_{2})$, and furthermore,
except for the singularity when $\gamma_{1}=\gamma_{2}$, which happens at
$\phi=0$, also the fraction (26) is a continuous function of $\phi$.
Disregarding this singularity, as $N\rightarrow\infty$ there will be
infinitely many $\phi_{j}$’s that fulfill the condition above, see Fig. 6 (a)
and (b). The spectrum at these “quasi equidistant points” becomes highly
degenerate, but it is only perfectly equidistant for one of them,
$\phi_{1}=\pi/2$. This is the reason why for the other of these points we call
the spectrum quasi equidistant, since here there exists more than a single
energy difference. For a given $\phi_{j}$, the energy differences between the
nearby energies can be ordered $\Delta_{j}$, $2\Delta_{j}$, $3\Delta_{j}$, …,
$j\Delta_{j}$. Note that the basic energy difference $\Delta_{j}$ is different
for the various equidistant points. Solving (26) for $\phi_{j}$ leads to the
following condition

$\tan\left(\frac{\phi_{j}}{3}\right)=\frac{\sqrt{3}}{2j+1},\quad\quad
j\in\mathbb{N}.$ (28)

which we have reached using (27), and introducing a new pair
$n,n^{\prime}\in\mathbb{N}$ where $j+j^{\prime}=n^{\prime}$ and
$j^{\prime}-j=n$. As this holds for any $q$ we may choose $n=1$, which implies
that $n^{\prime}=2j+1$ and leads to the expression above.

The equidistant points given by Eq. (28) occur denser and denser the larger
values $j$. This manifests as a fractal structure of the spectrum, see Fig. 6
showing both the spectrum (a) and its nearby energy differences $\delta E_{\bf
n}$ (b). The results of the energy differences demonstrate the concept of the
equidistant points. As we will demonstrate in the next section, picking $\phi$
according to one of these points will result in novel dynamics.

###  IV.4 The four-mode and higher-mode models

To fractal structure of the spectrum is rooted in the existence of special
fluxes fulfilling the condition (26). The $3\times 3$ matrix $h$ gives three
energies, and hence two energy difference, and the single condition deriving
from comparing these two. For the four-mode model (21) we find four
eigenvalues $\gamma_{k}$ (with $k=1,2,3,4$) of the corresponding matrix $h$.
Following the same arguing that led to the condition (26), for the four-mode
model we are left we three conditions

$\frac{\gamma_{3}-\gamma_{0}}{\gamma_{2}-\gamma_{0}}=q_{1},\qquad\frac{\gamma_{3}-\gamma_{0}}{\gamma_{1}-\gamma_{0}}=q_{2},\qquad\frac{\gamma_{2}-\gamma_{0}}{\gamma_{1}-\gamma_{0}}=q_{3},$
(29)

where $q_{1}$, $q_{2}$, and $q_{3}$ are all rational numbers. We can
generalize this to a $K$-mode model in which the number of conditions becomes
$(K-1)(K-2)/2$. It should be clear that it becomes harder and harder, the
larger $K$ is, to find fluxes $\phi_{j}$ fulfilling all conditions
simultaneously. As long as the conditions cannot be met for single fluxes the
spectrum cannot be fractal as for the three-mode model. This we demonstrate by
displaying the spectrum and energy difference $\delta E_{\bf n}$ for the four-
mode model in Fig. 26 (c) and (d). There is, however, possible to find fractal
spectra also for higher mode models (four and six modes) as discussed next.

####  IV.4.1 Ring models

The models defined via Eq. (21) are “fully connected”, i.e. each mode can
exchange excitations with every other mode. We may lift this by introducing
selections among the tunneling processes, and in particular we consider
bosonic chains with periodic boundary conditions

$\hat{H}=\sum_{i=1}^{K}\hat{a}^{\dagger}_{i}\hat{a}_{i+1}e^{i\phi/K}+\mathrm{h.c},$
(30)

where $\hat{a}_{K+1}=\hat{a}_{1}$. Since we impose periodic boundary
conditions, and thereby the lattice has a ring structure, the phase $\phi$
becomes non-trivial in the sense that we cannot gauge it away. In the boson
chain it mimics a magnetic flux through the lattice ring. For the fully
connected four mode case, the FSL takes the form of a tetrahedron lattice.
When we instead consider the four mode ring model, some of the links between
sites in the tetrahedron lattice are cut; the number of sites each site
couples to goes from 12 to eight. The resulting FSL has a rhombic structure,
i.e. any 2D lattice plane forms a rhombic lattice. The characteristic equation
for the corresponding matrix $h$ becomes
$\gamma^{4}-4\gamma^{2}+4\sin^{2}\left(\frac{\phi}{2}\right)=0$, and the four
eigenvalues (roots) can be expressed as

$\gamma_{k\pm}=\pm\cos\left(\frac{\phi+\pi\cdot k}{2}\right),\qquad
k\in\left\\{0,1\right\\}.$ (31)

This leads to the following condition for the flux $\phi=2\arctan(q)$ with
$q\in\mathcal{Q}$. Similarly, for the six-mode ring model we find the
eigenvalues

$\gamma_{k}=\pm 2\cos\left(\frac{\phi+2\pi\cdot k}{6}\right)\qquad
k\in\left\\{0,1,2\right\\},$ (32)

giving the following condition for the quasi equidistant points
$\phi=6\arctan\left(\frac{q}{\sqrt{3}}\right)$. In Fig. 7 we present the
spectra and energy differences for the two ring models. Clearly, the two
models’ spectra display a similar fractal structure as for the three-mode
model of Figs. 6 (a) and (b). Interestingly, the three-mode model is also a
ring model just like those of Fig. 7. However, we have numerically verified
that the five- and eight-mode models do not reproduce fractal spectra, i.e. it
is not a property of the ring-type coupling. We have noticed, though, that the
spectrum $E_{\bf n}(\phi)$ is symmetric around $E=0$ for all fluxes $\phi$
whenever one considers an even number of boson modes. This suggest an
additional chiral symmetry that holds only for ring models with an even number
of modes.

Figure 7: The same as Fig. 6, but for the four, upper two plots (a) and (b),
and six mode, lower two plots (c) and (d), ring models. The number of bosons
are $N=20$ and $N=10$ for the two cases respectively. By disregarding some of
the terms in the Hamiltonian (21) to obtain the ring Hamiltonian (30), the
number of conditions for having a quasi equidistant spectrum is reduced and it
is possible to fulfill them simultaneously, and thereby the reappearance of
fractal spectra for these two models.

##  V Evolution of initially localized states in the three-mode boson model

To further study manifestations of the fractal spectrum, we turn to study the
evolution over time. While our focus is on the three-mode case, let us start
by a remark on the two-mode model, which produces a 1D FSL with tunneling
rates $t\sqrt{(n+1)(N-n)}$. The Hamiltonian is easily diagonalized by a
Bogoliubov transformation, but we use instead the Schwinger spin-boson mapping
in which we identify a spin operator
$\hat{S}_{x}=\frac{1}{2}\left(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a}\right)$,
such that the Hamiltonian is simply $\hat{H}_{\mathrm{int}}=2\hat{S}_{x}$,
with the total spin $S=N/2$. The spectrum is clearly equidistant, and any
state returns to its initial state after the revival time $T_{R}=2\pi$. This
model was studied by Christandl et al. in terms of perfect state transfer in
1D chains with spatially varying tunneling rates cdel , e.g. after half the
revival time $T_{R}/2$ any initial state at one end of the chain will have
traversed the chain to populate the other end. Recently the model was
experimentally explored with the particular Hamiltonian realized by coupling
two photonic modes via a beam splitter FSsimulation . To see the effect of the
flux $\phi$ we need, however, to go up by one dimension.

As a demonstrating example we consider a system of $N=15$ bosons. Hence, the
Fock states have a form $\ket{n_{a},n_{b},n_{c}}$ where $n_{i}\in\mathbb{N}$
and $n_{a}+n_{b}+n_{c}=15$. $N=15$ results in $\mathcal{S}=136$ different Fock
states that make up the corresponding FSL. We shall look at the evolution of
the two different initial states, taken to be $\ket{15,0,0}$ and
$\ket{11,4,0}$, which both are initially completely localized on a single site
located somewhere in the FSL. In order to characterize the extent of the
states we consider the inverse partition ratio

$IPR(t)=\sum_{\bf n}|\langle{\bf n}|\psi(t)\rangle|^{4},$ (33)

which gives the degree of localization in the FSL. For a fully localized state
$IPR=1$, while for a maximally delocalized one $IPR=1/\mathcal{S}$, where
$\mathcal{S}$ is, as before, the number of sites of the FSL. When it comes to
the choice of fluxes we shall consider the time-evolution for the fluxes
$\phi_{1}$ and $\phi_{2}$ according to Eq. (28). In addition we pick a phase
not coinciding with any equidistant point, namely $\tilde{\phi}=\pi/100$.

Figure 8: (Color online) The $IPR$ for two different initial states, and two
different fluxes corresponding to equidistant points. For both fluxes the
states display perfect revivals, even though the revival times differ between
the fluxes (but not between the states!). As a comparison, we also give the
minimum possible $IPR$ (yellow line). How the $IPR$ is defined it is easy to
get fooled that a small value necessarily implies an extended distribution. In
fact, the distributions remain fairly localized throughout, see Fig. 10.

The numerically extracted $IPR$s are depicted in Fig. 8. As long as the flux
is taken to be one of the $\phi_{j}$ values we see perfect recurrences of the
$IPR$ to its initial value $IPR=1$. The state, however, has not necessarily
returned to its initial state at these instances, but it might have completely
localized on another lattice site as will be clear below. In between the
strong localization the state is spread over several lattice sites, but the
$IPR$ never drops close to its minimum value (a small $IPR$ can still
represent a localized distribution since it scales as $1/$occupied sites).
Hence, even in-between the maximum localization the distribution remains
fairly localized in the lattice, see further Fig. 10. For $\tilde{\phi}$ (not
shown in the figure), when the spectrum is not quasi equidistant, we see no
clear indicators of localization, even though the state is far from populating
the whole lattice.

Figure 9: (Color online) Examples of the evolution of the occupation vector
$\vec{n}(t)=\left(n_{a}(t),n_{b}(t),n_{c}(t)\right)$, with
$n_{i}(t)=\langle\hat{n}_{i}\rangle$ in the FSL for the three-mode model. The
initial states are $\ket{0,0,15}$ (blue curves), and $\ket{0,4,11}$ (red
curves), and the fluxes for the different curves are given in the inset. Note
how the trajectories close whenever the flux is chosen as one of the
equidistant points (signaling perfect revivals), while this is no longer true
for a general flux (thin curves). In addition, the trajectories are found to
be ‘shape invariant’; by which we mean that the two different initial states
result in the same type of trajectories only rotated and scaled in size.

Recall that the constraint stemming from particle conservation removes one
degree-of-freedom, which ensures, for example, that the participation vector
$\vec{n}(\phi,t)=\left(n_{a}(\phi,t),n_{b}(\phi,t),n_{c}(\phi,t)\right)$, with
$n_{i}(\phi,t)=\langle\hat{n}_{i}\rangle$, lives in a 2D plane, i.e. the FSL.
Thus, the trajectory gives us an idea of how the full distributions move
around in the lattice, even though it only captures the mean occupation but no
higher moments of the distribution. In Fig. 9 we show examples of such
trajectories for the same initial states and parameters as those used for Fig.
8. The first to notice is how the trajectory $\vec{n}(\phi_{1},t)$ sets an
outer boundary for all the other trajectories.

As pointed out, when $\phi_{j}$ agrees with one of the equidistant points we
encounter perfect revivals at times

$T_{j}=\frac{2\pi}{\Delta_{j}}$ (34)

where $\Delta_{j}$ is the aforementioned smallest energy splitting for a given
$\phi_{j}$. The revivals manifest as closed trajectories of $\vec{n}(t)$ in
the FSL. These loops form a symmetric pattern in the FSL, and the number of
full 360∘ turns it makes before closing upon itself is found to be $2j-1$.
Thus, for $j=1$ (i.e. $\phi_{1}=\pi/2$) the distribution makes a single loop
around the lattice, which is found to follow closely the lattice edges, i.e.
there always remains a large population imbalance among the three modes. For
$\phi_{2}$, the distribution makes instead three loops before returning tho
its initial states. For a flux away from any $\phi_{j}$ we do not encounter
such revivals and hence no closing of the trajectories. For the closed loops
there are ‘clustering points’ in which many trajectories cross; the lattice
origin being the most noticeable. A closer look reveals that for any initial
Fock state, the resulting trajectories pass the origin (i.e.
$n_{a}(t)=n_{b}(t)=n_{c}(t)=N/3$) for all $\phi_{j}$ with
$j=2,3,5,6,8,9,\dots$. Upon comparing the trajectories originating from the
two different initial Fock states, $\ket{0,0,15}$ vs. $\ket{0,4,11}$, we find
that the shape and structure remain intact, however rotated and contracted. In
the App. A we derive semi-classical equations-of-motion that exactly reproduce
the trajectories.

The trajectories of Fig. 9 give the mean positions of the Fock state
distribution

$P({\bf n},t)=|\langle{\bf n}|\psi(t)\rangle|^{2},$ (35)

but they say nothing about its actual shape. This is shown in Fig. 10 by
plotting the distribution $P({\bf n},t)$ at different times for the case of
$\phi_{1}=\pi/2$. Clearly, the distribution remains fairly localized along the
FSL edges throughout. Furthermore, it propagates counter-clockwise around the
lattice which reflects the breaking of time-inversion as soon as the flux is
non-zero (a zero flux does not produce a favoured direction for the
propagating distribution). The propagation is reminiscent of a Hall current; a
perpendicular magnetic field applied to a 2D electronic gas implies a current
living on the edge of the system. In our case, the magnetic field is
staggered, but we notice that the direction of the current can be reversed by
flipping the sign of the fluxes.

Figure 10: (Color online) The upper light blue triangles display snapshots of
the Fock state distributions (35) for the initial state $A=\ket{0,0,15}$ and
for the flux $\phi_{1}=\pi/2$. Initially the distribution occupies a single
corner site, and after $T_{R}/3$ it has been transferred to another corner
site. In between it stays well localized and traverses along the lower edge of
the lattice. The lower plot gives the corresponding $IPR$ (33), with the
vertical lines indicating the time instances used for the distribution
snapshots. Note that for the center distribution at $T_{1}/6$ we have $IPR\ll
1$, but still the distribution is well localized according to the upper
figure. Figure 11: (Color online) The FSL (a) emerging from the central spin
model (37) with anisotropic coupling amplitudes $g_{\pm}$. The red/greed dots
mark Fock states with the spin-1/2 in the $|g\rangle/|e\rangle$-states
respectively. The parity symmetry of the model implies a similar FSL but
belonging to the other parity sector with the green/red dots swapped. In (b)
we show the exponentially localized $E=0$ energy eigenestate
$|\psi_{0}\rangle$ in an 101 site lattice with $g_{+}/g_{-}=0.9$. In the other
parity sector there is an equivalent edge state localized to the right edge
instead.

##  VI Towards the SSH model - and the appearance of edge states

It is intriguing to speculate whether there is an underlying topological
structure that causes the edge currents seen in Fig. 10. A more thorough
discussion regarding this topic for the FSLs can be found in Refs. jctop ;
FSsimulation . Since translational invariance is absent we encounter discrete
energies rather than bands, but it is interesting to note that for the chiral
case with $\phi_{1}=\pi/2$ we find $Z=N/2+1$ or $Z=N/2+1/2$ degenerate $E=0$
energies for even and odd particle numbers $N$ respectively. These are not
unique for the FSL, but appear also for the translationally invariant lattice
with a staggered field. In other words, the degenerate $E=0$ energies in the
FSL are the counterpart of a flat band.

Instead of exploring properties of these $E=0$ states for the three-mode
model, let us introduce another model also supporting $E=0$ eigenstates which
are found to be living on the edge of the FSL. In order to construct a finite
1D FSL we do not a boson mode, but rather a spin-$S$ subsystem via the mapping

$\begin{array}[]{ccc}\left(\hat{a}^{\dagger}+\hat{a}\right)&\rightarrow&\hat{S}_{x},\\\
i\left(\hat{a}^{\dagger}-\hat{a}\right)&\rightarrow&\hat{S}_{y},\\\
\hat{a}^{\dagger}\hat{a}&\rightarrow&\hat{S}_{z},\end{array}$ (36)

where the $\hat{S}^{\alpha}$’s are the SU(2) spin operators for a spin-S
particle. Contrary to the case of a boson mode, this introduces both an upper
and a lower bound/edge for the lattice. When this is applied to the JC model
one obtains a so called central spin model - which is a model where a central
spin-$1/2$ particle interacts identically with $N$ non-interacting spin-$1/2$
particles which form the large spin-$S$ (i.e. $S=N/2$). The central spin model
has served as a toy model for studying quantum criticality, as well as for
analyzing non-Markovian decay of a qubit centralspin . To achieve the desired
FSL we consider the central spin model on the form

$\begin{array}[]{lll}\hat{H}_{int}^{cs}&=&g_{x}\hat{\sigma}^{x}\hat{S}^{x}+g_{x}\hat{\sigma}^{y}\hat{S}^{y}\\\
\\\
&=&g_{+}\left(\hat{\sigma}^{+}\hat{S}^{-}+\hat{\sigma}^{-}\hat{S}^{+}\right)+g_{-}\left(\hat{\sigma}^{+}\hat{S}^{+}+\hat{\sigma}^{-}\hat{S}^{-}\right),\end{array}$
(37)

where $g_{\pm}=g_{x}\pm g_{y}$. This is the spin analogue of the anisotropic
quantum Rabi model (11), i.e. the two types of interaction terms couple with
different strengths; for $g_{x}=g_{y}$ we regain a JC type interaction, and
for $g_{y}=0$ we instead get a quantum Rabi type interaction. Like for the
quantum Rabi model, the central spin model supports a $\mathbb{Z}_{2}$-parity
symmetry for general couplings $g_{x}$ and $g_{y}$. This results in two
decoupled 1D chains, one for each parity. Furthermore, the chains have lengths
$2S+1$, and provided $g_{+}\neq g_{-}$ and $g_{+},\,g_{-}\neq 0$, the
tunneling rates will alternate between neighbouring sites (see Fig. 11 (a)),
similar to that of the SSH model. As for the SSH model, also for the central
spin model one finds two $E=0$ eigenstates that become exponentially localized
at either of the edges. More precisely, the probability to populate site $i$
for such a state scales as
$P_{i}\propto\left(\frac{g_{+}}{g_{-}}\right)^{S-i}$. Each localized state
belongs to either of the two parity sectors, and in Fig. 11 (b) we show one of
them for $g_{+}=0.9g_{-}$. For the two-mode JC model on resonance, one also
find localized edge states jctop , but these, however, are not exponentially
localized as the ones of the central spin model.

Model | Interaction Hamiltonian | Lattice type
---|---|---
Jaynes-Cummings | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}\hat{\sigma}_{-}+\hat{a}\hat{\sigma}_{+}\right)$ | Double-well, Fig. 4 (a)
Quantum Rabi | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}+\hat{a}\right)\hat{\sigma}_{x}$ | 1D chain, Fig. 4 (b)
Anisotropic quantum Rabi | $\hat{H}_{\mathrm{int}}=g_{\mathrm{jc}}\left(\hat{\sigma}_{+}\hat{a}+\hat{a}^{\dagger}\hat{\sigma}_{-}\right)+g_{\mathrm{ajc}}\left(\hat{\sigma}_{-}\hat{a}+\hat{a}^{\dagger}\hat{\sigma}_{+}\right)$ | Infinite SSH chain
$N$ atom Dicke | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}+\hat{a}\right)\hat{S}_{x}$ | Square lattice
$N$ atom Tavis-Cummings | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}\hat{S}^{-}+\hat{a}\hat{S}^{+}\right)$ |  $N$ potential-well
Central spin model | $\hat{H}_{\mathrm{int}}=g_{x}\hat{\sigma}_{x}\hat{S}_{x}+g_{y}\hat{\sigma}_{y}\hat{S}_{y}$ | Finite SSH chain, Fig. 11 (a)
Driven quantum Rabi | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}+\hat{a}\right)\hat{\sigma}_{x}+\eta\left(\hat{a}^{\dagger}+\hat{a}\right)$ | Creutz ladder, Fig. 4 (d)
Single mode $\Lambda$ | $\hat{H}_{\mathrm{int}}=g\left(\hat{a}^{\dagger}+\hat{a}\right)\left(\hat{\lambda}^{(1)}+\hat{\lambda}^{(6)}\right)$ | Lieb ladder
Two-mode JC | $\hat{H}_{\mathrm{int}}=g_{a}\left(\hat{a}^{\dagger}\hat{\sigma}_{-}+\hat{a}\hat{\sigma}_{+}\right)+g_{b}\left(\hat{b}^{\dagger}\hat{\sigma}_{-}+\hat{b}\hat{\sigma}_{+}\right)$ | SSH chain
Two-mode detuned JC | $\hat{H}_{\mathrm{int}}=t\left(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a}\right)$ | CDEL chain
Two-mode quantum Rabi | $\hat{H}_{\mathrm{int}}=g\left[\left(\hat{a}^{\dagger}+\hat{a}\right)+\left(\hat{b}^{\dagger}+\hat{b}\right)\right]\hat{\sigma}_{x}$ | (Layered) square lattice
Two mode $\Lambda$ | $\hat{H}_{\mathrm{int}}=g_{a}\left(\hat{a}^{\dagger}+\hat{a}\right)\hat{\sigma}_{12}+g_{b}\left(\hat{b}^{\dagger}+\hat{b}\right)\hat{\sigma}_{23}$ | 2D Lieb lattice
Three-mode JC | $\hat{H}_{\mathrm{int}}=g_{a}\left(\hat{a}^{\dagger}\hat{\sigma}_{-}+\hat{a}\hat{\sigma}_{+}\right)+g_{b}\left(\hat{b}^{\dagger}\hat{\sigma}_{-}+\hat{b}\hat{\sigma}_{+}\right)+g_{c}\left(\hat{c}^{\dagger}\hat{\sigma}_{-}+\hat{c}\hat{\sigma}_{+}\right)$ | Hexagonal lattice
Three-mode detuned JC | $\hat{H}_{\mathrm{int}}=t\left(\hat{a}^{\dagger}\hat{b}\,e^{i\varphi}+\hat{b}^{\dagger}\hat{c}+\hat{a}^{\dagger}\hat{c}+\text{h.c.}\right)$ | Triangular lattice, Fig. 5
Three-mode tripod | $\hat{H}_{\mathrm{int}}=g_{a}\left(\hat{a}^{\dagger}+\hat{a}\right)\hat{\sigma}_{12}+g_{b}\left(\hat{b}^{\dagger}+\hat{b}\right)\hat{\sigma}_{13}+g_{c}\left(\hat{c}^{\dagger}+\hat{c}\right)\hat{\sigma}_{14}$ | Perovskite lattice
Three-mode quantum Rabi | $\hat{H}_{\mathrm{int}}=g\left[\left(\hat{a}^{\dagger}+\hat{a}\right)+\left(\hat{b}^{\dagger}+\hat{b}\right)+\left(\hat{c}^{\dagger}+\hat{c}\right)\right]\hat{\sigma}_{x}$ | Cubic lattice
Four-mode detuned JC1 | $\hat{H}_{\mathrm{int}}=t\sum_{j,i=1}^{4}e^{i\phi_{ij}}\hat{a}^{\dagger}_{i}\hat{a}_{j}$ | Tetrahedral lattice
Four-mode detuned JC2 | $\hat{H}_{\mathrm{int}}=t\sum_{i=1}^{4}e^{i\phi_{i}}\hat{a}^{\dagger}_{i}\hat{a}_{i+1}+\text{h.c}.$ | Rhombic 3D lattice
Table 1: A collection of atom-light interaction Hamiltonians and their
corresponding FSLs. Many of these appear frequently in the quantum optics
community. Some of them we have discussed in some detail in the main text,
e.g. the two-mode and three-mode JC models.
$\hat{\sigma}_{ij}=|i\rangle\langle j|+|j\rangle\langle i|$.

##  VII A zoo of Fock state lattices

In quantum optics one encounters a range of models, comprising few degree-of-
freedoms, that describe the interaction between quantized light and matter.
Many of these are extensions and generalizations of the JC model towards
multi-mode and multi-level atoms Themis . They often generate interesting FSLs
with known, translationally invariant, counterparts in condensed matter
physics. In Tab. 1 some of these models are listed, with their respective
interaction Hamiltonians, and with information about the type of FSL that
represents them. Below we discuss a few of these models

  * •

Driven quantum Rabi model. The FSL belonging to the quantum Rabi model was
shown in Fig. 4 (b), where the $\mathbb{Z}_{2}$-parity symmetry splits the
lattice into two 1D chains, one for each parity. We may break this parity by
including a driving term, which induces tunneling along the legs of the FSL
ladder, as depicted in Fig. 4 (d). This is a special type of a Creutz ladder
Creutz . The more standard Creutz ladder appears if we also include an atomic
driving term that results in non-zero tunneling terms along the rungs of the
ladder. In such a ladder we have competing nearest and next nearest
neighbouring tunneling terms. A property of the Creutz ladder lattice is the
presence of a dispersionless flat band when the tunneling rates are properly
tuned. For our model, if we rotate the spin,
$\hat{\sigma}_{x}\rightarrow\hat{\sigma}_{z}$, the Hamiltonian takes the
simple form $\hat{H}_{\mathrm{int}}=(\eta\pm
g)\left(\hat{a}^{\dagger}+\hat{a}\right)$. For $g=\eta$, i.e. the same
amplitude of the tunneling rates, we find an infinite number of degenerate
$E=0$ eigenstates. These are the corresponding flat band of the Creutz ladder,
now manifested in the FSL.

  * •

Two-mode $\Lambda$ model. The Lieb lattice is another paradigm model hosting a
flat band. In 2D it is a “punctured” squared lattice where every fourth
lattice point has been excluded. This implies that the unit cell now contains
three sites, with two of them having non-zero tunneling elements to two
neighbours, and the third site couples to four neighbours. The lattice
geometry provides a destructive interference mechanism that may completely
hinder mobility, resulting in a $E=0$ flat band in the tight-binding limit
Lieb . The remaining two bands touch in a Dirac cone in the center of the
Brillouin zone. We can achieve this lattice type in a system of a three-level
$\Lambda$ atom (two lower states and one excited), with the transition of each
arm of the atom coupled to different boson modes. One finds that the
destructive interference survives also in the FSL, and one thereby derives an
infinite number of degenerate $E=0$ eigenstates.

  * •

Three-mode tripod model. The generalization of the three-level $\Lambda$ atom
to four levels consists in adding one more arm, which is called the tripod
atom. We can then couple each arm of the tripod to an individual mode such
that we have a three-mode model. The resulting 3D lattice, bares similarities
with the Lieb lattice, i.e. it is a punctured cubic lattice with a unit cell
containing four points. The translational invaraint version is called a
Perovskite type lattice, and as for the Lieb lattice it supports bands with
non-trivial topology, i.e. their Chern numbers are non-zero perov . The flat
bands in both the two-mode $\Lambda$ and three-mode tripod systems are related
to the existence of dark states in these models, i.e. $E=0$ eigenstates.

The list of Tab. 1 is by no means complete. As shown in Ref. jctop , many
topological features of translationally invariant lattices survive in FSLs
with the same geometry. Topological invariants, like the Chern number, are not
directly applicable for these FSLs, and instead alternative topological
numbers are studied jctop . Hence, given your favourite lattice with say some
interesting properties, e.g. non-trivial topology, you can work backwards and
ask for some interaction Hamiltonian that generates the corresponding FSL. For
example, bi-layered lattices may host novel properties, like bi-layered Lieb
lattices can produce high Chern numbers, and it is often straightforward to
introduce more layers in the FSL by including additional internal atomic
levels.

##  VIII Concluding remarks

In this work we have taken a new approach in order to analyze a range of
different spin-boson models known from quantum optics. This idea relying on
FSLs was first brought up in FSL1 , and later extended in jctop . In these
works, focus was on the multi-mode JC model. Part of the present work explores
the same models, and we find a set of new properties. The spectra for the
three-mode boson model, when a synthetic magnetic flux has been added, is
found to be fractal, akin to the Hofstadter butterfly. The flux of the FSL is,
however, staggered rather than uniform (a uniform flux would not render a
fractal spectrum for the triangular FSL). The fractal structure cannot be
explained by the same arguments used for the Hofstadter butterfly, as this
relies on varying sizes of the lattice’s unit cell. Instead, the emergence of
our fractal structure is understood in terms of certain fluxes $\phi_{j}$ that
causes a high degeneracy in the spectrum. For higher number of boson modes,
the fractal structure is in general lost, but it can be restored for the four-
and six-mode cases provided one adds some selection among the tunneling terms
between the modes.

The evolution of the system very much depends on the flux – at the quasi
equidistant points $\phi_{j}$ we have that the system shows perfect revivals.
In the FSL, given that we start in a Fock states, for these fluxes the
distribution $P({\bf n},t)$ remains well localized and follows closed loops.
Since the flux breaks time-reversal symmetry, the loops are traversed either
clockwise or anti-clockwise depending on the sign of $\phi_{j}$. One also
finds that these closed trajectories have a mirror symmetry with respect to
the axis intersecting the lattice origin and the site where the state is
initialized. Furthermore, the index $j$ determines how many twists the
distribution makes before it returns to the initial site. There is also an
inherent scale invariance of these trajectories; for a given $\phi_{j}$ the
trajectories are identical in structure irrespective of the specific initial
Fock state, i.e. they only differ in orientation and size.

We also demonstrated how fractal spectra appear in a four- and six-mode boson
model. The properties of the corresponding trajectories, i.e. how the state
evolves within the FSLs, was not consider, but we expect similar behaviour as
for the three-mode model. The FSL of the three-mode model is two dimensional,
the occupation vector $\vec{n}(t)$ evolves in a plane defined by a fixed total
boson number $N$. Likewise, for the four-mode model we have a four component
occupation vector $\vec{n}(t)$, but the conserved particle number implies that
it is projected down to a 3D FSL.

A boson degree-of-freedom results in infinite FSLs, unless some symmetry
restricts it, e.g. particle conservation. Spin degrees-of-freedom have finite
size Hilbert spaces and thereby finite size FSLs. We considered one such
example, the anisotropic central spin model, and showed how this model has a
topological property manifested in $E=0$ states exponentially localized to the
edges of the lattice. Without going into any details we listed other models
and mentioned which type of FSLs they generate. One interesting example being
systems with “flat bands” being a huge $E=0$ degeneracy. Of course, one can
consider other models beyond those mentioned in this paper. Another note we
may make, here we use the Fock states, i.e. bare states, in order to construct
the FSL. Shortly, the full Hamiltonian has the form
$\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{int}}$, where the Fock states are the
eigenstates of $\hat{H}_{0}$, and $\hat{H}_{\mathrm{int}}$ constitute the FSL.
We could, generalize this to a situation where $\hat{H}_{0}$ is not diagonal
in the Fock basis, but instead in some other (dressed) basis. Take, for
example, the driven JC model,
$\hat{H}=\hat{H}_{\mathrm{JC}}+\eta\left(\hat{a}^{\dagger}+\hat{a}\right)$,
and consider the JC dressed states, rather than the bare Fock states. The
drive term will then constitute a new state space lattice (dressed state
lattice), which in this example becomes a 1D chain. As another example, the
quantum Rabi model could be written as
$\hat{H}_{\mathrm{R}}=\hat{H}_{\mathrm{JC}}+g\left(\hat{a}^{\dagger}\hat{\sigma}^{+}+\hat{a}\hat{\sigma}^{-}\right)$,
such that it is expressed as a JC Hamiltonian plus the counter rotating terms
Themis . In the JC dressed state basis, the resulting state space lattice
formed by the counter rotating terms will be two decoupled 1D parity chains.
Naturally, experimentally the bare Fock states are typically more relevant
from a perspective of detection, but we mention this to stress the greater
flexibility of the method by allowing for other bases.

##  Appendix A Equations of motion for the mean values of the three-mode boson
model

The multi-mode bosonic model (21) is quadratic and thereby solvable via a
Bogoliubov transformation. Moreover, as a quadratic model it follows that an
initial Gaussian state, e.g. any coherent or squeezed states, remains Gaussian
at all times. Thus, a Gaussian quantum phase space distribution stay localized
and follows the classical phase space trajectories. It is convenient to work
with the classical canonical variables for position

$\begin{split}\hat{x}=\frac{\left(\hat{a}^{\dagger}+\hat{a}\right)}{\sqrt{2}}\end{split}$
(38)

and momentum

$\begin{split}\hat{p}_{x}=\frac{\left(\hat{a}^{\dagger}-\hat{a}\right)}{\sqrt{2}}.\end{split}$
(39)

Using that the Heisenberg and Hamilton’s equations have the same structure for
a quadratic bosonic Hamiltonian, for the three-mode case we get the equations
of motion

$\begin{array}[]{lll}\dot{x}_{a}&=&-p_{b}-p_{c}\\\
\dot{x}_{b}&=&-p_{a}-\cos(\varphi)p_{c}-\sin(\varphi)x_{c},\\\
\dot{x}_{c}&=&-p_{a}-\cos(\varphi)p_{b}+\sin(\varphi)x_{b},\end{array}$ (40)

and

$\begin{array}[]{lll}\dot{p}_{a}&=&x_{b}+x_{c},\\\
\dot{p}_{b}&=&x_{a}+\cos(\varphi)x_{c}-\sin(\varphi)p_{c},\\\
\dot{p}_{c}&=&x_{a}+\cos(\varphi)x_{b}+\sin(\varphi)p_{b},\end{array}$ (41)

such that
$\langle\hat{n}_{\alpha}\rangle_{t}=\left(p_{\alpha}^{2}(t)+x_{\alpha}^{2}(t)\right)/2$.
Since a Gaussian state (apart from the vacuum) does not contain a fixed number
of bosons $N$, it means for the FSLs that it will populate different
triangular FSLs – one for each boson number $N$. However, the shape and
structure of the trajectories in the FSLs are invariant for different
$N$-values; the particle number only provides a scaling of the trajectories as
demonstrated in Fig. 9.

For the above equations to reproduce the correct trajectories of figs. 9,
where the initial state is a Fock state $|N,0,0\rangle$, one should use the
initial conditions
$(x_{a},x_{b},x_{c},p_{a},p_{b},p_{c})=(\sqrt{2N},0,0,0,0,0)$. We have
numerically conformed that these trajectories agree with those obtained from
diagonalizing the full $N$-boson Hamiltonian.

###### Acknowledgements.

The authors thank Themistoklis Mavrogordatos for helpful feedback on the
manuscript. We acknowledge financial support from VR-Vetenskapsråset (The
Swedish Research Council), and KAW (The Knut and Alice Wallenberg foundation).

## References

  * (1) A. Altland and B. Simmons, Condensed Matter Field Theory, (Cambridge University Press, 2010).
  * (2) Lewenstein M, Sanpera A, Ahufinger V, Damski B, Sen(De) A and Sen U 2007 Adv. Phys. 56 243; Bloch I, Dalibard J and Zwerger W 2008 Rev. Mod. Phys. 80 885; Zoller P and Gardiner C W The Quantum World Of Ultra-cold Atoms And Light - Book III: Ultra-cold Atoms (World Scientific, 2017).
  * (3) T. J. Sturges, T. McDermott, A. Buraczewski, W. R. Clements, J. J. Renema, S. W. Nam, T. Gerrits, A. Lita, W. S. Kolthammer, A. Eckstein, I. A. Walmsley, and M. Stobiska, NPJ Quant. Inf. 7, 91 (2021).
  * (4) O. Boada, A. Celi, J. I. Latorre, and M. Lewenstein Phys. Rev. Lett. 108, 133001 (2012).
  * (5) H. M. Price, T. Ozawa, and N. Goldman, Phys. Rev. A 95, 023607 (2017).
  * (6) A. Perez-Leija, H. Moya-Cessa, A. Szameit, and D. N. Christodoulides, Opt. Lett. 35, 2409 (2010).
  * (7) K. Tschernig, R. de J. León-Montiel, A. Pérez-Leija, and K. Busch, Phot. Res. 8, 1161 (2020).
  * (8) D.-W. Wang, H. Cai, R.-B. Liu, and M. O. Scully, Phys. Rev. Lett. 116, 220502 (2016).
  * (9) H. Cai and D.-W. Wang, National Sci. Rev. 8, nwaa196 (2021).
  * (10) E. T. Jaynes and F. W. Cummings, Jan 1963, Comparison of quantum and semiclassical radiation theories with application to the beam maser in Proceedings of the IEEE, vol. 51, no. 1, pp. 89-109, doi: 10.1109/PROC.1963.1664.
  * (11) J. Larson and T. Mavrogordatos, The Jaynes–Cummings Model and Its Descendants , (IOP Publishing, 2022).
  * (12) J. Casanova, G. Romero, I. Lizuain, J. J. García-Ripoll, and E. Solano, Phys. Rev. Lett. 105, 263603 (2010).
  * (13) Q.-T. Xie, S. Cui, J.-P. Cao, L. Amico, and H. Fan, Phys. Rev. X 4, 021046 (2014).
  * (14) M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
  * (15) N. D. Mermin, Physica A: Stat. Mech. App. 177, 561 (1991); H.-P. Breuer, D. Burgarth, and F. Petruccione, Phys. Rev. B 70, 045323 (2004).
  * (16) W. P. Su, J. R. Schrieffer and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * (17) M. Creutz, Phys. Rev. Lett. 83, 2636 (1999).
  * (18) C. Gardiner and P. Zoller, Quantum Noise, (Springer, 2004).
  * (19) A. Nunnenkamp, J. Koch, and S. M. Girvin, New J. Phys. 13, 095008 (2011).
  * (20) D. R. Hofstadter, Phys. Rev. B 14m 2239 (1976).
  * (21) M. Creutz, Phys. Rev. Lett. 83, 2636 (1999).
  * (22) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * (23) C. Weeks and M. Franz, Phys. Rev. B 82, 085310 (2010).
  * (24) J. Larson, E. Sjöqvist, and P. Öhberg. Conical Intersections in Physics, (Springer, 2020).
  * (25) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698, Solitons in Polyacetylene