# Cavity Quantum Electrodynamics with Hyperbolic van der Waals Materials

Yuto Ashida<EMAIL_ADDRESS>Department of Physics, University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Institute for Physics of
Intelligence, University of Tokyo, 7-3-1 Hongo, Tokyo 113-0033, Japan  Ataç
$\dot{\mathrm{I}}$mamo$\breve{\mathrm{g}}$lu  Institute of Quantum
Electronics, ETH Zurich, CH-8093 Zürich, Switzerland  Eugene Demler  Institute
for Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland

###### Abstract

The ground-state properties and excitation energies of a quantum emitter can
be modified in the ultrastrong coupling regime of cavity quantum
electrodynamics (QED) where the light-matter interaction strength becomes
comparable to the cavity resonance frequency. Recent studies have started to
explore the possibility to control an electronic material by embedding it in a
cavity that confines electromagnetic fields in deep subwavelength scales.
Currently, there is a strong motivation to realize ultrastrong-coupling cavity
QED in the terahertz (THz) range, since most of the elementary excitations of
quantum materials are in this frequency window. We propose and analyze an
ideal platform to achieve this aim where a two-dimensional electronic material
is encapsulated by a planar cavity consisting of ultrathin polar van der Waals
crystals. As a concrete setup, we show that nanometer-thick hexagonal boron
nitride layers allow for reaching the ultrastrong coupling regime for single-
electron cyclotron resonance in a bilayer graphene. The proposed cavity
setting can be realized by a wide variety of thin dielectric materials with
hyperbolic dispersions. Consequently, van der Waals heterostructures could
provide an ideal playground for exploring the ultrastrong-coupling physics of
cavity QED materials.

Strong coupling regime of cavity quantum electrodynamics (QED), where the
emitter-cavity coupling strength exceeds decay rates, has played a central
role in quantum information science. For instance, cavity-mediated coherent
interaction between two distant qubits has allowed for implementing two-qubit
gates with high fidelity DiCarlo  _et al._ (2009); Blais  _et al._ (2021).
Moreover, cavity-mediated Raman transitions have the potential to realize all-
to-all coupling with a tunable range and strength, providing an indispensable
tool for quantum simulators Stute  _et al._ (2012); Vaidya  _et al._ (2018).
Recent studies have started to explore the possibility of further increasing
the coupling strength so that it becomes comparable to the elementary
excitation energies. Many of the common simplifications in cavity QED fail in
this regime, rendering the theoretical analysis challenging Hagenmüller  _et
al._ (2010); De Bernardis  _et al._ (2018); Le Boite (2020); Ashida  _et al._
(2021). A remarkable feature here is that ultrastrong light-matter coupling
can alter the ground-state electronic properties due to virtual processes
where both emitters and cavity photons are excited Garcia-Vidal  _et al._
(2021); Schlawin  _et al._ (2022); Bloch  _et al._ (2022). Consequently, a
natural question to address is if and when the ultrastrong coupling regime can
be attained and used to control material properties simply by cavity
confinement in the absence of an external driving.

A quick calculation for cavities supporting purely photonic excitations
suggests that the ultrastrong coupling regime is out of reach due to the
smallness of the fine structure constant Devoret  _et al._ (2007). However,
this limit can be overcome in structured electromagnetic environments because
hybridization with matter excitations enables one to control cavity
frequencies independently of wavelengths. For instance, in superconducting
microwave circuits, a large kinetic inductance allows for high-impedance
electromagnetic excitations, leading to the single-quantum ultrastrong
coupling in the microwave range Forn-Díaz  _et al._ (2017); Yoshihara  _et
al._ (2017); Forn-Díaz  _et al._ (2019); Frisk Kockum  _et al._ (2019).

Meanwhile, many of the elementary excitations in quantum materials are in the
terahertz (THz) regime. Recent experimental and theoretical studies have shown
the potential of utilizing the cavity confinement as a means to modify such
excitations Scalari  _et al._ (2014); Smolka  _et al._ (2014); Maissen  _et
al._ (2014); Zhang  _et al._ (2016); Bayer  _et al._ (2017); Keller  _et al._
(2020); Thomas  _et al._ (2021); Flick  _et al._ (2017); Rokaj  _et al._
(2019); Juraschek  _et al._ (2021); Román-Roche and Zueco (2022), which shows
great promise for controlling quantum many body states and endowing them with
new functionalities, including superconductivity Schlawin  _et al._ (2019);
Curtis  _et al._ (2019); Sentef  _et al._ (2018); Gao  _et al._ (2021),
ferroelectricity Ashida  _et al._ (2020a); Latini  _et al._ (2021); Lenk  _et
al._ (2022), magnetotransport Paravicini-Bagliani  _et al._ (2019); Appugliese
_et al._ (2022); Rokaj  _et al._ (2022) and topological properties Wang  _et
al._ (2019); Masuki and Ashida (2022). In view of the prospects of observing
these intriguing phenomena in cavity QED materials, there is a strong
motivation to examine the feasibility of attaining the single-quantum
ultrastrong coupling at THz frequencies.

Figure 1:  Schematic figure illustrating the proposed planar cavity setup
consisting of thin polar van der Waals crystals (green shaded) whose optical
axis is along $z$ direction. In the narrow air gap between the two slabs, a 2D
material (red shaded) is inserted and the electron there can ultrastrongly
couple to electromagnetic fields of tightly confined hyperbolic polaritons
(blue solid curve). The thickness (lateral extension) of surrounding materials
is $d$ ($L$). (b) Out-of-plane (solid curve) and in-plane (dashed curve)
permittivities $\epsilon_{z,t}$ of hyperbolic materials. The Restrahlen bands
(red shaded) appear above each of the out-of-plane and in-plane phonon
resonances $\Omega_{z,t}$.

The aim of this Letter is to propose and analyze a planar cavity setup
consisting of polar van der Waals (vdW) crystals as an ideal platform to
explore the ultrastrong-coupling physics of cavity QED materials (Fig. 1a). We
point out that one can attain a single-quantum ultrastrong coupling by
embedding a 2D material with electronic transitions in the THz regime into
ultrathin hexagonal boron nitride (h-BN) layers. A special feature here is
that electrons in the 2D material couple to the electric field component of
tightly confined hyperbolic phonon polaritons. The strong photon-phonon
hybridization together with the low frequencies of polaritons then results in
the significantly enhanced coupling strength over a broad range of momenta.
This should be contrasted to conventional polar dielectrics with an isotropic
dispersion, where sizable light-matter hybridization takes place in a limited
range of momenta since light and matter are almost decoupled at large momenta
due to the fast speed of light.

As a proof-of-concept demonstration, we consider the case of a 2D electron
with parabolic dispersion under the static magnetic field. We analyze the
hybridization between hyperbolic phonon polaritons in h-BN and the cyclotron
motion of the parabolic electrons. We show that the single-electron
ultrastrong coupling regime is within the reach provided that the thicknesses
of h-BN layers are chosen to be nanometer-scale. This consideration is
motivated by recent advances demonstrating ultrasmall mode volumes of
hyperbolic phonon polaritons in h-BN nanostructures Caldwell  _et al._ (2014);
Dai  _et al._ (2014); Giles  _et al._ (2018); Dai  _et al._ (2019); Ma  _et
al._ (2022); Sheinfux  _et al._ (2022). While we present quantitative
estimates for h-BN nanocavities for the sake of concreteness, the proposed
cavity scheme is applicable to the majority of polar vdW crystals Geim and
Grigorieva (2013), which exhibit hyperbolic polaritons originating from
distinct in- and out-of-plane infrared-active phonons.

Hyperbolic phonon polaritons.— We begin our analysis by reviewing the
properties of a planar cavity made out of layered thin polar vdW materials.
Due to the weakness of interlayer coupling, such materials naturally possess
two types of optical phonons corresponding to in-plane and out-of-plane ionic
oscillations in the THz or mid-infrared regimes. This leads to the uniaxial
anisotropy characterized by the out-of-plane (in-plane) dielectric
permittivities $\epsilon_{z}$ ($\epsilon_{t}$). In the frequency windows above
each of the phonon resonances, the two dielectric responses can have opposite
signs (Fig. 1b). This unique feature leads to the hyperbolic isofrequency
surfaces defined by the dispersion relation of the transverse magnetic modes,

