# Thermodynamics of correlated electrons in a magnetic field

Jixun K. Ding<EMAIL_ADDRESS>Department of Applied Physics, Stanford
University, CA 94305, USA  Wen O. Wang  Department of Applied Physics,
Stanford University, CA 94305, USA  Brian Moritz  Stanford Institute for
Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand
Hill Road, Menlo Park, CA 94025, USA  Yoni Schattner  Stanford Institute for
Materials and Energy Sciences, SLAC National Accelerator Laboratory, 2575 Sand
Hill Road, Menlo Park, CA 94025, USA  Department of Physics, Stanford
University, Stanford, CA 94305, USA  Edwin W. Huang  Department of Physics and
Institute of Condensed Matter Theory, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  Thomas P. Devereaux<EMAIL_ADDRESS>Stanford Institute for Materials and Energy Sciences, SLAC National
Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA
Department of Materials Science and Engineering, Stanford University,
Stanford, CA 94305, USA

###### Abstract

The Hofstadter-Hubbard model captures the physics of strongly correlated
electrons in an applied magnetic field, which is relevant to many recent
experiments on Moiré materials. Few large-scale, numerically exact simulations
exists for this model. In this work, we simulate the Hubbard-Hofstadter model
using the determinant quantum Monte Carlo (DQMC) algorithm. We report the
field and Hubbard interaction strength dependence of charge compressibility,
fermion sign, local moment, magnetic structure factor, and specific heat. The
gross structure of magnetic Bloch bands and band gaps determined by the non-
interacting Hofstadter spectrum is preserved in the presence of $U$.
Incompressible regions of the phase diagram have improved fermion sign. At
half filling and intermediate and larger couplings, a strong orbital magnetic
field delocalizes electrons and reduces the effect of Hubbard $U$ on
thermodynamic properties of the system.

##  I Introduction

Strong magnetic fields allow us to probe the phase diagram of strongly
correlated materials and uncover novel phases. For example, a magnetic field
induces charge/pair density wave order in cuprate superconductors [1, 2, 3,
4]; magnetic field is a convenient tuning parameter for accessing quantum
critical points [5, 6]; and field-induced reentrant superconductivity also has
been reported in uranium compounds [7, 8]. With the recent proliferation of
experimental evidence for fractional quantum Hall effect, superconductivity,
and other correlated electron phases [9, 10, 11, 12, 13, 14] in graphene Moiré
superlattices, there is renewed interest in studying the behavior of strongly
correlated electronic systems in strong magnetic fields.

Properties of non-interacting electrons in a two-dimensional periodic lattice
under the influence of a strong magnetic field are fairly well-understood. In
this system, the competition between lattice and magnetic length scales leads
to the fractal Hofstadter butterfly spectrum with recursive magnetic subband
structure [15, 16], which generalizes the idea of Landau levels in a free
electron gas. The Chern numbers associated with these magnetic subbands
provide an elegant explanation of the integer quantum Hall effect [17].
Experimentally, the Hofstadter Hamiltonian has been realized in ultra-cold
atoms loaded on optical lattices [18], and direct observation of the
Hofstadter spectrum has been reported in Moiré superlattices in graphene with
high resolution [19, 20].

The most natural framework for understanding the simultaneous influence of
magnetic field and Coulomb interaction on electrons in a periodic lattice is
to take the Hofstadter Hamiltonian and add to it a Hubbard interaction term.
In the literature, this is sometimes called the Hofstadter-Hubbard or Hubbard-
Hofstadter model. This model has been investigated using Hartree-Fock mean-
field theory [21, 22], exact diagonalization [23, 24], dynamical mean-field
theory [25, 26], and in the large $U$ limit via renormalized mean-field theory
[27]. Aside from exact diagonalization, which is limited to small system
sizes, all methods used to study the Hubbard-Hofstadter model have been
approximate and don’t capture the full extent of quantum fluctuations. It is
not conclusive, for example, whether interactions change or preserve the gap
structure of the Hofstadter butterfly [22, 24, 26].

Determinant quantum Monte Carlo (DQMC) [28, 29, 30] is an unbiased and
numerically exact algorithm for studying quantum systems at finite
temperature. It employs a discrete Hubbard–Stratonovich transformation to
reduce the quartic Hubbard interaction term to quadratic at the cost of
introducing a fluctuating auxiliary field. This auxiliary field is then
sampled using the Metropolis-Hastings algorithm. DQMC has been employed
successfully in the (zero-field) Hubbard model to study spin and charge
excitations [31, 32] and superconducting fluctuations [33], as well as find
evidence for fluctuating stripes [34] and $T$-linear resistivity [35]. The
DQMC method is especially powerful at half-filling in the absence of kinetic
frustration [36], where the fermion sign problem is absent due to particle-
hole symmetry, even in the presence of a magnetic field. This allows
simulations to be performed at much lower temperatures, providing access to
properties more reflective of the ground state.

At half filling, the Fermi surface of the non-interacting Hofstadter model
consists of a finite number of Dirac points at even-denominator rational
fractions of magnetic flux per plaquette [37]. The ground state of the
Hubbard-Hofstadter model is thus expected to remain a Dirac semi-metal up to
some finite coupling strength $U_{\mathrm{c}}$. At half a magnetic flux
quantum per plaquette, the model also is known as the $\mathrm{\pi}$-flux
model [38], and has been studied extensively numerically. As interactions are
turned on, the $\mathrm{\pi}$-flux model exhibits a quantum phase transition
of the chiral Heisenberg Gross-Neveu universality class at
$U_{\mathrm{c}}\approx 5.6t$ into an antiferromagnetic Mott insulator (AFMI)
[39, 40, 41, 42, 43]. Since the $\mathrm{\pi}$-flux model corresponds to the
Hubbard-Hofstadter Hamiltonian threaded with maximum possible flux, we may
think of the zero-field Hubbard model on a half-filled square lattice as the
“$0$-flux model”, which exibits a metal$-$AFMI transition with
$U_{\mathrm{c}}=0$ [29, 30]. Our simulations address intermediate field
strengths between the $0$-flux and $\mathrm{\pi}$-flux Hubbard model, which,
to the best of our knowledge, has not been studied via DQMC.

In this work, we study the Hubbard-Hofstadter model using DQMC and present the
evolution of thermodynamic properties of correlated electrons in an orbital
magnetic field $B$. We demonstrate that the gross structure of magnetic Bloch
bands and band gaps determined by the non-interacting Hofstadter spectrum is
preserved in the presence of $U$. Moreover, we determine that the many-body
fermion sign is directly connected to electronic charge compressibility.
Finally, focusing on the half-filled AFMI, we find that an orbital magnetic
field tends to delocalize electrons and thus effectively lower the influence
of Hubbard $U$.

##  II Methods

We study the single-band Hubbard-Hofstadter model on a two dimensional square
lattice

$H=-t\sum_{\langle
ij\rangle\sigma}\left\\{\exp\left[\mathrm{i}\varphi_{ij}\right]c_{i\sigma}^{\dagger}c_{j\sigma}+\mathrm{h.c.}\right\\}\\\
-\mu\sum_{i\sigma}n_{i\sigma}+U\sum_{i}\left(n_{i\uparrow}-1/2\right)\left(n_{i\downarrow}-1/2\right),$
(1)

where $t$ is the hopping integral between the nearest neighbor sites $\langle
ij\rangle$, $\mu$ is chemical potential, and $U$ is the on-site Coulomb
interaction strength. $c_{i\sigma}^{\dagger}$ ($c_{i\sigma}$) is the creation
(annihilation) operator for an electron on site $i$ with spin
$\sigma=\uparrow,\downarrow$ and
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ measures the number of
electrons of spin $\sigma$ on site $i$. As this model only has nearest-
neighbor hopping, it preserves particle-hole symmetry at half-filling with
$\mu=0$. A uniform, orbital magnetic field is introduced by the Peierls
substitution via the phase

$\varphi_{ij}=\dfrac{2\mathrm{\pi}}{\Phi_{0}}\int_{\mathbf{R}_{i}}^{\mathbf{R}_{j}}\mathbf{A}\cdot
d\bm{\ell},$ (2)

where the integral is taken over the shortest straight line path,
$\Phi_{0}=h/e$ is the magnetic flux quantum, and $\mathbf{R}_{i}$ is the
position of site $i$. We choose the symmetric gauge
$\mathbf{A}=(-y\mathbf{\hat{x}}+x\mathbf{\hat{y}})B/2$ and do not include any
Zeeman coupling terms.

We simulate the Hamiltonian in Eq. 1 on a finite cluster with lattice constant
$a=1$, and $N_{x}$ and $N_{y}$ sites in the $x$ and $y$ directions,
respectively. $N=N_{x}N_{y}$ denotes the total number of sites. We implement
modified periodic boundary conditions consistent with magnetic translation
symmetry [44]. Requiring that the wave function be single-valued on the torus
gives the flux quantization condition $\Phi/\Phi_{0}=n_{f}/N$, where
$\Phi=Ba^{2}$ is the flux through a plaquette and $n_{f}$ is an integer.

