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

TL;DR

This paper proposes a method to simulate frustrated quantum antiferromagnetism using ultracold Bose gases at negative temperatures, enabling exploration of quantum phase transitions and phase boundaries in complex lattice systems.

Contribution

It introduces a feasible protocol for analog quantum simulation of frustrated antiferromagnetism at negative temperatures with numerical validation.

Findings

Simulation of quantum phase transitions via negative-temperature states.

Prediction of phase boundary between superfluid and Mott-insulator phases.

Benchmark results for triangular lattice Bose-Hubbard model.

Abstract

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.

Equations

$$\hat{\mathcal{H}}=-\sum _{i,j}{J}_{ij}{\hat{b}}_i^{†}{\hat{b}}_j+\frac{U}{2}\sum _i{\hat{n}}_i({\hat{n}}_i-1)+\frac{V}{a^2}\sum _i|{\mathbf{r}}_i{|}^2{\hat{n}}_i$$

$$\displaystyle{\bf Q}=\left\{\begin{array}[]{ccc}(2\pi/a,0){\equiv{\bf Q}_{\rm M}}&{\rm for}&0\leq J_{1}/J_{2}\leq 0.5,\\ (2\arccos[\frac{-J_{2}}{2J_{1}}]/a,0)&{\rm for}&J_{1}/J_{2}>0.5.\end{array}\right.$$

$${\hat{\mathcal{H}}}_{\mathrm{spin}}=-\sum _{i,j}{J}_{ij}^{\mathrm{XY}}({\hat{S}}_i^x{\hat{S}}_j^x+{\hat{S}}_i^y{\hat{S}}_j^y)+D\sum _i({\hat{S}}_i^z{)}^2$$

$${\mathbf{E}}_{\mathrm{tot}}(\mathbf{r},t)=\sum _{i=1}^3{E}_i\cos ({\mathbf{k}}_i\cdot \mathbf{r}-{\omega }_{\mathrm{L}}t+{\varphi }_i){\mathbf{e}}_z.$$

$$V(\mathbf{r})$$

$$V_{1D}(x)=V^{\prime }{\cos }^2(k_{\mathrm{L}}^{\prime }x+{\varphi }^{\prime })$$

$${\hat{\mathcal{H}}}_i^{\mathrm{GA}}$$

$${\hat{\mathcal{H}}}_i^{\mathrm{GA}}={\epsilon }_{\mathbf{q}}\bar{\psi }\left({\hat{b}}_i^{†}+{\hat{b}}_i-\bar{\psi }\right)-\mu {\hat{n}}_i+\frac{U}{2}{\hat{n}}_i({\hat{n}}_i-1),$$

$$i\frac{\partial }{\partial t}|{\Psi }_i(t)⟩={\hat{\mathcal{H}}}_i^{\mathrm{GA}}|{\Psi }_i(t)⟩.$$

$$i{\dot{f}}_i^{(n)}(t)$$

$$i{\dot{\mathbf{f}}}_i(t)$$

$${\mathbf{f}}_i(t+\Delta t)$$

$${\tilde{\mathcal{H}}}_{\mathrm{C}}$$

$$\bar{\psi }=\frac{1}{N_{\mathrm{C}}}\sum _{i\in \mathrm{C}}⟨{\tilde{b}}_i{⟩}_{{\tilde{\mathcal{H}}}_{\mathrm{C}}(\bar{\psi })}$$

Full text

Frustration is a key concept to understand various emergent phenomena in modern many-body physics diep-04 ; moessner-06 . 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 balents-10 and as a challenge for numerical techniques to handle highly entangled ground states verstraete-08 ; kulagin-13 . Quantum simulation with the use of ultracold atomic gases in optical lattices lewenstein-07 ; bloch-12 ; gross-17 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 damski-05 ; ruostekoski-09 ; eckardt-10 ; huber-10 ; chen-10 ; yamamoto-19 .

One straightforward idea for creating frustration is the use of two-component Fermi gases with (pseudo)spins $\sigma=\uparrow,\downarrow$ greif-13 ; hart-15 ; greif-15 ; parsons-16 ; boll-16 ; cheuk-16 ; mazurenko-17 ; hilker-17 ; brown-17 loaded into a nonbipartite (e.g., triangular becker-10 and kagome jo-12 ) optical lattice. The second-order hopping process provides antiferromagnetic superexchange interactions between the (pseudo)spins $\sigma$ kuklov-03 ; duan-03 , 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 mazurenko-17 , 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 eckardt-05 ; lignier-07 . For ultracold Bose gases with sign-inverted hopping, the relative local phase of Bose-Einstein condensates (BECs) tends to be $\pi$ 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 struck-11 . 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 ramsey-56 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\geq 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\leq 0$. In the pioneer work of Ref. braun-13 , 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.

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 dobrek-99 ; burger-99 ; taie-15 and the statistics of negative absolute temperatures. Our proposal is based on the fact that a negative-temperature state of a system $\hat{\mathcal{H}}$ at $T<0$ realizes the corresponding equilibrium state of the sign-inverted Hamiltonian $-\hat{\mathcal{H}}$ 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 time-dependent Gutzwiller approach (TDGA) to demonstrate that quantum phases of the frustrated Bose-Hubbard model, including chiral superfluid (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 yamamoto-12-2 ; yamamoto-14 ; yamamoto-17 with two-dimensional (2D) density matrix renormalization group (DMRG) solver yamamoto-19 . 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 $\hbar=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:

[TABLE]

with hopping integral $J_{ij}$ (for $i\neq j$), chemical potential $J_{ii}\equiv\mu$, onsite interaction $U$, harmonic trap potential $V|{\bf 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 ${\bf a}_{1}=(a,0)$ direction and $J_{ij}=J_{2}$ in the ${\bf a}_{2}=(-a/2,\sqrt{3}a/2)$ and ${\bf a}_{3}=(-a/2,-\sqrt{3}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}\gg 1$, are reduced to square lattice and 1D chain, respectively.

Gutzwiller analysis. We first discuss the ground state for $V=0$ within the site-decoupling Gutzwiller approach (GA) GA1 ; GA2 ; buonsante09 to get a basic insight into the problem. The ground state in the weak-coupling regime ($U\ll|J_{ij}|$) is well described within the GA under the assumption of the BEC order $\langle b_{i}\rangle\equiv\psi_{i}=\bar{\psi}e^{i{\bf q}\cdot{\bf r}_{i}+\varphi}$ with momentum ${\bf q}$ and global phase $\varphi$. The kinetic energy of Eq. (1) is given by $-\sum_{i,j\neq i}J_{ij}\langle\hat{b}_{i}^{\dagger}\hat{b}_{j}\rangle\approx\varepsilon_{{\bf q}}\bar{\psi}^{2}M$. Here, $\varepsilon_{{\bf q}}\equiv-2(J_{1}\cos{\bf q}\cdot{\bf a}_{1}+J_{2}\cos{\bf q}\cdot{\bf a}_{2}+J_{2}\cos{\bf q}\cdot{\bf 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 ${\bf q}={\bf 0}$, leading to a uniform SF state. The maximum of the kinetic energy is achieved at ${\bf q}=\pm{\bf Q}$ with

[TABLE]

Therefore, if sign-inverted hopping ($J_{1},J_{2}<0$) is prepared, a frustrated CSF state with finite BEC momentum ${\bf q}=\pm{\bf Q}$ is realized. The choice of ${\bf q}={\bf Q}$ or $-{\bf Q}$ represents the degeneracy with respect to the chirality of vortex in unit triangles. Hereafter we suppose ${\bf q}={\bf Q}$ to be spontaneously selected. In the equilateral case ($J_{1}=J_{2}\equiv J$), the momentum is ${\bf Q}=(4\pi/3a,0)\equiv{\bf Q}_{K}$. Therefore, the CSF state forms a “three-color” arrangement of the local phase (${\rm Arg}[\psi_{i}]=0$, $2\pi/3$, and $4\pi/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 ${\bf Q}$. In the parameter range $0\leq J_{1}/J_{2}\leq 0.5$ and the 1D limit $J_{1}/J_{2}\gg 1$, a “two-color” ([math] and $\pi$) pattern is formed with no chiral degeneracy. The $J_{1}/J_{2}$ dependence of ${\bf Q}$ is analogous to the pitch vector of spin spiral states in spatially-anisotropic triangular antiferromagnets spinwave1 ; spinwave2 ; coldea-01 ; zvyagin-19 .

For large, repulsive interaction $U$, lattice bosons undergo a quantum phase transition to the MI state when the filling factor $\rho\equiv M^{-1}\sum_{i}\langle\hat{n}_{i}\rangle$ is an integer greiner-02 . Performing the GA decoupling $\hat{b}_{i}^{\dagger}\hat{b}_{j}\approx\psi_{j}\hat{b}_{i}^{\dagger}+\psi_{i}^{\ast}\hat{b}_{j}-\psi_{i}^{\ast}\psi_{j}$ in Eq. (1), we calculate the order parameter $\bar{\psi}$ for $\rho=1$ in the unfrustrated ($J_{1},J_{2}>0$: ${\bf q}={\bf 0}$) and frustrated ($J_{1},J_{2}<0$: ${\bf q}={\bf Q}$) ground states as a function of $U/J_{2}$ (see Methods). Figures 2b and 2c 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 $\bar{\psi}$ vanishes is given as $U^{({\rm GA})}_{\rm c}/|J|=17.5$ for the frustrated case, which is a half of $U^{({\rm GA})}_{\rm 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 $\rho$, the critical point is given by $U^{({\rm GA})}_{\rm c}=-\varepsilon_{{\bf q}}(\sqrt{\rho}+\sqrt{\rho+1})^{2}$ with ${\bf q}={\bf 0}$ (${\bf q}={\bf Q}$) for the unfrustrated (frustrated) case. The ratio $|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|$ equals to 1 only at $J_{1}/J_{2}=0$ or $J_{1}/J_{2}\rightarrow\infty$, indicating that the reduction effect due to frustration exists even in the “two-color” region of $0<J_{1}/J_{2}\leq 0.5$.

It is noted that the curves of $\bar{\psi}$ 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 $|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|$ (specifically, the factor 1/2 when $J_{1}=J_{2}$). This is the case for general values of $J_{1}/J_{2}$ and $\rho$. 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 $|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|U$ are identical within the GA, except for the spatial phase distribution ${\rm Arg}[\psi_{i}]={\bf Q}\cdot{\bf r}_{i}$. This is because the scale of the kinetic energy for the frustrated system is reduced by the factor $|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{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 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 buonsante09 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 ${\rm Arg}[\psi_{i}]$ is uniform in the unfrustrated SF state.

First, one needs to introduce the spatial phase distribution ${\rm Arg}[\psi_{i}]={\bf Q}\cdot{\bf r}_{i}$ to create the frustrated CSF state. To this end, here we suggest the use of the phase-imprinting techniques dobrek-99 ; burger-99 ; taie-15 . When a single-particle energy difference $\delta E$ is introduced between two sites, the relative phase on the two sites starts to develop with $\exp[i\delta Et]$ in time $t$. In a region deep inside the SF phase, the kinetic energy ($\propto\varepsilon_{{\bf q}}$) reaches the maximum when the BEC momentum takes ${{\bf q}}={{\bf Q}}$ given in Eq. (4). One can transfer ${\bf q}$ from ${\bf 0}$ to ${\bf Q}$ by introducing a temporary linear gradient potential $\hat{V}^{\rm ext}=\delta E\sum_{i}(x_{i}/a)\hat{n}_{i}$ for appropriate time $\delta t=Q_{x}a/\delta E$ [see Fig. 3b]. Such a temporary potential can be created, e.g., by a magnetic field gradient or by an extra 1D optical lattice satisfying ${\bf Q}={\bf Q}_{\rm M}$ or ${\bf Q}={\bf Q}_{\rm K}$ (see Methods). One has to perform the phase-imprinting operation in a much shorter time than the time scale of the hopping ($\delta t\ll|J_{ij}|^{-1}$) in order to affect only the local phases. Besides, the local energy difference $\delta E$ has to be set to a large enough value compared to $|U|$ and $|V|$ so that one can safely avoid the influence of the inhomogeneity of the density profile (see Methods). Since such well-controlled phase imprinting has been successfully made in previous experiments burger-99 ; taie-15 , 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^{i{\bf Q}\cdot{\bf r}_{i}}$ should of course be dynamically unstable, since 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, in order to realize a stable negative-temperature ground state (at $T\approx-0$) for $J_{ij}>0$ and $U,V<0$. The created negative-temperature state should be equivalent to the ground state (at $|T|\approx+0$) of the frustrated system with $J_{ij}<0$ and $U,V>0$ since the physics of the two systems obey the same factor $\exp[-\hat{\mathcal{H}}/k_{\rm B}T]$. In Ref. braun-13 for square (unfrustrated) lattice, this has been performed simply with $U\rightarrow-U$ and $V\rightarrow-V$ by using the Feshbach resonance chin-10 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 order-parameter ($|\psi_{i}|$) distributions shown in Fig. 3a after the phase imprinting ${\rm Arg}[\psi_{i}]={\bf Q}\cdot{\bf r}_{i}$ is expected to correspond to the ground state of the frustrated Hamiltonian with the interactions rescaled by the factor $|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|$, as explained above in the homogeneous case. Therefore, the values of $U$ and $V$ have to be changed as $U\rightarrow-|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|U$ and $V\rightarrow-|\varepsilon_{\bf{Q}}/\varepsilon_{\bf{0}}|V$ (e.g., $U\rightarrow-U/2$ and $V\rightarrow-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.

TDGA simulation.— In the framework of the TDGA damski-03 ; zakrzewski-05 , 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 $\sum_{n}f_{i}^{(n)}|n\rangle\rightarrow\sum_{n}e^{in{\bf Q}_{\rm K}\cdot{\bf r}_{i}}f_{i}^{(n)}|n\rangle$ on the local wave function at every site $i$. Here, $f_{i}^{(n)}$ is the coefficient of the Fock basis $|n\rangle$. In addition, we perform the sudden changes of $U=U_{0}\rightarrow-U_{0}/2$ and $V=0.001U_{0}\rightarrow-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<t\leq 200U_{0}^{-1}$ 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<t\leq 200U_{0}^{-1}$. 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\lesssim 2U_{0}^{-1}$ 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 $\delta{n}^{2}\equiv\overline{\langle\hat{n}_{i}^{2}\rangle-\langle\hat{n}_{i}\rangle^{2}}$ averaged over the center sites within $|{\bf r}_{i}|\leq 10a$. As shown in Fig. 3f, the value of $\delta n^{2}$ rapidly increases and then exhibits an irregular oscillation in the case without the phase imprinting, while 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 negative-temperature CSF state shown in Fig. 3c, we slowly increase the value of $|U/J|$ for $t>200U_{0}^{-1}$ to observe the CSF-to-MI transition. In experiments, the tuning of $|U/J|$ is performed by controlling the height of the optical lattice greiner-02 . Here, we assume that $|U/J|$ increases simply through $J$ decreasing linearly with $J=0.08U_{0}-0.0001(tU_{0}-200)U_{0}$ for $t>200U_{0}^{-1}$. As shown in Fig. 4a, when $|U/J|$ increases, the MI plateau is gradually formed and $|\psi_{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 (Figs. 4b and 4c).

To show the transition process more clearly, we plot in Fig. 4d the time evolution of the scaled value of the density fluctuation $\delta{n}^{2}(t)/\delta{n}^{2}(0)$. There we see that $\delta{n}^{2}(t)/\delta{n}^{2}(0)$ decreases with $|U/J|$ and almost vanishes at $|U/J|\approx|U^{({\rm GA})}_{\rm 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 gaunt-13 ; chomaz-15 , which is modeled by $V=0$ and the open boundary at $|{\bf 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\approx 0$) and strong-coupling ($U\rightarrow\infty$) limits, and is thus expected to reasonably interpolate the two limits in two or higher dimensions GA1 ; zwerger-03 ; lewenstein-07 . In fact, the TDGA simulation has been often used to describe the dynamics of Bose-Hubbard systems altman-05 ; snoek-07 ; snoek-12 ; saito-12 , 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) mun-07 . It is therefore expected that the above TDGA simulation gave a proper description for the dynamical stability of the negative-temperature 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, since 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 yamamoto-19 , to the present system. We consider the $N_{\rm 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, $\tilde{b}_{i}\equiv e^{-i{\bf Q}\cdot{\bf r}_{i}}\hat{b}_{i}$, with optimizing ${\bf Q}=(Q_{x},0)$ (see Methods). The CMF+S method permits the systematic inclusion of quantum intersite correlations by increasing $N_{\rm C}$, which connects between the GA ($N_{\rm C}=1$) and exactly quantum ($N_{\rm C}\rightarrow\infty$) results.

Figure 5a summarizes the $N_{\rm C}=10,15,21$ data and the CMF+S ($N_{\rm C}\rightarrow\infty$) result for the $\rho=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_{\rm c}/|J_{2}|$ for the frustrated case exhibits a nonmonotonic behavior with a dip around $J_{1}/J_{2}\approx 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_{\rm c}$ is strongly reduced from the GA prediction $U_{\rm c}^{({\rm GA})}$ due to the inclusion of the intersite quantum correlations. It should be noted that the relative difference $(U_{\rm c}^{({\rm GA})}-U_{\rm c})/U_{\rm c}^{({\rm GA})}$ between the GA and CMF+S values is much larger for the frustrated case ($\sim 40$-$50\%$) than the unfrustrated case ($\sim 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\pi/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 ${\bf Q}$ in the CSF state and the transition to the MI can be simply detected by time-of-flight (TOF) images of momentum distributions greiner-02 . A more precise determination of the critical point can be made by extracting the condensate fraction from the TOF images spielman-08 , by observing the critical velocity using a moving optical lattice mun-07 , or by measuring the density fluctuation $\delta{n}^{2}$ using the quantum-gas microscope mazurenko-17 ; parsons-16 ; boll-16 ; cheuk-16 ; hilker-17 ; brown-17 . The frustrated CSF-MI transition is easily distinguishable from the standard SF-MI transition thanks to the sharp difference in the value of $U_{\rm 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 superfluid 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-of-the-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 altman-02 : $\hat{b}_{i}\rightarrow\sqrt{\frac{\rho}{2}}\hat{S}_{i}^{-},~{}\hat{n}_{i}\rightarrow\hat{S}_{i}^{z}+\rho$, which is valid in the vicinity of the transition between the CSF and MI phases at integer fillings $\rho$. The Hamiltonian (1) is mapped onto the spin-1 XY model

[TABLE]

with the XY spin exchange $J^{\rm XY}_{ij}=-\rho 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^{\circ}$” magnetic order (corresponding to CSF) to nonmagnetic state (MI) in spin-1 easy-plane triangular antiferromagnets such as CsFeCl3 kurita-16 ; hayashida-18 , which has recently attracted considerable attention in the connection with the novel excitations near quantum criticality hayashida-19 . From the standpoint of fundamental statistical physics, the CSF-to-MI transition should be associated with the spontaneous U(1)$\times\mathbb{Z}_{2}$ symmetry breaking with respect to the global phase and the chirality determined by ${\bf q}={\bf Q}$ or $-{\bf 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, while the real materials have strong 3D couplings between the triangular layers kurita-16 ; hayashida-18 .

Moreover, it has been expected that adding long-range interatomic interaction to the present system may give rise to an exotic chiral MI state zaletel-14 in between possible separate U(1) and $\mathbb{Z}_{2}$ symmetry breakings. This is essentially equivalent to the so-called “chiral liquid” expected in a spin-1 frustrated magnet wang-17 . 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 $\rho={\rm half}~{}{\rm integers}$ balents-10 ; eckardt-10 . 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 ${\bf k}_{1}$ = $k_{\rm L}(1,0)$, ${\bf k}_{2}$ = $k_{\rm L}(-1/2,-\sqrt{3}/2)$, and ${\bf k}_{3}$ = $k_{\rm L}(-1/2,\sqrt{3}/2)$ and equal frequency $\omega_{\rm L}$ becker-10 . 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

[TABLE]

The dipole potential generated by the electric field is proportional to its squared amplitude,

[TABLE]

where ${\bf b}_{1}={\bf k}_{1}-{\bf k}_{2}$, ${\bf b}_{2}={\bf k}_{3}-{\bf k}_{2}$, $\phi_{ij}=\phi_{i}-\phi_{j}$, and $A(t)$ represents the terms dependent on time $t$. Since the frequency of light is quite large, only the time-averaged value of $|{\bf E}_{\rm tot}|^{2}$ can affect atoms. All the terms in $A(t)$ oscillate at frequency $2\omega_{\rm L}$, and thus can be dropped. Finally, we obtain a periodic dipole potential

[TABLE]

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 $\phi_{i}$ of the laser beams yield only a global shift of the lattice in position. The primitive lattice vectors ${\bf a}_{1}$ and ${\bf a}_{2}$ are given so that ${\bf a}_{i}\cdot{\bf b}_{j}=2\pi\delta_{ij}$. We define ${\bf a}_{3}=-{\bf a}_{1}-{\bf 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=|{\bf a}_{i}|=4\pi/3k_{\rm L}=2\lambda_{\rm L}/3$ with $\lambda_{\rm 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. 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 $\delta 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 ${\bf Q}={\bf Q}_{\rm M}\equiv(2\pi/a,0)$ (two-color) and ${\bf Q}={\bf Q}_{\rm K}\equiv(4\pi/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 ${\bf a}_{1}$ direction:

[TABLE]

with amplitude $V^{\prime}$, wave vector $k_{\rm L}^{\prime}=2\pi/\lambda_{\rm L}^{\prime}$, and phase $\phi^{\prime}$. Here, let us set $\phi_{i}=0$ ($i=1,2,3$) for the triangular optical lattice without loss of generality.

Let us first consider the range of $0\leq J_{1}/J_{2}\lesssim 0.5$ ($V_{1}\gg 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 ${\bf Q}={\bf Q}_{\rm 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 $\delta E$ only between the two-sublattice groups of sites for time $\delta t=2\pi/\delta E$. This can be achieved by using an additional 1D optical lattice of magic wavelength defined by $\lambda_{\rm L}^{\prime}=4\lambda_{\rm L}/3$ and phase shift $\phi^{\prime}=0$ or $\pi/2$ (mod $\pi$). Figures 6a and 6b show an example of the total potential $V({\bf r})+V_{\rm 1D}(x)$ with the parameters $2V_{1}=V_{2}=V_{3}=V^{\prime}$ and $\phi^{\prime}=0$. Note that the two options in $\phi^{\prime}$ correspond to the exchange of A and B.

In a similar way, using a 1D periodic potential with $\lambda_{\rm L}^{\prime}=2\lambda_{\rm L}$, one can also implement a temporary energy difference $\delta E$ between the three-sublattice groups of sites, say A, B, and C. The additional lasers are shined for time $\delta t=4\pi/3\delta E$ to imprint the three-color phase configuration illustrated in the middle panel of Fig. 2a to the initial SF state. Figures 6c and 6d show an example of the total potential $V({\bf r})+V_{\rm 1D}(x)$ with the parameters $V_{1}=V_{2}=V_{3}=V^{\prime}$ and $\phi^{\prime}=\pi/12$. Note that the phase shift $\phi^{\prime}$ has six options, $(2n-1)\pi/12$ $(n=1,2,\cdots,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

[TABLE]

as the result of the decoupling $\hat{b}_{i}^{\dagger}\hat{b}_{j}\approx\psi_{j}\hat{b}_{i}^{\dagger}+\psi_{i}^{\ast}\hat{b}_{j}-\psi_{i}^{\ast}\psi_{j}$ in the original Hamiltonian (1). The results displayed in Fig. 2b are calculated under the assumption of finite-momentum BEC, $\psi_{i}=\bar{\psi}e^{i{\bf q}\cdot{\bf r}_{i}+\varphi}$, for $V=0$. In this case, one has only to consider a certain single site $i$, e.g., the site at the origin ${\bf r}_{i}={\bf 0}$, and $\varphi=0$ without loss of generality. Equation (10) becomes

[TABLE]

for the origin site $i$. The minimization of the kinetic energy gives ${\bf q}={\bf 0}$ for $J_{1},J_{2}>0$ and ${\bf q}={\bf Q}$ given in Eq. (4) for $J_{1},J_{2}<0$. The Hamiltonian $\hat{\mathcal{H}}_{i}^{\rm GA}$ can be easily diagonalized on the Fock state basis for the local wave function, $|\Psi_{i}\rangle\equiv\sum_{n=0}^{n_{\rm max}}f_{i}^{(n)}|n\rangle$, in which the maximum one-site occupation number $n_{\rm max}$ must be sufficiently large (we take $n_{\rm max}=10$). The order parameter is obtained from $\bar{\psi}=\sum_{n}\sqrt{n}f_{i}^{(n-1)\ast}f_{i}^{(n)}$ with eigenvector $f_{i}^{(n)}$ in a self-consistent way.

Note that the effective one-body Hamiltonians $\hat{\mathcal{H}}_{i}^{\rm GA}$ in the unfrustrated and frustrated cases differ only in $\varepsilon_{{\bf q}}$ ($\varepsilon_{{\bf 0}}$ or $\varepsilon_{{\bf Q}}$). Therefore, if the values of all the other terms are multiplied by $|\varepsilon_{{\bf Q}}/\varepsilon_{{\bf 0}}|$, the results of the GA calculations, such as the value of $|\psi_{i}|$, for the unfrustrated system become exactly same as those for the frustrated system, except for the chiral phase distribution $e^{i{\bf Q}\cdot{\bf r}_{i}}$ and the overall energy scale (which is also multiplied by $|\varepsilon_{{\bf Q}}/\varepsilon_{{\bf 0}}|$).

In the presence of the trap potential $V\neq 0$, the mean field $\psi_{i}$ can no longer be assumed to have spatially uniform amplitude. Therefore, one has to deal with Eq. (10) on the entire lattice sites, each of which is connected to the six neighboring sites through the mean fields $\{\psi_{i\pm{\bf a}_{n}};n=1,2,3\}$. To prepare the initial state shown in Fig. 3a, we solve the set of self-consistent equations $\psi_{i}=\sum_{n}\sqrt{n}f_{i}^{(n-1)\ast}f_{i}^{(n)}$ for all sites within a cutoff length $l_{\rm c}$ from the trap center. We take $l_{\rm c}=36a$. The chemical potential $\mu$ is determined from the global number equation $\sum_{i}\langle\hat{n}_{i}\rangle=N$, which must be solved simultaneously with the equations for $\psi_{i}$. To efficiently achieve the convergence in the self-consistent calculations, we take the uniform solution for $\psi_{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

[TABLE]

From Eq. (10) with $|\Psi_{i}(t)\rangle=\sum_{n}f_{i}^{(n)}(t)|n\rangle$, Eq. (12) becomes

[TABLE]

with $\psi_{j}(t)=\langle\Psi_{j}(t)|\hat{b}_{j}|\Psi_{j}(t)\rangle=\sum_{n}\sqrt{n}{f}_{j}^{(n-1)\ast}{f}_{j}^{(n)}(t)$.

We numerically solve the set of equations of motion for $f_{i}^{(n)}(t)$ using the Crank-Nicolson method javanainen-99 . We first rewrite Eq. (13) in a matrix form:

[TABLE]

where ${\bf f}_{i}$ is the vector with the components $f_{i}^{(n)}$ ($0\leq n\leq n_{\rm max}$) and ${\bf H}^{\rm GA}_{i}$ is the $(n_{\rm max}+1)\times(n_{\rm max}+1)$ matrix form of the GA local Hamiltonian, which is time-dependent via the order parameters $\psi_{j}(t)$ of neighboring sites of site $i$. The explicit components of ${\bf H}^{\rm GA}_{i}$ are easily obtained from Eq. (13). The value of ${{\bf f}}_{i}$ after time evolution in a short time $\Delta t$ is given by

[TABLE]

where ${\bf I}$ is the identity matrix. Here $\Delta 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 $\Delta t$. In order to calculate ${\bf f}_{i}(t+\Delta t)$ at site $i$, the values of $\psi_{j}(t)$ on its neighboring sites are required. Therefore, one has to calculate $\psi_{i}(t)=\sum_{n}\sqrt{n}{f}_{j}^{(n-1)\ast}{f}_{j}^{(n)}(t)$ and update ${\bf f}_{i}(t)$ according to Eq. (15) in parallel for all sites (within $|{\bf r}_{i}|\leq l_{\rm c}$). The local particle number at time $t$ can be calculated by $n_{i}(t)=\sum_{n}n|f_{i}^{(n)}(t)|^{2}$.

The phase imprinting operation. In the section “TDGA simulation,” we assume that the phase imprinting can be theoretically implemented by the operation $\sum_{n}f_{i}^{(n)}|n\rangle\rightarrow\sum_{n}e^{in{\bf Q}_{\rm K}\cdot{\bf r}_{i}}f_{i}^{(n)}|n\rangle$ 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 $\hat{V}^{\rm ext}=\delta E\sum_{i}(x_{i}/a)\hat{n}_{i}$. As seen in Fig. 7a, the chiral structure in the local phase, ${\rm Arg}[\psi_{i}]={\bf Q}\cdot{\bf r}_{i}$, could be successfully imprinted by applying the temporary potential for $\delta t=Q_{x}a/\delta E$ ($=4\pi/3\delta 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 $\delta E$ exceeds $\sim 10U_{0}$. The corresponding imprinting time $\delta t\lesssim 0.4U_{0}^{-1}$ is much shorter than the typical time scale of the experiments on the SF(CSF)-MI transition (see Fig. 4).

The CMF+S analysis with 2D DMRG solver for bosons. In the CMF+S analysis, we consider the $N_{\rm C}$-site cluster Hamiltonian

[TABLE]

on a triangular-shaped cluster of $N_{\rm C}=10,15,21$ sites. The mean-field boundary condition on the cluster-edge sites $\partial{\rm C}$ is implemented by the third term, and the twisted frame, $\tilde{b}_{i}\equiv e^{-i{\bf q}\cdot{\bf r}_{i}}\hat{b}_{i}$, is adopted. The cluster Hamiltonian (16) is treated with 2D DMRG solver. Here we take the maximum one-site occupation $n_{\rm max}=4$, which is confirmed to be sufficient for the discussion near the $\rho=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. (16) yamamoto-19 . The dimension of the truncated matrix product states kept in the present DMRG calculations is typically $\sim 10^{3}$ to obtain numerically precise results. In order to solve

[TABLE]

in a self-consistent way, we iteratively perform the DMRG calculations until convergence.

Note that when we put a real number $\bar{\psi}$ as an input for $\tilde{\mathcal{H}}_{\rm C}(\bar{\psi})$ in the fixed global gauge, the output $N_{\rm C}^{-1}\sum_{i\in{\rm C}}\langle\tilde{b}_{i}\rangle_{\tilde{\mathcal{H}}_{\rm C}(\bar{\psi})}$ includes a small but finite imaginary component ($\lesssim 4\%$ for $N_{\rm C}=21$). This is due to a finite-size effect; the order with uniform amplitude $\bar{\psi}$ and spiral phase twist $\exp[{i{\bf q}\cdot{\bf 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_{\rm C}$ since it decreases with $N_{\rm C}$ and is expected to vanish at the limit of $N_{\rm C}\rightarrow\infty$.

The optimization of the spiral twist ${\bf q}$ is performed in the following way: For different values of ${\bf q}=(q_{x},0)$, the “provisional” critical point $U^{\ast}_{\rm c}/|J_{2}|$ (at which $\bar{\psi}=0^{+}$) is numerically determined (Fig. 8b). The maximum value of $U^{\ast}_{\rm c}(q_{x})/|J_{2}|$ and the corresponding $(q_{x},0)$ were adopted as the CMF+S prediction of the critical point $U_{\rm c}/|J_{2}|$ and the BEC momentum ${\bf Q}$ at the critical point, respectively. The slight variance of $Q_{x}$ from $4\pi/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_{\rm C}=10,15,21$. Figure 9 shows the extrapolation of the $N_{\rm C}=10,15,21$ data to $N_{\rm C}\rightarrow\infty$ ($\lambda\rightarrow 1$) for several values of $J_{1}/J_{2}$ with a linear function of the scaling parameter $\lambda\equiv N_{\rm B}/3N_{\rm C}$ yamamoto-12-2 ; yamamoto-14 . Here, $N_{\rm B}$ is the number of NN bonds treated exactly in the cluster ($N_{\rm B}=18,30,45$ for $N_{\rm 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_{\rm C}=10,15,21$ data.