# Polaritonic Hofstadter Butterfly and Cavity-Control of the Quantized Hall
Conductance

Vasil Rokaj<EMAIL_ADDRESS>Max Planck Institute for the
Structure and Dynamics of Matter, Center for Free Electron Laser Science,
22761 Hamburg, Germany  ITAMP, Harvard-Smithsonian Center for Astrophysics,
Cambridge, MA 02138, USA  Markus Penz  Department of Mathematics, University
of Innsbruck, Technikerstraße 13/7, A-6020 Innsbruck, Austria  Michael A.
Sentef  Max Planck Institute for the Structure and Dynamics of Matter, Center
for Free Electron Laser Science, 22761 Hamburg, Germany  Michael Ruggenthaler
Max Planck Institute for the Structure and Dynamics of Matter, Center for Free
Electron Laser Science, 22761 Hamburg, Germany  Angel Rubio
<EMAIL_ADDRESS>Max Planck Institute for the Structure and Dynamics of
Matter, Center for Free Electron Laser Science, 22761 Hamburg, Germany  Center
for Computational Quantum Physics (CCQ), Flatiron Institute, 162 Fifth Avenue,
New York NY 10010

###### Abstract

In a previous work [Phys. Rev. Lett. 123, 047202 (2019)] a translationally
invariant framework called quantum-electrodynamical Bloch (QED-Bloch) theory
was introduced for the description of periodic materials in homogeneous
magnetic fields and strongly coupled to the quantized photon field in the
optical limit. For such systems, we show that QED-Bloch theory predicts the
existence of fractal polaritonic spectra as a function of the cavity coupling
strength. In addition, for the energy spectrum as a function of the relative
magnetic flux we find that a terahertz cavity can modify the standard
Hofstadter butterfly. In the limit of no quantized photon field, QED-Bloch
theory captures the well-known fractal spectrum of the Hofstadter butterfly
and can be used for the description of 2D materials in strong magnetic fields,
which are of great experimental interest. As a further application, we
consider Landau levels under cavity confinement and show that the cavity
alters the quantized Hall conductance and that the Hall plateaus are modified
as $\sigma_{xy}=e^{2}\nu/h(1+\eta^{2})$ by the light-matter coupling $\eta$.
Most of the aforementioned phenomena should be experimentally accessible and
corresponding implications are discussed.

###### pacs:

Valid PACS appear here

##  I Introduction

The study of two-dimensional systems perpendicular to a strong homogeneous
magnetic field has given rise to a plethora of macroscopic quantum phenomena
known as the quantum Hall effects von Klitzing  _et al._ (2020). This
fundamental branch of condensed matter physics was ignited by the discovery of
the integer Hall effect in 1980 Klitzing  _et al._ (1980), where the
macroscopic Hall conductance $\sigma_{xy}$ exhibits quantized plateaus whose
value depend solely on the fundamental charge $e$, Planck’s constant $h$, and
the filling factor $\nu$ in the picture of non-interacting Landau levels,
$\sigma_{xy}=e^{2}\nu/h$. In the four decades after the fundamental discovery
of the integer Hall effect, a great number of related phenomena like the
fractional quantum Hall effect Stormer  _et al._ (1999); Laughlin (1983), the
quantum spin Hall effect Kane and Mele (2005), the quantum anomalous Hall
effect Chang  _et al._ (2013) and more recently the light-induced anomalous
Hall effect McIver  _et al._ (2020) have been observed and studied
theoretically, as they provide a unique platform to address strongly
correlated electronic phases, topology and strong light-matter phenomena. All
these exciting developments have been reviewed beautifully in an article
celebrating the 40 years anniversary of the quantum Hall effect von Klitzing
_et al._ (2020).

In this cornucopia of phenomena offered by the quantum Hall setting there is
another important phenomenon which stands out due to its fundamental nature,
aesthetic beauty, and the connections that it provides between different
branches of mathematics and physics, namely the Hofstadter butterfly
Hofstadter (1976). The Hofstadter butterfly is a fractal pattern that
describes the energies of electrons on a periodic lattice perpendicular to a
homogeneous magnetic field as a function of the relative magnetic flux
$\Phi/\Phi_{0}$, where $\Phi$ is the magnetic flux through the fundamental
unit cell of the lattice and $\Phi_{0}=h/e$ is the magnetic flux quantum. It
was predicted by Hofstadter in 1976 Hofstadter (1976) and it was proven to be
a fractal by Avila and Jitomirskaya Avila and Jitomirskaya (2006). In recent
years, due to the advent of Moiré materials Andrei  _et al._ (2021); Kennes
_et al._ (2021), it has become possible to probe it experimentally, and
signatures of the fractal spectrum have been observed in the magnetotransport
properties of Moiré systems Dean  _et al._ (2013); Wang  _et al._ (2015);
Forsythe  _et al._ (2018). Further, the physics of the Hofstadter butterfly
has also been realized with ultracold atoms in optical lattices Aidelsburger
_et al._ (2013); Miyake  _et al._ (2013).

Another pillar of modern quantum physics is quantum electrodynamics (QED),
which describes the interaction of charged particles with photons Cohen-
Tannoudji  _et al._ (1997); Weinberg (2005); Spohn (2004). In the last decade
there has been a great interest in the regime of strong and ultrastrong light-
matter interactions Kockum  _et al._ (2019), where light and matter lose their
individual character and form hybrid quasiparticles known as polaritons Basov
_et al._ (2021). Many different platforms and routes have been explored to
reach strong light-matter coupling and several unprecedented phenomena
involving polaritonic states have been observed. Modifications of chemical
properties and chemical reactions have been achieved through coupling to
vacuum fields in polaritonic chemistry Ebbesen (2016); Hutchison  _et al._
(2013, 2012); Orgiu  _et al._ (2015); Feist  _et al._ (2017); Galego  _et al._
(2016); Flick  _et al._ (2017); Schäfer  _et al._ (2019). Cavity control of
excitons has been studied Latini  _et al._ (2019); Förg  _et al._ (2019) and
exciton-polariton condensation has been achieved Kasprzak  _et al._ (2006);
Keeling and Kena-Cohen (2020). It has been suggested that light-matter
interactions modify the electron-phonon coupling and the critical temperature
of superconductors Schlawin  _et al._ (2019); Cotleţ  _et al._ (2016); Sentef
_et al._ (2018); Curtis  _et al._ (2019), with the first experimental evidence
already at hand Thomas  _et al._ (2019). Further, the implications of coupling
to chiral electromagnetic fields is currently investigated Hübener  _et al._
(2021); Petersen  _et al._ (2014); Zhang  _et al._ (2019); Lodahl  _et al._
(2017); Wang  _et al._ (2019), and the possibility of cavity-induced
ferroelectric phases has been proposed Latini  _et al._ (2021); Ashida  _et
al._ (2020).

As a synthesis of QED and the quantum Hall setting, quantum Hall systems under
cavity confinement have been studied both experimentally and theoretically, in
the integer Hagenmüller  _et al._ (2010); Rokaj  _et al._ (2019); Keller  _et
al._ (2020); Scalari  _et al._ (2012); Li  _et al._ (2018) and the correlated
fractional Ravets  _et al._ (2018); Smolka  _et al._ (2014) regime, and
ultrastrong coupling to the photon field and modifications of their transport
properties Paravicini-Bagliani  _et al._ (2019) have been demonstrated.
Recently, a theoretical mechanism for a cavity-mediated hopping in the integer
regime was also proposed Ciuti (2021) and the breakdown of the topological
protection of the integer Hall effect due to cavity vacuum fields was
demonstrated experimentally Appugliese  _et al._ (2022).

In this article, we focus on this emerging field of quantum Hall systems
strongly coupled to the quantized photon-field originating from a cavity (see
Fig. 1). To be more precise, we investigate the modification of two particular
phenomena due to strong vacuum fluctuations inside a cavity: (i) the
Hofstadter butterfly and (ii) the quantization of the Hall conductance in the
integer regime. To describe these quantum Hall systems and phenomena in the
cavity setting we employ the recently-introduced quantum-electrodynamical
Bloch (QED-Bloch) theory Rokaj  _et al._ (2019); Rokaj (2021). QED-Bloch
theory solves the long-standing problem of broken translational invariance due
to an external magnetic field and provides a translationally symmetric
framework for periodic systems in the presence of a homogeneous magnetic field
and strongly coupled to the quantized photon field in the optical limit (or
dipole approximation). Our main findings are:

  * •

Polaritonic Hofstadter Butterfly.—For a 2D periodic system perpendicular to a
homogeneous magnetic field and under cavity confinement we find that for the
energy spectrum a self-similar pattern emerges as a function of the light-
matter coupling $\eta$ as depicted in Fig. 2(d). This is an extension of the
standard Hofstadter butterfly Hofstadter (1976) to the polaritonic (light-
matter) setting of cavity QED which introduces the concept of polaritonic
fractal spectra. We call this phenomenon the polaritonic Hofstadter butterfly.
In addition, we compute the energy spectrum for a periodic material as a
function of the relative magnetic flux and find that for a terahertz cavity
Scalari  _et al._ (2012); Paravicini-Bagliani  _et al._ (2019); Li  _et al._
(2018) the standard Hofstadter butterfly gets modified due to the strong
vacuum fluctuations of the photon field (see Fig. 3(c)).

  * •

Periodic Materials in Homogeneous Magnetic Fields.—In the semi-classical limit
of no quantized field our QED-Bloch theory recovers the standard results of
condensed matter systems in strong magnetic fields like the Hofstadter
butterfly Hofstadter (1976), the Landau levels Landau and Lifshitz (1997) and
the quantization of the Hall conductance Klitzing  _et al._ (1980); Laughlin
(1981) and provides a first-principles framework for the description of
periodic materials in strong magnetic fields. In addition, a dual relation
between the minimal-coupling Hamiltonian and the tight-binding models with the
Peierls substitution is established (see Fig. 7).

  * •

Cavity Modification of the Integer Hall Effect.—For a 2D electron gas
consisting of Landau levels strongly coupled to the cavity field we find that
the quantized Hall conductance gets modified. The Hall plateaus inside the
cavity depend on the light-matter coupling constant $\eta$ as
$\sigma_{xy}=e^{2}\nu/h(1+\eta^{2})$. This modification is a consequence of
the formation of hybrid quasiparticle states between the Landau levels and the
cavity photons known as Landau polaritons Scalari  _et al._ (2012);
Paravicini-Bagliani  _et al._ (2019); Rokaj  _et al._ (2019); Hagenmüller  _et
al._ (2010).

Figure 1: Cartoon depiction of a two-dimensional periodic material confined
inside a cavity with mirrors of length $L$ and area $S=L^{2}$. The distance
between the cavity mirrors is $L_{z}$. The whole system is placed
perpendicular to a classical homogeneous magnetic field
$\mathbf{B}_{\textrm{ext}}$. We note that typically in experimental setups the
space between the 2D material and the cavity is filled with a dielectric
medium.

Outline of the Article.—In Section II we recapitulate the basic steps of the
construction of QED-Bloch theory and the restoration of translational
symmetry. In Section III we construct the QED-Bloch ansatz which we use for
the description of periodic materials in strong magnetic fields and strongly
coupled to the quantized photon field and we derive the QED-Bloch central
equation. In Section IV we show that for periodic materials in strong magnetic
fields and under cavity confinement there exists a polaritonic fractal
spectrum as a function of the light-matter coupling constant. In Section V we
demonstrate that there is a duality between the minimal-coupling Hamiltonian
and the tight-binding models with the Peierls phase. In Section VI we compute
the Hall conductance for non-interacting Landau levels coupled to the cavity
field and we show that the cavity modifies the quantization of the Hall
conductance. Finally, in Section VII we conclude and highlight the future
perspectives of this work.

##  II QED-Bloch Theory

The description of periodic materials in the presence of a homogeneous
magnetic field has been a long-standing problem for condensed matter physics.
The problem arises due to the fact that, although the magnetic field is
homogeneous throughout the whole material, in the minimally-coupled
Schrödinger equation the electrons couple to the spatially inhomogeneous
vector potential of the electromagnetic field. Thus, translational symmetry is
broken and Bloch theory is not applicable.

Recently, a translationally invariant, quantum-electrodynamical framework for
the description of periodic solids in homogeneous magnetic fields was
introduced by the authors Rokaj  _et al._ (2019). Due to the fact that in this
framework translational symmetry is restored in the higher-dimensional
configuration space of both electrons and photons, in which Bloch’s theorem
can still be applied, it was named quantum-electrodynamical Bloch (QED-Bloch)
theory. Before we proceed with the application of this framework we would like
to briefly recapitulate the basic steps in the construction of QED-Bloch
theory.

Our starting point is the Pauli-Fierz Hamiltonian for $N$ interacting
electrons in a periodic potential and in the presence of a classical,
homogeneous magnetic field and further coupled to the quantized photon field
Cohen-Tannoudji  _et al._ (1997); Spohn (2004); Rokaj  _et al._ (2018),

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum\limits^{N}_{j=1}\left[\frac{1}{2m_{e}}\left(\mathrm{i}\hbar\mathbf{\nabla}_{j}+e\hat{\mathbf{A}}(\mathbf{r}_{j})+e\mathbf{A}_{\textrm{ext}}(\mathbf{r}_{j})\right)^{2}+v_{\textrm{ext}}(\mathbf{r}_{j})\right]$
(1) $\displaystyle+$
$\displaystyle\frac{1}{4\pi\epsilon_{0}}\sum\limits^{N}_{j<k}\frac{e^{2}}{|\mathbf{r}_{j}-\mathbf{r}_{k}|}+\sum\limits_{\bm{\kappa},\lambda}\hbar\omega(\bm{\kappa})\left[\hat{a}^{\dagger}_{\bm{\kappa},\lambda}\hat{a}_{\bm{\kappa},\lambda}+\frac{1}{2}\right],$

where we neglected the Pauli (Stern-Gerlach) term
$\hat{\bm{\sigma}}\cdot\hat{\mathbf{B}}(\mathbf{r})$ as it is commonly done
for the description of the Hofstadter butterfly Hofstadter (1976) and the
Landau levels in the quantum Hall effect Laughlin (1981); Landau and Lifshitz
(1997). Here $\mathbf{A}_{\textrm{ext}}(\mathbf{r})$ is the external vector
potential which gives rise to the homogeneous magnetic field
$\mathbf{B}_{\textrm{ext}}=\nabla\times\mathbf{A}_{\textrm{ext}}(\mathbf{r})=\mathbf{e}_{z}B$
in $z$-direction and $\mathbf{A}_{\textrm{ext}}(\mathbf{r})$ is chosen to be
in the Landau gauge $\mathbf{A}_{\textrm{ext}}(\mathbf{r})=-\mathbf{e}_{x}By$
Landau and Lifshitz (1997). The quantized vector potential
$\hat{\mathbf{A}}(\mathbf{r})$ of the electromagnetic field in the Coulomb
gauge is Spohn (2004); Greiner and Reinhardt (1996)

$\hat{\mathbf{A}}(\mathbf{r})=\sum_{\bm{\kappa},\lambda}\sqrt{\frac{\hbar}{\epsilon_{0}V2\omega(\bm{\kappa})}}\left[\hat{a}_{\bm{\kappa},\lambda}\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})+\hat{a}^{\dagger}_{\bm{\kappa},\lambda}\mathbf{S}^{*}_{\bm{\kappa},\lambda}(\mathbf{r})\right],$
(2)

with $\bm{\kappa}=(\kappa_{x},\kappa_{y},\kappa_{z})$ the wave vectors of the
photon field, $\omega(\bm{\kappa})=c|\bm{\kappa}|$ the allowed frequencies in
the quantization volume $V=L^{2}L_{z}$, $\lambda=1,2$ the two transverse
polarization directions and $\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})$ the
vector-valued mode functions, chosen such that the Coulomb gauge is satisfied,
$\nabla\cdot\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})=0$ Spohn (2004);
Greiner and Reinhardt (1996). In order for the mode functions
$\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})$ to satisfy the boundary
conditions of the cavity, the wave vectors of the photon field can only take
the values $\bm{\kappa}=(\kappa_{x},\kappa_{y},\kappa_{z})=(2\pi n_{x}/L,2\pi
n_{y}/L,\pi n_{z}/L_{z})$ with
$\mathbf{n}=(n_{x},n_{y},n_{z})\in\mathbb{Z}^{3}$. The operators
$\hat{a}_{\bm{\kappa},\lambda}$ and $\hat{a}^{\dagger}_{\bm{\kappa},\lambda}$
are the annihilation and creation operators of the photon field and obey the
bosonic commutation relations
$[\hat{a}_{\bm{\kappa},\lambda},\hat{a}^{\dagger}_{\bm{\kappa}^{\prime},\lambda^{\prime}}]=\delta_{\bm{\kappa}\bm{\kappa}^{\prime}}\delta_{\lambda\lambda^{\prime}}$.
We note that the annihilation and creation operators can also be defined in
terms of the displacement coordinates $q_{\bm{\kappa},\lambda}$ and their
conjugate momenta $\partial/\partial q_{\bm{\kappa},\lambda}$ as
$\hat{a}_{\bm{\kappa},\lambda}=\frac{1}{\sqrt{2}}\left(q_{\bm{\kappa},\lambda}+\partial/\partial
q_{\bm{\kappa},\lambda}\right)\;\textrm{and}\;\hat{a}^{\dagger}_{\bm{\kappa},\lambda}=\frac{1}{\sqrt{2}}\left(q_{\bm{\kappa},\lambda}-\partial/\partial
q_{\bm{\kappa},\lambda}\right)$.

