Cold atoms in twisted bilayer optical potentials

TL;DR

This paper proposes a cold atom setup in optical lattices to simulate twisted bilayer systems, enabling flexible control over coupling and exploration of phenomena like band narrowing and unconventional emitter interactions related to Moiré patterns.

Contribution

It introduces a novel cold atom platform for simulating twisted bilayer physics with adjustable parameters, extending the study of correlated phenomena beyond solid-state systems.

Findings

Potential to observe band narrowing at larger twist angles

Prediction of unconventional radiation patterns from emitter interactions

Flexible control over inter/intra-layer coupling

Abstract

The possibility of creating crystal bilayers twisted with respect to each other has led to the discovery of a wide range of novel electron correlated phenomena whose full understanding is still under debate. Here we propose and analyze a method to simulate twisted bilayers using cold atoms in state-dependent optical lattices. Our proposed setup can be used as an alternative platform to explore twisted bilayers which allows one to control the inter/intra-layer coupling in a more flexible way than in the solid-state realizations. We focus on square geometries but also point how it can be extended to simulate other lattices which show Dirac-like physics. This setup opens a path to observe similar physics, e.g., band narrowing, with larger twist angles, to rule out some of the mechanisms to explain the observed strongly correlated effects, as well as to study other phenomena difficult to realize with crystals. As an example of the latter we explore the quantum optical consequences of letting emitters interact with twisted bilayer reservoirs, and predict the appearance of unconventional radiation patterns and emitter interactions following the emergent Moir\'e geometry.

Equations

$$|{\hat{e}}_{{\sigma }_{\pm }}\cdot \mathbf{E}(\mathbf{R}){|}^2$$

$$|{\hat{e}}_z\cdot \mathbf{E}(\mathbf{R}){|}^2$$

$$I_{{\sigma }_{\pm }}(\mathbf{R})=$$

$$I_z(\mathbf{R};\theta )=$$

$$+$$

$$V_{\alpha }(\mathbf{R})\approx -\frac{|{\tilde{\Omega }}_{\alpha }{|}^2}{{\Delta }_{\alpha }}{I}_{\alpha }(\mathbf{R}).$$

$$\frac{|{\Omega }_a{|}^2}{{\Delta }_a}=\frac{|{\Omega }_b{|}^2}{{\Delta }_b}\equiv V_D.$$

$${\epsilon }_{2ph}\sim {\left(\frac{V_D}{\delta }\right)}^2\ll 1,$$

$${\Gamma }_{\ast }\approx \frac{V_D}{\Delta }{\Gamma }_g,$$

$$V_a(\mathbf{R})\approx -\frac{|{\Omega }_a{|}^2}{{\Delta }_a}{I}_{\pi }(\mathbf{R})-\frac{|{\Omega }_b{|}^2}{{\Delta }_a+{\delta }_g}{I}_{{\sigma }_{+}}(\mathbf{R}),$$

$$V_b(\mathbf{R})\approx -\frac{|{\Omega }_b{|}^2}{{\Delta }_b}{I}_{{\sigma }_{-}}(\mathbf{R})-\frac{|{\Omega }_a{|}^2}{{\Delta }_b+{\delta }_g}{I}_{\pi }(\mathbf{R}),$$

$$V_a(\mathbf{R})\approx -V_D\left(I_{\pi }(\mathbf{R})+\frac{{\Delta }_b}{{\Delta }_a+{\delta }_g}{I}_{{\sigma }_{+}}(\mathbf{R})\right),$$

$$V_b(\mathbf{R})\approx -V_D\left(I_{{\sigma }_{-}}(\mathbf{R})+\frac{{\Delta }_a}{{\Delta }_b+{\delta }_g}{I}_{\pi }(\mathbf{R})\right),$$

$$V_{a,b}(\mathbf{R})\propto -\frac{|{\Omega }_p{|}^2}{{\Delta }_p^2}\frac{|{\tilde{\Omega }}_{a,b}{|}^2}{{\Delta }_{a,b}}{I}_{{\sigma }_{-},\pi }(\mathbf{R}),$$

$$V_D=\frac{|{\Omega }_p{|}^2}{{\Delta }_p^2}\frac{|{\tilde{\Omega }}_{a,b}{|}^2}{{\Delta }_{a,b}}.$$

$$H_B$$

$$+\sum _{\alpha }\frac{4\pi {\hslash }^2{a}_{s,\alpha }}{2m}\sum _{\alpha }\int d\mathbf{R}{\hat{\psi }}_{\alpha }^{†}(\mathbf{R}){\hat{\psi }}_{\alpha }^{†}(\mathbf{R}){\hat{\psi }}_{\alpha }(\mathbf{R}){\hat{\psi }}_{\alpha }(\mathbf{R})$$

$$+{\Omega }_{\perp }\int d\mathbf{R}{e}^{i{\mathbf{k}}_L\cdot \mathbf{R}-i{\omega }_{L}t}{\hat{\psi }}_a(\mathbf{R}){\hat{\psi }}_b^{†}(\mathbf{R})+H.c.,$$

$$w_{\alpha }(\mathbf{R}-{\mathbf{R}}_{\alpha ,\mathbf{n}})=\frac{1}{L_{\alpha }\sqrt{\pi }}{e}^{-\frac{x^2+y^2}{2L_{\alpha }^2}}=\frac{1}{L_{\alpha }\sqrt{\pi }}{e}^{-\frac{r^2}{2L_0^2}}.$$

$$H_{\alpha }=-\sum _{\mathbf{n},\mathbf{m}}{J}_{\alpha }(\mathbf{n}-\mathbf{m}){\alpha }_{\mathbf{n}}^{†}{\alpha }_{\mathbf{m}},$$

$$J_{\alpha }(\mathbf{n}-\mathbf{m})$$

$$w_{\alpha }(\mathbf{R}-{\mathbf{R}}_{\alpha ,\mathbf{m}})\propto E_R{e}^{-|{\mathbf{R}}_{\alpha ,\mathbf{n}}-{\mathbf{R}}_{\alpha ,\mathbf{m}}{|}^2/(4L_{\alpha }^2)},$$

$$H_U=\sum _{\alpha ,\mathbf{n}}\frac{U_{\alpha }}{2}{\alpha }_{\mathbf{n}}^{†}{\alpha }_{\mathbf{n}}({\alpha }_{\mathbf{n}}^{†}{\alpha }_{\mathbf{n}}-1),$$

$$H_{ab}=-\sum _{\mathbf{n},\mathbf{m}}\left(J_{\perp }({\mathbf{R}}_{a,\mathbf{n}}-{\mathbf{R}}_{b,\mathbf{m}})a_{\mathbf{n}}^{†}{b}_{\mathbf{m}}+H.c.\right),$$

$$J_{\perp }(\mathbf{n}-\mathbf{m})\propto {\Omega }_{\perp }{e}^{-|{\mathbf{R}}_{a,\mathbf{n}}-{\mathbf{R}}_{b,\mathbf{m}}{|}^2/(4L_0^2)},$$

$$H_B=\sum _{\alpha }{H}_{\alpha }+H_{U,\alpha }+H_{ab}.$$

$$\theta (m,n)=\arccos \left(\frac{2mn}{m^2+n^2}\right),n,m\in \mathbb{Z},$$

