# Atomic Bose-Einstein condensate in twisted-bilayer optical lattices

Zengming Meng  Liangwei Wang  Wei Han  Fangde Liu  Kai Wen  State Key
Laboratory of Quantum Optics and Quantum Optics Devices,
Institute of Opto-Electronics,
Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan,
Shanxi 030006, People’s Republic of China  Chao Gao  Department of Physics,
Zhejiang Normal University, Jinhua, 321004, People’s Republic of China
Pengjun Wang  State Key Laboratory of Quantum Optics and Quantum Optics
Devices,
Institute of Opto-Electronics,
Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan,
Shanxi 030006, People’s Republic of China  Cheng Chin  James Franck Institute,
Enrico Fermi Institute, Department of Physics, University of Chicago, Chicago,
Illinois 60637, USA  Jing Zhang<EMAIL_ADDRESS>State Key Laboratory of
Quantum Optics and Quantum Optics Devices,
Institute of Opto-Electronics,
Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan,
Shanxi 030006, People’s Republic of China  Hefei National Laboratory, Hefei,
Anhui 230088, P. R. China

###### Abstract

Observation of strong correlations and superconductivity in twisted-bilayer-
graphene Cao2018-1 ; Cao2018-2 ; Yankowitz1059 ; Lu2019  has stimulated
tremendous interest in fundamental and applied physics WANG2019100099 ;
Andrei2020 ; Balents2020 ; Kennes2021 . In this system, the superposition of
two twisted honeycomb lattices, generating a moir$\acute{\mathrm{e}}$ pattern,
is the key to the observed flat electronic bands, slow electron velocity and
large density of states PhysRevLett.99.256802 ; PhysRevB.81.161405 ;
Bistritzer12233 ; PhysRevB.85.195458 . Extension of the twisted-bilayer system
to new configurations is highly desired, which can open up exciting prospects
to investigate twistronics beyond bilayer graphene. Here we demonstrate a
quantum simulation of superfluid-to-Mott insulator transition in twisted-
bilayer square lattices based on atomic Bose-Einstein condensates loaded into
spin-dependent optical lattices. The lattices are made of two sets of laser
beams that independently address atoms in different spin states, which form
the synthetic dimension accommodating the two layers. The interlayer coupling
is highly controllable by microwave field, which enables the occurrence of
flat lowest band and novel correlated phases in the strong coupling limit. We
directly observe the spatial moir$\acute{\mathrm{e}}$ pattern and the momentum
diffraction, which confirm the presence of two forms of superfluid and a
modified superfluid to insulator transition in twisted-bilayer lattices. Our
scheme is generic and can be applied to different lattice geometries and for
both boson and fermion systems. This opens up a new direction for exploring
moir$\acute{\mathrm{e}}$ physics in ultracold atoms with highly controllable
optical lattices.

###### pacs:

34.20.Cf, 67.85.Hj, 03.75.Lm

Novel band structures in lattice systems often lead to new material functions
and discoveries. Twistronics, originating from the twisted-bilayer-graphene as
a tunable experimental platform Cao2018-1 ; Cao2018-2 ; Yankowitz1059 ; Lu2019
; WANG2019100099 ; Andrei2020 ; Balents2020 ; Kennes2021  has attracted broad
attention in recent years and launched intensive theoretical research. Here
overlaying two graphene layers with a small relative angle exhibit the rich
phase diagram, such as the coexistence of unconventional superconductivity and
correlated insulating phases Cao2018-2 ; Yankowitz1059 ; Lu2019 . In recent
years, many examples of twisted-bilayer are discovered with remarkable
physical properties not present in their untwisted counterparts. Recently,
photonic moir$\acute{\mathrm{e}}$ lattices are explored for their capabilities
in localizing and delocalizing light Wang2020 ; Huang2016 ; Fu2020  and
twisted-bilayer materials for engineering the photonic dispersion Hu2020 .

Ultracold atoms in optical lattices constitute an ideal platform to simulate
emerging many-body phenomena in condensed matter physics RevModPhys.80.885 ;
Lewenstein-book ; Windpassinger_2013 . Different optical lattice geometries
can be realized by interfering different sets of laser beams Soltan-Panahi2011
; Wirth2011 ; Tarruell2012 ; PhysRevLett.108.045305 ; Taiee1500854 ; Gall2021
. In particular, a scheme of simulating twisted-bilayer lattice is recently
proposed using two overlapping optical lattices PhysRevA.100.053604 ;
PhysRevLett.126.103201 . Other schemes for simulating bilayer heterostructures
have also been put forward Gra__2016 ; PhysRevLett.125.030504 . These schemes
are based on coherent coupling between spin states of atoms, which simulates
interlayer tunnelling along an artificial, synthetic dimension
PhysRevLett.108.133001 ; PhysRevLett.112.043001 ; Ozawa2019 .

Figure 1:  Simulation of twisted-bilayer systems based on atoms in spin-
dependent optical lattices. $\bf{a}$ Atoms are loaded into a single layer, 2D
pancake-like potential formed by a vertical optical lattice (green) in the
$z$-direction. Two sets of square optical lattices $V_{1}$ (purple) and
$V_{2}$ (blue) on the horizontal plane with a small relative angle
$\theta=5.21^{\circ}$ form a spin-dependent lattice potential and confine Rb
atoms in spin state $\left|1\right\rangle$ (up-arrows) and
$\left|2\right\rangle$ (down-arrows) independently. A magnetic field is
applied in the $x$-$y$ plane along the 45∘ diagonal of the $V_{2}$ lattice.
The lattice beams for $V_{1}$ and $V_{2}$ are set with opposite circular
polarization to generate the vector shift with the opposite sign. $\bf{b}$
Left panel: Sketch of the bilayer lattices in the synthetic dimension. The
interlayer tunnelling is controlled by an MW field. Right panel: Superimposed
lattice structure with the lattice constant $\lambda/2$ and much larger
moir$\acute{\mathrm{e}}$ length $\lambda_{mo}$. $\bf{c}$ Energy diagram of the
two ground Zeeman states $\left|1\right\rangle$ and $\left|2\right\rangle$ and
the associated lattice beams at the tune-out wavelengths $\lambda_{1}=790.02$
nm and $\lambda_{2}=788.28$ nm.

In this paper we demonstrate Bose-Einstein condensates (BEC) of Rubidium-87
(87Rb) atoms loaded into a pair of twisted-bilayer optical lattices. Two
overlapping lattices $V_{1}$ and $V_{2}$ are formed by interfering laser beams
at the specific “tune-out” wavelengths leblanc2007species ; arora2011tune ;
Wen2021  $\lambda_{1}$ and $\lambda_{2}$ with proper polarizations such that
atoms in spin state $\left|1\right\rangle\equiv\left|F=1,m_{F}=1\right\rangle$
and state $\left|2\right\rangle\equiv\left|F=2,m_{F}=0\right\rangle$ only
experience the lattice potential $V_{1}$ and $V_{2}$, respectively, see Fig.
1. Here $F$ and $m_{F}$ are the angular momentum and projection quantum
numbers in the 87Rb ground state manifold. Each set of the laser beams forms a
two-dimensional (2D) square lattice on the horizontal x-y plane and the twist
of the two lattices is realized by orienting the beams of different
wavelengths with a small relative angle $\theta=5.21^{\circ}$. The sample is
tightly confined in the vertical $z$-direction such that the sample is in the
quasi-2D regime (see Supplementary Material section I for details).

Figure 2:  Independent diffraction of atoms in different spin states by the
twisted-bilayer optical lattices. The optical lattice potential is applied to
the atomic BEC with a short duration of 4 $\mu$s. The images show diffraction
patterns of the atoms after 18 ms free space expansion. At the tune-out
wavelength $\lambda_{1}=790.02$ nm and $\lambda_{2}=788.28$ nm, atoms in state
$\left|1\right\rangle$ and $\left|2\right\rangle$, are diffracted by the
associated optical lattice $V_{1}$ and $V_{2}$, respectively.

The two spin states of 87Rb atoms constitute the synthetic dimension that
accommodates the two twisted layers of lattices $V_{1}$ and $V_{2}$. To
precisely determine the tune-out wavelengths $\lambda_{1}$ and $\lambda_{2}$
of the optical lattices $V_{1}$ and $V_{2}$, we measure the diffraction of
atoms by the optical lattices. The experimental sequence starts with an almost
pure BEC in a crossed-beam dipole trap. The atoms are prepared in one of the
two spin states and a short pulse of the lattice beams is applied. The lattice
potential induces Bragg diffraction of atoms to high momentum states. After
turning off the lattice beams, we image the diffracted atoms. The wavelengths
of the lattice beams are finely adjusted to the “tune-out” wavelengths such
that atoms in state $\left|1\right\rangle$ are only diffracted by the lattice
potential $V_{1}$ and not by the potential $V_{2}$ as shown in Fig. 2.
Similarly atoms in state $\left|2\right\rangle$ only experience the potential
$V_{2}$, but not $V_{1}$. By eliminating the cross-talks, we determined the
“tune-out” wavelengths to be $\lambda_{1}=790.02$ nm and $\lambda_{2}=788.28$
nm. Remarkably the lattice beams are circularly polarized to produce spatial
intensity modulation such that the lattice potentials are attractive to atoms
in both spin states (see Supplementary Material section II for details).

Figure 3:  Interlayer coupling in twisted-bilayer optical lattices. $\bf{a}$
MW spectrum of atoms in the twisted-bilayer optical lattices. Atoms in spin
state $\left|1\right\rangle$ are driven by MW to spin state
$\left|2\right\rangle$ in the presence of the lattice potential with depth of
0 (no lattice), 4 and 10$E_{r}$. Here $E_{r}=q_{r}^{2}/2m=h\times 3.67$ kHz is
the recoil energy, $q_{r}=\hbar k=h/\lambda$ is the recoil momentum, $m$ is
the atomic mass of 87Rb, and $\lambda$ is the wavelength of the lattice laser.
The MW pulse length of 530 $\mu$s corresponds to a $\pi$ pulse in the absence
of the lattice potential. $\bf{b}$ Lattice band structure for the two spin
states calculated with the lattice depth 4$E_{r}$. The MW field drives atoms
from the $s$-band of state $\left|1\right\rangle$, labelled as $|1,S\rangle$
to $s$-, $p$\- and $d$-bands of state $\left|2\right\rangle$ with different
detuning $\Delta$. $\bf{c}$ and $\bf{d}$ Starting with all atoms in
$|1,S\rangle$, population in state $\left|2\right\rangle$ is measured in the
twisted-bilayer lattices at 4$E_{r}$ after the MW pulse that drives the atoms
to $|2,S\rangle$ with the detuning $\Delta=-0.9$ kHz, or to $|2,P\rangle$ with
detuning $\Delta=15.08$ kHz. Fits in panel $\bf{c}$ show an interlayer
coupling frequency of $\Omega_{\mathrm{R}}=$2$\pi\times 893$ Hz and a decay
rate of 1200/s. Lines in panel $\bf{d}$ are guides to the eye. Each point is
based on three or more measurements and error bars show the standard
deviations of the mean.

Experimentally intralayer hoppings $t_{1}$ and $t_{2}$ between lattice sites
are controlled by the depth of the optical lattices $V_{1}$ and $V_{2}$;
interlayer hopping $\Omega_{\mathrm{R}}$, on the other hand, is independently
induced by microwaves (MW) that couple the two spin states. Starting with
atoms in state $\left|1\right\rangle$ in the dipole trap, for example, the MW
spectrum displays a single narrow peak when atoms are driven to state
$\left|2\right\rangle$. By loading the atoms into the twisted-bilayer optical
lattices, the spectrum displays multiple peaks. The peaks correspond to
transitions from atoms in the ground band of lattice $V_{1}$, which we label
$|1,S\rangle$, to different Bloch bands of lattice $V_{2}$, which we label
$|2,S\rangle$, $|2,P\rangle$, $|2,D\rangle$ and so on, see Figs. 3a and 3b.
The peak locations agree with the calculated energies of the $s$-, $p$\- and
$d$-bands in lattice $V_{2}$. Remarkably, the multi-peak structure supports
that atoms in different spin states are confined in different lattices. If
atoms are loaded into a spin-independent lattice, only a single narrow peak
shows up in the spectrum, which belongs to the $|1,S\rangle$ to $|2,S\rangle$
transition. This is because MW transitions between different Bloch bands are
negligible in spin-independent lattices. In the twisted optical lattice, the
transitions from $s$-band of state $\left|1\right\rangle$ to other bands of
state $\left|2\right\rangle$ are allowed. In the presence of the twisted
bilayer lattices, the transitions are broadened since the two spin states
experience different trap potentials, which induce fast dephasing. Moreover,
the on-site interactions increase in deeper lattice potential, resulting in
faster decay from high bands to lower bands and thus broader spectral lines.
Our observation supports MW as a versatile and powerful tool to induce
interlayer hopping between the two twisted layers in the synthetic (spin)
dimension.

