Quantum simulation of $(1+1)$-dimensional U(1) gauge-Higgs model on a lattice
by cold Bose gases

To understand physics of a given quantum many-body system in a coherent manner is not an easy matter. For example, the time development of these systems have not been clarified enough because high dimensionality of the state-vector space brings great difficulty for solving many-body Schrödinger equation Noack . The sign problem also prevents us from carrying out straightforward Monte Carlo simulations for models including fermionic matter fields Troyer .

In the last several years, as an alternative approach, quantum simulations of various quantum systems by using ultra-cold atomic gases on an optical lattice have attracted lots of interests coldatoms . This comes mainly from high controllability and versatility of such atomic systems. Quantum simulations of systems under a synthetic gauge field have already become feasible due to experimental developments mg_optical . Furthermore, recent studies have proposed several ideas for realizing simulators of quantum systems involving dynamical gauge fields WieseZohar . Gauge systems appear in various fields of physics as an effective model, and to understand their dynamics is one of the most important problems. Also, proposals for an experimentally feasible atomic simulator of gauge systems still remain as an open problem.

In the previous papers U1GHM ; NJP1 ; Future , some of the present authors discussed how to simulate lattice gauge-Higgs models (GHMs) Fradkin ; ichimaturev by bosonic atoms on a $D$-dimensional optical lattice ($D=3$ in Refs. U1GHM ; Future and $D=2$ in Ref. NJP1 ). The starting model in that approach is a cold atomic system described by the extended Bose-Hubbard model (EBHM) Baier ; Dutta ; a lattice model of bosons having both on-site and off-site interactions. The target gauge model, GHM, is a U(1) lattice gauge theory wilson involving a compact U(1) gauge field and a Higgs field in the London limit. Reflecting the nonrelativistic nature of the EBHM, interactions in the GHMs are asymmetric in the space and time directions; some couplings in the spatial directions have no partners in the time direction in contrast to the lattice gauge theory considered in high-energy physics wilson . However, the GHMs have a very interesting phase diagram as we shall see, and quantum simulation of the GHMs is expected to give us a very useful insight on the dynamics of the gauge theory.

In this approach U1GHM ; NJP1 ; Future , the phase degree of freedom $\hat{\theta}_{a}$ of boson operator $\hat{\psi}_{a}=\exp(i\hat{\theta}_{a})\sqrt{\hat{\rho}_{a}}$ on the site $a$ of the optical lattice plays the role of the compact U(1) gauge field $\hat{\theta}_{r,i}$ defined on the link $(r,r+\hat{i})$ of a $D$-dimensional (spatial) lattice of lattice gauge theory ($i=1,\cdots,D$ is the direction index). Compared with the alternative approach called quantum link model (gauge magnet) quantumlink , this approach has an advantage to express a compact U(1) gauge theory in terms of the genuine U(1) operator $\exp(i\hat{\theta}_{r,i})$ advantage . The Hamiltonian of the EBHM is known to contain an annoying term that breaks Gauss law of the pure gauge theory (i.e., theory containing no matter fields) and forces one to make a severe fine-tuning of its coefficient to generate a local gauge symmetry WieseZohar . However, inclusion of the Higgs field converts this breaking term to a term expressing the genuine Gauss-law with Higgs charge and relieves us from the fine-tuning.

In this paper, we consider theoretically the atomic quantum simulation of the GHM with $D=1$ , supposing a system of Bosonic atoms on a one-dimensional (1D) optical lattice. This work is motivated by the following reasons:

(i) From experimental point of views, quantum simulation of gauge theories in a 1D spatial lattice is expected to be easier than that in higher-dimensional lattices WieseZohar . The present approach is based on a set of single-component atoms with long-range interaction (e.g., dipole-dipole interaction) loaded on a 1D optical lattice. It has the following advantages advantage ; first, this set up is proved to be feasible experimentally 1Doptical_DDI . Second, a digital quantum simulation of the Schwinger model has been realized recently by using a 1D chain of several sites of trapped ions Martinez . In contrast, the present 1D optical lattice is scalable, i.e., hundreds of sites with large number of atoms can be prepared. Third, some approaches such as Ref. Bazavov use the quantum link model, whose realization requires two-component atoms on the optical lattice. The single-component atoms in the present approach is certainly more feasible in experiments.

(ii) We can expect nontrivial behavior of the phase transitions of the target GHM, not found in the higher-dimensional cases U1GHM ; NJP1 ; Future . As shown later, the phase transition is governed by the Berezinskii-Kosterlitz-Thouless (BKT) type transition BKT .

(iii) Non-equilibrium physics of the lattice gauge theories has not been deeply understood yet. This subject is closely related to the recent heavy-ion collision experiments in high-energy physics HIC . In our previous paper NJP1 ; Future , we applied the Gross-Pitaevskii (GP) model GP to study the real time dynamics of an electric flux, showing that the dynamics qualitatively realizes the behavior expected in each phase of the GHM. On the other hand, since the quantum fluctuations become more important in 1D system, we should include some quantum correction to the GP approach. Actually, real-time dynamics of a lattice gauge model has been studied by tensor network method TNN , which is mainly applicable for 1D quantum systems with low filling. Here, we use the truncated Wigner approximation TWA to study the effects of quantum fluctuations.