$$D(\omega )\propto \sum _{j=1}^{N_{\mathrm{M}}}\sum _{\mathbf{k}}\mathcal{N}(\omega -\Delta \omega <{\omega }_j(\mathbf{k})<\omega ),$$

$${\theta }_{\mathrm{hc}}(m,n)=\arccos \left(\frac{n^2+4nm+m^2}{2(m^2+n^2+nm)}\right),n,m\in \mathbb{Z},$$

$$H_c={\omega }_c\sum _i{c}_i^{†}{c}_i+\frac{U_c}{2}\sum _i{c}_i^{†}{c}_i(c_i^{†}{c}_i-1),$$

$$H_{\mathrm{int}}=g\sum _i\left(a_{{\mathbf{n}}_i}^{†}{c}_i+H.c\right),$$

$$|\Psi (t)⟩=\left(C_e(t)c^{†}+\sum _{\alpha ,\mathbf{R}}{C}_{\alpha ,\mathbf{n}}(t){\alpha }_{\mathbf{n}}^{†}\right)|\mathrm{vac}⟩,$$

$$H_{\mathrm{spin}}=\sum _{i,j}{J}_{i,j}{c}_i^{†}{c}_j,$$

Full text

I Introduction

The recent observation of unconventional superconducting Cao et al. (2018a) and correlated insulating behaviour Cao et al. (2018b) in twisted bilayer graphene has brought the study of twisted van-der Waal materials, that is, rotated weakly coupled layers, at the forefront of condensed-matter research. One of the initial motivations of the twisting was to shift graphene van-Hove singularities closer to half-filling to induce superconductivity by doping Nandkishore et al. (2012); Li et al. (2010); McChesney et al. (2010). However, the interest in those systems burst with the observation that the Fermi velocity, $v_{F}$, around the Dirac point renormalizes strongly for small rotations Lopes dos Santos et al. (2007); Shallcross et al. (2010); Trambly de Laissardière et al. (2010); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011), especially at the so-called magic angles Trambly de Laissardière et al. (2010); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011). At these angles, the interplay between the interlayer hopping, $\hbar J_{\perp}$ and the momentum mismatch energy, $\hbar v_{F}k_{\theta}$, leads to bands whose bandwidth is comparable to the effective Coulomb interactions and with large density of states, making them an ideal candidate to observe strongly correlated phenomena. This potential has been evidenced by the renowned transport experiments Cao et al. (2018a, b) as well as the observation of other exotic effects such as super-planckian dissipation Cao et al. (2019) or ferromagnetism Sharpe et al. (2019). Beyond these observations, the experimental efforts have been devoted to finding ways of tuning these systems, e.g., using pressure Yankowitz et al. (2019), electrical gating Chen et al. (2018), dynamical angle control Ribeiro-Palau et al. (2018), magnetic fields Kim et al. (2017); Sharpe et al. (2019), or temperature Cao et al. (2019), as well as alternative forms to probe the system beyond electron transport Kerelsky et al. (2018); Li et al. (2010); Kim et al. (2017); Choi et al. (2019); Yoo et al. (2018); Ni et al. (2015); Sunku et al. (2018).

Since the inter/intra-layer hopping ratio is small in these systems, the magic angles occur for very small rotations, $\theta\lesssim 1^{\circ}$. This has several consequences: first, it makes the structures more sensitive to structural relaxation Yoo et al. (2018); Yankowitz et al. (2019), to disorder, or the effect of the substrate Sharpe et al. (2019); second, the Moiré primitive cells are enormous ($\sim 10^{4}$ atoms per unit cell), making ab-initio calculations challenging. This is why one of the main theoretical research directions consists in finding accurate descriptions for inter/intra-layer hoppings in tight-binding Shallcross et al. (2010); Trambly de Laissardière et al. (2010); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011); Moon and Koshino (2012); Fang et al. (2015); Fang and Kaxiras (2016) or continuum descriptions Lopes dos Santos et al. (2007); San-Jose et al. (2012), including lattice relaxation effects Uchida et al. (2014); Nam and Koshino (2017); Zhang and Tadmor (2018); Wen et al. (2018); Carr et al. (2017, 2018a, 2018b); Koshino et al. (2018); Guinea and Walet (2018a); Walet and Guinea (2018), as well as searching for minimal models that capture the physics of the so-called active bands Kang and Vafek (2018a); Po et al. (2018a); Zou et al. (2018); Else et al. (2018); Song et al. (2018); Yuan and Fu (2018); de la Peña et al. (2017); Classen et al. (2019). Despite all these efforts, the origin of the emergent superconductivity Lian et al. (2018); Wu et al. (2018); Peltonen et al. (2018); Isobe et al. (2018); You and Vishwanath (2018); Lin and Nandkishore (2019); González and Stauber (2019); Po et al. (2018b); Venderbos and Fernandes (2018); Xu and Balents (2018); Dodaro et al. (2018); Kennes et al. (2018); Guinea and Walet (2018b); Xu et al. (2018); Kang and Vafek (2018b), the correlated insulating behaviour Po et al. (2018b); Pizarro et al. (2018); Padhi et al. (2018); Padhi and Phillips (2018), or even the magic angles Tarnopolsky et al. (2019), is still under debate. For example, several mechanisms have been proposed to explain the experimentally observed superconductivity such as electron-phonon interaction Lian et al. (2018); Wu et al. (2018); Peltonen et al. (2018), assisted correlated hopping Guinea and Walet (2018b); Xu et al. (2018), or the effective electron attraction induced by the presence of van-Hove singularities and nested Fermi-surfaces Isobe et al. (2018); You and Vishwanath (2018); Lin and Nandkishore (2019); González and Stauber (2019). Thus, it seems very timely to find other platforms where to study these effects and enlarge our understanding of this exciting system.

In this work, we propose a way to engineer state-dependent optical potentials for cold atoms such that their dynamics mimic that of twisted bilayers. Cold atoms in optical lattices stand nowadays as one of the most mature implementations to simulate condensed-matter problems Lewenstein et al. (2007); Bloch et al. (2012), including the physics of graphene-like materials Tarruell et al. (2012) and even their modification when placed on a substrate Grass et al. . The key ingredient of our proposal is the creation of two orthogonally-polarized standing waves with an angle between them such that when judiciously coupled to an atom they generate twisted optical potentials which trap independently two long-lived atomic states. The localized atomic excitations in each potential can hop within the lattice through the hybridization of their atomic wavefunctions simulating the intralayer hoppings. On top of that, we add an additional global field to connect the two internal atomic states, e.g., via a two-photon Raman transition through an auxiliary level, to obtain the interlayer hopping. Thus, the inter/intra-layer couplings can be controlled independently by optical means. This control opens up the possibility of observing similar phenomena than in solid-state implementations, such as the emergence of narrow bands or the movement (or appearance) of van-Hove singularities, but for larger rotation angles. In this manuscript we focus on square Moiré geometries, but it can potentially be extended to brick-wall Tarruell et al. (2012) or hexagonal lattices Soltan-Panahi et al. (2011), which would allow the exploration of Dirac-like physics. Finally, we also show how our proposed setup can be adapted to explore physics beyond the one observed with crystals such as the one of emitters interacting with twisted bilayer reservoirs. For that configuration, we predict the appearance of strongly non-Markovian dynamics controlled by the twisting angle, and the generation of emission patterns and interactions following the emergent Moiré geometry.