Allowing the hopping integral to carry a complex phase requires us to modify
the standard DQMC algorithm to use complex numbers, which increases the run-
time of our algorithm $\sim 3$-fold. The complexified DQMC algorithm retains
the same $O(M^{3}L)$ scaling as the real DQMC algorithm, where $M=N_{x}=N_{y}$
is the linear size of the lattice, and $L$ is the number of imaginary time
discretization steps. Unless otherwise specified, all Hubbard-Hofstadter DQMC
simulations are performed on a $N_{x}=N_{y}=8$ square cluster. Error bars in
DQMC results, when shown, denote $\pm 1$ standard error of the mean, estimated
by jackknife resampling. Detailed simulation parameters are listed in
Supplementary Note 1. Non-interacting results, where shown, are obtained from
diagonalizing the Hofstadter model on a $40\times 40$ cluster in order to
minimize finite size effects.

##  III Results and Discussion

###  III.1 Interacting Gap Structure

In Fig. 1, we show the electron density $\langle n\rangle$ vs. chemical
potential $\mu$ at different field strengths for $U/t=0-8$. In Fig. 1(a) we
observe an electron density plateau where there are energy gaps between
magnetic Bloch bands, $\langle n\rangle=2B\nu/\Phi_{0},\nu\in\mathbb{Z}$. As
$U$ increases in Fig. 1(b)-(e), these plateaus are weakened and pushed outward
in chemical potential, but inflections of the $\langle n\rangle$ vs. $\mu$
curves are still visible at the same density values, indicating that
degeneracy of Landau levels is not modified by $U$. Since we are at relatively
high temperature, only the most prominent band gaps (with Chern number $C=\pm
1$) $\langle n\rangle=2B/\Phi_{0}$ remain visible at larger $U$ values.
Additionally, at half filling, as $U$ increases, a Mott gap appears and widens
for all values of magnetic field. In Fig. 1(d)-(e), for $U/t\geq 6$, when the
Mott gap is well-defined, it decreases monotonically as the magnetic field
increases, consistent with previous exact diagonalization results [24]. The
same trend can be seen in a “correlated Hofstadter butterfly” plot, as shown
in Supplementary Fig. S1 and described in Supplementary Note 2. We will
discuss later, and in more detail, the behavior of the Mott gap.

Figure 1: Electron density $\langle n\rangle$ vs. chemical potential $\mu$ (a)
in the non-interacting system, and (b)-(e) with Hubbard $U/t=2-8$. Curves with
the same color have the same magnetic field strength $\Phi/\Phi_{0}$ across
all panels.Each curve is plotted with an offset $2B/\Phi_{0}$ in order to
improve visibility of inflection points. Error bars, corresponding to $\pm 1$
standard error of the mean, estimated by jackknife resampling, are smaller
than the size of data points.All plots have inverse temperature $\beta=4/t$.

It is instructive to plot our data as Wannier diagrams [45], i.e. color
intensity plots of charge compressibility $\chi=\partial\langle
n\rangle/\partial\mu$ as a function of electron density $\langle n\rangle$ and
magnetic field strength $B$. Charge compressibility, or thermodynamic density
of states, is directly measurable in experiments [9, 13]. In a non-interacting
system, at zero temperature, charge compressibility is equivalent to the
single-particle density of states. We measure charge compressibility in DQMC
simulations as

$\chi=\frac{\beta}{N}\sum_{ij}\left[\langle n_{i}n_{j}\rangle-\langle
n_{i}\rangle\langle n_{j}\rangle\right],$ (3)

where $n_{i}=n_{i\uparrow}+n_{i\downarrow}$. In Fig. 2, we show Wannier
diagrams for $U/t=0-6$. For all values of $U$, we observe local minima of
$\chi$ (indicating incompresssible states) along straight lines satisfying the
Diophantine equation

$\frac{\langle n\rangle}{n_{0}}=r\left(\frac{\Phi}{\Phi_{0}}\right)+s,$ (4)

where $r$ and $s$ are integers, and $n_{0}=2$ is the electron density of the
completely filled system. This is consistent with what we expect from the
Hofstadter spectrum in the non-interacting system [45]. The most prominent
incompressible state with $r=1,s=0$ remains clearly visible up to $U/t=8$.
Less prominent incompressible states with $r=2$ and $r=3$ persist to $U/t=4$
and $U/t=2$, respectively. These results show that the integer quantum Hall
states for $r\leq 3,s=0$ have no weak coupling instabilities with respect to
Hubbard repulsion; the $r=1,s=0$ state remains stable up to large $U$. At half
filling, the vertical compressibility minima indicative of the Mott gap
becomes visible for $U/t\gtrsim 4$. Thus, we argue that the gross structure of
magnetic Bloch bands and band gaps determined by the non-interacting
Hofstadter spectrum is preserved in the presence of $U$, with the Mott gap at
half-filling when $U/t\gtrsim 4$ superimposed as an additional feature. Due to
the sign problem, our DQMC simulations are restricted to relatively high
temperature $\beta t\leq 5$, so we cannot resolve conclusively how much $U$
changes the fine structure of the Hofstadter spectrum.

Figure 2: Wannier diagrams (a) in the non-interacting system and (b)-(d) with
Hubbard $U/t=2-6$. Grey regions in (b)-(d) are parameter regions where we
don’t have simulation data. All plots have inverse temperature $\beta=5/t$.
Magnetic field strength is displayed as $\Phi/\Phi_{0}$, and electron density
is shown as $\langle n\rangle/n_{0}$, where $n_{0}$ is the electron density of
a completely filled system, which in our case is $2$. The system is particle-
hole symmetric, so only the range $\langle n\rangle/n_{0}\in[0,0.5]$ is shown.
Dotted cyan lines indicate where we expect local minima of charge
compressibility to occur from the non-interacting model.

###  III.2 Fermion Sign

An important quantity in QMC simulations of interacting fermions is the
fermion sign. The Hubbard-Hofstadter model is sign-problem-free at half
filling on a bipartite lattice. But the fermion sign problem [46]
fundamentally prevents us from obtaining high quality simulation data at low
temperatures and away from half filling. Thus, any insight into factors
affecting the severity of the sign problem is valuable. Since the fermion sign
problem is NP-hard [47], we do not expect a general solution to the fermion
sign problem to exist. Nevertheless, as the sign problem is representation-
dependent, it is possible to reduce or completely remove the sign problem for
specific classes of non-generic Hamiltonians [48].

Figure 3: (a) Charge compressibility $\chi=\partial\langle
n\rangle/\partial\mu$ and (b) fermion sign $\langle s\rangle$ plotted against
electron density $\langle n\rangle$ at magnetic field strength
$\Phi/\Phi_{0}=12/64$, inverse temperature $\beta=6/t$ and Hubbard interaction
$U/t=4$. This system is particle-hole symmetric, so we only plot the density
range $\langle n\rangle\in[0,1]$. Dashed lines indicate electron densities at
which $\chi$ reaches local minina and $\langle s\rangle$ reaches local maxima.
Error bars, corresponding to $\pm 1$ standard error of the mean, estimated by
jackknife resampling, are smaller than the size of data points.

In this work, we find a correlation between the fermion sign and the charge
compressibility. In Fig. 3, we show the fermion sign $\langle s\rangle$ and
charge compressibility $\chi$, both plotted against $\langle n\rangle$, for
one representative set of parameters. Local minima of charge compressibility
in this interacting system exactly correspond to local maxima of the fermion
sign. At these local maxima, the fermion sign may be an order of magnitude
improved over its value at other electron densities and that of the standard
zero-field Hubbard model. For an extended figure demonstrating that this
correspondence is general across our parameter space and not a finite size
artifact, see Supplementary Fig. S2 and Supplementary Note 3. Our results may
mean that although the Hubbard model in general suffers from a sign problem,
it is possible to obtain good results when we are precisely located on an
integer quantum Hall plateau. Since similar sign-compressibility
correspondence has been reported [30, 49, 50, 51], it appears that the
improvement of fermion sign in insulating phases is quite general, consistent
with our intuition that fermionic statistics become less important in
localized states. Our results also relate to recent work [52, 53, 54]
suggesting that the fermion sign is not merely a coincidental barrier to
accessing low-temperature physics, but may be reflective of intrinsic physics
of model Hamiltonians.

###  III.3 Half-Filling

Finally, we focus on half filling, where we believe interesting interplay
between Hofstadter physics and Hubbard physics occurs. In the absence of a
magnetic field, the ground state of the half-filled Hubbard model is an AFMI
at any nonzero value of $U$, i.e. ($U_{\mathrm{c}}=0$) [29, 30], due to
perfect nesting of the Fermi surface and a logarithmically divergent single-
particle density of states. In the limit of strong interactions $U\gg t$, the
half-filled Hubbard model maps to the Heisenberg model with antiferromagnetic
nearest neighbor spin exchange energy $J=4t^{2}/U$ [55].