To quantify the interlayer hopping energy, we measure the time evolution of
the population in state $\left|2\right\rangle$. We observe a coherent
oscillation at detuning $\Delta=-0.9$ kHz, which corresponds to the transition
from $|1,S\rangle$ to $|2,S\rangle$, see Fig. 3c. The interlayer coupling
strength can be determined from the oscillation frequencies. In our
experiment, the coupling strength is tunable up to 1$E_{r}$, which exceeds
that in typical twisted-bilayer graphene systems. On the other hand, coupling
to the $p$-band $|2,P\rangle$ leads to faster decay likely due to collisional
relaxation to the lower $s$-band, see Fig. 3d. In the following, we will focus
on atoms in the twisted-bilayer optical lattices with MW-induced coupling
between the $s$-bands of the two layers.

A key signature of atoms in the twisted-bilayer optical lattice is the
moir$\acute{\mathrm{e}}$ lattice with a period

$\displaystyle\lambda_{mo}=\frac{a}{2\sin\theta/2},$ (1)

which, for the lattice constant $a=395$ nm and twist angle
$\theta=5.21^{\circ}$, amounts to $\lambda_{mo}=4.35$ $\mu$m. The large
moir$\acute{\mathrm{e}}$ period gives rise to a mini-Brillouin zone in the
momentum space, which is expected to generate the flat bands and strongly
correlated states Cao2018-1 ; Cao2018-2 ; Yankowitz1059 ; Lu2019 ;
WANG2019100099 ; Andrei2020 ; Balents2020 ; Kennes2021 ; PhysRevLett.99.256802
; PhysRevB.81.161405 ; Bistritzer12233 ; PhysRevB.85.195458 . Notably _in-
situ_ moir$\acute{\mathrm{e}}$ pattern is also observed in 1D lattices with
two lattice constants PhysRevX.9.021001 . To identify the
moir$\acute{\mathrm{e}}$ length scale in our system, we employ _in situ_
absorption imaging to visualize the moir$\acute{\mathrm{e}}$ pattern, see Fig.
4. Here we first load the atoms in state $\left|1\right\rangle$ into the
lowest $s$-band of lattice $V_{1}$ and then ramp up the MW field with detuning
$\Delta=-0.9$ kHz to drive the transition from $|1,S\rangle$ to $|2,S\rangle$.
We then _in situ_ image the atoms in state $\left|2\right\rangle$.
Moir$\acute{\mathrm{e}}$ patterns in one and two dimensions are observed, and
the moir$\acute{\mathrm{e}}$ period is measured to be 4.35 $\mu$m consistent
with expectation, see Figs. 4a-f. Note that the primary optical lattice
spacing $a=395$ nm is indiscernible with our imaging optics.

We also examine the quantum state of atoms in the bilayer twisted lattices by
analyzing their momentum space distribution. After loading a BEC into the
bilayer lattice of 4$E_{r}$ in the presence of resonant MW transition, we hold
for some time and then perform the time-of-flight measurement, see Figs. 4g
and 4h. Two sets of diffractions manifest, which correspond to the primary
lattice momentum $\pi/a$ and the much smaller moir$\acute{\mathrm{e}}$
momentum $\pi/\lambda_{mo}$. The high contrast of both sets of diffraction
pattern suggests that the atoms remain in the superfluid phase with phase
coherence extending beyond the moir$\acute{\mathrm{e}}$ length scale. In
particular, the contrasts of the moir$\acute{\mathrm{e}}$ pattern in real and
momentum space persist over 40 ms, see Fig. 4i, from which we conclude that
the atoms maintain in the superfluid phase in the twisted-bilayer lattices.

Figure 4:  Moir$\acute{\mathrm{e}}$ pattern and superfluid ground state in
twisted-bilayer optical lattices. $\bf{a}$ and $\bf{c}$
Moir$\acute{\mathrm{e}}$ pattern of atoms in one dimensional bilayer optical
lattices in the $x$\- and $y$\- direction. $\bf{b}$ and $\bf{d}$ Plot of the
corresponding optical lattice potential. Note that the primary optical lattice
$\lambda/2$ is indiscernible in the _in situ_ images. In all experiments, an
MW field is applied to couple the $s$-band of state $\left|1\right\rangle$ and
$s$-band of state $\left|2\right\rangle$. The atoms in state
$\left|2\right\rangle$ are measured. The elliptical shape of the atomic cloud
is induced by a little asymmetric harmonic trap potential in x-y plane.
$\bf{e}$ and $\bf{f}$ Experimental and theoretical moir$\acute{\mathrm{e}}$
pattern of atoms in the presence of optical lattices in both directions.
$\bf{g}$ and $\bf{h}$ Time-of-flight (TOF) images with 18 ms from the
experiment and calculation show diffraction peaks associated with the primary
lattice constant $\lambda/2$ and moir$\acute{\mathrm{e}}$ length
$\lambda_{mo}$. The enlarged picture of the center part of $\bf{g}$ is
obtained by the imaging system with higher magnification. $\bf{i}$ The
contrast of moir$\acute{\mathrm{e}}$ pattern in real space and the contrast of
diffraction pattern for different hold times. Here, the lattice depths are
$V_{1x}=V_{2x}=V_{1y}=V_{2y}=4E_{r}$ ($U/t=1.67$,
$\Omega_{\mathrm{R}}/t=2.07$, $U$ is the on-site interaction). The real-space
and momentum-space distributions of $\bf{f}$ and $\bf{h}$ are theoretically
calculated by solving the mean-field ground states according to the Gross-
Pitaevskii equations (see Supplementary Material section IV for details). The
contrast in the real space is defined as
$(S_{max}-S_{min})/(S_{max}+S_{min})$, where $S_{max}$ and $S_{min}$ are the
maximum and minimum atomic density of the moir$\acute{\mathrm{e}}$ fringes.
The contrast in the momentum space is defined as
$(P^{av}_{max}-P_{min})/(P^{av}_{max}+P_{min})$, where,
$P^{av}_{max}=(P^{0}_{max}+P^{1}_{max})/2$ and $P_{min}$ are the average
maximum and minimum density of the diffraction pattern. Here $P^{0}_{max}$ is
zero momentum component and $P^{1}_{max}$ is the moir$\acute{\mathrm{e}}$
component near the zero momentum. Error bars show the standard deviation of
the mean.

Theoretically, depending on the twist angle $\theta$, the superimposed
twisted-bilayer lattice can yield either a periodic potential with supercells
that supports a delocalized ground state or a quasi-periodic one that supports
a localized ground state in the absence of interactions. In fact, only
specific twist angles give rise to periodic lattice potentials. For square
lattices, the twist angles that lead to commensurate superlattice should
satisfy $\theta=2\arctan(\bar{m}/\bar{n})$, where $\bar{m}$ and $\bar{n}$ are
integers Wang2020 . The twist angle $\theta=5.21^{\circ}$ used in our work is
close to the commensurate angle $\theta=2\arctan(1/22)\approx 5.205^{\circ}$,
and the period of the supercell is given by $22a=2\lambda_{mo}$ (see
Supplementary Material section III for details). While our twist angle does
not exactly match the commensurate angle $\theta=2\arctan(1/22)\approx
5.205^{\circ}$, the small difference can not be distinguished in a finite size
sample due to repulsive interactions (see Supplementary Material section IV
for details). In particular, the spatial moir$\acute{\mathrm{e}}$ period
remains a clear observable in our experiment because of the finite chemical
potential of our atomic superfluid. The persistence of the spatial and
momentum periodicity of the sample in the twisted-bilayer lattice supports the
superfluid as the ground state of the system.

Compared with electronic materials, where the flat band is investigated
frequently near the Fermi surface, we can also explore flatband physics with
bosons condensed in the lowest band. In our system, when interlayer coupling
increases, the long-wavelength moir$\acute{\mathrm{e}}$ potential becomes
deeper, so atoms in the lowest band are isolated at a larger spatial scale
(moir$\acute{\mathrm{e}}$ wavelength), which flattens the ground band and
enhances the localization of the atoms. In the large interlayer coupling
limit, the system can be regarded as a single-layer (single-component)
experiencing a twisted optical lattice (see Supplementary Material section III
for details). The single-layer system with a twisted optical lattice admits a
flatband structure in the ground band, which has also been studied
experimentally in photonic systems Wang2020 ; Huang2016 ; Fu2020 . The easily-
tuned intra- and interlayer couplings in our system offer the added advantages
of seeking novel quantum phases and phase transitions with cold atoms.

By varying the depth of optical lattices and interlayer coupling, we find
several distinct quantum phases, including superfluid (SF), superfluid with
only short-range coherence (SF-II), Mott Insulator (MI) and insulator (I),
(see Fig. 5a and Supplementary Material section IV). These phases can be
distinguished by the phase coherence and real-space density correlations.
Remarkably the SF-II phase emerges with finite interlayer coupling around the
transition from a regular SF to an insulator. The spatial range of phase
coherence is the key to distinguishing the two SF phases: while an SF supports
long-range phase coherence PhysRevA.72.053606 , the SF-II phase maintains the
coherence only up to the moir$\acute{\mathrm{e}}$ length scale. In addition,
the SF-II phase supports the moir$\acute{\mathrm{e}}$ pattern in the real-
space. Theoretically, SF-II is the phase with superfluid domains embedded in a
gapped insulator, induced by the interlayer coupling. Finally, the insulator
phases I and MI can be identified by the disappearance of spatial coherence at
all scales and integer fillings of all the sites. While the MI has uniform
atom density with weak interlayer coupling, the I phase features
moir$\acute{\mathrm{e}}$ pattern due to stronger interlayer coupling.

Figure 5:  Phase transition for the twisted-bilayer optical lattice. $\bf{a}$
Phase diagram (see Supplementary Material section IV for details), where SF,
SF-II, MI, and I refer to superfluid, superfluid only with short-range
coherence, Mott insulator, and insulator. The solid curves denote the
calculated phase boundaries with mean-field at zero temperature. The dots are
experimental measurements of the phase boundaries. The pink circles and blue
squares denote the loss of the coherence at the moir$\acute{\mathrm{e}}$
length scale and the length scale of the primary lattice respectively. The red
diamonds denote the appearance of the moir$\acute{\mathrm{e}}$ pattern in the
real-space (in situ images) as increases the interlayer coupling. The error
bars indicate the experimental uncertainties in determining the phase
transition. $\bf{b}$ Visibility curves for the moir$\acute{\mathrm{e}}$ and
primary lattice momentum components as a function of lattice depth (Path I in
$\bf{a}$). A sequential loss of phase coherence appears at the
moir$\acute{\mathrm{e}}$ momentum and the primary lattice momentum. The
intermediate regime indicates SF-II phase. The interlayer coupling frequency
is $\Omega_{\mathrm{R}}=0.24E_{r}$. $\bf{c}$ Visibility curves for the
moir$\acute{\mathrm{e}}$ momentum component as a function of interlayer
coupling with the lattice depth $8E_{r}$ (Path II in $\bf{a}$) and $14E_{r}$
(Path III in $\bf{a}$) respectively. The solid lines in $\bf{b}$ and $\bf{c}$
are fitted from the experimental data and only guide the eye.

