Real time dynamics and proposal for feasible experiments of lattice gauge-Higgs model simulated by cold atoms

Lattice gauge theory has provided a crucial non-perturbative method in studying canonical models in high-energy physics such as quantum chromodynamics. Among other models of lattice gauge theory, the lattice gauge-Higgs model is a quite important one because it describes wide variety of phenomena/models related to the Anderson-Higgs mechanism such as superconductivity, the standard model of particle physics, and inflation process of the early universe. In this paper, we first show that atomic description of the lattice gauge model allows us to explore real time dynamics of the gauge variables by using the Gross-Pitaevskii equations. Numerical simulations of the time development of an electric flux reveal some interesting characteristics of dynamical aspect of the model and determine its phase diagram. Next, to realize a quantum simulator of the U(1) lattice gauge-Higgs model on an optical lattice filled by cold atoms, we propose two feasible methods: (i) Wannier states in the excited bands and (ii) dipolar atoms in a multilayer optical lattice. We pay attentions to respect the constraint of Gauss's law and avoid nonlocal gauge interactions.


I. INTRODUCTION
Cold atoms in an optical lattice have been used as versatile quantum simulators for various many-body quantum systems and some important results were obtained [1]. Recently, there appeared several proposals to simulate models of lattice gauge theory (LGT) [2][3][4]. Since its introduction, LGT has been an indispensable tool to studying the non-perturbative aspect of quantum models in the high-energy physics (HEP) such as confinement of quarks, the spontaneous chiral-symmetry breaking, etc. Atomic simulations, if realized, shall certainly clarify the dynamics, i.e, the time evolution of lattice gauge models, which is far beyond the present theoretical standard.
At present, the proposals are classified into two approaches. In the first approach [5][6][7][8][9][10][11], atoms with spin degrees of freedom are put on links of the optical lattice. Such a lattice gauge model is called the "quantum link model" or "gauge magnet" [12][13][14]. Although the Gauss's law (divergence of the electric field is just the charge density of matter fields) [15] is assured as the conservation of "angular momentum" [16], the Hilbert space of this model itself truncates the full Hilbert space of the original U (1) LGT studied in HEP.
In the second approach [17], one considers the Bose-Einstein condensate (BEC) of atoms put on each link of the two or three dimensional optical lattice. The U(1) phase variable of the complex amplitude of BEC plays a role of dynamical gauge field on the links, and its density fluctuation corresponds to the electric field, the conjugate variable of the gauge field [17,18]. Adopting the phase variable of the atomic field as the gauge field assures us that we deal with a U(1) field as in HEP in contrast with the first approach. However, to keep the Gauss's law and the short-range gauge interactions simultaneously, a complex design of the system and also a fine tuning of interaction parameters are generally needed in the experimental setups [17,18].
Recently, we proposed a perspective to overcome the difficulty to respect the local gauge invariance in the second approach [19]. Namely, the cold atomic simulator of the LGT without exact local gauge invariance due to untuned interaction parameters (the Gauss's law is not satisfied) can be a simulator for a lattice "gauge-Higgs" (GH) model [20] with exact local gauge invariance, where the unwelcome interactions that violate the Gauss's law are viewed as gauge-invariant couplings of the gauge field to a Higgs field.
Needless to say, the GH model is a canonical model of Anderson-Higgs mechanism and plays a very important role in various fields of modern physics. Its list includes mass generation in the standard model of HEP, phase transition and the vortex dynamics [21] in superconductivity, time-evolution of the early universe crossing a Higgs phase transition [22,23]. Atomic quantum simulation of this model is certainly welcome because it simulates the real time evolution of the above exciting phenomena.
In what follows, we consider the second approach and propose two realistic experimental setups to the atomic quantum simulator of the U(1) lattice GH model by using cold atoms in a two-dimensional (2D) optical lattice. Our target GH model has a nontrivial phase structure, i.e., existence of the phase boundary between confinement and Higgs phases, and this phase boundary is to be observed by cold-atom experiments. As a reference to such experiments, we make numerical simulations of the time-dependent Gross-Pitaevskii (GP) equation and observe the real-time dynamics of the atomic simulators. In particular, we study the dynamical stability of a single electric flux connecting two charges with opposite signs, corresponding to a density hump and dip for the atomic simulators. The obtained phase boundary is discussed and compared with that of the Monte-Carlo simulations.
which describes the bosons in a single band of a 2D optical lattice. The bosonic atomic fieldsψ a = exp(iθ a ) √ρ a are put on the site a of the square optical lattice. The summation is taken over the unit cell k (yellow region in Fig. 1) and, in each unit cell, over the the site a(b) ∈ 1-6. We confine ourselves to the nearest-neighbor (NN) and next-nearest-neighbor (NNN) couplings for the site pairs (a, b) in the 1st and 3rd terms. The parameters J ab (= J ba ), V 0 , and V ab (= V ba ) are the coefficients of the hopping, the on-site interaction, and the intersite interaction, respectively, and calculable by using the Wannier functions in a certain band. The intersite terms V ab may arise when the atoms have a long-range dipole-dipole interaction (DDI) [25], or when the atoms are populated in the excited bands of the optical lattice [26].
To map the BH model onto the hamiltonian of LGT, we consider the diagonal lattice whose sites r = (r 1 , r 2 ) are positioned on the centers of the colored squares in Fig. 1. Then, the original sites can be viewed as links of the diagonal lattice. The links are labeled as (r, i) with the direction index i = 1, 2. To derive the hamiltonian of the target GH model, we consider the case such that J ab and V ab take values according to the following three groups (i)-(iii) for pairs (a, b) of sites as shown in Table I [17,19]. We note that Table I breaks the translational symmetry of atomic interactions, e.g., V 24 = 0 while V 15 = 0. Next, we assume that the equilibrium atomic density is uniform and sufficiently large ρ 0 ≡ ρ r,a 1. Then, we expand the density operator asρ r,i = ρ 0 +η r,i , Relation between two 2D lattices. The dashed red lines indicate the 2D optical lattice with the square geometry, and cold atoms reside on its sites denoted by black crosses. Its unit cell consists of a pair of white and blue squares (yellow region). The filled black lines indicate the 2D gauge lattice on which the U(1) lattice GH model is defined. Then the cold atoms are viewed to sit on each link of the gauge lattice to play the role of gauge field. The relevant sites of the original lattice (links among the site r of the gauge lattice) in Table I are numbered as a = 1 ∼ 6, as a = 1 = (r, 1), 2 = (r − 2, 2), 3 = (r − 1, 1), 4 = (r, 2), 5 = (r + 2, 1), 6 = (r + 1, 2), where (r, i) represents the link of the gauge lattice emanating from the site r into the positive i (= 1, 2)-th direction. The gauge lattice plays an important role also in the study of the BH model with V ab = 0 in a different context [24]. and keep terms up to O(η 2 ) to obtain where V 0 ≡ V 0 − 2γ −2 > 0, (r, δ) represents the NN links of (r, i), and1 = 2,2 = 1. The first order term O(η) is absent due to the stability condition for ρ 0 = [µ + 4J + 2(J + J )]/(V 0 + 8γ −2 ) with the chemical potential µ. In the atomic simulators of LGT, the phaseθ r,i plays a role of a gauge variable on the link (r, i) and its conjugate momentumη r,i is the electric field −Ê r,i [17][18][19]. By replacingη r,i → (−) rη r,i and θ r,i → (−) rθ r,i with (−) r ≡ (−) r1+r2 , the first term in the rhs of Eq. (2) describes the "Gauss's law" as Table I are necessary to generate the (∇ · E) 2 term without nonlocal interaction among E r,i . If these conditions are not fulfilled, a productÊ r,iÊr ,i over the different links appears additionally, and it gives rise to long-range interactions among the gauge field θ r,i in the target GH model. Although such a model still respects gauge symmetry, we reject it here in favor of short-range GH models.
In Ref. [19], it was shown that the partition function Z = Tr exp(−βĤ) of the atomic model of Eq. (2) is equivalent to that of the GH model. The GH model is the U(1) lattice gauge model on the (2+1)D lattice, and its partition function is given by Here, x = (x 0 , r 1 , r 2 ) is the site index of the (2+1)D lattice with the discrete imaginary time τ = x 0 × ∆τ [x 0 = 0, · · · , N 0 , N 0 ∆τ = β ≡ (k B T ) −1 ] and the 2D spatial coordinate r 1 , r 2 . µ and ν (= 0, 1, 2) are direction indices. The U(1) gauge variables U x,µ ≡ exp(iθ x,µ ) are defined on the link (x, µ). θ x,µ (µ = 1, 2) corresponds to the eigenvalue of the phase of atomic operatorψ ra through θ x,µ =θ r,i . The complex field φ x defined on site x is a bosonic matter field, referred to as "Higgs field" in the London limit, taking the form φ x = exp(iϕ x ) with frozen radial fluctuations. The integration [dU ][dφ] is over the angles θ x,µ , ϕ x ∈ [0, 2π). The coefficients c 1 ∼ c 5 are real dimensionless parameters for interactions among gauge fields. Each term of the action, hence the action itself, and the integration measure are invariant under the local U(1) gauge transformation, According to Ref. [19], the atomic simulator of the GH model in a 2D system corresponds to the following case of parameters for c 1µ , c 2µν : In terms of the atomic system, the c 1 -and c 2 -terms describe the sum of the self coupling and the neighboring correlations of densities of atoms, and the c 3 -and c 4,5terms describe the NN and the NNN hopping terms respectively. The relations among the parameters of Eq. (2) and Eq. (3) are In experiments, we expect low T ( 10 nK set by the parameters ofĤ), and the quantum phase transitions may be explored in a multi-dimensional space parametrized by the dimensionless and ∆τ -independent combinations such as c 1 /c 2 = γ 2 V 0 , c 3 c 2 = 2Jρ 0 /V 0 , etc.

III. IMPLEMENTATION WITH COLD ATOMS
In this section we present two methods to realize V ab as shown in Table I. A major way to prepare intersite interactions in BH systems is to use DDI between atoms or molecules [25,27,28]. In usual experiments, dipoles of an atomic cloud are uniformly polarized along a certain direction, and one may easily check that uniformly oriented dipoles generate V ab different from the configuration of V ab in Table I. This is partly because we consider a square lattice, and the similar requirement for V ab is satisfied on the triangular or Kagome lattice [18]. Although an individual control of the polarization of a dipole at each site may achieve V ab in Table I, its actual fulfillment is difficult (some discussions can be seen in the system of polar molecules [29]), and importantly the hopping process between sites with different dipole orientations are prohibited or reduced due to the conservation of the atomic spin. Recently, there is an interesting proposal to realize V ab in Table I by using the Rydberg p-states [30], but it also suffers from the same problem about the hopping.
In Sec. III A, we discuss the possibility to realize the values of V ab in Table I by using the excited bands of an optical lattice, which is an alternative route to get intersite interactions [26]. In Sec. III B, we discuss a system of multi-layer 2D optical lattices [31] to realize tunable DDI between of atoms. These proposals are within reach in current experimental techniques.

A. Method A: Using excited bands of an optical lattice
The Wannier functions in excited bands have extended anisotropic orbitals compared with the lowest s-orbital Table I  band. Thus, we expect the significant intersite densitydensity interaction without introducing DDI between atoms [26]. To implement this scheme, we assume the following optical lattice potential: which can be created in a current experimental setup. For V B /V A > 0.5, the potential forms a checkerboard lattice (line graph of a square lattice [24]) and its minima are characterized by anisotropic harmonic form as shown in Fig. 2(a). This anisotropy is necessary to prevent the intraband mixing dynamics. Excitation to higher orbitals can be achieved by stimulated Raman transition [32] or nonadiabatic control of the optical lattice [33,34].
The intersite density-density interaction is proportional to the overlap integral O ab = dr|w a | 2 |w b | 2 , where w a is the Wannier function at the link (r, a) and we assume a negligibly small DDI. For the horizontal links, by approximating a minimum of the optical lattice as a quadratic form mω 2 ho (α 2 x 2 +y 2 )/2, w a can be represented by the harmonic oscillator basis Φ n (r), where ω ho = 2 E R (V A + 2V B ) with the recoil energy E R = 2 k 2 /2m of the optical lattice and α = ( The band index takes n = (0, 0), (1, 0), and (2, 0) for the s-, p-, and d-orbitals. For the vertical links, the role of (x, y) is just exchanged by (y, x).
The conditions in Table I represent the parameter domain satisfying this condition for α andR 0 (defined as R 0 /a ho using a harmonic oscillator length a ho = /mω ho ). These parameter domains are experimentally realistic becausẽ R 0 ∼ 10 can be obtained for V A(B) ∼ 10 − 100E R and α is also a precisely tunable parameter. Even for the s-orbitals the condition is satisfied in a small but finite region for a highly asymmetric potential with α 1. Using the p-or d-orbitals allows us to get the condition of Table I more easily in the experimentally feasible condition. When the excited orbitals are used, we have significant hopping amplitudes J ab not only for the NN (J) but also the 1st half of NNN (J ); the 2nd half of NNN (J ) would be small because of higher potential height between the link of group (iii) as seen in Fig. 2(a).
Finally, we admit that, for actual parameter estimation, one should also try other more realistic Wannier functions such as ∼ x c exp(−h|x|) [35], although the qualitative feature captured here needs no modifications.

B. Method B: Using dipolar atoms in a multilayer optical lattice
The idea for the second method is to introduce new subsidiary 2D lattices and treat the DDI between atoms in the original 2D lattice and atoms in the subsidiary lattices by the second-order perturbation theory to obtain V ab effectively. For illustrative purpose, we explicitly describe the idea by using a triple-layer system consisting of three 2D square optical lattices (layer L A , L B , L C ) as seen in Fig. 3(a). Here, we neglect the contribution of short-range interaction for the intersite interaction. The scheme may be reduced to a double-layer system by approaching the distance between two layers, e.g. L A and L B , to zero, which is discussed for realistic parameter estimation at the end of Appendix C.
The boson system on the layer L A (we call them Abosons) is a playground of the (2+1)D U(1) GH model, which is sandwiched by B-bosons on L B and C-bosons on L C . The B-and C-bosons are trapped in deep optical lattices with negligible hopping. Each layer has different basis vectors of the lattice structure as shown in Fig. 3(b) and (c). Each species of bosons is assumed to have a dipole, perpendicular to the plane of the layer. By treating the DDI between A-boson and B-boson as a perturbation, the second-order perturbation theory generates an effective intersite interaction between the A-bosons. So is the DDI between A-and C-bosons, which generates another intersite interactions between the A-bosons. respectively. For the layers LB and LA, we take the NN DDI between atoms on the sites k in LB and k ± δx(δy) in LA into account. For LC and LA, we take the NN DDI between atoms on the sites in LC and +δ+(δ−) in LA into account.
These two kinds of interactions may be tuned to realize V ab as given in Table I. We omit the DDI between the Band C-bosons because of the large separation.
Let us focus on Fig. 3(d). When one projects the sites of L B onto L A , their image locates on the center of each plaquette of the L A lattice. Similarly, the image of sites of L C locates on the middle of NN pairs of the L A sites. In L A , the A-bosons at different sites have the repulsive DDI. Furthermore, the A-and B(C)-bosons are coupled through the NN attractive DDI given by where ρ A,k and n B(C),k are boson densities at the site k and V AB(C) < 0 is the DDI, which is tunable by controlling the interlayer separation. Our strategy is to trace out B-and C-bosons to get the effective attractive intersite interactions between the A-bosons themselves. According to the usual second-order perturbation theory with H AB(C) as perturbation, the effective attractive interaction between A-bosons may be estimated as ∼ −V 2 AB k,δ,δ ρ A,k+δ ρ A,k+δ and −V 2 AC l,δ ρ A,l+δ ρ A,l−δ .
They are due to density fluctuations of B-and C-bosons, respectively. The former term contributes a constant to V ab for (a,b) of the groups (i,iii) of Table I, while the latter contributes a constant only for the group (i). Then one may fulfill the condition of V ab in Table I. The detailed calculation of the effective interaction and the experimental feasibility are described in Appendix C.

IV. REAL-TIME DYNAMICS OF SIMULATORS: STABILITY OF AN ELECTRIC FLUX
In actual experiments observing the nonequilibrium time-evolution of a quantum simulator, the results globally reflect the phase structure of the target model. The (2+1)D GH model supports the confinement phase and the Higgs phase (see Appendix A). The confinement phase is characterized by the strong phase fluctuation; when static two point charges, such as density defects created by the focused potentials, are put on, they are connected by an almost straight electric flux (linearly-rising confinement potential). In contrast, the Higgs phase possesses the phase coherence over the system and the system can be regarded as a superfluid phase; the density wave can propagate around the charges [19].
To get some insight on the time-evolution as a guide to experiments, we study the dynamical features of the simulators through numerical simulations under the meanfield approximation of the two quantum hamiltonians: the base BH model Eq. (1) and the target GH model Eq. (2). The time-dependent equations used in Sec. IV can be derived from the real-time path-integral formulation under the saddle-point approximation (we put = 1). The operators of the original hamiltonian are replaced by the c-number fields. We confine ourselves to the models with only NN hopping J = 0 and J = J = 0 for simplicity. We note in advance that the mean-field equations necessarily underestimate quantum fluctuations, and their results should be taken as a guide to practical and future experiments which are expected to reveal the real dynamics of quantum systems.
The equation of motion for ψ in the BH model of Eq. (1) can be derived from the Lagrangian L = − r i=1,2 iψ * r,i (dψ r,i /dt) − H. It is the discretized version of the GP equation called the discrete nonlinear Schorödinger equation [36] and given by i ∂ψ r,i ∂t = −J(ψ r,ī + ψ r−ī,ī + ψ r+i,ī + ψ r+i−ī,ī ) where i = 1, 2 and1 ≡ 2,2 ≡ 1. The uniform stationary solution can be obtained by substituting ψ r,i = ψ 0 e −iµt as |ψ 0 | 2 = (µ + 4J)/(V 0 + 8γ −2 ), where µ is the chemical potential. Since an important quantity to observe the (2) In terms of the optical lattice, the summation over j of Eq. (9) implies to take over the four atomic sites which are NN to the atomic site (r, i) (i = 1, 2). In terms of the gauge lattice, given an atomic link (r, i), (r, j) takes (r,ī), (r−ī,ī), (r+i,ī), (r+i−ī,ī). Equations (9) and (10) can be also derived by linearizing Eq. (8) with respect to the density ρ r,i (t) = ρ 0 + η r,i (t). The constraint of the Gauss's law requires the replacement η r,i → (−1) r η r,i and θ r,i → (−1) r θ r,i . We make a dimensionless form of Eq. (8) and (10)  As an explicit example to apply the dynamical equations, we consider the dynamical stability of a single straight flux connecting two external charges, which is prepared as an initial condition. The details of the numerical simulations are described in Methods. In the confinement phase, a set flux string should be stable. The free parameters of this system are (γ 2 , V 0 , J), related to (c 1 , c 2 , c 3 ). By using the ∆τ -independent parameters, we expect the confinement (Higgs) phase for small (large) values of c 1 /c 2 = γ 2 V 0 and c 3 c 2 = 2Jρ 0 /V 0 (see Appendix A). Figure 4(b) represents the time evolution of the density distribution η 2 calculated by the above two models. For a certain value of J, both models show similar behaviors for small values of γ 2 , where the placed density flux is stable and does not spread out. This captures the characteristics of the confinement phase. With increasing γ 2 , i.e., the Higgs coupling, the structure of the density flux is gradually lost by emitting the density waves from the charge. This phonon emission is a characteristic of the superfluid phase, i.e, Higgs phase. The density waves are generated in a different way: successively in the BH model and intermittently in the GH model, propagating concentrically around the point charges with the sound velocity ∼ Jρ 0 V 0 for γ 2 1. To judge whether the system is in confinement-or Higgs-regime by dynamical simulations, we calculate the remnants of the flux σ(t) defined by where the sum is taken over the sites on which the density flux line is set initially. The flux is stable when σ is kept small during the time evolution. Figure 4 (a) shows the dynamical phase diagram obtained by the behavior of σ shown in Fig. 4(c). The rapid oscillation of σ reflects in the periodic vanish-revival cycle of the density flux. In the BH model, we calculate the time-average σ t and determine the phase boundary by finding the point at which σ t almost vanishes. In the GH model, the boundary is determined by the appearance of rapid growth of σ(t) due to the intermittent density-wave emission. The dynamical phase boundaries of both models are qualitatively consistent with the result of the Monte Carlo simulations (see Fig. 5 and Appendix A). The difference of the two models can be observed in the amplitude fluctuation of the simulating gauge field. Because the GH model is obtained by expanding ρ = ρ 0 + η around the constant density ρ 0 1, the BH model can approximately reproduce the GH model when the Thomas-Fermi limit is satisfied; note that the boundary of the BH model in Fig. 4(a) is obtained for the particular value ρ 0 = 1. In addition, near the situation µ = 0 (dotted curve in Fig. 4(a)), the density fluctuation is accidentally frozen and the dynamics of the BH model is similar to the GH model. This is a reason of the decrease of σ t around γ 2 V 0 = 2.5.  3)) is always in the confinement phase [37], in which the phase θ x,µ is unstable by strong fluctuation. With increasing matter couplings (c 1,3,4,5 ) from zero, the system enters into the "Higgs" phase (see Fig. 5), where both θ x,µ and ϕ x are stable.
To identify the location of the transitions, we measure the internal energy U = A and the specific heat C = A 2 − A 2 by using the standard Metropolis algorithm in Monte Carlo simulation with the periodic boundary condition for the cubic lattice of size V = L 3 with L up to 40. The typical number of sweeps is 100000 + 10000 × 10, where the first number is for thermalization and the second one is for measurement. The errors of U and C are estimated by the standard deviation over 10 samples. Acceptance ratios in updating variables are controlled to be 0.6 ∼ 0.8.
Explicitly, we confine ourselves to the case c 4 = c 5 = 0 and obtain the phase diagram in the c 1 −c 3 plane for several values of c 2 . The result is presented in Fig. 5. There are two phases: the Higgs phase in the large c 3 region (upper region) and the confinement phase in the small c 3 region (lower region). The confinement-Higgs transition here should correspond to various phase transitions such as the superconducting transition, the mass generation in the standard model, and the believed one to take place in the early universe [22,23]. In contrast to the phase diagram of the (3+1)D model for c 4 = c 5 = 0, c 1 = c 3 , and c 2 ≥ 0 [19], the Coulomb phase is missing due to the low dimensionality.
To understand Fig. 5, let us consider some limiting cases. First, after choosing the unitary gauge φ x = 1, let us consider the limit c 1 → ∞. Then the c 1 term makes θ x,0 = 0 [mod(2π)], and the action becomes up to constant. This is viewed as a 3D XY spin model with asymmetric couplings, where θ x,i (i = 1, 2) on the link (x, x + i) is the XY spin angleθx. In fact, the c 3 term is their NN coupling in the 12 plane and the c 2 term is their NN coupling along the µ = 0 axis. The region of sufficiently large c 2 and c 3 is the ordered phase of this XY spins, and corresponds to the Higgs phase with small gauge-field (θ x,µ ) fluctuations. As a check of Fig. 5, let us consider the case c 2 = c 3 of Eq. (A1), which reduces to the symmetric 3D XY spin model of A 3DXY = c XY x,µ cos(θx +µ −θx). It is known to have a genuine second-order phase transition at c XY 0.45. Therefore the transition line in Fig. 5 (b) should approach to c 2 · c 3 → 0.45 2 0.20 as c 1 /c 2 → ∞ as it shows.
Next, let us consider the case c 2 = 0. Then, each variable θ x,0 appears only through the c 1 term without couplings to other variables (we take the unitary gauge as before). Then the dynamics is controlled by the c 3 term. Again, this term is viewed as the energy of the XY spins θ xi (i = 1, 2). However they have no coupling along the µ = 0 direction, and therefore the system is a collection of decoupled 2D XY spin models. 2D XY spin model is known to exhibit Kosterlitz-Thouless (KT) transition which is infinitely continuous. Thus, although it is not drawn in Fig. 5, there should be added a horizontal line (independent of c 1 ) for c 2 = 0 consisting of a collections of KT transitions at around c 3 ∼ 0.96. We understand that the crossover points appearing in the smaller c 1 part in each curve for three c 2 drawn in Fig. 5 are the rem- which separate the Higgs phase (above) and the confinement phase (below). The transition points are located at the peak of C as a function of c3 for fixed c1. They are (i) secondorder (no marks) where the peak of C develops as the size L increases and U exhibits no hysteresis or (ii) cross-over or Kosterlitz-Thouless transition (both marked with co) where the peak does not develop. (b) The transition points on (a) are arranged in the (c1/c2) − (c3 · c2) plane (γ 2 V 0 − 2Jρ0/V 0 plane) . They almost sit on a single curve (gray line) (See the comment right after Eq. (5)).
nants of these KT transitions. They have a chance to be a genuine KT transition, although we called them crossover here. Another support of this interpretation is to consider the case c 1 = 0. Then there is no source term for θ x,0 and θ x,0 should determine their dynamics only through the c 2 term. Thus, even θ x,i could be set constant with no fluctuations, θ x,0 has the NN coupling in each 12 plane. However, two dimensions is not enough to stabilize θ x,0 . In turn, the c 2 term is not enough to sustain the coupling between θ x,i along the µ = 0 direction. The dynamics of θ x,i is essentially from the c 3 term, which is the 2D XY model as explained. Therefore, the transition, if any, for c 1 = 0 may be a KT transition. No genuine second-order one is possible.
The last case is c 2 = ∞. Then θ x,µ is frozen to be a pure gauge configuration, θ x,µ = λ x+µ − λ x . By plugging this into the c 1 and c 3 term, we obtain which belongs again to the class of 3D XY spin models, where λ x is the XY spin angles on the site x. So we should have a second-order transition at c 2 = ∞ as long as both c 1 and c 3 are nonvanishing. This is consistent with Fig. 5. Let us finally comment on the transition line of Fig.  5(b) and the boundaries of Fig. 4(a) calculated by dynamical simulation in Sec. IV. Their behaviors in the (c 1 /c 2 )-(c 3 · c 2 ) plane are qualitatively consistent but different in quantitative comparison. We understand that there are no inconsistency in these results because the two methods, MC and GP, are different in nature: MC is static and GP is dynamical, they treat fluctuations in contrasting manners, and especially, the dynamical simulations necessarily exhibit various properties of the system according to their setup and probes, etc. This certainly motivates exact quantum and dynamical simulation of the BH model in experiments.

Appendix B: Calculation of overlap integrals
In this section, we describe the calculation of the overlap integrals discussed in Sec. III A. For the horizontal links of the potential V OL of Eq. (6), the minimum is approximated by the harmonic oscillator V hx = mω 2 ho (α 2 x 2 + y 2 )/2. The basis function of V hx is given by where H n is the Hermite polynomial, A n is the normalization factor, and the length scale is normalized by the harmonic oscillator length a ho = /mω ho ; the dimensionless coordinates are denoted by putting tildes. The s, p, and d orbitals for these links correspond to n = (n x , n y ) = (0, 0), (1, 0), and (2, 0). As the Wannier function w a (r) at the link (r, a), we use Φ n withx,ỹ measured from the center of the link. For the vertical links, the minimum is written as V hy = mω 2 ho (x 2 + α 2 y 2 )/2 and the basis function is Φ n (x, √ αỹ) = AH nx (x)H ny ( √ αỹ)e −(x 2 +αỹ 2 )/2 . The s, p, and d orbitals for these links correspond to n = (n x , n y ) = (0, 0), (0, 1), and (0, 2). w a (r) are given as follows; whereR 0 represents the lattice spacing in unit of a ho and we shift the origin of the coordinate to (R 0 /2,R 0 /2) of Fig. 2(a). The intersite interaction strength V ab is proportional to the overlap integrals O ab = dr|w a | 2 |w b | 2 . It is sufficient to calculate only the integrals for the link pairs O 15 with decreasing α further as shown in Fig. 2(b). The advantage of the use of p-or d-orbitals is that the region of parameters for O 12 ≈ O 13 O 15 becomes wider because of the extended amplitude profile of the wave functions. These behaviors are clearly seen in Figs. 2(c) and (d).
Appendix C: Effective intersite interaction in the triple-layer system of Sec. III B In this section, we apply the second-order perturbation theory to the triple-layer system in Sec. III B to derive the effective intersite interaction of A-bosons, and estimate the possible values of involved parameters to realize V ab of Table I. After that, we briefly explain a double-layer system in which magnitude of the intersite interaction of A-bosons is controlled in a similar way.
We first confine ourselves to the subsystem of the Aand B-bosons (two layers L A and L B ) which has the NN DDI, H AB of Eq. (7). It implies that the B-boson on the site k interacts with the four NN A-bosons on the sites k ± δ x(y) as seen in Fig. 3(d). V AB in H AB is expressed as where r(r ) is the position of A(B)-boson, w A(B) (r) is their Wannier function, and C = µ 0μAμB ; µ 0 is the magnetic permeability of the vacuum andμ A(B) is the magnetic moment of A(B)-atoms.
We assume that the B-bosons of L B have a chemical potential µ B (> 0), a negligibly small NN hopping amplitude due to a deep trapping potential, an on-site repulsion U B (> 0), and the NN DDI with A-bosons V AB . One may forget the DDI between B-bosons, because it is a constant due to negligible NN hopping. Then, the HamiltonianĤ B and the partition function Z B = Tr exp(−βĤ B ) for the subsystem of B-bosons are written by using the B-boson density operatorn k at the site k aŝ Then we have The first-order terms ∝ W k are renormalized to the chemical potential of A-bosons and the second-order terms define the effective density-density interaction HamiltonianĤ ABA of A-bosons induced by B-boson density fluctuation (∆n) 2 , The DDI between A-atoms and C-atoms can be analyzed in the same way, and we obtain another effective density-density interaction for the A-bosons,Ĥ ACA = −(β/2)(∆n ) 2 V 2 AC ,δ,δ ρ +δρ +δ , where (∆n ) 2 is obtained by replacing µ B , U B by µ C , U C in (∆n) 2 .
The sumĤ ABA +Ĥ ACA contributes to the coefficients V ab of the intersite density-density interactions for Abosons as follows; For NN links (group (i) in Table I), Each x-y plane is separated by the distance z and has the 2D lattice structure of the superposition of the three layers LA,B,C shown in Fig. 3. Then we supply the A,B,C atoms so that every x − y layer is viewed as the triple-layer system of Fig. 3 with AB = AC = 0. (ii) We blow almost all the atoms off in such a way that only the A and B atoms in one x − y layer and the C atoms in the next layer left. The resulting double-layer system is viewed as the triple-layer system of Fig.  3 with AB = 0 and AC = z .
For NNN links, where V and V are the direct DDI for NN and NNN link pairs, respectively. The condition for V ab in Table I can be established by adjusting two inter-layer distances AB and AC and density fluctuations (∆n) 2 (µ B , V B , β) and (∆n ) 2 (µ C , V C , β) as Let us present some brief account for an example and estimation of the experimental parameters that satisfy the tuning relations Eq.(C7). We shall report detailed discussion on this example and related topics in a future publication.
For bosons loaded in each layer we consider 52 Cr atoms [38] as A bosons, 87 Rb atoms as B bosons, and 168 Er atoms [39] as C bosons. They have the permanent magnetic moments 6µ BM , µ BM , and 7µ BM (µ BM is a Bohr magneton), respectively. Then we are interested in the effective double-layer system, which is obtained from the triple-layer system explained above by choosing AB = 0. The reason for using the double-layer system is to make the intersite interaction as large as possible because the magnetic moment of 87 Rb atom is small.
The method to make such a double-layer system is sketched in Fig. 7. First, one prepares the 3D layer system as shown in Fig. 7(i) by emitting three standing waves with the wavelengths satisfying 2λ 1 = √ 2λ 2 = λ 3 (e.g., λ 1 = 410nm, λ 2 = 580nm and λ 3 = 820nm) in eight appropriate directions in the x-y plane, each being separated by 45 degrees. In addition, we emit another standing-wave laser in the z-direction with the wavelength λ z to establish the 3D structure. Because 52 Cr, 87 Rb and 168 Er exhibit the specific strong absorptions of photon with wave length 425nm, 780nm and 401nm, respectively, above standing waves load these atoms to the sites of corresponding layer L A,B,C [40]. This completes the step (i) in Fig. 7.
In the second step (ii) on Fig. 7, one needs to remove almost all the atoms except for those in two adjacent x-y layers. This can be experimentally realized by using the technique of a position-dependent microwave transfer in a magnetic field gradient perpendicular to the layers [41] successively. This achieve to make an effective doublelayer system with AB = 0, AC = z .
Finally, let us estimate the parameters to satisfy the tuning relation Eq. (C7). By making a straightforward calculation using DDI, we find that the following is a typical example of parameters : The average densities per site are ∼ 560 for B bosons and ∼ 8 for C bosons. The ratio |V AB(C) /µ B(C) | is ∼ 0.192(∼ 0.585), which seems to validate the perturbation theory.