Let us explain the basic steps of the present paper. First, by applying the method in our previous papers U1GHM ; NJP1 ; Future to the 1D EBHM for large fillings of atoms (large average occupation numbers of atoms per site) we obtain the (1+1)D GHM with short-range interactions. Here the additional +1D in the GHMs represents the imaginary-time axis in path-integral quantization for which we also use a lattice regularization. Next, we calculate the phase diagram of the derived GHM, which contains the Higgs phase and the confinement phase separated by the BKT transition. This phase diagram is compared with that of the EBHM to discuss the interesting relations between the strongly-correlated bosons and the gauge systems; namely, the Higgs phase corresponds to the superfluid phase and the confinement phase to Mott-insulator phase. This also provides an estimation of the relevant parameter region to simulate the phase transitions in experiments. Finally, we study the real-time dynamics of the GHM by using the GP equation GP and the truncated Wigner approximation TWA . This mean-field approach naturally arises because of the assumption of the large fillings of atoms. Numerical results show that a prepared electric flux exhibits quite different behaviors depending on the Higgs or confinement phase. A phenomenon similar to the Schwinger mechanism is observed numerically.

The paper is organized as follows. In Sec. II, we present a derivation from the experimentally realized model, the EBHM on the 1D optical lattice, to the target simulated model, the (1+1)D U(1) GHMs which follows our previous works U1GHM ; Future . In Sec. III we discuss the phase diagram of the GHM calculated by the Monte Carlo simulations. The obtained phase diagram is compared with that of the EBHM, whose details of the calculation is described in Appendix A. In Sec. IV, we study non-equilibrium physics of the GHM, especially focusing on real-time dynamics of an electric flux. We present conclusion in Sec. V

We start with the EBHM defined on a 1D optical lattice (see Fig. 1). The hamiltonian $H_{\rm EBH}$ is given as

where $\hat{\psi}^{\dagger}_{a}$ and $\hat{\psi}_{a}$ are creation and annihilation operators of bosonic atoms on the site $a$, respectively, satisfying $[\hat{\psi}_{a},\hat{\psi}^{\dagger}_{a^{\prime}}]=\delta_{aa^{\prime}}$, etc. By representing $\hat{\psi}_{\alpha}=e^{i\theta_{\alpha}}\sqrt{\hat{\rho}_{\alpha}}$, we have the phase operator $\hat{\theta}_{\alpha}$ and the density operator $\hat{\rho}_{\alpha}$. The coefficient $J$ represents the hopping strength, $U$ the on-site interaction, and $V(>0)$ the nearest-neighbor (NN) replusive interaction generated by, e.g., a dipole-dipole interaction DDI . The above on-site repulsion $U(>0)$ represents the sum of $s$-wave scattering interaction $U_{s}$ and on-site dipole-dipole interaction $U_{d}$; $U=U_{s}+U_{d}$. The $s$-wave scattering amplitude $U_{s}$ is highly controllable by Feshbach resonance Feshbach and the hopping $J$ is controlled by the strength of the laser that makes the optical lattice. On the other hand, the NN repulsion $V$ is determined by the kind of the atoms loaded on the lattice and also depends on the lattice spacing. In particular, ratios $J/U$ and $V/U$ can become highly-controllable parameters in real experiments by changing lattice depth coldatoms and $s$-wave scattering lengths Feshbach , respectively.

The procedure U1GHM ; NJP1 ; Future to derive the GHM with local interactions from the EBHM for homogeneous density average at large fillings is based on a couple of assumption; (i) the average density is homogeneous and large, $\rho_{0a}(=\langle\hat{\rho}_{a}\rangle)=\rho_{0}\gg 1$ (large filling) and (ii) the density fluctuation around $\rho_{0}$ is small. The point (ii) implies that we focus on the sector of low-energy excitations. Explicitly, we write the density operator as

and expand $H_{\rm EBH}$ with respect to the fluctuation operator $\hat{\eta}_{a}$ up to the second-order. This kind of expansion is widely used in atomic physics coldatoms .

To regard this expanded Hamiltonian as that of the GHM, we introduce a gauge lattice and a set of operators in gauge theory defined on it. The gauge lattice is the dual lattice of the optical lattice (see Fig. 1); the site $r$ of the gauge lattice sits on the middle point of the link $(a,a+1)$ of the optical lattice. We define the 1st component of the vector potential $\hat{\theta}_{r,1}$ and its conjugate momentum, namely the electric field $\hat{E}_{r}$, on the gauge link $(r,r+1)$ as

where we use the same character $\theta$ both for atomic and gauge variables (They can be distinguished by whether it carries the direction suffix 1). Strictly speaking, the electric field should carry the direction suffix 1 as $\hat{E}_{r,1}$, but we omit it in this paper because there are no other components in a 1D system. On the other hand, we keep the suffix 1 in $\hat{\theta}_{r,1}$ because we shall introduce the zeroth component of vector potential later [see $\theta_{x,0}$ in Eq. (10)]. These pairs satisfy the canonical commutation relations wilson ,

which come from the commutation relations of $\hat{\psi}_{a}$. The factor $(-)^{r}$ in Eq. (3) may look curious but is necessary to generate the Gauss law as we shall see. Then the expanded Hamiltonian up to $O(\eta^{2})=O(E^{2})$ is expressed as