In the experiment, we measure the phase coherence from the momentum-space
diffraction peaks in the TOF images and directly probe the
moir$\acute{\mathrm{e}}$ pattern by _in situ_ imaging following the
measurement method as shown in Fig. 4. The measurement of the phase boundaries
is shown in Fig. 5a. We employ three independent paths to study these phases.
Path I, we fix the interlayer coupling strength at a small value
$\Omega_{R}=0.24E_{r}$ and increase the lattice depth. The phase transition
from SF to MI and across SF-II is shown in the time-of-flight images, see Fig.
5b. Here the diffraction peaks at the moir$\acute{\mathrm{e}}$ momenta
disappear first before the disappearance of the primary lattice. The
intermediate regime indicates the SF-II phase, where the
moir$\acute{\mathrm{e}}$-scale long-range correlation is destroyed while a
short-range coherence remains, at the same time, the density correlations of
moir$\acute{\mathrm{e}}$ pattern appears in the real-space. Path II, we fix
lattice depth in the SF region and increase the interlayer coupling. The
diffraction peaks at the moir$\acute{\mathrm{e}}$ momenta persist with high
contrasts. However, in Path III, when the depth of optical lattices is fixed
at the MI region and the interlayer coupling increases, the visibility at the
moir$\acute{\mathrm{e}}$ momenta presents the threshold behavior and emerges
at $\Omega_{R}>0.5E_{r}$, see Fig. 5c. These observations are qualitatively
consistent with the theoretical expectation and demonstrate that the
interlayer coupling can induce a reentrant transition from MI to SF across SF-
II. One may understand such rich transitions from the fact that the interplay
between the interlayer coupling and interactions tends to localize the bosons,
primarily in the moir$\acute{\mathrm{e}}$ length scale.

This work provides a preliminary physical insight into the quantum phase
transition between SF and SF-II ( MI and SF-II or MI and I) and offers the
possibility to study the complex phases due to the presence of quasi-disorder
induced by large interlayer coupling and strong interaction, such as Bose
glass insulator, resembling that in disordered bosonic systems PhysRevB.40.546
; PhysRevLett.67.2307 ; PhysRevB.53.2691 . These complex phases are worth
further investigating in the future.

The present work focuses on the realization and the ground state properties of
atoms in the twisted-bilayer optical square lattice. Our success in loading a
superfluid into the bilayer lattice demonstrates a new versatile platform to
explore moir$\acute{\mathrm{e}}$ physics and the associated superfluidity in a
quantum many-body system. Beyond the tunable twist angle, the cold atom
platform offers remarkable controls such as different lattice depths and
interlayer coupling in different layers.

Furthermore, the twisted-bilayer square lattice closely connects to the
physics of heterostructures of twisted atomically thin semiconductors
PhysRevResearch.1.033076 ; Kennes2020 ; Kennes2021 . At the same time, our
experiment can in principle be extended to multi-layer lattice where the
interlayer couplings can be independently induced by MW and radio-frequencies.
Replacing the MW by optical Raman transitions, the interlayer coupling can be
spatial dependence, which can support topological ground states. Finally, our
optical lattice scheme can be applied to confine fermionic atoms in bilayer
hexagonal lattice, which faithfully simulates electrons in a bilayer graphene,
and may offer insight into the emergence of superconductivity in the strongly
correlated, flatband regime.

## DATA AVAILABILITY

All data generated or analysed during this study are included in this
published article. Additional data are also available from the corresponding
authors upon reasonable request.

###### Acknowledgements.

This research is supported by Innovation Program for Quantum Science and
Technology (Grant No. 2021ZD0302003), National Key Research and Development
Program of China (2016YFA0301602, 2018YFA0307601, 2022YFA1404101), NSFC (Grant
No. 12034011, 12022406, 11804203), the Fund for Shanxi “1331 Project” Key
Subjects Construction and Tencent (Xplorer Prize). CC acknowledges support by
the National Science Foundation (Grant No. PHY-2103542), and the Army Research
Office STIR (Grant W911NF2110108).

## Author contributions

J.Z. conceived the idea and performed the experimental designs. L.W., Z.M.,
F.L., K.W., P.W. and J.Z. performed the experiments. C.C., Z.M., L.W., F.L.,
W.H. and J.Z. analysed the data and all authors discussed the results. W.H.,
C.G. and J.Z. performed the simulation. Z.M. plotted the figures. J.Z. and
C.C. wrote the manuscript. All authors interpreted the results and reviewed
the manuscript. J.Z. designed and supervised the project.

## Competing interests

The authors declare no competing financial interests.

## References

  * (1) Cao, Y. _et al._ Correlated insulator behaviour at half-filling in magic-angle graphene superlattices.  _Nature_ 556, 80–84 (2018).
  * (2) Cao, Y. _et al._ Unconventional superconductivity in magic-angle graphene superlattices.  _Nature_ 556, 43–50 (2018).
  * (3) Yankowitz, M. _et al._ Tuning superconductivity in twisted bilayer graphene 363, 1059–1064 (2019).
  * (4) Lu, X. _et al._ Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene.  _Nature_ 574, 653–657 (2019).
  * (5) Wang, J., Mu, X., Wang, L. & Sun, M.  Properties and applications of new superlattice: twisted bilayer graphene.  _Mater. Today Phys._ 9, 100099 (2019).
  * (6) Andrei, E. Y. & MacDonald, A. H.  Graphene bilayers with a twist.  _Nat. Mater._ 19, 1265–1275 (2020).
  * (7) Balents, L., Dean, C. R., Efetov, D. K. & Young, A. F.  Superconductivity and strong correlations in moir$\acute{\mathrm{e}}$ flat bands.  _Nat. Phys._ 16, 725–733 (2020).
  * (8) Kennes, D. M. _et al._ Moir$\acute{\mathrm{e}}$ heterostructures as a condensed-matter quantum simulator.  _Nat. Phys._ 17, 155–163 (2021).
  * (9) Lopes dos Santos, J. M. B., Peres, N. M. R. & Castro Neto, A. H.  Graphene bilayer with a twist: electronic structure.  _Phys. Rev. Lett._ 99, 256802 (2007).
  * (10) Mele, E. J.  Commensuration and interlayer coherence in twisted bilayer graphene.  _Phys. Rev. B_ 81, 161405 (2010).
  * (11) Bistritzer, R. & MacDonald, A. H.  Moiré bands in twisted double-layer graphene 108, 12233–12237 (2011).
  * (12) Moon, P. & Koshino, M.  Energy spectrum and quantum hall effect in twisted bilayer graphene.  _Phys. Rev. B_ 85, 195458 (2012).
  * (13) Wang, P. _et al._ Localization and delocalization of light in photonic moiré lattices.  _Nature_ 577, 42–46 (2020).
  * (14) Huang, C. _et al._ Localization-delocalization wavepacket transition in pythagorean aperiodic potentials.  _Sci. Rep._ 6, 32546 (2016).
  * (15) Fu, Q. _et al._ Optical soliton formation controlled by angle twisting in photonic moiré lattices.  _Nat. Photonics_ 14, 663–668 (2020).
  * (16) Hu, G. _et al._ Topological polaritons and photonic magic angles in twisted $\alpha$-moo3 bilayers.  _Nature_ 582, 209–213 (2020).
  * (17) Bloch, I., Dalibard, J. & Zwerger, W.  Many-body physics with ultracold gases.  _Rev. Mod. Phys._ 80, 885–964 (2008).
  * (18) Lewenstein, M., Sanpera, A. & Ahufinger, V.  Ultracold atoms in optical lattices: simulating quantum many-body systems (Oxford University Press, 2012).
  * (19) Windpassinger, P. & Sengstock, K.  Engineering novel optical lattices.  _Rep. Prog. Phys._ 76, 086401 (2013).
  * (20) Soltan-Panahi, P. _et al._ Multi-component quantum gases in spin-dependent hexagonal lattices.  _Nat. Phys._ 7, 434–440 (2011).
  * (21) Wirth, G., $\ddot{O}$lschl$\ddot{a}$ger, M. & Hemmerich, A.  Evidence for orbital superfluidity in the p-band of a bipartite optical square lattice.  _Nat. Phys._ 7, 147–153 (2011).
  * (22) Tarruell, L., Greif, D., Uehlinger, T., Jotzu, G. & Esslinger, T.  Creating, moving and merging dirac points with a fermi gas in a tunable honeycomb lattice.  _Nature_ 483, 302–305 (2012).
  * (23) Jo, G.-B. _et al._ Ultracold atoms in a tunable optical kagome lattice.  _Phys. Rev. Lett._ 108, 045305 (2012).
  * (24) Taie, S. _et al._ Coherent driving and freezing of bosonic matter wave in an optical lieb lattice 1, e1500854 (2015).
  * (25) Gall, M., Wurz, N., Samland, J., Chan, C. F. & K$\ddot{o}$hl, M.  Competing magnetic orders in a bilayer hubbard model with ultracold atoms.  _Nature_ 589, 40 (2021).
  * (26) González-Tudela, A. & Cirac, J. I.  Cold atoms in twisted-bilayer optical potentials.  _Phys. Rev. A_ 100, 053604 (2019).
  * (27) Luo, X.-W. & Zhang, C.  Spin-twisted optical lattices: tunable flat bands and larkin-ovchinnikov superfluids.  _Phys. Rev. Lett._ 126, 103201 (2021).
  * (28) Graß, T., Chhajlany, R. W., Tarruell, L., Pellegrini, V. & Lewenstein, M.  Proximity effects in cold atom artificial graphene.  _2D Mater._ 4, 015039 (2016).
  * (29) Salamon, T. _et al._ Simulating twistronics without a twist.  _Phys. Rev. Lett._ 125, 030504 (2020).
  * (30) Boada, O., Celi, A., Latorre, J. I. & Lewenstein, M.  Quantum simulation of an extra dimension.  _Phys. Rev. Lett._ 108, 133001 (2012).
  * (31) Celi, A. _et al._ Synthetic gauge fields in synthetic dimensions.  _Phys. Rev. Lett._ 112, 043001 (2014).
  * (32) Ozawa, T. & Price, H. M.  Topological quantum matter in synthetic dimensions.  _Nat. Rev. Phys._ 1, 349–357 (2019).
  * (33) LeBlanc, L. J. & Thywissen, J. H.  Species-specific optical lattices.  _Phys. Rev. A_ 75, 053612 (2007).
  * (34) Arora, B., Safronova, M. S. & Clark, C. W.  Tune-out wavelengths of alkali-metal atoms and their applications.  _Phys. Rev. A_ 84, 043401 (2011).
  * (35) Wen, K. _et al._ Experimental study of tune-out wavelengths for spin-dependent optical lattice in 87rb bose-einstein condensation.  _J. Opt. Soc. Am. B_ 38, 3269 (2021).
  * (36) McDonald, M., Trisnadi, J., Yao, K.-X. & Chin, C.  Superresolution microscopy of cold atoms in an optical lattice.  _Phys. Rev. X_ 9, 021001 (2019).
  * (37) Gerbier, F. _et al._ Interference pattern and visibility of a mott insulator.  _Phys. Rev. A_ 72, 053606 (2005).
  * (38) Fisher, M. P. A., Weichman, P. B., Grinstein, G. & Fisher, D. S.  Boson localization and the superfluid-insulator transition.  _Phys. Rev. B_ 40, 546–570 (1989).
  * (39) Krauth, W., Trivedi, N. & Ceperley, D.  Superfluid-insulator transition in disordered boson systems.  _Phys. Rev. Lett._ 67, 2307–2310 (1991).
  * (40) Freericks, J. K. & Monien, H.  Strong-coupling expansions for the pure and disordered bose-hubbard model.  _Phys. Rev. B_ 53, 2691–2700 (1996).
  * (41) Kariyado, T. & Vishwanath, A.  Flat band in twisted bilayer bravais lattices.  _Phys. Rev. Res._ 1, 033076 (2019).
  * (42) Kennes, D. M., Xian, L., Claassen, M. & Rubio, A.  One-dimensional flat bands in twisted bilayer germanium selenide.  _Nat. Commun._ 11, 31124 (2020).

### Supplementary Material

I. Experimental setup

In our experiment, the ultracold 87Rb atoms in the $|F=2,m_{F}=2\rangle$ state
are prepared in the crossed optical dipole trap Xiong:10 (44). Forced
evaporation in the optical trap creates the BEC with up to $5\times 10^{5}$
atoms. The atoms can be transferred to the $|F=1,m_{F}=1\rangle$ state via a
rapid adiabatic passage induced by microwave transition. To load the atoms
into the 2D trap, a 532 nm laser beam is deflected by an acousto-optic
deflector (AOD) and then split into two beams with variable spacing adjusted
by the AOD. The two beams are focused onto the atoms with a 150 mm aspherical
lens. These beams interfere to form a standing wave in the vertical direction
with variable separation (Accordion lattice). This separation can be varied
from 12 $\mu m$ down to 3 $\mu m$. The advantage of variable spacing is that
we can load a 3D shaped cloud into a single layer of the 2D pancakes at
maximum separation and then compress the pancake adiabatically to reach a deep
2D regime. The maximum vertical confinement can reach more than 20 kHz and we
optimize it at 1 kHz to observe moir$\acute{\mathrm{e}}$ pattern and
superfluid of ultracold atoms.

