Frustrated Quantum Magnetism with Bose Gases in Triangular Optical Lattices at Negative Absolute Temperatures

We propose an experimental protocol to perform analog quantum simulation of frustrated antiferromagnetism with strong quantum fluctuations by using ultracold Bose gases in triangular optical lattices with isosceles anisotoropy of hoppings. Specifically, we combine a phase-imprinting scheme with sudden sign inversion of the interatomic interaction and the trap potential to prepare a chiral superfluid state with a negative absolute temperature. In the framework of the time-dependent Gutzwiller approach, we compute the time evolution of the state subjected to a slow sweep of the hopping energy. We show that in this process the system simulates a state near equilibrium of the Bose-Hubbard model with sign-inverted hoppings. By means of the cluster mean-field method with a cluster size scaling, we quantitatively predict the quantum critical point of the superfluid-Mott insulator transition as a function of the spatial anisotropy parameter, which serves as a benchmark for quantum simulation.

Introduction.-Frustration is a key concept to understand various emergent phenomena in the modern manybody 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, the effort in realizing frustrated quantum systems with cold atoms has not yet paid off even though 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 [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 neighbor-ing 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 not easy again to reach a quantum regime of low temperature and density because the lattice shaking can be a source of heating.
In this Letter, we propose and examine a realistic route to create frustrated Bose gases in a quantum regime using the combination of the phase-imprinting techniques [31][32][33] and the statistics of negative absolute temperatures [34]. A state in thermal equilibrium is usually described by a statistical ensemble in which the lower-energy states are more occupied than higherenergy 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 the pioneer work of Ref. [35], Braun et al. 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.
Our present proposal is based on the fact that a negative-temperature state of a systemĤ at T < 0 realizes the corresponding equilibrium state of the signinverted 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 to the lattice shaking.
Supposing Bose gases in a triangular lattice, we simulate the dynamics along the protocol within the timedependent Gutzwiller approach (TDGA) to demonstrate that quantum phases of the frustrated Bose-Hubbard model, including chiral superfluid (CSF), could indeed be realized. Moreover, considering general cases with isosceles-type anisotropy of hoppings, 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][37][38] with two-dimensional (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. See also the Supplementary Material (SM) [39] for the technical details of the calculations.
The Bose-Hubbard model on triangular lattice.-A system of Bose gases in a deep optical lattice is described by the Bose-Hubbard model: 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 hoppings of isosceles type, which can be created by tuning the intensity of one of the three lasers to be different from the others [39]. The spatial anisotropy 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, √ 3a/2) and a 3 = (−a/2, − √ 3a/2) directions [Figs. 1(a) and 1(b)].
We first discuss the ground state for V = 0 within the site-decoupling Gutzwiller approach (GA) [40,41]. Under the assumption of the BEC order b i ≡ ψ i = ψe iq·ri+ϕ with momentum q and global phase ϕ, the kinetic energy of Eq. (1) is given by 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. 1(b). 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 "twocolor" (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].
For large, repulsive interaction U , lattice bosons undergo a quantum phase transition to the MI state when the filling factor ρ ≡ 1/M i n i is an integer [45]. Per- , 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 [see of the anisotropy J 1 /J 2 and ρ, the critical point is given by U 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.
Negative absolute temperature.-To experimentally create the frustrated quantum states, it is required to prepare sign-inverted hoppings J i,j < 0 with avoiding serious heating of the system. Below, we explain the details of the protocol through 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. 2(a) 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). First, we suddenly change U and V to be attractive (U < 0) and anti-trapping (V < 0), and at the same time, achieve the maximum interaction and potential energies. In Ref. [35] for square (unfrustrated) lattice, this has been performed simply with U → −U and V → −V by using the Feshbach resonance [46] and blue-detuned lasers. In the present case with frustration, one requires special attention to achieve the maximum interaction and potential energies; 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 ) in order to adjust the energy scale in consideration of the reduction of the kinetic energy due to the frustration.
By making also the kinetic energy reach its maximum, one can realize a negative-temperature ground state (at T ≈ −0) for J ij > 0 and U, V < 0, which should be equivalent to the ground state (at |T | ≈ +0) of the system with J ij < 0 and U, V > 0 since the physics of the two systems obey the same factor exp[−Ĥ/k B T ]. To this end, here we suggest the use of the phase-imprinting techniques [31][32][33]. 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 potentialV ext = δE i (x i /a)n i for appropriate time δt = Q x a/δE [see Fig. 2(b)]. Such a temporary potential can be created, e.g., by a magnetic field gradient or by an extra 1D optical lattice satisfying Q = Q M or TDGA simulation.-In the framework of the TDGA [39, 47,48], we simulate the time evolution of the initial state in Fig. 2(a) after suddenly performing the three operations shown in Fig. 2(b). In the simulation, we implement the phase imprinting on the initial state by the operation n f As shown in Fig. 2(c), the created negativetemperature state is predicted to be indeed dynamically stable for 0 < t ≤ 200U −1 0 . Note that the state collapses with increasing |U/J|, the MI plateau is gradually formed and |ψ i | in the trap center decreases and eventually almost vanishes. Until the transition point, the three-color phase profile in the CSF state is properly kept within a global phase shift [ Fig. 2(d)]. To show the transition process more clearly, we plot in Fig. 3(a) the time evolution of the density fluctuation δn 2 ≡ n 2 i − n i 2 averaged over the center sites within |r i | ≤ 10a. The value of δn 2 (t)/δn 2 (0) decreases as |U/J| and almost vanishes at |U/J| ≈ |U (GA) c /J| = 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 [49,50], 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. 3(a) for N = 4692). We also plot in Fig. 3(b) the unfrustrated case of the standard SF-MI transition as a reference, which shows e sharp difference from the frustrated case in the values of |U/J| of the transition region.
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 inset of Fig. 4(a) for N C = 10] under the mean-field boundary condition. We work in the twisted frame,b i ≡ e iQ·rib i , with optimizing Q = (Q x , 0) [39]. 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 4(a) 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, while 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 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. 4(b). This indicates that the quantum effects are strongly enhanced by the interplay with frustration. Figure 4(c) 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 [39].
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 [45]. A more precise determination of the critical point can be made by extracting the condensate fraction from the TOF images [52], by observing the critical velocity using a moving optical lattice [53], 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 |.
Conclusions.-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. Simulatng the time evolution, we demonstrated that a dynamically-stable superfluid state with chiral symmetry breaking was indeed realized, and underwent the quantum transition to the MI as the hopping amplitude decreased. Moreover, we performed state-of-the-art numerical calculations on the quantum critical point as a function of the spatial hopping anisotoropy, which predicted a significant interplay of frustration and quantum fluctuations.
It has been expected that adding long-range interatomic interaction to the present system may give rise to an exotic chiral MI state [54], which is essentially equivalent to the so-called "chiral liquid" expected in a spin-1 frustrated magnet [55]. 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. In this Supplementary Material, we provide more details on the preparation of triangular optical lattice with spatially-anisotropic hoppings, the phase imprinting with an extra 1D optical lattice, and the calculation techniques used in the main text.

