# A momentum dependent optical lattice induced by artificial gauge potential

Zekai Chen<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Rochester, Rochester, New York 14627, USA  Center for Coherence
and Quantum Optics, University of Rochester, Rochester, New York 14627, USA
Hepeng Yao  Department of Quantum Matter Physics, University of Geneva, 24
Quai Ernest-Ansermet, CH-1211 Geneva, Switzerland  Elisha Haber  Nicholas P.
Bigelow<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Rochester, Rochester, New York 14627, USA  Center for Coherence
and Quantum Optics, University of Rochester, Rochester, New York 14627, USA

(August 27, 2024)

###### Abstract

We propose an experimentally feasible method to generate a one-dimensional
optical lattice potential in an ultracold Bose gas system that depends on the
transverse momentum of the atoms. The optical lattice is induced by the
artificial gauge potential generated by a periodically driven multi-laser
Raman process, which depends on the transverse momentum of the atoms. We study
the many-body Bose-Hubbard model in an effective 1D case and show that the
superfluid–Mott-insulator transition can be controlled via tuning the
transverse momentum of the atomic gas. We examined our prediction via a
strong-coupling expansion to an effective 1D Bose-Hubbard model and a quantum
Monte Carlo calculation, and discuss possible applications of our system.

### Introduction.–

In the past two decades, artificial gauge potentials have become a powerful
tool for Hamiltonian engineering in cold atom systems. They enable physicists
to use the cold atom ensemble for quantum simulation of complicated condensed
matter systems, and lead to interesting physics such as creating novel
topological defects Sugawa  _et al._ (2018), the spin Hall effect (SHE) Zhu
_et al._ (2006); Liu  _et al._ (2007); LeBlanc  _et al._ (2012); Beeler  _et
al._ (2013) and spin-orbit coupling (SOC) Dalibard  _et al._ (2011); Goldman
_et al._ (2014); Zhai (2015); Zhang and Liu (2018). The SHE in a cold atom
ensemble—a phenomenon that arises from a spin-dependent Lorentz-like force
acting on moving particles from a transverse direction—has drawn significant
attention, as it is closely related to quantum Hall physics and spintronics.
To date, much of the work studying the SHE in cold atom systems is focused on
weakly-interacting single-particle physics such as single-particle spin-
dependent trajectories and the SHE induced spin current. On the other hand,
SOC—a type of interaction that couples the spin of a particle with its
external degree of freedom—is also of extensive interest. It leads to the
novel energy dispersion Huang  _et al._ (2016); Hamner  _et al._ (2015) and
the topological order Lin  _et al._ (2011); Liu  _et al._ (2013); Wu  _et al._
(2016); et al. (2021) in both continuous gas and optical lattice systems.

Meanwhile, optical lattices with ultracold atoms are well developed platforms
for studying quantum phenomena, especially many-body effects of condensed
matter systems Bloch  _et al._ (2008). A familiar aspect of this many-body
physics is the superfluid–Mott-insulator (SF-MI) transition, which has been
well studied by many groups both theoretically Jaksch  _et al._ (1998); Boéris
_et al._ (2016) and experimentally Greiner  _et al._ (2002); Haller  _et al._
(2010); Soltan-Panahi  _et al._ (2011); Boéris  _et al._ (2016). For optical
lattice systems, a natural question arises: what many-body effects in an
optical lattice system involve an artificial gauge potential? There are a few
studies about the many-body effects in a lattice system with a gauge potential
Cai  _et al._ (2012); Cole  _et al._ (2012), but generally the many-body
effect caused by either SOC or the SHE in an optical lattice system demand
more exploration. In particular, a system where the lattice potential is
determined by the motional state of the atoms is of special interest.

In this letter, we propose to realize a novel transverse-momentum-dependent
optical lattice (TMDOL) where the longitudinal lattice potential depends on
the motional state of the atoms in the transverse dimension. The momentum
dependence derives from the same mechanism as the SHE and SOC, which makes it
possible to tune the many-body phase of the system via the transverse momentum
of the ultracold gas. We create such an artificial gauge potential with
translational symmetry in an ultracold pseudospin-1/2 gaseous system via a
periodically driven Raman process. To simplify the problem, we project the
Hilbert space to a single spin component and explore the many-body physics of
the TMDOL with the help of a strong transverse confinement. We construct the
effective 1D Bose-Hubbard model, and investigate the SF-MI phase diagram using
a strong coupling expansion. We then examine our Bose-Hubbard model result
with a quantum Monte Carlo (QMC) calculation. Both calculations show that the
SF-MI phase transition can be induced by changing the average transverse
momentum in such a TMDOL, which arises as a novel many-body effect caused by
the artificial gauge potential. In addition, we propose a feasible
experimental implementation of the TMDOL and discuss the possible
applications.

### Theoretical formalism.–

We first describe the theoretical construction of the optical lattice for a
pseudospin-1/2 system with a periodically driven Hamiltonian. We start by
considering the Hamiltonian

$H=\frac{\vec{p}^{2}}{2M}+V_{ext}(\vec{r})+H_{d},$ (1)

where $M$ is the mass of the particle, $V_{ext}(\vec{r})$ is the spin-
independent external trapping potential and $H_{d}$ is a periodically driven
Hamiltonian
$H_{d}(\vec{r},t)=\frac{\hbar}{2}\Omega_{0}\vec{n}\cdot\vec{\sigma}\cos{\omega
t}$, with $\vec{n}=(\cos{k_{L}x}\cos{k_{L}z},\cos{k_{L}x}\sin{k_{L}z},0)^{T}$
the direction vector, $k_{L}$ the wave vector, $\Omega_{0}$ the amplitude of
the laser field,
$\vec{\sigma}=(\hat{\sigma}_{1},\hat{\sigma}_{2},\hat{\sigma}_{3})^{T}$ the
vector of Pauli matrices, and $\omega$ the frequency of the periodic driving.
According to Floquet theory, we can apply a gauge transformation to the
Hamiltonian described by the micromotion operator $U=e^{-i\sin{\omega
t}\vec{\Omega}(r)\cdot\vec{\sigma}/2\omega}$ Račkauskas  _et al._ (2019);
Novičenko and Juzeliūnas (2019); Chen  _et al._ (2020), where
$\vec{\Omega}=\Omega_{0}\vec{n}$ and the Hamiltonian in the Floquet basis can
be written as

$H_{F}(r,t)=\frac{\left[\vec{p}-\vec{A}(r,t)\right]^{2}}{2M}+V_{ext}(\vec{r}),$
(2)

where $\vec{A}(r,t)=i\hbar U^{\dagger}\nabla U$ is the non-Abelian gauge
potential, with the specific form of the component given in the supplemental
material sup . We can ignore all terms higher than the zeroth-order in the
Fourier series if the adiabatic condition,
$|\bra{\vec{k}}H_{F}^{(n)}\ket{\vec{k}^{{}^{\prime}}}|\ll\hbar\omega(n\neq
0)$, is satisfied Račkauskas  _et al._ (2019); Novičenko and Juzeliūnas
(2019); sup , where
$H_{F}^{(n)}=\frac{1}{\mathcal{T}}\int_{0}^{\mathcal{T}}H_{F}(t)e^{-in\omega
t}dt$ is the $n$-th Fourier component of the Hamiltonian in the Floquet basis,
$\mathcal{T}=2\pi/\omega$, and $\ket{\vec{k}}$ is the momentum state with
momentum $\hbar\vec{k}$.

The gauge field exerts a spatially periodic ‘Lorentz’ force on atoms along the
$x$-axis due to the motion of atoms along the $z$-axis, which makes the
lattice potential $p_{z}$-dependent. The non-zero component of the zeroth
Fourier order of the gauge potential takes the form $A_{z}^{(0)}=\hbar
k_{L}[J_{0}(a)-1]\sigma_{3}/2$, where $J_{0}(a)$ is the zeroth-order Bessel
function of the first kind and $a=a_{0}|\cos{k_{L}x}|$ with
$a_{0}=\Omega_{0}/\omega=4$ throughout this paper. Also we get
$(\vec{A}^{2})^{(0)}=\hbar^{2}k_{L}^{2}\\{\frac{1}{16}a_{0}^{2}(1-\cos{2k_{L}x})+\frac{1}{2}[1-J_{0}(a)]\\}$.
Keeping only the zeroth Fourier component of the gauge potential in eqn.(2)
and projecting to spin up subspace, we get the effective Hamiltonian in the
Floquet basis as

