Observation of Chiral-Mode Domains in a Frustrated XY Model on Optical Triangular Lattices

In a tight-binding approximation, the renormalized dispersion relation in the periodically driven triangular lattice is

where $d=\frac{2}{3}\lambda$ is the lattice constant. The parabolic dispersion relation term along the direction perpendicular to the triangular lattice plane is omitted. According to the dispersion relation Eq. S1, the minima of energy $\epsilon(J,J^{\prime},\mathbf{k})$ are given by

where we assigned the same names as in Struck et al. (2011) to the spin configurations that are realized with the different spin-spin couplings $J,J^{\prime}$. The corresponding phase $\theta_{i}^{gs}$ at each lattice site $i$ is given by:

where $\mathbf{R}_{i}$ is a lattice vector. Equation S6 does not include the cases for $J^{\prime}=0$ where the system forms decoupled chains and the dispersion relation along $k_{x}$ becomes degenerate. Energy of dispersion relation at the position of minima $\epsilon(J,J^{\prime},\mathbf{q}_{gs})$ is equivalent to the ground state energy of the classical XY model per particle $E_{XY}^{gs}(J,J^{\prime})/N$:

Figure S1(a) shows a zero-temperature phase diagram of the classical XY model in a triangular lattice based on Eq. S10.

The expressions of the phase modulation signals for the triangular lattice beams can be written as

where $\Omega/2\pi$ is the driving frequency, and $\phi_{x},\phi_{y}$ are the modulation amplitudes. Thus, the effective tunneling processes are expressed as

where $\mathcal{J}_{0}$ is the zeroth order Bessel function of the first kind. Figure S1(b) plots $J,J^{\prime}$ as a function of $\phi_{x},\phi_{y}$. Figure S1(c) shows experimentally observed interference patterns of unfrustrated phases. Figure S1(d) shows the effective band dispersion relations of F, Sp1, and Sp2 based on Eq. S1 using $J,J^{\prime}$ in Fig. S1(a),(b).

Figure S2 shows the experimental sequence from the holding time of the crossed FORTs after evaporative cooling to in-plane TOF. After 500 ms of holding atoms in the crossed FORTs for thermalization, we ramped up the optical triangular lattice to $V_{0}=1.0~{}E_{R}$ in 100 ms and simultaneously ramped down the crossed FORTs to a moderate depth so that the Flquet evaporation could efficiently work (see the S.4 Flquet evaporation). In the 2nd ramp of 100 ms, the optical triangular lattice was increased to $V_{0}=3.0~{}E_{R}$. Subsequently, the phase modulation amplitudes to drive the positions of the optical triangular lattice were ramped up with various times $\tau_{\rm U}$ as in Fig. 2 of the main text. Before the in-plane TOF, the vertical lattice depth increased to $V_{\rm ver}/k_{B}=285$ nK. Simultaneously, one of the phase modulation signals was decreased to zero to eliminate the micromotion effect that inherently arises in Floquet systems Goldman and Dalibard (2014); Eckardt (2017); Guo et al. (2019). Figure S2(b) displays the in-plane TOF images in Sp1 with several ramp-down times, $\tau_{\rm D}$. In the case of $\tau_{\rm D}<2\pi/\Omega=0.833$ ms, the atoms undergo a momentum kick (i.e., a micromotion) at the end of elliptical modulation of the lattice position. Consequently, the centers of the envelope in the momentum distributions shift. In contrast, if $\tau_{D}\geq 2\pi/\Omega$, the elliptical modulation reduces to a liner modulation, where the micromotion can be removed as long as the final phase of the remaining modulation signal is set to fulfill arg$[\phi_{y}(\tau_{\rm U}+\tau_{\rm D}){\rm cos}(\Omega(\tau_{\rm U}+\tau_{\rm D}))]=\pm\pi/2$. A too long $\tau_{\rm D}$, however, leads to relaxation from Sp1 to Sp2 via Staggered 1D chains (see Fig. S1(b)). We set $\tau_{\rm D}=1.0$ ms in our experiments.

In quantum systems with Flquet engineering, the heating effect during the periodic drive is the most serious problem. In the case of driven optical lattices, the major cause of this heating is multi-photon excitation between time-averaged energy bands Arlinghaus and Holthaus (2010); Weinberg et al. (2015). The periodic modulation with frequency $\Omega/2\pi$ causes the interband coupling processes that conserve quasimomentum but change the energy by integer multiples of the photon energy $\hbar\Omega$. An $n_{p}^{\rm th}$-order multi-photon transition is assumed to occur when the following resonance condition is fulfilled:

where $\Delta E_{\beta}^{\rm eff}(\mathbf{q})=\epsilon_{\beta}^{\rm eff}(\mathbf{q})-\epsilon_{0}^{\rm eff}(\mathbf{q})$, and $\epsilon_{\beta}^{\rm eff}(\mathbf{q})$ is the time-averaged single-particle energy of $\beta$-th band with quasimomenta $\mathbf{q}$.