The manuscript is structured as follows: first, in Section II, we discuss how to generate twisted optical potentials for cold atoms using Alkali-Earth atoms Ludlow et al. (2015). Then, in Section III, we derive the effective atom dynamics in these twisted optical lattices (in III.1), and calculate its associated band structure and density of states (in III.2). Afterwards, in Section IV, we explain how to simulate quantum optical effects in these setups by using an additional atomic state, and calculate the spontaneous emission dynamics and radiation patterns of an emitter coupled to such twisted bilayer reservoirs. Finally, we summarize our findings in Section V, and provide an outlook of future directions of our work.

II Twisted optical potentials

In this section we describe in several steps how to obtain twisted independent optical potentials for cold atoms to simulate the physics of twisted crystal bilayers. First, in II.1 we explain how to generate a set of orthogonally polarized standing waves with a square geometry and twisted with respect to each other. Second, in II.2 we describe the basic atomic configuration required to harness these standing waves to generate independent optical potentials for two atomic states. Finally, in II.3 we discuss several specific realizations of that configuration based on Alkali-Earth atomic level structure Ludlow et al. (2015).

II.1 Twisted orthogonally polarized standing waves

The laser configuration is sketched in Fig. 1(a): one initially starts with two retro-reflected lasers with in-plane polarization in orthogonal directions ($x$ and $y$) which create a two-dimensional square-like potential $\propto\sin^{2}(kx)+\sin^{2}(ky)$ (in red). To create additional standing waves tilted an angle $\theta$ with respect to each of the beams, one can tilt one of the mirrors by an angle $\theta/2$, add a half-wave plate that changes the in-plane polarization $\hat{\epsilon}_{x}$ into $\hat{\epsilon}_{z}$ (and viceversa), and include an additional mirror that reflects the light back to create a standing wave. Denoting the light-paths as $1$-$4$ like in Fig. 1(a), the momenta and polarization of the lasers along this closed path read $\mathbf{k}_{1,2}=k(0,\pm 1,0)$, $\mathbf{k}_{3,4}=k(\pm\sin(\theta),\mp\cos(\theta),0)$, $\hat{\varepsilon}_{1,2}=\pm\hat{e}_{x}$, and $\hat{\varepsilon}_{3,4}=\pm\hat{e}_{z}$, where $\lambda=2\pi/k$ is the wavelength of the laser. Since the polarizations in the paths $1$-$2$ and $3$-$4$ are orthogonal, the two standing waves will not interfere even if they have the same frequency. In the orthogonal path, denoted by $5$-$8$ in the Fig. 1(a), we have something similar with $\mathbf{k}_{5,6}=k^{\prime}(\pm 1,0,0)$, $\mathbf{k}_{7/8}=k^{\prime}(\mp\cos(\theta),\mp\sin(\theta),0)$, $\hat{\varepsilon}_{5,6}=\pm\hat{e}_{y}$, and $\hat{\varepsilon}_{7,8}=\pm\hat{e}_{z}$. To avoid interference between the paths $3,4,7,8$ (and $1,2,5,6$) that have the same polarization, we use the standard procedure Grimm et al. (2000) of sending slightly detuned lasers with $\lambda\approx\lambda^{\prime}$ (so that the periodicities of the potentials are the same) but whose frequency $\omega_{L}(^{\prime})=2\pi/\lambda(^{\prime})$ are sufficiently far such that their cross-talk average out in the atomic timescales (of the order tens of kHz). Taking into account all these considerations, the time averaged intensity felt by a linear/circularly polarized atomic transition read:

[TABLE]

where $|E_{0}|^{2}$ is the intensity of the lasers generating the standing waves, and the functions $I_{\alpha}$ are the (adimensional) electric field profile intensity that read:

[TABLE]

where we have defined $z$ as the quantization axis, and written the in-plane electric field as a combination of circular polarized light: $\hat{e}_{x}=\frac{1}{\sqrt{2}}\left(\hat{e}_{\sigma_{-}}+\hat{e}_{\sigma_{+}}\right)$ and $\hat{e}_{y}=\frac{i}{\sqrt{2}}\left(\hat{e}_{\sigma_{-}}-\hat{e}_{\sigma_{+}}\right)$. Notice the factor $1/2$ which appears in the circular component intensity in Eq. 1 because the in-plane intensity $\hat{e}_{x/y}$ divides between the two circular polarizations. The vertical confinement can be obtained adding another pair of retro-reflected lasers in the vertical direction, $k_{9/10}=k_{z}(0,0,\pm 1)$, that will generate a standing wave $\propto\sin^{2}(k_{z}z)$. Since we do not want the atoms to hop in that direction we can choose a different frequency from the one generating the twisted standing waves such that no matter the polarization of $\hat{\varepsilon}_{9/10}$ it does not interfere with them.

II.2 Basic atomic scheme

Let us now explain how by adequately choosing the atomic level structure, this set of standing waves forms a different optical potential for two different long-lived atomic states labeled as $a,b$. We consider here the configuration depicted in Fig. 1(b) where these two states connect to an additional one $g$ (with decay rate $\Gamma_{g}$) through an orthogonal polarization of light, e.g., $\sigma_{+},\pi$ (we show in the next section how). Consequently, the effective Rabi amplitude of these optical transitions, $\Omega_{a,b}(\mathbf{R})$, inherit the spatial dependence of the intensity profiles of Eqs. 1-2. This can be explicitly shown by writing $\Omega_{a,b}(\mathbf{R})=\tilde{\Omega}_{a,b}I_{\sigma_{+},\pi}(\mathbf{R})$, where in $\tilde{\Omega}_{\alpha}$ we include the overall transition strength including the contribution of the Clebsch-Gordan coefficients as well as the $1/2$ factor of Eq. 1. Next, we assume that these long-lived states are split in energy by $\delta=\omega_{b}-\omega_{a}$, e.g., by a magnetic field, such that their optical transitions to $g$ have different detunings $\Delta_{\alpha}=\omega_{g}-\omega_{L}-\omega_{\alpha}$, for $\alpha=a,b$.

Then, in the limit where $|\tilde{\Omega}_{\alpha}|\ll\Delta_{\alpha}$ and $|\tilde{\Omega}_{a}\tilde{\Omega}_{b}|\ll\Delta_{\alpha}\delta$, the dressing of the states $\alpha$ with $g$ leads to the following state-dependent optical potential Grimm et al. (2000):

[TABLE]

The role of the energy splitting $\delta$ is two-fold:

• •

It allows to compensate for the different coupling amplitudes of the optical transitions generating the potentials. For example, if we want to obtain the same potential depth $V_{D}$ for the $a/b$ levels we must impose:

[TABLE]

As we will see afterwards, this guarantees obtaining the same hopping rate in each optical potential, which is what simulates the intralayer hopping $J$ for both atomic excitations.

• •

It suppresses the unwanted two-photon transitions between the states. The two-photon amplitude between two states scales with $\Omega_{\mathrm{2ph}}\approx\tilde{\Omega}_{a}\tilde{\Omega}_{b}^{*}/\Delta_{a}$, which is of the order of the potential depth ($V_{D}$), whereas its effective detuning is $\delta$. Thus, for the two-photon transition probability ($\varepsilon_{\mathrm{2ph}}$) to be small we require:

[TABLE]

Apart from these trapping lasers, to complete the proposal we need a mechanism that couples the $a,b$ states directly. Depending on the particular choice for the $a,b$ levels (see next section) this can be done in several ways as depicted in green in Fig. 1(b), either with a global microwave (optical) field which couples directly these states, or with a two-photon Raman process through an auxiliary excited state. Any of these mechanisms provides an effective coupling strength, $\Omega_{\perp}$, between the $a/b$ states that can be controlled independently from the lasers generating the twisted optical potentials. Thus, we have a knob to change the interlayer coupling between the states, $J_{\perp}$, which is independent from the one generating the intralayer hopping, $J$, as schematically depicted in Fig. 1(c).

Before giving a particular realization of this atomic level structure, let us discuss more conditions values for the proposal:

• •

The trapping depths must be big enough such that the atoms are trapped within their motional ground state. The rule-of-the thumb condition is that it should be bigger than the recoil energies Bloch et al. (2008) of the lattice $V_{D}>E_{R}=\frac{\hbar^{2}k^{2}}{2m}$, of the order of $1-10$ kHz for standard cold atoms experiments. Thus, $V_{D}$ should be $\sim 10$-$100$ kHZ.

• •

If $g$ is an unstable state with a decay rate $\Gamma_{g}$, the $a,b$ levels acquire a finite decoherence rate which scales as:

[TABLE]

and that will set a limitation on the coherence timescales of our simulation. In practice, $\Gamma_{*}$ should be much smaller that the effective atom dynamics that we want to simulate, e.g., the band-width of the emergent Moiré bands. This suggests using states for which $\Gamma_{g}$ is as small as possible, as we discuss in the next section.

II.3 Specific realizations using Alkali-Earth atoms

The ideal level structure considered in Fig. 1(b) can be easily engineered in an optical transition between a ground state with total electron angular momentum $J=1$ and an excited state with $J^{\prime}=0$. These optical transitions appear in rare-earth atoms like Samarium Nottelmann et al. (1993). Unfortunately, these atoms are difficult to cool down to their ground state. Other ways of generating such state-dependent optical traps consists in using ground and excited hyperfine levels in Alkali atoms Rubio-Abadal et al. ; Krinner et al. (2018) for the $a/b$ and $g$ states, respectively. However, these potentials typically suffer from photon loss decoherence due to the small hyperfine energy splitting so that the conditions required above would be hardly fulfilled.

For the above reasons we just provide several possible realizations of the scheme of Fig. 1(b) using Alkali-Earth atoms Ludlow et al. (2015), such as Calcium (Ca), Strontium (Sr), or Ytterbium (Yb). All these atoms share an electronic configuration with two electrons in an outer $s$-shell, generating a level structure (see Fig. 2) with several remarkable features:

• •

The two electrons can align in singlets ($S=0$) and triplets ($S=1$). Strong transitions (depicted in red in Fig. 2) occur between singlet-triplet states, e.g., ${}^{1}P_{1}$-${}^{1}S_{0}$ with decay rates $\Gamma/(2\pi)\sim 30$ MHz for Sr or Yb Ludlow et al. (2015). These transitions are well separated spectrally and can be used for trapping, state detection, or cooling. In fact, for all these atoms there exist already refined laser cooling techniques which can bring them to their ground state Ludlow et al. (2015); Takasu et al. (2003); Kraft et al. (2009); Stellmer et al. (2009); Snigirev et al. (2017); Riegger et al. (2018).

• •

Weaker decay rates occur for spin-forbidden transitions. For example, the transition ${}^{3}P_{1}$-${}^{1}S_{0}$ (in blue) is a simply spin forbidden transition leading to lifetimes of the order of $300$ Hz for Ca Kraft et al. (2009), $7.5$ kHz for Sr Ludlow et al. (2015), and $180$ kHz for Yb Ludlow et al. (2015). The clock transition ${}^{3}P_{0}$-${}^{1}S_{0}$ is doubly forbidden (in green) leading to narrow linewidths of $\sim 1$-$10$ mHz for Ytterbium and Strontium, respectively Ludlow et al. (2015), whereas the state ${}^{3}P_{2}$ has even longer predicted lifetime.

• •

All these atoms have bosonic (40Ca, 88Sr, 174Yb) and fermionic (43Ca, 87Sr, 171,173Yb) stable isotopes. The fermionic ones have non-zero nuclear spin $I\neq 0$ depending on their atomic mass, e.g., $I=1/2$ and $5/2$ for 171Yb and 173Yb, respectively. This nuclear spin can play an important role in the context of quantum simulation and information because it leads to degenerate multiplet hyperfine levels in the clock states Daley et al. (2008); Gorshkov et al. (2010). Moreover, using a magnetic field these levels can be split with differential Zeeman shifts for the ${}^{1}S_{0}$ and ${}^{3}P_{0}$ of the order of $1$ MHz/T Ludlow et al. (2015), which allows to address them independently Daley et al. (2008).

After reviewing the main spectral properties of these atoms, let us propose three independent ways of generating the twisted optical potentials as summarized in Fig. 3: the first two (in panels a-b) are inspired in the scheme discussed in the previous Section and relies on a judicious Zeeman splitting of hyperfine or fine structure levels using magnetic fields. On the other hand, the third one (panel c) does not rely on polarization or magnetic fields but on sending one of the lasers with different frequency and a certain out-of-plane angle with respect to the plane of the layer, as proposed in the context of quantum information in Ref. Daley et al. (2008).

II.3.1 Using hyperfine levels of the clock state ${}^{3}P_{0}$

One possible realization consists in using the hyperfine levels of the state ${}^{3}P_{0}$ of a fermionic isotope as our $a/b$ levels. In particular, we can take $\left|a\right\rangle=\left|{}^{3}P_{0},I,m_{I}=-I\right\rangle$ and $\left|{}^{3}P_{0},I,m_{I}=-I+1\right\rangle$. Differently from our ideal scheme, the off-resonant lasers $\Omega_{a/b}$ connect now these states to two different hyperfine levels of the ground state, namely, $\left|g_{1}\right\rangle=\left|{}^{1}S_{0},I,m_{I}=-I\right\rangle$ and $\left|g_{2}\right\rangle=\left|{}^{1}S_{0},I,m_{I}=-I+1\right\rangle$ (see Fig. 3(a)). Thus, both the $a,b$ states feel the linear and circularly polarized intensities, but with different detunings such that they yield the following optical potentials:

[TABLE]

where we have defined $\Delta_{a/b}=\omega_{g_{1}}-\omega_{a/b}-\omega_{L}$, and $\delta_{g}=\omega_{g_{2}}-\omega_{g_{1}}$ 111Notice that here for simplicity we assume that the Clebsch Gordan coefficients for the transitions $a\rightarrow g_{2}$ and $b\rightarrow g_{1}$ (as well as for $a\rightarrow g_{1}$ and $b\rightarrow a_{1}$) are the same. Even if they are different this will become irrelevant in the regime where one can neglect the non-ideal contributions.. Imposing Eq. 6 such that the potential depth $V_{D}$ for both potentials is the same, we can rewrite the potentials:

[TABLE]