where we omitted hat symbols for operators. There appear no terms linear in $E_{r}$ because $\rho_{0}$ is chosen such that it gives rise to a stationary point of the energy. $C_{\rho_{0}}$ is an irrelevant constant. The hopping $J$ term of Eq. (1) generates both the last term of $H_{\rm GH}$ and the term $\Delta H$. Their coefficients are $J\rho_{0}$ and $J/\rho_{0}$ respectively, and the latter is suppressed by an extra factor $\rho_{0}^{-2}$ for large fillings. Therefore we neglect $\Delta H$ and keep $H_{\rm GH}$ as the Hamiltonian of the GHM for large fillings.

The first term of $H_{\text{GH}}$ leads to Gaussian with respect to $\nabla_{1}E_{r}\equiv E_{r+1}-E_{r}$, the divergence of the electric field, and gives rise to $\langle(\nabla_{1}E_{r})^{2}\rangle\simeq V^{-1}$. This becomes the Gauss-law constraint $\nabla_{1}E_{r}=0$ for pure gauge theory in the limit $V\to\infty$, and this limit is just the fine tuning required for the pure gauge theory in Refs. WieseZohar . After introduction of a Higgs field in the London limit, this term represents the genuine Gauss law ($\nabla_{1}E_{r}=Q_{\rm Higgs}$) with the Higgs charge $Q_{\rm Higgs}$ (see Ref. Future for details).

The second term represents the energy of electric field. In the ordinary lattice gauge theory used in high-energy physics wilson , its coefficient $g^{2}$ is the square of so called real gauge coupling constant $g$ and $g^{2}$ is treated as positive definite. In our case, $g^{2}=U-2V$ may take negative value (although we used a suggestive notation $g^{2}$ that it may be positive). Here we can see that a homogeneous configuration of $\rho_{0a}=\rho_{0}$ is realized for $g^{2}>0$ by using a very simple type of the mean-field theory, i.e., two site mean-field theory. We take only two nearest-neighbor site of the Hamiltonian $E_{\rm{EBH}}$, called even-site and odd-site. By putting a mean field ansatz $\psi_{\rm even}\rightarrow\sqrt{\rho_{\rm even}}$ and $\psi_{\rm odd}\rightarrow\sqrt{\rho_{\rm odd}}$ into the two sites, the EBHM mean field energy $E_{\rm{EBH}}$, which depends on the two mean density is given as

Then by detecting the energy minimum under a canonical ensemble constraint, $\langle\rho_{a}\rangle=\rho_{0}$ $(\rho_{\rm even}+\rho_{\rm odd}=2\rho_{0})$, we can show the homogeneity of the atom density. Fig. 2 indicates the homogeneity of the atom density in ($U$-$V$)-plane. Here, we found that, $\rho_{0a}$ is homogeneous for $g^{2}>0$, i.e., $\rho_{0a}=\rho_{0}$, while for $g^{2}\ll 0$, $\rho_{0a}$ takes an inhomogeneous DW configuration such as $\rho_{0a}=\rho_{0}+(-)^{a}\delta\rho/2$ (see Fig. 2). Result in Fig. 2 also indicates the existence of a homogeneous ‘phase’ even for $g^{2}<0$. This comes from the fact that the ‘critical line’ $g^{2}=U-2V=0$ is renormalized by the hopping $J$-term in $H_{\rm EBH}$ in Eq. (1). In Appendix A, we give rather detailed discussion on this homogeneous state and obtain the global phase diagram of the EBHM.

The third $J\rho_{0}$ term of $H_{\rm GH}$ expresses the NN coupling of the U(1) gauge field $U_{r,1}=\exp(i\theta_{r,1})$ as $J\rho_{0}U_{r+1,1}U_{r,1}$+ H.c., but it breaks the gauge symmetry under $\theta_{r,1}\to\theta_{r,1}+\lambda_{r+1}-\lambda_{r}$ ($\lambda_{r}$ is an arbitrary $r$-dependent real parameter) wilson . In Refs. U1GHM ; NJP1 ; Future we introduced a U(1) Higgs field $\phi_{r}=\exp(i\varphi_{r})$ in the London limit (its amplitude is frozen to unity) on the site $r$, and regard this term as a gauge-fixed version of the Higgs coupling to the unitary gauge $\phi_{r}=1$ as

This Higgs coupling and the remaining terms in $H_{\rm GH}$ is invariant under the local U(1) gauge transformation wilson

where the bar symbol in $\bar{V}_{r}$ denotes the complex conjugate. In short, $\tilde{H}_{\rm GH}(\phi,U,E)\equiv H_{\rm GH}(\phi U\phi^{\dagger},E)$ is a gauge-invariant Hamiltonian which we formally work with and the gauge invariant observables are calculable by its gauge fixed version $H_{\rm GH}(U,E)$ of Eq. (5).

Here we comment that, in the present 1D system, the Gauss law term of gauge theory is to be supplied only by the $U$- and $V$-terms in $H_{\rm EBH}$. This is in contrast with the higher dimensional [(2+1)D and (3+1)D] systems in which one must take additional next-NN interactions into account NJP1 ; Future . This point shows that the 1D atomic system is more feasible to realize a quantum simulator in experiments than the higher dimensional systems.