It is clear that in the Pauli-Fierz Hamiltonian both the quantized field
$\hat{\mathbf{A}}(\mathbf{r})$ defined in Eq. (2) as well as the external
vector potential $\mathbf{A}_{\textrm{ext}}(\mathbf{r})$ that induces the
perpendicular, homogeneous magnetic field break translational symmetry because
they are spatially inhomogeneous. This implies that Bloch’s theorem cannot be
applied for the general Pauli-Fierz Hamiltonian. However, a way to restore
translational symmetry in the full electron-photon configuration space was
found after performing the long-wavelength limit (or dipole approximation) for
the quantized field $\hat{\mathbf{A}}(\mathbf{r})$ Rokaj  _et al._ (2019). In
the long-wavelength limit Rokaj  _et al._ (2018); Faisal (1987), which has
been proven adequate for cavity QED systems Schäfer  _et al._ (2019);
Ruggenthaler  _et al._ (2018), the mode functions
$\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})$ become spatially independent
vectors
$\mathbf{S}_{\bm{\kappa},\lambda}(\mathbf{r})=\bm{\varepsilon}_{\lambda}(\bm{\kappa})$
and satisfy the condition
$\bm{\varepsilon}_{\lambda}(\bm{\kappa})\cdot\bm{\varepsilon}_{\lambda^{\prime}}(\bm{\kappa})=\delta_{\lambda\lambda^{\prime}}$.
In addition, we keep only a single-mode of the quantized field and the
polarization of the dipolar quantized field is chosen to be parallel to the
external vector potential. Under these assumptions the Pauli-Fierz Hamiltonian
simplifies to

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum\limits^{N}_{j=1}\left[\frac{1}{2m_{e}}\left(\mathrm{i}\hbar\mathbf{\nabla}_{j}+e\hat{\mathbf{A}}+e\mathbf{A}_{\textrm{ext}}(\mathbf{r}_{j})\right)^{2}+v_{\textrm{ext}}(\mathbf{r}_{j})\right]$
(3) $\displaystyle+$
$\displaystyle\frac{1}{4\pi\epsilon_{0}}\sum\limits^{N}_{j<k}\frac{e^{2}}{|\mathbf{r}_{j}-\mathbf{r}_{k}|}+\hbar\omega\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right)$

and respectively the quantized field in the dipole approximation is Spohn
(2004)

$\hat{\mathbf{A}}=\left(\frac{\hbar}{\epsilon_{0}V}\right)^{\frac{1}{2}}\frac{\mathbf{e}_{x}}{\sqrt{2\omega}}\left(\hat{a}+\hat{a}^{\dagger}\right).$
(4)

This single-mode quantized field describes the cavity mode which is strongly
coupled to the matter system. However, in Appendix C we will also take into
account the effect of many modes for the integer Hall effect inside a cavity.

Translational symmetry can be restored for the combined electron-photon system
in the optical limit (or dipole approximation) $\omega\to 0$ Rokaj  _et al._
(2019).111We note that the limit $\omega\rightarrow 0$ does not imply that the
photonic Hamiltonian $\hat{H}_{p}$ is identically zero. To be able to perform
consistently this limit one needs to treat exactly and non-perturbatively the
$\hat{\mathbf{A}}^{2}$ term of the quantized mode. To do so, we isolate the
purely photonic part of the Pauli-Fierz Hamiltonian, namely the part that
depends only on the photonic annihilation and creation operators $\hat{a}$ and
$\hat{a}^{\dagger}$,

$\hat{H}_{p}=\hbar\omega\left(\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right)+\frac{Ne^{2}}{2m_{e}}\hat{\mathbf{A}}^{2}.$
(5)

The purely photonic part can be brought into the form of the standard harmonic
oscillator as follows. The diamagnetic $\hat{\mathbf{A}}^{2}$ term
renormalizes the photon frequency $\omega$ by introducing the diamagnetic
shift $\omega_{p}=\sqrt{e^{2}n_{e}/m_{e}\epsilon_{0}}$ which depends on the
electron density $n_{e}$. Then we introduce the scaled coordinate
$u=q\sqrt{\widetilde{\omega}/\omega}$ where the dressed photon frequency
$\widetilde{\omega}$ is defined by
$\widetilde{\omega}^{2}=\omega^{2}+\omega^{2}_{p}$. In terms of the new
coordinate $u$ and its conjugate momentum $\partial_{u}$ the photonic part
$\hat{H}_{p}$ takes the form of a simple harmonic oscillator with frequency
$\widetilde{\omega}$,

$\hat{H}_{p}=\frac{\hbar\widetilde{\omega}}{2}\left(-\frac{\partial^{2}}{\partial
u^{2}}+u^{2}\right).$ (6)

The quantized photon field in terms of $u$ is

$\hat{\mathbf{A}}=\sqrt{\frac{\hbar}{\epsilon_{0}V\widetilde{\omega}}}u\mathbf{e}_{x}$
(7)

and performing now the optical limit $\omega\rightarrow 0$ we find that the
dressed frequency $\widetilde{\omega}$ simply goes to the diamagnetic
frequency $\omega_{p}$. Substituting the expressions for the purely photonic
part
$\hat{H}_{p}=\left(\hbar\omega_{p}/2\right)\left(-\partial^{2}_{u}+u^{2}\right)$
and the vector potential

$\hat{\mathbf{A}}=\sqrt{\frac{\hbar}{\epsilon_{0}V\omega_{p}}}u\mathbf{e}_{x},$
(8)

back into (3) we obtain the Pauli-Fierz Hamiltonian in the optical limit

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum\limits^{N}_{j=1}\left[\frac{1}{2m_{e}}\left(\mathrm{i}\hbar\mathbf{\nabla}_{j}+e\hat{\mathbf{A}}+e\mathbf{A}_{\textrm{ext}}(\mathbf{r}_{j})\right)^{2}+v_{\textrm{ext}}(\mathbf{r}_{j})\right]$
(9) $\displaystyle+$
$\displaystyle\frac{1}{4\pi\epsilon_{0}}\sum\limits^{N}_{j<k}\frac{e^{2}}{|\mathbf{r}_{j}-\mathbf{r}_{k}|}-\frac{\hbar\omega_{p}}{2}\frac{\partial^{2}}{\partial
u^{2}}\;.$

We note that in order to obtain the above expression for the Hamiltonian we
used that $\hbar\omega_{p}u^{2}/2=Ne^{2}\hat{\mathbf{A}}^{2}/2m_{e}$.

Let us now check the translational properties of the Hamiltonian given by Eq.
(9). For simplicity and to avoid specifying a certain lattice geometry, we
will consider the case of no external potential,
$v_{\textrm{ext}}(\mathbf{r})=0$. The Hamiltonian is not periodic in the
electronic coordinates because $\mathbf{A}_{\textrm{ext}}(\mathbf{r})$ is
linear in $y$. However, the Hamiltonian in the optical limit is periodic under
a generalized translation in the full electronic plus photonic configuration
space,

$(\mathbf{r}_{j},u)\longrightarrow\left(\mathbf{r}_{j}+\mathbf{a},u+Ba_{y}\sqrt{\epsilon_{0}V\omega_{p}/\hbar}\right),$
(10)

where $\mathbf{a}=(a_{x},a_{y},a_{z})\in\mathbb{R}^{3}$ arbitrary. This is
true because the total vector potential
$\hat{\mathbf{A}}_{\textrm{tot}}=\hat{\mathbf{A}}+\mathbf{A}_{\textrm{ext}}(\mathbf{r}_{j})$
is invariant under the above generalized translation and obviously the kinetic
terms of the quantized mode and of the electrons are also invariant, as well
as the Coulomb interaction. The fact that the Hamiltonian in Eq. (9) is
invariant under the continuous translations of Eq. (10) implies that $\hat{H}$
will also be invariant under Bravais lattice translations in the case of a
periodic potential.

##  III Effective Hamiltonian & QED-Bloch Expansion

Having restored translational symmetry, our goal now is to go one step further
and construct a Bloch-type ansatz in the polaritonic (electronic plus
photonic) configuration space and to derive a Bloch-type central equation for
the description of solids in a classical, homogeneous magnetic field coupled
to a quantized electromagnetic field.

To make the problem tractable, instead of treating the unfeasible many-body
interacting Hamiltonian of Eq. (9), we will employ the independent electron
approximation which is commonly used in condensed matter physics. We note that
this independent-electron approach is consistent with Bloch theory, which is
not a theory of a single electron in a periodic potential but rather of many
non-interacting electrons.

To incorporate the fact that the charged particles couple collectively to the
photon field, we will use an effective electron density $n_{e}$ and we will
scale the strength of the quantized field $\hat{\mathbf{A}}$ by the square
root of the number of charges $\sqrt{N}$,

$\hat{\mathbf{A}}\longrightarrow\hat{\mathbf{\mathcal{A}}}=\sqrt{N}\hat{\mathbf{A}}=\sqrt{\frac{\hbar\omega_{p}m_{e}}{e^{2}}}u\mathbf{e}_{x}.$
(11)

Keeping an effective electron density and scaling the field by $\sqrt{N}$
allows us to capture the back-reaction of matter to the photon field and to
describe the emergence of novel quasiparticle excitations known as Landau
polaritons Rokaj  _et al._ (2019). We would also like to mention that the
scaling of light-matter interaction by $\sqrt{N}$ is the standard argument for
the description of collective coupling in the few-level models of quantum
optics Ribeiro  _et al._ (2021); Dicke (1954); Garraway (2011); Kockum  _et
al._ (2019). For the inclusion of any further effects, like exchange and
correlation effects, one would need the addition of effective fields as
introduced in quantum-electrodynamical density functional theory (QEDFT)
Ruggenthaler  _et al._ (2014); Ruggenthaler (2015); Tokatly (2013); Schäfer
_et al._ (2021).

Upon these assumptions, we obtain the following single-particle effective
Hamiltonian,

$\hat{H}_{\mathrm{eff}}=\frac{1}{2m_{e}}\left(\mathrm{i}\hbar\mathbf{\nabla}+e\hat{\mathbf{\mathcal{A}}}+e\mathbf{A}_{\textrm{ext}}(\mathbf{r})\right)^{2}-\frac{\hbar\omega_{p}}{2}\frac{\partial^{2}}{\partial
u^{2}}+v_{\textrm{ext}}(\mathbf{r}).$ (12)

which was already proposed by the authors in Ref. Rokaj  _et al._ (2019) and
was applied successfully to the description of Landau polariton systems.
Before we continue, we would like to specify the geometries in which we are
interested in this article.

###  III.1 Setting the Geometry

Our aim here is to treat all possible 2D geometries of periodic structures.
The external potential in a solid is assumed periodic
$v_{\textrm{ext}}(\mathbf{r})=v_{\textrm{ext}}(\mathbf{r}+\mathbf{R}_{\mathbf{n}})$
where $\mathbf{R}_{\mathbf{n}}$ is a Bravais lattice vector with
$\mathbf{n}=(n,m)\in\mathbb{Z}^{2}$. The Bravais lattice vectors in general
are $\mathbf{R}_{\mathbf{n}}=n\mathbf{a}_{1}+m\mathbf{a}$, where
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ are the primitive vectors which lie in
different directions and span the 2D lattice. Without loss of generality we
can choose the vector $\mathbf{a}_{1}$ to be in the $x$-direction
$\mathbf{a}_{1}=a_{1}\mathbf{e}_{x}$. The second primitive vector in this case
is
$\mathbf{a}_{2}=a_{2}\cos\theta\mathbf{e}_{x}+a_{2}\sin\theta\mathbf{e}_{y}$,
where $\theta$ is the angle between the vectors $\mathbf{a}_{1}$ and
$\mathbf{a}_{2}$. Thus, the Bravais lattice vectors are

$\mathbf{R}_{\mathbf{n}}=\left(na_{1}+ma_{2}\cos\theta\right)\mathbf{e}_{x}+ma_{2}\sin\theta\mathbf{e}_{y}.$
(13)

Then, the reciprocal lattice vectors are
$\mathbf{G}_{\mathbf{n}}=n\mathbf{b}_{1}+m\mathbf{b}_{2}$ with
$\mathbf{n}=(n,m)\in\mathbb{Z}^{2}$. The defining relation for the vectors
$\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ is
$\mathbf{b}_{i}\cdot\mathbf{a}_{j}=2\pi\delta_{ij}\;\;\textrm{with}\;\;i,j=1,2$
Aschroft and Mermin (1976); Callaway (1964). With the given choice of
primitive vectors the reciprocal primitive vectors are

$\mathbf{b}_{1}=\frac{2\pi}{a_{1}}\mathbf{e}_{x}-\frac{2\pi\cos\theta}{a_{1}\sin\theta}\mathbf{e}_{y}\;\;\textrm{and}\;\;\mathbf{b}_{2}=\frac{2\pi}{a_{2}\sin\theta}\mathbf{e}_{y}.$
(14)

Thus, the reciprocal lattice vectors are

$\mathbf{G}_{\mathbf{n}}=\frac{2\pi n}{a_{1}}\mathbf{e}_{x}+\left(\frac{2\pi
m}{a_{2}\sin\theta}-\frac{2\pi
n\cos\theta}{a_{1}\sin\theta}\right)\mathbf{e}_{y},$ (15)

which for convenience we will write as

$\displaystyle\mathbf{G}_{\mathbf{n}}=\left(G^{x}_{n},G_{\mathbf{n}}\right)\;\;\textrm{with}\;\;G_{\mathbf{n}}=-\frac{G^{x}_{n}\cos\theta}{\sin\theta}+\frac{G^{y}_{m}}{\sin\theta}$
(16) $\displaystyle\textrm{and}\;\;G^{x}_{n}=\frac{2\pi
n}{a_{1}},\;\;\;G^{y}_{m}=\frac{2\pi m}{a_{2}}.$

With these choices we have defined our geometrical setting and we have made
clear how all possible 2D Bravais lattices can be described.

###  III.2 Polaritonic Coordinates

To continue, we introduce the cyclotron frequency $\omega_{c}=eB/m_{e}$ and we
expand the covariant kinetic term of the effective Hamiltonian into

$\displaystyle\hat{H}_{\mathrm{eff}}=$ $\displaystyle-$
$\displaystyle\frac{\hbar^{2}}{2m_{e}}\nabla^{2}+\mathrm{i}\hbar\mathbf{e}_{x}\left(u\sqrt{\hbar\omega_{p}/m_{e}}-y\omega_{c}\right)\cdot\nabla$
$\displaystyle+$ $\displaystyle
v_{\textrm{ext}}(\mathbf{r})+\frac{m_{e}}{2}\left(u\sqrt{\hbar\omega_{p}/m_{e}}-y\omega_{c}\right)^{2}-\frac{\hbar\omega_{p}}{2}\frac{\partial^{2}}{\partial
u^{2}}.$

This effective Hamiltonian is invariant under a translation that acts on both
the electronic and photonic coordinates for any 2D Bravais lattice vector
$\mathbf{R}_{\mathbf{n}}$,

$(\mathbf{r},u)\longrightarrow\left(\mathbf{r}+\mathbf{R}_{\mathbf{n}},u+ma_{2}\sin\theta\omega_{c}\sqrt{m_{e}/\hbar\omega_{p}}\right).$
(18)

To describe properly this symmetry in the polaritonic space, we will switch to
a new set of coordinates. For that purpose, we introduce the relative distance
and center-of-mass coordinates between rescaled versions of $u$ and $y$,

$w=\frac{m_{p}u\sqrt{\hbar\omega_{p}/m_{e}}+m_{c}\omega_{c}y}{\sqrt{2}M},\;v=\frac{u\sqrt{\hbar\omega_{p}/m_{e}}-\omega_{c}y}{\sqrt{2}},$
(19)

and the Hamiltonian $\hat{H}_{\mathrm{eff}}$ simplifies to

$\displaystyle\hat{H}_{\mathrm{eff}}$
$\displaystyle=-\frac{\hbar^{2}}{2m_{e}}\frac{\partial^{2}}{\partial
x^{2}}-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
w^{2}}+v_{\textrm{ext}}(\mathbf{r})$ (20)
$\displaystyle+\textrm{i}\hbar\sqrt{2}v\frac{\partial}{\partial
x}+m_{e}v^{2}-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial v^{2}},$

with the mass parameters $M,\mu,m_{p}$ and $m_{c}$ being

$\displaystyle m_{p}$
$\displaystyle=\frac{m_{e}}{\omega^{2}_{p}},\;\;\;m_{c}=\frac{m_{e}}{\omega^{2}_{c}},$
(21) $\displaystyle M$
$\displaystyle=\frac{m_{p}+m_{c}}{2}\;\;\;\textrm{and}\;\;\;\mu=\frac{m_{p}m_{c}}{M}.$