In the presence of an orbital magnetic field, the non-interacting density of
states at half filling is modified significantly. As can be seen in Fig. 2(a),
the density of states/charge compressibility does not change monotonically
with field, but instead shows prominent minima at $\Phi/\Phi_{0}=p/q$, where
$p$ and $q$ are co-prime and $q$ is even, corresponding to a non-interacting
ground state with $q$ inequivalent Dirac cones [37]. The large-field limit
corresponds to the $\mathrm{\pi}$-flux model, in which a semimetal$-$AFMI
transition occurs at $U_{\mathrm{c}}\approx 5.6t$. As the orbital magnetic
field significantly changes the non-interacting density of states at half
filling, we expect that critical $U_{\mathrm{c}}$ should exhibit $B$ field
dependence. It would be interesting to investigate if $U_{\mathrm{c}}$ changes
monotonically with $B$, or if it exhibits non-monotonicity commensurate with
the oscillatory behavior of the density of states. We defer the mapping of
this $U$-$B$ phase diagram to future work.

For the remainder of this section, we focus on the parameter region
$U/t\in[6,10]$. Here, the system is safely an AFMI at all field strengths. We
examine the evolution of local magnetic moment $\langle m_{z}^{2}\rangle$,
antiferromagnetic structure factor $S(\mathrm{\pi},\mathrm{\pi})$ and specific
heat $c_{v}=\partial\langle E\rangle/\partial T$ with magnetic field strength,
and see that these thermodynamic quantities all consistently show that a
strong orbital magnetic field tends to modify the AFMI by delocalizing
electrons and thereby reducing the effect of $U$ on the low-energy properties
of the Mott insulating phase. For finite-size analysis of thermodynamic
observables at half filling, see Supplementary Fig. S3 and Supplementary Note
3.

Figure 4: Temperature and field dependence of local moment $\langle
m_{z}^{2}\rangle$ at half filling. Insets display the field dependence of
magnetic structure factor $S(\mathrm{\pi},\mathrm{\pi})$ at $\beta t=16$ for
magnetic field strength $\Phi/\Phi_{0}\in[0,0.5]$. (a)-(c) Correspond to
Hubbard interaction strength $U/t=6-10$, respectively. Curves with the same
color and marker type have the same magnetic field strength across all panels.
Error bars, corresponding to $\pm 1$ standard error of the mean, estimated by
jackknife resampling, are smaller than the size of data points. Figure 5:
Temperature and field dependence of specific heat $c_{v}$ at half filling.
Insets display the field dependence of kinetic energy $\langle K\rangle$ at
$\beta t=16$ for magnetic field strength $\Phi/\Phi_{0}\in[0,0.5]$. (a)-(c)
Correspond to Hubbard interaction strength $U/t=6-10$, respectively. Values of
$J=4t^{2}/U$ are also indicated for each $U$. Curves with the same color and
marker type have the same magnetic field strength across all panels. Error
bars denote $\pm 1$ standard error of the mean, estimated by jackknife
resampling.

In Fig. 4, we show the temperature and field dependence of the local moment.
The local moment or sublattice magnetization

$\langle
m_{z}^{2}\rangle=\frac{1}{N}\sum_{i}\langle(n_{i\uparrow}-n_{i\downarrow})^{2}\rangle$
(5)

measures the degree of spin localization. It is $0.5$ in the non-interacting
system and approaches $1$ in the $U/t\rightarrow\infty$ limit. In the zero-
field Hubbard model, $\langle m_{z}^{2}\rangle$ has features at $T\sim U$
associated with the formation of local moments, and at $T\sim J$, associated
with short- or long-range ordering of local moments [56, 36]. We see that at
fixed $U$, increasing magnetic field strength reduces the local moment
monotonically at all temperatures, with the effect largest below temperatures
$T\sim J$. The zero-field and $\mathrm{\pi}$-flux limit of local moment data
are consistent with previous work [43].

The magnetic structure factor is the Fourier transform of the real-space spin-
spin correlation function

$S(\mathbf{Q})=\frac{1}{N}\sum_{ij}e^{\mathrm{i}\mathbf{Q}\cdot(\mathbf{R}_{i}-\mathbf{R}_{j})}\langle(n_{i\uparrow}-n_{i\downarrow})(n_{j\uparrow}-n_{j\downarrow})\rangle.$
(6)

When the system has long-range antiferromagnetic order, the structure factor
is strongly peaked at the ordering wave vector
$\mathbf{Q}=(\mathrm{\pi},\mathrm{\pi})$, with peak height scaling linearly
with lattice size [36, 57]. In our simulations, we find that at all $B$ and
$U$ values, the magnetic structure factor is sharply peaked at
$(\mathrm{\pi},\mathrm{\pi})$, consistent with the system being in the AFMI
phase. Insets to Fig. 4 show that at all $U$, the magnetic field monotonically
reduces $S(\mathrm{\pi},\mathrm{\pi})$, indicating that the magnetic field
reduces the AFMI ordering tendencies.

Fig. 5 shows the evolution of the low temperature peak of $c_{v}$ with
magnetic field and Hubbard $U$. We calculate $c_{v}$ numerically by measuring
energy as a function of temperature $\langle E(T)\rangle$ and taking the
finite difference $\Delta\langle E\rangle/\Delta T$. In the zero-field half-
filled Hubbard model, at large $U$, the specific heat has a “two peak”
structure, with a broad high temperature peak at $T\sim U$ associated with
charge fluctuations and and a narrow low temperature peak at $T\sim J$
associated with spin fluctuations [56, 58]. When a magnetic field is turned
on, as shown in Fig. 5(b)-(c), $U$ is large enough that the two peaks remain
well-separated. The high temperature peak doesn’t move, while the low
temperature peak shifts to higher temperatures. In Fig. 5(a), the two peaks
are initially well-separated at low fields, but at $\Phi/\Phi_{0}\gtrsim 1/4$,
the low temperature peak shifts upwards and merges with the high-temperature
peak, complicating our interpretation. This is likely due to $U/t=6$ being low
enough for the system to not simply map to the Heisenberg model, and for the
system to be close to the AFMI phase transition at $\mathrm{\pi}$-flux.

Since the low temperature peak in specific heat is associated with the spin
exchange energy $J$, we are tempted to say that the orbital magnetic field
increases $J$. However, this interpretation may be overly naive. We believe
the more accurate statement is that the orbital magnetic field tends to
delocalize electrons, and thus, effectively lower the influence of $U$ on low
energy properties of the system. Insets to Fig. 5 show that a magnetic field
increases (in magnitude) kinetic energy in the insulating phase, which
supports this interpretation. This runs contrary to our usual intuition that
an orbital magnetic field localizes electrons by winding them up into Landau
orbits. However, here our starting point is a correlated insulator, rather
than free electrons (or a Fermi liquid). Our results in Figs. 4 and 5, along
with the decreasing width of the Mott gap in Fig. 1(d)-(e), suggest that in
the AFMI phase, the orbital magnetic field tends to delocalize electrons,
increase kinetic energy, and lower the effective influence of $U$. We observe
that the influence of magnetic field is suppressed as $U$ increases. As
$U/t\rightarrow\infty$, the influence from the $B$ field will diminish and
become negligible, since in the atomic limit, no hopping exists, and the
orbital magnetic field cannot have an influence on the system.

##  IV Conclusions

In this work, we implemented DQMC to simulate the Hubbard-Hofstadter model and
directly investigate field dependent thermodynamic properties of correlated
electrons, specifically focusing on the charge compressibility, local moment,
magnetic structure factor, and specific heat. By examining charge
compressibility, we find that magnetic Bloch bands are smeared out by the
Hubbard interaction, but the non-interacting band gaps away from half-filling
persist at temperatures accessible to DQMC in the presence of local
correlations. At half filling, we find that the orbital magnetic field and
Hubbard potential act antagonistically. At intermediate to strong coupling
$U/t\in[6,10]$, strong orbital magnetic fields reduce the apparent width of
the Mott gap, reduce the magnitude of the local moment and magnetic structure
factor, increase kinetic energy, and shift the low-$T$ peak of specific heat
to higher temperatures. Together, these phenomena indicate that an orbital
magnetic field tends to delocalize electrons and reduce the effect of $U$.
From the algorithmic perspective, we find that the fermion sign in DQMC
simulations is improved significantly when the physical system is
incompressible.

For future work, we are interested in mapping out the $U$-$B$ phase diagram
for $U/t\in[0,6]$ and in studying the evolution of $U_{\mathrm{c}}$ with
magnetic field strength. At smaller values of $U$, we can achieve lower
temperatures in DQMC simulations, but finite size effects also become more
significant. Pinning down the precise location of the phase transition will
require simulations on much larger lattices, as well as careful finite size
scaling analysis.