Notice, that $V_{a,b}(\mathbf{R})$ contains now two contributions: the desired one (left hand side of the bracket) plus a correction (right-hand side) that might spoil the twisted bilayer behaviour. Thus, in order for the right-hand term of these equations to be negligible, we require that $\frac{\Delta_{\alpha}}{\Delta_{\beta}+\delta_{g}}\ll 1$. Since $\delta_{g}\sim 2$ MHz for $B\sim 1$ T, this restricts $\Delta_{\alpha}/(2\pi)\sim 0.2$ MHz. The advantage here is that since the dressing is done with an stable state (${}^{1}S_{0}$), ones does not have to worry about canceling $\Gamma_{g}$. Thus, $\Omega_{\alpha}$ can be a significant fraction of $\Delta_{\alpha}$, e.g., $\Omega_{\alpha}\approx 0.25\Delta_{\alpha}$, yielding $V_{D}/(2\pi)\sim 10$ kHz.

II.3.2 Using fine structure levels of the ${}^{3}P_{2}$ through a two-photon transition

One of the challenges of the previous realization is that the use of hyperfine states limits the energy splitting $\delta_{g}$ required to cancel the non-desirable terms in Eqs. 11-12. An alternative is to use fine structure levels as our target levels, e.g., $\left|a,b\right\rangle=\left|{}^{3}P_{2},m_{J}=-1,0\right\rangle$, whose Zeeman splitting can be several orders of magnitude larger for similar magnetic fields ($\sim 10$ GHz for $\sim T$ fields). The problem here is that these states do not connect directly to the ground state ${}^{1}S_{0}$. Thus, we have to do it through a two-photon transition via the states $\left|g_{1,2}\right\rangle=\left|{}^{3}P_{1},m_{J}=-1,0\right\rangle$ as depicted in Fig. 3(b): the twisted standing wave lasers of section II.1 couple the ground state ${}^{1}S_{0}$ to the intermediate $g_{1,2}$ states through the spatially dependent $\Omega_{a,b}(\mathbf{R})$ lasers. This generates a state-dependent Stark-shift on the $g_{1,2}$ levels, which translates to the ${}^{3}P_{2}$ through a global laser field with linear polarization, Rabi amplitude $\Omega_{p}$, and frequency $\omega_{L,p}$. This laser couples the ${}^{3}P_{1}$ and ${}^{3}P_{2}$ states off resonantly with an effective detuning which in this case is approximately equal for both transitions $\Delta_{p}=\omega_{g_{1}}-\omega_{a}-\omega_{L,p}\approx\omega_{g_{2}}-\omega_{b}-\omega_{L,p}$ 222This is a good approximation because the ${}^{3}P_{J}$ states have $S=L=1$ such that the Landé factor $g_{J}\approx 3/2$ independent on the total angular momentum, whereas the global two-photon transition detuning for each path reads $\Delta_{a,b}=\omega_{s}-\omega_{L}-\omega_{L.p}-\omega_{a,b}$. This two photon transition can be shown to yield a spatially dependent potential for the $a,b$ states which reads:

[TABLE]

when $|\tilde{\Omega}_{a,b}|^{2}\ll\Delta_{a,b}\Delta_{p}$. Notice that we require that the laser $\Omega_{p}$ does not introduce any additional spatial dependence which can be achieved with a running-wave laser configuration. Thus, fixing $\Delta_{a,b}$ like in Eq. 6, the overall trapping depth scales in this case as:

[TABLE]

Since the intermediate states $g_{1,2}$ have a finite lifetime $\Gamma_{g}^{-1}$, the $a/b$ states acquire a decay rate which scales with $\Gamma_{*}\sim\frac{|\Omega_{p}|^{2}}{\Delta_{p}^{2}}\Gamma_{g}$. Thus, to obtain coherence times of the order of seconds, we must impose $|\Omega_{p}/\Delta_{p}|\sim 10^{-1},10^{-2}$ for Ca or Sr atoms whose $\Gamma_{g}/(2\pi)=0.3,7.5$ kHz, respectively. Since $\Delta_{a,b}$ can be $\sim 10$ GHz in this case, and $|\tilde{\Omega}_{a,b}|\ll\Delta_{a,b}$, this leads to $V_{D}/(2\pi)\sim 10$-$1000$ KHz, depending on the ratio $\Omega_{p}/\Delta_{p}$ required to cancel $\Gamma_{g}$.

The main technical difficulty of this proposal is to drive directly the transition ${}^{3}P_{1}$-${}^{3}P_{2}$ which requires a laser in the $10$ THz range. An alternative way might be driving indirectly with a two-photon Raman transition with optical lasers through the ${}^{3}S_{1}$ state. In that case, one must use a detuning large enough to cancel the decay rate of the ${}^{3}S_{1}$ state.

II.3.3 Turnout wavelengths with angled beams

Finally, let us mention one last possibility which does not rely on the generation of the orthogonally polarized standing waves explained in section II.1. It is inspired on the existence of the so-called turnout wavelengths ($\lambda_{1,2}$) of Alkaline-Earth atoms Daley et al. (2008) in which the polarizabilities of the two-long-lived states $\left|{}^{1}S_{0}\right\rangle=\left|a\right\rangle$ and $\left|{}^{3}P_{0}\right\rangle=\left|b\right\rangle$ vanish independently, that is, $\alpha_{a,b}(\lambda_{1,2})=0$ and $\alpha_{a,b}(\lambda_{2,1})\neq 0$. This allows one to create two independent potentials for the $a/b$ states using two lasers with these wavelengths. The only problem for our purposes is that these wavelengths differ significantly, e.g., $\lambda_{1,2}=689,627$ nm for Strontium Ludlow et al. (2015), such that if both lasers are sent in-plane as in Fig. 1(a) their potentials will have different periodicities. However, this can be compensated by sending the shorter wavelength laser with an out-of-plane angle $\alpha$ Daley et al. (2008), as schematically depicted in Fig. 3(c), such that the in-plane momenta is reduced by a factor $\cos(\alpha)=\lambda_{2}/\lambda_{1}$. For the vertical trapping, one can send an additional pair of retro-reflected lasers in the $z$ direction with the so-called magic-wavelength ($\lambda_{m}$) that leads to the same polarizability for both states $\alpha_{{}^{1}S_{0}}(\lambda_{m})=\alpha_{{}^{3}P_{0}}(\lambda_{m})$ Ludlow et al. (2015), and with a large intensity such that it provides a much deeper potential than the one in the XY plane. The coupling between the $a,b$ levels can be obtained by resonantly driving the clock transition or through an intermediate state, e.g., ${}^{3}S_{1}$, to avoid the small dipole moment of the clock transition.

Let us finally provide some estimates of the coherence times and trapping depths of this configuration. Since the metastable states $a,b$ are dressed by the excited states ${}^{1}P_{1}$ and ${}^{3}S_{1}$, respectively, with lifetimes $\Gamma_{e}/(2\pi)\sim 10$-$100$ MHz, for the coherence times of $\Gamma_{a/b}/(2\pi)\lesssim 1$ Hz the dressing with the excited states have to be kept small $\tilde{|\Omega}_{a,b}/\Delta_{a,b}|\sim 10^{-3}$-$10^{-4}$. Thus, for $\Delta_{a,b}$’s of the order of $100$ GHz, one can obtain $V_{D}/(2\pi)\sim 1$-$100$ kHz.

