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

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–4]. Since its introduction, LGT has been an indispensable tool for studying the non-perturbative aspect of quantum models in 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–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–14]. Although 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 in the 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 Gauss's law and the short-range gauge interactions simultaneously, a complex design of the system and a fine-tuning of interaction parameters are generally needed in the experimental setups [17, 18].

Recently, we proposed a perspective to overcome the difficulty with respect to 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 (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 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 the 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 such as the dynamics of Higgs phase transition and the related problems of topological defects, uniformity, etc [22, 23]. Atomic quantum simulation of this model is certainly welcome because it simulates the real-time evolution of the previously mentioned exciting phenomena.

In what follows, we consider the second approach to discuss 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. In the experiments, each phase could be generally studied through the non-equilibrium dynamics of the system, which are detected by, e.g., the density distribution of the time-of-flight imaging after the system is perturbed. 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 dynamic stability of a single electric flux connecting two charges with opposite signs, corresponding to a density hump and dip for the atomic simulators. We stress that this dynamic simulation in an interacting atomic system gives a new theoretical tool for the analysis of lattice gauge models far beyond the present standards of the theoretical study on the LGT using the 'classical' Monte Carlo simulations, the strong coupling expansion, etc. The obtained phase boundary is discussed and compared with that of the Monte Carlo simulations. Next, we propose two realistic experimental setups for the quantum simulators. To respect the constraint of Gauss's law and avoid nonlocal gauge interactions, it is necessary to tune suitably the intersite density-density interaction of the hamiltonian. We give two ideas: (i) using Wannier states in the excited bands and (ii) using dipolar atoms in a multilayer optical lattice, both of which are reachable under current experimental techniques.

The paper is organized as follows. In section 2, we introduce our target hamiltonian of the U(1) lattice GH model starting from the extended Bose–Hubbard (BH) model with the intersite density-density interaction. The exact correspondence of the atomic system to the LGT has been discussed in [19]. In section 3, we present the results of the dynamic simulations by means of the GP equation, which is obtained under the saddle-point approximation of the real-time path integral of the two quantum hamiltonians, i.e., the original BH model and the target GH model. The obtained dynamic phase diagram is compared with the result of the Monte Carlo simulations. We give two proposals for experiments to construct realistic atomic simulators for the lattice GH model in section 4. Method A in section 4.1 relies on extended orbits of the Wannier functions in excited bands of an optical lattice. Method B utilizes the long-range interaction between atoms in different layers of 2D optical lattices. Both methods may be used to tune the intersite density-density interactions in the 2D BH system for our purpose. Section 5 is devoted to our conclusion and an outlook on future directions.

Corresponding to the simplest realistic experimental situation of the quantum simulator, we focus on the boson system defined on a 2D square lattice. We start from a generalized BH hamiltonian [1]

$$\begin{eqnarray}&&\hat{H}=-\mathop{\sum }\limits_{k,a\ne b}{{J}_{ab}}\hat{\psi }_{a}^{\dagger }{{\hat{\psi }}_{b}}+\frac{{{V}_{0}}}{4}\mathop{\sum }\limits_{k,a}{{\hat{\rho }}_{a}}({{\hat{\rho }}_{a}}-1)+\mathop{\sum }\limits_{k,a\ne b}\frac{{{V}_{ab}}}{2}{{\hat{\rho }}_{a}}{{\hat{\rho }}_{b}},\end{eqnarray} \tag{ 1 }$$

which describes the bosons in a single band of a 2D optical lattice. The bosonic atomic fields ${{\hat{\psi }}_{a}}={\rm exp} (i{{\hat{\theta }}_{a}})\sqrt{{{{\hat{\rho }}}_{a}}}$ are put on the site a of the square optical lattice. The summation is taken over the unit cell k (yellow region in figure 1) and, in each unit cell, over the the site $a(b)\in 1$–6. We confine ourselves to the nearest-neighbor (NN) and next-nearest-neighbor (NNN) couplings for the site pairs $(a,b)$ in the first and third terms. The parameters ${{J}_{ab}}(={{J}_{ba}})$, V₀, and ${{V}_{ab}}(={{V}_{ba}})$ are the coefficients of the hopping, the on-site interaction, and the intersite interaction, respectively, and are 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].

Zoom In Zoom Out Reset image size

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 figure 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 1 [17, 19]. We note that table 1 breaks the translational symmetry of atomic interactions, e.g., ${{V}_{24}}\ne 0$ while ${{V}_{15}}=0$. Next, we assume that the equilibrium atomic density is uniform and sufficiently large ${{\rho }_{0}}\equiv \langle {{\hat{\rho }}_{r,a}}\rangle \gg 1$. Then, we expand the density operator as ${{\hat{\rho }}_{r,i}}={{\rho }_{0}}+{{\hat{\eta }}_{r,i}}$, and keep terms up to $O({{\hat{\eta }}^{2}})$ to obtain

$$\begin{eqnarray} {\hat{H}} & \simeq & \frac{1}{2{{\gamma }^{2}}}\mathop{\sum }\limits_{r}{{\left[ \mathop{\sum }\limits_{i}({{{\hat{\eta }}}_{r,i}}+{{{\hat{\eta }}}_{r-i,i}}) \right]}^{2}}+\frac{V_{0}^{\prime }}{2}\mathop{\sum }\limits_{r,i}\hat{\eta }_{r,i}^{2} \\ {} & {} & -{{\rho }_{0}}J\mathop{\sum }\limits_{r,i,\delta }{\rm cos} ({{{\hat{\theta }}}_{r,i}}-{{{\hat{\theta }}}_{r,\delta }})-2{{\rho }_{0}}J^{\prime} \mathop{\sum }\limits_{r,i}{\rm cos} ({{{\hat{\theta }}}_{r,i}}-{{{\hat{\theta }}}_{r-i,i}})-2{{\rho }_{0}}J^{\prime\prime} \mathop{\sum }\limits_{r,i}{\rm cos} ({{{\hat{\theta }}}_{r,i}}-{{{\hat{\theta }}}_{r-\bar{i},i}}), \end{eqnarray} \tag{ 2 }$$

where $V_{0}^{\prime }\equiv {{V}_{0}}-2{{\gamma }^{-2}}\gt 0$, $(r,\delta )$ represents the NN links of $(r,i)$, and $\bar{1}\equiv 2$, $\bar{2}\equiv 1$. The first-order term $O(\hat{\eta })$ is absent due to the stability condition for ${{\rho }_{0}}=[\mu +4J+2(J^{\prime} +J^{\prime\prime} )]/(V_{0}^{\prime }+8{{\gamma }^{-2}})$ with the chemical potential μ. In the atomic simulators of LGT, the phase ${{\hat{\theta }}_{r,i}}$ plays a role of a gauge variable on the link $(r,i)$ and its conjugate momentum ${{\hat{\eta }}_{r,i}}$ is the electric field $-{{\hat{E}}_{r,i}}$ [17–19]. By replacing ${{\hat{\eta }}_{r,i}}\to {{(-)}^{r}}{{\hat{\eta }}_{r,i}}$ and ${{\hat{\theta }}_{r,i}}\to {{(-)}^{r}}{{\hat{\theta }}_{r,i}}$ with ${{(-)}^{r}}\equiv {{(-)}^{{{r}_{1}}+{{r}_{2}}}}$, the first term in the rhs of equation (2) describes 'Gauss's law' as ${{(2{{\gamma }^{2}})}^{-1}}{{\sum }_{r}}{{[{{\sum }_{i}}({{\hat{\eta }}_{r,i}}-{{\hat{\eta }}_{r-i,i}})]}^{2}}\simeq {{(2{{\gamma }^{2}})}^{-1}}{{\sum }_{r}}{{(\nabla \cdot {\bf E})}^{2}}$. The two conditions ${{V}_{({\rm i})}}={{V}_{({\rm ii})}}(={{\gamma }^{-2}})$ and ${{V}_{({\rm iii})}}=0$ in table 1 are necessary to generate the ${{(\nabla \cdot {\bf E})}^{2}}$ term without non-local interaction among ${{E}_{r,i}}$. If these conditions are not fulfilled, a product ${{\hat{E}}_{r,i}}{{\hat{E}}_{r^{\prime} ,i^{\prime} }}$ over the different links appears additionally, and it gives rise to long-range interactions among the gauge field ${{\theta }_{r,i}}$ in the target GH model. Although such a model still respects gauge symmetry, we reject it here because all LGTs relevant to HEP are generally models with local interaction.