Furthermore, by performing a square completion, the effective Hamiltonian can
be written in the compact form

$\displaystyle\hat{H}_{\mathrm{eff}}=$
$\displaystyle-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
w^{2}}-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial v^{2}}$ (22)
$\displaystyle+\frac{\mu\Omega^{2}}{2}\left(v+\frac{\textrm{i}\hbar}{\sqrt{2}m_{e}}\frac{\partial}{\partial
x}\right)^{2}+v_{\textrm{ext}}(\mathbf{r}),$

where the dressed frequency $\Omega$ is defined by

$\Omega^{2}=\frac{2m_{e}}{\mu}=\omega^{2}_{p}+\omega^{2}_{c}.$ (23)

The original electronic vector $\mathbf{r}=(x,y)$ in the new polaritonic
coordinate system is

$\mathbf{r}=(x,y)=\left(x,\frac{w}{\sqrt{2}\omega_{c}}-\frac{m_{p}v}{\sqrt{2}M\omega_{c}}\right).$
(24)

It is important to note that the coordinates $v$ and $w$ are independent,
because their respective position and momentum operators commute. Moreover, we
would like to emphasize that since the external potential is periodic, it can
be written in terms of a Fourier series, which in terms of the polaritonic
coordinates is

$v_{\textrm{ext}}(\mathbf{r})=\sum_{\mathbf{n}}V_{\mathbf{n}}\mathrm{e}^{\textrm{i}\mathbf{G}_{\mathbf{n}}\cdot\mathbf{r}}=\sum_{\mathbf{n}}V_{\mathbf{n}}\mathrm{e}^{\textrm{i}\mathbf{G}^{w}_{\mathbf{n}}\cdot\mathbf{r}_{w}}\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}}v},$
(25)

where

$\displaystyle
G^{v}_{\mathbf{n}}=\frac{m_{p}G_{\mathbf{n}}}{\sqrt{2}M\omega_{c}},\quad\mathbf{r}_{w}=(x,w)\;\;\textrm{and}$
(26)
$\displaystyle\mathbf{G}^{w}_{\mathbf{n}}=(G^{x}_{n},G^{w}_{\mathbf{n}})=(G^{x}_{n},G_{\mathbf{n}}/\sqrt{2}\omega_{c}).$

###  III.3 QED-Bloch Expansion

The Hamiltonian $\hat{H}_{\mathrm{eff}}$ of Eq. (22) is invariant under
translations in the polaritonic configuration space

$(x,w)\longrightarrow(x+na_{1}+ma_{2}\cos\theta,w+m\sqrt{2}\omega_{c}a_{2}\sin\theta).$
(27)

This implies that we can use Bloch’s theorem in the $(x,w)$ plane.
Consequently, the eigenfunctions of $\hat{H}_{\mathrm{eff}}$ can be written
with the ansatz

$\Psi_{\mathbf{k}}(\mathbf{r}_{w},v)=\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}_{w}}U^{\mathbf{k}}(\mathbf{r}_{w},v)$
(28)

where $\mathbf{r}_{w}=(x,w)$. Here the function
$U^{\mathbf{k}}(\mathbf{r}_{w},v)$ is periodic under the translations in the
polaritonic space defined in Eq. (27). The crystal momentum
$\mathbf{k}=(k_{x},k_{w})$ corresponds to $\mathbf{r}_{w}$ and $k_{w}$ is a
polaritonic quantum number. Note that the polaritonic unit cell in the
$w$-direction scales linearly with the strength of the magnetic field. The
same feature appears also for the usual magnetic unit cell, but in this case
only field strengths which generate a rational magnetic flux through a unit
cell are allowed Kohmoto (1985). Contrary to that, the polaritonic unit cell
puts no restrictions on the strength of the magnetic field.

Since the function $U^{\mathbf{k}}(\mathbf{r}_{w},v)$ is periodic in
$\mathbf{r}_{w}$ we expand it in a Fourier series in $\mathbf{r}_{w}$, while
for the $v$-dependent part we consider a generic function,

$\Psi_{\mathbf{k}}(\mathbf{r}_{w},v)=\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}_{w}}\sum_{\mathbf{n}}U^{\mathbf{k}}_{\mathbf{n}}\mathrm{e}^{\mathrm{i}\mathbf{G}^{w}_{\mathbf{n}}\cdot\mathbf{r}_{w}}\phi^{\mathbf{k}}_{\mathbf{n}}(v),$
(29)

where $\mathbf{G}^{w}_{\mathbf{n}}=(G^{x}_{n},G^{w}_{\mathbf{n}})$ are the
reciprocal lattice vectors in the $(x,w)$-space. We substitute the above
ansatz wavefunction into the Schrödinger equation with
$\hat{H}_{\mathrm{eff}}$ and we get

$\sum_{\mathbf{n}}U^{\mathbf{k}}_{\mathbf{n}}\mathrm{e}^{\textrm{i}(\mathbf{k}+\mathbf{G}^{w}_{\mathbf{n}})\cdot\mathbf{r}_{w}}\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+v_{\textrm{ext}}(\mathbf{r})-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega}{2}\left(v-\frac{\hbar(k_{x}+G^{x}_{n})}{\sqrt{2}m_{e}}\right)^{2}-E_{\mathbf{k}}\Bigg{]}\phi^{\mathbf{k}}_{\mathbf{n}}(v)=0.$
(30)

The Hamiltonian includes a harmonic-oscillator part $\hat{H}_{v}$ which is
shifted by the momentum in the $x$-direction

$A^{k_{x}}_{n}=\frac{\hbar(k_{x}+G^{x}_{n})}{\sqrt{2}m_{e}},$ (31)

written out as

$\hat{H}_{v}=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega^{2}}{2}\left(v-A^{k_{x}}_{n}\right)^{2}.$ (32)

The eigenfunctions of this operator are the shifted Hermite functions
$\phi_{j}(v-A^{k_{x}}_{n})$ with eigenenergies

$\displaystyle\mathcal{E}_{j}$
$\displaystyle=\hbar\Omega\left(j+\frac{1}{2}\right)=\hbar\sqrt{\omega^{2}_{p}+\omega^{2}_{c}}\left(j+\frac{1}{2}\right)$
(33)
$\displaystyle=\hbar\omega_{c}\sqrt{1+\eta^{2}}\left(j+\frac{1}{2}\right)$

which are degenerate with respect to the momentum in $x$-direction. The
eigenstates $\phi_{j}\left(v-A^{k_{x}}_{n}\right)$, as it was shown in Ref.
Rokaj  _et al._ (2019), correspond to Landau polaritons Keller  _et al._
(2020) and have many structural similarities to the well-known Landau levels
Landau and Lifshitz (1997). Further, the Landau polariton energy levels can be
understood as Landau levels modified by the coupling to the cavity

$\eta=\frac{\omega_{p}}{\omega_{c}},$ (34)

which in our setting is defined as the ratio between the two fundamental
scales in the light-matter coupled system, namely the diamagnetic frequency
$\omega_{p}$ and the cyclotron frequency $\omega_{c}$. For $\eta=0$ the Landau
polariton energy levels reduce to the standard Landau levels Landau and
Lifshitz (1997).

We will now make use of the Landau polariton eigenfunctions by expanding
$\phi^{\mathbf{k}}_{\mathbf{n}}(v)$ in terms of this basis. Then, the
polaritonic Bloch ansatz takes the form

$\Psi_{\mathbf{k}}(\mathbf{r}_{w},v)=\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}_{w}}\sum_{\mathbf{n},j}U^{\mathbf{k}}_{\mathbf{n},j}\mathrm{e}^{\mathrm{i}\mathbf{G}^{w}_{\mathbf{n}}\cdot\mathbf{r}_{w}}\phi_{j}(v-A^{k_{x}}_{n}).$
(35)

With respect to our previous work Rokaj  _et al._ (2019) here we use shifted
Hermite functions $\phi_{j}(v-A^{k_{x}}_{n})$ for the construction of our
ansatz, instead of unshifted ones. The shifted basis helps to incorporate the
degeneracy of the Landau polaritons with respect to $k_{x}$.

Substituting the polaritonic Bloch ansatz again into our Schrödinger equation
and making use of the Fourier expansion of the external potential given in Eq.
(25) we have

$\displaystyle\sum_{\mathbf{n},j}U^{\mathbf{k}}_{\mathbf{n},j}\mathrm{e}^{\textrm{i}\mathbf{G}^{w}_{\mathbf{n}}\cdot\mathbf{r}_{w}}\phi_{j}\left(v-A^{k_{x}}_{n}\right)\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+\mathcal{E}_{j}-E_{\mathbf{k}}\Bigg{]}$
$\displaystyle+\sum_{\mathbf{n},\mathbf{n}^{\prime},j}V_{\mathbf{n}^{\prime}}U^{\mathbf{k}}_{\mathbf{n},j}\mathrm{e}^{\textrm{i}\mathbf{G}^{w}_{\mathbf{n}+\mathbf{n}^{\prime}}\cdot\mathbf{r}_{w}}\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}^{\prime}}v}\phi_{j}(v-A^{k_{x}}_{n})=0.$
(36)

To eliminate the plane waves depending on $\mathbf{r}_{w}$ we multiply the
above expression by
$\mathrm{e}^{-\textrm{i}\mathbf{G}^{w}_{\mathbf{q}}\cdot\mathbf{r}_{w}}$ and
integrate over $\mathbf{r}_{w}$,

$\displaystyle\sum_{j}U^{\mathbf{k}}_{\mathbf{n},j}\phi_{j}\left(v-A^{k_{x}}_{n}\right)\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+\mathcal{E}_{j}-E_{\mathbf{k}}\Bigg{]}$
$\displaystyle+\sum_{\mathbf{n}^{\prime},j}V_{\mathbf{n}-\mathbf{n}^{\prime}}U^{\mathbf{k}}_{\mathbf{n}^{\prime},j}\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}\phi_{j}(v-A^{k_{x}}_{n^{\prime}})=0.$
(37)

We note that after the integration we exchanged the index $\mathbf{q}$ with
$\mathbf{n}$ again. Next, we apply the bra222We note that the standard bra and
ket notation does not include the coordinate. Here, we kept the coordinate as
it will be convenient to perform shift transformations on these states to
obtain the matrix representation of the displacement operators. See Appendix A
for this. $\langle\phi_{i}(v-A^{k_{x}}_{n})|$ from the left,

$\displaystyle
0=U^{\mathbf{k}}_{\mathbf{n},i}\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+\mathcal{E}_{i}-E_{\mathbf{k}}\Bigg{]}+$
(38)
$\displaystyle\sum_{\mathbf{n}^{\prime},j}V_{\mathbf{n}-\mathbf{n}^{\prime}}U^{\mathbf{k}}_{\mathbf{n}^{\prime},j}\;\langle\phi_{i}(v-A^{k_{x}}_{n})|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}|\phi_{j}(v-A^{k_{x}}_{n^{\prime}})\rangle.$

The only thing left to be computed in order to obtain our QED-Bloch central
equation are the matrix elements

$\langle\phi_{i}(v-A^{k_{x}}_{n})|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}|\phi_{j}(v-A^{k_{x}}_{n^{\prime}})\rangle.$
(39)

This can be performed with the use of displacement operators and their
corresponding algebra Cahill and Glauber (1969). We present this derivation in
Appendix A and we find

$\displaystyle\langle\phi_{i}(v-A^{k_{x}}_{n})|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}|\phi_{j}(v-A^{k_{x}}_{n^{\prime}})\rangle$
(40)
$\displaystyle=\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle$

where $\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})$ is a displacement
operator with the shift given by

$\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}=-\sqrt{\frac{\mu\Omega}{2\hbar}}A^{0}_{n-n^{\prime}}-\textrm{i}\sqrt{\frac{\hbar}{2\mu\Omega}}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}.$
(41)

Further, the matrix representation of this displacement operator in the basis
$\\{\phi_{i}\\}$ is Cahill and Glauber (1969)

$\displaystyle\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=\sqrt{\frac{j!}{i!}}\alpha^{i-j}_{\mathbf{n}-\mathbf{n}^{\prime}}\mathrm{e}^{-\frac{|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2}}{2}}L^{(i-j)}_{j}(|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2})$

where $i\geq j$ and
$L^{(i-j)}_{j}(|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2})$ are the
associated Laguerre polynomials. We note that for $j>i$ one needs to take

$\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=(-1)^{j-i}\langle\phi_{j}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{i}\rangle^{*}$
(43)

because $\hat{D}^{\dagger}(\alpha)=\hat{D}(-\alpha)$. Finally, substituting
Eqs. (40) and (III.3) for the matrix representation of the displacement
operator into Eq. (38), we obtain the QED-Bloch central equation

$U^{\mathbf{k}}_{\mathbf{n},i}\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+\mathcal{E}_{i}-E_{\mathbf{k}}\Bigg{]}+\sum_{\mathbf{n}^{\prime},j}V_{\mathbf{n}-\mathbf{n}^{\prime}}U^{\mathbf{k}}_{\mathbf{n}^{\prime},j}\;\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=0.$
(44)

The equation above is the main result of QED-Bloch theory. The QED-Bloch
central equation provides a unified framework for the description of periodic
materials in the presence of homogeneous magnetic fields coupled to a
quantized electromagnetic field. It is important to mention that the QED-Bloch
central equation is also applicable in the case where there is no quantized
field, i.e., where the frequency $\omega_{p}$ is equal to zero. This implies
that QED-Bloch theory and the central equation we derived can be used also for
the description of periodic materials solely under the influence of the
homogeneous magnetic field. The semi-classical limit of our central equation,
where the quantized field is equal to zero, is performed in detail in Appendix
B. In addition, we will also discuss this limit in the context of the Harper
equation Harper (1955) and the Hofstadter butterfly Hofstadter (1976) in
Section V.

##  IV Polaritonic Hofstadter Butterfly

As a first application of our QED-Bloch theory we consider the case where we
have a 2D periodic system perpendicular to a homogeneous external magnetic
field and coupled to the quantized mode originating from a cavity as depicted
in Fig. 1. It is important to mention that in such a cavity setting the
diamagnetic frequency $\omega_{p}$ can be defined in terms of the 2D electron
density $n_{\textrm{2D}}=N/S$, where $S$ is the area of the 2D material, and
the fundamental cavity frequency $\omega_{\textrm{cav}}=\pi c/L_{z}$ as Rokaj
_et al._ (2019, 2022); Eckhardt  _et al._ (2021)

$\omega_{p}=\sqrt{\frac{e^{2}n_{e}}{m_{e}\epsilon_{0}}}=\sqrt{\frac{e^{2}n_{\textrm{2D}}\omega_{\textrm{cav}}}{m_{e}\epsilon_{0}\pi
c}}.$ (45)

###  IV.1 Polaritonic Butterfly in Square Lattice & Polaritonic Harper
Equation

(a) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=0.1$.

(b) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=0.2$.

(c) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=1$.

(d) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=1$.

Figure 2: Energy spectra as a function of the light-matter coupling $\eta$ for
the square cosine potential and for different values of the relative magnetic
flux $\Phi/\Phi_{0}$. The inset in Fig. 2(a) displays the spatial profile of
the applied square cosine lattice potential.

In this subsection, as a first example, we will look into a 2D square cosine
potential. For a square lattice potential the angle between the lattice
vectors is $\theta=\pi/2$ and the two lattice constants are equal
$a_{1}=a_{2}=a$. We emphasize that in order to achieve substantial fractions
of the magnetic flux quantum $\Phi_{0}=h/e$ in the magnetic flux
$\Phi=Ba^{2}$, we have to choose the lattice constant of the potential to be
50 times larger than the typical lattice constants in standard materials,
which are of the order of a few Ångström. Thus, for our lattice potential the
lattice constant is $a=50\times 3\textrm{\AA}=15\textrm{nm}$. Such enlarged
lattice periodicities can be achieved with Moiré materials Andrei  _et al._
(2021); Kennes  _et al._ (2021), and it is within such setups that the
experimental demonstration of the Hofstadter butterfly has been achieved Dean
_et al._ (2013); Wang  _et al._ (2015); Forsythe  _et al._ (2018).

For the square cosine potential the only non-zero Fourier components are
$V_{\pm 1,0}=V_{0,\pm 1}=V_{0}$, where $V_{0}$ determines the strength of the
potential that we choose to be $V_{0}=3\textrm{eV}$. Consequently, our 2D
square lattice potential takes the form depicted as an inset in Fig. 2(a). For
the reciprocal lattice vectors of the square cosine potential it holds
$G^{v}_{\mathbf{n}}=m_{p}G^{y}_{m}/\sqrt{2}M\omega_{c}$ and
$G^{w}_{\mathbf{n}}=G^{y}_{m}/\sqrt{2}\omega_{c}$, see for this Eqs. (16) and
(25). To further simplify our considerations, we also consider the case where
only the lowest Landau polariton $\phi_{0}(v-A^{k_{x}}_{n})$ is occupied. As
we will see later, restricting ourselves to the lowest level will help us to
connect to the well-known Harper equation Harper (1955) and the Hofstadter
butterfly Hofstadter (1976). Under these assumptions, the QED-Bloch central
equation takes the simplified form