$H_{L}=\frac{\vec{p}^{2}}{2M}+V(x,p_{z})+V_{ext}(\vec{r}),$ (3)

where $E_{r}=\hbar^{2}k_{L}^{2}/2M$ is the recoil energy.
$V(x,p_{z})=\\{(p_{z}+1/2)\left[1-J_{0}(a)\right]+a_{0}^{2}(1-\cos{2k_{L}x})/16\\}E_{r}$
is an optical lattice potential with periodicity $a_{L}=\pi/k_{L}$. In this
work, we only consider the case with $a_{0}=4$. Eqn.(3) shows the
$p_{z}$-dependence of the lattice potential, which results in
$p_{z}$-dependent ground states sup .

### Many-body Hamiltonian.–

We desire to study the many-body effects in our system with a non-negligible
interaction strength. Therefore, we introduce a deep two-dimensional optical
lattice confinement along the $y$ and $z$-axes so that the tunneling along $y$
and $z$ is negligible and the system can be regarded as an array of effective
1D Bose gases. With a tight 2D transverse confinement along the $y$ and
$z$-axes, the ultracold gas will be trapped in an array of effective 1D tubes.
Throughout this paper, we consider the case where the transverse confinement
is strong enough that the tunneling between nearest sites in the $y$ and $z$
is negligible Stöferle  _et al._ (2004); Anisimovas  _et al._ (2016) (e.g. a
lattice depth $V_{\perp}=30E_{r}$). Furthermore, we consider a homogeneous
lattice along the $x$-axis for simplicity. The transverse trapping potential
for each 1D tube can be well approximated by a harmonic trap with the trapping
frequency $\omega_{\perp}=2\sqrt{V_{\perp}E_{r,\perp}}/\hbar\gg\mu$, where
$\mu$ is the chemical potential, and $V_{\perp}$ and
$E_{r,\perp}=\hbar^{2}k_{\perp}/2M$ are the lattice depth and the recoil
energy of the 2D lattice. One can separate the dimensions and the ground state
wavefunctions along the $y$ and $z$-axes can be well approximated by the
harmonic ground state wave packets
$\varphi(y)=(\sqrt{\pi}a_{\perp})^{-\frac{1}{2}}exp\\{-\frac{y^{2}}{2a_{\perp}^{2}}\\}$
and
$\varphi(z)=(\sqrt{\pi}a_{\perp})^{-\frac{1}{2}}exp\\{-(\frac{z^{2}}{2a_{\perp}^{2}}+i\bar{k}_{z}z)\\}$,
where $a_{\perp}=\sqrt{\hbar/M\omega_{\perp}}$ is the characteristic length of
the tight harmonic trap and $\bar{k}_{z}$ corresponds to the average momentum
along $z$-axis, $\langle p_{z}\rangle=\hbar\bar{k}_{z}$.

Due to the uncertainty principle, the tight transverse confinement gives rise
to a non-negligible spreading in momentum space of the atoms along $y$ and
$z$-axes. By integrating along the $y$ and $z$-axes, one gets the single-
particle effective 1D Hamiltonian Barbiero and Salasnich (2014)

$\displaystyle H_{1D}(x,\bar{k}_{z})$ $\displaystyle=$ $\displaystyle\int\int
dydz\varphi^{\ast}(y)\varphi^{\ast}(z)H_{L}\varphi(z)\varphi(y)$ (4)
$\displaystyle=$
$\displaystyle\frac{p_{x}^{2}}{2M}+\sum_{k}|\tilde{\varphi}(k-\bar{k}_{z})|^{2}V(x,p_{z}),$

where $\tilde{\varphi}(k-\bar{k}_{z})$ is the Fourier transform of
$\varphi(z)$,
$V_{1D}(x,\bar{k}_{z})=\sum_{k_{1}}|\tilde{\varphi}(k_{1}-\bar{k}_{z})|^{2}V(x,p_{z})$
is the effective 1D potential and can be written in the form of

$V_{1D}(x,\bar{k}_{z})=\sum_{m=0}^{\infty}V_{m}(\bar{k}_{z})\cos{2mk_{L}x},$
(5)

where $V_{m}(\bar{k}_{z})$ is the $m$-th order coefficient that depends on
$\bar{k}_{z}$. We find truncating at $m=4$ sufficient with $a_{0}=4$, as the
coefficients $V_{m}$ with $m\geq 3$ are much smaller than the leading orders.
Notice that we ignored the constant energy offset,
$\hbar\omega_{\perp}+\hbar^{2}\bar{k}_{z}^{2}/2M$.

Next, we consider the 3D many-body Hamiltonian
$\mathcal{H}_{3D}=H_{0}+H_{int}$, where $H_{0}=\int
d^{3}r\hat{\psi}^{\dagger}(r)H_{L}\hat{\psi}(r)$ and $H_{int}=\frac{g}{2}\int
d^{3}r\int
d^{3}r^{\prime}\hat{\psi}^{\dagger}(r)\hat{\psi}^{\dagger}(r^{\prime})\delta(r-r^{\prime})\hat{\psi}(r^{\prime})\hat{\psi}(r)$.
The bosonic annihilation operator at $j$th lattice site is defined as
$\hat{\psi}(r)=\sum_{j}W_{j}(x)\varphi(y)\varphi(z)\hat{a}_{j}$ and
$g=4\pi\hbar^{2}a_{s}/M$, where $a_{s}$ is the $s$-wave scattering length.
$W_{j}(x)$ and $\hat{a}_{j}$ is the Wannier function and dimensionless
annihilation operator of $j$th site for the effective 1D Hamiltonian,
$H_{1D}$. We use the Gaussian approximation of the Wannier functions of each
site for the lattice potential, $V_{1D}(x,\bar{k}_{z})$. Combining eqn.(4)
with the 3D many-body Hamiltonian, we obtain the many-body effective 1D
Hamiltonian. Noticing that the next-nearest-neighbor tunneling is negligible
in a certain regime of $\langle p_{z}\rangle=\hbar\bar{k}_{z}$, we can apply
the tight-binding approximation and obtain the effective 1D Bose-Hubbard
Hamiltonian,

$\displaystyle H_{BH}=$
$\displaystyle-\sum_{j}J_{1D}(\hat{a}_{j+1}^{\dagger}\hat{a}_{j}+\hat{a}_{j+1}\hat{a}_{j}^{\dagger})$
(6)
$\displaystyle+\frac{U_{1D}}{2}\sum_{j}\hat{n}_{j}(\hat{n}_{j}-1)-\sum_{j}\tilde{\mu}\hat{n}_{j},$

where $\hat{n}_{j}=\hat{a}_{j}^{\dagger}\hat{a}_{j}$ is the number operator of
$j$th site and $\tilde{\mu}_{j}=\mu-\varepsilon_{j}$. The effective tunneling,
$J_{1D}$, on-site interaction, $U_{1D}$, and energy offset, $\varepsilon_{j}$,
take the forms

$\displaystyle J_{1D}(\bar{k}_{z})=-\int
dxW_{j}^{\ast}(x)H_{1D}(x,\bar{k}_{z})W_{j+1}(x),$ $\displaystyle
U_{1D}(\bar{k}_{z})=g_{1D}\int dx|W_{j}(x)|^{4},$ (7)
$\displaystyle\varepsilon_{j}(\bar{k}_{z})=-\int
dxW_{j}^{\ast}(x)H_{1D}(x,\bar{k}_{z})W_{j}(x),$

where $g_{1D}=2\hbar\omega_{\perp}a_{s}/(1-Aa_{s}/a_{\perp})$ with $A=1.036$
Bloch  _et al._ (2008). We can drop the site index $j$ so that
$\tilde{\mu}_{j}$ becomes $\tilde{\mu}$, since the lattice is assumed
isotropic.