The twisted-bilayer optical lattices are created by two sets of 2D square
lattice $V_{1}$ and $V_{2}$. A twisted angle of $\theta=5.21^{\circ}$ is set
between the two lattice potentials, namely, $V_{2}(\bf{r})$ $=V_{1}(S\bf{r})$,
$S=\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\\ \sin\theta&\cos\theta\\\
\end{array}\right)$. The optical lattices $V_{1}$ and $V_{2}$ are derived from
two CW Ti:Sapphire single frequency lasers (M Squared lasers SolsTiS and
Coherent MBR-110) respectively. Two lattice beams $V_{1x}$ and $V_{1y}$ of
$V_{1}$ are frequency-shifted +80 MHz and +95 MHz by two single-pass acousto-
optic modulators (AOMs), respectively. The same applies to the two lattice
beams $V_{2x}$ and $V_{2y}$ of lattice $V_{2}$. The four lattice beams are
coupled into polarization-maintaining single-mode fibers in order to improve
the stability of the beam pointing and achieve better beam-profile quality.
After the fibers, each lattice beam is focused by a lens and retroreflected by
a concave mirror. In order to generate the vector light shift, we use the same
circular polarization for two lattice beams to produce spatial intensity
modulation. In the experiment, we can determine and calibrate this angle by
measuring the intersection angle between two lasers and the
moir$\acute{\mathrm{e}}$ period from the _in situ_ images. The estimated
uncertainty of the two methods is about 0.05∘.

We use the MW field to couple the two spin states for manipulating the
interlayer coupling. The 6.8 GHz MW signal is amplified by a 10 W solid state
amplifier (Kuhne Electronic, KU PA 640720-10A). We place a circulator on the
output of the amplifier to reduce reflected power coming back to the
amplifier. The MW is emitted out to the atoms by a sawed-off waveguide, which
is placed outside of the high vacuum glass cell. We use MW cables to transfer
MW from the amplifier to the waveguide. With this MW power amplifier, we can
reach the maximum interlayer coupling strength of about 1.0$E_{r}$. It is
feasible to increase the interlayer coupling strength to about several $E_{r}$
by using an available higher power amplifier.

Our image system consists of an objective with a numerical aperture of
NA=0.69, working distance of 11 mm and effective focal length of 18 mm. A 900
mm lens after the objective leads to a magnification of 50 for _in situ_
imaging with an EMCCD (Andor iXon Ultra 897). We also employ a 200 mm (400 mm)
lens after the objective leads to a magnification of 11 (22) for the time-of-
flight absorption imaging with 18 ms. The atoms are detected by state-
selective absorptive imaging. Since we choose two ground hyperfine Zeeman
states of 87Rb $\left|F=2,m_{F}=0\right\rangle$ of the $F=2$, and
$\left|1,1\right\rangle$ of the $F=1$ hyperfine manifold as the two internal
spin states, we can fully resolve the population in each individual state. For
$\left|F=2,m_{F}=0\right\rangle$ state, a 50 $\mu$s long imaging pulse of
resonant light on the $F=2\rightarrow F^{\prime}=3$ D2 cycling transition is
used to detect the $\left|2\right\rangle$ atoms. In order to detect
$\left|F=1,m_{F}=1\right\rangle$ state, a resonant light pulse on the
$F=2\rightarrow F^{\prime}=3$ cycling transition is firstly used to remove the
$\left|2\right\rangle$ atoms and then a 50 $\mu$s long imaging pulse of
resonant light on the $F=2\rightarrow F^{\prime}=3$ is applied at the same
time with a repump light (resonant light $F=1\rightarrow F^{\prime}=2$) to
detect the $\left|1\right\rangle$ atoms.

When studying the superfluid-to-Mott insulator transition, we use the standard
method of interference pattern contrast (visibility) to reveal this transition
PhysRevA.72.053606 ; Greiner2002 (45). We first load the atoms in state
$\left|1\right\rangle$ into the lowest $s$-band of lattice $V_{1}$ by ramping
up $V_{1}$ and $V_{2}$ simultaneously with 30 ms, and then ramp up the MW
field with 10 ms to drive the transition from $|1,S\rangle$ to $|2,S\rangle$.
The atoms are detected by state-selective absorptive imaging with time-of-
flight of 18 ms after switching off all lattices and trapping light. In
experiment, we first check that BEC is still kept (see Fig. S1) as ramp up the
lattice $V_{1}$ (or $V_{2}$) to the higher lattice depth than $24E_{r}$ and
then ramp down again, which makes sure to perform the phase transition from SF
to MI successfully for the lattice V1 (or V2). When adding the interlayer
coupling between two spin states, and at the same time a quasi-disorder is
introduced, there are two more mechanics to make the system not completely
reversible. One is the finite coherent time between two spin states. When the
system is prepared initially in the spin down, the system will become the spin
mixture after the interlayer coupling is ramped back down. We define this
process as irreversibility. The other is that adiabaticity is broken down by a
quasi-period or disordered lattice, which induces not to completely remain in
the zero-momentum state after ramping the lattices back down.

Figure S1:  Coherence in the SF-MI transition. $\bf{a}$ The initial BEC in 2D
pancake-like potential. $\bf{b}$ Absorption image after atoms are released
abruptly from an optical lattice potential $V_{1}$ (or $V_{2}$) with a
potential depth $24E_{r}$. $\bf{c}$ Absorption image when the lattice is
ramped up to the lattice depth $24E_{r}$ and then ramp down to zero. The
images are obtained after 18 ms free space expansion.

II. Tune-out wavelength for twisted-bilayer optical lattices

Figure S2:  Determination of tune-out wavelengths. $\bf{a}$-$\bf{b}$ The
lattice depth $V_{1x}$ (blue) and $V_{1y}$ (red) as a function of wavelength
$\lambda$ for the two different hyperfine states $|F=1,{{m}_{F}}=1\rangle$ and
$|F=2,{{m}_{F}}=0\rangle$. The angles between $V_{1x}$, $V_{1y}$ and $B_{0}$
are $39.79^{\circ}$ and $50.21^{\circ}$ respectively. $\bf{c}$-$\bf{d}$ The
potential depth $V_{2x}$ (blue) and $V_{2y}$ (red) as a function of wavelength
$\lambda$ for the two different hyperfine states $|F=1,{{m}_{F}}=1\rangle$ and
$|F=2,{{m}_{F}}=0\rangle$. $\bf{e}$ Theoretical light shift of $V_{1x}$,
$V_{1y}$ for $|1,1\rangle$ and $|2,0\rangle$. $\bf{f}$ Theoretical lattice
depth of $V_{2x}$, $V_{2y}$ for $|1,1\rangle$ and $|2,0\rangle$. The bias
magnetic field of 10 Gauss is applied along the $45^{\circ}$ diagonal line of
the square lattice $V_{2}$.

AC Stark shift, or light shift is a light-induced change of energy level. For
alkali-metal atoms, the total ac Stark shift can be expressed in the
irreducible components (including scalar, vector, tensor components) of the
polarizability steck2007quantum (46):

$\displaystyle\Delta U=$ $\displaystyle\Delta
U\left(F,{{m}_{F}};\omega\right)$ (2) $\displaystyle=$
$\displaystyle-I\Big{[}\alpha_{~{}}^{\left(0\right)}\left(\omega\right)+{{\alpha}^{\left(1\right)}}\left(\omega\right)\left(\xi\hat{e}_{k}\cdot\hat{e}_{B}\right)\frac{{m}_{F}}{F}$
$\displaystyle+{{\alpha}^{\left(2\right)}}\left(\omega\right)\frac{3\cos^{2}\phi-1}{2}\frac{3{m}_{F}^{2}-F(F+1)}{F(2F-1)}\Big{]},$

where $F$ is the total atomic angular momentum, $m_{F}$ is the magnetic
quantum number, $\omega$ is the laser frequency, $I$ is the laser field
intensity, $\xi$ is light ellipticity, $\hat{e}_{k}$ and $\hat{e}_{B}$ are
unit vectors along the light wave-vector and magnetic field quantization axis
respectively, $\phi$ is the intersection angle between the linearly polarized
component of light field and $\hat{e}_{B}$. This formula comes from the
perturbation expansion. Note that the range of values of light ellipticity is
$\xi\in[-1,1]$, $\xi=\pm 1$ denotes left and right circular polarization.
${\alpha}^{\left(0\right)}\left(\omega\right),{\alpha}^{\left(1\right)}\left(\omega\right),{\alpha}^{\left(2\right)}\left(\omega\right)$
are the scalar, vector, and tensor polarizability, respectively. Scalar shift
is spin independent. Vector shift acts like an effective magnetic field to
generate the linear Zeeman splitting (light shift proportional to $m_{F}$),
which depends on the ellipticity of the light and the intersection angle
between the laser beam wave vector and magnetic field quantization axis
$\hat{e}_{B}$. So there are two methods to control the vector shift, rotating
bias magnetic field and changing light polarization. Tensor part is derived
from the linearly polarized light, and acts as an effective dc electric field.

For the first excited state of alkali-metal atoms, the fine structure
interaction induces the spectral lines of the D1 ($5^{2}S_{1/2}\rightarrow
5^{2}P_{1/2}$) and D2 ($5^{2}S_{1/2}\rightarrow 5^{2}P_{3/2}$) lines. The
coefficients of the scalar, vector and tensor shifts of the ground states
$5^{2}S_{1/2}$ of 87Rb atoms in Eq. 2 are given by

$\displaystyle{{{\alpha}^{\left(0\right)}}\left(\omega\right)\approx-\frac{\pi{{c}^{2}}{{\Gamma}_{{{D2}}}}}{2\omega^{3}_{0}}\left(\frac{2}{{{\delta}_{{{D2}}}}}+\frac{1}{{{\delta}_{{{D1}}}}}\right)},$
(3)
$\displaystyle{{\alpha}^{\left(1\right)}}\left(\omega\right)\approx-\frac{\pi{{c}^{2}}{{\Gamma}_{{{D2}}}}}{2\omega^{3}_{0}}\left(\frac{1}{{{\delta}_{{{D2}}}}}-\frac{1}{{{\delta}_{{{D1}}}}}\right){g}_{F}{F},$
$\displaystyle{{\alpha}^{\left(2\right)}}\left(\omega\right)\approx 0,$

where ${{\Gamma}_{{{D}_{2}}}}$ is the decay rate of the excited state for
$D_{2}$ line, $\omega_{0}=\frac{1}{3}\omega_{D1}+\frac{2}{3}\omega_{D2}$ is
the effective frequency,
${{\delta}_{{{D}_{1}}}}={{\omega}_{~{}}}-{{\omega}_{D1}}$,
${{\delta}_{{{D}_{2}}}}={{\omega}_{~{}}}-{{\omega}_{D2}}$ is the frequency
detuning of the laser. Therefore, according to Eq. 3 we only consider the
scalar and vector shift in this work. We employ tune-out wavelength for spin-
dependent optical lattice, in which ac Stark shifts cancel. Two internal spin
states have different tune-out wavelengths when the contributions of both the
scalar and vector shifts are included Wen2021 .

We choose two ground hyperfine Zeeman states of 87Rb
$\left|F=2,m_{F}=0\right\rangle$ of the $F=2$, and $\left|1,1\right\rangle$ of
the $F=1$ hyperfine manifold as the two internal spin states. A bias magnetic
field with 10 Gauss is applied along the $45^{\circ}$ diagonal line of the
square lattice $V_{2}$. We scan the wavelength of the optical lattice beams to
determine the tune-out wavelength precisely, as shown in Fig. S2. The tune-out
wavelength for $\left|1,1\right\rangle$ state is determined at 788.28 nm with
$\sigma_{-}$ circular polarization as shown in Fig. S2(c), which balances the
contribution of the scalar and vector shift. Thus we choose this wavelength
for the lattice $V_{2}$. Note that the tune-out wavelength for
$\left|1,1\right\rangle$ state is sensitive to the intersection angle between
the laser beam wave vector and magnetic field quantization axis, which
requires a careful alignment of the bias magnetic field carefully. The spin
state $\left|2,0\right\rangle$ only experiences the square lattice $V_{2}$
with the red-detuning AC stark shift (which is only from scalar shift) as
shown in Fig. S2(d) and (f), in contrast, the spin state
$\left|1,1\right\rangle$ experiences no shift.