$\displaystyle\frac{|{\boldsymbol{q}}|^{2}}{\epsilon_{z}}+\frac{|{\boldsymbol{\kappa}}|^{2}}{\epsilon_{t}}=\frac{\omega^{2}}{\epsilon_{0}c^{2}},$
(1)

where $\boldsymbol{q}$ ($\boldsymbol{\kappa}$) is the in-plane (out-of-plane)
wavevector, $\epsilon_{0}$ is the vacuum permittivity, and $c$ is the speed of
light. The opposite signs of $\epsilon_{z,t}$ allow for excitations to have
low frequencies even at large momenta. The resulting hybridized elementary
excitations of photons and optical phonons are known as hyperbolic phonon
polaritons. The corresponding dispersions in this range of frequencies are
called the Reststrahlen bands of either Type I when
$\epsilon_{z}<0,\epsilon_{t}>0$ or Type II when
$\epsilon_{z}>0,\epsilon_{t}<0$.

As a representative hyperbolic material, h-BN has a strong crystalline
anisotropy leading to two spectrally distinct Reststrahlen bands with ultralow
loss Caldwell  _et al._ (2014); Dai  _et al._ (2014); Giles  _et al._ (2018);
Dai  _et al._ (2019); Ma  _et al._ (2022); Sheinfux  _et al._ (2022). Below we
focus on the electromagnetic couplings with Type I h-BN hyperbolic modes
$\omega_{\boldsymbol{q}n}$ lying in the narrow frequency window above the out-
of-plane phonon frequency $\Omega_{z}=41.1\,{\rm THz}$ (cf. top panel in Fig.
2). As we discuss below, these modes exhibit several advantages for the
purpose of attaining the ultrastrong couplings thanks to their relatively low
frequencies and sizable photon components over a broad range of momenta $q$.

Figure 2:  (Top) Dispersions $\omega_{{\boldsymbol{q}}n}$ of the hyperbolic
phonon polaritons in the planar cavity setting. The inset shows the spatial
profiles of the in-plane electric fields along $z$ direction for each mode at
$qd=4$. For the sake of visibility, only the three lowest modes are plotted.
(Bottom) Inverse of the square root of the effective dimensionless confinement
length, $\sqrt{d/d_{{\boldsymbol{q}}n}}$, characterizing the dimensionless
single-electron coupling strength for each mode. The maximum value in the
principal branch $n=1$ is reached at $q=q^{*}$.

Single-electron ultrastrong coupling with hyperbolic materials.— We consider
the planar cavity setup where a 2D material (e.g., a bilayer graphene) is
inserted into the narrow air gap between two thin h-BN slabs whose thickness
and lateral extensions are denoted by $d$ and $L$, respectively (Fig. 1a). Our
starting point is the following cavity QED Hamiltonian of the 2D parabolic
electron interacting with hyperbolic phonon polaritons:

$\displaystyle\hat{H}=\frac{\left[\hat{\boldsymbol{p}}+e\boldsymbol{A}_{s}(\hat{\boldsymbol{r}})+e\hat{\boldsymbol{A}}(\hat{\boldsymbol{r}})\right]^{2}}{2m}+\hat{H}_{\rm
pol},$ (2)

which can be derived from an effective theory of uniaxial polar dielectrics
coupled to quantized electromagnetic fields in the Coulomb gauge Cohen-
Tannoudji  _et al._ (1989); SM (1). Here, $e$ is the elementary charge, $m$ is
the electron mass, $\hat{\boldsymbol{r}}$ and $\hat{\boldsymbol{p}}$ are the
electron position and momentum operators in the 2D lateral directions,
respectively, $\boldsymbol{A}_{s}$ represents an arbitrary static field, and
$\hat{H}_{\rm pol}$ is the free polariton Hamiltonian,

$\displaystyle\hat{H}_{{\rm
pol}}=\sum_{\boldsymbol{q}n}\hbar\omega_{\boldsymbol{q}n}\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}\hat{\gamma}_{\boldsymbol{q}n},$
(3)

where $\hat{\gamma}_{\boldsymbol{q}n}$
($\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}$) annihilates (creates) a
hyperbolic polariton with the in-plane wavevector $\boldsymbol{q}$ in the
branch $n\in{\mathbb{N}}$. The 2D vector field
$\hat{\boldsymbol{A}}({\boldsymbol{r}})$ is obtained by projecting the vector
potential onto the 2D plane where the electron is placed, and can be expanded
in terms of the polariton operators by

$\displaystyle\hat{\boldsymbol{A}}(\hat{\boldsymbol{r}})=\sum_{\boldsymbol{q}n}{\cal
A}_{\boldsymbol{q}n}\boldsymbol{e}_{\boldsymbol{q}}\left(\hat{\gamma}_{\boldsymbol{q}n}e^{i\boldsymbol{q}\cdot\hat{\boldsymbol{r}}}+\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}e^{-i\boldsymbol{q}\cdot\hat{\boldsymbol{r}}}\right),$
(4)

where ${\cal
A}_{\boldsymbol{q}n}\simeq\sqrt{{\hbar}/({2L^{2}\epsilon_{0}\omega_{\boldsymbol{q}n}d_{\boldsymbol{q}n}})}$
is the amplitude with the effective confinement length $d_{\boldsymbol{q}n}$
whose value is characterized by the polariton mode function. The effective
polarization vector becomes
$\boldsymbol{e}_{\boldsymbol{q}}\equiv\boldsymbol{q}/|{\boldsymbol{q}}|$
because, for the symmetric arrangement we assumed, the electric fields of the
polaritons only have in-plane components along the propagation direction in
the 2D plane where the electronic material is located. In general, while the
vector potential is originally a 3D transverse vector field in the Coulomb
gauge, it can effectively acquire longitudinal components when projected onto
the 2D tangential plane.

Figure 2 shows the results for Type I h-BN hyperbolic modes. The dispersions
(top panel) start from the longitudinal phonon frequency and saturate to
$\Omega_{z}$ at high $q\equiv|{\boldsymbol{q}}|$. Naturally, these hybridized
modes are almost purely longitudinal (transverse) phonon excitations in the
limit $q\to 0$ ($q\to\infty$), which do not allow for a strong coupling to
electrons in the air gap. Crucially, however, except for these two limits,
there still exist the nonvanishing photon contributions leading to the sizable
electric-dipole couplings at intermediate momenta (bottom panel). This is
because the in-plane component, which couples to the 2D electron, has almost
equal photonic and phononic contents while the out-of-plane component is
predominantly phonon-like SM (1).

Figure 3:  Magnetoabsorption spectra $A(\omega)$ of the 2D electron in the
cavity plotted against a cyclotron resonance $\omega_{c}=eB/m$ at different
cavity thicknesses $d$. The blue solid curve in each panel corresponds to the
bound-state energy of the Landau polariton. We impose the periodic boundary
conditions on the lateral directions and fix the aspect ratio $L/d=\pi$. We
use the h-BN parameters in Ref. Caldwell  _et al._ (2014) and set the in-plane
momentum cutoff $\Lambda=2{\rm\,nm}^{-1}$.

To further proceed, we use the unitary transformation
$\hat{U}=e^{-i\hat{\boldsymbol{r}}\cdot\hat{\boldsymbol{P}}_{b}/\hbar}$ with
$\hat{\boldsymbol{P}}_{b}=\sum_{\boldsymbol{q}n}\hbar\boldsymbol{q}\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}\hat{\gamma}_{\boldsymbol{q}n}$
Lee  _et al._ (1953), resulting in the Hamiltonian
$\hat{H}_{U}\equiv\hat{U}^{\dagger}\hat{H}\hat{U}$ given by

$\displaystyle\hat{H}_{U}=\frac{\left[\hat{\boldsymbol{p}}\\!+\\!e\boldsymbol{A}_{s}(\hat{\boldsymbol{r}})\\!+\\!e\hat{\boldsymbol{A}}(\boldsymbol{0})\\!-\\!\hat{\boldsymbol{P}}_{b}\right]^{2}}{2m}+\hat{H}_{\rm
pol}.$ (5)