Figure 1: Effective 1D on-site interaction, $U_{1D}$ (blue line), and
effective nearest-neighbor tunneling, $J_{1D}$ (red line), with respect to the
average transverse momentum, $\hbar\bar{k}_{z}$.

Eqn.(Many-body Hamiltonian.–) shows that both $J_{1D}$ and $U_{1D}$ depend on
the average transverse momentum, $\hbar\bar{k}_{z}$. In Fig.1, we show the
plot of $J_{1D}$ and $U_{1D}$ as functions of average transverse momentum.

### Phase diagram.–

In the effective 1D regime, the on-site interaction strength is much larger
than the nearest-neighbor tunneling so that we can use the strong coupling
expansion to calculate the SF-MI phase diagram of the Bose-Hubbard model
Kühner and Monien (1998); Ejima  _et al._ (2012); et al. (2018) and get a
$\langle p_{z}\rangle$-dependent phase diagram. The solid line in Fig.2 shows
the $\bar{n}=1$ Mott lobe of the phase diagram. The phase diagram displays the
tunability of the SF-MI phase transition upon the transverse average momentum
$\langle p_{z}\rangle$ in TMDOL. The detailed parameters we use are chosen for
87Rb atoms. Specifically, $k_{L}=2\pi/\lambda_{R}$, $\lambda_{R}\approx
791\textrm{nm}=2a_{L}$, $k_{\perp}=2\pi/\lambda_{\perp}$, and
$\lambda_{\perp}=852\textrm{nm}$ is the wavelength for the 2D optical lattice
along transverse directions.

Next, we use path-integral Monte Carlo calculations in continuous space in the
grand-canonical ensemble Ceperley (1995) to obtain the QMC phase diagram. We
consider the effective 1D many-body Hamiltonian in the continuous space,

$H_{QMC}=\sum_{1\leq j\leq
N}\Big{[}-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x_{j}^{2}}+V_{1D}(x_{j},\bar{k}_{z})\Big{]}+g_{1D}\sum_{j<\ell}\delta(x_{j}-x_{\ell}),$
(8)

where $V_{1D}(x,\bar{k}_{z})$ is the effective potential with the form of
eqn.(5). For a given temperature, $T$, chemical potential, $\mu$, and external
potential, $V_{1D}(x,\bar{k}_{z})$, we can generate the imaginary-time Feynman
diagram configurations with a certain probability distribution. Then, the
particle density, $\bar{n}=N/L$, is obtained from the statistics of the closed
worldlines, where $N$ denotes the atom number in the ensemble and $L$ the
system size. Thanks to the worm algorithm implementations Boninsegni  _et al._
(2006a, b), we can further compute the superfluid density, $\rho_{s}$,
efficiently. These two quantities help us to distinguish the MI from the SF
phase. We use the same QMC algorithm as in Refs. Boéris  _et al._ (2016); Yao
_et al._ (2018, 2020); Gautier  _et al._ (2021), where details of the
implementations are discussed. In the actual QMC calculations, we take a large
enough size, $L=50a_{L}$, and small enough temperature, $k_{B}T=0.004E_{r}$,
to make sure our system is in thermodynamic and zero-temperature limits. We
confirm that there is no further finite-size and finite-temperature effects at
larger sizes and smaller temperatures.

Figure 2: Phase diagram of the effective 1D TMDOL. The pink shaded area with a
black boundary is the $\bar{n}=1$ Mott lobe calculated from the strong
coupling expansion of the effective 1D Bose-Hubbard model. Diamonds with
errorbars are QMC results with system size $L=50a_{L}$ and temperature
$k_{B}T=0.004E_{r}$. Note that all errorbars are $\pm 0.005E_{r}$. The inset
plot shows the effective 1D lattice potential for $\bar{k}_{z}=-4.73k_{L}$
(green line) and $\bar{k}_{z}=-2.49k_{L}$ (violet line), which correspond to
the potentials for the green diamond and the violet diamond data points,
respectively. Blue diamonds with errorbars are QMC results with the other
momentum.

Our final QMC results are shown in Fig.2 as the diamonds with errorbars, which
indicates the SF-MI transition boundaries for corresponding parameters. We
generally find good agreements between the strong coupling expansion and the
QMC calculation for the $\bar{n}=1$ Mott lobe. Despite the small difference
for small $|k_{z}|$, both methods confirms a significant Mott lobe for
$\bar{n}=1$, which verifies that there exists a TMD SF-MI transition. As
$\bar{k}_{z}$ gets closer to zero, the strong coupling expansion and the QMC
results start to deviate from each other. As shown by the inset plot of Fig.2,
the lattice potential becomes shallower when $|\bar{k}_{z}|$ decreases, which
leads to the increase of $J_{1D}/U_{1D}$, and the next-nearest neighbor
tunneling becomes non-negligible. Therefore, the Bose-Hubbard model becomes
less accurate and the deviation between two methods arises. Moreover, the
$\bar{n}=2$ Mott lobe is shown in the supplemental material sup , and we
believe that the bigger deviation between two methods is due to the density-
dependent ground state Lühmann  _et al._ (2012); Dutta  _et al._ (2011);
Bissbort  _et al._ (2012). Despite the slight differences, generally the TMD
SF-MI transition is still confirmed by both approaches, which displays the
many-body effect caused by the artificial gauge potential in TMDOL.

Figure 3: (a) Level diagram of the Raman process. Red solid, blue solid and
red dashed arrows indicate lasers $a$, $b$, and $c$, respectively. (b) Laser
configuration of the Raman process. The laser colors are the same as (a), with
their polarizations indicated by the small arrows. The bias magnetic field is
in the $+z$ direction, and the four black arrows represent the two-dimensional
trapping lattice beams in $y-z$ plane.

### Experimental implementation.–

Experimentally, the periodically driven Hamiltonian in eqn.(1) can be
constructed by coupling the atoms to laser fields. Specifically, we consider
the dipole interaction of 87Rb atom with continuous weak (CW) laser fields
Wright  _et al._ (2008). The pseudospin-1/2 system is formed by two Zeeman
states in the $5^{2}S_{1/2},F=1$ ground state manifold, namely
$\ket{\uparrow}=\ket{F=1,m_{F}=-1}$, and $\ket{\downarrow}=\ket{F=1,m_{F}=0}$.
In the presence of a bias magnetic field, $\vec{B}$, which defines the
$z$-axis of the system, the degeneracy between the two states is lifted. Both
states are coupled via a Raman process by applying CW Raman lasers $a$, $b$
and $c$, where $a$ and $c$ are standing waves along the $x$-axis with
$\pi$-polarization and $b$ propagates in the $+z$-direction with
$\sigma^{-}$-polarization. The Raman lasers are detuned far enough from the
$D_{1}$ and $D_{2}$ transitions to avoid any excitations over the timescale of
the experiment. The states we considered in our calculations and the laser
configurations are shown in Fig.3. Laser coupling to all other states will be
negligible.

The sinusoidal driving of the Hamiltonian can be achieved by using a elctro-
optical modulator (EOM) to drive the intensity of Raman lasers $a$ and $c$
with $|\cos{\omega t}|$, and at the same time modulating the relative phase,
$\eta$, between lasers $a$ and laser $b$ such that
$\eta(t)=\pi\left[1-sgn(\cos{\omega t})\right]$, where $sgn(\cdot)$ denotes
the signum function. The electric fields of the Raman lasers take the form

$\displaystyle\vec{E}_{a}=E_{a}\hat{\varepsilon}_{0}\cos{k_{L}x}|\cos{\omega
t}|e^{-i\left[\omega_{a}t+\eta(t)\right]}$
$\displaystyle\vec{E}_{b}=E_{b}\hat{\varepsilon}_{-}e^{-i(k_{L}z+\omega_{b}t)}$
$\displaystyle\vec{E}_{c}=E_{c}\hat{\varepsilon}_{0}\cos{k_{L}x}|\cos{\omega
t}|e^{-i\omega_{c}t},$

