Frustrated quantum magnetism with Bose gases in triangular optical lattices at negative absolute temperatures

Quantum antiferromagnets with geometrical frustration exhibit rich many-body physics but are hard to simulate by means of classical computers. Although quantum-simulation studies for analyzing such systems are thus desirable, they are still limited to high-temperature regions, where interesting quantum effects are smeared out. Here we propose a feasible protocol to perform analog quantum simulation of frustrated antiferromagnetism with strong quantum fluctuations by using ultracold Bose gases in optical lattices at negative absolute temperatures. Specifically, we show from numerical simulations that the time evolution of a negative-temperature state subjected to a slow sweep of the hopping energy simulates quantum phase transitions of a frustrated Bose–Hubbard model with sign-inverted hoppings. Moreover, we quantitatively predict the phase boundary between the frustrated superfluid and Mott-insulator phases for triangular lattices with hopping anisotropy, which serves as a benchmark for quantum simulation. Classical computer simulations of quantum antiferromagnet exhibiting geometrical frustration are very demanding, and quantum simulation allows accessing high-temperature regimes where quantum effects are less relevant. By using a protocol for ultracold bosonic gases in optical lattices, the authors show that it is possible to achieve a regime of negative absolute temperature at which to study the physics of a frustrated Bose-Hubbard model.

F rustration is a key concept to understand various emergent phenomena in modern many-body physics 1,2 . When different interactions among particles strongly compete with each other, e.g., for a geometric reason, the system is "frustrated" in determining the true ground state. The study on the interplay of the frustration and strong quantum fluctuations has been one of the core challenges of quantum many-body physics, presenting many open problems in connection with nontrivial magnetic states including quantum spin liquids 3 and as a challenge for numerical techniques to handle highly entangled ground states 4,5 . Quantum simulation with the use of ultracold atomic gases in optical lattices [6][7][8] has been discussed as a promising approach to make a critical breakthrough in this research area. However, there still remain many challenges that have to be overcome in realizing and controlling frustrated quantum systems with cold atoms, whereas many theoretical proposals have been made [9][10][11][12][13][14] .
One straightforward idea for creating frustration is the use of two-component Fermi gases with (pseudo)spins σ = ↑, ↓ [15][16][17][18][19][20][21][22][23] loaded into a nonbipartite (e.g., triangular 24 and kagome 25 ) optical lattice. The second-order hopping process provides antiferromagnetic superexchange interactions between the (pseudo) spins σ 26,27 , which result in a frustrated situation, because complete staggered spin configuration is not allowed by the lattice geometry. Although long-range magnetic correlation over a distance comparable to the system size has recently been observed in a square optical lattice 21 , a further technical breakthrough is required to realize far lower temperatures to study frustrated quantum magnetism. Another interesting idea is a fast shaking of optical lattice, which can effectively invert the sign of the hopping integral from the natural one 28,29 . For ultracold Bose gases with sign-inverted hopping, the relative local phase of Bose-Einstein condensates (BECs) tends to be π on neighboring sites, analogous to antiferromagnetic spin coupling, which induces geometric frustration in nonbipartite lattices. A frustrated classical XY model has been successfully simulated with this technique 30 . However, it is challenging to reach a quantum regime of low temperature and density, because the lattice shaking can be a source of heating.
Recently, it has become realistic to create a well-controlled system at negative absolute temperatures 31 in laboratory. A state in thermal equilibrium is usually described by a statistical ensemble in which the lower-energy states are more occupied than higher-energy ones, obeying the probability proportional to the Boltzmann factor with temperature T ≥ 0. However, if the system has an upper energy bound, the opposite distribution with the largest occupation of the highest energy could also manage to be prepared. Such a state is characterized by a negative absolute temperature T ≤ 0. In a pioneer work, Braun et al. 32 have created a thermodynamically stable negative-temperature state of Bose gases in a square optical lattice by achieving the maximum interaction and potential energies in the regime of negligible kinetic energy. There the absolute temperature remains so low that the quantum phase transition from the Mott insulator (MI) to the superfluid (SF) has been observed.
Here we propose and examine a realistic route to create frustrated Bose gases in a quantum regime using the combination of the phase-imprinting techniques [33][34][35] and the statistics of negative absolute temperatures. Our proposal is based on the fact that a negative-temperature state of a systemĤ at T < 0 realizes the corresponding equilibrium state of the sign-inverted Hamiltonian ÀĤ at |T| > 0. Using this, one can achieve the same effect as inverting the sign of hopping, instead, by inverting the other factors, namely interatomic interaction and trap potential, and then preparing a negative-temperature state. To this end, we propose a phase-imprinting scheme combined with sudden sign inversion of the interaction and potential, which causes much less heating compared with the lattice shaking. Supposing Bose gases in a triangular lattice, we simulate the dynamics along the protocol within the time-dependent Gutzwiller approach (TDGA) to demonstrate that quantum phases of the frustrated Bose-Hubbard model, including chiral SF (CSF), could indeed be realized. Moreover, considering more general hoppings with spatial anisotropy, in which the hopping amplitude in one direction can be different from those in the other two directions, we give more quantitative analysis on the quantum phase transition between frustrated CSF and MI by means of the cluster mean-field plus scaling (CMF + S) method 36-38 with a twodimensional (2D) density matrix renormalization group (DMRG) solver 14 . This enables us to discuss the interplay of frustration and quantum fluctuations, which is a critical factor for producing various exotic frustrated states. The theoretical predictions provide a solid guidepost for future experiments to confirm that the interplay effects are properly captured in the quantum simulator. We set ℏ = 1, except in the figures.