Triangular optical lattice with spatially-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, − √ 3/2), and k 3 = k L (−1/2, √ 3/2) and equal frequency ω L [1]. 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 (S1) The dipole potential generated by the electric field is proportional to its squared amplitude, , and A(t) represents the terms dependent on time t. Since 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 . In the main text, 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 E 1 = 1.6E 2 = 1.6E 3 (1.6V 1 = V 2 = V 3 ) yields the potential landscape shown in Fig.  1(a), 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 in the main text, one has to introduce a singleparticle 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 , wave vector k L = 2π/λ L , and phase φ . 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 Fig. 1(b). 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 two-sublattice groups of sites for time δt = 2π/δE. This can be achieved by using an additional 1D optical lattice of magic wavelength defined by λ L = 4λ L /3 and phase shift φ = 0 or π/2 (mod π). Figure S1(a) shows an example of the total potential V (r) + V 1D (x) with the parameteres 2V 1 = V 2 = V 3 = V and φ = 0. Note that the two options in φ correspond to the exchange of A and B.
In a similar way, using a 1D periodic potential with λ 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 Fig. 1(b) to the initial SF state. Figure S1(b) shows an example of the total potential V (r) + V 1D (x) with the parameteres V 1 = V 2 = V 3 = V and φ = π/12. Note that the phase shift φ has six options, (2n − 1)π/12 (n = 1, 2, · · · , 6), reflecting the possible permutation of A, B, and C.  Fig. 3(a).

The TDGA equation in a trap potential
The TDGA equation is given by for each site i. We numetically solve the equations for all sites (within |r i | ≤ l c ) in parallel. In Figs. 2(c) and 3(a), we show the time evolution of the initial state in Fig. 2(a) after the sudden changes U → −U/2 and V → −V /2 at t = 0 as well as the phase-imprinting Fig. S2, as a reference for the comparison, we present the results in the case with the same settings but without the phase-imprinting operation. It can be seen in the case without the phase-imprinting that the order parameter decreases within a short time, although the density profile is kept. The density fluctuation δn 2 rapidly increases and then exhibits an irregular oscillation. These results indicate that the state immediately collapses due to the dynamical instability, which is attributed to the fact that the conditions for the negative-temperature ground state -the maximum kinetic, interaction, and potential energies -are not satisfied.

The CMF+S analysis with 2D DMRG solver for bosons
In the CMF+S analysis, we consider the N C -site cluster Hamiltoniañ 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·rib i , is adopted. The cluster Hamiltonian (S8) 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. S3). The DMRG calculation is performed in the standard way but with the mean-field terms in Eq. (S8) [2]. The dimension of the truncated matrix product states kept in the present DMRG calculations is typically ∼ 10 3 to obtain numerically precise results. In order to solvē in a self-consistent way, we iteratively perform the DMRG calculations until convergence.
Note that when we put a real numberψ as an input forH C (ψ) in the fixed global gauge, the output (1/N C ) i∈C b i H 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 since it decreases as 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 /|J 2 | (at whichψ = 0 + ) is numerically determined [Fig. S4]. The maximum value of U * c (q x )/|J 2 | and (q x , 0) at which it occurs 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. 4(c)] is thought to stem from the same finite cluster-size effect mentioned above.
The CMF+S curves in Fig. 4(a) are obtained from the size scaling of the phase boundaries for N C = 10, 15, 21. Figure (S5) 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 [3,4]. 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 estimated from the variation in the linear fittings for different pairs of the N C = 10, 15, 21 data are smaller than the symbol size in Fig. 4(a).