where $E_{j}$ and $\omega_{j}$ ($j=a,b,c$) is the amplitude and the angular
frequency of electric field of each laser. $\hat{\varepsilon}_{0}$,
$\hat{\varepsilon}_{+}$ and $\hat{\varepsilon}_{-}$ are $\pi$, $\sigma^{+}$
and $\sigma^{-}$ polarization vectors, respectively. Using second order
perturbation theoryCohen-Tannoudji  _et al._ (1998); Wright  _et al._ (2008);
sup , the Raman coupling Hamiltonian can be made into the desired form. In our
calculations, we apply $\hbar\omega=20E_{r}$ with Raman laser wavelengths
close to $791\textrm{nm}$ and $806\textrm{nm}$, which are achievable in the
laboratory. These parameters result in a total laser power of around
$0.8\textrm{W}$ with Gaussian laser beams with beam waists of $200\upmu$m, and
a scattering limited lifetime of $\tau_{sc}\approx 133\textrm{ms}$, which is
enough for experimental applications.

To load the system into the ground state of the TMDOL, we can start from a
cigar-shaped Bose–Einstein condensate with the desired atom number in an
optical dipole trap with the long axis of the trap aligned with $x$-axis in
the lab frame. Then, we adiabatically ramp up the two-dimensional optical
lattice along the $y$ and $z$-axes to avoid excitation to higher bands. Next,
applying the shortcut adiabatic loading technique Torrontegui  _et al._
(2011); Corgier  _et al._ (2018); Guéry-Odelin  _et al._ (2019); Ding  _et
al._ (2020) to load the Bose gas into the ground state of a moving 2D optical
lattice. Due to the relative motion between the lattice along the $z$-axis and
the lab frame, the atoms will carry a transverse momentum with expectation
value $\langle p_{z}\rangle=-\hbar\bar{k}_{z}$ in the moving lattice frame.
Then, adiabatically turn on the Raman lasers and the system will have been
loaded into the ground state of the TMDOL with a certain $\langle
p_{z}\rangle$. To make observations, one can suddenly shut off all lasers at
the end of the periodic driving cycle when the micromotion operator, $U$,
becomes the identity operator so that the system is automatically projected
from the Floquet basis to the Zeeman sublevel basisChen  _et al._ (2020).
Then, a time-of-flight (TOF) image can be taken with the imaging beam aligned
along the $z$-axis. The interference pattern contrast, $\mathcal{I}$, can be
measured from the TOF image Gerbier  _et al._ (2005); Soltan-Panahi  _et al._
(2011), which signals SF and MI phases.

### Conclusion.–

To conclude, we proposed an experimentally feasible method to create an
effective 1D TMD optical lattice induced by an artificial gauge potential in a
periodically driven system with Raman couplings. We explored the many-body
physics of the system induced by a periodic artificial gauge potential by
constructing the TMD effective 1D Bose-Hubbard model with a deep 2D transverse
confinement. Using the strong coupling expansion and QMC methods, we
calculated the phase diagram of an effective 1D TMDOL, and showed that the SF-
MI transition can be tuned by changing the average transverse momentum of the
Bose gas. Then, we proposed a possible experimental implementation that
achieves the desired periodically driven Hamiltonian and includes loading and
measurement procedures.

Our work provides a new platform for studying the many-body physics of a cold
atom system with an artificial gauge potential. For instance, the TMDOL makes
it possible to create two spatially overlapping BECs of the same species in
different quantum phases, which gives rise to opportunities to study if the
two BECs, one SF and one MI, will exchange their phases during a collision in
which they exchange their momenta. Additionally, more complicated TMDOL
potentials may be considered by engineering the transverse momentum
distribution of the Bose gas and introducing Raman lasers along other
directions to achieve a 2D TMDOL. The transverse momentum in a TMDOL can also
be regarded as a synthetic dimension, which is analogous to spin or harmonic
trap eigenstates in previous work Celi  _et al._ (2014); Taddia  _et al._
(2017); Price  _et al._ (2017). Finally, for fermionic systems, one could also
implement a TMDOL with appropriate Raman couplings, and such a system provides
a broad playground for studying TMD spin dynamics and many-body phases.

###### Acknowledgements.

We thank Thierry Giamarchi for insightful discussions. This work is supported
by NSF grant PHY 1708008, NASA/JPL RSA 1656126 and the Swiss National Science
Foundation under Division II. Numerical calculations make use of the ALPS
scheduler library and statistical analysis tools Troyer  _et al._ (1998); et
al. (2007, 2011).

