Dipolar bosons in a twisted bilayer geometry

When two-dimensional periodic structures are overlapped and rotated with respect to each other, new periodic structures, moiré patterns, emerge. In recent years, condensed matter and material scientists have shown that this mechanism can be used to non-trivially modify properties of two-dimensional systems. At certain twist angles (magic angles), properties of bilayer graphene dramatically change Cao et al. (2018). What was a weakly-correlated Fermi liquid becomes a strongly-correlated system which supports superconductivity, quantized anomalous Hall states, and more Andrei and MacDonald (2020). The quantum anomalous Hall effect has also been observed with moiré hetero-bilayers Serlin et al. (2020); Li et al. (2021). The field of moiré bilayers remains active Xie et al. (2022); Zhao et al. (2024a); Tao et al. (2024); Zhao et al. (2024b).

Moiré bilayers have also been engineered with ultracold gases. The unprecedented level of control in these setups provides an ideal setting to simulate condensed matter systems in a highly controlled environment Bloch et al. (2012); Schäfer et al. (2020). A few years ago, a theoretical proposal to simulate twisted bilayers with cold atoms trapped in state-dependent optical lattices was put forward González-Tudela and Cirac (2019). In a recent experiment Meng et al. (2023), twisted bilayers at the magic angle $\theta=5.21^{\circ}$ have been artificially realized using bosonic atoms in two different spin states trapped in spin-dependent square optical lattices. This introduces yet another facet to the cold atoms toolbox – an ideal platform to investigate the intricate interplay between lattice geometry, interactions, and quantum fluctuations. The tunability of these systems, achieved by adjusting the twist angle between layers, provides a unique opportunity to manipulate quantum many-body systems at the microscopic level and explore the emergence of exotic states of matter. These studies are driven not only by a fundamental curiosity about the behavior of matter in these artificial structures but also by the potential applications in quantum simulation and quantum information processing.

Lattice bosons interacting via long-range, anisotropic dipolar interaction can stabilize various solid and supersolid phases Zhang et al. (2022, 2021, 2015); Suthar and Ng (2022); Chomaz et al. (2022). Recently, various solid phases have been experimentally observed with cold magnetic atoms Su et al. (2023). Dipolar bosons in layered geometries can also induce pairing leading to paired solid, paired superfluid, and paired supersolid phases Capogrosso-Sansone and Kuklov (2011); Safavi-Naini et al. (2013, 2014). In this work, we consider dipolar lattice bosons trapped in a twisted bilayer geometry. Layers are placed at a fixed distance $d_{z}$. Particles cannot hop between layers and a fixed dipolar interaction provides coupling between particles on different layers. We are interested in investigating the effects on the phase diagram of the bilayer geometry when a finite, small twist angle between the layers is introduced. Our findings can be summarized as follows. At small enough twist angle $\theta\sim 0.1^{\circ}$, the quantum phases and phase transitions remain consistent with those observed at zero twist angle Safavi-Naini et al. (2013). That is, paired superfluid, paired checkerboard solid, and paired supersolid phases are stabilized. However, a slight increase in the twist angle to $\theta=0.2^{\circ}$ disrupts the paired superfluid and paired supersolid phases in favor of phase separation between the checkerboard solid and superfluid phases. At the magic angle $\theta=5.21^{\circ}$, the local occupation number follows the moiré pattern of the underlying moiré bilayers so that a periodic structure of insulating islands is formed. These insulating islands are surrounded by a superfluid.

This paper is organized as follows: In Sec. II, we introduce the Hamiltonian of the system. In Sec. III, we discuss various phases, the corresponding order parameters and correlators. In Sec. IV, we present the results of the system for three small fixed values of twisted angles $0.1^{\circ},0.2^{\circ}$, and magic angle $5.21^{\circ}$; we outline our conclusions in Sec. V.

In the single band approximation, the system can be described by the Hamiltonian Safavi-Naini et al. (2013):