To derive the gauge model, we start with the partition function $Z_{\rm GH}={\rm Tr}\exp(-\beta\tilde{H}_{\rm GH}(\phi,U,E))$ defined from $H_{\rm GH}$ of Eq. (5) with the inverse temperature $\beta=1/(k_{\rm B}T)$ expressed in path-integral in canonical formalism (sum over canonical pair of coordinate and momentum). Then all the fields are defined on the (1+1)D lattice with its site denoted by $x=(x_{0},x_{1})\equiv(\tau,r)$ as $\phi_{x},U_{x,1},E_{x}$, where $\tau=0,\cdots,N-1$ is the site index along the imaginary-time axis discretized by the spacing $\Delta\tau\equiv\beta/N$ and we take the limit $N\to\infty$ formally. In this paper we omit the hat symbol in $\hat{\mu}\ (\mu=0,1)$ to denote the unit vector in the $\mu$-th direction. Therefore $x+0$ and $x+1$ imply the NN site $(\tau+1,r)$ and $(\tau,r+1)$ respectively. We use $\nabla_{\mu}$ as the forward difference operator $\nabla_{\mu}f_{x}\equiv f_{x+\mu}-f_{x}$.

We introduce the imaginary-time ($\tau$) component of the gauge field $\theta_{x,0}$ on the link $(x,x+0)$ to make the first Gauss-law term of $H_{\rm GH}$ a linear form in $\nabla_{1}E_{x}$, i.e.,

Because every step in the approach of Refs. U1GHM ; NJP1 ; Future is applied in the present (1+1)D system in a straightforward manner, we present the expressions for $Z_{\rm GH}$ as (For details see Sec. II of Ref. Future ),

The action $A_{\rm GH}$ may be expressed in terms of the compact U(1) variables $U_{x,\mu}=\exp(i\theta_{x,\mu})$ on the link $(x,x+\mu),\mu=0,1$ $[x=(x_{0},x_{1})=(\tau,r)]$ and $\phi_{x}=\exp(i\varphi_{x})$ as

The parameters in $A_{\rm GH}$ are expressed in terms of those in $H_{\rm GH}$ and $\Delta\tau$ as

The action $A_{\rm GH}$ in Eq. (11) contains only the short-range interaction, and each term is obviously invariant under a time-dependent local gauge transformation [given by Eq. (8) with the replacement $r\to x=(\tau,r)$], and the term $A_{I}$ is nothing but the gauge-invariant kinetic term of the Higgs field $\phi_{x}$ Fradkin . The present action (11) is derived from the nonrelativistic Hamiltonian (1) and has no relativistic invariance under change of directions $\mu=0\leftrightarrow\mu=1$, which is in contrast with the conventional models wilson in lattice gauge theory. However, this model supports both the confinement and Higgs phases (see Sec. III) and therefore is interesting enough for simulation of the gauge theory.

In this section we calculate the phase diagram of the GHM defined by Eqs.(10) and (11) by MC simulations. For ordinary second-order phase transitions, one may use the specific heat $C_{\rm GH}$,

to locate the transition point by the peak of $C_{\rm GH}$. However, because the present GHM is a (1+1)D system with a global U(1) symmetry as we shall see, the general theorem allows neither second-order transition nor long-range order. One may expect, if any, a phase transition of the same type as the BKT transition that is well known in the 2D XY spin model BKT . This expectation is plausible because the original 1D EBHM itself is known to have a phase transition at zero temperature between the Mott-insulator and superfluid, where the superfluid has a quasi-long range order due to the low dimensionality of the space.

Let us summarize the 2D XY spin model and its BKT transition. Its action is given by

where $\sigma_{x}\in[0,2\pi)$ is the XY spin angle defined on $x$ and $\vec{S}_{x}\equiv(\cos\sigma_{x},\sin\sigma_{x})$. The BKT transition of the 2D XY model is of infinite order, signaled by divergence of the susceptibility $\chi$,

at the transition point $J=J_{c}\simeq 1.12$. This value deviates from the round peak of specific heat by $\sim 20\%$. The spin-spin correlation function $G(x)\equiv\langle\vec{S}_{x}\vec{S}_{0}\rangle$ for large $|x|$ behaves exponentially $\exp(-|x|/\xi)$ in the spin disordered ($J<J_{c}$) phase and algebraically $|x|^{-\eta}$ in the spin quasi-ordered ($J>J_{c}$) phase. The correlation length $\xi$ and the critical exponent $\eta$ depend on $J$ as $\xi(J)\propto\exp(b/\sqrt{J_{c}-J})$ near $J_{c}$ and $\eta(J_{c})=1/4$.

To locate the phase transition point, one may observe the scaling behavior of $\chi$ of Eq. (15) for lattice of size $L\times L$ with the periodic boundary condition. In the quasi-ordered phase, power decay of $G(x)$ gives rise to $\chi\propto L^{2-\eta(J)}$ for large $L$. One may define the scaled susceptibility density as $\bar{\chi}(J,\eta^{*})\equiv L^{\eta^{*}-2}\chi(J)$ with arbitrary parameter $\eta^{*}$. Then one has

Therefore, at $\eta^{*}=\eta(J)$, $\bar{\chi}(J,\eta(J))\propto L^{0}$ becomes scale invariant (no $L$-dependence). We note that this scale invariance holds for every point of $J(>J_{c})$. Among possible multiple set of scale-invariant points $(J,\eta(J))$, the smallest value of $J$ is just the critical point $J_{c}$. On the other hand, in the disordered phase, behavior of $\chi(J)$ depends on the relation between $L$ and $\xi$. For $L\gg\xi$ (i.e., $J$ is well below $J_{c}$), the exponential decay of $G(x)$ gives rise to $\chi\propto\xi^{2}$ and $\bar{\chi}(J,\eta^{*})\propto L^{\eta^{*}-2}\xi(J)^{2}$, which shows no scale-invariant points are possible (for $\eta^{*}\neq 2$). For other case ($L\lesssim\xi$), no explicit formula of $\chi(J)$ seems available. Although there is no analytic proof that $\bar{\chi}(J,\eta^{*})$ has no scale invariant points throughout the disordered phase, we confirmed that it is the case by numerical analysis applying the same method to obtain Fig. 3(d).