III Emergent band structure and density of states

Here, in Section III.1 we write the Hamiltonian describing the dynamics of the atoms moving in the twisted optical potentials $V_{a/b}(\mathbf{R})$. Then, in Section III.2 we calculate the associated density of states of the emergent band structure. For completeness, in Section III.3 we discuss the possibility to go to other geometries to obtain Dirac-like physics.

III.1 Effective hopping dynamics in the twisted optical potentials: intra- and interlayer couplings

The dynamics of atoms subject to the state-dependent optical potential $V_{a,b}(\mathbf{R})$ plus the coupling between the internal atomic states $a,b$ through $\Omega_{\perp}$ can be described by the following Hamiltonian Jaksch et al. (1998):

[TABLE]

where $\hat{\psi}_{\alpha}(\mathbf{R})$ is the atom field operator of the state $\alpha$, that will be fermionic or bosonic depending on the isotope we are considering. Let us now explain in detail the different parts of the Hamiltonian:

• •

Hopping dynamics (Eq. 15). This term corresponds to the kinetic energy of the atoms plus the optical potential in each internal atomic state. When the trapping depth is big enough ($V_{D}>E_{R}$), the atoms localize within the minima of the potential, that we denote by $\mathbf{R}_{\alpha,\mathbf{n}}$ with $\mathbf{n}\in\mathbb{Z}^{2}$ and $\alpha=a,b$, and which can be described as parabolic isotropic traps 333To obtain accurate descriptions beyond this limit one needs to calculate the exact band-structure using Bloch theorem and obtain the maximally localized Wannier functions Marzari and Vanderbilt (1997). In any case, the Hamiltonian parameters, e.g., tunneling, can also be controlled through the optical potential features.. Notice that here we are already considering that the trapping in the $z$ direction is much deeper than in the $x$-$y$ plane such that the atoms can only hop in these directions. Within that regime, we can expand the atomic field operators in a separable Wannier basis in the $x$-$y$ direction, and keep only the lowest ground state levels such that $\hat{\psi}_{\alpha}(\mathbf{R})=\sum_{\mathbf{n}}w_{\alpha}(\mathbf{R}-\mathbf{R}_{\alpha,\mathbf{n}})\alpha_{\mathbf{n}}$, with

[TABLE]

with $L_{\alpha}=\sqrt{\frac{\hbar}{m\omega_{t,\alpha}}}$ being the ground state wavefunction size, $\hbar\omega_{t,\alpha}$ the trapping frequency. The $\alpha_{\mathbf{n}}$ ($\alpha^{\dagger}_{\mathbf{n}}$) represent the annihilation (creation) of an atomic excitations at site $\mathbf{R}_{\alpha,\mathbf{n}}$. Projecting this part of the Hamiltonian in the Wannier basis we arrive to the following form:

[TABLE]

where the tunnelings between the different sites are given by the following overlap integrals Bloch et al. (2008):

[TABLE]

while for $\mathbf{n}=\mathbf{m}$ is equal to $\omega_{t,\alpha}$. When the trapping depth of the potentials $V_{a,b}(\mathbf{R})$ is the same, so are the trapping frequencies ($\omega_{t,\alpha}\equiv\omega_{t}$) and ground state sizes ($L_{\alpha}\equiv L_{0}$), which ultimately determine the strength and range of the hoppings, respectively. Thus, the hopping within each simulated layer will be the same $J_{\alpha}(\mathbf{n})\equiv J(\mathbf{n})$.

• •

On-site interactions (Eq. 16). In the low energy limit these interactions are well approximated by $s$-wave scattering whose strength is determined by so-called scattering length $a_{s,\alpha}$ Bloch et al. (2008). In the Wannier basis, this part of the Hamiltonian reads:

[TABLE]

with $U_{\alpha}=\frac{4\pi\hbar^{2}a_{s,\alpha}}{m}\int d\mathbf{R}|w_{\alpha}(\mathbf{R})|^{4}$. One of the advantage of our setup is that the scattering length, $a_{s,\alpha}$ (and consequently $U_{\alpha}$) can be controlled through Feshbach resonances Bloch et al. (2008). Thus, one can tune from having purely non-interacting models with $U_{\alpha}\approx 0$ to strongly interacting ones of repulsive or attractive character. Furthermore, with the proper choice of the atomic levels one can also obtain crossed on-site interactions between the different internal atomic states $\propto U_{ab}a^{\dagger}_{\mathbf{n}}a_{\mathbf{n}}b^{\dagger}_{\mathbf{n}}b_{\mathbf{n}}$.

• •

Direct $a,b$ coupling (Eq. 17). Here, $\Omega_{\perp}$ is the overall strength of the field connecting the states $a$ and $b$, $\mathbf{k}_{L}$ is the in-plane projection of the momentum of the field mediating the transition, and $\omega_{L}$ its frequency. For simplicity, we can set $\omega_{L}\equiv 0$ and choose $\mathbf{k}_{L}=(0,0)$ (assuming it is coming perpendicular to the layers). With these simplifications, and expanding in the Wannier basis, the $H_{ab}$ Hamiltonian reads de Vega et al. (2008); Navarrete-Benlloch et al. (2011):

[TABLE]

whose spatial dependence will be same than for the intralayer couplings:

[TABLE]

since it comes from the overlap of the Wannier states between two sites, but whose strength can be tuned independently through $\Omega_{\perp}$.

Summing up, the global Hamiltonian describing the atom dynamics reads:

[TABLE]

where we can control: i) the inter/intra-layer hopping ratios ($J_{\perp}/J$) with $V_{D}$ and $\Omega_{\perp}$ as well as their range (through $L_{0}$); ii) the strength of the on-site atomic interactions $U_{\alpha},U_{ab}$ through Feshbach resonances; iii) the statistics of the particles with the atomic isotope. Even though we will not consider it for this manuscript, the range of the atomic interactions can be extended by dressing the $a,b$ states with a Rydberg level Gallagher (2005). Thus, we have a very unique platform to explore and simulate the physics of fermionic and bosonic twisted bilayers.

In the next section, we illustrate the potential of the proposed setup to simulate similar phenomena than in twisted van der Waals materials, in terms of observation of narrow bands or higher-order van-Hove singularities. For that, we calculate the emergent band structure and density of states for the twisting angles in which the two layers have a crystalline structure (Moiré patterns).

III.2 Moiré patterns and associated density of states

As it occurs with the honeycomb geometry Lopes dos Santos et al. (2007); Shallcross et al. (2010); Trambly de Laissardière et al. (2010); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011), only certain angles give rise to a crystal like configuration with a new periodicity (Moiré patterns). For square geometries the conmensuration angles can be written as:

[TABLE]

where we can restrict to $\theta(m,n)<\pi/4$ due to the symmetry of the layers.

Since the band structure and density of states are properties of the non-interacting part of the Hamiltonian we only need to consider $H_{a}+H_{b}+H_{ab}$ for the calculation. Furthermore, in this manuscript we restrict to the strong confinement limit ($L_{0}\ll\lambda/2$) such that we only include nearest neighbour hopping for the intralayer hopping with strength $J$, and on-site inter-layer hopping, $J_{\perp}(\mathbf{R})=J_{\perp}\delta_{\mathbf{R},0}$. For the large angles considered along this manuscript this minimal model will be a good approximation with small quantitative differences when longer range hoppings are included.