Here, $\alpha,\beta=1,2$ and $i,j$ label the layers and the lattice sites in each layer respectively, $a_{i,\alpha}(a_{i,\alpha}^{\dagger})$ are the bosonic annihilation (creation) operators obeying the hard-core constraint $a_{i,\alpha}^{\dagger 2}=0$, and $n_{i,\alpha}=a_{i,\alpha}^{\dagger}a_{i,\alpha}$ is the number operator for particles. The brackets $\langle\cdots\rangle$ denote summation over nearest neighbors only. The first term in equation 1 describes the kinetic energy with in-plane hopping rate $t$. The second term is the dipole-dipole interaction given by $V_{i\alpha;j\beta}=C_{\text{dd}}(1-3\text{cos}^{2}\gamma_{i_{\alpha};j_{\beta}})/4\pi|r_{i_{\alpha};j_{\beta}}|^{3}$, where $\gamma_{i_{\alpha};j_{\beta}}$ characterizes the angle between the direction of the polarization and the relative position of the two particles given by $r_{i_{\alpha};j_{\beta}}$, and $C_{\text{dd}}=d^{2}/\epsilon_{0}(C_{\text{dd}}=\mu_{0}d^{2})$ for electric (magnetic) dipoles of strength $d$. The repulsive nearest neighbor intra-layer interaction is $V_{\text{dd}}=C_{\text{dd}}/(4\pi a^{3})$, with $a$ the in-plane lattice constant. The attractive nearest neighbor inter-layer interaction is $V_{\text{dd}}^{\perp}=-2C_{\text{dd}}/4\pi d_{z}^{3}$, with $d_{z}$ the inter-layer distance. The quantity $\mu_{\alpha}$ is the chemical potential which sets the number of particles in each layer. Here, we fix $\mu_{1}=\mu_{2}$, i. e. $N_{1}=N_{2}$. The twist angle between these two layers is $\theta$.

In the following, we present unbiased results based on path-integral quantum Monte Carlo using the two-worm algorithm Söyler et al. (2009); Lingua et al. (2018). We have performed the simulations with system size $L\times L=N_{site}$ for each layer, with $L=20,22,33$ (we choose the lattice constant $a$ to be our unit of length). The filling factor is $n=N/N_{site}$ with $N=N_{1}=N_{2}$ the particle number on each layer. The inverse temperature $\beta$ is set to $\beta=L$. As for boundary conditions, each layer is separately mapped into a torus. As a result, lattice sites on the twisted layer which “fall outside” the region covered by the untwisted layer are considered within the untwisted region, and particle position and relative particle distances entering the hamiltonian are calculated accordingly. In other words, we set the simulation box in each direction to be $(1/2,L+1/2]$, see sketch in Fig. 1(b). With this approach, we maximize homogeneity of the system and reduce boundary effects so that our results highlight the effect of the twisted geometry on system properties. Due to computational limitations, our simulations are restricted to a system size of $L=33$. Due to this limitation, with open boundary conditions, we would expect boundary effects to be conflated with the effects of the twisted geometry. By using the torus-mapped boundary, we can avoid this issue and isolate the impact of the twisted geometry.

To discern superfluid phase (SF) from paired superfluid phase (PSF), we measure the superfluid stiffness $\rho_{s}$ for each layer and the superfluid stiffness for pairs $\rho_{\text{PSF}}$. The superfluid stiffness is $\rho_{s}=\langle\mathbf{W}^{2}\rangle/dL^{d-2}\beta$ Ceperley and Pollock (1989), where $\langle\mathbf{W}^{2}\rangle=\sum_{i=1}^{d}\langle W_{i}^{2}\rangle$ represents the expectation value of the winding number squared, $d$ is the dimension of the system (in our case, $d=2$), $L$ is the linear system size, and $\beta$ is the inverse temperature. The superfluid stiffness for pairs is directly associated with a pair condensate and is given by $\rho_{\text{PSF}}=\langle\mathbf{W}_{+}^{2}\rangle/dL^{d-2}\beta$, where $\mathbf{W}_{+}=\mathbf{W}_{1}+\mathbf{W}_{2}$ represents the sum of winding numbers in layers 1 and 2 and can be calculated within the two-worm path-integral quantum Monte Carlo. We note that, due to pairing across the layers, in the PSF phase the fluctuation of difference in winding numbers is zero, $\langle(\mathbf{W}_{1}-\mathbf{W}_{2})^{2}\rangle=0$. Diagonal long-range order, instead, corresponds to finite structure factor, $S(\mathbf{k})=\sum_{\mathbf{r},\mathbf{r^{\prime}}}\exp[i\mathbf{k}(\mathbf{r}-\mathbf{r^{\prime}})\langle n_{r}n_{r^{\prime}}\rangle]/N,$ with $N$ the particle number. $\mathbf{k}$ is the reciprocal lattice vector. Here, $\mathbf{k}=(\pi,\pi)$ to identify a checkerboard (CB) density pattern.