Below we shall see a close relation between the GHM and the 2D XY model, and explore the possibility of a BKT-type phase transition of the GHM. Let us first consider the limit of $c_{1}\to\infty$. Then $A_{\rm I}$ of Eq. (11) forces $U_{x,0}=1\ [\theta_{x,0}=0$(mod $2\pi$)] up to gauge rotation. By introducing the gauge-invariant angle variable $\sigma_{x}$ (we use the same letter as the above XY spin angle) on the spatial link $(x,x+1)$ [i.e., between the site $(x_{0},x_{1})$ and $(x_{0},x_{1}+1)$] as

$A_{\rm GH}$ reduces up to a constant to

This is just the action of the 2D XY spin model with the asymmetric couplings $c_{2}$ and $c_{3}$. It has the global U(1) symmetry under $\sigma_{x}\to\sigma_{x}+\alpha$. This limiting model exhibits BKT transitions on a certain critical line in the $c_{2}$-$c_{3}$ plane asymXY .

As $c_{1}$ is reduced from infinity, $\theta_{x,0}$ is allowed to fluctuate more. However, it is plausible to expect that the above critical line in the $c_{2}$-$c_{3}$ plane survive for finite but large enough $c_{1}$, and the BKT-type scaling relation (16) holds. Then, we determine the transition point $c_{1}=c_{1c}$ for a given set of $(c_{2},c_{3})$ using the scaled susceptibility density $\bar{\chi}$ of Eq. (16) calculated by $A_{\rm GHM}$ and Eq. (17). Explicit procedures are as follows;
(i) Fix $c_{2}$ and $c_{3}$ and measure the susceptibility $\chi(c_{1},L)$ along $c_{1}$ for three $L$’s ($L_{1},L_{2},L_{3}$).
(ii) Obtain interpolating curves $y(c_{1},L)$ for $\chi(c_{1},L)$ in third-order polynomials in $c_{1}$;
(iii) For $\bar{\chi}(c_{1},\eta,L)=L^{\eta-2}y(c_{1},L)$ (below we use $\eta$ instead of $\eta^{*}$ because no confusions may arise) , set the conditions of a triple intersection; $\bar{\chi}(c_{1},\eta,L_{1})=\bar{\chi}(c_{1},\eta,L_{2})=\bar{\chi}(c_{1},\eta,L_{3})$, and solve these two equalities numerically.
From Eq. (16), one expects a continuous (multiple) set $(c_{1},\eta)$ of solutions for every value of $c_{1}$ starting at and larger than the critical value $c_{1}=c_{1c}$.

In Fig. 3, we show (a) $C_{\rm GH}$ and (b) susceptibility density $\chi/L^{2}$ for $c_{2}=c_{3}=2.0$ and $L_{1}=16$, $L_{2}=32$, $L_{3}=64$. Instead of multiple set of solutions, we find a unique triple intersection at $c_{1}=c_{1c}\simeq 1.22$ for $\eta\simeq 0.68$. In Fig. 3 (c) we show three $\bar{\chi}$ with this value. To understand the absence of multiple solutions, we introduce the distance $\Delta\eta$ from the triple interaction defined as follows; Each pair of three curves $\bar{\chi}(c_{1},\eta,L)$ (e.g., $L_{a}$ and $L_{b}$) with fixed $c_{1}$ has an intersection at $\eta=\eta_{ab}$. We define $\Delta\eta\equiv$ Max of three distances, $|\eta_{ab}-\eta_{a^{\prime}b^{\prime}}|$. In Fig. 3 (d) we show $\Delta\eta$ together with the average of $\eta$, $\bar{\eta}\equiv(\eta_{12}+\eta_{23}+\eta_{31})/3$. The triple intersection at $c_{1}=c_{1c}$ locates at the unique minimum of $\Delta\eta\ [\Delta\eta(c_{1c})=0]$ and separates the large $\Delta\eta$-region ($c_{1}<c_{1c})$ and the small $\Delta\eta$-region ($c_{1c}<c_{1}$). By regarding $\Delta\eta$ in $c_{1c}<c_{1}$ almost vanishing, Eq. (16) tells us that the point $c_{1}=c_{1c}$ is an approximate critical point of BKT-type phase transition separating the disordered phase ($c_{1}<c_{1c}$) and the quasi-ordered phase $(c_{1c}\leq c_{1}$). Slow but steady development of $\Delta\eta$ from zero with increasing $c_{1}(>c_{1c})$ shows that there exist some corrections to the power-law decay of spin correlations. One may account for it by finite-size effect generally; in particular, by the logarithmic corrections janke that may become relatively important for smaller $\eta$. Then, we expect that $\Delta\eta$ for $c_{1}>c_{1c}$ tends to small if we use the date obtained for systems with larger $L$. To calculate the critical point for general $(c_{2},c_{3})$ below, we repeat to identify a unique triple intersection (scale-invariant point) as the BKT-type critical point. This manipulation is supported because a quite similar behavior of $\Delta\eta(J)$ is observed for the 2D XY model giving rise to $J_{c}\sim 1.11(1.12),\eta(J_{c})\sim 0.23(0.25)$ in good agreement with the known results in parentheses for larger $L$.