Table 1.  Atomic parameters J_ab and V_ab in equation (1). Those not shown below are set to zero to avoid double counting ((1, 6), etc.) or due to longer-ranges ((3, 5), etc).

Group	Range	(a, b)	Jab	Vab
(i)	NN	(1, 2), (2, 3), (3, 4), (1, 4)	J	${{\gamma }^{-2}}$
(ii)	1st half of NNN	(1, 3), (2, 4)	J'	${{\gamma }^{-2}}$
(iii)	2nd half of NNN	(1, 5), (4, 6)	J''	0

In [19], it was shown that the partition function $Z={\rm Trexp}(-\beta \hat{H})$ of the atomic model of equation (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

$$\begin{eqnarray} {{Z}_{{\rm GH}}} & = & \int [dU][d\phi ]{\rm exp} ({{A}_{{\rm I}}}+{{A}_{{\rm P}}}+{{A}_{{\rm L}}}+{{A}_{2{\rm I}}}+{{A}_{{{{\rm I}}^{2}}}}), \\ {{A}_{{\rm I}}} & = & \frac{1}{2}\mathop{\sum }\limits_{x,\mu }{{c}_{1\mu }}({{{\bar{\phi }}}_{x+\mu }}{{U}_{x,\mu }}{{\phi }_{x}}+{\rm c}{\rm .c}{\rm .}), \\ {{A}_{{\rm P}}} & = & \frac{1}{2}\mathop{\sum }\limits_{x,\mu \lt \nu }{{c}_{2\mu \nu }}({{{\bar{U}}}_{x,\nu }}{{{\bar{U}}}_{x+\nu ,\mu }}{{U}_{x+\mu ,\nu }}{{U}_{x,\mu }}+{\rm c}{\rm .c}{\rm .}), \\ {{A}_{{\rm L}}} & = & \frac{{{c}_{3}}}{2}\mathop{\sum }\limits_{x}\left[ {{{\bar{\phi }}}_{x+1+2}}{{U}_{x+1,2}}{{U}_{x,1}}{{\phi }_{x}} \right.+{{{\bar{\phi }}}_{x+2}}{{{\bar{U}}}_{x+2,1}}{{U}_{x+1,2}}{{\phi }_{x+1}} \\ {} & {} & +{{{\bar{\phi }}}_{x}}{{{\bar{U}}}_{x,2}}{{{\bar{U}}}_{x+2,1}}{{\phi }_{x+1+2}}+\left. {{{\bar{\phi }}}_{x+1}}{{U}_{x,1}}{{{\bar{U}}}_{x,2}}{{\phi }_{x+2}}+{\rm c}{\rm .c}{\rm .} \right], \\ {{A}_{2{\rm I}}} & = & \frac{{{c}_{4}}}{2}\mathop{\sum }\limits_{x}\mathop{\sum }\limits_{i=1,2}\left[ {{{\bar{\phi }}}_{x+i}}{{U}_{x,i}}{{U}_{x-i,i}}{{\phi }_{x-i}}+{\rm c}{\rm .c}{\rm .} \right], \\ {{A}_{{{{\rm I}}^{2}}}} & = & \frac{{{c}_{5}}}{2}\mathop{\sum }\limits_{x}\left[ {{{\bar{\phi }}}_{x+2+1}}{{U}_{x+2,1}}{{\phi }_{x+2}}\cdot {{{\bar{\phi }}}_{x+1}}{{U}_{x,1}}{{\phi }_{x}} \right.+(1\leftrightarrow 2)+{\rm c}{\rm .c}{\rm .}]. \\ \end{eqnarray} \tag{ 3 }$$

Here, $x=({{x}_{0}},{{r}_{1}},{{r}_{2}})$ is the site index of the (2+1)D lattice with the discrete imaginary time $\tau ={{x}_{0}}\times \Delta \tau \ [{{x}_{0}}=0,\;\cdots ,{{N}_{0}},{{N}_{0}}\Delta \tau =\beta \equiv {{({{k}_{{\rm 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,\mu }}\equiv {\rm exp} (i{{\theta }_{x,\mu }})$ are defined on the link $(x,\mu )$. ${{\theta }_{x,\mu }}\ (\mu =1,2)$ corresponds to the eigenvalue of the phase of atomic operator ${{\hat{\psi }}_{r,a}}$ through ${{\theta }_{x,\mu }}={{\hat{\theta }}_{r,i}}$. The complex field ${{\phi }_{x}}$ defined on site x is a bosonic matter field, referred to as 'Higgs field' in the London limit, taking the form ${{\phi }_{x}}={\rm exp} (i{{\varphi }_{x}})$ with frozen radial fluctuations. The integration $\int [dU][d\phi ]$ is over the angles ${{\theta }_{x,\mu }},{{\varphi }_{x}}\in [0,2\pi )$. The coefficients ${{c}_{1}}\sim {{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, ${{\theta }_{x,\mu }}\to {{\theta }_{x,\mu }}+{{\lambda }_{x+\mu }}-{{\lambda }_{x}},\ {{\varphi }_{x}}\to {{\varphi }_{x}}+{{\lambda }_{x}}$.

According to [19], the atomic simulator of the GH model in a 2D system corresponds to the following case of parameters for ${{c}_{1\mu }},{{c}_{2\mu \nu }}$:

$$\begin{eqnarray} {{c}_{10}} & = & {{c}_{1}},\quad {{c}_{11}}={{c}_{12}}=0, \\ {{c}_{201}} & = & {{c}_{202}}={{c}_{2}},\quad {{c}_{212}}=0. \\ \end{eqnarray} \tag{ 4 }$$

In terms of the atomic system, the c₁ and c₂ terms describe the sum of the self-coupling and the neighboring correlations of densities of atoms, and the c₃ and ${{c}_{4,5}}$ terms describe the NN and the NNN hopping terms, respectively. The relations among the parameters of equations (2) and (3) are

$$\begin{eqnarray} {{c}_{1}} & = & \frac{{{\gamma }^{2}}}{\Delta \tau },\quad {{c}_{2}}=\frac{1}{\Delta \tau V_{0}^{\prime }},\quad {{c}_{3}}=2\;J{{\rho }_{0}}\Delta \tau ,{{c}_{4}} & = & 2J^{\prime} {{\rho }_{0}}\Delta \tau ,\quad {{c}_{5}}=2J^{\prime\prime} {{\rho }_{0}}\Delta \tau , \\ \end{eqnarray} \tag{ 5 }$$

In experiments, we expect low T(≲10 nK set by the parameters of $\hat{H}$), and the quantum phase transitions may be explored in a multi-dimensional space parameterized by the dimensionless and $\Delta \tau$-independent combinations such as ${{c}_{1}}/{{c}_{2}}={{\gamma }^{2}}V_{0}^{\prime },{{c}_{3}}{{c}_{2}}=2J{{\rho }_{0}}/V_{0}^{\prime },$ etc.

In actual experiments observing the non-equilibrium 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 into the time evolution of the system, we study the dynamic features of the simulators through numerical simulations under the mean-field approximation of the two quantum hamiltonians: the base BH model equation (1) and the target GH model equation (2). The time-dependent equations can be derived from the real-time path-integral formulation under the saddle-point approximation (we put $\hbar =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\ne 0$ and $J^{\prime} =J^{\prime\prime} =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 equation (1) can be derived from the Lagrangian $L=-{{\sum }_{r}}{{\sum }_{i=1,2}}i\psi _{r,i}^{*}(d{{\psi }_{r,i}}/dt)-H$. It is the discretized version of the GP equation called the discrete nonlinear Schorödinger equation [27] and given by

$$\begin{eqnarray} i\frac{\partial {{\psi }_{r,i}}}{\partial t} & = & -J({{\psi }_{r,\bar{i}}}+{{\psi }_{r-\bar{i},\bar{i}}}+{{\psi }_{r+i,\bar{i}}}+{{\psi }_{r+i-\bar{i},\bar{i}}}) \\ {} & {} & +\left[ \left( V_{0}^{\prime }+\frac{2}{{{\gamma }^{2}}} \right)|{{\psi }_{r,i}}{{|}^{2}}+\frac{1}{{{\gamma }^{2}}}\left( |{{\psi }_{r,\bar{i}}}{{|}^{2}}+|{{\psi }_{r-\bar{i},\bar{i}}}{{|}^{2}} \right. \right. \\ {} & {} & +|{{\psi }_{r+i,\bar{i}}}{{|}^{2}}+|{{\psi }_{r+i-\bar{i},\bar{i}}}{{|}^{2}}+\left. \left. |{{\psi }_{r-i,i}}{{|}^{2}}+|{{\psi }_{r+i,i}}{{|}^{2}} \right) \right]{{\psi }_{r,i}}, \\ \end{eqnarray} \tag{ 6 }$$

where $i=1,2$ and $\bar{1}\equiv 2,\bar{2}\equiv 1$. The uniform stationary solution can be obtained by substituting ${{\psi }_{r,i}}={{\psi }_{0}}{{e}^{-i\mu t}}$ as $|{{\psi }_{0}}{{|}^{2}}=(\mu +4J)/(V_{0}^{\prime }+8{{\gamma }^{-2}})$, where μ is the chemical potential. Since an important quantity to observe the dynamics of electric fluxes is the density fluctuation, we give the equilibrium density $|{{\psi }_{0}}{{|}^{2}}={{\rho }_{0}}$ by controlling the chemical potential as $\mu ={{\rho }_{0}}(V_{0}^{\prime }+8{{\gamma }^{-2}})-4J$ and see the evolution of the density fluctuation $\eta =\rho -{{\rho }_{0}}$.

The time-dependent equation of motion for η and θ in the GH model of equation (2) is derived in the similar way from $L=-{{\sum }_{r,i}}{{\eta }_{r,i}}(d{{\theta }_{r,i}}/dt)-H$ as

$$\begin{eqnarray}&&\frac{d{{\eta }_{r,i}}}{dt}=2J{{\rho }_{0}}\mathop{\sum }\limits_{j}{\rm sin} ({{\theta }_{r,i}}-{{\theta }_{r,j}}),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} \frac{d{{\theta }_{r,i}}}{dt} & = & -V_{0}^{\prime }{{\eta }_{r,i}}-\frac{1}{{{\gamma }^{2}}}({{\eta }_{r,i}}+{{\eta }_{r-i,i}}+{{\eta }_{r,\bar{i}}}+{{\eta }_{r-\bar{i},\bar{i}}})-\frac{1}{{{\gamma }^{2}}}({{\eta }_{r+i,i}}+{{\eta }_{r,i}}+{{\eta }_{r+i,\bar{i}}}+{{\eta }_{r+i-\bar{i},\bar{i}}}). \\ \end{eqnarray} \tag{ 8 }$$

In terms of the optical lattice, the summation over j of equation (7) implies the takeover of 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,\bar{i}),(r-\bar{i},\bar{i}),(r+i,\bar{i}),(r+i-\bar{i},\bar{i})$. Equations (7) and (8) can be also derived by linearizing equation (6) with respect to the density ${{\rho }_{r,i}}(t)={{\rho }_{0}}+{{\eta }_{r,i}}(t)$. The constraint of Gauss's law requires the replacement ${{\eta }_{r,i}}\to {{(-1)}^{r}}{{\eta }_{r,i}}$ and ${{\theta }_{r,i}}\to {{(-1)}^{r}}{{\theta }_{r,i}}$. We make a dimensionless form of equations (6)–(8) by using the energy scale $V_{0}^{\prime }$. In solving both sets of equations of motion, we use the $200\times 200$ discretized space and the time step $\Delta t={{10}^{-4}}$.

As an explicit example to apply the dynamic equations, we consider the dynamic stability of a single straight flux connecting two external charges, which is prepared as an initial condition. In the confinement phase, a set flux string should be stable. To see the stability of the flux configuration, we put the density modulation ${{\eta }_{r,1}}={{(-)}^{r}}0.1{{\rho }_{0}}$ for $-R\leqslant {{r}_{1}}\leqslant R-1$ in the background initial density ${{\psi }_{0}}=1$, in which the length of the flux is R = 10. The presence of point charges is taken into account by fixing ${{\eta }_{R,1}}=0.1{{\rho }_{0}}$ and ${{\eta }_{R-1,1}}=-0.1{{\rho }_{0}}$ through the time evolution. The free parameters of this system are $({{\gamma }^{2}},V_{0}^{\prime },J)$, related to $({{c}_{1}},{{c}_{2}},{{c}_{3}})$. By using the $\Delta \tau$-independent parameters, we expect the confinement (Higgs) phase for small (large) values of ${{c}_{1}}/{{c}_{2}}={{\gamma }^{2}}V_{0}^{\prime }$ and ${{c}_{3}}{{c}_{2}}=2\;J{{\rho }_{0}}/V_{0}^{\prime }$ (see appendix A).

Figures 2(b) and (c) represent the time evolution of the density distribution ${{\eta }^{2}}$ calculated by the previous two models. For a certain value of J, both models show similar behaviors for small values of ${{\gamma }^{2}}$, where the placed density flux is stable and does not spread out. This captures the characteristics of the confinement phase with strong phase fluctuation, where the density fluctuation can be localized by the mechanism similar to the self-trapping effects as observed in a cold atom experiment [28]. However, the underlying physics are slightly different because the system in [28] possesses only on-site interaction, without a long-range one. With increasing ${{\gamma }^{2}}$, i.e., the Higgs coupling, the structure of the density flux is gradually lost by emitting the density waves from the charge. This emission is a characteristic of the superfluid phase, i.e, Higgs phase, where the phase coherence can generate a long-wavelength phonon. 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 $\sim \sqrt{J{{\rho }_{0}}V_{0}^{\prime }}$ for ${{\gamma }^{2}}\gg 1$.

Zoom In Zoom Out Reset image size

To judge whether the system is in confinement or Higgs regime by dynamic simulations, we calculate the remnants of the flux $\sigma (t)$ defined by

$$\begin{eqnarray}&&\sigma (t)=\mathop{\sum }\limits_{\ell \in {\rm initial}\,{\rm flux}\,{\rm line}}{{[\eta _{\ell }^{2}(t)-\eta _{\ell }^{2}(0)]}^{2}},\end{eqnarray} \tag{ 9 }$$

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 2(a) shows the dynamical phase diagram obtained by the behavior of σ shown in figures 2(d) and (e). The rapid oscillation of σ reflects the periodic vanish-revival cycle of the density flux. In the BH model, we calculate the time average ${{\langle \sigma \rangle }_{t}}$ and determine the phase boundary by finding the point at which ${{\langle \sigma \rangle }_{t}}$ almost vanishes (below 0.001; see figure 2(d)). In the GH model, the boundary is determined by the appearance of rapid growth of $\sigma (t)$ due to the intermittent density-wave emission as seen in figure 2(e).

It is important to note that our dynamic approach can give a new method to explore the phase structure of the LGT. The validity of our approach exactly stems from the correspondence of the LGT to the theoretical description of the atomic systems in section 2. Although the dynamic results are obtained under the mean field approximation and are only applicable to the GH model with the unitary gauge of the Higgs field [19], the dynamic phase boundaries of both models are qualitatively in good agreement with the result of the Monte Carlo simulations of the full GH model of equation (3). (See figure 5 and appendix A.)

Zoom In Zoom Out Reset image size

The dynamic difference of the BH and GH models can be observed in the amplitude fluctuation of the simulating gauge field. Because the GH model is obtained by expanding $\rho ={{\rho }_{0}}+\eta$ around the constant density ${{\rho }_{0}}\gg 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 figure 2(a) is obtained for the particular value ${{\rho }_{0}}=1$. In addition, near the situation $\mu =0$ represented by a dotted curve in figure 2(a), the density fluctuation is accidentally frozen because the development of the homogeneous wave function is driven as ${{\psi }_{0}}{{e}^{-i\mu t}}$. Then, the dynamics of the BH model are similar to the GH model. This is a reason of the decrease of ${{\langle \sigma \rangle }_{t}}$ around ${{\gamma }^{2}}V_{0}^{\prime }=2.5$. Another point is that the amplitude fluctuation in the BH model can give rise to a similar effect of the fluctuation of the Higgs coupling. When the Higgs field moves away from the London limit, the Higgs-confinement transition may become first order and its boundary can be sharp [29]. Since our GH model corresponds to the London limit, in which the amplitude fluctuation of the Higgs field is absent, the phase boundary becomes less clear because the two phases connect with each other through crossover. The significant amplitude fluctuation in the BH model can lead to the stabilization of the Higgs phase as seen in figure 2(a).

In this section we present two methods to realize V_ab as shown in table 1. A major way to prepare intersite interactions in BH systems is to use DDI between atoms or molecules [25, 30, 31]. 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 1. 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 1, its actual fulfillment is difficult (some discussions can be seen in the system of polar molecules [32]), and, importantly, the hopping process between sites with different dipole orientations are prohibited or reduced due to the conservation of the atomic spin. We note that the bipartite structures of the nanoscale ferromagnetic islands have been proposed for realizing the right Gauss law constraint using dipolar interactions [33]. Recently, there is an interesting proposal to realize V_ab in table 1 by using the Rydberg p-states of cold atoms [34].

In section 4.1, we discuss the possibility of realizing the values of V_ab in table 1 by using the excited bands of an optical lattice, which is an alternative route to get intersite interactions [26]. In section 4.2, we discuss a system of multi-layer 2D optical lattices [35] to realize tunable DDI between atoms. The difference from the proposal in [33] is that the long-range interaction of dipoles between different layers is controlled by tuning the height of the two layers and the length of dipoles in [33], while in our case, the long-range interaction in the same layer is controlled through the mediation of atomic interaction in different layers. These proposals are within reach in current experimental techniques.

4.1. 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 band. Thus, we expect the significant intersite density-density interaction without introducing DDI between atoms [26]. To implement this scheme, we assume the following optical lattice potential:

$$\begin{eqnarray} {{V}_{{\rm OL}}} & = & {{V}_{{\rm A}}}\left( {{{\rm cos} }^{2}}kx+{{{\rm cos} }^{2}}ky \right)+{{V}_{{\rm B}}}\left[ {{{\rm cos} }^{2}}k(x-y)+{{{\rm cos} }^{2}}k(x+y) \right],\quad {{V}_{{\rm A}}},{{V}_{{\rm B}}}\geqslant 0, \\ \end{eqnarray} \tag{ 10 }$$

which can be created in a current experimental setup. For ${{V}_{{\rm B}}}/{{V}_{{\rm A}}}\gt 0.5$, the potential forms a checkerboard lattice (line graph of a square lattice [24]) and its minima are characterized by an anisotropic harmonic form as shown in figure 3(a). This anisotropy is necessary to prevent the intraband mixing dynamics. Excitation to higher orbitals can be achieved by stimulated Raman transition [36] or nonadiabatic control of the optical lattice [37, 38].

Zoom In Zoom Out Reset image size

The intersite density-density interaction is proportional to the overlap integral ${{O}_{ab}}=\int d{\bf r}|{{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\omega _{{\rm ho}}^{2}({{\alpha }^{2}}{{x}^{2}}+{{y}^{2}})/2$, w_a can be represented by the harmonic oscillator basis ${{\Phi }_{{\bf n}}}({\bf r})$, where $\hbar {{\omega }_{{\rm ho}}}=2\sqrt{{{E}_{R}}({{V}_{{\rm A}}}+2{{V}_{{\rm B}}})}$ with the recoil energy ${{E}_{R}}={{\hbar }^{2}}{{k}^{2}}/2m$ of the optical lattice and $\alpha =\sqrt{(2{{V}_{{\rm B}}}-{{V}_{{\rm A}}})/(2{{V}_{{\rm B}}}+{{V}_{{\rm A}}})}$. The band index takes ${\bf 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 1 read ${{O}_{12}}\approx {{O}_{13}}\gg {{O}_{15}}$. Figures 3(b)–(d) represent the parameter domain satisfying this condition with respect to V_A and ${{V}_{{\rm B}}}/{{V}_{{\rm A}}}$, where the amplitudes V_A and V_B of the optical lattice are precisely tunable parameters. Because of the characteristics of the potential equation (10), we can have significant overlap of the Wannier functions even for the high potential height such as ${{V}_{{\rm A},{\rm B}}}=100{{E}_{R}}$ for ${{V}_{{\rm B}}}/{{V}_{{\rm A}}}\geqslant 0.5$; see appendix B for more details. For the s orbitals the domain is limited to a narrow region (${{V}_{{\rm B}}}\sim 0.55{{V}_{{\rm A}}}$) of the parameter space. Using the p- or d-orbitals allows us to get the condition of table 1 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 first half of NNN $(J^{\prime} )$; the second half of NNN (J'') is small because of the higher potential height between the link of group (iii) as seen in figure 3(a).

Finally, we admit that for actual parameter estimation, one should also try other more realistic Wannier functions such as $\sim {{x}^{c}}{\rm exp} (-h|x|)$ [39], although the qualitative feature captured here needs no modifications.

4.2. 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 purposes we explicitly describe the idea by using a triple-layer system consisting of three 2D square optical lattices (layer ${{{\rm L}}_{{\rm A}}}$, ${{{\rm L}}_{{\rm B}}}$, ${{{\rm L}}_{{\rm C}}}$) as seen in figure 4(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., ${{{\rm L}}_{{\rm A}}}$ and ${{{\rm L}}_{{\rm B}}}$, to zero, which is discussed for realistic parameter estimation at the end of appendix C.

Zoom In Zoom Out Reset image size

The boson system on the layer ${{{\rm L}}_{{\rm A}}}$ (we call them A-bosons) is a playground of the (2+1)D U(1) GH model, which is sandwiched by B-bosons on ${{{\rm L}}_{{\rm B}}}$ and C-bosons on ${{{\rm L}}_{{\rm 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 figures 4(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-bosons and B-bosons 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 interaction between the A-bosons. These two kinds of interactions may be tuned to realize V_ab as given in table 1. We omit the DDI between the B- and C-bosons because of the large separation.

Let us focus on figure 4(d). When one projects the sites of ${{{\rm L}}_{{\rm B}}}$ onto ${{{\rm L}}_{{\rm A}}}$, their image is located in the center of each plaquette of the ${{{\rm L}}_{{\rm A}}}$ lattice. Similarly, the image of sites of ${{{\rm L}}_{{\rm C}}}$ is located in the middle of NN pairs of the ${{{\rm L}}_{{\rm A}}}$ sites. In ${{{\rm L}}_{{\rm 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

$$\begin{eqnarray} {{H}_{{\rm AB}}} & = & {{V}_{{\rm AB}}}\mathop{\sum }\limits_{k,\delta }{{\rho }_{A,k+\delta }}{{n}_{B,k}}, \\ {{H}_{{\rm AC}}} & = & {{V}_{{\rm AC}}}\mathop{\sum }\limits_{l,\tilde{\delta }}{{\rho }_{{\rm A},l+\tilde{\delta }}}{{n}_{{\rm C},l}}, \\ \end{eqnarray} \tag{ 11 }$$

where ${{\rho }_{{\rm A},k}}$ and ${{n}_{{\rm B}({\rm C}),k}}$ are boson densities at the site k and ${{V}_{{\rm AB}({\rm C})}}\lt 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}_{{\rm AB}({\rm C})}}$ as perturbation, the effective attractive interaction between A-bosons may be estimated as $\sim -V_{{\rm AB}}^{2}{{\sum }_{k,\delta ,{{\delta }^{\prime }}}}{{\rho }_{{\rm A},k+\delta }}{{\rho }_{{\rm A},k+{{\delta }^{\prime }}}}$ and $-V_{{\rm AC}}^{2}{{\sum }_{l,\tilde{\delta }}}{{\rho }_{{\rm A},l+\tilde{\delta }}}{{\rho }_{{\rm A},l-\tilde{\delta }}}$. 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 1, while the latter contributes a constant only for the group (i). Then one may fulfill the condition of V_ab in table 1. The detailed calculation of the effective interaction and the experimental feasibility are described in appendix C. Although there is a small contribution of long-range interaction beyond the NNN links due to the power-law tail of r³, this correction may suppress the density fluctuation and result in the enhancement of the confinement phase.

In conclusion, realization of the quantum simulator of the U(1) lattice GH model provides a significant innovation to tackle unresolved problems such as the inflation Universe, being able to be constructed by the cold atomic architecture. The phase structure of the atomic simulators may be explored by the non-equilibrium dynamics, where the electric flux dynamics can be observed from the behavior of the density fluctuation. We proposed two experimentally feasible schemes (Methods A and B) to respect the constraint of Gauss's law and locality of the gauge interaction in the atomic simulators.

Many works have been devoted to the dynamic properties of phase defects, namely quantized vortices, by analyzing the GP equation [40]. In terms of the gauge theory, these phase defects correspond to the magnetic fluxes. Our work focuses on the density fluxes, corresponding to electric fluxes, whose dynamics are under constraint by Gauss's law. Such a density flux in the GP model has not been discussed before, and this point of view could open the door for a new avenue of the GP dynamics, such as dynamic features of various configurations of an electric flux or many fluxes. These non-equilibrium dynamics are interesting themselves, although they could also give references as a guide not only to the atomic simulator experiments, but also to the LGT. The dynamic equations can be derived and give some insights for various models of the LGT.

The other problems for future study include the clarification of the global phase diagram of equation (3) for the general sets of parameters and of how to implement the general terms in equation (3) experimentally. It has been proposed in [19] that the Higgs coupling (${{c}_{1i}}$-term) in the spatial dimension can be implemented by using an idea of [41]. An idea to generate the spatial plaquette (${{c}_{2ij}}$-) term is discussed in [42]. There is still insufficient discussion on how to combine these schemes toward the quantum simulation of the full GH model, which is a subject for future study. Fine-tuning of the intersite density-density interaction is also an important task, and we believe that the method in section 4.1 is the most feasible scheme in actual experiments. Our method in section 4.2 provides a new scheme for tuning the intersite atom-atom interactions, and more elaborate discussion using concrete atomic species, optical lattice structures, etc, remains to be studied. All of these issues will be reported in future publications.

This work was supported by KAKENHI from JSPS (Grant Nos. 26400371, 25220711, 26400246 and 26400412).

Let us explain the phase structure of the GH model defined by equation (3) with asymmetric couplings ${{c}_{1\mu }},{{c}_{2\mu \nu }}$ given by equation (4). First, we note that the (2+1)D version of the standard 4D U(1) GH theory [20], which is considered in HEP and has the symmetric couplings (${{c}_{1\mu }}={{c}_{1}}\geqslant 0,{{c}_{2\mu \nu }}={{c}_{2}}\geqslant 0,{{c}_{3,4,5}}=0$ in equation (3)), is always in the confinement phase [43], in which the phase ${{\theta }_{x,\mu }}$ is unstable by strong fluctuation. In our model, inclusion of sufficient c₃ in addition to the asymmetric couplings ${{c}_{1\mu }}$ and ${{c}_{2\mu \nu }}$ lets the system enter into the 'Higgs' phase, where both ${{\theta }_{x,\mu }}$ and ${{\varphi }_{x}}$ are stable (see figure 5).

To identify the location of the transitions, we measure the internal energy $U=\langle A\rangle$ and the specific heat $C=\langle {{A}^{2}}\rangle -{{\langle A\rangle }^{2}}$ by using the standard Metropolis algorithm in Monte Carlo (MC) 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\times 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\sim 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₂. The result is presented in figure 5. There are two phases: the Higgs phase in the large c₃ region (upper region) and the confinement phase in the small c₃ 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 one believed 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}}\geqslant 0$ [19], the Coulomb phase is missing due to the low dimensionality.

To understand figure 5, let us consider some limiting cases. First, after choosing the unitary gauge ${{\phi }_{x}}=1$, let us consider the limit ${{c}_{1}}\to \infty$. Then the c₁ term makes ${{\theta }_{x,0}}=0$ [mod(2π)], and the action becomes

$$\begin{eqnarray} {{A}_{{{{\rm c}}_{1}}=\infty }} & = & {{c}_{2}}\mathop{\sum }\limits_{x}\mathop{\sum }\limits_{i=1}^{2}{\rm cos} ({{\theta }_{x+0,i}}-{{\theta }_{x,i}}) \\ {} & {} & +{{c}_{3}}\mathop{\sum }\limits_{x}\left[ {\rm cos} ({{\theta }_{x,1}}-{{\theta }_{x,2}})+{\rm cos} ({{\theta }_{x,1}}+{{\theta }_{x+1,2}}) \right. \\ {} & {} & +\left. {\rm cos} ({{\theta }_{x+1,2}}-{{\theta }_{x+2,1}})+{\rm cos} ({{\theta }_{x,2}}+{{\theta }_{x+2,1}}) \right], \\ \end{eqnarray} \tag{ A.1 }$$

up to constant. This is viewed as a 3D XY spin model with asymmetric couplings, where ${{\theta }_{x,i}}\ (i=1,2)$ on the link $(x,x+i)$ is the XY spin angle ${{\tilde{\theta }}_{\tilde{x}}}$. In fact, the c₃ term is their NN coupling in the 12 plane and the c₂ term is their NN coupling along the $\mu =0$ axis. The region of sufficiently large c₂ and c₃ is the ordered phase of this XY spin and corresponds to the Higgs phase with small gauge-field (${{\theta }_{x,\mu }}$) fluctuations. As a check of figure 5, let us consider the case ${{c}_{2}}={{c}_{3}}$ of equation (A.1), which reduces to the symmetric 3D XY spin model of ${{A}_{3{\rm DXY}}}={{c}_{{\rm XY}}}{{\sum }_{\tilde{x},\mu }}{\rm cos} ({{\tilde{\theta }}_{\tilde{x}+\mu }}-{{\tilde{\theta }}_{\tilde{x}}})$. It is known to have a genuine second-order phase transition at ${{c}_{{\rm XY}}}\simeq 0.45$. Therefore the transition line in figure 5(b) should approach to ${{c}_{2}}\cdot {{c}_{3}}\to {{0.45}^{2}}\simeq 0.20$ as ${{c}_{1}}/{{c}_{2}}\to \infty$ as it shows.

Next, let us consider the case ${{c}_{2}}=0$. Then, each variable ${{\theta }_{x,0}}$ appears only through the c₁ term without couplings to other variables (we take the unitary gauge as before). Then the dynamics is controlled by the c₃ term. Again, this term is viewed as the energy of the XY spins ${{\theta }_{xi}}\ (i=1,2)$. However, they have no coupling along the $\mu =0$ direction, and therefore the system is a collection of decoupled 2D XY spin models. The 2D XY spin model is known to exhibit KT transition, which is infinitely continuous. Thus, although it is not drawn in figure 5, there should be added a horizontal line (independent of c₁) for ${{c}_{2}}=0$ consisting of a collections of KT transitions at around ${{c}_{3}}\sim 0.96$. We understand that the crossover points appearing in the smaller c₁ part in each curve for three c₂ drawn in figure 5 are the remnants 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 ${{\theta }_{x,0}}$ and ${{\theta }_{x,0}}$ should determine their dynamics only through the c₂ term. Thus, even ${{\theta }_{x,i}}$ could be set constant with no fluctuations, ${{\theta }_{x,0}}$ has the NN coupling in each 12 plane. However, two dimensions is not enough to stabilize ${{\theta }_{x,0}}$. In turn, the c₂ term is not enough to sustain the coupling between ${{\theta }_{x,i}}$ along the $\mu =0$ direction. The dynamics of ${{\theta }_{x,i}}$ is essentially from the c₃ 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}}=\infty$. Then ${{\theta }_{x,\mu }}$ is frozen to be a pure gauge configuration, ${{\theta }_{x,\mu }}={{\lambda }_{x+\mu }}-{{\lambda }_{x}}$. By plugging this into the c₁ and c₃ term, we obtain

$$\begin{eqnarray} {{A}_{{{{\rm c}}_{2}}=\infty }} & = & {{c}_{1}}\mathop{\sum }\limits_{x}{\rm cos} ({{\lambda }_{x+3}}-{{\lambda }_{x}})+2{{c}_{3}}\mathop{\sum }\limits_{x}\left[ {\rm cos} ({{\lambda }_{x+1}}-{{\lambda }_{x+2}})+{\rm cos} ({{\lambda }_{x}}-{{\lambda }_{x+1+2}}) \right], \\ \end{eqnarray} \tag{ A.2 }$$

which belongs again to the class of 3D XY spin models, where ${{\lambda }_{x}}$ is the XY spin angles on the site x. So we should have a second-order transition at ${{c}_{2}}=\infty$ as long as both c₁ and c₃ are nonvanishing. This is consistent with figure 5.

Let us finally comment on the transition line of figure 5(b) and the boundaries of figure 2(a) calculated by dynamical simulation in section 3. Their behaviors in the $({{c}_{1}}/{{c}_{2}})$–$({{c}_{3}}\cdot {{c}_{2}})$ plane are qualitatively consistent but different in quantitative comparison. We understand that there is no inconsistency in these results because the two methods, MC and GP, are different in nature: MC is static and GP is dynamic, they treat fluctuations in contrasting manners, and especially, the dynamic simulations necessarily exhibit various properties of the system according to their setup and probes, etc This certainly motivates exact quantum and dynamic simulation of the BH model in experiments.

In this section, we describe the calculation of the overlap integrals discussed in section 3.1. For the horizontal links of the potential V_OL of equation (10), the minimum is approximated by the harmonic oscillator ${{V}_{hx}}=m\omega _{{\rm ho}}^{2}({{\alpha }^{2}}{{x}^{2}}+{{y}^{2}})/2$. The basis function of V_hx is given by

$$\begin{eqnarray}&&{{\Phi }_{{\bf n}}}(\sqrt{\alpha }x,y)={{A}_{{\bf n}}}{{H}_{{{n}_{x}}}}\left( \sqrt{\alpha }\frac{x}{{{a}_{{\rm ho}}}} \right){{H}_{{{n}_{y}}}}\left( \frac{y}{{{a}_{{\rm ho}}}} \right){{e}^{-(\alpha {{x}^{2}}+{{y}^{2}})/2a_{{\rm ho}}^{2}}},\end{eqnarray} \tag{ B.1 }$$

where H_n is the Hermite polynomial, ${{A}_{{\bf n}}}$ the normalization factor, and ${{a}_{{\rm ho}}}=\sqrt{\hbar /m{{\omega }_{{\rm ho}}}}$ the harmonic oscillator length.

The s, p, and d orbitals for these links correspond to ${\bf n}=({{n}_{x}},{{n}_{y}})=(0,0)$, $(1,0)$, and $(2,0)$. As the Wannier function ${{w}_{a}}({\bf r})$ at the link $(r,a)$, we use ${{\Phi }_{{\bf n}}}$ with $(x,y)$ measured from the center of the link. For the vertical links, the minimum is also approximated as ${{V}_{hy}}=m\omega _{{\rm ho}}^{2}({{x}^{2}}+{{\alpha }^{2}}{{y}^{2}})/2$ and the basis function is ${{\Phi }_{{\bf n}}}(x,\sqrt{\alpha }y)={{A}_{n}}{{H}_{{{n}_{x}}}}(x/{{a}_{{\rm ho}}}){{H}_{{{n}_{y}}}}(\sqrt{\alpha }y/{{a}_{{\rm ho}}}){{e}^{-({{x}^{2}}+\alpha {{y}^{2}})/2a_{{\rm ho}}^{2}}}$. The s, p, and d orbitals for these links correspond to ${\bf n}=({{n}_{x}},{{n}_{y}})=(0,0)$, $(0,1)$, and $(0,2)$. Then, the Wannier functions ${{w}_{a}}({\bf r})$ relevant to the following calculations are given as follows:

$$\begin{eqnarray} {{w}_{1}}({\bf r}) & = & {{\Phi }_{{\bf n}}}(\sqrt{\alpha }{{R}_{0}}(\tilde{x}-1/2),{{R}_{0}}\tilde{y}), \\ {{w}_{2}}({\bf r}) & = & {{\Phi }_{{\bf n}}}({{R}_{0}}\tilde{x},\sqrt{\alpha }{{R}_{0}}(\tilde{y}+1/2)), \\ {{w}_{3}}({\bf r}) & = & {{\Phi }_{{\bf n}}}(\sqrt{\alpha }{{R}_{0}}(\tilde{x}+1/2),{{R}_{0}}\tilde{y}), \\ {{w}_{5}}({\bf r}) & = & {{\Phi }_{{\bf n}}}(\sqrt{\alpha }{{R}_{0}}(\tilde{x}-1/2),{{R}_{0}}(\tilde{y}-1), \\ \end{eqnarray} \tag{ B.2 }$$

where R₀ represents the lattice spacing and we shift the origin of the coordinate to $({{R}_{0}}/2,{{R}_{0}}/2)$ of figure 3(a). The length scale of the coordinate is normalized by R₀ and the dimensionless coordinates are denoted by putting tildes.

The intersite interaction strength V_ab is proportional to the overlap integrals ${{O}_{ab}}=\int d{\bf r}|{{w}_{a}}{{|}^{2}}|{{w}_{b}}{{|}^{2}}$.

It is sufficient to calculate only the integrals for the link pairs $(a,b)=(1,2)$, $(1,3)$, and $(1,5)$, because ${{O}_{12}}={{O}_{23}}={{O}_{34}}={{O}_{41}}$, ${{O}_{13}}={{O}_{24}}$, and ${{O}_{15}}={{O}_{26}}$ due to the lattice symmetry.

The typical results of the overlap integrals for the three orbitals are shown in figure B1 for ${{V}_{{\rm B}}}/{{V}_{{\rm A}}}=0.6$ as a function of V_A. We also show the integral for the on-site contribution ${{O}_{11}}=\int d{\bf r}|{{w}_{a}}{{|}^{4}}$ and the hopping integrals $J=\int d{\bf r}{{w}_{1}}(-{{\hbar }^{2}}{{\nabla }^{2}}/2m+{{V}_{{\rm OL}}}){{w}_{2}}$, $J^{\prime} =\int d{\bf r}{{w}_{1}}(-{{\hbar }^{2}}{{\nabla }^{2}}/2m+{{V}_{{\rm OL}}}){{w}_{3}}$, and $J^{\prime\prime} =\int d{\bf r}{{w}_{1}}(-{{\hbar }^{2}}{{\nabla }^{2}}/2m+{{V}_{{\rm OL}}}){{w}_{5}}$. In any case, ${{O}_{11}}$ is monotonically increased with V_A, and J'' and ${{O}_{15}}$ (not seen in figure B1) are negligibly small. In the case of the s-orbital, ${{O}_{12}}$ and ${{O}_{13}}$ are also monotonically decreasing functions, so that the range satisfying ${{O}_{12}}\approx {{O}_{13}}$ is only limited by a narrow range or a point with respect to V_A. On the other hand, for the p- and d-orbitals ${{O}_{12}}$ and ${{O}_{13}}$ change non-monotonically because of the node structure and the extended amplitude profile of the wave functions. This fact extends the range of ${{O}_{12}}\approx {{O}_{13}}$ as seen in figures B1(b) and (c).

Zoom In Zoom Out Reset image size

Note that the hopping integrals J and J' are of $\mathcal{O}(1)$ even for ${{V}_{{\rm A}}}=100{{E}_{R}}$. This is because the energy barrier between links of groups (i) and (ii) in table 1 is the sub-maximum with the height $2{{V}_{{\rm B}}}-{{V}_{{\rm A}}}$ at $({{R}_{0}}/2,{{R}_{0}}/2)$ in figure 3(a). Since the value of $J(J^{\prime} )$ is bigger than that of ${{O}_{12(13)}}$ by two orders of magnitude, one needs to increase considerably the s-wave scattering length via a Feshbach resonance to get the exact Gauss's law constraint, namely, ${{V}_{ab}}\gg J$.

In this section, we apply the second-order perturbation theory to the triple-layer system in section 4.2 to derive the effective intersite interaction of A-bosons, and estimate the possible values of involved parameters to realize V_ab of table 1. 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 A- and B-bosons (two layers ${{{\rm L}}_{{\rm A}}}$ and ${{{\rm L}}_{{\rm B}}}$), which has the NN DDI, ${{H}_{{\rm AB}}}$ of equation (11). It implies that the B-boson on the site k interacts with the four NN A-bosons on the sites $k\pm {{\delta }_{x(y)}}$ as seen in figure 4(d). ${{V}_{{\rm AB}}}$ in ${{H}_{{\rm AB}}}$ is expressed as

$$\begin{eqnarray} {{V}_{{\rm AB}}} & = & \int d{\bf r}d{{{\bf r}}^{\prime }}{{U}_{{\rm DD}}}({\bf r},{{{\bf r}}^{\prime }})|{{w}_{{\rm A}}}({\bf r}){{|}^{2}}|{{w}_{{\rm B}}}({{{\bf r}}^{\prime }}){{|}^{2}}, \\ {{U}_{{\rm DD}}}({\bf r},{{{\bf r}}^{\prime }}) & = & \frac{C}{4\pi |{\bf r}-{{{\bf r}}^{\prime }}{{|}^{3}}}\left( 1-\frac{3\ell _{{\rm AB}}^{2}}{|{\bf r}-{{{\bf r}}^{\prime }}{{|}^{2}}} \right), \\ \end{eqnarray} \tag{ C.1 }$$

where ${\bf r}({{{\bf r}}^{\prime }})$ is the position of A(B)-boson, ${{w}_{{\rm A}({\rm B})}}({\bf r})$ is their Wannier function, and $C={{\mu }_{0}}{{\tilde{\mu }}_{{\rm A}}}{{\tilde{\mu }}_{{\rm B}}}$; ${{\mu }_{0}}$ is the magnetic permeability of the vacuum and ${{\tilde{\mu }}_{{\rm A}({\rm B})}}$ is the magnetic moment of A(B)-atoms.

We assume that the B-bosons of ${{{\rm L}}_{{\rm B}}}$ have a chemical potential ${{\mu }_{{\rm B}}}$(>0), a negligibly small NN hopping amplitude due to a deep trapping potential, an on-site repulsion ${{U}_{{\rm B}}}(\gt 0)$, and the NN DDI with A-bosons ${{V}_{{\rm AB}}}$. One may forget the DDI between B-bosons, because it is a constant due to negligible NN hopping. Then, the Hamiltonian ${{\hat{H}}_{{\rm B}}}$ and the partition function ${{Z}_{{\rm B}}}={\rm Trexp}(-\beta {{\hat{H}}_{{\rm B}}})$ for the subsystem of B-bosons are written by using the B-boson density operator ${{\hat{n}}_{k}}$ at the site k as

$$\begin{eqnarray} {{{\hat{H}}}_{{\rm B}}} & = & \mathop{\sum }\limits_{k}\left[ -{{\mu }_{{\rm B}}}{{{\hat{n}}}_{k}}+{{U}_{{\rm B}}}{{{\hat{n}}}_{k}}({{{\hat{n}}}_{k}}-1)+{{V}_{{\rm AB}}}\mathop{\sum }\limits_{\delta =\pm {{\delta }_{x,y}}}{{{\hat{n}}}_{k}}{{{\hat{\rho }}}_{k+\delta }} \right], \\ {{Z}_{{\rm B}}} & = & \mathop{\prod }\limits_{k}{{z}_{{\rm B},k}},\ {{z}_{{\rm B},k}}=\mathop{\sum }\limits_{n=0}^{\infty }{\rm exp} \left[ -\beta ({{E}_{n}}+n{{V}_{{\rm AB}}}\mathop{\sum }\limits_{\delta }{{{\hat{\rho }}}_{k+\delta }}) \right], \\ {{E}_{n}} & = & -{{\mu }_{{\rm B}}}n+{{U}_{{\rm B}}}n(n-1). \\ \end{eqnarray} \tag{ C.2 }$$

By assuming ${{\mu }_{{\rm B}}},{{U}_{{\rm B}}}\gg {{V}_{{\rm AB}}}$, we expand ${{z}_{{\rm B},k}}$ up to $O(V_{{\rm AB}}^{2})$,

$$\begin{eqnarray*} {{z}_{{\rm B},k}} & = & \mathop{\sum }\limits_{n=0}^{\infty }{{e}^{-\beta {{E}_{n}}}}\left[ 1-n{{W}_{k}}+\frac{1}{2}{{n}^{2}}{{({{W}_{k}})}^{2}}\;\cdots \right] \\ {} & = & {{F}_{0}}\left[ 1-\frac{{{F}_{1}}}{{{F}_{0}}}{{W}_{k}}+\frac{1}{2}\frac{{{F}_{2}}}{{{F}_{0}}}{{({{W}_{k}})}^{2}}\;\cdots \right], \\ \end{eqnarray*}$$

$$\begin{eqnarray} {{W}_{k}} & = & \beta {{V}_{{\rm AB}}}\mathop{\sum }\limits_{\delta }{{{\hat{\rho }}}_{k+\delta }},\quad {{F}_{m}}\equiv \mathop{\sum }\limits_{n=0}^{\infty }{{n}^{m}}{{e}^{-\beta {{E}_{n}}}}. \\ \end{eqnarray} \tag{ C.3 }$$

Then we have

$$\begin{eqnarray} {{Z}_{{\rm B}}} & = & {{({{F}_{0}})}^{{{L}^{2}}}}{\rm exp} \left\{ \mathop{\sum }\limits_{k}\left[ {{\langle -n\rangle }_{0}}{{W}_{k}}+\frac{1}{2}{{(\Delta n)}^{2}}{{({{W}_{k}})}^{2}}+\;\cdots \right] \right\}, \\ {{\langle {{n}^{m}}\rangle }_{0}} & \equiv & \frac{{{F}_{m}}}{{{F}_{0}}},\quad {{(\Delta n)}^{2}}\equiv {{\langle {{n}^{2}}\rangle }_{0}}-\langle n\rangle _{0}^{2}. \\ \end{eqnarray} \tag{ C.4 }$$

The first-order terms $\propto \;{{W}_{k}}$ are renormalized to the chemical potential of A-bosons, and the second-order terms define the effective density-density interaction Hamiltonian ${{\hat{H}}_{{\rm ABA}}}$ of A-bosons induced by B-boson density fluctuation ${{(\Delta n)}^{2}}$,

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm ABA}}}=-\frac{\beta }{2}{{(\Delta n)}^{2}}V_{{\rm AB}}^{2}\mathop{\sum }\limits_{k}\mathop{\sum }\limits_{\delta ,{{\delta }^{\prime }}}{{\hat{\rho }}_{k,\delta }}{{\hat{\rho }}_{k,{{\delta }^{\prime }}}}.\end{eqnarray} \tag{ C.5 }$$

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, ${{\hat{H}}_{{\rm ACA}}}=-(\beta /2){{(\Delta n^{\prime} )}^{2}}V_{{\rm AC}}^{2}{{\sum }_{l,\bar{\delta },{{{\bar{\delta }}}^{\prime }}}}{{\hat{\rho }}_{l+\bar{\delta }}}{{\hat{\rho }}_{l+{{{\bar{\delta }}}^{\prime }}}},$ where ${{(\Delta n^{\prime} )}^{2}}$ is obtained by replacing ${{\mu }_{{\rm B}}},{{U}_{{\rm B}}}$ by ${{\mu }_{{\rm C}}},{{U}_{{\rm C}}}$ in ${{(\Delta n)}^{2}}$.

The sum ${{\hat{H}}_{{\rm ABA}}}+{{\hat{H}}_{{\rm ACA}}}$ contributes to the coefficients V_ab of the intersite density-density interactions for A-bosons as follows:

$$\begin{eqnarray*}&&\begin{array}{ccccccccccccccc} {\rm For}\;{\rm NN}\;{\rm links}\;{\rm (group}\;{\rm (i)}\;{\rm in}\;{\rm table}\;{\rm 1)}, \\ {{V}_{ab}}=V-\beta {{(\Delta n)}^{2}}V_{{\rm AB}}^{2}-\beta {{(\Delta n^{\prime} )}^{2}}V_{{\rm AC}}^{2}, \\ \end{array}\end{eqnarray*}$$

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} {\rm For}\;{\rm NNN}\;{\rm links}, \\ {{V}_{ab}}=\left\{ \begin{array}{ccccccccccccccc} V^{\prime} & {\rm for}\;{\rm group}\;{\rm (ii)} \\ V^{\prime} -\beta {{(\Delta n)}^{2}}V_{{\rm AB}}^{2} & {\rm for}\;{\rm group}\;{\rm (iii)}, \\ \end{array} \right. \\ \end{array}\end{eqnarray} \tag{ C.6 }$$

where V and V' are the direct DDI for NN and NNN link pairs, respectively. The condition for V_ab in table 1 can be established by adjusting two inter-layer distances ${{\ell }_{{\rm AB}}}$ and ${{\ell }_{{\rm AC}}}$ and density fluctuations ${{(\Delta n)}^{2}}({{\mu }_{{\rm B}}},{{V}_{{\rm B}}},\beta )$ and ${{(\Delta n^{\prime} )}^{2}}({{\mu }_{{\rm C}}},{{V}_{{\rm C}}},\beta )$ as

$$\begin{eqnarray} \beta {{(\Delta n^{\prime} )}^{2}}V_{{\rm AC}}^{2} & = & V-2V^{\prime} , \\ \beta {{(\Delta n)}^{2}}V_{{\rm AB}}^{2} & = & V^{\prime} (={{\gamma }^{-2}}). \\ \end{eqnarray} \tag{ C.7 }$$

Let us present some brief account for an example and estimation of the experimental parameters that satisfy the tuning relations equation (C.7). We shall report detailed discussion on this example and related topics in a future publication.

For bosons loaded in each layer, we consider ⁵²Cr atoms [44] as A bosons, $^{87}{\rm Rb}$ atoms as B bosons, and $^{168}{\rm Er}$ atoms [45] as C bosons. They have the permanent magnetic moments $6{{\mu }_{{\rm BM}}},{{\mu }_{{\rm BM}}}$, and $7{{\mu }_{{\rm BM}}}$ (${{\mu }_{{\rm 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 earlier by choosing ${{\ell }_{{\rm 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 the ⁸⁷Rb atom is small.

The method to make such a double-layer system is sketched in figure C1 . First, one prepares the 3D layer system as shown in figure C1(i) by emitting three standing waves with the wavelengths satisfying $2{{\lambda }_{1}}=\sqrt{2}{{\lambda }_{2}}={{\lambda }_{3}}$ (e.g., ${{\lambda }_{1}}=410$ nm, ${{\lambda }_{2}}=580\;{\rm nm}$ and ${{\lambda }_{3}}=820$ nm) 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 ${{\lambda }_{z}}$ to establish the 3D structure. Because ⁵²Cr,⁸⁷Rb and ¹⁶⁸Er exhibit the specific strong absorptions of photons with wave length 425 nm, 780 nm and 401 nm, respectively, above standing waves load these atoms to the sites of the corresponding layer ${{{\rm L}}_{{\rm A},{\rm B},{\rm C}}}$ [46]. This completes the step (i) in figure C1.

Zoom In Zoom Out Reset image size

In the second step (ii) in figure C1, 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 [47] successively. This achieves an effective double-layer system with ${{\ell }_{{\rm AB}}}=0,{{\ell }_{{\rm AC}}}={{\ell }_{z}}$.

Finally, let us estimate the parameters to satisfy the tuning relation equation (C.7). By making a straightforward calculation using DDI, we find that the following is a typical example of the parameters:

$$\begin{eqnarray} V & \sim & \frac{36{{\mu }_{0}}\mu _{{\rm BM}}^{2}}{4\pi \lambda _{2}^{3}},\quad {{\ell }_{{\rm AC}}}={{\ell }_{z}}\sim 580\ [{\rm nm}], \\ \beta & \sim & \frac{1}{2V},\quad {{U}_{{\rm B}}}\sim 0.3V,\quad {{\mu }_{{\rm B}}}\sim 2.5V, \\ {{U}_{{\rm C}}} & \sim & 1.3V,\quad {{\mu }_{{\rm C}}}\sim 2V. \\ \end{eqnarray} \tag{ C.8 }$$

The average densities per site are ∼560 for B-bosons and ∼8 for C-bosons. The ratio $|{{V}_{{\rm AB}({\rm C})}}/{{\mu }_{{\rm B}({\rm C})}}|$ is ∼0.192(∼0.585), which seems to validate the perturbation theory.