$\displaystyle
U^{k_{x},k_{w}}_{n,m}\Bigg{[}\frac{\hbar^{2}\left(k_{w}+\frac{G^{y}_{m}}{\sqrt{2}\omega_{c}}\right)^{2}}{2M}+\mathcal{E}_{0}-E_{k_{x},k_{w}}\Bigg{]}+V_{0}U^{k_{x},k_{w}}_{n-1,m}\mathrm{e}^{-\frac{|\alpha_{1,0}|^{2}}{2}}+V_{0}U^{k_{x},k_{w}}_{n+1,m}\mathrm{e}^{-\frac{|\alpha_{-1,0}|^{2}}{2}}$
(46)
$\displaystyle+V_{0}U^{k_{x},k_{w}}_{n,m-1}\exp\left(\frac{-\textrm{i}m_{p}G^{y}_{1}A^{k_{x}}_{n}}{\sqrt{2}M\omega_{c}}\right)\mathrm{e}^{-\frac{|\alpha_{0,1}|^{2}}{2}}+V_{0}U^{k_{x},k_{w}}_{n,m+1}\exp\left(\frac{-\textrm{i}m_{p}G^{y}_{-1}A^{k_{x}}_{n}}{\sqrt{2}M\omega_{c}}\right)\mathrm{e}^{-\frac{|\alpha_{0,-1}|^{2}}{2}}=0.$

For $\theta=\pi/2$ and $a_{1}=a_{2}=a$ the $\alpha$-matrix defined in Eq. (41)
is

$\alpha_{\mathbf{n}}=-\frac{2\pi}{a}\sqrt{\frac{\hbar}{2m_{e}\Omega}}\left(n+\textrm{i}\frac{\omega_{c}}{\Omega}m\right).$
(47)

Using the above expression we find for the four components of the $\alpha$
matrix entering our equation

$\displaystyle|\alpha_{1,0}|^{2}=|\alpha_{-1,0}|^{2}=\frac{4\pi^{2}}{a^{2}}\frac{\hbar}{2m_{e}\Omega}=\frac{\pi\Phi_{0}}{\Phi\left(1+\eta^{2}\right)^{1/2}}\;\;\;\textrm{and}$
$\displaystyle|\alpha_{0,1}|^{2}=|\alpha_{0,-1}|^{2}=\frac{\pi\Phi_{0}}{\Phi\left(1+\eta^{2}\right)^{3/2}}.$
(48)

To obtain the above results we used the definition for $\mu$ and $\Omega$ and
the coupling constant $\eta$ given in Eqs. (21), (23) and (34) respectively.
In addition, we use the definitions for $m_{p}$, $M$, $\Omega$, and
$A^{k_{x}}_{n}$ given respectively in Eqs. (21), (23) and (31) and we find

$\frac{m_{p}}{\sqrt{2}M\omega_{c}}G^{y}_{\pm 1}A^{k_{x}}_{n}=\frac{\pm
1}{1+\eta^{2}}\frac{2\pi\Phi_{0}}{\Phi}\left(\frac{ak_{x}}{2\pi}+n\right).$
(49)

After these manipulations we obtain

$\displaystyle\left[\frac{\hbar^{2}\left(\sqrt{2}\omega_{c}k_{w}+G^{y}_{m}\right)^{2}}{2m_{e}(1+\eta^{-2})}+\mathcal{E}_{0}-E_{k_{x},k_{w}}\right]U^{k_{x},k_{w}}_{n,m}+t_{1}(\Phi,\eta)\left(U^{k_{x},k_{w}}_{n-1,m}+U^{k_{x},k_{w}}_{n+1,m}\right)$
(50)
$\displaystyle+t_{2}(\Phi,\eta)\left[U^{k_{x},k_{w}}_{n,m-1}\exp\left(\frac{-\textrm{i}2\pi\Phi_{0}}{\Phi(1+\eta^{2})}\left(\frac{ak_{x}}{2\pi}+n\right)\right)+U^{k_{x},k_{w}}_{n,m+1}\exp\left(\frac{\textrm{i}2\pi\Phi_{0}}{\Phi(1+\eta^{2})}\left(\frac{ak_{x}}{2\pi}+n\right)\right)\right]=0,$

where we defined the functions $t_{1}(\Phi,\eta)$ and $t_{2}(\Phi,\eta)$,
shown below, which play a similar role as the hopping matrix elements in a
tight-binding description.

$\displaystyle t_{1}(\Phi,\eta)$ $\displaystyle=$ $\displaystyle
V_{0}\exp\left(-\frac{\pi\Phi_{0}}{2\Phi(1+\eta^{2})^{1/2}}\right)$
$\displaystyle t_{2}(\Phi,\eta)$ $\displaystyle=$ $\displaystyle
V_{0}\exp\left(-\frac{\pi\Phi_{0}}{2\Phi(1+\eta^{2})^{3/2}}\right)$ (51)

The hopping functions above depend on the relative magnetic flux and the
light-matter coupling $\eta$. Equation (50) is a polaritonic extension of the
Harper equation. This can be understood from the fact that in the limit
$\eta\rightarrow 0$ the polaritonic Harper equation (50) reduces to the
standard Harper equation Harper (1955). We will see how this limit can be
performed and discuss this important point in detail in Section V.

However, the polaritonic Harper equation (50) has several important
differences to the standard Harper equation. First of all, Eq. (50) does not
describe simply electrons on a lattice under the influence of a magnetic
field, but it describes Landau polaritons on a lattice. Further, there is an
additional degree of freedom $k_{w}$ corresponding to the polaritonic Bloch
wave in the $w$ direction. Most importantly, the polaritonic Harper equation
does not only depend parametrically on the relative magnetic flux
$\Phi/\Phi_{0}$, but also on the light-matter coupling constant
$\eta=\omega_{p}/\omega_{c}$. This opens the possibility of not only having a
fractal/self-similar spectrum as a function of the relative flux, but also a
fractal as a function of the light-matter coupling constant $\eta$. The
coupling constant $\eta$ can be tuned experimentally either by varying the
strength of the external magnetic field (i.e., changing the cyclotron
frequency $\omega_{c}$) or by varying the diamagnetic frequency $\omega_{p}$
via the 2D electron density and the fundamental cavity frequency, or by
shaping the cavity environment in order to achieve a smaller effective volume
Paravicini-Bagliani  _et al._ (2019); Appugliese  _et al._ (2022).

To test the existence of this polaritonic fractal, we plot the energy spectrum
of the polaritonic Harper equation (50) as a function of the light-matter
coupling $\eta$ for different values of the relative magnetic flux
$\Phi/\Phi_{0}$. We note that in all the computations performed in this and
the following section the momenta $k_{x}$ and $k_{w}$ are taken to be equal to
zero. This is done for computational simplicity and because we found that the
inclusion of different momenta throughout the Brillouin zone has very little
influence on the energy spectra and fractal patterns.

First, we start with computing the energy spectrum as a function of $\eta$ for
a relatively small magnetic flux $\Phi/\Phi_{0}=0.1$. In Fig. 2(a) the energy
spectrum consists of well separated energy levels without much overlap between
them. For small $\eta$ the gaps between the energy levels are small and as the
light-matter coupling $\eta$ increases the energies fan out and the gaps
increase without any significant pattern emerging.

Subsequently, we double the magnetic flux and show in Fig. 2(b) the respective
energy spectrum. For $\Phi/\Phi_{0}=0.2$ we see that the energies broaden and
each energy band demonstrates an internal oscillatory behavior as a function
of $\eta$. However, the energy bands are still well-separated and there is not
much overlap between them.

In Fig. 2(c) we increase even further the relative magnetic flux and we plot
the polaritonic energies for magnetic flux equal to the flux quantum,
$\Phi/\Phi_{0}=1$, and varying $\eta$. In this case we see that a self-similar
pattern emerges as a function of the light-matter coupling, which is similar
to well-known fractal pattern of the Hofstadter butterfly Hofstadter (1976).

In addition, we compute and plot in Fig. 2(d) the polaritonic energy spectrum
for $\Phi/\Phi_{0}=1$ but now for $\eta$ ranging from $0$ to $5$. With respect
to $\eta=1$ we see that on the left and right there is self-similarity but
clearly the pattern is different on the two sides of the plot.

From these computations of the energy spectrum for different magnetic fluxes
and over different regimes of light-matter interaction we conclude that for 2D
periodic materials strongly coupled to the quantized cavity field and placed
perpendicular to a homogeneous magnetic field there is not only a fractal
spectrum emerging as a function of the magnetic flux but there is also a novel
fractal pattern showing up as a function of the light-matter coupling $\eta$.
This implies that fractal structures do not only appear due to the magnetic
field but also due to the quantized cavity field and the interaction of the
Landau polariton states with the periodic potential of the material. Thus,
what we have presented here introduces the novel concept of polaritonic
fractals or fractal polaritons. To the best of our knowledge such a phenomenon
has not been reported before.

###  IV.2 Cavity Fluctuations in the Hofstadter Butterfly

(a) Energy spectrum as a function of the relative magnetic flux
$\Phi/\Phi_{0}$ for $\omega_{p}=10^{-3}\textrm{THz}$ (weak coupling).

(b) Energy spectrum as a function of the relative magnetic flux
$\Phi/\Phi_{0}$ for $\omega_{p}=10^{-1}\textrm{THz}$ (intermediate coupling).

(c) Energy spectrum as a function of the relative magnetic flux
$\Phi/\Phi_{0}$ for $\omega_{p}=1\textrm{THz}$ (strong coupling).

Figure 3: Energy spectra for the square cosine potential as a function of the
relative magnetic flux $\Phi/\Phi_{0}$ for different values of the diamagnetic
frequency $\omega_{p}$.

So far, we computed the energy spectrum for a fixed value of the relative
magnetic flux $\Phi/\Phi_{0}$ as a function of the light-matter coupling
$\eta$ and we demonstrated that a fractal polaritonic butterfly spectrum
arises. Now we will do the opposite, namely fix the value for the strength of
the quantized photon field determined by the diamagnetic frequency
$\omega_{p}$ (see Eq. (11)) and then plot the energy spectrum as a function of
the relative magnetic flux $\Phi/\Phi_{0}$, as it was done by Hofstadter for
the purely electronic problem Hofstadter (1976). For our computations we will
use again the square cosine potential, whose energy spectrum is given by the
polaritonic Harper equation 50.

First, we compute the energy spectrum for a relatively small value of
$\omega_{p}=10^{-3}\textrm{THz}$, which is a thousand times smaller than the
standard terahertz cavity modes used in the setting of Landau polaritons
Paravicini-Bagliani  _et al._ (2019); Rokaj  _et al._ (2019); Li  _et al._
(2018). In Fig. 3(a) we show the energy spectrum as a function of the relative
magnetic flux. In the low flux regime we see the energy of the lowest Landau
level without any splitting. As the flux increases the energy level broadens
and mini-gaps start showing up in the region $0.25\leq\Phi/\Phi_{0}\leq 0.5$.
For $\Phi/\Phi_{0}>0.5$ the gaps become much larger and the fractal pattern of
the Hofstadter butterfly emerges. Due to the small value of $\omega_{p}$ no
modification of the butterfly spectrum is visible and we recover the standard
butterfly pattern. Only at $\Phi/\Phi_{0}\approx 0$ we see a little spike
deviating from the linear Landau-level dispersion, which is due to the fact
that for small fluxes the diamagnetic frequency $\omega_{p}$ is larger than
the cyclotron frequency $\omega_{c}$.

Next, in Fig. 3(b) we increase the diamagnetic frequency by two orders of
magnitude to $\omega_{p}=10^{-1}\textrm{THz}$. In this case, around
$\Phi/\Phi_{0}=0$ the energies of the system spread over a wide range and the
deviation from the linear dispersion persists over a larger region of magnetic
fluxes because the diamagnetic frequency is comparable to $\omega_{c}$ in a
larger region. However, as the magnetic flux increases, we get the linear
dispersion of the Landau level which then evolves into the Hofstadter
butterfly spectrum just like in Fig. 3(a).

Finally, we choose the diamagnetic frequency to be of the order of one
terahertz, $\omega_{p}=1\textrm{THz}$. As we showed in Rokaj  _et al._ (2019),
such a value for $\omega_{p}$ is within experimental reach in Landau polariton
platforms Paravicini-Bagliani  _et al._ (2019); Scalari  _et al._ (2012); Li
_et al._ (2018) where the fundamental cavity frequency $\omega_{\textrm{cav}}$
is also in the terahertz and the 2D electron densities are of the order
$n_{\textrm{2D}}\sim 10^{12}\textrm{cm}^{-2}$ . In Fig. 3(c) we plot the
energy spectrum as a function of the relative magnetic flux. In this case we
see that in the region $0\leq\Phi/\Phi_{0}\leq 0.5$ the linear dispersion of
the Landau level no longer shows up. In this region the energy spectrum has
been completely deformed due to the formation of Landau polariton states as a
consequence of the large value of $\omega_{p}$ Rokaj  _et al._ (2019). The
energy spectrum in this regime consists of a set of energy levels with finite
width which spread like tentacles over the energies from $-6\textrm{eV}$ to
$+6\textrm{eV}$. Then, as the magnetic flux increases the energy levels merge
together and recombine to form the fractal spectrum of the Hofstadter
butterfly. Figure 3(c) shows clearly that the fractal spectrum of the
Hofstadter butterfly gets strongly modified by the interaction with the cavity
photons. This phenomenon is a novel prediction of QED-Bloch theory.

###  IV.3 Polaritonic Butterfly in Hexagonal Lattice

(a) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=0.2$.

(b) Energy spectrum as a function of the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$ for magnetic flux ratio $\Phi/\Phi_{0}=1$.

Figure 4: Energy spectra as a function of the light-matter coupling $\eta$ for
the hexagonal cosine potential for different values of the relative magnetic
flux $\Phi/\Phi_{0}$. The inset displays the spatial profile of the applied
hexagonal cosine lattice potential.

In addition to the results that we presented for the square-lattice cosine
potential we now also apply the QED-Bloch central equation (44) for the case
of a hexagonal cosine potential (see the inset in Fig. 4(a)). We note that in
the case of the hexagonal lattice the lattice constants $a_{1}$ and $a_{2}$
are equal and we choose them to be the same with the square-lattice case,
$a_{1}=a_{2}=15\textrm{nm}$. Further, the parameter $V_{0}$ which controls the
strength of the lattice potential is the same as with the square-cosine
potential before, $V_{0}=3\textrm{eV}$.

First, we compute the energy spectrum of the combined electron-photon system
as a function of $\eta$, for a relatively small value of the relative magnetic
flux $\Phi/\Phi_{0}=0.2$. In Fig. 4(a) we plot the energy spectrum and we see
that for small values of $\eta$ the energies levels are well concentrated but
as the light-matter coupling increases the energy levels fan out (as in Fig.
2(a) for the square lattice) and the energy gaps increase. However, no
significant pattern emerges due to the small value of the magnetic flux.

Subsequently, we increase the value for the relative magnetic flux to
$\Phi/\Phi_{0}=1$. In Fig. 4(b) we plot the energy spectrum as a function of
the light-matter coupling $\eta$ and we clearly see that a self-similar
pattern shows up. In comparison to the energy spectrum for the square lattice
depicted in Fig. 2(c), we notice that for the hexagonal lattice the energy
spectrum it is not vertically symmetric. This is typical for hexagonal
lattices and was also reported for the Hofstadter butterfly in the hexagonal
lattice case Claro and Wannier (1979).

###  IV.4 Connection to Floquet Engineering

In the recent years, another topic of interest in the field of 2D materials in
homogeneous magnetic fields has been the Floquet driving of the Hofstadter
butterfly Du  _et al._ (2018); Zhao  _et al._ (2020); Kooi  _et al._ (2018).
Also the connection between Floquet driving and cavity engineering has
attracted considerable attention and has been explored extensively Sentef  _et
al._ (2020); Schäfer  _et al._ (2018); Hübener  _et al._ (2021); Eckhardt  _et
al._ (2021); Lloyd-Hughes  _et al._ (2021). In what follows, we will try to
compare some of the results coming from our QED-Bloch theory to the ones
obtained with Floquet driving.

A generic feature of Floquet theory is that the Floquet field produces
photonic copies, above and below the original bare electronic bands Hübener
_et al._ (2021); Lloyd-Hughes  _et al._ (2021). The photonic Floquet copies
show up because the Floquet field generates states in the time-domain which
dress the electronic states. Then, the quasienergy spectrum can become dense,
as it was shown in Ref. Hone  _et al._ (1997).

Such photonic copies are also a feature of QED. However, there is a
fundamental difference between non-relativistic QED and Floquet theory. The
Floquet Hamiltonian is unbounded from below and produces states which have an
arbitrarily negative energy. On the other hand, the Pauli-Fierz Hamiltonian of
non-relativistic QED is bounded from below and the photons do not produce
states with arbitrarily negative energies Spohn (2004). This means that
electron-photon systems in QED have a stable ground state Spohn (2004); Rokaj
_et al._ (2018); Schäfer  _et al._ (2019). This is an advantage of QED
compared to semi-classical Floquet theory. Thus, in non-relativistic QED one
cannot expect to obtain photonic copies of the energy bands below the original
bare electronic bands, yet above such copies occur.