In Fig. 4, we show the phase diagram in the ($c_{1}c_{3}$-$c_{2}c_{3}$) plane obtained in the same procedure as in Fig. 3 (c). The four transition lines correspond to $c_{3}=0.5,1.0,2.0$ and $3.0$. The two coordinates $c_{1}c_{3}$ and $c_{2}c_{3}$ of this plane are chosen so that they are dimensionless and $\Delta\tau$ independent combinations (see Sec. III of Ref. Future for detailed discussion on this point). We note that there are no transitions in the region of $c_{1}c_{3}\lesssim 1.4$ because the data fitting in the form of Eq. (16) failed there. This is consistent with our assumption that finite but sufficiently large $c_{1}$ is a necessary condition for BKT-type transition.

Because the disordered phase necessarily realizes for small $c_{1},c_{2},c_{3}$’s, the lower (smaller $c_{2}c_{3}$) region in Fig. 4 corresponds to the disordered phase, and the higher (larger $c_{2}c_{3}$) region corresponds to the quasi-ordered phase. In gauge theory, typical and well-known phases are Higgs, confinement, and Coulomb phases. The magnitude $\Delta\theta$ of fluctuations of gauge field $\theta_{x,\mu}$ increases in this order wilson ; Kogut ; Fradkin . The magnitude $\Delta E$ of fluctuations of conjugate electric field $E_{x}$ decreases in this order due to the uncertainty principle. Because we have two phases in Fig. 4, we naturally identify the spin-disorder phase as the confinement phase and the quasi-ordered phase as the Higgs phase. This is supported by the phase diagrams obtained for related GHMs derived from the 2D and 3D EBHM U1GHM ; NJP1 ; Future .

In order to identify the relevant parameter regions in experiments, it is important to compare the phase diagram of Fig. 4 to that of the original EBHM Eq. (1). In Fig. 5, we show the phase diagram of the EBHM in small $V/J$ regime, which is obtained by the MC simulations in Appendix A. There are two phases. In the Mott insulator (MI) phase, the strong on-site repulsion $U$ stabilizes the uniform distribution of atoms without the phase coherence. For the large hopping $J$, the superfluid (SF) forms as a result of the Bose-Einstein condensation. In Fig. 5, the red dotted line shows the phase boundary of the EBHM between the MI and the SF, which was determined by a scaling method. For details, see Appendix A and Fig. A.1. To compare the phase diagrams of BHM and GHMs, we transport critical lines shown in Fig. 4 to those in the ($U/J\rho_{0}$-$V/J\rho_{0}$) plane of Fig. 5, by using the relations between the parameters such as $c_{1}c_{3}=2(V/(J\rho_{0}))^{-1}$ and $c_{2}c_{3}=2(J\rho_{0})/(U-2V)$. The phase boundaries of the two models are fairly in good agreement for $g^{2}>0$ and this result certifies the existence the confinement-Higgs phase transition.

In other words, our target phases are the MI phase and the SF phase in relatively small $V$ regime; the two phases in cold atomic system has explicit meaning of the lattice gauge theory. Fig. 5 also shows suitable parameter regions for the quantum-simulation experiment on the GHM. This fact certainly supports that the confinement-Higgs phase transition can be observed by experiments of ultra-cold atoms in a 1D optical lattice.

The phase diagram of the BHM in large $V$ regime itself has various interesting phases DMRG1 ; Kawaki . As the NN repulsion is getting larger compared to the on-site repulsion, the assumption of uniform density configuration $\hat{\rho}_{a}\to\rho_{0}$ is broken, and the density wave phase appears in the phase diagram, in which an imbalance in the density at the even and odd sites exists. The result of the MC simulation for large $V$ is also given in Appendix A.

In the high-energy physics, study of non-equilibrium dynamics of the gauge theories in Minkowski space-time has been a challenging problem. Its detailed study is going to be even more important inspired by recent heavy-ion collision experiments HIC . In this section, we study real-time dynamics of an electric flux in the Higgs and confinement phases observed by the path-integral MC simulations in Sec. III. To this end, we employ the GP equation, which is a reliable method to simulate dynamical behavior of a condensed fluid in condensed matter physics especially for the superfluid behavior GP2 . In general the GP equation is suitable for the SF (Higgs) phase, however in this section we extend the limitation, i.e., apply the GP equation to the MI (confinement) phase and focus on investigating qualitative behavior of both phases. A similar problem of real-time dynamics of an electric flux was addressed in an atomic quantum simulation of the Schwinger model, the (1+1)D quantum electrodynamics (QED), in Refs. Schw ; TNN , where the tensor network simulation was used. Interestingly enough, the calculations in Refs. Schw ; TNN report qualitatively different behavior of the electric flux depending on values of the gauge coupling and mass of the electron.

From Sec. III, the present $(1+1)$D GHMs have two different phases, the Higgs and the confinement phases, therefore study of dynamical behavior of the electric flux gives important and interesting insight into the gauge dynamics in each phase. A similar problem for the higher-dimensional cases has been already addressed in our previous work NJP1 ; Future , where we have obtained qualitative difference of flux string behavior for each phases. In the present $(1+1)$D case, the direct comparison between the GHM in Eq. (5) and the Schwinger model is possible for the string dynamics.