Figure S3(a) shows the effective energy bands of the triangular lattice with $V_{0}=1.0~{}E_{R}$ in the case of Sp1. The red arrows indicate the two-photon transition from the lowest to the first excited band. These interband transitions caused by periodic driving significantly reduce the contrast $C$ defined in Fig. 2(a). The decrease in $C$ can be attributed to two distinct processes. First, the atoms excited into higher bands are no longer trapped in the optical lattice. Thus, the atomic losses reduce the bosonic enhancement $n|J^{(\prime)}|$. Second, an interacting BEC in the excited bands rapidly decays owing to the scattering processes, which results in a decrease in the degree of coherence in the system Martikainen (2011); Paul and Tiesinga (2013).

Figure S3(b) shows the results of a spectroscopic study of these multi-photon transitions in the driven triangular lattice. The excitation spectra were taken with various lattice depths at fixed forcing amplitudes $F_{A},F^{\prime}_{A}$. The resonant driving frequencies in Eq. S18 with different $n_{p}$ are plotted together with the excitation spectra. These $ab\;initio$ calculations agree with the experimental data, and multi-photon excitations up to the fourth-order are visible, which is consistent with a previous study Weinberg et al. (2015). We note that to keep $F_{A},F^{\prime}_{A}$ constant throughout the spectroscopy, the final amplitudes of phase modulation $\phi_{x}(\tau_{\rm U}),\phi_{y}(\tau_{\rm U})$ were adjusted according to Eqs. S16, S17. The maximum $C$ was reached within a narrow parameter region between $n_{p}=4$ and $n_{p}=5$. In the main text, we use $V_{0}=3.0~{}E_{R}$ and $\Omega/2\pi=1.2$ kHz for this reason.

As explained in the previous section, the triangular lattice depth $V_{0}$ and driving frequency $\Omega/2\pi$ are important parameters to avoid multi-photon interband transitions. Another salient factor that should be considered to reduce heating in Floquet-engineered quantum gases is Floquet evaporation Reitter et al. (2017). As in Fig. S4(a), we found that the crossed FORT depth during the Floquet drive should be decreased such that the resonantly scattered atoms can leave the trap before dissipating their energy into the system. We examined the effect of Floquet evaporation by applying a B-field gradient along the z-axis (Fig. S4(b)). Without the B-field gradient, heated atoms can evaporate because of gravity. When a B-field gradient was applied to cancel the gravity, the atomic ensemble suffered from heating and exhibited no interference peaks. When the B-field gradient was stronger than gravity, the visibility recovered, as expected.

Figure S4(c) shows the dependence of $C$ and $|\chi|$ for Sp1 on the vertical lattice depth $V_{\rm ver}$ in the driven triangular lattice. As $V_{\rm ver}$ increases, both $C$ and $|\chi|$ decrease, which indicates that Floquet evaporation along the z-axis does not work well in the presence of the vertical lattice because it prevents thermal atoms from evaporating along the gravitational direction. This could be a crucial problem for the experimental realization of the strong interaction regime ($U/n|J|\gg 1$) in the two-dimensional driven lattice systems.

We systematically investigated dependence on the crossed FORT depth and LS ramp-up time. Figure S5(a) shows the results. The contrast and chiral contrast of the middle column in Fig. S5(a) correspond to Fig. 2(c) in the main text. Figure S5(b) summarizes the fitting results. These results and Fig. S4(a) ensure that although the crossed FORT depth affects the maximum values of $C$ and $|\chi|$, it does not change the fact that the relaxation from F to Sp1 is much faster than that from F to Sp2.

FL count and $C$ in Fig. S5(a) are fitted by single or double exponential functions, while $|\chi|$ is fitted by a linear function with offset. Ideally, $|\chi|$ might also be fitted by an exponential-like function in the same way as the contrast $C$. In reality, however, the loss and heating of atoms cause both $C$ and $|\chi|$ to decrease in long ramp-up times. Since these loss and heating mechanisms are too complex, we excluded the experimental data of the long ramp-up times from fitting. Unlike $C$, the plateau of $|\chi|$, which is important if one wants to make a fit to the data with an exponential-like function, can hardly be seen. Therefore, we approximated the fitting function to a linear one in the case of $|\chi|$.

In Fig. 2(c) of the main text, we observed relaxation from F to Sp1 and Sp2. We found that the rise time of contrast $C$ differs between Sp1 and Sp2, which we attribute to the length of the path and the effective band structure. To ensure this scenario, we checked relaxation from R to Sp1 and Sp2 as seen in Fig. S6(a), which is the opposite case of relaxation from F to Sp1 and Sp2 used in the main text. Therefore, the relaxation time to Sp2 should be faster than that to Sp1.