To demonstrate these two basic features, (i) the photonic copies and (ii)
denseness of the spectrum, we plot in Fig.5(a) the energy spectrum for the
square lattice, in the small flux regime and for small diamagnetic frequency
$\omega_{p}=10^{-3}$THz with four Landau polariton states included. As it is
shown in Fig. 5(a), we have indeed the higher states introducing copies of the
lowest polariton band. The copies here are well-separated without much
interference, and the energy separation is given by the Landau polariton
excitation $\hbar\Omega$. It is important to note that we do not work with a
tensor-product basis between the electrons and photons, like in Floquet
theory, and as consequence the higher states are polaritonic and the energy
separation is the one of the Landau polariton. In addition, we clearly see
that the energy spectrum in the small flux regime becomes dense, as expected
Hone  _et al._ (1997). Thus, our QED-Bloch theory recovers these two basic
features of Floquet theory. It is important to note that for $\omega_{p}\ll
10^{-3}$THz the states between Landau polaritons which produce the dense
spectrum do not show up, because we practically have no photonic contribution,
and the spectrum is no longer dense, as expected in the case of no photons.

Furthermore, in the studies of the Hofstadter butterfly under Floquet driving
Du  _et al._ (2018); Zhao  _et al._ (2020); Kooi  _et al._ (2018), another
interesting phenomenon has been identified: In the regime where the frequency
of the driving field is small, such that the energy of the photon is not much
larger than the band width, the Floquet copies overlap with the original bare
electronic Hofstadter butterfly and modify the spectrum substantially Du  _et
al._ (2018); Zhao  _et al._ (2020); Kooi  _et al._ (2018).

This regime, from the perspective of cavity QED can be understood as a strong-
coupling scenario where the electronic and the photonic states mix strongly.
To compare our QED-Bloch theory to the Floquet results we choose the
diamagnetic frequency to be large $\omega_{p}=1$THz and we plot the energy
spectrum as a function of $\Phi/\Phi_{0}$ for the square cosine potential with
two Landau polariton states in Fig. 5(b). As previously, the photon field
introduces a copy of the Hofstadter butterfly, but in this case the copies
overlap and interfere. This leads to a substantial modification of the result
we obtained in the case where only the lowest Landau polariton was taken into
account (see Fig. 3(c)). As it is shown in Fig. 5(b), the two butterflies
merge together and form a new fractal pattern. This result is in qualitative
agreement with the Floquet driving of the Hofstadter butterfly. To conclude,
the results that we presented in this subsection show that our QED-Bloch
theory recovers the common features of Floquet theory and of the Floquet
driving of the Hofstadter butterfly.

(a) Energy spectrum as a function of $\Phi/\Phi_{0}$ for a small diamagnetic
frequency $\omega_{p}=10^{-3}$THz with four Landau polariton states.

(b) Energy spectrum as a function of $\Phi/\Phi_{0}$ for a large diamagnetic
frequency $\omega_{p}=1$THz with two Landau polariton states.

Figure 5: Energy spectra for the square cosine potential as a function of the
relative magnetic flux $\Phi/\Phi_{0}$ with higher Landau polariton states,
for two different values of the diamagnetic frequency $\omega_{p}$ which
correspond to weak (left) and strong (right) coupling to the cavity photons.

##  V Harper Equation & Minimal Coupling-Peierls Phase Duality

Our aim in this section is to connect the polaritonic Harper equation that we
derived in Eq. (50) to the well-known Harper equation Harper (1955) and the
fractal spectrum of the Hofstadter butterfly Hofstadter (1976). Both, the
Harper equation and the butterfly spectrum were derived from a one-band tight-
binding model with next-neighbor hopping in real space on a square lattice,
with the electrons coupled to the magnetic field via the Peierls substitution.
We note that using the Peierls substitution relies on the assumption that the
external magnetic field is weak enough such that no mixing between different
bands occurs Harper (1955); Luttinger (1951).

The polaritonic Harper equation we derived is for a square lattice potential
as well, and we chose the cosine potential because it introduces only next-
neighbor hopping in Fourier space, as the only non-zero Fourier components of
the potential are $V_{\pm 1,0}=V_{0,\pm 1}=V_{0}$. With the extra assumption
of only the lowest Landau polariton being occupied we obtained Eq. (50). To
connect now the polaritonic Harper equation to the original Harper equation we
need to take the limit of the light-matter coupling to zero, $\eta\rightarrow
0$. Then the kinetic term which depends on the polaritonic momentum $k_{w}$
vanishes because $\eta^{-2}$ goes to infinity and the equation becomes
completely independent of $k_{w}$ and as a consequence also the Fourier index
$m$ can be dropped. In addition, for $\eta=0$, the hopping functions
$t_{1}(\Phi,\eta)$ and $t_{2}(\Phi,\eta)$ defined in Eq. (IV.1) become equal,

$t_{1}(\Phi,\eta)=t_{2}(\Phi,\eta)=V_{0}\mathrm{e}^{-\frac{\pi\Phi_{0}}{2\Phi}}\equiv
t(\Phi).$ (52)

Also, in the described limit the energy of the lowest Landau polariton
$\mathcal{E}_{0}=\hbar\Omega/2$ for zero light-matter coupling becomes equal
to the energy of the lowest Landau level $\mathcal{E}_{0}=\hbar\omega_{c}/2$.
This means that for $\eta\rightarrow 0$ the polaritonic Harper equation does
not describe Landau polaritons on a lattice but actually Landau levels on a
lattice and it takes the form

$U^{k_{x}}_{n}\left(\frac{\hbar\omega_{c}}{2}-E_{k_{x}}\right)+t(\Phi)\left[U^{k_{x}}_{n-1}+U^{k_{x}}_{n+1}+2U^{k_{x}}_{n}\cos\left(\frac{2\pi\Phi_{0}}{\Phi}\left(\frac{ak_{x}}{2\pi}+n\right)\right)\right]=0.$
(53)

In the work of Hofstadter Hofstadter (1976) the fractal spectrum appears not
for the energy itself, $E$, but for the unitless scaled energy,
$\mathcal{E}=E/t$, where $E$ is divided by the constant hopping parameter $t$.
Of course this does not make a difference within the tight-binding model
because the hopping parameter $t$ is a constant anyway. On the contrary, for
the minimal-coupling Hamiltonian the magnetic field is part of the covariant
(physical) momentum of the electron. Thus, the kinetic energy of the electrons
naturally depends on the magnetic field and as a consequence the hopping
(which represents the kinetic energy in the tight-binding approach) should be
a function of the magnetic field as well. In our setting we have the flux-
dependent hopping parameter $t(\Phi)$ defined in Eq. (52) and as a consequence
the dimensionless scaled energies can be defined as

$\mathcal{E}_{k_{x}}=\frac{1}{t(\Phi)}\left(E_{k_{x}}-\frac{\hbar\omega_{c}}{2}\right)=\frac{\mathrm{e}^{\frac{\pi\Phi_{0}}{2\Phi}}}{V_{0}}\left(E_{k_{x}}-\frac{\hbar\omega_{c}}{2}\right).$
(54)

With this definition, we find the following equation for the scaled,
dimensionless energies of the system,

$\mathcal{E}_{k_{x}}U^{k_{x}}_{n}=U^{k_{x}}_{n-1}+U^{k_{x}}_{n+1}+2U^{k_{x}}_{n}\cos\left(\frac{2\pi\Phi_{0}}{\Phi}\left(\frac{ak_{x}}{2\pi}+n\right)\right).$
(55)

The equation above is precisely the usual Harper equation Harper (1955) and
plotting the eigenenergies of this equation we obtain the fractal spectrum of
the Hofstadter butterfly Hofstadter (1976), which is depicted in Fig 6(a).
However, there is one important difference: In the original Harper equation
the energy spectrum is a function of the relative magnetic flux
$\Phi/\Phi_{0}$ and the corresponding butterfly spectrum also appears as a
function of the relative flux $\Phi/\Phi_{0}$. In our case the energies are a
function of the reciprocal relative flux $\Phi_{0}/\Phi$ and the fractal
spectrum appears with respect to the reciprocal flux $\Phi_{0}/\Phi$. This
fact, that starting from the minimal-coupling Hamiltonian one obtains the
Hofstadter butterfly as a function of the reciprocal flux $\Phi_{0}/\Phi$
comprises a fundamental difference to tight-binding models with the Peierls
phase, as has been shown and discussed in several publications Pfannkuche and
Gerhardts (1992); Langbein (1969); Rauh (1975); Petschel and Geisel (1993).

Already Langbein Langbein (1969) discussed and compared the two approaches
while it was noted by Wannier Wannier (1978) that it is highly surprising for
the energy spectrum to be periodic either in the relative magnetic flux
$\Phi/\Phi_{0}$ or in the reciprocal magnetic flux $\Phi_{0}/\Phi$, because
the minimal-coupling Hamiltonian has a linear and a quadratic dependence with
respect to the magnetic field. Both periodic behaviors are in fact an
artefact. As it was pointed out by Wannier Wannier (1978) this periodicity
cannot be physical, and the tight-binding model with the Peierls phase cannot
be trusted for magnetic fields beyond one flux quantum, because for such
strong magnetic fields the original assumption that each band can be treated
separately and that there is no mixing between them will be invalid.

In the case of the minimal-coupling Hamiltonian the exact periodicity is due
to the redefinition (the scaling) of the energies in Eq. (54). This
redefinition of the energy spectrum cuts out the linear dependence of the
energy of the lowest Landau level on $\omega_{c}=eB/m_{e}$ and the exponential
increase due to the flux dependent hopping parameter $t(\Phi)$. To understand
what the actual dependence of the energy spectrum is, we plot the unscaled,
dimensionful energies of our system given by Eq. (53) in Fig. 6(b).

(a) Scaled dimensionless energy spectrum obtained from the Harper equation
(55) as a function of the reciprocal magnetic flux $\Phi_{0}/\Phi$.

(b) Energy spectrum as function of the relative magnetic flux $\Phi/\Phi_{0}$
as given by Eq. (53).

Figure 6: Comparison between the scaled and the unscaled spectrum for the
square cosine potential. The potential strength is $V_{0}=3$eV and the lattice
constant $a=2$Å. We use a small lattice constant to make the linear dispersion
of the Landau level in Fig. 6(b) clearly visible.

We see that for small fluxes we have the linear dispersion coming from the
energy of the lowest Landau level. As the magnetic field increases, the Landau
level starts to split and gaps show up. Then, for $\Phi/\Phi_{0}>1/2$ the
fractal nature of the spectrum shows up and the Hofstadter butterfly becomes
clearly visible. The Hofstadter butterfly however is not periodic, but it
actually spreads out due to the flux-dependent hopping $t(\Phi)$. Figure 6(b)
reconciles the two fundamental properties of the minimal-coupling Hamiltonian:
(i) the energy has to increase as a function of the magnetic field and (ii)
due to the reduced symmetry manifest in the magnetic translation group a
splitting of the energy bands for every fractional value of the relative
magnetic flux $\Phi/\Phi_{0}=p/q$ needs to occur, which subsequently leads to
the formation of the fractal Brown (1964); Zak (1964a, b).

Dual Descriptions.—The question that finally arises is: Why do the energy
spectra of the minimal-coupling Hamiltonian and the tight-binding model with
the Peierls phase differ in such a fundamental way?

In many cases, these two descriptions for electrons in periodic structures
coupled to electromagnetic fields are indeed equivalent descriptions and match
at least to some certain accuracy. Typically, this is true for electromagnetic
fields slowly varying within the unit cell of the solid. However, the problem
with our particular system is that although the magnetic field is constant,
the vector potential that actually couples to the electrons is linear in space
$\mathbf{A}_{\textrm{ext}}=-\mathbf{e}_{x}By$. As it was pointed out by
Luttinger, the spatial variation of the vector potential makes the Peierls
substitution questionable for large magnetic fields Luttinger (1951). This was
one of the first papers deriving the tight-binding model with Peierls phase
starting from the minimal-coupling Hamiltonian. To arrive at this model,
Luttinger had to drop a term from the Hamiltonian which consequently breaks
the actual relation between the minimal-coupling Hamiltonian and the tight-
binding model with the Peierls phase.

However, this does not mean that the two descriptions are completely
disconnected. They both yield the Harper equation and the Hofstadter
butterfly, with the difference that in the one case it shows up as a function
of the magnetic flux $\Phi/\Phi_{0}$ while in the other as a function of the
reciprocal flux $\Phi_{0}/\Phi$. This means that the minimal-coupling
Hamiltonian and the tight-binding model with the Peierls phase are not
equivalent but they are actually dual. It is important to mention that this
duality holds only in the lowest Landau level and for a single band,
respectively.

This dual (or reciprocal) relation between the two models can be understood on
the one hand from the fact that in the tight-binding model we have electrons
with next-neighbor hopping on a lattice in real space, while in the minimal-
coupling Hamiltonian we have next-neighbor hopping in k-space. Further, the
dual relation manifests itself as a weak-to-strong duality with respect to the
strength of the magnetic field, in the sense that the tight-binding model with
the Peierls phase mainly describes the regime of not too large fluxes and the
butterfly appears as a function of the flux, where the Peierls substitution is
still applicable. Instead the minimal-coupling Hamiltonian, restricted in the
lowest Landau level, works best for not too small fluxes and the butterfly
shows up as a function of the inverse flux. The two approaches match in the
regime where $\Phi/\Phi_{0}\approx\Phi_{0}/\Phi$ which is the region around
one flux quantum. We note that by the inclusion of higher Landau levels also
small magnitudes of the magnetic flux can be correctly described within the
minimal-coupling framework. The duality between the two approaches and the
steps to obtain the respective Hofstadter butterflies are summarized in Fig.
7.

Figure 7: Schematic illustration of the duality between the minimal-coupling
Hamiltonian and the tight-binding model with Peierls substitution. In both
models the Hofstadter butterfly emerges but in a dual fashion. Namely, in the
one case as a function of the reciprocal flux while in the other as a function
of flux.

Having established this duality between the minimal-coupling Hamiltonian and
the tight-binding model with the Peierls phase in the semi-classical setting,
we would now like to make a connection to the tight-binding model with Peierls
phase in the case where the electrons are coupled also to the quantized cavity
field. In subsection IV.2 in Fig. 3(c) we showed how the butterfly pattern
coming from the minimal-coupling Hamiltonian gets modified due to strong
coupling to a terahertz cavity. In the semi-classical case, to connect the two
approaches we applied the scaling transformation defined in Eq. (54) to the
energies of the minimal-coupling Hamiltonian. In the light-matter setting the
analogous scaling transformation is

$\mathcal{E}_{k_{x},k_{w}}=\frac{E_{k_{x},k_{w}}-\hbar\Omega/2}{t_{1}(\Phi,\eta)+t_{2}(\Phi,\eta)},$
(56)

where the hopping functions $t_{1}(\Phi,\eta)$ and $t_{2}(\Phi,\eta)$ were
defined in Eq. (IV.1) and $\hbar\Omega/2$ is the lowest Landau polariton
energy. Using this scaling transformation on the polaritonic energies of Fig.
3(c) we obtain the scaled energy spectrum depicted in Fig. 8. The scaled
energy spectrum provides a hint on how the original Hofstadter butterfly
(coming from the tight-binding model) might look when the 2D material is
coupled also to the cavity photon-field in the terahertz regime. In the region
of small and intermediate fluxes, where the cavity field dominates, we see
that the butterfly pattern is dissolved into distinct energy bands with an
internal oscillatory behavior. As the flux becomes larger than half of a flux
quantum we see that the butterfly pattern emerges. Still at values around one
flux quantum there are modifications of the spectra and the effect of the
cavity is noticeable. There is a clear similarity between the butterfly
patterns in and outside the cavity around one, but in contrast to the
perfectly periodic spectrum (of the Hofstadter butterfly coming from the
tight-binding model) the cavity induces a significant distortion that should
be visible experimentally.

Figure 8: Scaled dimensionless energy spectrum for the square cosine lattice
potential as a function of the magnetic flux $\Phi/\Phi_{0}$ for
$\omega_{p}=1\textrm{THz}$. The unscaled physical spectrum is shown in Fig.
3(c).

##  VI Cavity Modification of the Integer Hall Conductance

QED-Bloch theory is also applicable to the case of no external potential,
which for a homogeneous magnetic field (without a quantized field) constitutes
the usual setting of non-interacting Landau levels that describe the integer
quantum Hall effect Klitzing  _et al._ (1980); Laughlin (1981).

The aim of this section is to study non-interacting Landau levels under strong
coupling to the cavity field, in the presence of an external classical
electric field (as shown in Fig. 9), and to demonstrate that the cavity field
modifies the plateaus of the Hall conductance in the integer regime. Such a
cavity modification has also been measured recently experimentally Appugliese
_et al._ (2022).

As it was shown in Ref. Rokaj  _et al._ (2019), in the case where the external
potential is zero, $v_{\textrm{ext}}(\mathbf{r})=0$, the effective Hamiltonian