The GP equation can be derived from the coherent path-integral of the EBHM as a saddle point equation. When the Lagrangian $L$ is described as $L=-{{\sum}_{a}}i\psi_{a}^{*}(d{{\psi}_{a}}/dt)-H_{\rm EBH}$, the Euler-Lagrangian equation $\delta L/\delta\psi^{*}_{a}=0$ becomes the GP equation. From this procedure, the GP equation for the EBHM is given by

where we have added chemical potential term $\mu^{\prime}\psi_{a}$ by which the average density $\langle|\psi_{0}|^{2}\rangle=(1/N)\sum_{a}|\psi_{a}|^{2}=\rho_{0}$ is controlled. In what follows, we set the equilibrium density $\rho_{0}=1$ by putting $\mu^{\prime}=\rho_{0}(U+2V)-2J$. Here it should be remarked that this is a kind of the normalization and results for other densities can be obtained by a simple rescaling as $\psi_{a}\rightarrow\sqrt{\rho_{0}}\psi_{a}$ with $(U,V)\rightarrow{1\over\rho_{0}}(U,V)$. A system of $200$ spatial lattice sites with the free boundary condition is used and the time step $\Delta t$ for the numerical calculation is set as $\Delta t=10^{-4}$. The time scale $t=1$ is $\sim$ 0.32 [msec] by using the energy scale $U/h\sim$ 500[Hz] in typical experiments. In Eq. (19), the dissipation term is not included, and therefore the total energy is conserved during the time evolution.

As an example of a non-equilibrium state of the system, we consider the state with a single electric flux with a length $R$, i.e., we set the electric flux string in the initial state and solve the GP equation numerically. Hereafter we use the gauge lattice label, $a\rightarrow(r,1)$. Explicitly, we set the initial configuration as $|\psi_{r,1}|^{2}=[1+(-1)^{r}(0.1)]\rho_{0}$ for $r_{1}-R\leq r\leq r_{1}+R-1$ ($r_{1}$ is the center of the string). On the other sites, $|\psi_{r,1}|^{2}=\rho_{0}$. In experiments the above initial flux is created by the density modulation Wurtz . In our numerical simulation, we set the length of the electric flux $R=60$. In the confinement phase, this length is expected to be long enough for observing the typical quark-confinement phenomena Kogut ; Smit .

We first briefly consider the case of sufficiently small $J$, and verify that the electric flux string is stable by solving the GP equation (19). In our numerical simulation, the electric field $E_{r}$ are defined as

This corresponds to a density fluctuation from mean density. The typical numerical results is shown in Fig. 6. When $J$ is very small, the effects of the Higgs field is very small and the system has properties of the pure compact U(1) gauge theory.

In Fig. 7 (b)-(d), we show the time development of an electric flux along line ${g^{2}}/J=0.2$ in the phase diagram. The results show that from the Higgs phase to the confinement phase, the stability of the electric flux is increased. In fact, there is a striking difference between the Higgs phase and the confinement phase; in the Higgs phase the electric flux spreads out. On the other hand in the confinement phase, the form of the initial flux does not change. In addition, we find that in both phases an oscillation of the electric field in the flux takes place. This behavior of the flux string means that a pair production of the Higgs particle and the resultant string breaking takes place. But in the confinement phase the restoration of the flux string readily occurs as a result of the strong gauge coupling. By contrast, in the Higgs phase the initial electric flux is unstable, i.e., bits of electric flux appear and they spread out from the initial region $r_{1}-R\leq r\leq r_{1}+R-1$. This behavior comes from the condensation of the Higgs field and the resultant spontaneous symmetry breaking of the global U(1) gauge symmetry, i.e., the charge is not conserved in the Higgs phase.

In order to evaluate the spread of the electric field from the initial string, we measure the mean value of electric field strength $\langle E^{2}\rangle_{\rm out}$ outside of the initial place of the electric flux,

Fig. 8 shows the $\langle E^{2}\rangle_{\rm out}$ for some typical values of the Gauss-law coupling $V$ along the line ${g^{2}}/J=0.2$. It is obvious that as the value of $V$ decreases, $\langle E^{2}\rangle_{\rm out}$ increases rapidly as a function of time, i.e., the initial configuration of the electric flux is unstable, and electric-field excitations propagate out of the initial electric flux.

In addition, Fig. 7 exhibits another remarkable behavior of the electric field as non-equilibrium phenomena. The results show that on the sites, where the initial electric flux exists, the electric fields are oscillating and changing their sign in both the Higgs and the confinement phases. The data in Fig. 7 (e) and (f) shows the time evolution of the electric field averaged over the central ten sites $\langle E\rangle_{\rm center}$. Similar behavior was observed in all three cases; the central electric field in all cases oscillates in time and takes even negative values. The result is similar to behavior of an electric flux in Schwinger model with a small electron mass as reported in Refs. Schw ; TNN , and the phenomenon is called the string-antistring oscillation. Quantum link version of QED was also studied from this point of view Banerjee . All the above studies indicate that the pair creation takes place in the central region of the electric flux as the Schwinger mechanism tells Schwinger . The two initial static charges residing on the edges of the initial electric flux tend to be shielded by particles of opposite charge. The initial electric flux vanishes as a result. After that, the reverse process happens, i.e., the pair creation recurs and the electric flux in the opposite direction is created. Iterating this process causes consequently the string-antistring oscillation. Our results in Figs.7 (b)-(f) for the U(1) GHM suggest that a similar string-antistring oscillation process takes place as in QED. Fig. 9 illustrates schematically the string-antistring oscillation process.