As this work is partially motivated by understanding the large variety of
exotic states in Moiré materials [9, 10, 12, 11, 13, 14], another potential
direction for future work is to use methods developed here to study the
Hubbard-Hofstadter model on the honeycomb or triangular lattice. It has been
suggested that the single-band Hubbard model on a triangular lattice is
directly applicable to transition metal dichalcogenide heterobilayers [59,
60]. However, due to the linear dispersion relation of Dirac fermions, the
Coulomb interaction is poorly screened in graphene. It is likely that some
multi-orbital, extended Hubbard model may be required to capture the full
effect of electron correlations, for example, in twisted bilayer graphene [61,
62].

##  V Data Availability

Aggregated numerical data and analysis routines required to reproduce the
figures can be found at https://doi.org/10.5281/zenodo.6383764. Raw simulation
data that support the findings of this study are stored on the Sherlock
cluster at Stanford University and are available from the corresponding author
upon reasonable request.

##  VI Code Availability

The most up-to-date version of our DQMC simulation code can be accessed at
https://github.com/edwnh/dqmc.

##  VII Author Contributions

EWH wrote the DQMC simulation code. JKD performed simulations and data
analysis. TPD conceptualized the work. All authors (JKD, WOW, BM, YS, EWH,
TPD) participated in discussions and manuscript writing.

##  VIII Competing Interests

The authors declare no competing interests.

##  IX Acknowledgements

We acknowledge helpful discussions with Richard Scalettar, Philip W. Phillips,
Allan.H. MacDonald, Benjamin E. Feldman, Young S. Lee and Jiachen Yu.

This work was supported by the U.S. Department of Energy (DOE), Office of
Basic Energy Sciences, Division of Materials Sciences and Engineering. EWH was
supported by the Gordon and Betty Moore Foundation’s EPiQS Initiative through
grants GBMF 4305 and GBMF 8691. YS was supported by the Gordon and Betty Moore
Foundation’s EPiQS Initiative through grants GBMF 4302 and GBMF 8686.
Computational work was performed on the Sherlock cluster at Stanford
University and on resources of the National Energy Research Scientific
Computing Center, supported by the U.S. DOE, Office of Science, under Contract
no. DE-AC02-05CH11231.