## References

  * Sugawa  _et al._ (2018) S. Sugawa, F. Salces-Carcoba, A. R. Perry, Y. Yue, and I. B. Spielman, Science 360, 1429–1434 (2018).
  * Zhu  _et al._ (2006) S.-L. Zhu, H. Fu, C.-J. Wu, S.-C. Zhang, and L.-M. Duan, Phys. Rev. Lett. 97, 240401 (2006).
  * Liu  _et al._ (2007) X.-J. Liu, X. Liu, L. C. Kwek, and C. H. Oh, Phys. Rev. Lett. 98, 026602 (2007).
  * LeBlanc  _et al._ (2012) L. J. LeBlanc, K. Jiménez-García, R. A. Williams, M. C. Beeler, A. R. Perry, W. D. Phillips, and I. B. Spielman, Proceedings of the National Academy of Sciences 109, 10811–10814 (2012).
  * Beeler  _et al._ (2013) M. C. Beeler, R. A. Williams, K. Jimenez-Garcia, L. J. LeBlanc, A. R. Perry, and I. B. Spielman, Nature 498, 201–204 (2013).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Goldman  _et al._ (2014) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Reports on Progress in Physics 77, 126401 (2014).
  * Zhai (2015) H. Zhai, Reports on Progress in Physics 78, 026001 (2015).
  * Zhang and Liu (2018) L. Zhang and X.-J. Liu, in  _Synthetic Spin-Orbit Coupling in Cold Atoms_ (World Scientific, 2018) pp. 1–87.
  * Huang  _et al._ (2016) L. Huang, Z. Meng, P. Wang, P. Peng, S.-L. Zhang, L. Chen, D. Li, Q. Zhou, and J. Zhang, Nature Physics 12, 540–544 (2016).
  * Hamner  _et al._ (2015) C. Hamner, Y. Zhang, M. A. Khamehchi, M. J. Davis, and P. Engels, Phys. Rev. Lett. 114, 070401 (2015).
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83–86 (2011).
  * Liu  _et al._ (2013) X.-J. Liu, Z.-X. Liu, and M. Cheng, Phys. Rev. Lett. 110, 076401 (2013).
  * Wu  _et al._ (2016) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Science 354, 83–88 (2016).
  * et al. (2021) Z.-Y. Wang et al., Science 372, 271–276 (2021).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
  * Boéris  _et al._ (2016) G. Boéris, L. Gori, M. D. Hoogerland, A. Kumar, E. Lucioni, L. Tanzi, M. Inguscio, T. Giamarchi, C. D’Errico, G. Carleo, G. Modugno, and L. Sanchez-Palencia, Phys. Rev. A 93, 011601(R) (2016).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415, 39–44 (2002).
  * Haller  _et al._ (2010) E. Haller, R. Hart, M. J. Mark, J. G. Danzl, L. Reichsöllner, M. Gustavsson, M. Dalmonte, G. Pupillo, and H.-C. Nägerl, Nature 466, 597–600 (2010).
  * Soltan-Panahi  _et al._ (2011) P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nature Physics 7, 434–440 (2011).
  * Cai  _et al._ (2012) Z. Cai, X. Zhou, and C. Wu, Phys. Rev. A 85, 061605 (2012).
  * Cole  _et al._ (2012) W. S. Cole, S. Zhang, A. Paramekanti, and N. Trivedi, Phys. Rev. Lett. 109, 085302 (2012).
  * Račkauskas  _et al._ (2019) P. Račkauskas, V. Novičenko, H. Pu, and G. Juzeliūnas, Phys. Rev. A 100, 063616 (2019).
  * Novičenko and Juzeliūnas (2019) V. Novičenko and G. Juzeliūnas, Phys. Rev. A 100, 012127 (2019).
  * Chen  _et al._ (2020) Z. Chen, J. D. Murphree, and N. P. Bigelow, Phys. Rev. A 101, 013606 (2020).
  * (27) See supplemental material. We discuss calculation of Floquet Hamiltonian and the adiabatic conditions, the single particle ground state of the TMDOL, the $\bar{n}=2$ Mott lobe, and the detailed construction of the periodically driven Hamiltonian.
  * Stöferle  _et al._ (2004) T. Stöferle, H. Moritz, C. Schori, M. Köhl, and T. Esslinger, Phys. Rev. Lett. 92, 130403 (2004).
  * Anisimovas  _et al._ (2016) E. Anisimovas, M. Račiūnas, C. Sträter, A. Eckardt, I. B. Spielman, and G. Juzeliūnas, Phys. Rev. A 94, 063632 (2016).
  * Barbiero and Salasnich (2014) L. Barbiero and L. Salasnich, Phys. Rev. A 89, 063605 (2014).
  * Kühner and Monien (1998) T. D. Kühner and H. Monien, Phys. Rev. B 58, R14741 (1998).
  * Ejima  _et al._ (2012) S. Ejima, H. Fehske, F. Gebhard, K. zu Münster, M. Knap, E. Arrigoni, and W. von der Linden, Phys. Rev. A 85, 053644 (2012).
  * et al. (2018) T. Wang et al., Phys. Rev. B 98, 245107 (2018).
  * Ceperley (1995) D. M. Ceperley, Rev. Mod. Phys. 67, 279–355 (1995).
  * Boninsegni  _et al._ (2006a) M. Boninsegni, N. Prokof’ev, and B. Svistunov, Phys. Rev. Lett. 96, 070601 (2006a).
  * Boninsegni  _et al._ (2006b) M. Boninsegni, N. V. Prokof’ev, and B. V. Svistunov, Phys. Rev. E 74, 036701 (2006b).
  * Yao  _et al._ (2018) H. Yao, D. Clément, A. Minguzzi, P. Vignolo, and L. Sanchez-Palencia, Phys. Rev. Lett. 121, 220402 (2018).
  * Yao  _et al._ (2020) H. Yao, T. Giamarchi, and L. Sanchez-Palencia, Phys. Rev. Lett. 125, 060401 (2020).
  * Gautier  _et al._ (2021) R. Gautier, H. Yao, and L. Sanchez-Palencia, Phys. Rev. Lett. 126, 110401 (2021).
  * Lühmann  _et al._ (2012) D.-S. Lühmann, O. Jürgensen, and K. Sengstock, New Journal of Physics 14, 033021 (2012).
  * Dutta  _et al._ (2011) O. Dutta, A. Eckardt, P. Hauke, B. Malomed, and M. Lewenstein, New Journal of Physics 13, 023019 (2011).
  * Bissbort  _et al._ (2012) U. Bissbort, F. Deuretzbacher, and W. Hofstetter, Phys. Rev. A 86, 023617 (2012).
  * Wright  _et al._ (2008) K. C. Wright, L. S. Leslie, and N. P. Bigelow, Phys. Rev. A 78, 053412 (2008).
  * Cohen-Tannoudji  _et al._ (1998) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  _Atom-photon interactions: basic processes and applications_ (1998).
  * Torrontegui  _et al._ (2011) E. Torrontegui, S. Ibáñez, X. Chen, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. A 83, 013415 (2011).
  * Corgier  _et al._ (2018) R. Corgier, S. Amri, W. Herr, H. Ahlers, J. Rudolph, D. Guéry-Odelin, E. M. Rasel, E. Charron, and N. Gaaloul, New Journal of Physics 20, 055002 (2018).
  * Guéry-Odelin  _et al._ (2019) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).
  * Ding  _et al._ (2020) Y. Ding, T.-Y. Huang, K. Paul, M. Hao, and X. Chen, Phys. Rev. A 101, 063410 (2020).
  * Gerbier  _et al._ (2005) F. Gerbier, A. Widera, S. Fölling, O. Mandel, T. Gericke, and I. Bloch, Phys. Rev. A 72, 053606 (2005).
  * Celi  _et al._ (2014) A. Celi, P. Massignan, J. Ruseckas, N. Goldman, I. B. Spielman, G. Juzeliūnas, and M. Lewenstein, Phys. Rev. Lett. 112, 043001 (2014).
  * Taddia  _et al._ (2017) L. Taddia, E. Cornfeld, D. Rossini, L. Mazza, E. Sela, and R. Fazio, Phys. Rev. Lett. 118, 230402 (2017).
  * Price  _et al._ (2017) H. M. Price, T. Ozawa, and N. Goldman, Phys. Rev. A 95, 023607 (2017).
  * Troyer  _et al._ (1998) M. Troyer, B. Ammon, and E. Heeb, Lect. Notes Comput. Sci. 1505, 191 (1998).
  * et al. (2007) A. F. Albuquerque et al., J. Magn. Magn. Mater. 310, 1187 (2007).
  * et al. (2011) B. Bauer et al., J. Stat. Mech.: Th. Exp. 05, P05001 (2011).
  * Grimm  _et al._ (2000) R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov, Advances in atomic, molecular, and optical physics 42, 95–170 (2000).

Supplemental Material for

A momentum dependent optical lattice induced by artificial gauge potential

In this supplemental material, we provide details about zeroth-order component
of the Floquet Hamiltonian and the adiabatic condition calculations (Sec. S1),
the single particle ground state of the TMD optical lattice (Sec. S2), the
$\bar{n}=2$ Mott lobe calculation result (Sec. S3), and the detailed
construction of the periodically driven Hamiltonian via a Raman process (Sec.
S4).

##  S1 Zeroth-order component of the Floquet Hamiltonian and the adiabatic
condition

In this work, the system we consider satisfies the Floquet adiabatic condition
in the theoretical formalism of the main text. We can rewrite eqn.(2) in the
main text in the form

$H_{F}=\frac{\vec{p}^{2}}{2M}-\frac{1}{2M}(\vec{A}\cdot\vec{p}+\vec{p}\cdot\vec{A})+\frac{\vec{A}^{2}}{2M}+V_{ext}(\vec{r}).$
(S1)

The adiabatic condition can be rewritten as

$|\langle\vec{k}|[\frac{1}{2M}(\vec{A}\cdot\vec{p}+\vec{p}\cdot\vec{A})+\frac{\vec{A}^{2}}{2M}]^{(n)}|\vec{k}^{{}^{\prime}}\rangle|\ll\hbar\omega(n\neq
0),$ (S2)

where $[\cdots]^{(n)}$ is the $n$-th Fourier component. From eqn.(S2), we can
see that the adiabatic condition depends not only on the ratio
$\Omega_{0}/\omega$, but also on the momentum of the particle. A rigorous
solution of the adiabatic condition can be lengthy, so we focus on a
operational condition.

We first work on the term $\vec{A}^{2}/2M$. Writing down the analytic form of
the vector potential, $\vec{A}$ Račkauskas  _et al._ (2019); Novičenko and
Juzeliūnas (2019); Chen  _et al._ (2020):

$A_{\xi}=d_{1\xi}+d_{2\xi}+d_{3\xi},$ (S3)

where

$\displaystyle d_{1\xi}=\hbar a\sin{\omega
t}\frac{(\vec{\Omega}\cdot\partial_{\xi}\vec{\Omega})(\vec{\Omega}\cdot\vec{\sigma})}{2|\Omega|^{3}},$
$\displaystyle d_{2\xi}=\hbar\sin{(a\sin{\omega
t})}\frac{[(\vec{\Omega}\times\partial_{\xi}\vec{\Omega})\times\vec{\Omega}]\cdot\vec{\sigma}}{2|\Omega|^{3}},$
(S4) $\displaystyle d_{3\xi}=\hbar[\cos{(a\sin{\omega
t})}-1]\frac{\vec{\Omega}\times\partial_{\xi}\vec{\Omega}\cdot\vec{\sigma}}{2|\Omega|^{2}},$

are the $\xi$-th component of vectors $\vec{d}_{1}$, $\vec{d}_{2}$ and
$\vec{d}_{3}$. Then,