In addition, we also measure two-point correlation functions for each species $f_{\alpha}(r_{i_{\alpha};j_{\alpha}})\propto\langle a_{i\alpha}a^{\dagger}_{j\alpha}\rangle$, where the vector $r_{i_{\alpha};j_{\alpha}}=(x_{i_{\alpha}}-x_{j_{\alpha}},y_{i_{\alpha}}-y_{j_{\alpha}})=(x_{ij},y_{ij})$ is the relative position of the annihilation and creation operators at layer $\alpha$ (to ease the notation, we dropped the label $\alpha$ in the coordinates representation). Notice that, here, to calculate relative distances between annihilation and creation operators on each layer, we consider the relative position vector to originate at the site of $a^{\dagger}$ so that the maximum distance in $x$ and $y$ directions is $L-1$. Here, $\langle\rangle$ denotes a quantum and thermal average. $f_{\alpha}(r_{i_{\alpha};j_{\alpha}})$ is normalized so that its maximum value is 1. This correlator is expected to be long-ranged in a SF phase and short-ranged in an insulating phase.

Finally, we also produce density maps and condensate maps. The first is a map of the average of the local occupation number $\langle n_{i_{\alpha}}\rangle$ of a single typical Monte Carlo configuration. The latter is a map of the probability distribution for lattice sites to be visited by $a_{i\alpha}$ and $a^{\dagger}_{j\alpha}$ when their distance is larger than a chosen cutoff $R_{c}$. This map informs us on the lattice regions where off-diagonal long-range order exists.

Let us start with a description of the phase diagram in the absence of twist and at fixed $d_{z}=0.36$. The phase diagram is comprised of pair checkerboard solid (PCB), pair supersolid (PSS), paired superfluid (PSF), and, at lower dipolar interaction, independent superfluids (2SF) Safavi-Naini et al. (2013). In the PCB phase, reached at half filling and large enough dipolar interaction, atoms across the layers are strongly paired due to attractive interlayer interactions, i.e., particles sit on top of each other. The system can be thus envisioned as a solid of pairs. Upon doping the PCB solid with particles or holes, the PSS phase is stabilized. In this phase, CB density-density correlations coexist with a finite superfluid stiffness associated with pairs, i.e. $\rho_{\text{PSF}}\neq 0$ while $\langle(\mathbf{W}_{1}-\mathbf{W}_{2})^{2}\rangle=0$. For large enough doping, density-density correlation are lost (via a 2+1 Ising transition), and the system becomes a PSF. At lower values of dipolar interaction, a 2SF phase can be stabilized.

In the following, we fix the distance between the two layers to $d_{z}=0.36$, and the intra-layer dipolar interaction to $V_{\text{dd}}/t=0.26$. We denote the repulsive (attractive) nearest neighbor intra-layer (inter-layer) interaction by $V_{\text{dd}}=C_{\text{dd}}/(4\pi a^{3})(V_{\text{dd}}^{\perp}=-2C_{\text{dd}}/4\pi d_{z}^{3})$, with $a$ the in-plane lattice constant. These choices are made to compare with existing results on bilayer systems with no twist between the layers Safavi-Naini et al. (2013). We consider twist angles $\theta=0.1^{\circ}$, $0.2^{\circ}$, and magic angle $5.21^{\circ}$.

In conclusion, our investigation into dipolar bosons confined within a twisted bilayer geometry has uncovered a wealth of intriguing phenomena. We have elucidated the profound impact of twisted geometries on the quantum phases exhibited by the system. At a small twist angle of $\theta=0.1^{\circ}$, the system behavior remains consistent with that of the untwisted case. However, a slight increase in the twist angle to $\theta=0.2^{\circ}$ induces significant changes, disrupting paired phases and leading to the emergence of phase separation phenomena. This observation highlights the sensitivity of the system to slight alterations in the twisted geometry. At larger twist angle $\theta=5.21^{\circ}$, we have observed that particle density follows the pattern of the underlying moiré bilayer resulting in the emergence of periodic Mott insulator islands surrounded by superfluid regions, a phenomenon not observed in untwisted systems. This discovery highlights the novel and unexpected behaviors that can arise from the interplay between twisted geometries and dipolar interactions.

Our findings offer valuable insights into the interplay between geometry, interactions, and quantum fluctuations in correlated bosonic systems with potential implications for cold-atom experiments. By harnessing the precision and control afforded by cold-atom setups, the phenomena observed theoretically can be further explored and validated, advancing our understanding of complex quantum systems.