In the above, we studied the real-time dynamics of the electric flux string by using the GP equation. Results concerning to the Higgs phase are reliable as the Higgs phase is nothing but the SF and the order parameter $\psi_{a}$ has only small fluctuations there. On the other hand, one may think that use of the GP equation for the MI regime is not appropriate even in the region close to the phase boundary to the SF. However, we expect that our study for the weak MI regime captures at least qualitatively correct picture of the string dynamics. In fact, the numerical results presented above are very similar to those obtained by other numerical results for closely related gauge models in Refs. Schw ; TNN .

Anyway, a more detailed study of the dynamics of the flux string in the confinement-MI regime is important and it gives some useful remarks on the experimental setup. To this end, the idea of the truncated Wigner approximation (TWA) is useful TWA , that is, we simulate the real-time dynamics of the flux string with different initial conditions. In the above, on solving the GP equation, we employed the initial condition such that the phase of the wave function $\theta_{r,1}=0$. However in the MI state, the phase of the ‘condensate’ fluctuates as a result of the uncertainty relation between the particle number and the phase. In other words, the fluctuations of the phase cannot be avoided in the MI even if a prominent experimental technique to control the phase of the ‘condensate’ is used. This fact arises questions about the reliability of the results obtained by the GP equation.

Here, to study the effect of the phase fluctuation, we employ the initial condition such that $\theta_{r,1}=\lambda_{r}\pi$ and solve the GP equation, where the variables $\{\lambda_{r}\}$ are homogeneously distributed random numbers between $[-\Delta_{m},\Delta_{m}]$ with $\Delta_{m}>0$. To calculate the physical quantities like the density at each site, we take 100 samples, which have different initial conditions produced by the random numbers $\{\lambda_{r}\}$, and take average of the calculated quantities over the samples. In this way, we obtain results of the real-time dynamics of the electric flux. This sampling method essentially corresponds to the TWA TWA .

Fig. 10 shows the real-time dynamics of the flux string for the parameters $V/J=1$, $U/J=2.2$, which locates in the weak MI region in the phase diagram, and for two cases of the random distribution of the phase, $\Delta_{m}=0.01$ and $0.05$. The other conditions are the same with those in the previous calculation. For $\Delta_{m}=0.01$, the stable evolution of the string is observed, whereas for $\Delta_{m}=0.05$, only small portion of the initial string configuration remains. From the above result and the uncertainty relation of the number and phase, $\Delta n\cdot\Delta\theta\geq{1\over 2}$, the fluctuation of the atomic number $\Delta n$ at the initial state is required as $\Delta n\gtrsim 1/(0.05\cdot 2\pi)\sim 3.2$. To satisfy this requirement in experiments, the average particle number must be sufficiently large and also a sufficient number modulation in the flux string is needed. If this condition is satisfied, the string-antistring oscillation in Fig. 7 is to be observed. More detailed study on this problem by using the full TWA and its extension TWA will be reported in a future publication.

Finally let us comment on the real-time dynamics of the flux string in the confinement-MI region apart from the phase boundary to the SF. In this region, the GP equation is not applicable, as the phase of the ‘condensate’ fluctuates rather strongly. However, this region corresponds to the strong-coupling region of the GHM, and therefore the string tension of the electric flux is strong. Then it is naturally expected that the string is quite stable and its fluctuations in time are small, and as a result its electric field is kept positive, $E_{r}>0$. In fact, this phenomenon was observed in the numerical study of the gauge models in Refs. Schw ; TNN .

We explained how to realize the atomic quantum simulation of the (1+1)D GHM by using cold atomic gases on the optical lattice. By means of the path-integral MC simulations, we obtained the phase diagrams of the GHM and also the EBHM at large filling and verified that these two phase diagrams are in good agreement with each other for the parameter region of the positive gauge coupling. We examined the phases of the EBHM from the gauge-theoretical point of view and the result is summarized in Table 1.

By using the GP equations for the EBHM, we investigated dynamical properties of the electric flux put on a line as an initial state. In the confinement region of the GHM, the string breaking and restoration of the electric flux is observed, which has similar properties to the Schwinger mechanism in QED. On the other hand in the Higgs region, the flux spreads out and many bits of electric flux develop. Finally we discussed the reliability of the results obtained by the GP equations by using the idea of the TWA, and gave some remarks for the experimental setup to observe the string-antistring oscillation. We hope that the above phenomena will be observed experimentally in the near future. In particular, dynamical behavior of the transition from the confinement to Higgs phases is quite interesting and expected to be observed by the strong controllability and the versatility of the cold atomic system.

Finally, let us note that our approach identifies the phase and the amplitude of atomic variable as U(1) vector potential and the electric field respectively, and introduce the Higgs field to relax the crucial fine tuning to realize the Gauss’ law. These general and universal procedures make us possible to define corresponding GHMs of wide variety both for high and low fillings and also for homogeneous and inhomogeneous atomic distributions. Although this GHM may contain nonlocal interactions in some cases and look quite different from the conventional lattice gauge theory, it is important that GHM is of course a “gauge theory having gauge symmetry”. It obviously opens the way to study gauge theory in more general view points. Also such gauge-theoretical view points may shed some lights backward on studying the original atomic systems as partly shown in Sec. A.3.