$A_{\xi}^{2}=d_{1\xi}^{2}+d_{2\xi}^{2}+d_{3\xi}^{2}+(d_{1\xi}d_{2\xi}+d_{2\xi}d_{1\xi})+(d_{1\xi}d_{3\xi}+d_{3\xi}d_{1\xi})+(d_{2\xi}d_{3\xi}+d_{3\xi}d_{2\xi}).$
(S5)

We can simplify eqn.(S5) by noticing that the last three terms in the
parenthesis vanish:

$\displaystyle(d_{1\xi}d_{2\xi}+d_{2\xi}d_{1\xi})$ $\displaystyle\propto$
$\displaystyle\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\hat{\sigma}_{j}\epsilon_{klm}C_{\xi
k}\Omega_{l}\hat{\sigma}_{m}+\epsilon_{klm}C_{\xi
k}\Omega_{l}\hat{\sigma}_{m}\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\hat{\sigma}_{j}$
$\displaystyle=$
$\displaystyle\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\epsilon_{klm}C_{\xi
k}\Omega_{l}(\hat{\sigma}_{j}\hat{\sigma}_{m}+\hat{\sigma}_{m}\hat{\sigma}_{j})$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle(d_{2\xi}d_{3\xi}+d_{3\xi}d_{2\xi})$ $\displaystyle\propto$
$\displaystyle\epsilon_{ijk}C_{\xi i}\Omega_{j}\hat{\sigma}_{k}C_{\xi
l}\hat{\sigma}_{l}+C_{\xi l}\hat{\sigma}_{l}\epsilon_{ijk}C_{\xi
i}\Omega_{j}\hat{\sigma}_{k}$ (S6) $\displaystyle=$
$\displaystyle\epsilon_{ijk}C_{\xi i}\Omega_{j}C_{\xi
l}(\hat{\sigma}_{k}\hat{\sigma}_{l}+\hat{\sigma}_{l}\hat{\sigma}_{k})$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle(d_{1\xi}d_{3\xi}+d_{3\xi}d_{1\xi})$ $\displaystyle\propto$
$\displaystyle\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\hat{\sigma}_{j}C_{\xi
k}\hat{\sigma}_{k}+C_{\xi
k}\hat{\sigma}_{k}\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\hat{\sigma}_{j}$
$\displaystyle=$
$\displaystyle\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}C_{\xi
k}(\hat{\sigma}_{j}\hat{\sigma}_{k}+\hat{\sigma}_{k}\hat{\sigma}_{j})$
$\displaystyle=$ $\displaystyle
2\Omega_{i}\partial_{\xi}\Omega_{i}\Omega_{j}\epsilon_{jlm}\Omega_{l}\partial_{\xi}\Omega_{m}$
$\displaystyle=$ $\displaystyle 0,$

where we let $C_{\xi
i}=\\{\vec{\Omega}\times\partial_{\xi}\vec{\Omega}\\}_{i}=\epsilon_{ijk}\Omega_{j}\partial_{\xi}\Omega_{k}$
and we used the property
$\hat{\sigma}_{i}\hat{\sigma}_{j}+\hat{\sigma}_{j}\hat{\sigma}_{i}=2\delta_{ij}\mathbbm{1}$.

After the simplification above, eqn.(S5) becomes

$A_{\xi}^{2}=\hbar^{2}a^{2}\sin^{2}{\omega
t}\frac{(\vec{\Omega}\cdot\partial_{\xi}\vec{\Omega})^{2}}{4|\Omega|^{4}}+2\hbar^{2}[1-\cos{(a\sin{\omega
t)}}]\frac{|\Omega|^{2}(\partial_{\xi}\vec{\Omega})^{2}-(\vec{\Omega}\cdot\partial_{\xi}\vec{\Omega})^{2}}{4|\Omega|^{4}}.$
(S7)

We can write $\cos{(a\sin{\omega t})}$ as an expansion of Bessel functions of
the first kind,

$\cos{(a\sin{\omega t})}=J_{0}(a)+2\sum_{m=1}^{\infty}J_{2m}(a)\cos{(2m\omega
t)}.$ (S8)

The zeroth order Fourier term of $\vec{A}^{2}/2M$ is

$\left[\frac{\vec{A}^{2}}{2M}\right]^{(0)}=\\{\frac{1}{16}a_{0}^{2}(1-\cos{2k_{L}x})+\frac{1}{2}\left[1-J_{0}(a)\right]\\}E_{r},$
(S9)

and the second order ($\pm 2$) harmonic terms are

$\left[\frac{\vec{A}^{2}}{2M}\right]^{(\pm
2)}=[\frac{1}{32}a_{0}^{2}\cos{2k_{L}x}\pm\frac{i}{2}J_{2}(a)]E_{r}.$ (S10)

Next, we calculate the contribution of the
$\frac{1}{2M}(\vec{A}\cdot\vec{p}+\vec{p}\cdot\vec{A})$ term. From eqn.(S1),
we get

$\displaystyle\vec{A}$ $\displaystyle=$ $\displaystyle\frac{1}{2}\hbar
k_{L}\\{[-a_{0}\sin{k_{L}x}e^{-ik_{L}z}\sin{\omega
t}\hat{x}+ie^{ik_{L}z}\sin{(a\sin{\omega
t})}sgn(a)\hat{z}]\hat{\sigma}_{+}+h.c.\\}$ (S11) $\displaystyle+$
$\displaystyle\frac{1}{2}\hbar k_{L}\left[\cos{(a\sin{\omega
t})}-1\right]\hat{\sigma}_{3}\hat{z},$

where $sgn(\cdot)$ is the signum function,
$\hat{\sigma}_{+}=(\hat{\sigma}_{1}-i\hat{\sigma}_{2})/\sqrt{2}$, $\hat{x}$
and $\hat{z}$ are unit vector along the $x$ and $z$-axes. Then, we can write
down the absolute value of the non-zero matrix element of the Floquet
Hamiltonian that is proportional to the first Fourier order,

$\displaystyle\bra{k_{x},k_{z},\uparrow}H_{F}^{(\pm
1)}\ket{k_{x}^{{}^{\prime}},k_{z}^{{}^{\prime}},\downarrow}$ $\displaystyle=$
$\displaystyle\pm\frac{a_{0}}{8}E_{r}\\{\bra{k_{x},k_{z},\uparrow}[(e^{ik_{L}x}-e^{-ik_{L}x})e^{-ik_{L}z}\frac{k_{x}^{{}^{\prime}}}{k_{L}}]\ket{k_{x}^{{}^{\prime}},k_{z}^{{}^{\prime}},\downarrow}$
$\displaystyle+$
$\displaystyle\bra{k_{x},k_{z},\uparrow}[\frac{k_{x}}{k_{L}}(e^{ik_{L}x}-e^{-ik_{L}x})e^{ik_{L}z}]\ket{k_{x}^{{}^{\prime}},k_{z}^{{}^{\prime}},\downarrow}\\}$
$\displaystyle\pm$
$\displaystyle\frac{E_{r}}{2}\bra{k_{x},k_{z},\uparrow}J_{1}(a)sgn(a)(e^{ik_{L}z}\frac{k_{z}}{k_{L}}\mp
e^{-ik_{L}z}\frac{k_{z}^{{}^{\prime}}}{k_{L}})\ket{k_{x}^{{}^{\prime}},k_{z}^{{}^{\prime}},\downarrow},$

where $\ket{\uparrow(\downarrow)}$ indicates two spin states. Notice that the
above equation can be split into several matrix elements between different
momentum states. After collecting terms and simplifying, we get the adiabatic
condition for the first Fourier order as

$\\{\frac{a_{0}}{8}|\frac{k_{x}}{k_{L}}|+\frac{1}{2}|\langle
k_{L}|sgn(a)J_{1}(a)|k_{L}\pm
1\rangle||\frac{k_{z}}{k_{L}}|\\}E_{r}=\frac{1}{2}\\{|\frac{k_{x}}{k_{L}}|+0.19|\frac{k_{z}}{k_{L}}|\\}E_{r}\ll\hbar\omega,$
(S12)