On the other hand, there is only the contribution of the scalar shift for the
spin state $\left|2,0\right\rangle$, the tune-out wavelength for
$\left|2,0\right\rangle$ state is 790.02 nm as shown in Fig. S2(a), which is
well known and studied experimentally arora2011tune ; Wen2021 . We choose this
tune-out wavelength of 790.02 nm with $\sigma_{+}$ circular polarization as
the wavelength of the lattice $V_{1}$. Thus the spin state
$\left|1,1\right\rangle$ experiences the square lattice $V_{1}$ with the red-
detuned AC stark shift. In contrast, the spin state $\left|2,0\right\rangle$
see zero light shift. Note that the tune-out wavelength for
$\left|2,0\right\rangle$ state is insensitive to the intersection angle
between the laser beam wave vector and magnetic field quantization axis. The
spin state $\left|1,1\right\rangle$, however, has a different lattice depth in
two orthogonal directions of the lattice $V_{1}$ respectively, and feels the
lattice $V_{1}$ with the red-detuning AC stark shift (which is only from
vector shift at the wavelength of 790.02 nm) as shown in Fig. S2(b) and (e).

Moir$\acute{\mathrm{e}}$ superlattice can be generated by a small difference
in lattice constant or orientation. Since two different wavelengths are used
for twisted-bilayer lattices in this work, there is a large-period
superlattice with $\Delta\lambda=179$ $\mu$m, much larger than the size of
atomic cloud. Therefore, we can adjust the retroreflected concave mirror to
load atoms into the lower potential well of the long-period superlattice and
neglect the influence on the measurement of moir$\acute{\mathrm{e}}$ pattern.
In the future, we can correct this effect of two different wavelengths by
using a slight angle lattice beam for $V_{2}$ to ensure the same lattice
constant for two lattice potentials.

III. Band structures and flat band

Figure S3: Band structure of the twisted-bilayer optical lattices. The twist
angle of the commensurate optical lattice is $\theta=2\arctan(1/22)$, whose
band structure is regarded as an approximation of the experimental case
$\theta=5.21^{\circ}$. $\bf{a}$, $\bf{b}$ and $\bf{c}$ show the band
structures for the interlayer coupling strength $\Omega_{\mathrm{R}}=0E_{r}$,
0.1$E_{r}$ and 1$E_{r}$ respectively. $\bf{a}$ also gives the band structure
without the interlayer coupling in the form of the superlattice minibands
within the same reduced Brillouin zone. $\bf{d}$, $\bf{e}$ and $\bf{f}$ are
the enlargement of the lowest bands of $\bf{a}$, $\bf{b}$, and $\bf{c}$,
respectively. $\bf{g}$, $\bf{h}$ and $\bf{i}$ are the further enlargement of
the lowest bands of $\bf{d}$, $\bf{e}$ and $\bf{f}$, respectively. Here, the
MW detuning is $\Delta=0$ and $E_{0}$ corresponds to the energy of the lowest
band.

For square lattices, the commensurate angles $\theta$ satisfy
$\tan(\theta/2)=\bar{m}/\bar{n}$, where $\bar{m}$ and $\bar{n}$ are integers.
An equivalent condition is $\cos\theta=\bar{a}/\bar{c}$ and
$\sin\theta=\bar{b}/\bar{c}$, which can be defined by Pythagorean triples
($\bar{a}^{2}+\bar{b}^{2}=\bar{c}^{2}$, where $(\bar{a},\bar{b},\bar{c})\in N$
are positive integers) Wang2020 . The relationship between
($\bar{m}$,$\bar{n}$) and $(\bar{a},\bar{b},\bar{c})$ is
$(\bar{m}+i\bar{n})^{2}=(\bar{a}+i\bar{b})$ when $(\bar{m}+\bar{n})\in odd$
and $(\bar{m}+i\bar{n})^{2}=2(\bar{a}+i\bar{b})$ when $(\bar{m}+\bar{n})\in
even$. For the commensurate optical lattice, the period of its supercell is
given by $\lambda_{sc}=\bar{m}a/\sin(\theta/2)=2\bar{m}\lambda_{mo}$ when
$(\bar{m}+\bar{n})\in odd$ and
$\lambda_{sc}=\bar{m}a/\sin(\theta/2)=\sqrt{2}\bar{m}\lambda_{mo}$ when
$(\bar{m}+\bar{n})\in even$.

Here, we choose the band structure of the commensurate optical lattice with
the commensurate angle $\theta=2\arctan(1/22)$ as an approximation of the
experimental case. If getting a better approximation of band structure for the
experimental case, we can choose the larger supercell to calculate the energy
band structure, whose commensurate angle is closer to the experimental case.
The band structure $E(\bf k)$ of the commensurate optical lattice can be
obtained by solving the stationary Schrödinger equation, $H\Psi=E\Psi$, with
the Bloch function, $\Psi({\bf r})=\exp^{i{\bf k}\cdot{\bf r}}u(\bf r)$. Here
the Hamiltonian $H$ is given as

$H=\left(\begin{array}[]{cc}-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{1}+\frac{\Delta}{2}&\Omega_{\mathrm{R}}\\\
\Omega_{\mathrm{R}}&-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{2}-\frac{\Delta}{2}\\\
\end{array}\right),$ (4)

$u(\bf r)$ is a periodic function with the same periodicity as the coupled
lattice. The spin-dependent square optical lattice with a twist angle $\theta$
can be described by the potentials

$\displaystyle
V_{1}=V_{0}\Big{[}\sin^{2}\big{(}kx\cos\frac{\theta}{2}-ky\sin\frac{\theta}{2}\big{)}$
$\displaystyle+\sin^{2}\big{(}ky\cos\frac{\theta}{2}+kx\sin\frac{\theta}{2}\big{)}\Big{]},$
(5a) $\displaystyle
V_{2}=V_{0}\Big{[}\sin^{2}\big{(}kx\cos\frac{\theta}{2}+ky\sin\frac{\theta}{2}\big{)}$
$\displaystyle+\sin^{2}\big{(}ky\cos\frac{\theta}{2}-kx\sin\frac{\theta}{2}\big{)}\Big{]},$
(5b)

where $k=2\pi/\lambda$ is the wave number of lasers for the lattice and
$V_{0}$ describes the lattice depth. In numerics we first discretize the unit
supercell of area $\sqrt{\bar{c}}a\times\sqrt{\bar{c}}a$ in real space
($\bar{c}$ is the largest value in the Pythagorean triple) into $l\times l$
grids, and then diagonalize the effective Hamiltonian for $u(\bf r)$. As shown
in Fig. S3, the band structure approaches the flat band when increasing the
interlayer coupling.

Since our system allows for flexible control of the interlayer couplings, the
flatband in the lowest energy band can be realized. The Hamiltonian Eq. 4 can
be formally diagonalized as

$H=\left(\begin{array}[]{cc}H_{\text{eff}}^{+}&0\\\ 0&H_{\text{eff}}^{-}\\\
\end{array}\right),$ (6)

where

$H_{\text{eff}}^{\pm}=h_{0}\pm h_{1},$ (7)

with $h_{0}=-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{V_{1}+V_{2}}{2}$, and
$h_{1}=\sqrt{\Omega_{\mathrm{R}}^{2}+\frac{(V_{1}-V_{2}+\Delta)^{2}}{4}}$. In
the large interlayer coupling limit, $\Omega_{\mathrm{R}}\gg V_{0},\Delta$,
the low-energy band structure is encoded in the effective Hamiltonian
$H_{\text{eff}}^{-}$ in the lower-right block, which can be further
approximated as

$H_{\text{eff}}^{-}\approx-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{V_{1}+V_{2}}{2}-\frac{(V_{1}-V_{2}+\Delta)^{2}}{8\Omega_{\mathrm{R}}}-\Omega_{\mathrm{R}},$
(8)

or in a more rough way

$H_{\text{eff}}^{-}\approx-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{V_{1}+V_{2}}{2}-\Omega_{\mathrm{R}},$
(9)

The approximated effective Hamiltonians correspond to some effective lattices
for a single-layer (single-component) system, separately,
$V=\frac{V_{1}+V_{2}}{2}-\frac{(V_{1}-V_{2}+\Delta)^{2}}{8\Omega_{\mathrm{R}}}$
for Eq. 8, and $V=\frac{V_{1}+V_{2}}{2}$ for Eq. 9, with certain global energy
shift. Specifically Eq. 9 implies that the system becomes the single-layer
(single-component) experiencing a twisted optical lattice.

When increasing the interlayer coupling into the strong region, the long-
wavelength moir$\acute{\mathrm{e}}$ potential becomes deeper, so atoms in the
lowest band are isolated at a larger spatial scale (moir$\acute{\mathrm{e}}$
wavelength), which enhances the wavefunction localization and contributes to
the creating of flatband. The single-layer system with a twisted optical
lattice (approximation at the strong interlayer coupling limit) admits a
flatband structure in the lowest band, which has been studied experimentally
in photonic system Wang2020 ; Huang2016 ; Fu2020 . The
moir$\acute{\mathrm{e}}$ flat bands have several advantages. First, the flat
bands quench kinetic energy scales (wavefunction localization), thereby
drastically enhance the role of interactions and amplify the effects of
interactions. Second, the moir$\acute{\mathrm{e}}$ superlattice leads to the
emergence of minibands within a reduced Brillouin zone. The small Brillouin
zone means that low atomic densities are sufficient for full filling or
depletion of the superlattice bands, which is easily controlled in experiment.

IV. Theoretical calculation of the modified superfluid to insulator transition

In the mean-field approximation, the system for the superfluid phase can be
well described by the coupled Gross-Pitaevskii (GP) equations

$\displaystyle i\hbar\frac{\partial\psi_{1}}{\partial
t}=\Big{[}-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{1}{2}m\omega^{2}_{\perp}(x^{2}+y^{2})+V_{1}+$
$\displaystyle\eta g_{11}|\psi_{1}|^{2}+\eta
g_{12}|\psi_{2}|^{2}\Big{]}\psi_{1}+\hbar\Omega_{\mathrm{R}}\psi_{2},$ (10)
$\displaystyle i\hbar\frac{\partial\psi_{2}}{\partial
t}=\Big{[}-\frac{\hbar^{2}}{2m}\nabla^{2}+\frac{1}{2}m\omega^{2}_{\perp}(x^{2}+y^{2})+V_{2}+$
$\displaystyle\eta g_{12}|\psi_{1}|^{2}+\eta
g_{22}|\psi_{2}|^{2}\Big{]}\psi_{2}+\hbar\Omega_{\mathrm{R}}\psi_{1},$ (11)

where the MW detuning is $\Delta=0$ and the wave function is normalized as
$\sum_{i}\int|\psi_{i}|^{2}d\mathbf{r}=N$, with $N$ the total atom number. The
strong confinement along the $z$ axis gives rise to the quasi-2D interaction
strengths represented by a reduction coefficient $\eta$ multiplied by
$g_{ij}=4\pi\hbar^{2}a_{ij}/m$, where $\eta^{-1}=\sqrt{h/m\omega_{z}}$ defines
the characteristic length along the z axis, and $a_{ij}$ is the 3D s-wave
scattering length. In the experiment, the trapping frequency
$\omega_{z}\approx 2\pi\times 1\mathrm{kHz}$, and the scattering length for
the 87Rb atoms is about $a_{ij}\approx 100a_{B}$ with $a_{B}$ the Bohr radius.
This implies that even though the system is thermodynamically 2D, the
collisions still keep their 3D character with $\eta^{-1}\gg a_{ij}$.
Considering the similarity in scattering lengths $a_{11}$, $a_{22}$ and
$a_{12}$ for the 87Rb atoms, in the calculation we focus on the SU(2)
symmetric interaction with $g=g_{11}=g_{22}=g_{12}$. In addition to the
intercomponent atomic interaction, the two components are also coupled by a
microwave pulse, which causes Rabi oscillations with frequency
$\Omega_{\mathrm{R}}$.