In this way, the electron-coordinate dependence in the quantized vector
potential can be eliminated at the expense of generating the polariton
momentum $\hat{\boldsymbol{P}}_{b}$, leading to the nonlocal interaction among
polaritons mediated by the electron. The dimensionful coupling strength
between the electron and each of the dynamical quantized electromagnetic modes
is given by

$\displaystyle g_{\boldsymbol{q}n}$ $\displaystyle=$ $\displaystyle e{\cal
A}_{\boldsymbol{q}n}\sqrt{\frac{\omega_{\boldsymbol{q}n}}{m\hbar}}.$ (6)

Since $g_{\boldsymbol{q}n}$ characterizes the single-electron coupling
strength rather than the collective one, it depends on the lateral size $L$
through $g_{\boldsymbol{q}n}\propto L^{-1}$. Consequently, a natural measure
for the effective coupling strength between the 2D electron and the continuum
of polariton modes is given by $g_{{\rm
eff}}\equiv\sqrt{\sum_{\boldsymbol{q}n}g_{\boldsymbol{q}n}^{2}}$, which scales
as $g_{\rm eff}\propto{\rm O}(L^{0})$ Ashida  _et al._ (2022).

Previously, the deep strong coupling regime has been experimentally realized
in superconducting circuits, where $g_{\rm eff}$ becomes comparable to the
microwave photon frequency. In the present setting, the use of ultrathin h-BN
slabs with nanometer-scale thicknesses enables one to reach the deep or
extremely strong coupling regimes, where $g_{{\rm eff}}$ becomes comparable to
or even exceeds a THz cavity frequency. For instance, using the h-BN
parameters Caldwell  _et al._ (2014), we estimate coupling strengths of order
$g_{\boldsymbol{q}n}/\Omega_{z}=1.7\times\sqrt{(10{\rm\,nm})^{3}/(L^{2}d_{\boldsymbol{q}n})}$,
which, together with the results of the effective confinement length
$d_{{\boldsymbol{q}}n}$ in Fig. 2b, can lead to $g_{\rm eff}\sim\Omega_{z}$ in
nanoscale heterostructures. More specifically, one can attain $g_{\rm
eff}\simeq 2.0\,\Omega_{z}$ for the $n=1$ principal branch when the cavity
thickness (in-plane momentum cutoff) is set to be $d=5{\rm\,nm}$
($\Lambda=2{\rm\,nm}^{-1}$). As demonstrated below, one important consequence
of attaining $g_{\rm eff}\sim\Omega_{z}$ is that the electron mixes the
otherwise independent cavity modes and creates a localized state at a
frequency below the cavity resonance. We emphasize that this key feature is
largely insensitive to a choice of the lateral size $L$ and thus remains even
in the 2D thermodynamic limit $L/d\to\infty$.

Application to a 2D electron under the magnetic field.— As a proof-of-concept
demonstration, we now focus on a prototypical setting of ultrastrong-coupling
physics, namely, a 2D parabolic electron subject to a static perpendicular
magnetic field. From now on, we consider the $n=1$ principal hyperbolic mode
that most strongly couples to the cyclotron motion of the electron, and
abbreviate the subscript $n$ for the sake of notational simplicity. To
simplify the Hamiltonian, we choose the symmetric gauge,
$\boldsymbol{A}_{s}(\boldsymbol{r})=(-By/2,Bx/2)^{{\rm T}}$, with the magnetic
field $B$ and introduce the annihilation operator of Landau levels by

$\displaystyle\hat{a}=\frac{1}{\sqrt{2}}\left[\frac{l_{B}}{\hbar}\left(\hat{p}_{x}-i\hat{p}_{y}\right)-\frac{i}{2l_{B}}\left(\hat{x}-i\hat{y}\right)\right],$
(7)

where $l_{B}=\sqrt{\hbar/(eB)}$ is the magnetic length. Using the cyclotron
frequency $\omega_{c}=eB/m$ and the operator
$\hat{\boldsymbol{\pi}}=\frac{1}{\sqrt{2}}\left(\hat{a}+\hat{a}^{\dagger},i(\hat{a}-\hat{a}^{\dagger})\right)^{{\rm
T}}$, we obtain

$\displaystyle\hat{H}_{U}$ $\displaystyle\\!=\\!$
$\displaystyle\frac{\hbar\omega_{c}}{2}\left(\hat{\boldsymbol{\pi}}\\!+\\!\sum_{\boldsymbol{q}}\boldsymbol{e}_{\boldsymbol{q}}\left[c_{\boldsymbol{q}}\left(\hat{\gamma}_{\boldsymbol{q}}\\!+\\!\hat{\gamma}_{\boldsymbol{q}}^{\dagger}\right)\\!-\\!ql_{B}\hat{\gamma}_{\boldsymbol{q}}^{\dagger}\hat{\gamma}_{\boldsymbol{q}}\right]\right)^{2}\\!\\!+\\!\hat{H}_{\rm
pol},$

where
$c_{\boldsymbol{q}}=g_{\boldsymbol{q}}/\sqrt{\omega_{\boldsymbol{q}}\omega_{c}}$
is the dimensionless coefficient characterizing the coupling strength of the
cyclotron motion to the hyperbolic polaritons with momentum $\boldsymbol{q}$.
In the long-wavelength limit $ql_{B}\to 0$, the polariton-polariton
interaction in Eq. (Cavity Quantum Electrodynamics with Hyperbolic van der
Waals Materials) disappears and the problem reduces to the quadratic one. As
demonstrated below, however, this interaction term can in general contribute
to the dynamics and affect the absorption spectrum especially in the
ultrastrong coupling regimes. This is because the electron couples most
strongly with the electromagnetic modes at a finite momentum around $q\sim
q^{*}$ (cf. Fig. 2b) whose corresponding length scale $1/q^{*}$ is comparable
to the magnetic length $l_{B}$ near the cyclotron resonance. More generally,
such spatial dependence of the vector potential must be taken into account
when $1/q^{*}$ is comparable to the characteristic length scale of the
electron motion.

The low-energy excitations can be studied by analyzing the magnetoabsorption
spectrum

$\displaystyle A(\omega)={\rm Re}\left[\int_{0}^{\infty}e^{i\omega
t}\langle{\rm GS}|\hat{a}e^{-i\hat{H}t}\hat{a}^{\dagger}|{\rm
GS}\rangle\right],$ (9)

where $|{\rm GS}\rangle$ is the ground state. To reveal its qualitative
features, we perform a simple variational analysis as follows. Specifically,
we first determine the variational ground state of Eq. (Cavity Quantum
Electrodynamics with Hyperbolic van der Waals Materials) in the form of a
product of coherent states. We then expand the Hamiltonian around this state,
obtain the fluctuations up to the quadratic terms, and determine the
excitation spectrum via the exact diagonalization of the effective quadratic
Hamiltonian SM (1).

Figure 3 shows the obtained magnetoabsorption spectra, where the cavity
thickness $d$ is varied while the aspect ratio $L/d$ is kept constant. The
blue solid curve in each panel shows the bound-state energy of the Landau
polariton. As the thickness is decreased, the spectrum starts to exhibit the
anticrossings around the cyclotron resonance. In particular, when the cavity
length becomes a few nanometers, the large separations between the branches
that are comparable to the elementary excitation energies themselves emerge;
this is a hallmark of the ultrastrong coupling regime.

As discussed before, a key feature here is the formation of the dressed bound
state consisting of the electron and the localized phonon polaritons, which
manifests itself as the anticrossed lower branch in the spectrum (cf. the blue
curve in Fig. 3). Importantly, this feature remains independently of the
lateral size $L$, while the effect of increasing $L/d$ is the appearance of a
continuum of cavity modes above the lower branch SM (1). It is worthwhile to
note that the positions of the anticrossings are shifted above the bare
resonances $\omega_{c}/\Omega_{z}\simeq 1$. This upward shift originates from
the renormalization of the effective polariton energies due to the repulsive
polariton-polariton interaction. Since the latter comes from the spatial
dependence of the vector potential, this effect is absent in the long-
wavelength limit. We also remark that the appearance of the multiple
anticrossed branches in Fig. 3 are due to the discretized in-plane momentum
$q$, whose value is set by the periodic boundary conditions in the lateral
directions.