Results
The Bose-Hubbard model on triangular lattice. A system of Bose gases in a deep optical lattice is described by the Bose-Hubbard model: Fig. 1 Triangular optical lattice with anisotropic hoppings. a Typical potential landscape of isosceles triangular optical lattice created by three laser beams in the case where one of the three lasers (in the a 1 = (a, 0) direction with a being the lattice constant) has a larger intensity than the other two (in the a 2 ¼ ðÀa=2; ffiffiffi 3 p a=2Þ and a 3 ¼ ðÀa=2; À ffiffiffi 3 p a=2Þ directions). b Schematic figures of the modeling with Bose-Hubbard parameters: the hoppings J 1 in the a 1 direction and J 2 in the a 2 and a 3 directions, and the onsite interparticle interaction U. The size of the arrows schematically shows the difference of the hopping amplitudes depending on the potential height between the neighboring lattice sites in each direction.
with hopping integral J ij (for i ≠ j), chemical potential J ii ≡ μ, onsite interaction U, harmonic trap potential V|r i | 2 ∕a 2 , and lattice constant a. Here we consider a triangular optical lattice with spatially anisotropic hopping of "isosceles type," which is parameterized by J ij = J 1 for nearest-neighbor (NN) sites (i, j) in the a 1 = (a, 0) direction and J ij = J 2 in the a 2 ¼ ðÀa=2; ffiffi ffi 3 p a=2Þ and a 3 ¼ ðÀa=2; À ffiffi ffi 3 p a=2Þ directions (see Fig. 1). The spatial anisotropy can be created by tuning the intensity of one of the three lasers to be different from the others (see Methods). The two extreme limits, J 1 /J 2 = 0 and J 1 /J 2 ≫ 1, are reduced to square lattice and one-dimensional (1D) chain, respectively.
Gutzwiller analysis. We first discuss the ground state for V = 0 within the site-decoupling Gutzwiller approach (GA) [39][40][41] to get a basic insight into the problem. The ground state in the weakcoupling regime (U ≪ |J ij |) is well described within the GA under the assumption of the BEC order hb i i ψ i ¼ ψe iqÁr i þφ with momentum q and global phase φ. The kinetic energy of Eq. (1) is Here, ε q À2ðJ 1 cos q Á a 1 þ J 2 cos q Á a 2 þ J 2 cos q Á a 3 Þ and M is the number of lattice sites. For natural-sign hoppings J 1 , J 2 > 0, the minimum kinetic energy is obtained at q = 0, leading to a uniform SF state. The maximum of the kinetic energy is achieved at q = ±Q with Therefore, if sign-inverted hopping (J 1 , J 2 < 0) is prepared, a frustrated CSF state with finite BEC momentum q = ±Q is realized. The choice of q = Q or −Q represents the degeneracy with respect to the chirality of vortex in unit triangles. Hereafter, we suppose q = Q to be spontaneously selected. In the equilateral case (J 1 = J 2 ≡ J), the momentum is Q = (4π/3a, 0) ≡ Q K . Therefore, the CSF state forms a "three-color" arrangement of the local phase (Arg½ψ i ¼ 0, 2π/3, and 4π/3 within a global phase shift) as shown in Fig. 2a. This can be understood as a compromise solution for the frustration in the bond energy minimization. For generic J 1 /J 2 > 0.5, the phase factor changes spatially with incommensurate pitch vector Q. In the parameter range 0 ≤ J 1 / J 2 ≤ 0.5 and the 1D limit J 1 /J 2 ≫ 1, a "two-color" (0 and π) pattern is formed with no chiral degeneracy. The J 1 /J 2 dependence of Q is analogous to the pitch vector of spin spiral states in spatially anisotropic triangular antiferromagnets [42][43][44][45] . For large, repulsive interaction U, lattice bosons undergo a quantum phase transition to the MI state when the filling factor ρ M À1 P i hn i i is an integer 46 (1), we calculate the order parameter ψ for ρ = 1 in the unfrustrated (J 1 , J 2 > 0: q = 0) and frustrated (J 1 , J 2 < 0: q = Q) ground states as a function of U/J 2 (see Methods). Figure 2b, c show the GA results for the SF-to-MI and CSF-to-MI quantum phase transitions, respectively. At J 1 = J 2 = J, the critical point at which ψ vanishes is given as U ðGAÞ c =jJj ¼ 17:5 for the frustrated case, which is a half of U ðGAÞ c =J ¼ 35:0 for the unfrustrated case. The strong reduction is attributed by the fact that the CSF state is less stable due to the frustrated local-phase arrangement in which the NN bonds are not fully satisfied in minimizing the local kinetic energy. For general values of the anisotropy J 1 /J 2 and ρ, the critical point is given by U ðGAÞ for the unfrustrated (frustrated) case. The ratio |ε Q /ε 0 | equals to 1 only at J 1 /J 2 = 0 or J 1 /J 2 → ∞, indicating that the reduction effect due to frustration exists even in the "two-color" region of 0 < J 1 /J 2 ≤ 0.5. It is noted that the curves of ψ for the unfrustrated system in Fig. 2b exactly overlap those for the frustrated system in Fig. 2c, respectively, for each J 1 /J 2 , by changing the scale of the interaction U by the factor |ε Q /ε 0 | (specifically, the factor 1/2 when J 1 = J 2 ). This is the case for general values of J 1 /J 2 and ρ. Indeed, the ground-state properties, including the density and order parameter, for the unfrustrated system with certain J ij , U > 0 and those for the frustrated system with negative hopping −J ij and repulsion |ε Q /ε 0 |U are identical within the GA, except for the spatial phase distribution Arg½ψ i ¼ Q Á r i . This is because the scale of the kinetic energy for the frustrated system is reduced by the factor |ε Q ∕ ε 0 | due to the frustrated phase configuration (see Methods).
Negative absolute temperature. To experimentally create the frustrated quantum states, it is required to prepare sign-inverted hoppings J ij < 0 with avoiding serious heating of the system. Below, we explain the details of the protocol through the use of negative-temperature statistics. Let us suppose an initial state in which N particles are distributed in a triangular lattice and a harmonic potential, which realizes the SF ground state of the standard Bose-Hubbard model (1) with J ij > 0 and U, V > 0. See a typical example in Fig. 3a obtained within the GA 41 for N = 1400, J 1 = J 2 = J = 0.08U 0 , U = U 0 , and V = 0.001U 0 (with U 0 > 0 being the energy unit). The phase Arg½ψ i is uniform in the unfrustrated SF state.
First, one needs to introduce the spatial phase distribution Arg½ψ i ¼ Q Á r i to create the frustrated CSF state. To this end, here we suggest the use of the phase-imprinting techniques [33][34][35] . When a single-particle energy difference δE is introduced between two sites, the relative phase on the two sites starts to develop with exp½iδEt in time t. In a region deep inside the SF phase, the kinetic energy (∝ε q ) reaches the maximum when the BEC momentum takes q = Q given in Eq. (2). One can transfer q from 0 to Q by introducing a temporary linear gradient potential Fig. 3b). Such a temporary potential can be created, e.g., by a Fig. 2 Predictions from the Gutzwiller approach. a The J 1 /J 2 dependence of the local phase pattern in the frustrated case J 1 , J 2 < 0 is illustrated. The colors indicate the sublattice structure in the local phase and the arrows depict the chirality. b, c The order parameter ψ as a function of |U/J 2 | (being the ratio of the interaction and the hopping in the a 2 and a 3 directions) at the filling factor ρ = 1 for b J 1 , J 2 > 0 and c J 1 , J 2 < 0. The value of ψ becomes zero at the transition from the superfluid (SF) or chiral superfluid (CSF) to Mott insulator (MI) phase. magnetic field gradient or by an extra 1D optical lattice satisfying Q = Q M or Q = Q K (see Methods). One has to perform the phase-imprinting operation in a much shorter time than the time scale of the hopping (δt ≪ |J ij | −1 ), to affect only the local phases. Besides, the local energy difference δE has to be set to a large enough value compared with |U| and |V| so that one can safely avoid the influence of the inhomogeneity of the density profile (see Methods). As such well-controlled phase imprinting has been successfully made in previous experiments 34, 35 , we will assume in the theoretical simulations presented below that a perfect phase imprinting can be achieved. The created CSF state with the "forced" phase distribution e iQÁr i should of course be dynamically unstable, as it has the maximum kinetic energy.
By changing U and V to be attractive (U < 0) and anti-trapping (V < 0), we make the interaction and potential energies also reach their maximum, to realize a stable negative-temperature ground state (at T ≈ −0) for J ij > 0 and U, V < 0. The created negativetemperature state should be equivalent to the ground state (at |T| ≈ +0) of the frustrated system with J ij < 0 and U, V > 0, as the physics of the two systems obey the same factor exp½ÀĤ=k B T. In ref. 32 for square (unfrustrated) lattice, this has been performed simply with U → −U and V → −V by using the Feshbach resonance 47 and blue-detuned lasers. In the present case with frustration, one has to pay special attention to the change of the kinetic energy scale after the sign inversion of U and V. Specifically, the initial state with the density (n i ) and orderparameter (|ψ i |) distributions shown in Fig. 3a after the phase imprinting Arg½ψ i ¼ Q Á r i is expected to correspond to the ground state of the frustrated Hamiltonian with the interactions rescaled by the factor |ε Q /ε 0 |, as explained above in the homogeneous case. Therefore, the values of U and V have to be changed as U → −|ε Q /ε 0 |U and V → −|ε Q /ε 0 |V (e.g., U → −U/2 and V → −V/2 for J 1 = J 2 ), to adjust the energy scale in consideration of the reduction of the kinetic energy due to the frustration.
TDGA simulation. In the framework of the TDGA 48,49 , we simulate the time evolution of the initial state in Fig. 3a after suddenly performing the three operations shown in Fig. 3b. In the simulation, we implement the phase imprinting on the initial state by the operation is the coefficient of the Fock basis n j i. In addition, we perform the sudden changes of U = U 0 → −U 0 /2 and V = 0.001U 0 → −0.0005U 0 . After those three operations are performed at t = 0, we calculate the time evolution of the state fixing the value of |U/J| for 0 ≤ 200U À1 0 to see the stability of the created negative-temperature CSF state.
As shown in Fig. 3c, the created negative-temperature state is predicted to be indeed dynamically stable for a long time 0 ≤ 200U À1 0 . As a reference for the comparison, we also simulate the case with the same settings but without the phase-imprinting operation (Fig. 3d). In this case, as shown in Fig. 3e, the state collapses immediately within t ≲ 2U À1 0 due to the dynamical instability; the order parameter rapidly decreases, although the density profile is kept. The difference between the two cases (with and without the phase imprinting) can be clearly seen in the time evolution of the density fluctuation δn 2 hn 2 i i À hn i i 2 averaged over the center sites within |r i | ≤ 10a. As shown in Fig. 3f, the value of δn 2 rapidly increases and then exhibits an irregular oscillation in the case without the phase imprinting, whereas it is almost constant in the case that the negative-temperature state is properly created by simultaneously achieving the maximum kinetic, interaction, and potential energies.
The CSF-to-MI phase transition. For the stable negativetemperature CSF state shown in Fig. 3c, we slowly increase the value of |U/J| for t > 200U À1 0 to observe the CSF-to-MI transition. In experiments, the tuning of |U/J| is performed by controlling the height of the optical lattice 46 . Here we assume that |U/J| increases simply through J decreasing linearly with The three operations to make frustrated negative-temperature states (for spatially isotropic hoppings J 1 = J 2 = J). The single-particle spectrum ε q in units of J is plotted as a function of the wave vector q in units of πa −1 . The illustrations depict a tilting of the optical lattice by energy difference δE for the phase imprinting, the inverting from replusive (U > 0) to attractive (−U/2 < 0) interaction, and the inverting from confinement (V > 0) to anti-confinement (−V/2 < 0) trap potential, respectively. c Negative-temperature CSF state after time 200 in units of the inverse of the initial interaction strength U À1 0 , which shows its stability for a sufficiently long time. The case without the phase imprinting operation (d) is shown in e for a comparison. f The evolution of the density fluctuation δn 2 of the initial state with and without the phase imprinting.
0 . As shown in Fig. 4a, when |U/J| increases, the MI plateau is gradually formed and |ψ i | in the trap center decreases towards zero. Until the transition point, the three-color phase profile in the CSF state is properly kept within a global phase shift (Fig. 4b, c).
To show the transition process more clearly, we plot in Fig. 4d the time evolution of the scaled value of the density fluctuation δn 2 (t) ∕ δn 2 (0). There we see that δn 2 (t) ∕ δn 2 (0) decreases with |U/J| and almost vanishes at jU=Jj % jU ðGAÞ c =Jj ¼ 17:5, which is consistent with the GA prediction of the critical point for the frustrated system. In the case of a box-shaped trap potential 50,51 , which is modeled by V = 0 and the open boundary at |r i | = 36a, the transition is more sharply observed (the blue dashed line in Fig. 4d for N = 4692). We also plot in Fig. 4e the unfrustrated case of the standard SF-MI transition as a reference, which shows a sharp difference from the frustrated case in the values of |U/J| of the transition region.
At the end of this section, let us briefly mention the validity and limitation of the TDGA method with respect to the calculations presented above. The mean-field analysis with the site-decoupling approximation reproduces the exact wave function in the weak-coupling (U ≈ 0) and strong-coupling (U → ∞) limits, and is thus expected to reasonably interpolate the two limits in two or higher dimensions 6,39,52 . In fact, the TDGA simulation has been often used to describe the dynamics of Bose-Hubbard systems [53][54][55][56] , which should be fairly reliable, at least, in the absence of strong quantum correlations, e.g., in the deep-SF regime or in three dimensions (3D) 57 . It is therefore expected that the above TDGA simulation gave a proper description for the dynamical stability of the negativetemperature CSF state after the quench in the deep-SF regime. On the other hand, the site-decoupling treatment could fail to capture some quantitative features in the vicinity of the CSF-MI transition, including the location of the critical point and the values of critical exponents, as the intersite quantum correlations become important.
Quantitative analysis by CMF + S with DMRG. We provide more quantitative estimation of the CSF-MI critical point for V = 0 beyond the site-decoupling approximations (GA and TDGA) to predict the interplay effect of frustration and quantum fluctuations as a guideline necessary for experiments. The quantum effects on frustrated Bose gases in 2D lattices have been poorly studied due to the lack of effective computational methods. Here we generalize and apply the numerical CMF + S method with 2D DMRG solver, which recently established in studies on frustrated quantum spins 14 , to the present system. We consider the N C -site cluster Hamiltonian on a triangular-shaped cluster (see the illustrations in Fig. 5a) under the mean-field boundary condition. We work in the twisted frame,b i e ÀiQÁr ib i , with optimizing Q = (Q x , 0) (see Methods). The CMF + S method permits the systematic inclusion of quantum intersite correlations by increasing N C , which connects between the GA (N C = 1) and exactly quantum (N C → ∞) results. Figure 5a summarizes the N C = 10, 15, 21 data and the CMF + S (N C → ∞) result for the ρ = 1 CSF-MI (SF-MI) critical point in the frustrated (unfrustrated) case with J 1 , J 2 < 0 (J 1 , J 2 > 0). The value of U c /|J 2 | for the frustrated case exhibits a nonmonotonic behavior with a dip around J 1 /J 2 ≈ 0.8, whereas the unfrustrated one simply increases as the total hopping J 1 + 2J 2 increases. This is indeed the frustration effect, which destabilizes the CSF state. Besides, the value of U c is strongly reduced from the GA prediction U ðGAÞ c due to the inclusion of the intersite quantum correlations. It should be noted that the relative difference ðU ðGAÞ c À U c Þ=U ðGAÞ c between the GA and CMF + S values is much larger for the frustrated case (~40-50%) than the unfrustrated case (~20%) as shown in Fig. 5b. This indicates that the quantum effects are strongly enhanced by the interplay with frustration. Figure 5c shows that the BEC momentum Q x in the CSF state is little affected by the inclusion of quantum correlations. The slight variance of Q x from 4π/3 at J 2 /J 1 = 1 is thought to be due to the finite cluster-size effect (see Methods).
Detection method. The BEC momentum Q in the CSF state and the transition to the MI can be simply detected by time-of-flight (TOF) images of momentum distributions 46 . A more precise determination of the critical point can be made by extracting the condensate fraction from the TOF images 58 , by observing the critical velocity using a moving optical lattice 57 , or by measuring the density fluctuation δn 2 using the quantum-gas microscope [18][19][20][21][22][23] . The frustrated CSF-MI transition is easily distinguishable from the standard SF-MI transition thanks to the sharp difference in the value of U c /|J 2 |.