To calculate the associated density of states we define periodic boundary conditions, and diagonalize the Hamiltonian in momentum space for a discretization of $\mathbf{k}$ space of $N\times N$ points. The density of states is then given by:

[TABLE]

where $N_{\mathrm{M}}=2(n^{2}+m^{2})$ is the total number of $a/b$ lattice sites per unit cell (which is also the number of bands $\omega_{i}(\mathbf{k})$), and $\mathcal{N}(x)$ is a function that counts the number of states in a given frequency interval of width $\Delta\omega$, that we set to be $2\pi/N$ for the figures of this manuscript.

In Fig. 4(a) we plot the density of states at a fixed angle $\theta(2,1)\approx 36.9^{\circ}$, whose Moiré pattern is plotted in Fig. 4(b), for several $J_{\perp}/J\in[0,4]$. We observe several phenomena: i) the central van-Hove singularity splits for increasing $J_{\perp}/J$; ii) there appear extra divergences in the middle of the band; iii) at a critical $J_{\perp}/J\approx 1.7$ (for this angle) two of them separate as isolated bands in the lower/upper band part of the spectrum. These isolated bands get flatter as the angle decreases. This is clear in Fig. 4(c) where we zoom around these bands for fixed $J_{\perp}/J=4$, and plot the density of states for several angles ranging from $36.9^{\circ}$ to $12.1^{\circ}$. Though not shown, the critical ratio $J_{\perp}/J$ where these isolated upper and lower bands split from the rest also decreases with smaller angles. Thus the possibility to tune $J_{\perp}/J$ provide us a knob to observe narrow bands or the splitting van-Hove singularities with larger rotation angles.

To understand better the origin of the different contributions appearing in the density of states, we calculate the projected band-structure for $J_{\perp}/J=4$ and $\theta(2,1)$ along the square symmetry points, and plot it together with the associated density of states in Fig. 4(d). There, we can identify that the band-structure of the narrow upper and lower bands is similar to a square tight-binding model, e.g., with a saddle-point at the $X$ point, but with reduced band-width. As it happens in twisted bilayer graphene Ni et al. (2015); Sunku et al. (2018), the band-width is reduced because the hopping occurs mostly along the diagonals of the new primitive cell. Even more interesting is what happens in the middle of the band, where the accidental degeneracy between several bands give rise to different type of van-Hove singularities. For example, in $\omega\approx J(2.2J)$, six (two) bands touch at the $M(X)$-point leading to an asymmetric van-Hove singularity in the density of states. Similar asymmetric singularities have also been reported in twisted bilayer graphene Kerelsky et al. (2018), and explained in terms of high-order van-Hove singularities Yuan et al. (2019).

III.3 Considerations about other geometries

Along this manuscript we have focused on square-like potentials whose physics is dominated by the Van-Hove singularities explored in Fig. 4. For completeness, we want to note that it is possible to adapt some of the schemes discussed in II, e.g., the one using turnout wavelengths, together with existing experimental ideas Tarruell et al. (2012); Soltan-Panahi et al. (2011) to simulate Dirac like physics. The turnout wavelength scheme is convenient because it does not rely on having orthogonal polarization for the twisted standing waves, but rather on different frequencies (and angled beams), providing more freedom to choose the polarization for the rotated and unrotated potentials without having to worry about cross interference between them.

One possibility to obtain Dirac-like physics consists in keeping the square geometry behaviour but with staggered hoppings, forming the so-called brick-wall lattice. This configuration has already been realized experimentally for a single atomic state in Ref. Tarruell et al. (2012). For each layer it requires two retro-flected lasers in the X and Y directions with out-of-plane linear polarization and the same frequency, plus an additional retro-reflected laser in the X direction with slight different frequency such that its contribution adds up incoherently to the potential of the first two. With an appropriate choice of intensities (see Ref. Tarruell et al. (2012)), this potential gives rise to an energy dispersion with two Dirac points. However, since the symmetry of the layer is different than that of graphene, the underlying physics emerging under rotation will be different.

Another possibility consists in directly creating honeycomb lattices, as experimentally done in Ref. Soltan-Panahi et al. (2011) for a single atomic state. This can be done by sending three counter-propagating lasers within the same plane with $60^{\circ}$ angles between them and linear polarization for each of the turnout wavelengths, using angled beams for the one with larger momentum. In this case, the Moiré patterns occur for the same angles than in graphene Trambly de Laissardière et al. (2010):

[TABLE]

Since one of the features of this platform is the possibility to tune the inter/intra-layer hopping ratio, in Fig. 5 we calculate the change in the density of states (panel a), and band structure (panel b) for several $J_{\perp}/J$ and a honeycomb bilayer with fixed angle $\theta_{\mathrm{hc}}(2,1)$ (zooming around the Dirac point present when $J_{\perp}/J=0$). As in the square lattice model, we use a minimal hopping model restricting to nearest neighbour intra-layer hoppings and local inter-layer ones. Despite the simplicity of the model, we already observe how increasing the interlayer hopping can also be used to shift Van-Hove singularities close to Dirac points and to narrow the bands around the Dirac point for a given angle.

IV Probing twisted bilayer systems by spontaneous emission

Finally, to illustrate the potential of the proposed setup to explore physics beyond the one observed in solid state platforms, we place it in the context of the recent experiments simulating quantum optical phenomena with matter-waves Krinner et al. (2018). There, one studies the physics of one or several (simulated) emitters interacting with the matter-waves propagating in an structured optical potential. As shown in recent works Douglas et al. (2015); González-Tudela et al. (2015); González-Tudela and Cirac (2017, 2018); González-Tudela and Galve (2018); González-Tudela and Cirac (2018), these structured ”photons” can give rise to unconventional emitter dynamics and interactions with no analogue in standard photonic environments, and which can be instrumental, for example, for quantum simulation of long-range interacting spin systems Douglas et al. (2015); González-Tudela et al. (2015) or quantum chemistry problems Argüello-Luengo et al. (2018).

Apart from the states playing the role of matter-waves (the $a/b$ levels), in these experiments one requires an additional long-lived atomic state, $c$, to play the role of the emitter. Depending on the realization considered in section II.3, this can be either the ground state ${}^{1}S_{0}$ (for the ones relying in orthogonally polarized standing waves), or the other metastable excited state ${}^{3}P_{2}$ for the one exploiting turnout wavelengths. The advantage is that the potential for this state does not need to have the same periodicity than $V_{a,b}(\mathbf{R})$, such that the laser can have a different frequency. Another possibility is to use optical tweezers to trap them at controlled positions instead of retro-reflected lasers. In all cases, the potential depth must be deep enough such that once the atom is in the $c$ state in a given position it can not hop anywhere, but only emit into matter-waves. The Hamiltonian of these emitter-like atoms can be written as:

[TABLE]

where $U_{c}$ is the on-site interaction of this atomic state. As explained in section III, these interactions can be controlled, such that one can tune between having bosonic emitters when $U_{c}=0$ to purely two level ones $c_{i}\approx\sigma_{i}$ when $U_{C}\rightarrow\infty$.