We used a two-step ramp-up of the Floquet drive. The first ramp-up is to prepare R as an initial state. Figure S6(b) shows relaxation from F to R with various ramp-up times. We set the first ramp-up time at 60 ms, where atoms already become an R state. In the second ramp-up, atoms relax from R to Sp1 and Sp2. Figure S6(c) shows the results. As expected, the rise time of $C$ in Sp1 $t_{\rm rise}^{\rm Sp1}=24.0$ ms is slower than that in Sp2 $t_{\rm rise}^{\rm Sp2}=10.0$ ms, which is the opposite case of Fig. 2(c) in the main text. In making a fit to $C$, the definition of the empirical function in Eq. 2 is reversed such as

The reason why a double-exponential function is used for Sp1 in Eq. 2 and for Sp2 in Eq. S21 is because atoms go through a non-condensate state around the center of the phase diagram. There are two distinct processes; First, the distribution of atoms swiftly spreads in momentum space. Second, the atoms condensate again at a true ground state relatively slowly. FL count in the second ramp-up is fitted by a single exponential function instead of a double exponential function because the initial fast decay has already occurred in the first ramp-up. $|\chi|$ is fitted by a linear function with offset, the same as in the main text.

In Fig. 3(e) and Fig. 4(e) of the main text, the angle in $in\mathchar 45\relax situ$ images and FFT images and the widths in FFT images are plotted. In this section, we describe how the information was obtained.

As in Fig. S7, the $in\mathchar 45\relax situ$ images were fitted by a two-dimensional (2D) elliptic Gaussian function with tilt and offset to know the angles and widths of the atomic distribution before in-plane TOF. The fitting results were used for the numerical simulation of spin domains, the details of which are described in the following section S7.

As for the FFT signals, we extracted the information about angles and widths by using the following triple 2D Gaussian function;

where $A_{1},A_{2},A_{3},\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4},\theta,B$ are fitting parameters. The first 2D Gaussian (Gauss1) is for the inherent interference pattern of the triangular lattice, the second 2D Gaussian (Gauss2) for the fringes originating from spin domains, and the third 2D Gaussian (Gauss3) for noises in the experimental data. Gauss1 and Gauss2 have the same form. Because of this, although different initial guesses were used for $(\sigma_{0},\sigma_{1})$ and $(\sigma_{2},\sigma_{3})$ in making a fit, they were sometimes swapped. So, they were sorted after the fit. Before the fit, we applied a Gaussian filter with $\sigma=0.53~{}d^{-1}(=2~{}{\rm pixels})$ to the raw FFT signals, as in Fig. S7(b), to smooth out the spikes of the inherent interference patterns, whose spacing in FFT signal is $1.18~{}d^{-1}$. Figure S7(c) shows the Gaussian-filtered FFT signal of Sp2 together with the fitting results accumulated along the x-axis or y-axis. Figure S7(d) shows the FFT signals with Gauss3 subtracted to clarify that the experimental data is well fitted by Gauss1 and Gauss2.

Figure 4 in the main text shows the angle $\theta$ and widths $(\sigma_{2},\sigma_{3})$ of Gauss2 in Eq. S22 using the data of Sp2. We conducted the same analysis on the data of Sp1 and confirmed that the experimental data agree with the numerical simulation, as in Fig. S7(e).

To support our spin domains scenario, we numerically simulated the TOF of samples with spin domains. To simplify the numerical simulation, we assumed that Gaussian-shaped wavefunctions were located at each lattice site. The width of the wavefunction was extracted from the Gaussian fit to the Wannier function in the $ab\;initio$ calculation with a lattice depth $V_{0}=3.0~{}E_{R}$. The phases of Sp1 or Sp2 configurations for a homogeneous system were attached to the wavefunctions. We set domain walls perpendicular to the long axis of the elongated initial atom distribution and the opposite chiral modes for each region (Fig. S8(a)). The weight of the wavefunctions was determined by the initial density profile approximated by the 2D elliptic Gaussian function with the tilt angle (Fig. S7(a)). We obtained the density profiles after the TOF of 5 ms by summing all wavefunctions after free expansion. Here, we neglect the effect of residual trapping confinement owing to the vertical lattice.

The simulated density distributions with the number of domain walls from 1 to 3 are presented in Fig. S8(b). In the case of a single domain wall, one can see the effect of the interference between the two chiral modes. When multiple domain walls are formed, on the other hand, interference between the same chiral modes also occurs. We also checked the fast Fourier transformation of the density distributions, as we conducted for the experimental data in the main text, and observed that high-frequency components appeared along the long axis of the initial atom distributions (Fig. S8(c)). To quantify the FFT signals, we made a fit to the data in the same way as is described in the previous section S6. Since the numerical simulation is free from noise, we omitted Gauss3 in Eq. S22. Figure S8(d) shows the extracted widths in FFT space with various numbers of domain walls. The results tell us that the larger the number of domain walls, the larger the aspect of the FFT signals.