where we evaluated $|\langle k_{L}|sgn(a)J_{1}(a)|k_{L}\pm 1\rangle|$
numerically to be 0.19 by setting $a_{0}=4$. In our calculations, we use
$\hbar\omega/E_{r}=20$, and we can get the adiabatic condition for the first
Fourier order by combining the two conditions in eqn.(S12) to be
$|k_{x}/k_{L}|+0.19|k_{z}/k_{L}|\ll 40$. We also evaluated the matrix elements
between higher order momentum states and found them to be much smaller than
the matrix element in eqn.(S12).

Next, we consider the second order Fourier components. Both the
$\frac{1}{2M}(\vec{A}\cdot\vec{p}+\vec{p}\cdot\vec{A})$ and $A^{2}$ terms
contribute to the even order $\sin{2m\omega t}$ terms. Specifically, for the
second order, from eqn.(S10) and eqn.(S11) we can get the matrix element of
second Fourier component,

$\displaystyle\langle k_{x},k_{z},\uparrow(\downarrow)|H_{F}^{(\pm
2)}|k_{x}^{{}^{\prime}},k_{z},\uparrow(\downarrow)\rangle$ $\displaystyle=$
$\displaystyle\langle
k_{x},k_{z},\uparrow(\downarrow)|\left[\frac{1}{64}a_{0}^{2}\pm\frac{i}{2}J_{2}(a)\right]E_{r}|k_{x}^{{}^{\prime}}\pm
2k_{L},k_{z},\uparrow(\downarrow)\rangle$ $\displaystyle+$
$\displaystyle\langle
k_{x},k_{z},\uparrow(\downarrow)|2\frac{k_{z}}{k_{L}}J_{2}(a)E_{r}|k_{x}^{{}^{\prime}}\pm
2k_{L},k_{z},\uparrow(\downarrow)\rangle$ $\displaystyle\approx$
$\displaystyle(\frac{1}{64}a_{0}^{2}\pm
0.1i-0.4\frac{k_{z}}{k_{L}})E_{r}\delta(k_{x}-k_{x}^{{}^{\prime}}\pm 2k_{L}).$

where we used $\langle k_{x}|J_{2}(a)|k_{x}\pm 2k_{L}\rangle=-0.20$ from
numerical evaluation using the same parameters as before.The adiabatic
condition for the second Fourier order becomes $|(\frac{1}{64}a_{0}^{2}\pm
0.1i-0.4\frac{k_{z}}{k_{L}})E_{r}|\ll\hbar\omega$, and we get
$|k_{z}/k_{L}|\ll 49$. Combine this result with
$|k_{x}/k_{L}|+0.19|k_{z}/k_{L}|\ll 40$, we get $|k_{x}/k_{L}|\ll 31$.

For higher order Fourier terms, one can do the calculation and see that the
adiabatic conditions are all weaker than the first two Fourier terms. In our
calculations, the above adiabatic conditions, $|k_{x}/k_{L}|\ll 31$ and
$|k_{z}/k_{L}|\ll 49$, can easily be satisfied.

##  S2 Single particle ground state of the TMD optical lattice

In the main text, the single particle Hamiltonian takes the form

$H_{L}=\frac{\vec{p}^{2}}{2M}+V(x,p_{z})+V_{ext}(\vec{r}),$

where
$V(x,p_{z})=\\{(p_{z}+1/2)\left[1-J_{0}(a)\right]+a_{0}^{2}(1-\cos{2k_{L}x})/16\\}E_{r}$.
We can expand Hamiltonian $H_{L}$ in the momentum basis, $\ket{\pm 2l\hbar
k_{L}}$ ($l=0,\pm 1,\pm 2,\cdots$), and diagonalize the Hamiltonian to get the
single-particle band structure and eigenstates. In all our calculations above,
we have used the expansion form of the Bessel function of the first kind,
$J_{q}(a)=\sum_{m=0}^{\infty}(-1)^{m}(a)^{2m+q}/2^{2m+q}\Gamma(m+1)\Gamma(m+q+1)$,
where $q$ indicates the order of Bessel function. In practice, we truncate the
expansion at the order of $l=\pm 10$, which is accurate enough with our choice
of parameters.

Figure S1: $s$-band relative population of each momentum state along the
$x$-axis with different $p_{z}$ at zero quasi-momentum. The inset plot is the
corresponding $p_{z}$-dependent 1D optical lattice potential.

Consider the simplest case where $V_{ext}(\vec{r})=0$. When $p_{z}$ is low the
on-site interaction in such a 1D optical lattice will be negligible. In such a
weak interaction regime, we can study the physics of the system by solving the
single-particle Hamiltonian, $H_{L}$. As shown in Fig.S1, the 1D lattice
potential is determined by $p_{z}$, which can be experimentally observed by
measuring the population of each momentum state along the $x$-direction.

##  S3 Phase diagram around $\bar{n}=2$ Mott lobe

Figure S2: $\bar{n}=2$ Mott lobe. The pink shaded area with black solid line
boundaries is calculated from the occupation-independent Bose-Hubbard model,
blue diamonds with errorbars ($\pm 0.005Er$) are the QMC calculation results,
and black dashed lines are Bose-Hubbard model calculation results with a band
mixing ground state $|\psi_{g}\rangle\approx 0.994|s\rangle+0.109|d\rangle$.

In Fig.2 of the main text, we showed the phase diagram of our effective 1D
optical lattice and found good agreement on the $\bar{n}=1$ Mott lobe between
our effective 1D Bose-Hubbard model and the QMC calculation. In Fig.S2 we show
the $\bar{n}=2$ Mott lobe. There is a significant difference between the Bose-
Hubbard model (pink shaded area with black solid line boundaries) and the QMC
calculations (blue diamonds with errorbars). We believe that the difference is
caused by the occupation-dependence of the tunneling and interactions terms,
which makes the occupation-independent Bose-Hubbard model less accurate. Our
QMC calculation is more accurate because we used a continuous potential in the
QMC calculations instead of the Bose-Hubbard model. In previous work,
occupation-dependent tunneling and onsite interactions in a modified Bose-
Hubbard model have been studied Lühmann  _et al._ (2012); Dutta  _et al._
(2011); Bissbort  _et al._ (2012), and the true many-body ground state was
found to be a band mixing state. Phenomenologically, in our calculations we
find a mixing state with $|\psi_{g}\rangle\approx
0.994|s\rangle+0.109|d\rangle$ gives a better agreement between a Bose-Hubbard
model with modified parameters and the QMC results (black dashed lines), where
$|s\rangle$ and $|d\rangle$ are the single-particle eigenstates of the $s$ and
$d$-bands, respectively. Our result agrees with the finding in Ref.Lühmann
_et al._ (2012) that the population on higher orbitals is smaller than $1\%$.
Although for more rigorous analysis one would need to numerically diagonalize
the many-body Hamiltonian Lühmann  _et al._ (2012); Dutta  _et al._ (2011);
Bissbort  _et al._ (2012), we believe that our qualitative calculations are
good enough to understand the disagreement on the $\bar{n}=2$ Mott lobe.

##  S4 Detailed construction of the periodically driven Hamiltonian via Raman
process and lifetime estimation

Using second order perturbation theory, we can calculate the effective
Hamiltonian of the ground state manifold of the atom from the multi-laser
Raman process Hamiltonian Cohen-Tannoudji  _et al._ (1998); Wright  _et al._
(2008). Here, we consider 87Rb atoms in the $|5S_{\frac{1}{2}},F=1\rangle$
manifold to be our system. More specifically, we consider
$|1\rangle=|F=1,m_{F}=-1\rangle$, $|2\rangle=|F=1,m_{F}=0\rangle$ and
$|3\rangle=|F=1,m_{F}=1\rangle$. By applying a bias magnetic field, $\vec{B}$,
we can define the $+z$-direction. We then choose our Raman laser $a$ to be a
standing wave along $x$-axis with $\pi$-polarization, laser $b$ propagates
along $+z$-direction with $\sigma^{-}$-polarization, and laser $c$ is chosen
to be a standing wave along the $x$-axis with $\pi$-polarization that is
oppositely detuned compared to laser $a$ to cancel the spatially dependent AC
stark shift. Specifically, lasers $a$ and $b$ are blue detuned from the 87Rb
$D_{1}$ transition and laser $c$ is red detuned from $D_{1}$ transition. All
the wavelengths of the Raman lasers considered here can be achieved in the
laboratory (e.g.with a Ti:Sapphire laser). All one-photon detunings are much
larger than Rabi frequencies so that there is negligible transitions to
excited states. The level diagram is shown in Fig.3 in the main text. We only
consider $D_{1}$ and $D_{2}$ lines and all other atomic states are too far
from the laser frequency we consider so that the laser couplings to those
states are negligible.