$\hat{H}_{\mathrm{eff}}=\frac{1}{2m_{e}}\left(\mathrm{i}\hbar\mathbf{\nabla}+e\hat{\mathbf{\mathcal{A}}}+e\mathbf{A}_{\textrm{ext}}(\mathbf{r})\right)^{2}-\frac{\hbar\omega_{p}}{2}\frac{\partial^{2}}{\partial
u^{2}}$ (57)

is analytically solvable. As a reminder for the reader
$\mathbf{A}_{\textrm{ext}}(\mathbf{r})=-\mathbf{e}_{x}By$ is the vector
potential describing the external magnetic field and
$\hat{\mathbf{\mathcal{A}}}$ is the quantized vector potential of the cavity
defined in Eq. (11)

The energy spectrum of this system can be directly obtained from the QED-Bloch
central equation (44) by simply setting the external potential to zero,
$v_{\textrm{ext}}(\mathbf{r})=0$. Since there is no external potential, the
reciprocal lattice vectors have no role and have to be taken equal to zero.
This is done by simply setting $n=m=0$ and we get

$U^{\mathbf{k}}_{i}\left(\frac{\hbar^{2}k^{2}_{w}}{2M}+\mathcal{E}_{i}-E_{\mathbf{k}}\right)=0.$
(58)

Then, from the above equation it is clear that the eigenspectrum of the 2D
Landau levels coupled to the cavity is

$E_{\mathbf{k},i}=\frac{\hbar^{2}k^{2}_{w}}{2M}+\hbar\Omega\left(i+\frac{1}{2}\right).$
(59)

Further, the components of $U^{\mathbf{k}}_{i}$ from the QED-Bloch ansatz
defined in Eq. (35) become trivial and we obtain the full set of eigenfuctions
corresponding to the 2D Landau levels in the cavity,

$\Psi_{\mathbf{k},i}(\mathbf{r}_{w},v)=\mathrm{e}^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}_{w}}\phi_{i}\left(v-\frac{\hbar
k_{x}}{\sqrt{2}m_{e}}\right).$ (60)

This is the analytic solution for the non-interacting 2D Landau levels in the
cavity that was found in Ref. Rokaj  _et al._ (2019). The eigenfunctions above
are plane waves in the directions $x$ and $w$ because in these directions we
have translational invariance. The $v$-dependent eigenfunctions are Hermite
funcitons Griffiths (1995) with argument $v-\hbar k_{x}/\sqrt{2}m_{e}$. The
eigenfunctions in Eq. (60), since they are functions of the combined
polaritonic coordinates $w$ and $v$ defined in Eq. (19), describe
quasiparticles formed between the Landau levels and the photons, which are
known as Landau polaritons. Such Landau polariton states have been studied
theoretically Hagenmüller  _et al._ (2010) and have been also observed
experimentally Scalari  _et al._ (2012); Keller  _et al._ (2020); Li  _et al._
(2018).

Figure 9: Cartoon depiction of a 2D electron gas (material in black) confined
inside a cavity. The whole system is placed perpendicular to a classical
homogeneous magnetic field $\mathbf{B}_{\textrm{ext}}$, and an external
constant electric field $\mathbf{E}_{\textrm{ext}}$ is applied to the 2D
material.

###  VI.1 Integer Hall Conductance in the Cavity

Having briefly reviewed how the Landau polariton quasiparticles emerge in QED-
Bloch theory, our aim is now to study the implications of these polaritonic
quasiparticles for the quantum Hall effect. To compute the Hall conductance of
the Landau polaritons we need to add an external potential
$V_{E}(\mathbf{r})=e\phi_{E}(\mathbf{r})=-e\mathbf{E}_{\textrm{ext}}\cdot\mathbf{r}=-eEy$
to our effective Hamiltonian, which generates an electric field in the $y$
direction,
$\mathbf{E}_{\textrm{ext}}=-\nabla\phi_{E}(\mathbf{r})=E\mathbf{e}_{y}$. This
is what was done by Laughlin Laughlin (1981) for the description of the
integer quantum Hall effect in terms of non-interacting Landau levels. We note
that in his description the spin degrees of freedom of the electrons were
neglected and we do the same here.

Then, with the addition of the external potential $V_{E}(\mathbf{r})$, the
effective Hamiltonian for the 2D Landau levels coupled to the cavity is
$\hat{H}_{E}=\hat{H}_{\mathrm{eff}}+V_{E}(\mathbf{r})=\hat{H}_{\mathrm{eff}}-eEy$.
In terms of the polaritonic coordinates $w$ and $v$ introduced in Eq. (19)
$\hat{H}_{E}$ takes the simple form

$\displaystyle\hat{H}_{E}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
w^{2}}-\frac{eE}{\sqrt{2}\omega_{c}}w$ $\displaystyle-$
$\displaystyle\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega^{2}}{2}\left(v+\frac{\textrm{i}\hbar}{\sqrt{2}m_{e}}\frac{\partial}{\partial
x}\right)^{2}+\frac{eEm_{p}}{\sqrt{2}M\omega_{c}}v.$

The expression above can be easily derived from Eq. (22) and from the
expression for the $y$ coordinate in terms of the polaritonic coordinates $w$
and $v$,

$y=\frac{w}{\sqrt{2}\omega_{c}}-\frac{m_{p}v}{\sqrt{2}M\omega_{c}}.$ (62)

Due to translational invariance it is clear that the eigenfunctions with
respect to $x$ are plane waves of the form
$f(x)=\mathrm{e}^{\textrm{i}k_{x}x}$. Applying $\hat{H}_{E}$ on
$f(x)=\mathrm{e}^{\textrm{i}k_{x}x}$ and then dividing by $f(x)$ we obtain

$\displaystyle\hat{H}_{E}[k_{x}]$ $\displaystyle\equiv$
$\displaystyle\frac{1}{f(x)}\hat{H}_{E}f(x)=-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial
w^{2}}-\frac{eE}{\sqrt{2}\omega_{c}}w$ $\displaystyle-$
$\displaystyle\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega^{2}}{2}\left(v-\frac{\hbar
k_{x}}{\sqrt{2}m_{e}}\right)^{2}+\frac{eEm_{p}}{\sqrt{2}M\omega_{c}}v.$

The $w$-dependent part of the Hamiltonian is independent of the electronic
momentum $k_{x}$ and as a consequence it cannot give any contribution to the
induced current along the $x$ direction which we are interested in for the
computation of the Hall conductance $\sigma_{xy}$. Thus, for our purpose, we
can safely eliminate the $w$-dependent terms from the Hamiltonian and remain
with

$\hat{H}_{E}[k_{x}]=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega^{2}}{2}\left(v-\frac{\hbar
k_{x}}{\sqrt{2}m_{e}}\right)^{2}+\frac{\sqrt{2}m_{e}E}{B(1+\eta^{2})}v.$ (64)

In the previous equation we substituted the definitions for the cyclotron
frequency $\omega_{c}=eB/m_{e}$ and the mass parameters $M$ and $m_{p}$ given
by Eq. (21), and we introduced the light-matter coupling
$\eta=\omega_{p}/\omega_{c}$. Next we perform a square completion and the
Hamiltonian takes the form of a shifted harmonic oscillator

$\hat{H}_{E}[k_{x}]=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
v^{2}}+\frac{\mu\Omega^{2}}{2}\left(v-\frac{\hbar
k_{x}}{\sqrt{2}m_{e}}+\frac{E}{\sqrt{2}B(1+\eta^{2})}\right)^{2}+\frac{\hbar
k_{x}E}{B(1+\eta^{2})}-\frac{m_{e}}{2}\left(\frac{E}{B(1+\eta^{2})}\right)^{2}.$
(65)

In order to obtain the above expression we also used the definitions for $\mu$
and $\Omega$ given by Eqs. (21) and (23). The eigenfunctions of the above
Hamiltonian are Hermite functions $\phi_{n}(Z)$ which depend on the variable

$Z=v-\frac{\hbar k_{x}}{\sqrt{2}m_{e}}+\frac{E}{\sqrt{2}B(1+\eta^{2})}$ (66)

with eigenenergies

$E_{n,k_{x}}=\hbar\Omega\left(n+\frac{1}{2}\right)+\frac{\hbar
k_{x}E}{B(1+\eta^{2})}-\frac{m_{e}}{2}\left(\frac{E}{B(1+\eta^{2})}\right)^{2}.$
(67)

We note that the eigenstates $\phi_{n}(Z)$ are Landau polariton states shifted
by the external electric field.

From the energy spectrum it is clear that the degeneracy with respect to
$k_{x}$ is now lifted due to the electric field. But here the strength of the
electric field $E$ is considered to be much smaller than the strength of the
magnetic field $B$, which implies that $E/B\approx 0$. As a consequence, the
degeneracy with respect to $k_{x}$ remains effectively the same. Further,
having the expression for the energies we can straightforwardly compute the
group velocity of the Landau polariton states in the $x$ direction Aschroft
and Mermin (1976)

$v_{x}=\frac{1}{\hbar}\frac{\partial E_{n,k_{x}}}{\partial
k_{x}}=\frac{E}{B(1+\eta^{2})}.$ (68)

The total current of the system in the $x$ direction $J_{x}$ is obtained by
summing the velocity $v_{x}$ over all occupied Landau-polariton states,

$J_{x}=e\sum^{\nu-1}_{n=0}\frac{L}{2\pi}\int^{\frac{eBL}{\hbar}}_{0}v_{x}dk_{x}.$
(69)

Here we assumed that we have $\nu$ Landau-polariton states occupied, from
$n=0$ to $n=\nu-1$. Further, in order to specify the region of integration for
the electronic momentum $k_{x}$, we followed the standard procedure used in
the Landau-level setting, where the system is considered to have a finite size
in 2D with area $S=L^{2}$ and the momenta $k_{x}$ range from $0$ to
$eBL/\hbar$ Tong (2016); Peierls (1955). After integrating over $k_{x}$,
summing over $n$, and dividing by the area of the system $S=L^{2}$, we find
that the total current density $j_{x}=J_{x}/S$ induced by the electric field
is

$j_{x}=\frac{e^{2}\nu}{h(1+\eta^{2})}E.$ (70)

From the above expression we deduce that the Hall conductance $\sigma_{xy}$ of
the 2D Landau levels coupled to the cavity field is still quantized and is
given by the expression

$\sigma_{xy}=\frac{e^{2}}{h(1+\eta^{2})}\nu\;\;\;\textrm{with}\;\;\nu\in\mathbb{N}.$
(71)

The Hall conductance depends on the electron charge $e$, on Planck’s constant
$h$, on the amount of occupied Landau-polariton states $\nu$ and on the light-
matter coupling between the Landau levels and the cavity photons
$\eta=\omega_{p}/\omega_{c}$. This result makes clear that the cavity modifies
the plateaus of the Hall conductance in the integer regime. This is a
significant result as it demonstrates that the cavity field, due to its long-
range nature, can circumvent the topological protection of edge states which
carry the Hall current and modify the fundamental Hall conductance. In the
limit of the light-matter coupling going to zero, $\eta\to 0$, we recover the
standard result for the quantization of the Hall conductance
$\sigma_{xy}=e^{2}\nu/h$ as it is shown in Fig. 10. This is a beautiful
consistency check of our QED-Bloch theory. However, as the light-matter
coupling increases, the value of the quantum Hall plateaus decrease as it is
shown in Fig. 10. This decrease of the Hall plateaus can be understood as a
renormalization or screening effect due to the strong vacuum fluctuations
induced by the cavity field.

Figure 10: Depiction of the first quantum Hall plateau as a function of the
light-matter coupling constant $\eta$. For $\eta=0$ we obtain the standard
value of the Hall conductance. As the light-matter coupling increases the Hall
conductance decreases. This effect can be understood as a screening or
renormalization effect due to the cavity. The shaded area indicates the regime
in which typically experiments are performed.

To gain more understanding on the behavior of the Hall conductance inside the
cavity we would also like to plot the cavity-modified Hall conductance given
by Eq. (71) as a function of the strength of the external magnetic field $B$.
To do so, we rewrite the dimensionless light-matter coupling constant
$\eta=\omega_{p}/\omega_{c}$ as a function of the magnetic field $B$ by simply
using the definition for the cyclotron frequency, $\omega_{c}=eB/m_{e}$. Upon
this substitution $\eta$ takes the form

$\eta=\frac{m_{e}\omega_{p}}{eB}=\frac{B_{\textrm{cav}}}{B}\;\;\;\textrm{where}\;\;\;B_{\textrm{cav}}=\frac{m_{e}\omega_{p}}{e}.$
(72)

We note that the quantity $B_{\textrm{cav}}$ has magnetic-field dimensions and
it describes the magnetic field strength that corresponds to the diamagnetic
frequency $\omega_{p}$ of the cavity. For $\omega_{p}$ at the order of one
THz, as in the setting of Landau polaritons Paravicini-Bagliani  _et al._
(2019); Scalari  _et al._ (2012); Li  _et al._ (2018), the strength of the
cavity field is on the order of one Tesla.

Substituting now the expression for $\eta$ as a function of the magnetic
field, the Hall conductance takes the form

$\sigma_{xy}=\frac{e^{2}\nu}{h}\frac{1}{1+B^{2}_{\textrm{cav}}/B^{2}}.$ (73)

In Fig. 11 we plot the Hall conductance as a function of the relative magnetic
field $B/B_{\textrm{cav}}$. Figure 11 shows that in the regime where the
strength of the cavity field $B_{\textrm{cav}}$ is larger than the external
magnetic field $B$, the cavity-modified Hall conductance deviates strongly
from the value of the Hall plateau outside the cavity. As the strength of the
external magnetic field increases and becomes larger than the cavity field,
$B>B_{\textrm{cav}}$ (this is the common regime in experiments), the Hall
conductance approaches the value outside the cavity and demonstrates the
standard Hall plateau. This result provides an alternative understanding, from
a QED point of view, of why strong magnetic fields are required for the
formation of an exact Hall plateau. Namely, that the external magnetic field
should reach such a value that it becomes stronger and dominates the vacuum
fluctuations of the quantized electromagnetic field.

Figure 11: First quantized Hall plateau as a function of the ratio between the
strength of the external magnetic field and the strength of the cavity
magnetic field $B/B_{\textrm{cav}}$. In the region where the cavity strength
is larger, the Hall conductance deviates strongly from its standard value
$e^{2}/h$. As the external magnetic field increases and becomes larger than
the strength of the cavity we obtain the standard quantized Hall plateau. The
shaded area indicates the regime in which typically experiments are performed.

Finally, we note that in the case where many modes contribute to the light-
matter coupling, the strength of the cavity magnetic field gets enhanced and
renormalizes the light-matter coupling and the Hall conductance as we show in
Appendix C.

##  VII Summary & Outlook

In this article we studied two-dimensional non-interacting quantum Hall
systems strongly coupled to the quantized photon field of a cavity in the
optical limit (or dipole approximation). To do so, we employed the recently
proposed QED-Bloch theory in which the broken translational symmetry due to an
external homogeneous magnetic field is restored by taking into account the
quantum fluctuations of the photon field Rokaj  _et al._ (2019) and which thus
provides a first-principles framework for the description of periodic
materials in homogeneous magnetic fields and strongly coupled to the photon
field. We reviewed the basic steps in the construction of QED-Bloch theory and
we introduced the single-particle effective Hamiltonian which has been
successful in the description of Landau polaritons Rokaj  _et al._ (2019). In
this general framework we constructed a QED-generalization of the Bloch ansatz
(the QED-Bloch ansatz) and we derived the respective QED-Bloch central
equation.

Subsequently, we applied the effective Hamiltonian to the study of two
particular quantum Hall systems under strong coupling to the quantized cavity
field: (i) the Hofstadter butterfly and (ii) the integer quantum Hall effect
in the regime of non-interacting Landau levels.

Polaritonic Hofstadter Butterfly.—In the case of a non-interacting 2D periodic
material under cavity confinement and in the presence of a homogeneous
magnetic field we found that for the energy spectrum of the system as a
function of the dimensionless light-matter coupling constant a self-similar
pattern emerges which we call the polaritonic Hofstadter butterfly (see Fig.
2(d)). This polaritonic fractal is an extension of the standard Hofstadter
butterfly Hofstadter (1976) in cavity QED. With the advent of Moiré materials
the Hofstadter butterfly has become now experimentally accessible through
magnetotransport measurements of the Wannier diagram (integrated density of
states) Wannier (1978); Dean  _et al._ (2013); Wang  _et al._ (2015); Barrier
_et al._ (2020); Forsythe  _et al._ (2018). We believe that our prediction of
the existence of fractal polaritonic spectra due to strong light-matter
coupling can be observed in such Moiré systems under cavity confinement via
transport measurements. Further, the prediction of the polaritonic butterfly
opens a new avenue for the exploration of fractal physics in the field of
cavity QED. This novel polaritonic fractal could potentially be interesting
also from a mathematical point of view, like the original Hofstadter butterfly
Hofstadter (1976); Avila and Jitomirskaya (2006), and might provide new
connections between mathematics and physics Ikeda (2018).