Discussions.— The proposed platform for ultrastrong-coupling cavity QED
materials can be realized by various polar vdW materials exhibiting hyperbolic
phonon polaritons, including Bi2Se3, Bi2Te3, MoS2 and MoO3 as well as h-BN
Esslinger  _et al._ (2014); Geim and Grigorieva (2013); Sun  _et al._ (2014);
Ma  _et al._ (2018); Zheng  _et al._ (2018, 2019). Thus, by confining
materials in the cavity consisting of these vdW structures, one can strongly
couple electronic excitations to quantized electromagnetic modes in a wide
spectral range from mid- or far-infrared to THz regimes. Moreover, the layered
nature of vdW crystals should readily allow one to tune the cavity coupling
strengths by controlling the thickness of the surrounding crystals.

While polaritons in these materials are known to exhibit ultralow loss
Caldwell  _et al._ (2014); Dai  _et al._ (2014); Giles  _et al._ (2018); Dai
_et al._ (2019); Ma  _et al._ (2022); Sheinfux  _et al._ (2022); Ma  _et al._
(2018); Zheng  _et al._ (2018, 2019), it merits further study to examine a
role of dissipation in the present cavity setup. In particular, it is
intriguing to explore if and when a confinement of quantum materials in the
cavity with dissipative hyperbolic polaritons could lead to phase transitions
or dynamics unique to open systems Ashida  _et al._ (2020b); Mivehvar  _et
al._ (2021). We also note that in general the contribution from the
instantaneous Coulomb energy could be important when the 2D material is placed
near the surfaces. In the present planar setup, this contribution might lead
to sizable effects in many-body regimes due to dielectric screening of the
electron-electron Coulomb interaction Li  _et al._ (2015).

In summary, we showed that the planar cavity consisting of vdW
heterostructures provides a promising platform to attain the single-electron
ultrastrong coupling in the THz or mid-infrared regions. As a proof-of-concept
demonstration, we presented the analysis of the magnetoabsorption spectrum for
the 2D electron confined in the h-BN cavity, where the coupling can be
ultrastrong provided that the thicknesses are judiciously controlled. We
expect that our results open a way to study wealth of recent predictions in
the emerging field of cavity QED materials. We hope that our work simulates
further studies in these directions.

###### Acknowledgements.

We are grateful to Iliya Esin, Ilya Esterlis, Tim Kaxiras, Igor Khanonkin,
Frank Koppens, Kanta Masuki, Gil Refael, Tao Shi, and Kenji Yasuda for
fruitful discussions. Y.A. acknowledges support from the Japan Society for the
Promotion of Science through Grant No. JP19K23424. A.I. was supported by the
Swiss National Science Foundation (SNSF) under Grant Number 200020_207520.
E.D. acknowledges support from the ARO grant “Control of Many-Body States
Using Strong Coherent Light-Matter Coupling in Terahertz Cavities” and the
Swiss National Science Foundation under Division II.