The coupling between this internal atomic state and the photonic modes can be done in a similar way than the inter-layer hopping one (see Eqs. 17,22): either via a two-photon Raman transition through an auxiliary state, or with a direct microwave/optical transition depending on the level structure considered. Assuming we couple the $c$ state only to the $a$-modes, the simulated light-matter Hamiltonian in the strong confinement limit can be written as:

[TABLE]

where $\mathbf{n}_{i}$ denotes the lattice position where the $i$-th emitter-like atom couples to, and $g$ is the coupling strength which we assume to be real. For more details we refer the reader to the original Refs. de Vega et al. (2008); Navarrete-Benlloch et al. (2011) where all the details of this type of simulation were laid out.

To illustrate the physics that can emerge in these bilayer reservoirs we concentrate in the spontaneous emission dynamics of a single emitter. This means we assume to have a single emitter-like atom that is initially excited, $\left|\Psi(0)\right\rangle=c^{\dagger}\left|\mathrm{vac}\right\rangle$, and study its evolution under the global light-matter Hamiltonian given by: $H=H_{c}+H_{\mathrm{int}}+H_{B}$ as a function of time. Since $H$ conserves the number of excitations, the wavefunction at any time $\left|\Psi(t)\right\rangle=e^{-iHt}\left|\Psi(0)\right\rangle$ reads:

[TABLE]

Notice that since we have a single excitation in the system at any time, the interactions $U_{\alpha}$ or the statistics of the atoms play no role in the dynamics, which only depends on the system-bath coupling $g$, the inter/intra-layer coupling $J_{\perp},J$, and the detuning between the emitter-like and the bath frequencies $\Delta=\omega_{c}-\omega_{t}$. Perturbative treatments, like Markov approximation Cohen-Tannoudji et al. (1992), predict an exponential decay of the population, i.e., $|C_{e}(t)|^{2}\approx e^{-\Gamma_{M}t}$, with $\Gamma_{M}\propto g^{2}D(\omega_{c})$. However, since the density of states of these systems has a non-trivial structure, as shown in Fig. 4, the dynamics can differ significantly from this prediction.

To exemplify it, in Figs. 6(a-b) we plot the exact spontaneous emission dynamics of a single emitter coupled locally with $g=0.1J$ to a twisted square bilayer photonic reservoir with $J_{\perp}/J=4$. We consider two different situations: when the emitter’s frequency exactly matches the peak of the upper separated band (that we call the resonant case) in Figs. 6(a-b), and when the atomic frequency lies in the bandgap but close to lower edge of this separated band (that we call the off-resonant case) in Figs 6(c-d). The exact parameters are depicted in the legend and caption of Fig. 6. For the resonant configuration (see Fig. 6(a)) we observe that for the larger angle $\theta(2,1)$ the dynamics is mostly irreversible, following approximately an exponential decay. However, as the angle decreases the dynamics becomes more and more reversible (non-Markovian). The underlying reason is that the effective bath tunnelling rate is being reduced with the twisting angle and starts being comparable with $g$. For very small angles, the emitter exchanges excitations mostly with the localized mode that emerges in the Moiré supercell, leading the dynamics observed in the figure. This transition between Markovian and non-Markovian dynamics also occurs when the emitter’s frequency lies in the band-gap, shown in Fig. 6(c). In that case, the emitter goes from the no-decay situation predicted by Markov to a fractional decay one John and Quang (1994) in which the emitter decays partially and oscillates around a constant value. Thus, the twisting angle can be used as a knob to boost non-Markovian behaviour in such systems without increasing system-bath couplings.

Another interesting feature to consider is the behaviour of the matter-waves created in both the resonant and off-resonant situations, plotted in Figs. 6(b) and (d), respectively. When the emitter is resonant with the middle-band peak, the matter-wave pattern is highly anisotropic emitting mostly in $2$ directions which correspond to the diagonal directions of the emergent Moiré primitive cell. In the band-gap case, the bath excitations become exponentially localized around the emitter, forming a bound-state around them Bykov (1975); John and Wang (1990); Kurizki (1990). However, this bound state is anisotropic, being more extended again in the diagonal direction of the new primitive cell. We believe that the preferential emission (or localization) in the (off-)resonant cases along these new Moiré directions has a similar origin than the coupled-wire picture Wu et al. (2019) used to explain the results in twisted bilayer graphene, or the response of the material when it is probed by local plasmons Ni et al. (2015); Sunku et al. (2018) or preparing a localized electron wavepacket Do et al. (2019). From the quantum optical perspective, the main interest of such non-trivial matter-wave behaviour is that when putting more than one emitter, the induced collective decays ($\gamma_{ij}$) and dipole couplings ($J_{ij}$) will follow the same patterns Dicke (1954); Lehmberg (1970a, b) giving rise to very exotic equilibrium and out-of equilibrium many-body dynamics. For example, when the emitter’s frequency lie in the band-gap, the dynamics will be fully coherent and governed only by the dipole-dipole couplings:

[TABLE]

which in the blockaded limit ($U_{c}\rightarrow\infty$) corresponds to an effective spin model with Moiré-like interactions.

V Conclusions and outlook

Summing up, we have shown how to create two rotated state-dependent optical potentials to simulate the physics of twisted bilayers using cold-atoms. The possibility of tuning independently the inter/intra-layer hopping ratio opens the door to observe narrow bands or to move and split van-Hove singularities for larger rotation angles than in the solid state realizations. The potential of the proposed setup is three-fold: first, since the simulated particles interact locally and do not couple to phonons this setup can be used to discard some of the explanations introduced to explain the strongly correlated phenomena observed in experiments. Second, the possibility of observing similar physics with larger rotation angles simplifies ab-initio calculation and can be used as a platform to benchmark some of the effective descriptions made for the small rotation angles. Finally, as we have shown in the last part of the manuscript, the setup has the potential to explore different physics from solid-state platforms, e.g., quantum optical phenomena with structured Moiré photonic baths. To illustrate it, we have studied the spontaneous decay of a single emitter coupled to such type of bath, and predicted how the twisting angle can be used as a knob to enhance non-Markovianity, or to induce non-trivial emission patters and emitter-emitter couplings which follow the emergent Moiré geometry.

Taking this work as a basis, we believe there are still many open questions that will trigger more studies on the subject. For example, from the implementation point of view it will be interesting to find simpler atomic/laser schemes for the simulation of Dirac like physics as well as other ones that allow the simulation multi-layered materials that also show flat-bands. On the fundamental side, there are several relevant questions to answer such as the role of longer range hoppings when we tune away from the strong confinement situation, or whether we can observe the magic-angle behaviour (with a complete flattening of the band-structure) for larger angles by increasing $J_{\perp}/J$. From the perspective of the quantum optical experiments of matter-waves, it will be interesting to explore what kind of dynamics emerges from the higher-order van-Hove singularities that we found in the square twisted bilayers, and which also appear in the graphene case Yuan et al. (2019).

Acknowledgements

AGT acknowledges support from CSIC Research Platform PTI-001 and from the national project PGC2018-094792-B-I00 from Ministerio de Ciencia e Innovacion. JIC acknowledges the ERC Advanced Grant QENOCOBA under the EU Horizon 2020 program (grant agreement 742102). The authors acknowledge very useful discussions with S. Blatt, S. Fölling, I. Bloch, J. Knörzer, T. Udem, G. Giedke, and J. González.