By using the imaginary time evolution method, one can solve the GP equations
numerically for the ground states in the harmonic trap. Theoretically the non-
commensurate twist angle $\theta=5.21^{\circ}$ should be a localized single
particle ground state while the commensurate angle
$\theta=2\arctan{\frac{1}{22}}$ gives rise to extended ground states in the
absence of interactions. Experimentally the interatomic interaction is
dominant, and always leads to extended many-body states with the aperiodic and
periodic bilayer lattices becoming almost indistinguishable.

The phase transition from superfluid to Mott-insulator can be well described
by the Bose-Hubbard model in the tight-binding approximation. For simplicity
we consider the interlayer coupling as a quasi-periodically perturbed
potential, which leads to a site-dependent energy offset

$\displaystyle
M_{i}=M_{R}\big{[}\sin^{2}(i_{x}\pi\cos\theta+i_{y}\pi\sin\theta)$
$\displaystyle+\sin^{2}(i_{y}\pi\cos\theta-i_{x}\pi\sin\theta)\big{]},$ (12)

where the subindex $i_{x}$ and $i_{y}$ label the position of the $i$-th site
in the two-dimensional space. The tight-binding Hamiltonian for one layer then
is given by

$H=-t\sum_{\langle
i,j\rangle}b^{\dagger}_{i}b_{j}+\frac{U}{2}\sum_{i}\hat{n}_{i}(\hat{n}_{i}-1)+\sum_{i}(M_{i}-\mu)\hat{n}_{i},$
(13)

where the first term describes the nearest-neighbor tunnelling with
$b^{\dagger}$ and $b$ being the creation and annihilation operators, and the
second term represents the on-site interaction. The hopping amplitude $t$ is
considered to be site-independent for weak interlayer coupling and can be
estimated by
$t=\frac{4}{\sqrt{\pi}}E_{r}\left(V_{0}/E_{r}\right)^{3/4}e^{-2\left(V_{0}/E_{r}\right)^{1/2}}$.
The local repulsion $U$ depends on the depth of the optical lattice, and is
given by $U=\sqrt{8/\pi}ka_{s}E_{r}\left(V_{0}/E_{r}\right)^{3/4}$
Zwerger_2003 (47). The chemical potential $\mu$ controls the average number of
atoms in the moiré lattice.

The mean-field phase diagram (Fig. 5a in the Main text) can be mapped by using
the Gutzwiller method, which expands the local state $|\psi_{i}\rangle$ at
site $i$ in the Fock basis PhysRevB.45.3137 (48, 49). When the interlayer
coupling $M_{i}=0$, the system is reduced to the standard Bose-Hubbard model
PhysRevB.40.546 , which includes two phases, the superfluid phase and the Mott
insulator phase PhysRevLett.81.3108 (50, 45). While the superfluid phase is
identified by the superfluid order parameter $\langle\hat{b}_{i}\rangle\neq 0$
and an arbitrary filling of the atoms on the site, the Mott insulator phase
emerges with an integer number of atoms per site with
$\langle\hat{b}_{i}\rangle=0$. When the interlayer coupling $M_{i}\neq 0$, the
persistent coherence of the moir$\acute{\mathrm{e}}$ and primary lattice
length scale, as well as density distribution in real space can be used to
distinguish the phases, which is determined by the order parameter
$\langle\hat{b}_{i}\rangle$ and the filling of the atoms on the site $n$ as
shown in Fig. S4. The chemical potential $\mu/U=1$ is considered in this
calculation.

Figure S4:  Characteristics of the different phases. $\bf{a}$ Phase diagram,
where SF, SF-II, MI, and I refer to superfluid, superfluid only with short-
range coherence, Mott insulator, and insulator. $\bf{b}$ Table shows the
features of the different phases. $\bf{c}$ Plots of the order parameter
$\langle\hat{b}_{i}\rangle$ and the filling of the atoms on the site $n$ for
the different phases. Parameters ($V/E_{r}$, $\Omega_{\mathrm{R}}/E_{r})$ are
(10,0.6), (15,0.6), (23,0.3) and (23,1.1) for the plots from left to right
respectively. The chemical potential $\mu/U=1$ is considered.

Phase diagram with zero temperature and the homogeneous system is predicted
theoretically as shown in Fig. S4a, in which four phases of superfluid (SF),
superfluid II (SF-II), Mott insulator (MI), and insulator (I) are included.
The SF-II phase is a state with superfluid domains embedded in a gapped
insulate state, which is caused by interlayer coupling. So SF-II phase can be
identified by checking the disappearance of the moir$\acute{\mathrm{e}}$-scale
long-range correlation with vanishing moir$\acute{\mathrm{e}}$ lattice
momentum but the remaining of short-range coherence with residual primary
lattice momentum. At the same time, the SF-II phase supports the
moir$\acute{\mathrm{e}}$ pattern in the real-space. Since MI is an
incompressible insulator for integer filling factor with a gap $U$ for
particle–hole excitations induced by the onsite interaction $U$, the
moir$\acute{\mathrm{e}}$ pattern in the real-space appears only when the
interlayer coupling strength is larger than $U$ to break this gap. Therefore,
the I phase supports the moir$\acute{\mathrm{e}}$ pattern in the real-space
and no spatial coherence at all scales, which approaches the MI phase without
the moir$\acute{\mathrm{e}}$ pattern in the real-space in the limit of weak
interlayer coupling. Here, the phase transition between the SF to the MI phase
should have an intermediate phase SF-II. Obviously, the coherence is lost
almost simultaneously in all length scales at the critical point at very weak
interlayer coupling only for zero temperature as shown in Fig. 5a. In
contrast, at finite temperature, the thermodynamic quantities behave smoothly
near the critical point and the long-range moir$\acute{\mathrm{e}}$ coherence
is lost more before the short-range primary lattice length. Therefore, there
exists an SF-MI critical regime, which seems likely to be a thermodynamic
phase and is not predicted theoretically at zero temperature. The SF-MI
critical regime has the same coherence as SF-II without the
moir$\acute{\mathrm{e}}$-scale long-range correlation but with the short-range
coherence. However, the SF-MI critical regime does not support the
moir$\acute{\mathrm{e}}$ pattern in the real-space which is distinguished from
SF-II. In fact, the insulator phase presents special characteristics due to
the strong interaction and quasi-disorder induced by the larger interlayer
coupling, such as similar to Bose glass insulator from the model of disordered
strongly interacting bosonic system PhysRevB.40.546 ; PhysRevLett.67.2307 ;
PhysRevB.53.2691 .

## References

  * (1) Xiong, D., Wang, P., Fu, Z., Chai, S. $\&$ Zhang, J.  Evaporative cooling of 87Rb atoms into Bose-Einstein condensate in an optical dipole trap.  _Chin. Opt. Lett._ 8, 627–629 (2010).
  * (2) Greiner, M., Mandel, O., Esslinger, T., H$\ddot{a}$nsch, T. W. $\&$ Bloch, I.  Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms.  _Nature_ 415, 39 (2002).
  * (3) Steck, D. A.  Quantum and atom optics  (2007).
  * (4) Zwerger, W.  Mott Hubbard transition of cold atoms in optical lattices.  _J. Opt. B: Quantum Semiclass. Opt._ 5, S9–S16 (2003).
  * (5) Krauth, W., Caffarel, M. $\&$ Bouchaud, J.-P.  Gutzwiller wave function for a model of strongly interacting bosons.  _Phys. Rev. B_ 45, 3137–3140 (1992).
  * (6) Sheshadri, K., Krishnamurthy, H. R., Pandit, R. $\&$ Ramakrishnan, T. V.  Superfluid and Insulating Phases in an Interacting-Boson Model: Mean-Field Theory and the RPA.  _Europhys. Lett._ 22, 257–263 (1993).
  * (7) Jaksch, D., Bruder, C., Cirac, J. I., Gardiner, C. W. $\&$ Zoller, P.  Cold Bosonic Atoms in Optical Lattices.  _Phys. Rev. Lett._ 81, 3108–3111 (1998).


# Surface induced electronic Berry curvature in Berry curvature free bulk
materials

Dennis Wawrzik1 Jorge I. Facio1,2 Jeroen van den Brink1,3

1Institute for Theoretical Solid State Physics, IFW Dresden, Helmholtzstr. 20,
01069 Dresden, Germany
2Centro Atomico Bariloche, Instituto de Nanociencia y Nanotecnologia (CNEA-
CONICET) and Instituto Balseiro. Av. Bustillo 9500, Bariloche (8400),
Argentina
3Institute for Theoretical Physics and Cluster of Excellence ct.qmat, TU
Dresden, 01069 Dresden, Germany

###### Abstract

In recent years it has become clear that electronic Berry curvature (BC) is a
key concept to understand and predict physical properties of crystalline
materials. A wealth of interesting Hall-type responses in charge, spin and
heat transport are caused by the BC associated to electronic bands inside a
solid: anomalous Hall effects in magnetic materials[1, 2], and various
nonlinear Hall[3, 4, 5] and Nernst effects[6, 7] in non-magnetic systems that
lack inversion symmetry. However, for the largest class of known materials
–non-magnetic ones with inversion symmetry– electronic BC is strictly zero[8,
9]. Here we show that precisely for these bulk BC-free materials, a finite BC
can emerge at their surfaces and interfaces. This immediately activates
certain surfaces in producing Hall-type transport responses. We demonstrate
this by first principles calculations of the BC at bismuth, mercury-telluride
(HgTe) and rhodium surfaces of various symmetries, revealing the presence of a
surface Berry curvature dipole and associated quantum nonlinear Hall effects
at a number of these. This opens up a plethora of materials to explore and
harness the physical effects emerging from the electronic Berry curvature
associated exclusively to their boundaries.

Inside a solid the anomalous velocity of Bloch electrons caused by the
electronic Berry curvature (BC) leads to a sideways deflection of currents,
much as an external magnetic field[1, 8, 9]. Such a current deflection caused
by the intrinsic electronic structure of a material and its associated Hall-
type response must be allowed by the symmetries of the system. It is well-
known that an absence of magnetism, and consequently a presence of time-
reversal symmetry in the electronic band structure, leads to zero total BC,
integrated over Brillouin-zone momenta. This consequently forbids any linear
Hall response. But as long as at least lattice inversion symmetry is broken, a
finite BC is associated to the electronic structure, which then can allow,
e.g., a quantum nonlinear Hall effect where the observed transversal voltage
is proportional to the applied current squared. This nonlinear response is
related to the dipole moment of the Brillouin-zone integrated BC, the Berry
curvature dipole (BCD). For a bulk BCD to be present in a system with bulk BC,
the crystal must belong to certain point groups[3].

However, for materials with both time-reversal invariance and inversion, the
BC vanishes identically at all momenta in the 3D Brillouin-zone. So from the
BC point of view such materials seem to become uninteresting as their
anomalous electron velocities now appear to vanish as well. We demonstrate
that while this conclusion is valid for 3D bulk materials, at surfaces (and
interfaces) of these the opposite holds: at a general Miller index surfaces
Bloch-electrons do attain a finite anomalous velocity, also for materials with
bulk inversion and time-reversal symmetry. Indeed we find that as a general
rule of thumb the surface BC only vanishes at unreconstructed low Miller index
surfaces, whereas a finite surface BC is associated to any higher index
surface, including miscut and vicinal ones with surface steps that break
C2-rotation around the surface normal. If there is in addition not more than
one mirror plane containing the normal, a finite surface BCD forms. Very many
elementary surfaces fulfil this symmetry requirement.

To explicitly evaluate the surface Berry curvature (SBC) and the SBC dipole
(SBCD) for real materials we start from symmetry-conserving maximally-
projected Wannier functions from density-functional calculations using the
Full-Potential-Local-Orbital code, FPLO[10, 11], in order to construct bulk
Hamiltonians. We subsequently obtain the surface Green’s function for a semi-
infinite slab of material using an iterative construction (see Methods). After
extracting the surface Green’s function from the bulk Hamiltonian we evaluate
its associated SBC. The numerically efficiency of calculating the Green’s
function for a semi-infinite system sets this method apart from a direct real-
space slab approach, the numerical cost of which strongly increases with slab
thickness and which would in addition require disentangling the BC
contribution of different layers in the finite slab.