Cavity Engineering of the Hofstadter Butterfly.—In addition, we computed the
energy spectrum for a periodic material as a function of the relative magnetic
flux and we found that for a terahertz cavity Scalari  _et al._ (2012);
Paravicini-Bagliani  _et al._ (2019); Li  _et al._ (2018) the standard
Hofstadter butterfly gets modified due to the strong vacuum fluctuations of
the photon field (see Fig. 3(c)). This phenomenon is most prominent in the
intermediate regime of not exceedingly large magnetic fluxes, where the cavity
field dominates. The modification of the Hofstadter butterfly should again be
observable via transport measurements in the respective Wannier diagram
Wannier (1978). The Wannier diagram might provide a straightforward path to
the observation of cavity effects on the Hofstadter spectrum. Finally, we also
compared our results of the cavity engineering of the Hofstadter butterfly to
the Floquet driving of the butterfly Du  _et al._ (2018); Zhao  _et al._
(2020); Kooi  _et al._ (2018) and we found a basic consistency between QED-
Bloch theory and Floquet theory.

Semi-classical Limit of QED-Bloch Theory.—Further, we showed that our QED-
Bloch theory and the QED-Bloch central equation (44) are applicable also in
the semi-classical limit of no quantized field. In this limit our theory
recovers the standard Hofstadter butterfly fractal spectrum and provides a
first-principles framework for its description. We believe it can thus help to
understand recent experiments performed on Moiré systems Dean  _et al._
(2013); Wang  _et al._ (2015); Barrier  _et al._ (2020); Forsythe  _et al._
(2018). Moreover, in the semi-classical limit, the dual relation between the
minimal-coupling Hamiltonian and the tight-binding models with the Peierls
phase was described (see Fig. 7). This duality is of fundamental importance
for understanding the behavior of periodic materials perpendicular to
homogeneous magnetic fields and has been noted in experimental studies of the
Hofstadter butterfly Dean  _et al._ (2013).

Modification of the Integer Hall Effect.—As a further application of QED-Bloch
theory we considered the quantum Hall effect in the integer regime when
strongly coupled to the photon field of a cavity. In this case, our system
consists of non-interacting 2D Landau levels coupled to the cavity. Due to the
strong coupling between the cavity photons and the Landau levels hybrid
quasiparticles emerge, known as Landau polaritons Rokaj  _et al._ (2019);
Scalari  _et al._ (2012); Paravicini-Bagliani  _et al._ (2019). The formation
of the Landau polaritons modifies the plateaus of the Hall conductance,
$\sigma_{xy}=e^{2}\nu/h(1+\eta^{2})$, which now depend on the dimensionless
light-matter coupling $\eta$. This is a very important result because it
demonstrates that the long-range nature of the cavity photon-field circumvents
the topological protection of the integer Hall effect and actually modifies
this fundamental phenomenon of condensed matter physics. We believe that the
modification of the Hall conductance can be measured for a 2DEG inside a
cavity and could provide further insights on the recently observed
modification of the integer Hall effect Appugliese  _et al._ (2022). In
connection to these exciting, recent experiments, we would like to mention
that the modification of the Hall plateaus in Ref. Appugliese  _et al._ (2022)
is attributed to the presence of impurities and disorder in the 2DEG. The
cavity field benefits from the presence of impurities and mediates a hopping
mechanism which breaks the topological protection of the Hall effect Ciuti
(2021). In the absence of impurities, however, the proposed cavity-mediated
hopping vanishes. On the other hand, our work proposes a mechanism for the
modification of the Hall plateaus in the absence of impurities, for a clean
2DEG, without any disorder. Our prediction should be interpreted as a
renormalization effect (analogous to a Lamb shift Lamb and Retherford (1947))
of the Hall plateaus or equivalently as a screening effect on the external
magnetic field by the internal magnetic field of the cavity (see Fig. 11). It
is important to emphasize that the two approaches, i.e., the cavity-mediated
hopping and our renormalization effect, are not contradictory but rather
complementary as they apply to different physically relevant settings. In
conclusion, our findings provide new insights and pave the way for the
exploration of quantum Hall physics, in the integer and the fractional regime,
embedded in the field of cavity QED.

###### Acknowledgements.

We would like to thank J. Faist for useful discussions. This work was
supported by the European Research Council (ERC-2015-AdG694097), the Cluster
of Excellence “Advanced Imaging of Matter” (AIM), Grupos Consolidados
(IT1249-19), SFB925 “Light induced dynamics and control of correlated quantum
systems”, the Austrian Science Fund (J 4107-N27), and the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) via the Emmy Noether
program (SE 2558/2). The Flatiron Institute is a division of the Simons
Foundation.

##  Appendix A Displacement Operator Algebra

The aim of this appendix is to show how the matrix elements defined in Eq.
(39) can be computed and how the result in Eq. (40) was obtained. The matrix
elements that we are interested in computing are

$\langle\phi_{i}(v-A^{k_{x}}_{n})|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}|\phi_{j}(v-A^{k_{x}}_{n^{\prime}})\rangle.$
(74)

To calculate these matrix elements we will first perform a change of
coordinates, $s=v-A^{k_{x}}_{n}$, which will give us an overall phase
independent of the integration, and we have for the matrix elements of Eq.
(74)

$\mathrm{e}^{-\textrm{i}A^{k_{x}}_{n}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}}\langle\phi_{i}\left(s\right)|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}s}|\phi_{j}\left(s+A^{0}_{n-n^{\prime}}\right)\rangle.$
(75)

In order to compute the matrix elements above we will use the algebra of
displacement operators Cahill and Glauber (1969). The plane wave
$\exp(-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}s)$ can be written as a
displacement operator by using the expression for the coordinate $s$ in terms
of the annihilation and creation operators $\hat{b},\hat{b}^{\dagger}$
Griffiths (1995),

$s=\sqrt{\frac{\hbar}{2\mu\Omega}}\left(\hat{b}+\hat{b}^{\dagger}\right).$
(76)

Using the latter, we have for the plane wave in terms of the displacement
operator

$\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}s}=\hat{D}\left(-\textrm{i}\sqrt{\frac{\hbar}{2\mu\Omega}}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}\right).$
(77)

In addition, the wavefunction $\phi_{j}(s+A^{0}_{n-n^{\prime}})$ can be
written as

$\phi_{j}\left(s+A^{0}_{n-n^{\prime}}\right)=\hat{T}\left(A^{0}_{n-n^{\prime}}\right)\phi_{j}(s)$
(78)

using the translation operator, which is given by the expression Aschroft and
Mermin (1976)

$\hat{T}\left(A^{0}_{n-n^{\prime}}\right)=\exp\left(A^{0}_{n-n^{\prime}}\partial_{s}\right).$
(79)

The differential operator $\partial_{s}$ in terms of annihilation and creation
operators is

$\partial_{s}\equiv\frac{\partial}{\partial
s}=\sqrt{\frac{\mu\Omega}{2\hbar}}\left(\hat{b}-\hat{b}^{\dagger}\right).$
(80)

This implies that the translation operator can also be written as a
displacement operator Cahill and Glauber (1969),

$\hat{T}\left(A^{0}_{n-n^{\prime}}\right)=\hat{D}\left(-\sqrt{\frac{\mu\Omega}{2\hbar}}A^{0}_{n-n^{\prime}}\right).$
(81)

Using the expressions we derived in terms of the displacement operators we
obtain the following expression for the matrix elements in Eq. (74).

$\mathrm{e}^{-\textrm{i}A^{k_{x}}_{n}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}}\langle\phi_{i}|\hat{D}\left(-\frac{\textrm{i}\sqrt{\hbar}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}}{\sqrt{2\mu\Omega}}\right)\hat{D}\left(-\sqrt{\frac{\mu\Omega}{2\hbar}}A^{0}_{n-n^{\prime}}\right)|\phi_{j}\rangle$
(82)

We now use the formula from Cahill-Glauber Cahill and Glauber (1969)

$\hat{D}(\alpha)\hat{D}(\beta)=\hat{D}(\alpha+\beta)\exp((\alpha\beta^{*}-\alpha^{*}\beta)/2)$
(83)

and we obtain the following result for the product of displacement operators

$\displaystyle\hat{D}\left(-\textrm{i}\sqrt{\frac{\hbar}{2\mu\Omega}}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}\right)\hat{D}\left(-\sqrt{\frac{\mu\Omega}{2\hbar}}A^{0}_{n-n^{\prime}}\right)$
(84) $\displaystyle=$
$\displaystyle\hat{D}\left(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}\right)\mathrm{e}^{\frac{\textrm{i}}{2}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{0}_{n-n^{\prime}}},$

where the matrix elements $\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}$ are

$\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}=-\sqrt{\tfrac{\mu\Omega}{2\hbar}}A^{0}_{n-n^{\prime}}-\textrm{i}\sqrt{\tfrac{\hbar}{2\mu\Omega}}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}.$
(85)

We substitute the expression for the product of the displacement operators
into Eq. (82) and we have

$\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}\left(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}\right)|\phi_{j}\rangle.$
(86)

The matrix representation of this displacement operator in the basis
$\\{\phi_{i}(s)\\}$ is given by Cahill and Glauber (1969)

$\displaystyle\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=\sqrt{\frac{j!}{i!}}\alpha^{i-j}_{\mathbf{n}-\mathbf{n}^{\prime}}\mathrm{e}^{-\frac{|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2}}{2}}L^{(i-j)}_{j}(|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2})$

where $i\geq j$ and
$L^{(i-j)}_{j}(|\alpha_{\mathbf{n}-\mathbf{n}^{\prime}}|^{2})$ are the
associated Laguerre polynomials. We note that for $j>i$ one needs to take

$\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=(-1)^{j-i}\langle\phi_{j}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{i}\rangle^{*}$
(88)

because $\hat{D}^{\dagger}(\alpha)=\hat{D}(-\alpha)$ Cahill and Glauber
(1969). Finally, combining the result that we obtained in Eq. (86) with the
previous definitions, we obtain the expression for the matrix elements in Eq.
(74)

$\langle\phi_{i}(v-A^{k_{x}}_{n})|\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}v}|\phi_{j}(v-A^{k_{x}}_{n^{\prime}})\rangle=\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle.$
(89)

##  Appendix B Semi-classical Limit or No-quantized-field Limit

In this appendix we are interested in performing the semi-classical limit of
no quantized field for our QED-Bloch central equation (44). In this limit the
2D periodic material is only under the influence of the external magnetic
field $\mathbf{B}_{\textrm{ext}}$, while the quantized field
$\hat{\mathbf{\mathcal{A}}}$ goes to zero. Mathematically, this limit can be
performed by taking the diamagnetic frequency to zero, $\omega_{p}\rightarrow
0$, because the quantized field $\hat{\mathbf{\mathcal{A}}}$ is proportional
to $\omega_{p}$ (see Eq. (11)). This limit is equivalent to taking the light-
matter coupling $\eta$ to zero, $\eta\rightarrow 0$. Our QED-Bloch central
equation (44) was