Discussion
We made a proposal and provided the necessary theoretical analysis for analog quantum simulation of frustrated quantum magnetism by using ultracold Bose gases in triangular optical lattices. We proposed an experimental protocol to create a frustrated quantum state at negative absolute temperature by performing a phase imprinting together with sudden inversion of the interatomic interaction and the trap potential. Simulating the time evolution, we demonstrated that a dynamically stable SF state with chiral symmetry breaking was indeed realized and underwent the quantum phase transition to the MI as the hopping amplitude decreased. Moreover, we performed state-ofthe-art numerical calculations on the quantum critical point as a function of the spatial hopping anisotropy, which predicted a significant interplay of frustration and quantum fluctuations. The quantum dynamics in the presence of such enhanced fluctuations near the criticality is out of reach with the currently available numerical methods and thus gives a strong motivation for future experimental quantum simulations as an interesting and important subject to be investigated. A connection of the present synthetic system to real materials of frustrated antiferromagnets can be obtained by using the approximate mapping from the bosonic to spin-1 operators 59 which is valid in the vicinity of the transition between the CSF and MI phases at integer fillings ρ. The Hamiltonian (1) is mapped onto the spin-1 XY model with the XY spin exchange J XY ij ¼ ÀρJ ij and the single-ion anisotropy D = U/2. Therefore, the physics discussed in the present study is deeply related to the pressure-induced phase transition from the "120 ∘ " magnetic order (corresponding to CSF) to nonmagnetic state (MI) in spin-1 easy-plane triangular antiferromagnets such as CsFeCl 3 60,61 , which has recently attracted considerable attention in the connection with the novel excitations near quantum criticality 62 . From the standpoint of fundamental statistical physics, the CSF-to-MI transition should be associated with the spontaneous U(1) Z 2 symmetry breaking with respect to the global phase and the chirality determined by q = Q or −Q, whose quantum critical phenomena and universality class have not yet been established. The present bosonic system of synthetic antiferromagnets is advantageous for exploring exotic quantum critical phenomena in low dimensions, whereas the real materials have strong 3D couplings between the triangular layers 60,61 .
Moreover, it has been expected that adding long-range interatomic interaction to the present system may give rise to an exotic chiral MI state 63 in between possible separate U(1) and Z 2 symmetry breakings. This is essentially equivalent to the socalled "chiral liquid" expected in a spin-1 frustrated magnet 64 . Besides, our protocol for direct creation of a frustrated quantum state is advantageous for preparing the phases that are not neighboring to the MI phase, such as quantum spin liquids expected for ρ = half integers 3,11 . Thus, the present study gives a crucial guidepost for cold-atom quantum simulations of those exotic quantum frustrated physics.