We demonstrate the principles and method outlined above by considering three
different surfaces of inversion and time-reversal invariant elemental bismuth,
a strictly BC-free bulk crystal. The effects of (surface) symmetry can be
nicely illustrated on Bi crystals, which have a 3-fold rotoinversion axis (cf.
Fig. 1(a-b)) that cancels all bulk BC. The ($\bar{1}$10) surface[12] of Bi has
a two-dimensional center of inversion (which is equivalent to a C2-rotation
symmetry around the surface normal) and thus no SBC, whereas at the (111) and
(100) surfaces inversion is broken, see Fig. 1(c). Consequently Bi (111) and
(100) have a finite SBC associated to them, as the results of a direct
computation shown in Fig. 2(a-c) confirm. Since the bulk inversion relates top
and bottom surface, opposite surfaces will have opposite SBC. Consequently the
total BC for an inversion symmetric piece of bulk crystal including its
surfaces, vanishes again. We note that for Bi (100) in Fig. 2(c) the
projection of the bulk states gives a large contribution to the BC close to
$\vec{k}_{\parallel}=(0,\pi)$ while the states around $\Gamma$ and
$\vec{k}_{\parallel}=(\pi,\pi)$ are clearly separated from the bulk.

From the symmetry point of view the Bi (100) surface has only a single mirror
line, as depicted in Fig. 1(c). Thus it allows for a non-zero SBCD that in
principle is accessible experimentally: an electric field parallel to the SBCD
will result in a nonlinear current perpendicular to it, flowing along the
mirror line, see Fig. 2(c). The SBCD on other crystal facets with a single
mirror line will also give rise to a current in direction of the 3-fold
rotation axis, as depicted in Fig. 1(b). To quantitatively determine the BCD
associated to the Bi (100) surface states we evaluate the SBC dipole
${\vec{D}}^{S}$ for a k-mesh with $512^{2}$ k-points and an energy window
$\Delta\mu=0.04\,\text{eV}$, see Fig. 1(d). This yields
$D^{S}_{x}=0.03\,\text{\AA}$ for the state at $\vec{k}_{\parallel}=(\pi,\pi)$
and $D^{S}_{x}=-0.22\,\text{\AA}$ at $\Gamma$. Of course deeper layers will
also contribute to the total BC dipole associated to this surface, but it is
interesting to note that the top layer by itself already has a BCD that is of
the same order as the BCD of monolayer WTe2 [13].

The binary compound mercury-telluride, HgTe, has gained much interest as the
basis material for topologically insulating phases[14, 15]. Its zincblende
structure breaks inversion allowing for a finite bulk BC but has too many
symmetries for a bulk BC dipole. Similar to Bi we find a finite SBC on the
(111) surface but no SBCD due to the 3-fold rotation and the three mirror
lines, see Fig. 2(d). In contrast, HgTe possesses two-fold rotation axes which
is naturally preserved by surfaces perpendicular to it. Since a two-fold
rotation is identical to inversion in 2D, the HgTe (100) surface has zero SBC
(Fig. 2(f)). On the other hand, the ($\bar{1}$10) surface only has a single
mirror line yielding a clear SBCD (Fig. 2(e)). This implies that the natural
(110) cleaving plane[16] of HgTe, has a SBCD associated to it.

A feature setting three-dimensional bulk HgTe apart is it becoming a
topological insulator when tensile strain is applied and a gap opens up.
Straining HgTe along the xy-plane, the inversion on (100) and the mirror on
($\bar{1}$10) remain and a SBC is associated to the topological surface state
on ($\bar{1}$10). Strain breaks the 3-fold rotation and several mirror planes
and consequently the strained HgTe (111) surface, and in particular its
topological surface state, now also acquires a finite SBCD, as shown in Fig.
2(g-i).

Similar considerations as the ones above for Bi and HgTe apply to other
interfaces and materials. For example, a plethora of noble metal surfaces are
experimentally available with electronic properties that that are important
for catalysis [17]. From our analysis we conclude that (111) surfaces of hcp
and fcc noble metals (e.g., Au, Pt, Rh, Ir) have a finite SBC, as is
illustrated by a direct calculation for elemental Rh in Fig. 3. For such (111)
surfaces the SBCD vanishes but it can become finite when reducing the
symmetry, e.g. by applying uniaxial strains as shown above for HgTe (111).
Since these noble metals are used as substrates of other 2D materials[18],
e.g. graphene[19], a proper use of the latter can be another interesting path
to reduce the symmetry and induce a finite SBCD. Alternatively, higher Miller
index surfaces can be considered. In Fig. 3(d) we show the BC associated to Rh
(311), which is a stepped surface with a symmetry reduced to a single mirror
line, resulting in a SBCD. Similarly a SBCD is also associated to, e.g.,
chemically relevant Rh (210) [20] and Rh (310) [21] surfaces. The rather
complex Fermi surface of bulk Rh produces an intricate Rh (311) SBC
distribution in Fig. 3, certain detail of which depends on longer range-matrix
elements (see Methods), but it has a dstinct SBCD associated to it due to the
presence of just single mirror line.

Having established the presence of SBCs and SBCDs at the boundaries of a very
large class of materials, including the semimetal Bi, topological HgTe and the
transition metal Rh, it should be very interesting to corroborate these
theoretical findings by experimentally determining the momentum-resolved
surface Berry curvature by dichroic angle-resolved photoelectron
spectroscopy[22], and to observe the emerging surface Berry curvature dipole
that we predict here by the presence of a nonlinear Hall effect at relevant
surfaces and interfaces.

##  1 Methods

Berry Curvature from Surface Green’s function — Associated to the Hamiltonian
matrix $H_{\bm{k}}$ of a 3D system is the imaginary-frequency Green’s function
$G(i\omega,{\bm{k}})=(i\omega-H_{\bm{k}})^{-1}$ and a Berry curvature of all
occupied bands given by [23, 24]

$\Omega_{c}({\bm{k}})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\text{tr}\left(g_{c}(i\omega,{\bm{k}})\right)$
(1)

with

$g_{c}(i\omega,{\bm{k}})=\frac{\epsilon_{\mu\nu\rho
c}}{3!}~{}G\left(\partial_{\mu}G^{-1}\right)G\left(\partial_{\nu}G^{-1}\right)G\left(\partial_{\rho}G^{-1}\right)$
(2)

where $c\in\\{x,y,z\\}$ and $\mu,\nu,\rho\in\\{\omega,x,y,z\\}$. Whereas
evaluating this Green’s function expression is not a very efficient method to
calculate the BC for a 3D band structure, we rely on it to compute the SBC
after extracting the surface Green’s function from the bulk Hamiltonian. The
iterative algorithm in Ref. [25] is usually used to compute surface spectral
functions of semi-infinite systems. The method can efficiently calculate the
Green’s functions $G^{S}(i\omega,{\bm{k}}_{\parallel})$ of the top, bottom,
and central layers of a system of size $L$ for very large $L$ by doubling the
slab size in every iteration. The corresponding tight-binding Hamiltonian
$H(\vec{k}_{\parallel})$ must be of the form

$\displaystyle H=\begin{pmatrix}H_{0}&H_{t}&~{}&~{}&~{}\\\
H_{t}^{\dagger}&H_{0}&H_{t}&~{}&~{}\\\ ~{}&H_{t}^{\dagger}&H_{0}&\ddots&~{}\\\
~{}&~{}&\ddots&\ddots&H_{t}\\\
~{}&~{}&~{}&H_{t}^{\dagger}&H_{0}\end{pmatrix},$ (3)

i.e. it only allows for hoppings between nearest neighbors (perpendicular to
the surface of interest). For longer hoppings, the unit cell has to be
enlarged correspondingly. The algorithm for the surface Green’s function has
initial values

$\displaystyle\epsilon^{t}_{0}$
$\displaystyle=\epsilon^{b}_{0}=\epsilon^{c}_{0}=H_{0}$ (4)
$\displaystyle\alpha_{0}$ $\displaystyle=H_{t}$ (5) $\displaystyle\beta_{0}$
$\displaystyle=H_{t}^{\dagger}$ (6)

and the iterative step is:

$\displaystyle G$ $\displaystyle=\left(i\omega-\epsilon^{c}_{i}\right)^{-1}$
(7) $\displaystyle\epsilon^{t}_{i+1}$
$\displaystyle=\epsilon^{t}_{i}+\alpha_{i}G\beta_{i}$ (8)
$\displaystyle\epsilon^{b}_{i+1}$
$\displaystyle=\epsilon^{b}_{i}+\beta_{i}G\alpha_{i}$ (9)
$\displaystyle\epsilon^{c}_{i+1}$
$\displaystyle=\epsilon^{c}_{i}+\alpha_{i}G\beta_{i}+\beta_{i}G\alpha_{i}$
(10) $\displaystyle\alpha_{i+1}$ $\displaystyle=\alpha_{i}G\alpha_{i}$ (11)
$\displaystyle\beta_{i+1}$ $\displaystyle=\beta_{i}G\beta_{i}$ (12)

The iteration can be stopped if $\alpha_{i}$ and $\beta_{i}$ are sufficiently
small. The Green’s function of the top ($G^{t}$), bottom ($G^{b}$), and center
($G^{c}$) layer is given by:

$\displaystyle G^{t}$
$\displaystyle=\left(i\omega-\epsilon^{t}_{i}\right)^{-1}$ (13) $\displaystyle
G^{b}$ $\displaystyle=\left(i\omega-\epsilon^{b}_{i}\right)^{-1}$ (14)
$\displaystyle G^{c}$
$\displaystyle=\left(i\omega-\epsilon^{c}_{i}\right)^{-1}$ (15)

The numerically efficiency in reaching very large $L$ sets the Green’s
function method apart from a direct real-space slab approach, the numerical
cost of which strongly increases with $L$ and which would in addition require
disentangling the BC contribution of different layers in the finite slab. On
the other hand a real-space slab calculation might in principle accommodate
surface lattice reconstructions, even if slab sizes would be moderate, which
our present Green’s function approach does not. We use the iterative procedure
to access the Green’s function of a surface layer of thickness $N$, where $N$
is determined by the choice of unit cell size at the first iteration of the
method. The results presented below are typically for the choice $N=1$
(corresponding to the top primitive unit cell) unless stated otherwise.
$G^{S}$ has a finite momentum dependent self-energy and only depends on
momenta ${\bm{k}}_{\parallel}=(k_{x},k_{y})$ parallel to the surface. Based on
it, we define and calculate the associated surface Berry curvature
$\Omega^{S}_{z}({\bm{k}}_{\parallel})$ by a frequency integral over
$g^{S}_{z}(i\omega,{\bm{k}}_{\parallel})$ analogous to Eqs. (1) and (2).

In order to obtain the relevant quantity for transport measurements, the BCD,
we use the identity $\frac{\partial f}{\partial k_{j}}=\frac{\partial
f}{\partial\mu}\frac{\partial G^{-1}}{\partial k_{j}}$ where
$f(\mu-H({\bm{k}}))$ is the Fermi function, and we rewrite the BCD in terms of
(surface) Green’s functions as follows:

$\displaystyle D^{S}_{a}$
$\displaystyle=-\int\frac{\text{d}^{2}{\bm{k}}_{\parallel}}{(2\pi)^{2}}~{}\partial_{a}f(\mu,\vec{k}_{\parallel})~{}\Omega^{S}_{z}({\bm{k}_{\parallel}})$
$\displaystyle=-\int\frac{\text{d}^{2}{\bm{k}}_{\parallel}}{(2\pi)^{2}}~{}d^{S}_{a}(\mu,{\bm{k}}_{\parallel})$
(16) $\displaystyle d^{S}_{a}({\bm{k}}_{\parallel})$
$\displaystyle=\int_{-\infty}^{\infty}\frac{\text{d}\omega}{2\pi}~{}\text{tr}\Big{[}\partial_{a}\big{(}G^{S}(\mu,{\bm{k}}_{\parallel})\big{)}^{-1}~{}\partial_{\mu}g^{S}_{z}(i\omega+\mu,{\bm{k}}_{\parallel})\Big{]}$