$U^{\mathbf{k}}_{\mathbf{n},i}\Bigg{[}\frac{\hbar^{2}(k_{w}+G^{w}_{\mathbf{n}})^{2}}{2M}+\mathcal{E}_{i}-E_{\mathbf{k}}\Bigg{]}+\sum_{\mathbf{n}^{\prime},j}V_{\mathbf{n}-\mathbf{n}^{\prime}}U^{\mathbf{k}}_{\mathbf{n}^{\prime},j}\;\mathrm{e}^{-\textrm{i}G^{v}_{\mathbf{n}-\mathbf{n}^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}(\alpha_{\mathbf{n}-\mathbf{n}^{\prime}})|\phi_{j}\rangle=0$
(90)

and we will take the limit $\omega_{p}\rightarrow 0$ for all the parameters in
the above equation. First, we consider the limit $\omega_{p}\rightarrow 0$ for
the mass parameter $M$ defined in Eq. (21) and we find that
$\lim_{\omega_{p}\rightarrow 0}M=\infty$. This implies that the kinetic term
depending on $k_{w}$ in the central equation vanishes and the Fourier
components of our polaritonic Bloch wave no longer depend on $k_{w}$. Due to
the vanishing of the $w$ degree of freedom the index $m$ in the Bloch wave
becomes redundant,

$U^{k_{x},k_{w}}_{n,m,j}\longrightarrow U^{k_{x}}_{n,j}.$ (91)

Consequently, the central equation reduces to

$U^{k_{x}}_{n,i}\left(\mathcal{E}_{i}-E_{k_{x}}\right)+\sum_{n^{\prime},m^{\prime},j}V_{n-n^{\prime},m^{\prime}}U^{k_{x}}_{n^{\prime},j}\;\mathrm{e}^{-\textrm{i}G^{v}_{n-n^{\prime},m^{\prime}}A^{k_{x}}_{(n+n^{\prime})/2}}\langle\phi_{i}|\hat{D}(\alpha_{n-n^{\prime},m^{\prime}})|\phi_{j}\rangle=0.$
(92)

To obtain the result above we also relabelled the index
$-m^{\prime}\rightarrow m^{\prime}$. Now what is left to be done is to perform
the $\omega_{p}\rightarrow 0$ limit for the rest of the parameters in the
central equation which depend on $\omega_{p}$. By doing so we find

$\displaystyle\lim_{\omega_{p}\rightarrow
0}\Omega=\omega_{c},\;\;\lim_{\omega_{p}\rightarrow
0}\mu\Omega=\tfrac{2m_{e}}{\omega_{c}},$ (93)
$\displaystyle\lim_{\omega_{p}\rightarrow
0}G^{v}_{n-n^{\prime},m^{\prime}}=\tfrac{\sqrt{2}}{\omega_{c}}G_{n-n^{\prime},m^{\prime}}$
(94)

and

$\displaystyle\lim_{\omega_{p}\rightarrow 0}\alpha_{n-n^{\prime},m^{\prime}}$
$\displaystyle=\sqrt{\tfrac{\hbar}{2m_{e}\omega_{c}}}\left(-G^{x}_{n-n^{\prime}}-\textrm{i}G_{n-n^{\prime},m^{\prime}}\right)$
(95) $\displaystyle\equiv\beta_{n-n^{\prime},m^{\prime}}.$

Substituting all the above results and the definition for
$A^{k_{x}}_{(n+n^{\prime})/2}$ given by Eq. (31) we have

$U^{k_{x}}_{n,i}\left[\hbar\omega_{c}\left(i+\frac{1}{2}\right)-E_{k_{x}}\right]+\sum_{n^{\prime},m^{\prime},j}V_{n-n^{\prime},m^{\prime}}U^{k_{x}}_{n^{\prime},j}\mathrm{e}^{\frac{-\textrm{i}\hbar(k_{x}+\frac{1}{2}G^{x}_{n+n^{\prime}})G_{n-n^{\prime},m^{\prime}}}{m_{e}\omega_{c}}}\langle\phi_{i}|\hat{D}(\beta_{n-n^{\prime},m^{\prime}})|\phi_{j}\rangle=0.$
(96)

The central equation derived depends solely on electronic parameters like the
electronic crystal momentum $k_{x}$, the mass of the electron $m_{e}$ and the
cyclotron frequency $\omega_{c}=eB/m_{e}$, which is characteristic for
electrons in a constant magnetic field Landau and Lifshitz (1997). As a
consequence, the above central equation describes consistently the physics of
two-dimensional periodic systems in the presence of a perpendicular
homogeneous magnetic field. From this equation we can compute the energy bands
for such systems for all values of the magnetic field because our approach is
non-perturbative and does not rely on the magnetic translation group which
puts particular restrictions on the value of the magnetic field Brown (1964);
Zak (1964a, b).

For completeness, we would also like to give the expression of the polaritonic
Bloch ansatz defined in Eq. (35) in the limit of no quantized field. The QED-
Bloch ansatz depends on the polaritonic coordinates $w$ and $v$ defined in Eq.
(19). Taking the no-quantized-field limit $\omega_{p}\rightarrow 0$ the
coordinate $w$ vanishes while the coordinate $v$ becomes $v=-\omega_{c}y$.
Thus, we find that the polaritonic QED-Bloch ansatz in the limit of no
quantized field is

$\Psi_{k_{x}}(x,y)=\mathrm{e}^{\textrm{i}k_{x}}\sum_{n,j}U^{k_{x}}_{n,j}\mathrm{e}^{\textrm{i}G^{x}_{n}}\phi_{j}\left(-\frac{\omega_{c}y}{\sqrt{2}}-A^{k_{x}}_{n}\right).$
(97)

The above wavefunction corresponds to a correlated expansion between Bloch
waves in the $x$ coordinate and Landau levels
$\phi_{j}(-\omega_{c}y/\sqrt{2}-A^{k_{x}}_{n})$ in the $y$ coordinate. Such an
expansion has been used for the description of 2D materials in homogeneous
magnetic fields in several publications Pfannkuche and Gerhardts (1992);
Langbein (1969); Rauh (1975); Petschel and Geisel (1993) and central equations
analogous to Eq. (96) have been derived.

##  Appendix C The Effect of Many Modes

In this appendix we will look into the effect of many modes for the strength
of the light-matter coupling and the corresponding implication for the Hall
conductance. To do so, we will not try to take into account exactly a finite
amount of modes, but we will rather follow an effective approach in which the
dependence of the single-mode coupling constant $\eta$ on the photonic momenta
is introduced back and then the sum of the single-mode couplings over the
photonic momenta defines an effective many-mode coupling. This effective
approach was also followed in Ref. Rokaj  _et al._ (2022) and it was shown
that it captures the exact running of the light-matter coupling as a function
of the photonic upper cutoff and recovers well-known perturbative mass-
renormalization results of quantum field theory.

The single-mode coupling $\eta$ depends on the cavity frequency
$\omega_{\textrm{cav}}$ via the diamagnetic frequency $\omega_{p}$ defined in
Eq. (45). The cavity frequency itself is a function of the photonic momentum
$\kappa_{z}=\pi n_{z}/L_{z}$, $\omega_{\textrm{cav}}=c|\kappa_{z}|$, where we
take only the out-of-plane momenta into account. Substituting this expression
for $\omega_{\textrm{cav}}$ into the definition of $\eta$ we find
$\eta^{2}=e^{2}n_{\textrm{2D}}|\kappa_{z}|/m_{e}\pi\epsilon_{0}\omega^{2}_{c}$.
Then the effective many-mode light-matter coupling constant is

$\eta^{2}(\Lambda_{0})=\frac{e^{2}n_{\textrm{2D}}}{m_{e}\pi\epsilon_{0}\omega^{2}_{c}}\frac{\pi}{L_{z}}\sum^{\Lambda_{0}}_{n_{z}=-\Lambda_{0}}|n_{z}|=\frac{\omega^{2}_{p}}{\omega^{2}_{c}}\Lambda_{0}(\Lambda_{0}+1)$
(98)

where the photonic momenta $\kappa_{z}$ were summed up to the upper cutoff
$\Lambda=\Lambda_{0}\pi/L_{z}$, with $\Lambda_{0}\in\mathbb{N}$, which is
defined as an arbitrary multiple of the inverse of the cavity length $L_{z}$.
For large $\Lambda_{0}\gg 1$ we can approximately consider
$\Lambda_{0}+1\approx\Lambda_{0}$ and the many-mode effective coupling takes
the simple form

$\eta^{2}(\Lambda_{0})=\Lambda^{2}_{0}\frac{\omega^{2}_{p}}{\omega^{2}_{c}}.$
(99)

The many-mode effective coupling $\eta(\Lambda_{0})$ can be rewritten as the
the ratio between the strength of the external magnetic field $B$ and the
many-mode cavity magnetic field $B(\Lambda_{0})=\Lambda_{0}m_{e}\omega_{p}/e$
as $\eta(\Lambda_{0})=B(\Lambda_{0})/B$. We note that the many-mode cavity
field $B(\Lambda_{0})$ is merely a multiple of the of the single-mode cavity
magnetic field in the singe-mode case $B_{\textrm{cav}}=m_{e}\omega_{p}/e$.
This means that the inclusion of many modes acts as an amplifier for the
strength of the cavity field. Finally, replacing the single-mode coupling
$\eta$ in the formula for the cavity-modified Hall conductance in Eq. (71)
with the effective many-mode coupling $\eta(\Lambda_{0})$ one
straightforwardly obtains the effect of many modes for the modified Hall
conductance.

## References

  * von Klitzing  _et al._ (2020) K. von Klitzing, T. Chakraborty, P. Kim, V. Madhavan, X. Dai, J. McIver, Y. Tokura, L. Savary, D. Smirnova, A. M. Rey, C. Felser, J. Gooth, and X. Qi, “40 years of the quantum Hall effect,” Nat. Rev. Phys. 2, 38–41 (2020).
  * Klitzing  _et al._ (1980) K. v. Klitzing, G. Dorda, and M. Pepper, “New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance,” Phys. Rev. Lett. 45, 494–497 (1980).
  * Stormer  _et al._ (1999) Horst L. Stormer, Daniel C. Tsui, and Arthur C. Gossard, “The fractional quantum Hall effect,” Rev. Mod. Phys. 71, S298–S305 (1999).
  * Laughlin (1983) R. B. Laughlin, “Anomalous quantum Hall effect: An incompressible quantum fluid with fractionally charged excitations,” Phys. Rev. Lett. 50, 1395–1398 (1983).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, “Quantum spin Hall effect in graphene,” Phys. Rev. Lett. 95, 226801 (2005).
  * Chang  _et al._ (2013) C-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L-L. Wang, Z-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S-C. Zhang, K. He, Y. Wang, L. Lu, X-C. Ma, and Q-K. Xue, “Experimental observation of the quantum anomalous Hall effect in a magnetic topological insulator,” Science 340, 167–170 (2013).
  * McIver  _et al._ (2020) J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri, “Light-induced anomalous Hall effect in graphene,” Nat. Phys. 16, 38–41 (2020).
  * Hofstadter (1976) D. R. Hofstadter, “Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields,” Phys. Rev. B 14, 2239–2249 (1976).
  * Avila and Jitomirskaya (2006) A. Avila and S. Jitomirskaya, “Solving the ten martini problem,” in  _Mathematical Physics of Quantum Mechanics: Selected and Refereed Lectures from QMath9_, edited by Joachim Asch and Alain Joye (Springer Berlin Heidelberg, 2006) pp. 5–16.
  * Andrei  _et al._ (2021) E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, “The marvels of Moiré materials,” Nat. Rev. Mater. 6, 201–206 (2021).
  * Kennes  _et al._ (2021) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, “Moiré heterostructures as a condensed-matter quantum simulator,” Nat. Phys. 17, 155–163 (2021).
  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim, “Hofstadter’s butterfly and the fractal quantum Hall effect in Moiré superlattices,” Nature 497 (2013).
  * Wang  _et al._ (2015) L. Wang, Y. Gao, B. Wen, Z. Han, T. Taniguchi, K. Watanabe, M. Koshino, J. Hone, and C. R. Dean, “Evidence for a fractional fractal quantum Hall effect in graphene superlattices,” Science 350, 1231–1234 (2015).
  * Forsythe  _et al._ (2018) C. Forsythe, X. Zhou, K. Watanabe, Z. Han, T. Taniguchi, K. Watanabe, M. Koshino, T. Taniguchi, A. Pasupathy, P. Moon, M. Koshino, P. Kim, and C. R. Dean, “Band structure engineering of 2D materials using patterned dielectric superlattices,” Nature Nanotech. 13, 566–571 (2018).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, “Realization of the Hofstadter hamiltonian with ultracold atoms in optical lattices,” Phys. Rev. Lett. 111, 185301 (2013).
  * Miyake  _et al._ (2013) Hirokazu Miyake, Georgios A. Siviloglou, Colin J. Kennedy, William Cody Burton, and Wolfgang Ketterle, “Realizing the Harper hamiltonian with laser-assisted tunneling in optical lattices,” Phys. Rev. Lett. 111, 185302 (2013).
  * Cohen-Tannoudji  _et al._ (1997) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms-Introduction to Quantum Electrodynamics_ (Wiley-VCH, 1997).
  * Weinberg (2005) S. Weinberg,  _The Quantum Theory of Fields_ (Cambridge University Press, 2005).
  * Spohn (2004) H. Spohn,  _Dynamics of Charged Particles and their Radiation Field_ (Cambridge university press, 2004).
  * Kockum  _et al._ (2019) A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, “Ultrastrong coupling between light and matter,” Nat. Rev. Phys. 1, 19–40 (2019).
  * Basov  _et al._ (2021) D. N. Basov, Ana Asenjo-Garcia, P. James Schuck, Xiaoyang Zhu, and Angel Rubio, “Polariton panorama,” Nanophotonics 10, 549–577 (2021).
  * Ebbesen (2016) T. W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspective,” Acc. Chem. Res. 49, 2403–2412 (2016).
  * Hutchison  _et al._ (2013) J.A. Hutchison, A. Liscio, T. Schwartz, A. Canaguier-Durand, C. Genet, V. Palermo, P. Samorì, and T. W. Ebbesen, “Tuning the work-function via strong coupling,” Adv. Mater. 25, 2481–2485 (2013).
  * Hutchison  _et al._ (2012) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, “Modifying chemical landscapes by coupling to vacuum fields,” Angew. Chem. Int. Ed. 51, 1592–1596 (2012).
  * Orgiu  _et al._ (2015) E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samorì, and T. W. Ebbesen, “Conductivity in organic semiconductors hybridized with the vacuum field,” Nat. Mater. 14, 1123–1129 (2015).
  * Feist  _et al._ (2017) J. Feist, J. Galego, and F. J. Garcia-Vidal, “Polaritonic chemistry with organic molecules,” ACS Photonics 5, 205–216 (2017).
  * Galego  _et al._ (2016) J. Galego, F. J. Garcia-Vidal, and J. Feist, “Suppressing photochemical reactions with quantized light fields,” Nat. Commun. 7, 13841 (2016).
  * Flick  _et al._ (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry,” PNAS 114, 3026–3034 (2017).
  * Schäfer  _et al._ (2019) C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Modification of excitation and charge transfer in cavity quantum-electrodynamical chemistry,” PNAS 116, 4883–4892 (2019).
  * Latini  _et al._ (2019) S. Latini, E. Ronca, U. De Giovannini, H. Hübener, and A. Rubio, “Cavity control of excitons in two-dimensional materials,” Nano Lett. 19, 3473–3479 (2019).
  * Förg  _et al._ (2019) M. Förg, L. Colombier, R. K. Patel, J. Lindlau, A. D. Mohite, H. Yamaguchi, M. M. Glazov, D. Hunger, and A. Högele, “Cavity-control of interlayer excitons in Van der Waals heterostructures,” Nat. Commun. 10, 3697 (2019).
  * Kasprzak  _et al._ (2006) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and Le Si Dang, “Bose–Einstein condensation of exciton polaritons,” Nature 443, 409–414 (2006).
  * Keeling and Kena-Cohen (2020) J. Keeling and S. Kena-Cohen, “Bose–Einstein condensation of exciton-polaritons in organic microcavities,” Annual Review of Physical Chemistry 71, 435–459 (2020).
  * Schlawin  _et al._ (2019) F. Schlawin, A. Cavalleri, and D. Jaksch, “Cavity-mediated electron-photon superconductivity,” Phys. Rev. Lett. 122, 133602 (2019).
  * Cotleţ  _et al._ (2016) O. Cotleţ, S. Zeytinoǧlu, M. Sigrist, E. Demler, and A. Imamoǧlu, “Superconductivity and other collective phenomena in a hybrid bose-fermi mixture formed by a polariton condensate and an electron system in two dimensions,” Phys. Rev. B 93, 054510 (2016).
  * Sentef  _et al._ (2018) M. A. Sentef, M. Ruggenthaler, and A. Rubio, “Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity,” Sci. Adv. 4, eaau6969 (2018).
  * Curtis  _et al._ (2019) J. B. Curtis, Z. M. Raines, A. A. Allocca, M. Hafezi, and V. M. Galitski, “Cavity quantum Eliashberg enhancement of superconductivity,” Phys. Rev. Lett. 122, 167002 (2019).
  * Thomas  _et al._ (2019) A. Thomas, E. Devaux, K. Nagarajan, T. Chervy, M. Seidel, D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genet, G. Pupillo, and T. W. Ebbesen, “Exploring superconductivity under strong coupling with the vacuum electromagnetic field,” arXiv:1911.01459 [cond-mat.supr-con] (2019).
  * Hübener  _et al._ (2021) H. Hübener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, “Engineering quantum materials with chiral optical cavities,” Nat. Mater. 20, 438–442 (2021).
  * Petersen  _et al._ (2014) J. Petersen, J. Volz, and A Rauschenbeutel, “Chiral nanophotonic waveguide interface based on spin-orbit interaction of light,” Science 346, 67–71 (2014).
  * Zhang  _et al._ (2019) F. Zhang, J. Ren, L. Shan, X. Duan, Y. Li, T. Zhang, Q. Gong, and Y. Gu, “Chiral cavity quantum electrodynamics with coupled nanophotonic structures,” Phys. Rev. A 100, 053841 (2019).
  * Lodahl  _et al._ (2017) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, “Chiral quantum optics,” Nature 541, 473–480 (2017).
  * Wang  _et al._ (2019) X. Wang, E. Ronca, and M. A. Sentef, “Cavity quantum electrodynamical Chern insulator: Towards light-induced quantized anomalous Hall effect in graphene,” Phys. Rev. B 99, 235156 (2019).
  * Latini  _et al._ (2021) S. Latini, D. Shin, S. A. Sato, C. Schäfer, U. De Giovannini, H. Hübener, and A. Rubio, “The ferroelectric photo-groundstate of $\textrm{SrTiO}_{3}$: Cavity materials engineering,” arXiv:2101.11313 [cond-mat.mtrl-sci] (2021).
  * Ashida  _et al._ (2020) Y. Ashida, A. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, “Quantum electrodynamic control of matter: Cavity-enhanced ferroelectric phase transition,” Phys. Rev. X 10, 041027 (2020).
  * Hagenmüller  _et al._ (2010) D. Hagenmüller, S. De Liberato, and C. Ciuti, “Ultrastrong coupling between a cavity resonator and the cyclotron transition of a two-dimensional electron gas in the case of an integer filling factor,” Phys. Rev. B 81, 235303 (2010).
  * Rokaj  _et al._ (2019) V. Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, and A. Rubio, “Quantum electrodynamical Bloch theory with homogeneous magnetic fields,” Phys. Rev. Lett. 123, 047202 (2019).
  * Keller  _et al._ (2020) J. Keller, G. Scalari, F. Appugliese, S. Rajabali, M. Beck, J. Haase, C. A. Lehner, W. Wegscheider, M. Failla, M. Myronov, D. R. Leadley, J. Lloyd-Hughes, P. Nataf, and J. Faist, “Landau polaritons in highly nonparabolic two-dimensional gases in the ultrastrong coupling regime,” Phys. Rev. B 101, 075301 (2020).
  * Scalari  _et al._ (2012) G. Scalari, C. Maissen, D. Turčinková, D. Hagenmüller, S. De Liberato, C. Ciuti, C. Reichl, D. Schuh, W. Wegscheider, M. Beck, and J. Faist, “Ultrastrong coupling of the cyclotron transition of a 2D electron gas to a THz metamaterial,” Science 335, 1323–1326 (2012).
  * Li  _et al._ (2018) X. Li, M. Bamba, Q. Zhang, S. Fallahi, G. C. Gardner, W. Gao, M. Lou, K. Yoshioka, M. J. Manfra, and J. Kono, “Vacuum Bloch–Siegert shift in Landau polaritons with ultra-high cooperativity,” Nature Photon. 12, 324–329 (2018).
  * Ravets  _et al._ (2018) S. Ravets, P. Knüppel, S. Faelt, O. Cotlet, M. Kroner, W. Wegscheider, and A. Imamoglu, “Polaron polaritons in the integer and fractional quantum Hall regimes,” Phys. Rev. Lett. 120, 057401 (2018).
  * Smolka  _et al._ (2014) S. Smolka, W. Wuester, F. Haupt, S. Faelt, W. Wegscheider, and A. Imamoglu, “Cavity quantum electrodynamics with many-body states of a two-dimensional electron gas,” Science 346, 332–335 (2014).
  * Paravicini-Bagliani  _et al._ (2019) G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, S. Fallahi, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, “Magneto-transport controlled by Landau polariton states,” Nat. Phys. 15, 186–190 (2019).
  * Ciuti (2021) Cristiano Ciuti, “Cavity-mediated electron hopping in disordered quantum hall systems,” Phys. Rev. B 104, 155307 (2021).
  * Appugliese  _et al._ (2022) F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist, “Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect,” Science 375, 1030–1034 (2022).
  * Rokaj (2021) V. Rokaj, “Condensed matter systems in cavity quantum electrodynamics,” Ph.D. Thesis Published by the University of Hamburg (arXiv:2201.01331 [quant-ph]) (2021).
  * Landau and Lifshitz (1997) L. D. Landau and E. M. Lifshitz,  _Quantum Mechanics, Third Edition: Non-relativistic Theory_ (Pergamon Press, 1997).
  * Laughlin (1981) R. B. Laughlin, “Quantized Hall conductivity in two dimensions,” Phys. Rev. B 23, 5632–5633 (1981).
  * Rokaj  _et al._ (2018) V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, “Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy,” J. Phys. B: At. Mol. Opt. Phys 51, 034005 (2018).
  * Greiner and Reinhardt (1996) W. Greiner and J. Reinhardt,  _Field Quantization_ (Springer, 1996).
  * Faisal (1987) F. H. Faisal,  _Theory of Multiphoton Processes_ (Springer, Berlin, 1987).
  * Ruggenthaler  _et al._ (2018) M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, “From a quantum-electrodynamical light–matter description to novel spectroscopies,” Nat. Rev. Chem. 2, 0118 (2018).
  * Ribeiro  _et al._ (2021) R. F. Ribeiro, J. A. Campos-Gonzalez-Angulo, N. C. Giebink, W Xiong, and J. Yuen-Zhou, “Enhanced optical nonlinearities under collective strong light-matter coupling,” Phys. Rev. A 103, 063111 (2021).
  * Dicke (1954) R. H. Dicke, “Coherence in spontaneous radiation processes,” Phys. Rev. 93, 99–110 (1954).