Similar to the calculations in previous worksWright  _et al._ (2008); Chen
_et al._ (2020), we can write the reduced 3-level effective Hamiltonian for
the Raman process with the laser fields in eqn.(9) of the main text,

$\displaystyle H_{R}$
$\displaystyle=\sum_{\mathcal{D},F^{\prime}}\sum_{m_{F}=-1}^{1}\\{[\frac{(c_{\mathcal{D},F,F^{\prime}}^{m_{F},m_{F}})^{2}}{\Delta_{a,F,F^{\prime},\mathcal{D}}^{m_{F},m_{F}}}|\Omega_{a,\mathcal{D}}|^{2}+\frac{(c_{\mathcal{D},F,F^{\prime}}^{m_{F}-1,m_{F}})^{2}}{\Delta_{b,F,F^{\prime},\mathcal{D}}^{m_{F},m_{F}}}|\Omega_{b,\mathcal{D}}|^{2}+\frac{(c_{\mathcal{D},F,F^{\prime}}^{m_{F},m_{F}})^{2}}{\Delta_{c,F,F^{\prime},\mathcal{D}}^{m_{F},m_{F}}}|\Omega_{c,\mathcal{D}}|^{2}+\delta_{m_{F}}]|m_{F}\rangle\langle
m_{F}|$ (S14)
$\displaystyle+\frac{1}{2}(\frac{1}{\Delta_{a,F,F^{\prime},\mathcal{D}}^{m_{F},m_{F}}}+\frac{1}{\Delta_{b,F,F^{\prime},\mathcal{D}}^{m_{F}+1,m_{F}}})c_{\mathcal{D},F,F^{\prime}}^{m_{F},m_{F}}c_{\mathcal{D},F,F^{\prime}}^{m_{F}-1,m_{F}}\Omega_{a,\mathcal{D}}^{\ast}\Omega_{b,\mathcal{D}}|m_{F}\rangle\langle
m_{F}+1|+h.c.\\},$

where $F=1$, $\mathcal{D}=1,2$ indicates the transitions for $D_{1}$ and
$D_{2}$ lines, $F^{\prime}=1,2$ and $F^{\prime}=0,1,2$ for $\mathcal{D}=1$ and
$\mathcal{D}=2$, respectively.
$c_{\mathcal{D},F,F^{\prime}}^{m_{F},m_{F^{\prime}}}$ is the Clebsh-Gordon
coefficient between $|F,m_{F}\rangle$ and $|F^{\prime},m_{F^{\prime}}\rangle$
states for $D_{1}$ ($D_{2}$) line and
$\Delta_{\zeta,F,F^{\prime},\mathcal{D}}^{m_{F},m_{F^{\prime}}}=\omega_{\zeta}-(\omega_{F^{\prime},m_{F^{\prime}}}^{e}-\omega_{F,m_{F}}^{g})$
($\zeta=a,b,c$) is the one-photon detuning, where $\omega_{\zeta}$,
$\omega_{F,m_{F}}^{g}$ and $\omega_{F^{\prime},m_{F^{\prime}}^{e}}$ are
angular frequencies for laser $\zeta=a,b,c$ and the corresponding energy
levels of ground state and excited state, respectively.
$\Omega_{\zeta,\mathcal{D}}=-\vec{d}_{\mathcal{D}}\cdot\vec{E}_{\zeta}/\hbar$
($\zeta=a,b,c$) are the Rabi frequencies, where $\vec{d}_{\mathcal{D}}$
($\mathcal{D}=1,2$) are the effective dipole moment vector of $D_{1}$ and
$D_{2}$ transition, respectively. $|m_{F}\rangle=|F=1,m_{F}=-1,0,1\rangle$
correspond to $|1\rangle$, $|2\rangle$ and $|3\rangle$, respectively.
$\delta_{m_{F}}$ ($m_{F}=-1,0,1$) are defined as $\delta_{-1}=0$,
$\delta_{0}=\Delta_{a,F,F^{\prime},\mathcal{D}}^{-1,-1}-\Delta_{b,F,F^{\prime},\mathcal{D}}^{-1,0}$
and
$\delta_{1}=\Delta_{a,F,F^{\prime},\mathcal{D}}^{-1,-1}-\Delta_{b,F,F^{\prime},\mathcal{D}}^{-1,0}+\Delta_{a,F,F^{\prime},\mathcal{D}}^{0,0}-\Delta_{b,F,F^{\prime},\mathcal{D}}^{0,1}$.
Since laser $a$ and $b$ are blue detuned from $D_{1}$ line while laser $c$ is
red detuned from $D_{1}$ line, laser $c$ only contributes to the AC stark
shift terms in diagonal matrix elements. If we rewrite eqn.(S14) as

$H_{R}=H_{11}|1\rangle\langle 1|+H_{22}|2\rangle\langle
2|+H_{33}|3\rangle\langle 3|+H_{12}|1\rangle\langle 2|+H_{23}|2\rangle\langle
3|+h.c.,$

and if $H_{22}=H_{11}$ and $|H_{33}-H_{22}|\gg H_{23}$, then $|1\rangle$ and
$|2\rangle$ are well isolated from the third state, $|3\rangle$, and the
effective two-level Hamiltonian takes the form $H_{e}=H_{11}|1\rangle\langle
1|+H_{22}|2\rangle\langle 2|+H_{12}|1\rangle\langle 2|+h.c.$ The above
condition can be satisfied by introducing proper laser frequencies in our
Raman process. Specifically, in our calculations, with $a_{0}=4$, one group of
possible experimental parameters are $B=10\textrm{G}$, $\lambda_{a}\approx
791.56\textrm{nm}$, $\lambda_{b}\approx 791.56\textrm{nm}$ and
$\lambda_{c}\approx 805.96\textrm{nm}$. The two-photon detuning between
$\ket{1}$ and $\ket{2}$ transition is around $400\textrm{kHz}$. All the time-
dependent and spatially dependent AC Stark shifts in the diagonal elements are
negligibly small due to the opposite AC stark shift caused by laser $c$. The
corresponding Floquet driving frequency is $\omega=2\pi\times 72\textrm{kHz}$
and $\Omega_{0}=2\pi\times 288\textrm{kHz}$. The Raman laser intensities are
$I_{a}\approx 24.2\textrm{W}/\textrm{cm}^{2}$, $I_{b}\approx
605\textrm{W}/\textrm{cm}^{2}$ and $I_{c}\approx
9.8\textrm{W}/\textrm{cm}^{2}$. Ignoring the constant term proportional to the
identity matrix, we get the desired Hamiltonian, $H_{d}$, in the main text
from the effective two-level Hamiltonian, $H_{e}$.

The next issue we care about is the scattering limited lifetime of the Bose
gas. Using the well-developed theory in optical scattering in cold atom
ensemblesGrimm  _et al._ (2000), the scattering rate can be written as

$\gamma_{sc}\approx\frac{\pi
c^{2}\Gamma^{2}}{2\hbar\omega_{0}^{3}}(\frac{2}{\Delta_{D_{2}}}+\frac{1}{\Delta_{D_{1}}})(I_{a}+I_{b}+I_{c}),$
(S15)

where $c$ is the speed of light, $\Gamma\approx 2\pi\times 6\textrm{MHz}$ is
the natural linewidth, and $\Delta_{D_{1}}$ and $\Delta_{D_{2}}$ are the one-
photon detunings with respect to the $D_{1}$ and $D_{2}$ lines, respectively.
By plugging in our parameters and using a Gaussian beam waist of $200\mu m$,
the life time $\tau_{sc}=1/\gamma_{sc}\approx 133\textrm{ms}$, which is
sufficient for experimental applications.