Methods
Triangular optical lattice with anisotropic hoppings. A triangular optical lattice can be created by superposing three laser beams that intersect in the x-y plane with wave vectors k 1 = k L (1, 0), k 2 = k L ðÀ1=2; À ffiffi ffi 3 p =2Þ, and k 3 = k L ðÀ1=2; ffiffi ffi 3 p =2Þ, and equal frequency ω L 24 . All beams are linearly polarized orthogonal to the plane and each has field strength E i (i = 1, 2, 3). The total electric field is given by The dipole potential generated by the electric field is proportional to its squared amplitude, , and A(t) represents the terms dependent on time t. As the frequency of light is quite large, only the time-averaged value of |E tot | 2 can affect atoms. All the terms in A(t) oscillate at frequency 2ω L and thus can be dropped. Finally, we obtain a periodic dipole potential where V 1 , V 2 , and V 3 are proportional to E 2 E 3 , E 1 E 3 , and E 1 E 2 , respectively. A variation of the phases ϕ i of the laser beams yield only a global shift of the lattice in position. The primitive lattice vectors a 1 and a 2 are given so that a i ⋅ b j = 2πδ ij . We define a 3 = − a 1 − a 2 for convenience. For red-detuned lasers with V i > 0, the potential minima form a geometrically equilateral triangular lattice in lattice constant: a = |a i | = 4π/3k L = 2λ L ∕ 3 with λ L being the laser wavelength. The spatial anisotropy in the hopping amplitudes can be introduced by the difference in V 1 , V 2 , and V 3 through tuning the laser intensities E 1 , E 2 , and E 3 . For example, the set of the field strengths with the relation of yields the potential landscape shown in Fig. 1a, which gives anisotropy of "isosceles type" in the nearest-neighbor hopping amplitudes.
Tilting triangular optical lattice with an additional 1D optical lattice. To perform the phase-imprinting process in the protocol proposed here, one has to introduce a single-particle energy difference δE between adjacent two sites. This could be directly implemented by a magnetic field gradient, which introduces an extra linear gradient potential. Here, let us also provide another way to perform the phase-imprinting process for preparing the commensurate Q = Q M ≡ (2π/a, 0) (two-color) and Q = Q K ≡ (4π/3a, 0) (three-color) states by the use of an additional 1D optical lattice. We suppose that a potential of 1D optical lattice is created with additional laser beams in the a 1 direction: with amplitude V 0 , wave vector k 0 L ¼ 2π=λ 0 L , and phase ϕ 0 . Here, let us set ϕ i = 0 (i = 1, 2, 3) for the triangular optical lattice without loss of generality.
Let us first consider the range of 0 ≤ J 1 /J 2 ≲ 0.5 (V 1 ≫ V 2 = V 3 ), in which the configuration of the local phase factor is expected to form the two-sublattice (say A and B) structure with the pitch vector Q = Q M as illustrated in the left panel of Fig. 2a. For creating the local phase configuration by the phase imprinting, it is required to introduce a temporary single-particle energy difference δE only between the twosublattice groups of sites for time δt = 2π/δE. This can be achieved by using an additional 1D optical lattice of magic wavelength defined by λ 0 L ¼ 4λ L =3 and phase shift ϕ 0 ¼ 0 or π/2 (mod π). Figure 6a, b show an example of the total potential V(r) + V 1D (x) with the parameters 2V 1 ¼ V 2 ¼ V 3 ¼ V 0 and ϕ 0 ¼ 0. It is noteworthy that the two options in ϕ 0 correspond to the exchange of A and B.
In a similar way, using a 1D periodic potential with λ 0 L ¼ 2λ L , one can also implement a temporary energy difference δE between the three-sublattice groups of sites, say A, B, and C. The additional lasers are shined for time δt = 4π/3δE to imprint the three-color phase configuration illustrated in the middle panel of Fig. 2a to the initial SF state. Figure 6c, d show an example of the total potential V(r) + V 1D (x) with the parameters V 1 ¼ V 2 ¼ V 3 ¼ V 0 and ϕ 0 ¼ π=12. It is noteworthy that the phase shift ϕ 0 has six options, (2n − 1)π/12(n = 1, 2, ⋯, 6), reflecting the possible permutation of A, B, and C.
The GA analysis for finite-momentum BEC states. Within the site-decoupling mean-field approximation, known as the GA, the effective local Hamiltonian at site i is given byĤ as the result of the decouplingb y ibj % ψ jb y i þ ψ Ã ibj À ψ Ã i ψ j in the original Hamiltonian (1). The results displayed in Fig. 2b are calculated under the assumption of finite-momentum BEC, ψ i ¼ ψe iqÁr i þφ , for V = 0. In this case, one has only to consider a certain single site i, e.g., the site at the origin r i = 0, and φ = 0 without loss of generality. Equation (8) becomeŝ It is worth noting that the effective one-body HamiltoniansĤ GA i in the unfrustrated and frustrated cases differ only in ε q (ε 0 or ε Q ). Therefore, if the values of all the other terms are multiplied by |ε Q /ε 0 |, the results of the GA calculations, such as the value of |ψ i |, for the unfrustrated system become exactly the same as those for the frustrated system, except for the chiral phase distribution e iQÁr i and the overall energy scale (which is also multiplied by |ε Q /ε 0 |).
In the presence of the trap potential V ≠ 0, the mean field ψ i can no longer be assumed to have spatially uniform amplitude. Therefore, one has to deal with Eq. (8) on the entire lattice sites, each of which is connected to the six neighboring sites through the mean fields fψ i ± a n ; n ¼ 1; 2; 3g. To prepare the initial state shown in Fig. 3a, we solve the set of self-consistent equations ψ i ¼ P n ffiffiffi n p f ðnÀ1ÞÃ i f ðnÞ i for all sites within a cutoff length l c from the trap center. We take l c = 36a. The chemical potential μ is determined from the global number equation P i hn i i ¼ N, which must be solved simultaneously with the equations for ψ i . To efficiently achieve the convergence in the self-consistent calculations, we take the uniform solution for ψ i obtained in the absence of the trap potential as the initial condition and employ the Newton-Raphson method.
The TDGA simulation in a trap potential. The TDGA equation is given by  65 . We first rewrite Eq. (11) in a matrix form: where f i is the vector with the components f ðnÞ i (0 ≤ n ≤ n max ) and H GA i is the Fig. 6 Periodically tilted triangular optical lattice. a Potential landscapes of triangular optical lattice without and with an additional 1D periodic potential (whose amplutude and phase are V 0 and ϕ 0 , respectively) for 2V 1 The amplutudes of the three standing waves forming the triangular optical lattice are denoted by V 1 , V 2 , and V 3 . The cuts along r = r(1, 0) and r ¼ rð1=2; ffiffiffi 3 p =2Þ are shown in b. The spatial coordinate r = (x, y) and its norm r are measured in units of the lattice constant a. c, d Same as a and b, respectively, for ðn max þ 1Þ ðn max þ 1Þ matrix form of the GA local Hamiltonian, which is timedependent via the order parameters ψ j (t) of neighboring sites of site i. The explicit components of H GA i are easily obtained from Eq. (11). The value of f i after time evolution in a short time Δt is given by where I is the identity matrix. Here, Δt must be sufficiently shorter than the time scale of the physics considered. The time evolution of the system is numerically calculated, step by step, with short time interval Δt. To calculate f i (t + Δt) at site i, the values of ψ j (t) on its neighboring sites are required. Therefore, one has to calculate ψ j ðtÞ ¼ P The phase-imprinting operation. In the section "TDGA simulation", we assume that the phase imprinting can be theoretically implemented by the operation P n f ðnÞ i n j i ! P n e inQ K Ár i f ðnÞ i n j i on the local wave functions. Here we provide a brief discussion on the experimental conditions for achieving such perfect phase imprinting.
Let us consider the same initial state as shown in Fig. 3a (with the parameters N = 1400, J 1 = J 2 = J = 0.08U 0 , U = U 0 , and V = 0.001U 0 ) and its time evolution within the TDGA method in the presence of a temporary linear gradient potential V ext ¼ δE P i ðx i =aÞn i . As seen in Fig. 7a, the chiral structure in the local phase, Arg½ψ i ¼ Q Á r i , could be successfully imprinted by applying the temporary potential for δt = Q x a/δE (=4π/3δE for J 1 = J 2 ). Figure 7b indicates that the phase imprinting becomes almost perfect with no changes other than the local phase distribution when δE exceeds~10U 0 . The corresponding imprinting time δt ≲ 0.4U 0 −1 is much shorter than the typical time scale of the experiments on the SF (CSF)-MI transition (see Fig. 4). Fig. 7 Phase imprinting. a Time evolution of the phase difference between the neighboring sites in the a 1 direction, δφ ðArg½ψ i À Arg½ψ iþa 1 Þ, averaged over the center sites within |r i | ≤ 10a, in the presence of a temporary linear gradient potential. The colors correspond to different potential strengths δE/U 0 = 1, 2, 5, 10, and 20, although all the curves are almost overlapped. b The density (solid lines) and order-parameter (dashed lines) profiles in the final state after applying a gradient potential for time δt. Fig. 8 Cluster mean-field calculations with 2D density matrix renormalization group. a Mapping of the 2D cluster problem with the mean-field boundary onto an equivalent 1D chain with long-range interactions and mean fields. b Typical behavior of the provisional critical point U Ã c =jJ 2 j as a function of the given momentum q = (q x , 0) in units of πa −1 for the hopping anisotropy J 1 ∕J 2 = 0.7 and the cluster size N C = 21. Fig. 9 Scaling analysis in the cluster mean-field plus scaling analysis. Cluster-size scalings of the critical points for the transitions (a) between the superfluid and Mott insulator states in the unfrustrated (J 1 , J 2 > 0) case and b between the chiral superfluid and Mott insulator states in the frustrated (J 1 , J 2 < 0) case. The extrapolated (λ → 1) values are plotted as a function of the hopping anisotropy J 1 /J 2 in Fig. 5a.
The CMF + S analysis with 2D DMRG solver for bosons. In the CMF + S analysis, we consider the N C -site cluster Hamiltoniañ J ij e iqÁðr j Àr i Þb y ibj þ U 2 J ij e iqÁðr j Àr i Þb y i þ H:C: on a triangular-shaped cluster of N C = 10, 15, 21 sites. The mean-field boundary condition on the cluster-edge sites ∂C is implemented by the third term and the twisted frame,b i e ÀiqÁr ib i , is adopted. The cluster Hamiltonian (14) is treated with 2D DMRG solver. Here we take the maximum one-site occupation n max ¼ 4, which is confirmed to be sufficient for the discussion near the ρ = 1 SF-MI (CSF-MI) transition. For large-size clusters and especially for the large Hilbert space of bosons, the exact diagonalization is practically not realistic as a solver for the cluster problem. Therefore, we employ the DMRG on the equivalent 1D chain model with long-range hoppings and mean fields (see Fig. 8a). The DMRG calculation is performed in the standard way but with the mean-field terms in Eq. (14) 14 . The dimension of the truncated matrix product states kept in the present DMRG calculations is typically~10 3 to obtain numerically precise results. To solve in a self-consistent way, we iteratively perform the DMRG calculations until convergence.
It is noteworthy that when we put a real number ψ as an input forH C ðψÞ in the fixed global gauge, the output N À1 C P i2C hb i iH C ðψÞ includes a small but finite imaginary component (≲4% for N C = 21). This is due to a finite-size effect; the order with uniform amplitude ψ and spiral phase twist exp½iq Á r i is not fully commensurate with the shape of the finite-size clusters with the mean-field boundary. We just ignore the small imaginary component in the calculations for each N C , as it decreases with N C and is expected to vanish at the limit of N C → ∞.
The optimization of the spiral twist q is performed in the following way: For different values of q = (q x , 0), the "provisional" critical point U Ã c =jJ 2 j (at which ψ ¼ 0 þ ) is numerically determined (Fig. 8b). The maximum value of U Ã c ðq x Þ=jJ 2 j and the corresponding (q x , 0) were adopted as the CMF + S prediction of the critical point U c /|J 2 | and the BEC momentum Q at the critical point, respectively. The slight variance of Q x from 4π/3 at J 1 /J 2 = 1 (see Fig. 5c) is thought to stem from the same finite cluster-size effect mentioned above.
The CMF + S curves in Fig. 5a are obtained from the size scaling of the phase boundaries for N C = 10, 15, 21. Figure 9 shows the extrapolation of the N C = 10, 15, 21 data to N C → ∞ (λ → 1) for several values of J 1 /J 2 with a linear function of the scaling parameter λ ≡ N B ∕ 3N C 36,37 . Here, N B is the number of NN bonds treated exactly in the cluster (N B = 18, 30, 45 for N C = 10, 15, 21, respectively). The error bars in Fig. 5a are estimated from the variation in the linear fittings for different pairs of the N C = 10, 15, 21 data.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.