## References

  * DiCarlo  _et al._ (2009) L. DiCarlo, J. M. Chow, J. M. Gambetta, L. S. Bishop, B. R. Johnson, D. Schuster, J. Majer, A. Blais, L. Frunzio, S. Girvin, and R. J. Schoelkopf, Nature 460, 240 (2009).
  * Blais  _et al._ (2021) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Stute  _et al._ (2012) A. Stute, B. Casabone, P. Schindler, T. Monz, P. Schmidt, B. Brandstätter, T. Northup, and R. Blatt, Nature 485, 482 (2012).
  * Vaidya  _et al._ (2018) V. D. Vaidya, Y. Guo, R. M. Kroeze, K. E. Ballantine, A. J. Kollár, J. Keeling, and B. L. Lev, Phys. Rev. X 8, 011002 (2018).
  * Hagenmüller  _et al._ (2010) D. Hagenmüller, S. De Liberato, and C. Ciuti, Phys. Rev. B 81, 235303 (2010).
  * De Bernardis  _et al._ (2018) D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, Phys. Rev. A 98, 053819 (2018).
  * Le Boite (2020) A. Le Boite, Adv. Quant. Tech. 3, 1900140 (2020).
  * Ashida  _et al._ (2021) Y. Ashida, A. İmamoğlu, and E. Demler, Phys. Rev. Lett. 126, 153603 (2021).
  * Garcia-Vidal  _et al._ (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, Science 373, eabd0336 (2021).
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef, Appl. Phys. Rev. 9, 011312 (2022).
  * Bloch  _et al._ (2022) J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Nature 606, 41 (2022).
  * Devoret  _et al._ (2007) M. Devoret, S. Girvin, and R. Schoelkopf, Ann. Phys. 519, 767 (2007).
  * Forn-Díaz  _et al._ (2017) P. Forn-Díaz, J. J. García-Ripoll, B. Peropadre, J.-L. Orgiazzi, M. Yurtalan, R. Belyansky, C. M. Wilson, and A. Lupascu, Nat. Phys. 13, 39 (2017).
  * Yoshihara  _et al._ (2017) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Nat. Phys. 13, 44 (2017).
  * Forn-Díaz  _et al._ (2019) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * Frisk Kockum  _et al._ (2019) A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Nat. Rev. Phys. 1, 19 (2019).
  * Scalari  _et al._ (2014) G. Scalari, C. Maissen, S. Cibella, R. Leoni, P. Carelli, F. Valmorra, M. Beck, and J. Faist, New J. Phys. 16, 033005 (2014).
  * Smolka  _et al._ (2014) S. Smolka, W. Wuester, F. Haupt, S. Faelt, W. Wegscheider, and A. Imamoglu, Science 346, 332 (2014).
  * Maissen  _et al._ (2014) C. Maissen, G. Scalari, F. Valmorra, M. Beck, J. Faist, S. Cibella, R. Leoni, C. Reichl, C. Charpentier, and W. Wegscheider, Phys. Rev. B 90, 205309 (2014).
  * Zhang  _et al._ (2016) Q. Zhang, M. Lou, X. Li, J. L. Reno, W. Pan, J. D. Watson, M. J. Manfra, and J. Kono, Nature Phys 12, 1005 (2016).
  * Bayer  _et al._ (2017) A. Bayer, M. Pozimski, S. Schambeck, D. Schuh, R. Huber, D. Bougeard, and C. Lange, Nano Lett. 17, 6340 (2017).
  * 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, Phys. Rev. B 101, 075301 (2020).
  * Thomas  _et al._ (2021) A. Thomas, E. Devaux, K. Nagarajan, G. Rogez, M. Seidel, F. Richard, C. Genet, M. Drillon, and T. W. Ebbesen, Nano Lett. 21, 4365 (2021).
  * Flick  _et al._ (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, Proc. Natl. Acad. Sci. U.S.A. 114, 3026 (2017).
  * Rokaj  _et al._ (2019) V. Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, and A. Rubio, Phys. Rev. Lett. 123, 047202 (2019).
  * Juraschek  _et al._ (2021) D. M. Juraschek, T. Neuman, J. Flick, and P. Narang, Phys. Rev. Research 3, L032046 (2021).
  * Román-Roche and Zueco (2022) J. Román-Roche and D. Zueco, SciPost Phys. Lect. Notes , 50 (2022).
  * Schlawin  _et al._ (2019) F. Schlawin, A. Cavalleri, and D. Jaksch, Phys. Rev. Lett. 122, 133602 (2019).
  * Curtis  _et al._ (2019) J. B. Curtis, Z. M. Raines, A. A. Allocca, M. Hafezi, and V. M. Galitski, Phys. Rev. Lett. 122, 167002 (2019).
  * Sentef  _et al._ (2018) M. A. Sentef, M. Ruggenthaler, and A. Rubio, Sci. Adv. 4, eaau6969 (2018).
  * Gao  _et al._ (2021) H. Gao, F. Schlawin, and D. Jaksch, Phys. Rev. B 104, L140503 (2021).
  * Ashida  _et al._ (2020a) Y. Ashida, A. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, Phys. Rev. X 10, 041027 (2020a).
  * Latini  _et al._ (2021) S. Latini, D. Shin, S. A. Sato, C. Schäfer, U. De Giovannini, H. Hübener, and A. Rubio, Proc. Natl. Acad. Sci. U.S.A. 118, e2105618118 (2021).
  * Lenk  _et al._ (2022) K. Lenk, J. Li, P. Werner, and M. Eckstein, arXiv:2205.05559 (2022).
  * Paravicini-Bagliani  _et al._ (2019) G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, Nat. Phys. 15, 186 (2019).
  * 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, Science 375, 1030 (2022).
  * Rokaj  _et al._ (2022) V. Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, and A. Rubio, Phys. Rev. B 105, 205424 (2022).
  * Wang  _et al._ (2019) X. Wang, E. Ronca, and M. A. Sentef, Phys. Rev. B 99, 235156 (2019).
  * Masuki and Ashida (2022) K. Masuki and Y. Ashida, arXiv:2209.01363 (2022).
  * Caldwell  _et al._ (2014) J. D. Caldwell, A. V. Kretinin, Y. Chen, V. Giannini, M. M. Fogler, Y. Francescato, C. T. Ellis, J. G. Tischler, C. R. Woods, A. J. Giles, M. Hong, K. Watanabe, T. Taniguchi, S. A. Maier, and K. S. Novoselov, Nat. Commun. 5, 1 (2014).
  * Dai  _et al._ (2014) S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. McLeod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, Science 343, 1125 (2014).
  * Giles  _et al._ (2018) A. J. Giles, S. Dai, I. Vurgaftman, T. Hoffman, S. Liu, L. Lindsay, C. T. Ellis, N. Assefa, I. Chatzakis, T. L. Reinecke, J. G. Tischler, M. M. Fogler, J. H. Edgar, D. N. Basov, and J. D. Caldwell, Nat. Mater. 17, 134 (2018).
  * Dai  _et al._ (2019) S. Dai, W. Fang, N. Rivera, Y. Stehle, B.-Y. Jiang, J. Shen, R. Y. Tay, C. J. Ciccarino, Q. Ma, D. Rodan-Legrain, P. Jarillo-Herrero, E. H. T. Teo, M. M. Fogler, P. Narang, J. Kong, and D. N. Basov, Adv. Mater. 31, 1806603 (2019).
  * Ma  _et al._ (2022) E. Y. Ma, J. Hu, L. Waldecker, K. Watanabe, T. Taniguchi, F. Liu, and T. F. Heinz, Nano Lett. 22, 8389 (2022).
  * Sheinfux  _et al._ (2022) H. H. Sheinfux, L. Orsini, M. Jung, I. Torre, M. Ceccanti, R. Maniyara, D. B. Ruiz, A. Hötger, R. Bertini, S. Castilla, N. C. H. Hesp, E. Janzen, A. Holleitner, V. Pruneri, J. H. Edgar, G. Shvets, and F. H. L. Koppens, arXiv:2202.08611 (2022).
  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013).
  * Cohen-Tannoudji  _et al._ (1989) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Photons and Atoms_ (Wiley, New York, 1989).
  * SM (1) See Supplemental Materials for further details on the statements and the derivations.
  * Lee  _et al._ (1953) T. D. Lee, F. E. Low, and D. Pines, Phys. Rev. 90, 297 (1953).
  * Ashida  _et al._ (2022) Y. Ashida, T. Yokota, A. İmamoğlu, and E. Demler, Phys. Rev. Research 4, 023194 (2022).
  * Esslinger  _et al._ (2014) M. Esslinger, R. Vogelgesang, N. Talebi, W. Khunsin, P. Gehring, S. De Zuani, B. Gompf, and K. Kern, ACS Photonics 1, 1285 (2014).
  * Sun  _et al._ (2014) J. Sun, N. M. Litchinitser, and J. Zhou, ACS Photon. 1, 293 (2014).
  * Ma  _et al._ (2018) W. Ma, P. Alonso-González, S. Li, A. Y. Nikitin, J. Yuan, J. Martín-Sánchez, J. Taboada-Gutiérrez, I. Amenabar, P. Li, S. Vélez, _et al._ , Nature 562, 557 (2018).
  * Zheng  _et al._ (2018) Z. Zheng, J. Chen, Y. Wang, X. Wang, X. Chen, P. Liu, J. Xu, W. Xie, H. Chen, S. Deng, and N. Xu, Adv. Mater. 30, 1705318 (2018).
  * Zheng  _et al._ (2019) Z. Zheng, N. Xu, S. L. Oscurato, M. Tamagnone, F. Sun, Y. Jiang, Y. Ke, J. Chen, W. Huang, W. L. Wilson, A. Ambrosio, S. Deng, and H. Chen, Sci. Adv. 5, eaav8690 (2019).
  * Ashida  _et al._ (2020b) Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020b).
  * Mivehvar  _et al._ (2021) F. Mivehvar, F. Piazza, T. Donner, and H. Ritsch, Adv. Phys. 70, 1 (2021).
  * Li  _et al._ (2015) L. H. Li, E. J. Santos, T. Xing, E. Cappelluti, R. Roldán, Y. Chen, K. Watanabe, and T. Taniguchi, Nano Lett. 15, 218 (2015).

Supplementary Materials

##  .1 Derivation of the effective cavity QED Hamiltonian

We here derive the effective Hamiltonian of the proposed cavity QED system
discussed in the main text. To this end, we start from a general microscopic
model of the 2D electron coupled to quantized electromagnetic fields in
uniaxial polar dielectrics in the Coulomb gauge Cohen-Tannoudji  _et al._
(1989); Ashida  _et al._ (2020a):

$\displaystyle\hat{H}=\frac{\left[\hat{\boldsymbol{p}}+e\boldsymbol{A}_{s}(\hat{x},\hat{y})+e\hat{\boldsymbol{A}}(\hat{x},\hat{y},z=0)\right]^{2}}{2m}+\hat{H}_{0},\;\;\;\hat{\boldsymbol{p}}=\left(\begin{array}[]{c}\hat{p}_{x}\\\
\hat{p}_{y}\end{array}\right),$ (S3)

where $m$ is the electron mass, $e$ is the elementary charge, and
$\boldsymbol{A}_{s}$ is an arbitrary static external field. We note that the
2D vector field $\hat{\boldsymbol{A}}(x,y,z=0)$ is defined by the projection
of the following 3D quantized vector potential onto the 2D tangential plane at
$z=0$:

$\displaystyle\hat{\boldsymbol{A}}(\boldsymbol{r})=\sum_{\boldsymbol{k}\lambda}\sqrt{\frac{\hbar}{2V\epsilon_{0}ck}}\left[\hat{a}_{\boldsymbol{k}\lambda}\boldsymbol{\epsilon}_{\boldsymbol{k}\lambda}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}+{\rm
H.c.}\right],$ (S4)

where $V$ is the volume of the entire space, and
$\hat{a}_{\boldsymbol{k}\lambda}$
($\hat{a}_{\boldsymbol{k}\lambda}^{\dagger}$) is the annihilation (creation)
operator of photons with wavevector $\boldsymbol{k}$ and polarization
$\lambda$; note that these excitations correspond to bare photons defined in
the absence of dielectric materials. In the Coulomb gauge, the orthonormal
transverse polarization vectors
$\boldsymbol{\epsilon}_{\boldsymbol{k}\lambda}$ satisfy
$\boldsymbol{k}\cdot\boldsymbol{\epsilon}_{\boldsymbol{k}\lambda}=0$ and
$\boldsymbol{\epsilon}_{\boldsymbol{k}\lambda}^{\dagger}\boldsymbol{\epsilon}_{\boldsymbol{k}\nu}=\delta_{\lambda\nu}$.
The term $\hat{H}_{0}$ is the quadratic Hamiltonian consisting of the energies
of the quantized electromagnetic fields and the polar phonons in dielectrics:

$\displaystyle\hat{H}_{0}=\int
d^{3}r\left[\frac{\hat{\boldsymbol{\Pi}}^{2}}{2\boldsymbol{\epsilon}(z)}+\frac{\epsilon_{0}c^{2}}{2}\left(\nabla\times\hat{\boldsymbol{A}}\right)^{2}\right]+\int_{{\rm
D}}\frac{d^{3}r}{v}\sum_{i\in\\{x,y,z\\}}\left[\frac{\left(\hat{\pi}_{i}-Z^{*}e\hat{A}_{i}\right)^{2}}{2M_{i}}+\frac{1}{2}M_{i}\Omega_{i}^{2}\hat{\phi}_{i}^{2}+\frac{\left(Z^{*}e\hat{\phi}_{i}^{\parallel}\right)^{2}}{2\epsilon_{i}(z)v}\right],$
(S5)

where $\hat{\boldsymbol{\Pi}}(\boldsymbol{r})$ is the canonically conjugate
vector field of $\hat{\boldsymbol{A}}(\boldsymbol{r})$, $v$ is the unit-cell
volume of the dielectrics, $\hat{\boldsymbol{\phi}}(\boldsymbol{r})$ is a
vector field representing polar optical phonons,
$\hat{\boldsymbol{\pi}}(\boldsymbol{r})$ is its canonically conjugate vector
field, $Z^{*}e$ is the effective charge of polar phonons, and $\int_{{\rm
D}}d^{3}r$ represents the integration over the dielectric regions, $\int_{{\rm
D}}d^{3}r=\int
d^{2}r\left(\int_{-d/2}^{-\delta/2}+\int_{\delta/2}^{d/2}\right)dz$, with
$\delta$ being the thickness of the narrow air gap that is assumed to be much
shorter than the cavity thickness $d$. We also define the permittivity tensor
$\boldsymbol{\epsilon}(z)$ by

$\displaystyle\boldsymbol{\epsilon}(z)={\rm
diag}\left(\epsilon_{t\infty},\epsilon_{t\infty},\epsilon_{z\infty}\right)\vartheta\left(\frac{d}{2}-|z|\right)\vartheta\left(|z|-\frac{\delta}{2}\right)+\epsilon_{0}I\left[\vartheta\left(|z|-\frac{d}{2}\right)+\vartheta\left(\frac{\delta}{2}-|z|\right)\right],$
(S6)

where $\epsilon_{t(z)\infty}$ is the in-plane (out-of-plane) permittivity in
the infinite-frequency limit and $I$ is the identity matrix. We denote the in-
plane masses and phonon frequencies by $M_{x,y}=M_{t}$ and
$\Omega_{x,y}=\Omega_{t}$, respectively. We remark that the permittivity in
the air gap region can in general take a value different from the vacuum value
depending on a specific choice of the cavity geometry, while we take it to be
$\epsilon_{0}$ in the present work for the sake of simplicity.

To derive the effective Hamiltonian, we first need to identify a polariton
eigenmode that diagonalizes the quadratic part $\hat{H}_{0}$. This can be done
by solving the macroscopic Maxwell’s equations derived from the above model.
We remark that the influence of the static external field can be neglected
when constructing such polariton modes. Without loss of generality, we
consider eigenmodes propagating along the $x$ direction with in-plane momentum
$\boldsymbol{q}=(q,0,0)^{{\rm T}}$. Also, we recall that we are interested in
the hyperbolic eigenmodes that strongly couple to the 2D electron moving on
the $z=0$ plane. Such modes correspond to the P-polarized (or equivalently,
the transverse magnetic) modes with the even parity for the tangential
component, which satisfy the following conditions: $[\boldsymbol{B}]_{y}\neq
0,[\boldsymbol{B}]_{x,z}=0$,
$[\boldsymbol{E}]_{y}=0,[\boldsymbol{E}]_{x,z}\neq 0$, and
$E_{x}(z)=E_{x}(-z),E_{z}(z)=-E_{z}(-z)$. The resulting eigenequations that
characterize these modes are

$\displaystyle\frac{q^{2}}{\epsilon_{z}(\omega)}+\frac{\kappa^{2}}{\epsilon_{t}(\omega)}$
$\displaystyle=$ $\displaystyle\frac{\omega^{2}}{\epsilon_{0}c^{2}},$ (S7)
$\displaystyle\tan\left(\frac{\kappa d}{2}\right)$ $\displaystyle=$
$\displaystyle-\frac{\kappa}{\epsilon_{t}\sqrt{q^{2}-\frac{\omega^{2}}{c^{2}}}},$
(S8)

where $\kappa$ is the out-of-plane momentum and the in-plane and out-of-plane
permittivities $\epsilon_{t,z}(\omega)$ are given by

$\displaystyle\epsilon_{t}(\omega)=\epsilon_{t\infty}\left(1+\frac{g_{t}^{2}}{\Omega_{t}^{2}-\omega^{2}}\right),\;\;\epsilon_{z}(\omega)=\epsilon_{z\infty}\left(1+\frac{g_{z}^{2}}{\Omega_{z}^{2}-\omega^{2}}\right).$
(S9)

Here, we denote the coupling strengths of the dielectrics by

$\displaystyle
g_{t}=\sqrt{\frac{\left(Z^{*}e\right)^{2}}{\epsilon_{t\infty}M_{t}v}},\;\;\;g_{z}=\sqrt{\frac{\left(Z^{*}e\right)^{2}}{\epsilon_{z\infty}M_{z}v}}.$
(S10)

For the h-BN crystals considered in the main text, we use the parameters in
Ref. Caldwell  _et al._ (2014) where the out-of-plane and in-plane phonon
frequencies are $\Omega_{z}=23.3\;{\rm THz}$ and $\Omega_{t}=41.1\;{\rm THz}$,
the corresponding couplings are $g_{z}=0.37\Omega_{z}$ and
$g_{t}=0.61\Omega_{t}$, and the out-of-plane and in-plane infinite-frequency
permittivities are $\epsilon_{z\infty}=2.95$ and $\epsilon_{t\infty}=4.87$,
respectively.

Figure S1:  Mode profiles of the h-BN Type I hyperbolic phonon polaritons. The
top, middle, and bottom panels plot the spatial dependences of the mode
functions along the thickness direction for their in-plane components, out-of-
plane components, and root mean squares, respectively. The red solid (blue
dotted) curves correspond to the electric-field (phonon) mode functions
$\boldsymbol{u}_{\boldsymbol{q}n}$
($\boldsymbol{u}_{\boldsymbol{q}n}^{\phi}$).

Solving the eigenequations (S7) and (S8) with respect to $\kappa$ and $\omega$
for a given $q$, we obtain a hyperbolic dispersion $\omega_{\boldsymbol{q}n}$
associated with the out-of-plane momentum $\kappa_{\boldsymbol{q}n}$, which
are labeled by the discrete number $n$. The corresponding mode function,
$\boldsymbol{u}_{\boldsymbol{q}n}$, characterizes the spatial profile of the
electric field and is given by

$\displaystyle\boldsymbol{u}_{\boldsymbol{q}n}$ $\displaystyle=$
$\displaystyle
f_{\boldsymbol{q}n}\left[\cos(\kappa_{\boldsymbol{q}n}z)\boldsymbol{e}_{x}-\frac{i\epsilon_{t}q}{\epsilon_{z}\kappa_{\boldsymbol{q}n}}\sin(\kappa_{\boldsymbol{q}n}z)\boldsymbol{e}_{z}\right]e^{iqx}\vartheta\left(\frac{d}{2}-|z|\right)\vartheta\left(|z|-\frac{\delta}{2}\right)$
(S11) $\displaystyle+$ $\displaystyle f^{\prime}_{\boldsymbol{q}n}\sum_{s=\pm
1}\left[\boldsymbol{e}_{x}+is(q/\nu_{\boldsymbol{q}n})\boldsymbol{e}_{z}\right]e^{iqx-s\nu_{\boldsymbol{q}n}z}\vartheta\left(|z|-\frac{d}{2}\right)$
$\displaystyle+$ $\displaystyle
f^{\prime\prime}_{\boldsymbol{q}n}\left[\cosh(\nu_{\boldsymbol{q}n}z)\boldsymbol{e}_{x}-i(q/\nu_{\boldsymbol{q}n})\sinh(\nu_{\boldsymbol{q}n}z)\boldsymbol{e}_{z}\right]e^{iqx}\vartheta\left(\frac{\delta}{2}-|z|\right),$