## References

  * Wu  _et al._ [2011] T. Wu, H. Mayaffre, S. Krämer, M. Horvatić, C. Berthier, W. N. Hardy, R. Liang, D. A. Bonn, and M.-H. Julien, Magnetic-field-induced charge-stripe order in the high-temperature superconductor $\mathrm{YBa_{2}Cu_{3}O_{y}}$, Nature 477, 191 (2011).
  * Gerber  _et al._ [2015] S. Gerber, H. Jang, H. Nojiri, S. Matsuzawa, H. Yasumura, D. A. Bonn, R. Liang, W. N. Hardy, Z. Islam, A. Mehta, S. Song, M. Sikorski, D. Stefanescu, Y. Feng, S. A. Kivelson, T. P. Devereaux, Z.-X. Shen, C.-C. Kao, W.-S. Lee, D. Zhu, and J.-S. Lee, Three-dimensional charge density wave order in $\mathrm{YBa_{2}Cu_{3}O_{6.67}}$ at high magnetic fields, Science 350, 949 (2015).
  * Jang  _et al._ [2016] H. Jang, W.-S. Lee, H. Nojiri, S. Matsuzawa, H. Yasumura, L. Nie, A. V. Maharaj, S. Gerber, Y.-J. Liu, A. Mehta, D. A. Bonn, R. Liang, W. N. Hardy, C. A. Burns, Z. Islam, S. Song, J. Hastings, T. P. Devereaux, Z.-X. Shen, S. A. Kivelson, C.-C. Kao, D. Zhu, and J.-S. Lee, Ideal charge-density-wave order in the high-field state of superconducting $\mathrm{YBCO}$, Proceedings of the National Academy of Sciences 113, 14645 (2016).
  * Edkins  _et al._ [2019] S. D. Edkins, A. Kostin, K. Fujita, A. P. Mackenzie, H. Eisaki, S. Uchida, S. Sachdev, M. J. Lawler, E.-A. Kim, J. C. S. Davis, and M. H. Hamidian, Magnetic field-induced pair density wave state in the cuprate vortex halo, Science 364, 976 (2019).
  * Grigera  _et al._ [2001] S. A. Grigera, R. S. Perry, A. J. Schofield, M. Chiao, S. R. Julian, G. G. Lonzarich, S. I. Ikeda, Y. Maeno, A. J. Millis, and A. P. Mackenzie, Magnetic field-tuned quantum criticality in the metallic ruthenate $\mathrm{Sr_{3}Ru_{2}O_{7}}$, Science 294, 329 (2001).
  * Custers  _et al._ [2003] J. Custers, P. Gegenwart, H. Wilhelm, K. Neumaier, Y. Tokiwa, O. Trovarelli, C. Geibel, F. Steglich, C. Pépin, and P. Coleman, The break-up of heavy electrons at a quantum critical point, Nature 424, 524 (2003).
  * Lévy  _et al._ [2005] F. Lévy, I. Sheikin, B. Grenier, and A. D. Huxley, Magnetic field-induced superconductivity in the ferromagnet $\mathrm{URhGe}$, Science 309, 1343 (2005).
  * Ran  _et al._ [2019] S. Ran, I.-L. Liu, Y. S. Eo, D. J. Campbell, P. M. Neves, W. T. Fuhrman, S. R. Saha, C. Eckberg, H. Kim, D. Graf, F. Balakirev, J. Singleton, J. Paglione, and N. P. Butch, Extreme magnetic field-boosted superconductivity, Nature Physics 15, 1250 (2019).
  * Hunt  _et al._ [2013] B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, and R. C. Ashoori, Massive Dirac fermions and Hofstadter butterfly in a van der Waals heterostructure, Science 340, 1427 (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 (2015).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Spanton  _et al._ [2018] E. M. Spanton, A. A. Zibrov, H. Zhou, T. Taniguchi, K. Watanabe, M. P. Zaletel, and A. F. Young, Observation of fractional Chern insulators in a Van der Waals heterostructure, Science 360, 62 (2018).
  * Yu  _et al._ [2022] J. Yu, B. A. Foutty, Z. Han, M. E. Barber, Y. Schattner, K. Watanabe, T. Taniguchi, P. Phillips, Z.-X. Shen, S. A. Kivelson, and B. E. Feldman, Correlated hofstadter spectrum and flavour phase diagram in magic-angle twisted bilayer graphene, Nature Physics 18, 825 (2022).
  * Saito  _et al._ [2021] Y. Saito, J. Ge, L. Rademaker, K. Watanabe, T. Taniguchi, D. A. Abanin, and A. F. Young, Hofstadter subband ferromagnetism and symmetry-broken Chern insulators in twisted bilayer graphene, Nature Physics 17, 478 (2021).
  * Hofstadter [1976] D. R. Hofstadter, Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * MacDonald [1983] A. H. MacDonald, Landau-level subband structure of electrons on a square lattice, Phys. Rev. B 28, 6713 (1983).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * Miyake  _et al._ [2013] H. Miyake, G. A. Siviloglou, C. J. Kennedy, W. C. Burton, and W. Ketterle, Realizing the Harper hamiltonian with laser-assisted tunneling in optical lattices, Phys. Rev. Lett. 111, 185302 (2013).
  * 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, 598 (2013).
  * Ponomarenko  _et al._ [2013] L. A. Ponomarenko, R. V. Gorbachev, G. L. Yu, D. C. Elias, R. Jalil, A. A. Patel, A. Mishchenko, A. S. Mayorov, C. R. Woods, J. R. Wallbank, M. Mucha-Kruczynski, B. A. Piot, M. Potemski, I. V. Grigorieva, K. S. Novoselov, F. Guinea, V. I. Fal’ko, and A. K. Geim, Cloning of dirac fermions in graphene superlattices, Nature 497, 594 (2013).
  * Gudmundsson and Gerhardts [1995] V. Gudmundsson and R. R. Gerhardts, Effects of screening on the Hofstadter butterfly, Phys. Rev. B 52, 16744 (1995).
  * Doh and Salk [1998] H. Doh and S.-H. S. Salk, Effects of electron correlations on the Hofstadter spectrum, Phys. Rev. B 57, 1312 (1998).
  * Barelli  _et al._ [1996] A. Barelli, J. Bellissard, P. Jacquod, and D. L. Shepelyansky, Double butterfly spectrum for two interacting particles in the harper model, Phys. Rev. Lett. 77, 4752 (1996).
  * Czajka  _et al._ [2006] K. Czajka, A. Gorczyca, M. M. Maśka, and M. Mierzejewski, Hofstadter butterfly for a finite correlated system, Phys. Rev. B 74, 125116 (2006).
  * Acheche  _et al._ [2017] S. Acheche, L.-F. Arsenault, and A.-M. S. Tremblay, Orbital effect of the magnetic field in dynamical mean-field theory, Phys. Rev. B 96, 235135 (2017).
  * Markov  _et al._ [2019] A. A. Markov, G. Rohringer, and A. N. Rubtsov, Robustness of the topological quantization of the hall conductivity for correlated lattice electrons at finite temperatures, Phys. Rev. B 100, 115102 (2019).
  * Tu  _et al._ [2018] W.-L. Tu, F. Schindler, T. Neupert, and D. Poilblanc, Competing orders in the Hofstadter ${t}-{J}$ model, Phys. Rev. B 97, 035154 (2018).
  * Blankenbecler  _et al._ [1981] R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, Monte carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24, 2278 (1981).
  * Hirsch [1985] J. E. Hirsch, Two-dimensional hubbard model: Numerical simulation study, Phys. Rev. B 31, 4403 (1985).
  * White  _et al._ [1989] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, Numerical study of the two-dimensional hubbard model, Phys. Rev. B 40, 506 (1989).
  * Jia  _et al._ [2014] C. J. Jia, E. A. Nowadnick, K. Wohlfeld, Y. F. Kung, C.-C. Chen, S. Johnston, T. Tohyama, B. Moritz, and T. P. Devereaux, Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering, Nature Communications 5, 3314 (2014).
  * Kung  _et al._ [2015] Y. F. Kung, E. A. Nowadnick, C. J. Jia, S. Johnston, B. Moritz, R. T. Scalettar, and T. P. Devereaux, Doping evolution of spin and charge excitations in the hubbard model, Phys. Rev. B 92, 195108 (2015).
  * Khatami  _et al._ [2015] E. Khatami, R. T. Scalettar, and R. R. P. Singh, Finite-temperature superconducting correlations of the hubbard model, Phys. Rev. B 91, 241107 (2015).
  * Huang  _et al._ [2017] E. W. Huang, C. B. Mendl, S. Liu, S. Johnston, H.-C. Jiang, B. Moritz, and T. P. Devereaux, Numerical evidence of fluctuating stripes in the normal state of high-tc cuprate superconductors, Science 358, 1161 (2017).
  * Huang  _et al._ [2019a] E. W. Huang, R. Sheppard, B. Moritz, and T. P. Devereaux, Strange metallicity in the doped Hubbard model, Science 366, 987 (2019a).
  * Varney  _et al._ [2009] C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell, and R. T. Scalettar, Quantum monte carlo study of the two-dimensional fermion hubbard model, Phys. Rev. B 80, 075116 (2009).
  * Wen and Zee [1989] X. Wen and A. Zee, Winding number, family index theorem, and electron hopping in a magnetic field, Nuclear Physics B 316, 641 (1989).
  * Affleck and Marston [1988] I. Affleck and J. B. Marston, Large-n limit of the Heisenberg-Hubbard model: Implications for high-${T}_{c}$ superconductors, Phys. Rev. B 37, 3774 (1988).
  * Chang and Scalettar [2012] C.-C. Chang and R. T. Scalettar, Quantum disordered phase near the mott transition in the staggered-flux Hubbard model on a square lattice, Phys. Rev. Lett. 109, 026404 (2012).
  * Otsuka  _et al._ [2014] Y. Otsuka, S. Yunoki, and S. Sorella, Mott transition in the 2d hubbard model with $\pi$-flux JPS Conference Proceedings, 3, 10.7566/JPSCP.3.013021 (2014), 0.
  * Parisen Toldin  _et al._ [2015] F. Parisen Toldin, M. Hohenadler, F. F. Assaad, and I. F. Herbut, Fermionic quantum criticality in honeycomb and $\pi$-flux Hubbard models: Finite-size scaling of renormalization-group-invariant observables from quantum monte carlo, Phys. Rev. B 91, 165108 (2015).
  * Otsuka  _et al._ [2016] Y. Otsuka, S. Yunoki, and S. Sorella, Universal quantum criticality in the metal-insulator transition of two-dimensional interacting dirac electrons, Phys. Rev. X 6, 011029 (2016).
  * Guo  _et al._ [2018] H. Guo, E. Khatami, Y. Wang, T. P. Devereaux, R. R. P. Singh, and R. T. Scalettar, Unconventional pairing symmetry of interacting dirac fermions on a $\pi$-flux lattice, Phys. Rev. B 97, 155146 (2018).
  * Assaad [2002] F. F. Assaad, Depleted kondo lattices: Quantum monte carlo and mean-field calculations, Phys. Rev. B 65, 115104 (2002).
  * Wannier [1978] G. H. Wannier, A Result Not Dependent on Rationality for Bloch Electrons in a Magnetic Field, Physica Status Solidi B Basic Research 88, 757 (1978).
  * Loh  _et al._ [1990] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, Sign problem in the numerical simulation of many-electron systems, Phys. Rev. B 41, 9301 (1990).
  * Troyer and Wiese [2005] M. Troyer and U.-J. Wiese, Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * Li and Yao [2019] Z.-X. Li and H. Yao, Sign-problem-free fermionic quantum monte carlo: Developments and applications, Annual Review of Condensed Matter Physics 10, 337 (2019).
  * Mondaini  _et al._ [2012] R. Mondaini, K. Bouadim, T. Paiva, and R. R. dos Santos, Finite-size effects in transport data from quantum monte carlo simulations, Phys. Rev. B 85, 125127 (2012).
  * Kung  _et al._ [2016] Y. F. Kung, C.-C. Chen, Y. Wang, E. W. Huang, E. A. Nowadnick, B. Moritz, R. T. Scalettar, S. Johnston, and T. P. Devereaux, Characterizing the three-orbital hubbard model with determinant quantum monte carlo, Phys. Rev. B 93, 155166 (2016).
  * Huang  _et al._ [2019b] E. W. Huang, M.-S. Vaezi, Z. Nussinov, and A. Vaezi, Enhanced correlations and superconductivity in weakly interacting partially flat-band systems: A determinantal quantum monte carlo study, Phys. Rev. B 99, 235128 (2019b).
  * Mondaini  _et al._ [2022] R. Mondaini, S. Tarat, and R. T. Scalettar, Quantum critical points and the sign problem, Science 375, 418 (2022).
  * Wessel  _et al._ [2017] S. Wessel, B. Normand, F. Mila, and A. Honecker, Efficient Quantum Monte Carlo simulations of highly frustrated magnets: the frustrated spin-1/2 ladder, SciPost Phys. 3, 005 (2017).
  * Götz  _et al._ [2022] A. Götz, S. Beyl, M. Hohenadler, and F. F. Assaad, Valence-bond solid to antiferromagnet transition in the two-dimensional su-schrieffer-heeger model by langevin dynamics, Phys. Rev. B 105, 085151 (2022).
  * Fazekas [1999] P. Fazekas,  _Lecture Notes on Electron Correlation and Magnetism_ (World Scientific, 1999).
  * Paiva  _et al._ [2001] T. Paiva, R. T. Scalettar, C. Huscroft, and A. K. McMahan, Signatures of spin and charge energy scales in the local moment and specific heat of the half-filled two-dimensional hubbard model, Phys. Rev. B 63, 125116 (2001).
  * Huse [1988] D. A. Huse, Ground-state staggered magnetization of two-dimensional quantum heisenberg antiferromagnets, Phys. Rev. B 37, 2380 (1988).
  * Wang  _et al._ [2022] W. O. Wang, J. K. Ding, B. Moritz, E. W. Huang, and T. P. Devereaux, Magnon heat transport in a two-dimensional mott insulator, Phys. Rev. B 105, L161103 (2022).
  * Wu  _et al._ [2018] F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Hubbard model physics in transition metal dichalcogenide moiré bands, Phys. Rev. Lett. 121, 026402 (2018).
  * Tang  _et al._ [2020] Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Simulation of hubbard model physics in wse2/ws2 moiré superlattices, Nature 579, 353 (2020).
  * Andrews and Soluyanov [2020] B. Andrews and A. Soluyanov, Fractional quantum hall states for moiré superstructures in the Hofstadter regime, Phys. Rev. B 101, 235312 (2020).
  * Po  _et al._ [2019] H. C. Po, L. Zou, T. Senthil, and A. Vishwanath, Faithful tight-binding models and fragile topology of magic-angle bilayer graphene, Phys. Rev. B 99, 195455 (2019).


# Double layered transformation cavities: Crossover between local and global
transformations

Jung-Wan Ryu<EMAIL_ADDRESS>Center for Theoretical Physics of Complex
Systems, Institute for Basic Science (IBS), Daejeon 34126, Republic of Korea

###### Abstract

We study the resonant modes in double-layered transformation cavities
consisting of inner and outer layer boundaries where the mode intensity is
mainly located inside the inner layer. This is an intermediate design between
a transformation cavity and a cavity in a wholly transformed space of
transformation optics. We demonstrate the crossover between these two extreme
cases as the outer layer of the cavity varies and also explore the properties
of the resonant modes in the cavity. While the near-field patterns of the
resonant modes do not change as the outer layer becomes larger, the Q-factors
approach those of a cavity in the wholly transformed space of transformation
optics, and the far-field patterns are modified.

##  I Introduction

Transformation optics (TO) provides a general design method to manipulate
electromagnetic waves in a desired manner Leo06 ; Pen06 . It is based on the
concept that Maxwell’s equations can be form-invariant under coordinate
transformations, such that only the permittivity and permeability tensors are
modified Pos62 ; Jac75 . Consequently, TO has become a powerful tool to
control the propagation of electromagnetic waves Che10-2 ; Kil11 ; Liu12 .
Many novel optical devices and unprecedented applications have been
demonstrated with TO, including individual cloaks Pen06 ; Sch06 ; Cai07 ;
Lan13 , wave concentrators Jia08 ; Sad15 , optical black holes Gen09 ; Che10 ;
Deh15 , waveguide devices Rob08 ; Zha09 ; Tic10 , and negative refractive
index materials Ves68 ; She01 ; Zha16 .

Transformation optics can be applied not only to the control of the
propagation path of electromagnetic waves in a desired manner but also to the
manipulation of the resonant mode properties of two-dimensional dielectric
optical cavities Kim16 . One new kind of microcavity designed utilizing
conformal TO to manipulate the resonant mode properties is called a
transformation cavity (TC), which has a deformed boundary shape and
inhomogeneous profile of the refractive index inside the microcavity
associated with a given conformal mapping. The variation in the Q-factors and
emission directionality of resonant modes in the system parameter space was
systematically studied to find the optimized resonant modes in a limaçon-
shaped cavity Ryu19 . The applicability of TCs can be further extended to
arbitrary shapes by invoking a quasi-conformal mapping method Par19 . It has
been expected that metamaterials, e.g., subwavelength structures such as holes
and posts that can produce homogenized effective media, is one of the
promising candidates for the implementation of TCs, like in numerous
applications of TO theory Ves68 ; Pen99 ; Smi00 ; Smi04 ; Eng06 ; Val09 ;
Gab09 ; Vas10 ; Gao12 .

While space is globally transformed by a given conformal mapping in the case
of original TO, space is locally transformed by mapping only inside the cavity
in the case of TCs. Accordingly, the refractive indices inside and outside
cavities have inhomogeneous profiles in the case of original TO but have
inhomogeneous and homogeneous profiles, respectively, in the case of TCs. This
discrepancy results in significant modifications of the resonant modes,
thereby enabling us to independently control optical properties such as the
Q-factor and emission directionality of the conformal whispering gallery mode
(cWGM) in a TC. In this work, we propose a double-layered transformation
cavity (DLTC) consisting of inner and outer layers with gradually varying
refractive index profiles, which is an intermediate design between a TC and a
cavity in a wholly transformed (WT) space of TO. We study the optical
properties, namely the Q-factors and near- and far-field patterns, of the
resonant modes in the DLTC. In the cavity, as the outer layer boundary becomes
larger while the inner layer boundary does not change, the Q-factors of the
resonant modes approach those of a cavity in WT space, and the emission
directionalities change only by modifying the outer layer boundary shape.

##  II Systems: double layered transformation cavity

Deformed cavity and resonant modes - We consider a two-dimensional dielectric
cavity with a given boundary shape, where the optical modes are described by
the resonant modes. The modes are obtained by solving the following scalar
wave equation,

$[\nabla^{2}+n^{2}(x,y)k^{2}]\psi(x,y)=0,$ (1)

where $n(x,y)$ is the refractive index function at the position $(x,y)$. The
resonant modes should satisfy the outgoing-wave boundary condition,
$\psi(x,y)\sim h_{k}(\theta)e^{ikr}/\sqrt{r}$, where
$(x,y)=(r\cos\theta,r\sin\theta)$ and $h_{k}(\theta)$ is the far-field angular
distribution of the emission. In a deformed cavity with a homogeneous
refractive index profile, $n(x,y)$ is $n_{0}$ inside and $1$ outside the
cavity. The real and imaginary parts of the complex wave number $k$ are,
respectively, $\omega/c$ where $\omega$ is the frequency of the resonant mode
and $c$ is the speed of light, and $-1/(2c\tau)$ where $\tau$ is the lifetime
of the mode. The quality factor $Q$ of the mode is defined as
$Q=2\pi\tau/T=-\mathrm{Re}(k)/2\mathrm{Im}(k)$. Here, we focus on transverse
magnetic (TM) polarization, of which the wave function $\psi(x,y)$ corresponds
to $E_{z}$, the $z$ component of the electric field Jac75 . In the case of TM
polarization, both the wave function $\psi(x,y)$ and its normal derivative
$\partial_{\nu}\psi$ are continuous across the cavity boundary.

Figure 1: Boundary shapes (black) and refractive indices (scaled color) for
(a) a circular cavity in an original virtual space, (b) a limaçon-shaped
cavity in a wholly transformed space, (c) a circular cavity in a reciprocal
virtual space, and (d) a transformation cavity in a physical space. The black
and red arrows denote the conformal mapping of Eq. (2) and its inverse,
respectively. The white regions represent a refractive index of $n=1$.

Four spaces and the difference between WT and physical spaces - There are four
spaces to understand the relations between original TO theory and a TC in
depth: an original virtual (OV) space or a reference space [Fig. 1 (a)], a WT
space corresponding to the physical space of conventional TO theory [Fig. 1
(b)], a physical space for a TC [Fig. 1 (d)], and a reciprocal virtual (RV)
space [Fig. 1 (c)]. We first consider a circular cavity in an OV space. Here,
the cavity has a homogeneous refractive index, and the resonant modes are the
solution of Maxwell’s equations in dielectric media. Maxwell’s equations
remain invariant under a given conformal mapping, by which an OV space is
transformed into a WT space that corresponds to the physical space of
conventional TO theory. In a WT space, the cavity boundary is deformed and the
inhomogeneous profile of a refractive index in a whole space can be determined
by the ratio between the local length scales in both spaces. To make a TC in a
physical space, the refractive index outside the cavity is set to $1$, that
is, the conformal mapping is applied not only to the boundary deformations but
also to the inside of the cavity, unlike in WT space where the mapping is
applied to the whole space. These characteristic settings for TCs result in
novel properties of the resonant modes in TCs. Besides, an RV space is given
by an inverse conformal mapping from the physical space, in which refractive
indices have homogeneous and inhomogeneous profiles inside and outside the
cavity, respectively. Although RV space was originally introduced for
numerical calculations such as boundary element methods Kim16 ; BEM  and
Husimi functions Kim18 , which are used in this work, the space is useful to
understand the resonant mode properties since the interior of a cavity in RV
space is the same as that of a cavity in OV space. While the resonant modes
obtained from OV and corresponding WT spaces are intrinsically the same, modes
in WT and corresponding physical spaces are different in principle.

Figure 2: Boundary shapes (black) and refractive indices (scaled color) for
(a) a concentric circular layered cavity in a reciprocal virtual space and (b)
a double-layered transformation cavity in a physical space. The regions of the
inner and outer layers are $\Omega_{i}$ and $\Omega_{o}$, respectively, and
the refractive indices of the inner and outer layers are $n_{i}(x,y)$ and
$n_{o}(x,y)$, respectively.

Double-layered limaçon-shaped transformation cavity - In WT space, conformal
mapping is applied to the whole space [Fig. 1 (b)], but in physical space,
conformal mapping is applied only to the inside of the cavity [Fig. 1 (d)]. In
order to bridge the gap between these two limiting cases, we propose a second
layer outside of the inner layer. The resulting cavity has two sub-regions
with inhomogeneous refractive index profiles: inner and outer layers, denoted
by $\Omega_{i}$ and $\Omega_{o}$, respectively, as shown in Fig. 2 (b). As an
example, we consider a DLTC with boundaries given by a limaçon shape, which is
one of the widely studied shapes in the fields of quantum billiards Rob83 ,
deformed microcavities Wie08 , and transformation cavities Kim16 . The
corresponding conformal mapping from the unit circle to the limaçon shape is
given by

$\zeta=\beta(\eta+\epsilon\eta^{2}),$ (2)

where $\eta=u+iv$ and $\zeta=x+iy$ are complex variables denoting positions in
virtual spaces, i.e., OV and RV spaces, and physical spaces (including WT
spaces), respectively. $\beta$ is a positive size-scaling parameter and
$\epsilon$ is a deformation parameter. The inner and outer layer boundaries
are given by

$\displaystyle x$ $\displaystyle=$ $\displaystyle\beta(r\cos(\phi)+\epsilon
r^{2}\cos(2\phi)),$ (3) $\displaystyle y$ $\displaystyle=$
$\displaystyle\beta(r\sin(\phi)+\epsilon r^{2}\sin(2\phi)),$

where size scale factors $r$’s are $r_{i}$ and $r_{o}$ for the inner and outer
layer boundaries, which are equivalent to the radii of two circles of a
concentric circular layered cavity in an RV space [Fig. 2 (a)], respectively,
and $\beta=\beta_{max}=1/(1+2\epsilon)$. The maximal value $\beta_{max}$ of
$\beta$, that supports the cWGMs in a limaçon-shaped TC with $r=1.0$, can be
obtained from the condition $|d{\zeta}/d{\eta}|^{-1}\geq 1$ necessary for
total internal reflection in the TC with outside refractive index $n_{out}=1$.
The refractive index profiles are given by

$\displaystyle
n(x,y)=\begin{cases}n_{0}\left|\frac{d\zeta}{d\eta}\right|^{-1}=\frac{n_{0}}{\beta\left|\sqrt{1+4\epsilon\zeta/\beta}\right|},~{}&(x,y)\in\Omega_{i}\\\
\left|\frac{d\zeta}{d\eta}\right|^{-1}=\frac{1}{\beta\left|\sqrt{1+4\epsilon\zeta/\beta}\right|},~{}&(x,y)\in\Omega_{o}\\\
1,~{}&\mathrm{for~{}air}\end{cases}$ (4)

where $n_{0}$ is the refractive index of the disk cavity in the OV space. In
this DLTC, the boundary shapes of the inner and outer layers are limaçon, and
the refractive index ratio at the interface between them is equivalent to
$n_{0}$ irrespective of boundary positions. It is noted that we use the
parameter $r_{o}\leq 1.0$ for the outer layer boundary in order to avoid the
local refractive index being smaller than $1.0$.

##  III Results

Resonant modes - We obtained the resonant modes in a DLTC in which both inner
and outer layer boundary shapes are limaçon. The DLTC corresponds to a
concentric circular layered cavity with $r=r_{i}$ and $r=r_{o}$ when
$\epsilon=0$, where there are whispering gallery modes (WGMs) circling the
inner layer boundary supported by the total internal reflection at the
interface between the inner and outer layers. We start from the WGM with mode
numbers $(m,l)=(13,1)$, where $m$ is the azimuthal mode number and $l$ the
radial mode number, of which the Q-factor is $1230$, in a concentric circular
layered cavity with refractive indices $n_{i}=1.8$ and $n_{o}=1.0$,
respectively, when $r_{i}=0.6$ and $\epsilon=0$. This concentric circular
layered cavity is the same as a circular cavity with a radius of $r_{i}$ and a
refractive index of $n_{i}$. As we increase $\epsilon$, the WGM in the
concentric circular layered cavity changes into a cWGM in a limaçon-shaped
DLTC.

Figure 3: Q-factors of the resonant modes in double-layered limaçon-shaped
transformation cavities as a function of $r_{o}$ when $\epsilon=0.05$ (black
circles), $\epsilon=0.1$ (green diamonds), $\epsilon=0.15$ (blue triangles),
and $\epsilon=0.2$ (red squares). The insets display selected far-field
intensity patterns. The orange dashed line represents the Q-factor of the
whispering gallery mode with $(m,l)=(13,1)$ in a circular cavity with
$n_{0}=1.8$.

A physical space consists of a gradient index cavity and the air, while a WT
space should consist of a gradient index cavity and a gradient index of an
infinite whole space which cannot be realized. As the size of the outer layer
increases, resonant modes in the TC become closer to those in a WT space; for
example, a cWGM in a TC approaches the corresponding WGM in a gradient index
cavity in a WT space, which is equivalent to a WGM in a homogeneous circular
cavity in an OV space. It should be noted that this crossover between a
gradient index cavity in a WT space and a TC is not limited to the cWGM, i.e.,
it is not the mode properties but rather the system characteristics.

Q-factors as a function of $r_{o}$ - To show the crossover from a TC to a
gradient index cavity in a WT space, we fix the deformation of the cavity
$\epsilon$ and the size of the inner layer $r_{i}$ but change the outer layer
boundary. In the case of a TC without an outer layer, $r_{o}=0.6$, as
$\epsilon$ increases, the Q-factors of the corresponding modes increase in the
region of $\epsilon<0.2$ because of the effect of increasing the refractive
index inside the cavity due to the distortion of the length scale by the
conformal mapping Ryu21 . As we increase the size of outer layer $r_{o}$, the
Q-factors of the resonant modes monotonically decrease if $r_{o}\lesssim 0.8$,
as shown in Fig. 3, because the outer layer effectively reduces the effect of
higher refractive index ratios at the inner layer boundary originating from
the distortion of the length scale inside the cavity by conformal mapping. If
$r_{o}\gtrsim 0.8$, the Q-factors show fluctuations as $r_{o}$ increases since
a sufficiently large outer layer plays the role of an additional resonator.
This kind of transition between monotonic and oscillatory behaviors depending
on the size of an additional resonator where the resonances are coupled to
open channels is well known in coupled microcavities Ryu06 . If we can
infinitely increase $r_{o}$ like in the case of a WT space, the Q-factor
approaches that of the corresponding WGM in a circular cavity, such as the
dashed line in Fig. 3, because the two resonant modes in OV and WT spaces are
equivalent. It is noted that the Q-factors of resonant modes in slightly
deformed cases such as $\epsilon=0.05$ can be smaller than that of the
corresponding WGM in a circular cavity in an oscillatory region.

Far-field intensity patterns - For the different $r_{o}$, the near-field
patterns inside the inner layer are not distinct from each other [Fig. 4 (a)],
i.e., they are cWGMs in the inner layers. On the other hand, the far-field
patterns show clearly different emission directionalities, as shown in Fig. 3.
In the case of a single transformation cavity with sufficiently large
$\epsilon$, the far-field patterns of cWGMs show bidirectional emissions,
which are explained by the universal tunneling mechanism related to the
boundary shapes and the inhomogeneous refractive index profiles determined by
the given conformal mapping Kim16 . As $r_{o}$ increases, the emission
directionalities of the far-field patterns change from bidirectional to
unidirectional via mixed-direction for large $\epsilon$ cases. Additionally,
the nearly isotropic emissions of the far-field patterns change into
unidirectional emissions for small $\epsilon$ cases such as $\epsilon=0.05$.

Figure 4: (a) Near-field and (b) far-field intensity patterns of a resonant
mode in a double-layered transformation cavity when $r_{i}=0.6$, $r_{o}=0.95$,
and $\epsilon=0.2$. The corresponding Husimi functions are shown for (c) the
inner and (d) the outer layer boundary, where $s$ is the boundary coordinate
in the reciprocal virtual space Kim18 . The yellow lines represent the
critical lines for total internal reflection, i.e., the ratios of the
refractive indices at the dielectric interfaces of the inner layer boundaries.
In the inset, straight (dotted) arrows inside the cavity denote the boundary
positions and incident (reflection) angles corresponding to $s_{1}$ and
$s_{2}$. The arrows outside cavity denote the boundary positions and
refractive angles of the dominant emissions. The red and blue arrows represent
the counter-clockwise ($\sin{\chi}>0$) and clockwise ($\sin{\chi}<0$)
components, respectively.

Near-field intensity patterns and Husimi functions - The remarkable changes of
the far-field intensity patterns as $r_{o}$ increases are caused by the
modification of the outer layer boundary, not by the significant changes of
the near-field intensity patterns around the perimeter of the inner layer.
Figure 4 (c) and (d) show Husimi functions at the inner and the outer layer
boundary, respectively, of cWGMs forming along the inner layer boundary. The
Husimi function is a phase space representation of intracavity wave intensity
and is defined at the dielectric interfaces by the overlap of the boundary
wave function of a resonant mode with a Gaussian wave packet on the cavity
boundary in a phase space Hen03 ; Lee05 ; Kim18 . In a conventional single TC,
the cWGM has varying tunneling emission according to cavity boundary position
since the ratio of the refractive indices at the dielectric interfaces of the
cavity boundaries is given as a function of boundary position. Tunneling
emission mainly occurs at the position where the ratio is minimum, and light
is emitted along the tangent direction at this position. In the DLTC, however,
the cWGM exhibits a nearly isotropic tunneling emission at the inner layer
boundary since the ratio of the refractive indices at the dielectric interface
between the inner and outer layers is constant, $n_{0}$, independent of
boundary position as shown in Fig. 4 (c). In Fig. 4 (a), the modes are better
confined in the left part of the outer layer because of the large contrast
between the refractive indices of the outer layer and air. A localized high
intensity of the Husimi function of the resonant mode above the critical line
for total internal reflection around $s_{1}$ of Fig. 4 (d) confirms well-
confined light. In contrast, the modes are emitted through the right part of
the outer layer boundary because the ratio of the refractive index between the
outer layer and air at that region is small. The second highest intensity of
the Husimi function around $s_{2}$ mainly locates below the critical line for
total internal reflection, and this produces the refractive emission in the
right part of the outer layer boundary shown in the inset of Fig. 4 (d). This
mechanism achieves unidirectional emission in the DLTC.

Figure 5: Q-factors of the resonant modes in double-layered limaçon-shaped
transformation cavities as a function of $\epsilon$ for $r_{o}=0.6$ (black
circles, transformation cavity without an outer cavity), $0.7$ (red squares),
$0.8$ (green diamonds), $0.9$ (blue triangles), and $1.0$ (orange inverted
triangles) when $r_{i}=0.6$. The black dashed line represents the Q-factor of
the whispering gallery mode with $(m,l)=(13,1)$ in a circular cavity with
$n_{0}=1.8$. The insets depict the boundary shapes and refractive index
profiles of single and double-layered limaçon-shaped transformation cavities.
The color scales of the refractive index profiles are the same as in Fig. 2.

Q-factors as a function of $\epsilon$ - Finally, we study the Q-factors of the
modes as a function of $\epsilon$ for different $r_{o}$ when $r_{i}=0.6$ (Fig.
5). In the case of a TC without an outer cavity ($r_{o}=0.6$), as $\epsilon$
increases, the Q-factors of the corresponding modes increase as discussed
above. In some cases, if $\epsilon$ increases further, the Q-factors show a
decreasing trend Kim16 . If we can infinitely increase $r_{o}$, like in the
case of conventional TO, the Q-factor does not depend on the change of
$\epsilon$ such as the dashed line in Fig. 5 because the two resonant modes in
OV and WT spaces are equivalent. If we increase $r_{o}$, the rate of Q-factor
increase decreases, and the Q-factors are mostly maintained regardless of
$\epsilon$ when $r_{o}=1.0$.

##  IV Summary

As an alternative model of a transformation cavity having a gradually varying
refractive index profile in the interior region of the cavity, a double-
layered transformation cavity consisting of inner and outer layers with
gradually varying refractive index profiles has been proposed. We have studied
the Q-factors, near- and far-field intensity patterns, and Husimi functions of
the resonant modes in double-layered transformation cavities. As the outer
layer becomes larger, the Q-factors and near-field intensity patterns of the
resonant modes in the newly proposed cavities become closer to those in a
cavity in the wholly transformed space of transformation optics. As the outer
layer boundary is changed without modifications to the inner layer, the far-
field patterns of the modes in the cavities can be controlled, while the near-
field patterns in the inner layers do not change. We expect that such a
control scheme of the emission directionality of the resonant modes without
significant deformation of the mode pattern only by modification of the outer
layer without changing the inner layer paves the way to design a versatile
cavity through the suggested double-layered transformation cavity. As a
similar example, hybrid optical microcavities comprising an inner
transformation cavity and an outer homogeneous deformed cavity have recently
been proposed to extract unidirectional narrow-beam emission from high-Q
whispering gallery modes Lim21 .

## acknowledgments

The author thanks J. Cho, I. Kim, S. Rim, and M. Choi for helpful discussions.
We acknowledge financial support from the Institute for Basic Science in the
Republic of Korea through the project IBS-R024-D1.

## References

  * (1) U. Leonhardt, Science 312, 1777–1780 (2006).
  * (2) J. B. Pendry, D. Schurig, and D. R. Smith, Science 312, 1780–1782 (2006).
  * (3) E. J. Post, Formal Structure of Electromagnetics, (Dover books on physics, New York, 1962).
  * (4) J. D. Jackson, Classical Electrodynamics, 2nd Edition, (John Wiley and Sons, New York, 1975).
  * (5) H. Chen, C. T. Chan, and P. Sheng, Nature Mater. 9, 387–396 (2010).
  * (6) A. V. Kildishev and V. M. Shalaev, Phys. Usp. 54, 53 (2011).
  * (7) Y. Liu and X. Zhang, Nanoscale 4, 5277–5292 (2012).
  * (8) D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, Science 314, 977–980 (2006).
  * (9) W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, Nat. Photon. 1, 224–227 (2007).
  * (10) N. Landy and D. R. Smith, Nature Mater. 12, 25–28 (2013).
  * (11) W. X. Jiang, T. J. Cui, Q. Cheng, J. Y. Chin, X. M. Yang, R. Liu, and D. R. Smith, Appl. Phys. Lett. 92, 264101 (2008).
  * (12) M. M. Sadeghi, L. Xu, H. Nadgaran, and H. Chen, Sci. Rep. 5, 11015 (2015).
  * (13) D. A. Genov, S. Zhang, and X. Zhang, Nature Phys. 5, 687–692 (2009).
  * (14) H. Chen, R. X. Miao, and M. Li, Opt. Express 18, 15183–15188 (2010).
  * (15) S. Dehdashti, H. Wang, Y. Jiang, Z. Xu, H. Chen, Prog. Electromagn. Res. 154, 181–193 (2015).
  * (16) D. A. Roberts, M. Rahm, J. B. Pendry, and D. R. Smith, Appl. Phys. Lett. 93, 251111 (2008).
  * (17) J. Zhang, Y. Luo, H. Chen, J. Huangfu, B. I. Wu, L. Ran, and J. A. Kong, Opt. Express 17, 6203-6208 (2009).
  * (18) P. H. Tichit, S. N. Burokur, and A. de Lustrac, Opt. Express 18, 767–772 (2010).
  * (19) V. G. Veselago, Sov. Phys. Usp. 10, 509 (1968).
  * (20) R. A. Shelby, D. R. Smith, and S. Schultz, Science 292, 77 (2001).
  * (21) Y. Y. Zhao, Y. L. Zhang, M. L. Zheng, X. Z. Dong, X. M. Duan, and Z. S. Zhao, Laser Photon. Rev. 10, 665–672 (2016).
  * (22) Y. Kim, S.-Y. Lee, J.-W. Ryu, I. Kim, J.-H. Han, H.-S. Tae, M. Choi, and B. Min, Nat. Photonics 10, 647-652 (2016).
  * (23) J.-W. Ryu, J. Cho, I. Kim, and M. Choi, Sci. Rep. 9, 8506 (2019).
  * (24) S.-J. Park, I. Kim, J. Cho, Y. Kim, and M. Choi, Opt. Express 27, 16320–16328 (2019).
  * (25) J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, IEEE Trans. Microwave Theory Tech. 47, 2075–2084 (1999).
  * (26) D. R. Smith, J. W. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, Phys. Rev. Lett. 84, 4184–4187 (2000).
  * (27) N. Engheta and R. W. Ziolkowski, Metamaterials: Physics and Engineering Explorations. (John Wiley & Sons, 2006).
  * (28) D. R. Smith, J. B. Pendry, and M. C. K. Wiltshire, Science 305, 788-792 (2004).
  * (29) J. Valentine, J. Li, T. Zentgraf, G. Bartal, and X. Zhang, Nat. Mater. 8, 568–571 (2009).
  * (30) L. H. Gabrielli, J. Cardenas, C. B. Poitras, and M. Lipson, Nat. Photonics 3, 461–463 (2009).
  * (31) B. Vasić, G. Isić, R. Gajić, and K. Hingerl, Opt. Exp. 18, 20321–20333 (2010).
  * (32) H. Gao, B. Zhang, S. G. Johnson, and G. Barbastathis, Opt. Exp. 20, 1617-1628 (2012).
  * (33) J.-W. Ryu, J. Cho, S.-Y. Lee, Y. Kim, S.-J. Park, S. Rim, M. Choi, and I. Kim, Sci. Rep. 9, 19684 (2019).
  * (34) I. Kim, J. Cho, Y. Kim, B. Min, J.-W. Ryu, S. Rim, and M. Choi, Opt. Exp. 26, 6851-6859 (2018).
  * (35) M. Robnik, J. Phys. A 16, 3971 (1983).
  * (36) J. Wiersig and M. Hentschel, Phys. Rev. Lett. 100, 033901 (2008).
  * (37) J.-W. Ryu and J. Cho, Opt. Lett. 46, 1900-1903 (2021).
  * (38) J.-W. Ryu, S.-Y. Lee, C.-M. Kim, and Y.-J. Park, Phys. Rev. A 74, 013804 (2006).
  * (39) M. Hentschel, H. Schomerus, and R. Schubert, Europhys. Lett. 62, 636 (2003).
  * (40) S.-Y. Lee, J.-W. Ryu, T.-Y. Kwon, S. Rim, and C.-M. Kim, Phys. Rev. A 72, 061801(R) (2005).
  * (41) J.-H. Lim, Y.-H. Lee, I. Kim, J. Cho, S. Rim, S.-J. Park, and M. Choi, Opt. Express, 29, 14736-14744 (2021).