where $D^{S}_{a}$ is the SBC dipole, $d^{S}_{a}$ its density in momentum space
and we omit their dependence on $\mu$ for brevity. Since only states at the
Fermi surface contribute to the BCD, it is convenient for a numerical
integration over the Brillouin zone to broaden the derivative of
$g^{S}_{z}(i\omega+\mu,\vec{k})$ by approximating it with a finite symmetric
difference quotient with width $\Delta\mu$. In order to reduce numerical
errors in the $\omega$ integral, we determine the poles of $G^{S}$ and shift
the contour away from the poles if necessary and only calculate the velocity
of the bands at the Fermi energy by diagonalizing $G^{S}$. The integration
might require high precision for $G^{S}$ at the projection of bulk states
because $G^{S}$ diverges at $i\omega=0$ along a whole line. In this case, we
checked that leaving a small gap to avoid the line-divergence in the
integration contour at the real axis still provides a good approximation.

Density functional calculations for Bi, Rh and HgTe — Using the Full-
Potential-Local-Orbital code, FPLO [10], we constructed symmetry-conserving
maximally-projected Wannier functions which render an effective Hamiltonian
$H_{\bm{k}}$ of our 3D system depending on momenta
${\bm{k}}=(k_{x},k_{y},k_{z})$. This scheme can be regarded as a zeroth-order
approximation of the maximally localized approach in projecting a subset of
wave functions onto a set of suitably chosen local trial functions with
subsequent orthonormalization [11]. The particular nature of the local
orbitals in FPLO make them an efficient set of projectors, since they are
constructed to be a “chemical” basis, leading to highly localized Wannier
functions, which obey the space group symmetry of the crystal by construction.

Bismuth crystallizes in the space group R$\bar{3}$m and here we use lattice
parameters $a=4.609$ Å and $c=11.975$ Å. HgTe has a cubic zincblende
structure, space group F$\bar{4}$3m. In this work we consider both such
structure as well as one obtained by applying $2\%$ tensile strain that causes
opening of a gap and stabilizes a topological insulating phase [15]. For this,
we first relax the volume by minimization of the total energy obtaining the
equilibrium lattice parameter $a_{0}=6.699$ Å. We then enforce the above-
mentioned deformation and relax by similar means the lattice parameter $c$. To
exploit the remaining symmetries, in this case we work with the space group
I$\bar{4}m2$ with lattice parameters $a=b=4.832$ Å and $c=6.532$ Å. We work in
the generalized gradient approximation and include the spin-orbit coupling in
the four-component formalism. We use a linear tetrahedron method for Brillouin
zone integrations using a $k$-mesh with $12^{3}$ subdivisions in the case of
Bi and $16^{3}$ in the case of HgTe. For the Wannier models, for Bi we include
6p orbitals for every Bi atom and, for HgTe, 6s orbitals of Hg and 5p of Te.
For Rh the calculation is based on a Wannier Hamiltonian with 5s and 4d
orbitals and to simplify the calculation of the surface Green’s function, we
kept matrix elements coupling atoms closer than 12Å.

##  2 Acknowledgments

We thank Inti Sodemann and Jhih-Shih You for fruitful discussions, Ulrike
Nitzsche for technical assistance and the DFG for support through the
Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum
Matter, ct.qmat (EXC 2147, project-id 39085490) and through SFB 1143 (project-
id 247310070) project A5. J.I.F. would like to thank the support from the
Alexander von Humboldt Foundation during the part of his contribution to this
work done in Germany and ANPCyT grants PICT 2018/01509 and PICT 2019/00371.

Figure 1: Berry curvature at elemental bismuth surfaces.  (a) Lattice
structure and symmetries of bismuth. Mirror planes are marked in orange and
the 3-fold rotoinversion is denoted by the red arrow. (b) Surface Berry
curvature (SBC) of different facets of a Bi crystal. Blue facets have a SBC
dipole perpendicular to the mirror line while there is no SBC on black facets.
The grey surface possesses SBC but no SBC dipole. (c) Symmetries of certain Bi
surfaces. (d) SBC dipole density $\vec{d}^{S}$ of the (100) surface where the
mirror symmetry causes vanishing $\vec{D}^{S}_{y}$ and finite
$\vec{D}^{S}_{x}$. An electric field parallel to the SBC dipole,
$\vec{E}\parallel\vec{D}^{S}$, results in a nonlinear perpendicular in-plane
current $\vec{j}\perp\vec{D}^{S}$ along the mirror line.

.

Figure 2: Electronic surface Berry curvature. From left to right, data
corresponds to (111), ($\bar{1}$10) and (100) surfaces for (a)-(c) Bismuth,
(d)-(f) HgTe, and (g)-(i) strained HgTe. The Fermi surface is marked in black
and dashed green lines denote mirror symmetries. If a surface has a SBCD its
direction is marked by the black arrow.

.

Figure 3: Surface structure and Berry curvature of elemental fcc rhodium.
(a)-(b) (111) and (311) surface geometry and symmetries (c)-(d) surface Berry
curvature (SBC) of (111) and (311) top layers. Whereas bulk Rh is Berry
curvature free, Rh (111) surface supports a SBC and Rh (113) even a SBC
dipole. Legends as in Figs. 1,2.

## References

  * [1] Karplus, R. & Luttinger, J. M.  Hall effect in ferromagnetics.  _Phys. Rev._ 95, 1154–1160 (1954).  URL https://link.aps.org/doi/10.1103/PhysRev.95.1154.
  * [2] Chang, C.-Z. _et al._ Experimental Observation of the Quantum Anomalous Hall Effect in a Magnetic Topological Insulator.  _Science_ 340, 167–170 (2013).  URL https://www.science.org/doi/10.1126/science.1234414.
  * [3] Sodemann, I. & Fu, L.  Quantum Nonlinear Hall Effect Induced by Berry Curvature Dipole in Time-Reversal Invariant Materials.  _Phys. Rev. Lett._ 115, 216806 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevLett.115.216806.
  * [4] Moore, J. E. & Orenstein, J.  Confinement-Induced Berry Phase and Helicity-Dependent Photocurrents.  _Phys. Rev. Lett._ 105, 026805 (2010).  URL https://link.aps.org/doi/10.1103/PhysRevLett.105.026805.
  * [5] Ma, Q. _et al._ Observation of the nonlinear Hall effect under time-reversal-symmetric conditions.  _Nature_ 565, 337–342 (2018).  URL https://www.nature.com/articles/s41586-018-0807-6.
  * [6] Xiao, D., Yao, Y., Fang, Z. & Niu, Q.  Berry-phase effect in anomalous thermoelectric transport.  _Phys. Rev. Lett._ 97, 026603 (2006).  URL https://link.aps.org/doi/10.1103/PhysRevLett.97.026603.
  * [7] Yu, X.-Q., Zhu, Z.-G., You, J.-S., Low, T. & Su, G.  Topological nonlinear anomalous Nernst effect in strained transition metal dichalcogenides.  _Phys. Rev. B_ 99, 201410(R) (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.99.201410.
  * [8] Xiao, D., Chang, M.-C. & Niu, Q.  Berry phase effects on electronic properties.  _Rev. Mod. Phys._ 82, 1959–2007 (2010).  URL https://link.aps.org/doi/10.1103/RevModPhys.82.1959.
  * [9] Vanderbilt, D.  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018).
  * [10] Koepernik, K. & Eschrig, H.  Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme.  _Phys. Rev. B_ 59, 1743–1757 (1999).  URL https://link.aps.org/doi/10.1103/PhysRevB.59.1743.
  * [11] Koepernik, K., Janson, O., Sun, Y. & van den Brink, J.  Symmetry Conserving Maximally Projected Wannier Functions.  _arXiv preprint arXiv:2111.09652_ (2021).  URL https://arxiv.org/abs/2111.09652.
  * [12] For a better comparison with HgTe, the Miller indices refer to a simple cubic system, even though Bi slightly deviates from an fcc lattice. Surfaces (111), ($\bar{1}$10), and (100) correspond to Bravais-Miller indices (0001), (11$\bar{2}$0), and (01$\bar{1}$2).
  * [13] You, J.-S., Fang, S., Xu, S.-Y., Kaxiras, E. & Low, T.  Berry curvature dipole current in the transition metal dichalcogenides family.  _Phys. Rev. B_ 98, 121109(R) (2018).  URL https://link.aps.org/doi/10.1103/PhysRevB.98.121109.
  * [14] König, M. _et al._ Quantum spin hall insulator state in HgTe quantum wells.  _Science_ 318, 766 –770 (2007).  URL http://www.sciencemag.org/content/318/5851/766.abstract.
  * [15] Brüne, C. _et al._ Quantum Hall Effect from the Topological Surface States of Strained Bulk HgTe.  _Phys. Rev. Lett._ 106, 126803 (2011).  URL https://link.aps.org/doi/10.1103/PhysRevLett.106.126803.
  * [16] Liu, C. _et al._ Tunable spin helical dirac quasiparticles on the surface of three-dimensional hgte.  _Phys. Rev. B_ 92, 115436 (2015).  URL https://link.aps.org/doi/10.1103/PhysRevB.92.115436.
  * [17] Li, C. & Baek, J.-B.  Recent Advances in Noble Metal (Pt, Ru, and Ir)-Based Electrocatalysts for Efficient Hydrogen Evolution Reaction.  _ACS Omega_ 5, 31–40 (2020).  URL https://pubs.acs.org/doi/10.1021/acsomega.9b03550.
  * [18] Gao, Y., Liu, Y. & Liu, Z.  Controllable growth of two-dimensional materials on noble metal substrates.  _Iscience_ 24, 103432 (2021).  URL https://www.sciencedirect.com/science/article/pii/S2589004221014036.
  * [19] Varykhalov, A. _et al._ Ir(111) Surface State with Giant Rashba Splitting Persists under Graphene in Air.  _Phys. Rev. Lett._ 108, 066804 (2012).  URL https://link.aps.org/doi/10.1103/PhysRevLett.108.066804.
  * [20] Rebholz, M., Prins, R. & Kruse, N.  Subsurface oxygen and its influence on the chemisorption behaviour of CO on Rh (210).  _Surface science_ 269, 293–299 (1992).  URL https://www.sciencedirect.com/science/article/abs/pii/003960289291263B.
  * [21] Ren, D. M. & Liu, W.  A study of chemisorption behavior of carbon monoxide on rhodium surfaces.  _Surf. Sci._ 232, 316–322 (1990).  URL https://www.sciencedirect.com/science/article/pii/003960289090124Q.
  * [22] Schüler, M. _et al._ Local Berry curvature signatures in dichroic angle-resolved photoelectron spectroscopy from two-dimensional materials.  _Sci. Adv._ 6, eaay2730 (2020).  1905.09404.
  * [23] Ishikawa, K. & Matsuyama, T.  A microscopic theory of the quantum Hall effect.  _Nucl. Phys. B_ 280, 523–548 (1987).  URL https://www.sciencedirect.com/science/article/pii/055032138790160X.
  * [24] Witczak-Krempa, W., Knap, M. & Abanin, D.  Interacting Weyl Semimetals: Characterization via the Topological Hamiltonian and its Breakdown.  _Phys. Rev. Lett._ 113, 136402 (2014).  URL https://link.aps.org/doi/10.1103/PhysRevLett.113.136402.
  * [25] Sancho, M. P. L., Sancho, J. M. L. & Rubio, J.  Highly convergent schemes for the calculation of bulk and surface Green functions.  _J. Phys. F_ 15, 851–858 (1985).  URL https://doi.org/10.1088/0305-4608/15/4/009.


The interstellar absorption arise in the Galactic foreground and in the SMC.
Figure 16: Observed (blue) and synthetic (red) He ii lines that are sensitive
to the wind parameters. The sharp interstellar absorptions arise in the
Galactic foreground and in the SMC, and are also modeled with their respective
RV shifts. (a) He ii $\lambda 1640$ line in the COS spectrum (ID 2 in Table
1). (b) He ii $\lambda 4868$ line in the X-SHOOTER spectrum (ID 5 in Table 1).

###  B.4 Synthetic spectrum calculated with the stellar parameters from the
MESA model

Figure 17: Spectral model calculated with the stellar parameters of the most
favored binary evolutionary model. The plots and line styles are the same as
in Figs. 9 and 10. In this plot we use a distance of
$d=$49\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$ to compensate the used lower
luminosities of the binary components. This plot demonstrates that the
predicted surface abundances of the binary evolutionary models are unable to
reproduce the observed spectrum. The nitrogen abundance in this model is too
strong as it overestimates all optical N iv lines and also predicts unseen N v
absorption. The reduced hydrogen abundance leads to an increased helium
abundance. This enrichment would produce too deep He ii lines in the spectrum
of the primary.