where
$\nu_{\boldsymbol{q}n}=\sqrt{q^{2}-{\omega_{\boldsymbol{q}n}^{2}}/{c^{2}}}$
and the amplitude coefficients satisfy the following conditions:

$\displaystyle f_{\boldsymbol{q}n}/f^{\prime}_{\boldsymbol{q}n}$
$\displaystyle=$ $\displaystyle
e^{-\nu_{\boldsymbol{q}n}d/2}/\cos(\kappa_{\boldsymbol{q}n}d/2),$ (S12)
$\displaystyle f_{\boldsymbol{q}n}/f^{\prime\prime}_{\boldsymbol{q}n}$
$\displaystyle=$
$\displaystyle\cosh\left(\nu_{\boldsymbol{q}n}\delta/2\right)/\cos(\kappa_{\boldsymbol{q}n}\delta/2).$
(S13)

Without loss of generality, we choose $f_{\boldsymbol{q}n}>0$. The
corresponding phonon mode function is given by

$\displaystyle\boldsymbol{u}_{\boldsymbol{q}n}^{\phi}=f_{\boldsymbol{q}n}\left[\frac{g_{t}\Omega_{t}}{\Omega_{t}^{2}-\omega_{\boldsymbol{q}n}^{2}}\cos(\kappa_{\boldsymbol{q}n}z)\boldsymbol{e}_{x}-\frac{i\epsilon_{t}q}{\epsilon_{z}\kappa_{\boldsymbol{q}n}}\frac{g_{z}\Omega_{z}}{\Omega_{z}^{2}-\omega_{\boldsymbol{q}n}^{2}}\sin(\kappa_{\boldsymbol{q}n}z)\boldsymbol{e}_{z}\right]e^{iqx}\vartheta\left(\frac{d}{2}-|z|\right)\vartheta\left(|z|-\frac{\delta}{2}\right).$
(S14)

We then impose the normalization condition,

$\displaystyle\int_{-\infty}^{\infty}dz\left[\left|\boldsymbol{u}_{\boldsymbol{q}n}(z)\right|^{2}+\left|\boldsymbol{u}_{\boldsymbol{q}n}^{\phi}(z)\right|^{2}\right]=1,$
(S15)

which fixes the amplitude coefficients.

The results for the three lowest branches of the h-BN Type I hyperbolic modes
are shown in Fig. S1. We remark that, on the $z=0$ plane where the 2D material
is inserted, the electric fields have only the in-plane components along the
propagation direction. We also note that, in the air gap region, the electric
field is written as the curl of the magnetic field and thus it is purely
transverse, i.e.,
$\boldsymbol{u}_{\boldsymbol{q}n}\propto\nabla\times\left(\sinh(\nu_{\boldsymbol{q}n}z)e^{iqx}\boldsymbol{e}_{y}\right)$
for $|z|<\frac{\delta}{2}$. Since we assume the narrow gap
$\nu,\kappa\ll\delta^{-1}$, we can introduce the effective confinement length
$d_{\boldsymbol{q}n}$ through the relation
$|\boldsymbol{u}_{\boldsymbol{q}n}(z=0)|\simeq f_{\boldsymbol{q}n}\equiv
1/\sqrt{d_{\boldsymbol{q}n}}$, which characterizes the coupling strength
$g_{\boldsymbol{q}n}$ of the 2D electron on the $z=0$ plane as discussed in
the main text.

##  .2 Details about the mean-field analysis of the magnetoabsorption spectrum

Figure S2:  Magnetoabsorption spectra $A(\omega)$ at different aspect ratios
$L/d$. The cavity thickness $d$ is fixed to be $d=3{\rm\,nm}$. The in-plane
momentum $\boldsymbol{q}$ is discretized as in Eq. (S29). The results are
plotted for the h-BN parameters and the in-plane momentum cutoff
$\Lambda=2{\,{\rm nm}^{-1}}$.

We here provide technical details about the mean-field variational analysis of
the magnetoabsorption spectrum discussed in the main text. To this end, we
first perform an additional unitary transformation,

$\displaystyle\hat{V}=e^{\sum_{\boldsymbol{q}n}\frac{c_{\boldsymbol{q}n}}{l_{B}}\left(\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}-\hat{\gamma}_{\boldsymbol{q}n}\right)},\;\;\hat{V}^{\dagger}\hat{\gamma}_{\boldsymbol{q}n}\hat{V}=\hat{\gamma}_{\boldsymbol{q}n}+\frac{c_{\boldsymbol{q}n}}{l_{B}},$
(S16)

which displaces the polariton operators. Using the relation
$\sum_{\boldsymbol{q}n}\boldsymbol{q}c_{\boldsymbol{q}n}^{2}=\boldsymbol{0}$,
the Hamiltonian in Eq. (8) in the main text is transformed to

$\displaystyle\hat{H}_{VU}=\hat{V}^{\dagger}\hat{H}_{U}\hat{V}=\frac{\hbar\omega_{c}}{2}\left(\hat{\boldsymbol{\pi}}-\sum_{\boldsymbol{q}n}\boldsymbol{q}l_{B}\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}\hat{\gamma}_{\boldsymbol{q}n}\right)^{2}+\sum_{\boldsymbol{q}n}\hbar\omega_{\boldsymbol{q}n}\left[\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}\hat{\gamma}_{\boldsymbol{q}n}+\frac{c_{\boldsymbol{q}n}}{l_{B}}\left(\hat{\gamma}_{\boldsymbol{q}n}+\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}\right)\right].$
(S17)

It is useful to rewrite it in terms of the canonically conjugate variables as
follows:

$\displaystyle\frac{\hat{H}_{VU}}{\hbar\omega_{c}}=\frac{\hat{X}^{2}+\hat{P}^{2}}{2}+\sum_{\boldsymbol{q}n}\epsilon_{\boldsymbol{q}n}\frac{\hat{X}_{\boldsymbol{q}n}^{2}+\hat{P}_{\boldsymbol{q}n}^{2}}{2}-\hat{\boldsymbol{\pi}}\cdot\hat{\boldsymbol{\Xi}}_{b}+\frac{1}{2}:\hat{\boldsymbol{\Xi}}_{b}^{2}:+\sum_{\boldsymbol{q}n}\xi_{\boldsymbol{q}n}\hat{X}_{\boldsymbol{q}n},$
(S18)

where we denote the conjugate variables associated with the cyclotron motion
and the polaritons by

$\displaystyle\hat{\boldsymbol{\pi}}=\left(\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}},\frac{i(\hat{a}-\hat{a}^{\dagger})}{\sqrt{2}}\right)^{{\rm
T}}=\left(\hat{X},\hat{P}\right)^{{\rm
T}},\;\;\left(\frac{\hat{\gamma}_{\boldsymbol{q}n}+\hat{\gamma}_{\boldsymbol{q}n}^{\dagger}}{\sqrt{2}},\frac{i(\hat{\gamma}_{\boldsymbol{q}n}-\hat{\gamma}_{\boldsymbol{q}n}^{\dagger})}{\sqrt{2}}\right)^{{\rm
T}}=\left(\hat{X}_{\boldsymbol{q}n},\hat{P}_{\boldsymbol{q}n}\right)^{{\rm
T}},$ (S19)

and introduce

$\displaystyle\hat{\boldsymbol{\Xi}}_{b}=\sum_{\boldsymbol{q}n}\frac{l_{B}\boldsymbol{q}}{2}\left(\hat{X}_{\boldsymbol{q}n}^{2}+\hat{P}_{\boldsymbol{q}n}^{2}\right),\;\;\epsilon_{\boldsymbol{q}n}=\frac{\omega_{\boldsymbol{q}n}}{\omega_{c}}+\frac{l_{B}^{2}q^{2}}{2},\;\;\xi_{\boldsymbol{q}n}=\sqrt{\frac{2\omega_{\boldsymbol{q}n}}{\omega_{c}^{3}}}\frac{g_{\boldsymbol{q}n}}{ql_{B}},$
(S20)

