Simulating twistronics without a twist

We propose a scheme to emulate the essence of twisted bilayer graphene by exploiting ultracold atoms in an optical lattice. In our scheme, no bilayer nor twist are directly realized. Instead, two synthetic layers are produced exploiting coherently-coupled internal atomic states, and a supercell structure is generated via a spatially-dependent Raman coupling. We show that this system displays a band structure similar to that of magic angle twisted bilayer graphene, and explain its origin by deriving underlying effective Hamiltonians via perturbative approaches. Our proposal can be implemented using state-of-the-art experimental techniques, and opens the route towards the controlled study of strongly-correlated flat band accompanied by hybridization physics akin to magic angle bilayer graphene in cold atom quantum simulators.

The interesting correlated phenomenology is apparently related to Moiré patterns around small twist angles, the so-called magic angles, which lead to band flattening or effective mass reduction already at the single-particle level [22][23][24][25][26]. The geometrical Moiré patterns physically induce spatially varying interlayer couplings that are behind the strong modification of the band structure. As in artificial graphene systems [27], emulating this physics beyond materials research can allow identifying key minimal ingredients that give rise to the phenomenology of TwBLG while also providing additional microscopic control. Photonic systems are particularly suited to explore this physics at the single-particle level. Very recently, single-particle transport in tunable photonic Moiré lattices has been experimentally studied [28], where two dimensional localization of light and localization-delocalization have been experimentally demonstrated. Ultracold atoms in optical lattices [29,30] are the most promising platform to explore expeirmentally also the corresponding emerging manybody phenomena. The experimental realization of artificial graphene geometries [31,32], lattice geometries displaying flat bands like Kagome [33] and Lieb [34,35], or quasi-crystal structures [36][37][38], provides the building blocks for such exploration.
One obvious approach to study twisted bilayer graphene physics with ultracold atoms is to directly implement twisted bilayers using two intertwined optical lattices, as recently proposed in Ref. [39]. Schemes for simulating other bilayer heterostructures have also been put forward [40]. This direct strategy poses significant experimental challenges, as it is difficult to stabilize the two layers at relative small angles and simultaneously achieve a sufficiently large lattice containing several supercells. Here, we propose an alternative and highly flexible experimental approach, corresponding to effectively twisting the system without a twist. Our scheme builds on the concept of synthetic dimensions -reinterpreting the coherent coupling between spin states of an atom as tunneling along an artificial extra dimension [41][42][43]. As summarized in Fig. 1, we propose to realize a synthetic bilayer structure formed by two Raman-coupled atomic states with a spatially dependent synthetic tunneling Ω in analogy to the effects of Moiré patterns of twisted bilayer graphene. This leads to the creation of supercells with controllable sizes. By adjusting the strength, phase and spatial periodicity of Ω, this system allows engineering a broad range of band structures. In particular, magic values of the periodicity result in a band structure analogous to that of magic angle twisted bilayer graphene, including Dirac cones touching a set of quasi-flat bands. Our proposal could be realized with fermionic two-electron atoms, such as strontium or ytterbium, using available experimental techniques.
Concept. We consider a two-dimensional Fermi gas with four internal states, indicated here generically as spin states {m, σ} = ±1/2. The system is loaded into a spinindependent square optical lattice of lattice spacing d, which lies in the x − y plane and is characterized by a real tunneling amplitude t. We select two states to play the role of the electron spin σ = ↑, and the other two of electron spin σ = ↓. In addition, spin states corresponding to the same σ are coupled in pairs. We label them by the index m, and make them play the role of a synthetic layer dimension. Since m = ±1/2, we obtain a bilayer structure of synthetic layer tunneling given by the coherent coupling. In order to obtain a lattice geometry with a large supercell analogous to a Moiré pattern in twisted bilayer graphene, the amplitude of the synthetic layer tunneling amplitude is spatially modulated according to Ω(x, y) = Ω 0 [1 − α(1 + cos (2πx/l x ) cos (2πy/l y ))]. Here l x (l y ) is its periodicity along the x (y) axis. The synthetic tunneling also induces a Peierls phase γ · r, where γ = γx + γŷ and r = xx + yŷ. This mimics the effect of a magnetic flux that pierces the system perpendicularly to the synthetic layer dimension [42]. As depicted in Fig. 1(a), the complete scheme represents a synthetic spinfull bilayer structure subjected to a magnetic field, denoted as Θ(l x , l y ).
The Hamiltonian of the system is given by where we distinguish the in-layer and the inter-layer tunnelings.
To diagonalize it, we combine a gauge transformation and a Fourier transform such that a m,σ (r) = q exp (i(q · r + mγ · r)) a m,σ (q). Here, q is the momentum conjugated to r. The Hamiltonian can then be rewritten as H = q H q , where the dimension of H q is set by the spatial periodicity of the synthetic tunneling. Fig. 1(b) sketches the Brillouin zone of the bilayer system and a three-dimensional view of its energy spectrum for l x = l y = 4d, corresponding to Θ(4, 4), for γ = 0. In the vicinity of E = ±Ω 0 (1 − α), it features two quasi-flat bands and a Dirac point touching them (only one of them is represented in Fig. 1(b)). This band structure is reminiscent of that of magic angle twisted bilayer graphene.
Experimental proposal. For specificity, we focus on the realization of this scheme employing a subset of four states out of the large nuclear spin manifold I = 9/2 of 87 Sr. Note however that our proposal is directly transposable to 173 Yb (I = 5/2). Thanks to the SU (N ) invariant interactions characteristic of two-electron systems, collisional redistribution of the atoms among the different states is inhibited. We select two of them to play the role of the electron spin σ = ↑, and the other two of spin σ = ↓. All are subjected to a two-dimensional spin-independent optical lattice potential, created by two counter-propagating lattice beams. We choose λ L = 813 nm, which is commonly used because it corresponds to the magic wavelength of the clock transition 1 S 0 → 3 P 0 . We set a lattice depth 8 E L , which yields t/h = 107 Hz. Here, E L = 2 k 2 L /2m is the lattice recoil energy, k L = 2π/λ L , and d = λ L /2.
To create the synthetic layer tunneling, we exploit twophoton Raman transitions between spins m = ±1/2 [44]. We employ a pair of Raman beams of wavelength λ R = 689 nm near-resonant to the intercombination transition 1 S 0 → 3 P 1 , Each plane corresponds to one spin state m = ±1/2 (orange, green), which experiences a square lattice potential (tunneling t) and is connected to the other layer by a spatially-dependent and complex coupling Ω(x, y) (vertical red lines of variable width). A top view of the lattice indicating the unit cell of the system containing 2 × 8 sites for lx = ly = 4d is shown (black line). (b) Sketch of the first Brillouin zone, indicating the position of the high-symmetry points, and three dimensional view of the energy spectrum in the vicinity of E = −Ω0(1 − α) for Ω0α/h = 20t with α = 0.2 and γ = 0. It has two quasi-flat bands intersecting a Dirac point. Note that a simple square lattice supports neither flat-bands nor Dirac cones. (c) Proposed experimental realization. Top: Two retroreflected optical lattice beams (green) create the square lattice. Two Raman beams of opening angle θ (red) produce complex synthetic tunneling between the two layers. One "modulation laser" with a spatially varying intensity distribution (blue) modulates the amplitude of the Raman coupling. Bottom: laser beams involved in the synthetic bilayer coupling scheme. The single-photon detuning of the Raman beams (red arrows) is spatially modulated with respect to its initial value ∆0/2π ∼ 75 MHz using a laser beam blue detuned with respect to the 3 P1 → 3 S1 transition (blue arrow). It produces a light shift of maximal amplitude 2δ/2π ∼ 30 MHz.
which produce a coupling of amplitude Ω 0 = Ω 1 Ω 2 /∆ 0 . Here Ω 1 and Ω 2 are the individual coupling amplitudes of the Raman lasers and ∆ 0 the single-photon detuning. The Raman beams propagate in a plane perpendicular to the lattice potential, are aligned along its diagonal, and form an angle θ with the lattice plane (see Fig. 1(c)). This yields an in-plane momentum transfer per beam k R = ±2π cos θ/λ R , with projections k R / √ 2 along the lattice axes. Therefore, the phase of the synthetic tunneling is γ · r = γ(xx + yŷ), with γ = ±2π cos θλ L /( √ 2λ R ). The sign is determined by the relative detuning of the Raman lasers. Experimentally, the simplest choice is to use counterpropagating Raman beams (θ = 0 • ), which yields γ = 0.8 (mod 2π). However, other magnetic fluxes can be easily realized by adjusting the value of θ.
To implement a periodic modulation of the Raman coupling amplitude on the scale of several lattice sites, which is the key ingredient of our scheme, we propose to exploit a periodic potential created by a laser close-detuned from the excited state to excited state transition 3 P 1 → 3 S 1 (corresponding to 688 nm [46]). This results in a large light shift of the 3 P 1 excited state of amplitude δ, leading to a detuning of the Raman beams ∆(x, y) = ∆ 0 + δ(1 + cos (2πx/l x ) cos (2πy/l y )). Its effect is to modulate the Raman coupling amplitude [46,47]. We therefore name it "modulation laser". Band structures analogous to the one depicted in Fig. 1(b) are obtained for large values of αΩ 0 /h 20t = 10.7 kHz and spatial periodicities of the Raman coupling of several lattice sites [48]. The necessary patterns could therefore be projected by combining a spatial light modulator and an optical system of moderate optical resolution, ensuring a large flexibility.
Magic configuration and twisting analogy. The emerging band structures are sensitive to the spatial modulation and the strength of the laser coupling. Proper manipulations of the involved parameters can yield flat bands. Typically, a system with weak Raman coupling (Ω 0 α/h t) hosts a large number of extended hybridized bands. Enhanced coupling strength (Ω 0 α/h ≈ 10t) tends to foster band narrowing. Magic-angle TwBLG is characterized by flat-bands leading to a strong suppression of the Fermi velocity around the Dirac points. Analogous band structures can be realized for judicious adjustment of the Raman coupling periodicity, and as we discuss below, Θ(4, 4) turns out to be the configuration with smallest Moiré supercell (2 × 8) sites supporting a TBLG-like band structure. Fig. 1(b) illustrates spin degenerate bands around E/t = −Ω 0 (1 − α) for an exemplary case with Ω 0 α/h = 20t, α = 0.2 and γ = 0. The spectrum of the system is symmetric around E = 0, so that we only discuss the band structure at E < 0. The proposed experiment considers the effect of an additional magnetic flux, set by the parameter γ = 0.8. Narrow bands are formed for sufficiently large Ω 0 α at the energies −Ω 0 , −Ω 0 (1−α) and −Ω 0 (1−2α). The bands close to the energy −Ω 0 (1 − α) along the highsymmetry points are shown in Fig. 2 for (a) Ω 0 α/h = 2t and (b) Ω 0 α/h = 20t. In this case, two near-degenerate quasiflat bands at energy −Ω 0 (1 − α) and two Dirac cones appear. Such values of Ω 0 and α are reachable in our setup, leading to tunable bandwidth for the quasi-flat bands. Figs. 2 also shows the associated density of states (DOS), which is given by D(E) = L −d/2 i (E − E (k i )). The central peak at the energy −Ω 0 (1 − α) corresponds to a van-Hove singularity associated with the almost flat bands at sufficiently large Raman coupling.
Interestingly, band structures similar to the Θ(4, 4) case appear when l x = l y = 4νd, with ν integer. This can be explained by treating the intra-layer tunneling as a perturbation to the inter-layer tunneling. As explained in detail in the Supplemental Material [49], the nodal lines of the periodic modulation determine a bilayer Lieb lattice of sites. The two layers are energetically well separated with on-site energies ±Ω 0 (1 − α), respectively. The perturbation then induces tunnelings within the Lieb lattice topology, which at first order are composed of nearest neighbor tunneling matrix elements within a single layer. The Lieb lattice in its simplest form [50] is known to host a pair of Dirac cones intersecting at a single k point on a completely flat band, the Dirac point. The dispersion of the flat bands in the full model described by Eq. (1) originates from higher-order contributions in perturbation theory. More generally, for l x(y) = 4ν x(y) , where ν x(y) are posi- tive integers, a similar argument shows that the system can be effectively described by super-Lieb lattices with a supercell of 2(ν x + ν y ) − 1 sites. Moving away from these configurations leads generically to absence of Dirac-like physics similarly to moving away from magic angles in TwBLG [49].
Our proposal offers a powerful setup for engineering a wide range of band structures. By manipulating the periodicity and strength of the Raman coupling, and controlling the value of the chemical potential, it is possible to drive the system from insulating to semi-metallic and then metallic-type phases. An additional control parameter in our system is the artificial magnetic flux γ. It affects the band structure in a number of ways. As earlier, we focus on the case Θ M = Θ(4, 4) and the 6 bands closest to E/t = −Ω 0 (1 − α). Increasing γ leads to strong band narrowing. More interestingly, a non-zero γ opens a gap between the two quasi-flat bands at the Γ = (0, 0) point in the Brillouin zone. Moreover, the lower Dirac cone detaches from the lower quasi-flat band. This is reminiscent of the effect of a staggered chemical potential in the Lieb lattice [51]. However, in our model the upper quasi-flat band remains pinned exactly at the central energy (E/t = −Ω 0 (1 − α)) around the Γ point, and the upper Dirac cone remains gapless. Typical band configurations for two values of γ are shown in Fig. 3. The flux γ also controls the Dirac velocity of the cone which tends to zero as γ → π.
As already mentioned, the exact features of the spectrum can be traced via second order perturbation theory for large Ω 0 (1−α) [49]. The narrowing bandwidth and flattening of the Dirac cone are primarily driven by the same factor: the generically dominant (first order) nearest neighbour effective tunneling on the effective Lieb lattice decreases as −t cos (γ/2) when γ increases. On the other hand, the gap opening is driven by the tunneling modulation pattern breaking the C 4 (π/2 rotation) symmetry of the effective lattice for γ = 0.
Conclusions and Outlook. The basic element in the physics of TwBLG is the creation of large unit cells by rotating two layers with respect to each other. Around the magic angles, small rotations have a dramatic effect on the band structure of these systems. In this Letter, we have discussed a versatile method to create a new class of systems with controllable supercell structures for cold Fermi gases trapped in optical lattices. The size of the supercells is easily tunable and should allow addressing whether the physics of TwBLG is uniquely related to their macroscopic periodicity or indeed can be accessed for small unit cells. An inherent advantage of our optical-lattice-based construction is the possibility to modify over a wide range the interlayer coupling, which is controlled by a combination of optical Raman transitions and excited-state light shifts. We have shown that a square lattice synthetic bilayer displays a band structure that can be easily engineered by modifying the spatial periodicity, strength and Peierls phase imparted by the Raman lasers. It bears analogies to TwBLG in that it supports Dirac cones and quasi-flat bands in particular energy ranges and at certain magic periodicities.
The existence of identical scattering lengths parameterizing interactions between the atoms in the four internal states allows simulating the effect of both intra-and inter-layer interactions in the synthetic bilayer structure. The interacting Hamiltonian can be written as H I = U/2 r n(r)(n(r) − 1) where n(r) = m,σ a † mσ (r)a mσ (r), is the occupation of site (r) of the square optical lattice. The magnitude of U could be tuned by varying the transverse confinement. In particular, choosing a value of U smaller than Ω 0 α but much larger than the bandwidth of the quasi-flat band should allow achieving the strongly interacting regime in the latter. Projection of interactions onto the quasi-flat and hybridizing bands leads to extended Hubbard models with large on-site interactions as well as other terms, such as correlated tunneling. Probing such interacting systems at partial filling could potentially shed new light into theoretical debates on unconventional superconductivity [17,18,20,52], and topological order [21,53,54] in TwBLG. Finally, extending our approach to other lattice structures represents an exciting perspective for future studies. This describes a coupled bilayer system with intralayer tunneling − cos (γ/2) and interlayer tunneling −i sin (γ/2). Recall that the two layers actually correspond to different internal states on a single physical layer, so the mixing term is an effective spin-orbit coupling. The two layers have different on-site chemical potentials. We choose the lower layer to host the c +1/2,σ fermions and have potential −Ω(r).
The perturbation is highly effective if the potential Ω(r) displays equipotential connected regions, i.e. if the nodes of the periodic modulation (cosines) lie on the lattice. This certainly occurs when l x and l y are multiples of four. For concreteness, we discuss here the case l x = l y = 4 and 0 < α < 1. In this case the nodal lines of the periodic modulation determine a Lieb lattice (see Fig. S1) of sites with the same on-site energy. Furthermore, we focus on the layer with m = 1/2 and determine the spectrum centered around −(1 − α)Ω 0 [3]. At first order, the perturbation Eq. (2) partially lifts the degeneracy by inducing tunneling between sites on the Lieb lattice (black lines in Fig. S1).
Nearest neighbor tunneling on a Lieb lattice leads to a three band energy spectrum consisting of a completely flat band containing a Dirac point at which a pair of dispersive bands (with energy respectively higher and lower than the flat band) intersect. The Dirac point is located at the corner of the Brillouin zone (±π/2, ±π/2). This explains the origin of the band structure around energy −Ω 0 (1 − α) (and by analogy around Ω 0 (1 − α)) as shown in the main text). However, although the Lieb lattice has a unit cell of 3 sites (1 corner and two bridge sites), our full Hamiltonian has a unit cell containing six Lieb lattice sites (one choice is shown in Fig. S1). This full unit cell is recovered already in second order perturbation theory. The doubling of the Lieb lattice unit cell leads to folding of the Brillouin zone. Therefore, the two Dirac cones [4] in our system are at wave vector k = (0, 0). Moreover, the spectrum around −Ω 0 (1 − α) consists of six bands.
The detailed band structures, with quasi-flat bands, shown in the main text in Fig. 2 can be well recovered in second order perturbation theory for large Ω 0 (1−α), as seen in Fig. S3. The modification of the band structure is due to an additional periodic pattern of tunnelings in the Lieb lattice of the same periodicity of the supercell Θ(4, 4) generated by virtual tunneling to energy forbidden sites as shown in Fig. S2. The various terms described in Fig. S2 are gathered below.
Top view of the lattice for a modulation Ω(x, y) Ω0[(1 − α) − α cos (2πx/lx) cos (2πy/ly)], with lx, ly = 4. Black sites correspond to Ω(r) = Ω0(1−α), grey sites marked with ± correspond to sites with Ω(r) = Ω0(1 − α ± α). The latter are energy forbidden sites. Green boundary: unit cell of the lattice. Black lines denote the tunneling between the Lieb lattice sites generated. They are responsible for the main features of the band structures centered at ±Ω0(1 − α) described in the main text: a flat band intersected by the Dirac cones. The effects of the grey sites are only taken into account in second order perturbation theory (see Fig. S2).
Now, we comment on the evolution of the bands with the change of the flux γ. From Eqs. (3)-(11) we see that for γ = π the effective lattice consists of two inequivalent square lattices composed of the bridge and corner sites, respectively, that are completely decoupled from each other. The corner sites form a simple square lattice structure with lattice spacing 2d, while  Fig. S1) as well as from spin-orbit coupling interaction coupling to the other m = −1/2 lattice (grey sites in Fig. S1 as well as Lieb lattice sites. As a result for γ = 0, π, the chemical potential within a layer takes two different values µc on corner (black) and µ b on bridge sites (pink). Additionally a periodic pattern of next nearest and third neighbor tunneling is generated as shown by the colored lines (both full and dotted). The unit cell is that of the original lattice shown in Fig. S1. The values of the appropriate tunnelings are written in the text: t1(solid black lines, first order tunneling on the Lieb lattice), t Lieb (black and pink dotted lines), horizontal and vertical tunnelings t (2) ± between bridge sites over forbidden sites with Ω(r) = Ω0(1 − α ± α) (dotted orange (t    the bridge sites form a more complicated square lattice. However, neither of these support Dirac cones. These lattices have dispersive yet very narrow bands since the tunnel coupling is of the order of 0.01t. This is distinct from the opposite limit γ = 0, which is a Lieb lattice with additional tunnelings between bridge sites. Increasing γ decreases the total bandwidth of the 6-band system due to the dependence of the dominant coupling t 1 Eq. (3). The flat band and Dirac cone subsystem is modified as γ is increased. In particular, the upper Dirac -78 cone angle widens as a consequence of the decreasing bandwidth and therefore the Dirac velocity decreases. The upper flat band however exists pinned to the bare chemical potential −Ω 0 (1 − α) at the Brillouin zone center for all γ. We have found that the behaviour of the lower Dirac cone and the opening of a gap to the upper flat band is due to the C 4 symmetry breaking modulation in tunnelling between bridge sites. Indeed, on one hand this symmetry is explicitly unbroken where there is no gap. Moreover, we have also checked that setting the other possible factor, i.e. chemical potential staggering to zero does not influence the magnitude of the gap. The unimportance of the staggered chemical potential comes from the fact that the staggering in chemical potential is two orders of magnitude smaller than the staggering in tunneling.
Finally, in the context of twistronics, it is amusing to note that the two disentangled lattices at γ = π are rotated with respect to each other by π/4. By tuning γ towards 0, their coupling becomes stronger. This leads to a change in the band structure with the formation of Dirac cones that are attracted to each other at the center of the Brillouin zone towards the quasi-flat band.

Effect of phase shifts in the interlayer coupling
We briefly discuss the effect of allowing a phase in the spatial modulation of the synthetic layer tunneling amplitude, such that Ω(x, y) Ω 0 [(1 − α) − α cos (2πx/l x + φ x ) cos (2πy/l y + φ y )]. Small phases φ x and φ y displace the near-degenerate quasi-flat bands at energies ±(1 − α)Ω 0 away from each other. As a result the central peak of the DOS splits into a double peak structure for small values of these phases. Interestingly, a double peak structure in magic angle TwBLG has been reported in previous works [1,2]. Associated band calculations and DOS are shown in Fig. S4.

Non-magic configurations
As already mentioned before, the proposed scheme can be exploited to engineer a broad range of band structures by simply manipulating the periodicity of the spatial modulation of the Raman lasers in the square lattice under consideration, or by modifying the lattice structure in itself, along with other parameters, such as Raman coupling strength and magnetic flux. In order to demonstrate that the band structures obtained for Θ(4, 4) are not generic, we illustrate exemplary results for the configurations Θ(4, 7) and Θ (4,8) in Fig. S5(a-b). The chosen parameters have been selected for experimental convenience. We again focus in the vicinity of E/t = −Ω 0 (1 − α). While Θ(4, 7) supports a branch of isolated or hybridised flat bands, Θ(4, 8) supports semi-metallictype bands at E/t = −Ω 0 (1 − α) = −80.