which correspond to the dimensionless total polariton momentum, the normalized
effective polariton energies, and the dimensionless coupling strengths,
respectively. The mean-field energy $\overline{E}$ is given by taking the
expectation value of Eq. (S18) with respect to the product of coherent states
for these conjugate variables. The mean-field ground state is then determined
by finding the zero of the derivatives of $\overline{E}$, leading to

$\displaystyle\overline{\boldsymbol{\pi}}=\overline{\boldsymbol{\Xi}}_{b}=\sum_{\boldsymbol{q}n}\frac{l_{B}\boldsymbol{q}}{2}\left(\overline{X}_{\boldsymbol{q}n}^{2}+\overline{P}_{\boldsymbol{q}n}^{2}\right),\;\;\overline{X}_{\boldsymbol{q}n}=\frac{-\xi_{\boldsymbol{q}n}}{\epsilon_{\boldsymbol{q}n}},\;\;\overline{P}_{\boldsymbol{q}n}=0.$
(S21)

Figure S3:  Coupling strengths $g_{q}$ at (a) different cavity thicknesses $d$
with $L/d=\pi$ and (b) different aspect ratios $L/d$ with $d=3{{\rm\,nm}}$.
The values corresponding to each of the discrete in-plane momenta
$\boldsymbol{q}$ of hyperbolic modes are indicated by the filled red circles.
The coupling strengths indicated in panel (a) and (b) correspond to Fig. 3 in
the main text and Fig. S2, respectively. The results are plotted for the h-BN
parameters and the in-plane momentum cutoff $\Lambda=2{\,{\rm nm}^{-1}}$.

To obtain the low-energy excitation spectrum, we perform the projection
$\hat{P}$ onto the quadratic fluctuations on top of this mean-field ground
state. Specifically, we describe the low-energy excitations in terms of the
effective quadratic Hamiltonian given by

$\displaystyle\hat{H}_{{\rm
eff}}=\hat{P}\hat{D}^{\dagger}\hat{H}_{VU}\hat{D}\hat{P},$ (S22)

where $\hat{D}$ induces the mean-field displacements

$\displaystyle\hat{D}^{\dagger}\hat{\boldsymbol{\pi}}\hat{D}=\hat{\boldsymbol{\pi}}+\overline{\boldsymbol{\pi}},\;\;\hat{D}^{\dagger}\hat{X}_{\boldsymbol{q}n}\hat{D}=\hat{X}_{\boldsymbol{q}n}+\overline{X}_{\boldsymbol{q}n},\;\;\hat{D}^{\dagger}\hat{P}_{\boldsymbol{q}n}\hat{D}=\hat{P}_{\boldsymbol{q}n}.$
(S23)

The result is

$\displaystyle\frac{\hat{H}_{{\rm eff}}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle\frac{\hat{X}^{2}+\hat{P}^{2}}{2}+\sum_{\boldsymbol{q}n}\tilde{\epsilon}_{\boldsymbol{q}n}\frac{\hat{X}_{\boldsymbol{q}n}^{2}+\hat{P}_{\boldsymbol{q}n}^{2}}{2}+\frac{1}{2}\left[\sum_{\boldsymbol{q}n}\tilde{\xi}_{\boldsymbol{q}n}\boldsymbol{e}_{\boldsymbol{q}}\hat{X}_{\boldsymbol{q}n}\right]^{2}-\hat{\boldsymbol{\pi}}\cdot\sum_{\boldsymbol{q}n}\tilde{\xi}_{\boldsymbol{q}n}\boldsymbol{e}_{\boldsymbol{q}}\hat{X}_{\boldsymbol{q}n}$
(S24) $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\hat{\boldsymbol{\phi}}^{{\rm
T}}M\hat{\boldsymbol{\phi}},$

where we introduce the modified effective polariton energies and the effective
coupling strengths by

$\displaystyle\tilde{\epsilon}_{\boldsymbol{q}n}\equiv\frac{\omega_{\boldsymbol{q}n}}{\omega_{c}}+\frac{l_{B}^{2}q^{2}}{2}-\frac{1}{2}l_{B}\boldsymbol{q}\cdot\overline{\boldsymbol{\Xi}}_{b},\;\;\;\tilde{\xi}_{\boldsymbol{q}n}\equiv\frac{g_{\boldsymbol{q}n}}{\omega_{c}}\sqrt{\frac{2\omega_{\boldsymbol{q}n}}{\omega_{c}}}\frac{1}{\frac{\omega_{\boldsymbol{q}n}}{\omega_{c}}+\frac{l_{B}^{2}q^{2}}{2}},$
(S25)

respectively. In the last line, we rewrite the Hamiltonian by using the matrix
$M$ and the vector of the conjugate variables

$\displaystyle\hat{\boldsymbol{\phi}}=\left(\hat{X},\ldots,\hat{X}_{\boldsymbol{q}n},\ldots,\hat{P},\ldots,\hat{P}_{\boldsymbol{q}n},\ldots\right)^{{\rm
T}}.$ (S26)

The low-energy excitation spectrum $\left\\{\omega_{\lambda}\right\\}$ can
then be obtained from the Williamson eigenvalues of $M$ as follows:

$\displaystyle S^{{\rm T}}MS={\rm
diag}(\left\\{\omega_{\lambda}/\omega_{c}\right\\},\left\\{\omega_{\lambda}/\omega_{c}\right\\}),$
(S27)

where $S$ is a symplectic matrix. We note that this treatment is exact in the
long-wavelength limit $ql_{B}\to 0$ in which Eq. (S18) simplifies to the
quadratic Hamiltonian.

Finally, the absorption spectrum can be obtained from

$\displaystyle A(\omega)$ $\displaystyle\simeq$ $\displaystyle{\rm
Re}\left[\int_{0}^{\infty}e^{i\omega t}\langle
0|\hat{D}^{\dagger}\hat{a}e^{-i\hat{H}_{VU}t}\hat{a}^{\dagger}\hat{D}|0\rangle\right]$
(S28) $\displaystyle\simeq$ $\displaystyle{\rm
Re}\left[\int_{0}^{\infty}e^{i\omega t}\langle
0|\left(\hat{a}+\frac{\overline{\pi}_{x}-i\overline{\pi}_{y}}{\sqrt{2}}\right)e^{-i\hat{H}_{{\rm
eff}}t}\left(\hat{a}^{\dagger}+\frac{\overline{\pi}_{x}+i\overline{\pi}_{y}}{\sqrt{2}}\right)|0\rangle\right].$

The last line can easily be evaluated by using the diagonalized basis in Eq.
(S27). We note that the results in the main text are obtained for the
discretized in-plane wavenumbers as follows:

$\displaystyle
q_{i}\in\left\\{-\Lambda,-\frac{N-1}{N}\Lambda,\ldots,-\frac{\Lambda}{N},0,\frac{\Lambda}{N},\ldots,\frac{N-1}{N}\Lambda,\Lambda\right\\},\;\;i\in\\{x,y\\}.$
(S29)

Here, the integer number $N$ is related to the lateral system size $L$ via
$\left[\frac{\Lambda L}{2\pi}\right]=N$. The results obtained for the cavity
consisting of ultrathin h-BN materials with different aspect ratios $L/d$ are
plotted in Fig. S2. The coupling strengths to each hyperbolic phonon polariton
mode with discretized in-plane momentum are plotted for different thicknesses
and aspect ratios in Fig. S3. When the lateral size $L$ is increased, the
coupling strength $g_{q}$ to each mode decreases uniformly as $g_{q}\propto
L^{-1}$ while there appear a larger number of cavity modes. Consequently, a
key feature of the ultrastrong coupling regime in the spectrum, that is, the
formation of the localized Landau-polariton mode remains almost the same
independently of the lateral size $L$. Meanwhile, the increase of $L/d$ leads
to the appearance of dense anticrossed branches originating from the
hybridization with the continuum cavity modes above the lower Landau-polariton
mode.