Dynamics of first-order quantum phase transitions in extended Bose-Hubbard model: From density wave to superfluid and vice-versa

In this paper, we study the nonequilibrium dynamics of the Bose-Hubbard model with the nearest-neighbor repulsion by using time-dependent Gutzwiller (GW) methods. In particular, we vary the hopping parameters in the Hamiltonian as a function of time, and investigate the dynamics of the system from the density wave (DW) to the superfluid (SF) crossing a first-order phase transition and vice-versa. From the DW to SF, we find scaling laws for the correlation length and vortex density with respect to the quench time. This is a reminiscence of the Kibble-Zurek scaling for continuous phase transitions and contradicts the common expectation. We give a possible explanation for this observation. On the other hand from the SF to DW, the system evolution depends on the initial SF state. When the initial state is the ground-state obtained by the static GW methods, a coexisting state of the SF and DW domains forms after passing through the critical point. Coherence of the SF order parameter is lost as the system evolves. This is a phenomenon similar to the glass transition in classical systems. When the state starts from the SF with small local phase fluctuations, the system obtains a large-size DW-domain structure with thin domain walls.


Introduction
In recent years, dynamics of quantum-many-body systems is one of the most actively studied subjects in physics. Process in which a system approach to an equilibrium is of fundamental interests, and also evolution of system under a quench has attracted many physicists. Nowadays, ultra-cold atomic gas systems play a very important role for the study on these subjects because of their versatility, controllability and observability [1]. Theoretical ideas proposed to understand transient phenomena are to be tested by experiments on ultra-cold atomic systems. This is one of examples of so-called quantum simulations [2][3][4][5].
For the second-order thermal phase transition, time evolution of systems under a change in temperature has been studied extensively so far. From the viewpoint of cosmology, Kibble [6,7] claimed that the phase transitions lead to disparate local choices of the broken symmetry state and as a result, topological defects called cosmic strings are generated. Later, Zurek [8][9][10] pointed out that a similar phenomenon is realized in laboratory experiments on the condensed matter systems like the superfluid (SF) of 4 He. After the above seminal works, many theoretical and experimental studies on the Kibble-Zurek (KZ) mechanism have appeared [11].
Concerning to experiments on Bose-condensed ultra-cold atomic gases, the correlation length of the SF and the rate of topological defect formation were measured and the KZ scaling hypothesis was examined [12,13].

Extended Bose-Hubbard model and slow quench
We consider the EBHM whose Hamiltonian is given by [39] H J a a U n n V n n n h.c. = † , and μ is the chemical potential. J(>0) and U(>0) are the hopping amplitude and the on-site repulsion, respectively. We also add the NN repulsion with the coefficient V, which plays an important role in the present work.
In this study, we are interested in cases near the half filling, i.e., N n 1 1 2 s i i r º å á ñ » , where N s is the total number of the lattice sites, and we take N s =64×64 or 100×100 for the practical calculation. We set U=1 as the energy unit, and time t is measured in the unit ÿ/U. We investigated the system in equation (1) by using the static GW approximation and show obtained ground-state phase diagram in figure 1 for V/U=0.05. There exist three phases, i.e., the DW, SF and supersolid (SS) although the area of the SS in the phase diagram is small for V/U=0.05. We also show the system energy, particle density and amplitude of the SF order parameter, figure 2 for μ/U=0.1. From the results in figure 2, it is obvious that the system exists near the half filling ρ≈1/2, and a first-order phase transition between the DW and SF takes place at J c /U;0.022 as a finite jump in Y | |indicates. The existence of the first-order phase transition is quite plausible as the DW and SF have both the own long-range order. In recent paper [40], we studied the EBHM for V/U=0.375 and near the unit filling ρ≈1. There exists a substantially finite region of the SS in addition to the DW and SF. These three phases are separated by two second-order phase transitions. This result is in agreement with the quantum Monte-Carlo (MC) study [41].
In the following, we shall study dynamics of the system under 'slow quenchs'. To this end, we employ the tGW methods [42][43][44][45][46][47][48]. In the tGW approximation, the Hamiltonian of the EBHM in equation (1) is approximated by a single-site Hamiltonian H i , which is derived by introducing the expectation value where iNN denotes the NN sites of site i, and Hartree-Fock type approximation has been used for the hopping and NN repulsion. To solve the quantum system H GW in equation (2), we introduce GW wave function f t n n n n n , , 3 where n c is the maximum number of particle at each site, and we mostly take n c =6 in the present work. Some quantities are calculated with n c =10 to verify that n c =6 is large enough for the study of the half filling case. See figures 3 and 7.
In terms of f t n i { ( )}, the order parameter of the SF is given as where Q t is the quench time, which is a controllable parameter in experiments. We employed 10 samples as the initial state at t=−τ Q (i.e., J 0 , which have the DW order with small local density fluctuations from the perfect DW. Then, we solve equation (5) to obtain GW F ñ | . Physical quantities for which scaling lows are examined are obtained by averaging over samples. The linear quench in equation (6) is terminated at t=t f with J t J 0.044 c f = > ( ) ( ) in the numerical study. Subsequent behavior of the system is also observed to see how the system approaches to an equilibrium.
We show the typical behavior of Y | |as a function of t in figure 3 for τ Q =300. At t=0, the system crosses the critical point at J c /U;0.022. After crossing the critical point, Y | |remains vanishingly small for some period, and then it develops very rapidly. After the rapid increase, Y | |starts to fluctuate and coarsening of the and t eq ≈120, respectively. On the other hand, t ex ≈400, at which the oscillation of Y | | terminates. From t eq to t ex , coarsening process of the phase of Ψ i takes place in large scales [26]. (Lower panel) Calculation of Y | | in the n c =10 case is also shown. It is in good agreement with that of n c =6.
phase of the SF order parameter takes place there [26]. t in figure 3 is defined as t | ( )|, and t eq is the time at which the oscillation of Y | |starts. Similarly, t ex is the time at which that oscillation terminates. Similar behavior to the above was observed in the Mott to SF quench dynamics and examined carefully [26]. Compared with the Mott to SF dynamics, the SF amplitude Y | |is smaller, e.g., for t t eq > , 0.8 0.9 Y| | ( -)in the Mott to SF transition, whereas 0.5 Y| | in the present case. This difference simply comes from the difference of the mean particle density, i.e., ρ∼1 in the Mott to SF transition case. It is interesting to study the correlation length ξ of the SF order parameter and the vortex density N v as a function of the quench time τ Q . These quantities are defined as follows;

The DW order parameters
is the unit vector in the x (y) direction. For continuous second-order phase transitions, the KZ hypothesis predicts a scaling law such as Recently, applicability of the above KZ scaling law for second-order QPT has been discussed for several quantum systems. On the other hand for first-order phase transitions, it is commonly expected that such a scaling law does not hold as the relaxation-time cannot be defined properly. For a classical statistical model, another type of scaling law was proposed for first-order phase transitions [36]. It should be also noted that off-equilibrium dynamics of a quantum Ising ring was investigated recently and finite-size scaling laws for first-order phase transitions were proposed [49]. There, off-equilibrium scaling variables were given in terms of an energy gap and quench time, and physical quantities were obtained as a function of time.
To see if scaling law exists or not, we measured ξ and N v at t t =ˆand t t eq = . In the original KZ hypothesis for continuous phase transitions [11], t is the time at which the system re-enters an equilibrium after the freezing (or impulse) period. On the other hand, t eq is the time at which a coarsening process of the SF phase coherence starts [26]. . On the other hand at t t eq = , data at each Q t exhibits slightly large fluctuations but scaling laws for the correlation length, N v and t eq seem to Figure 5. Scaling laws observed for the correlation length ξ, vortex number N v at t t =ˆ, and t with respect to τ Q . Figure 6. Scaling lows observed for the correlation length, vortex number at t t eq = , and t eq with respect to Q t . exist for 20 Q t > . The above results indicate that besides the KZ mechanism, there exists another mechanism to generate the scaling laws. Possible explanation is given in section 4. It should be noted that after passing the critical point, Δ DW and Δ SF have even-odd site fluctuations, and therefore, the system is not homogeneous. We think that because of this inhomogeneity, the critical exponents of ξ and N v at t t =ˆdo not satisfy the expected relation such as b=d/2. On the other hand at t t eq = , the system is rather homogeneous, and therefore b∼d/2.
In the appendix, we consider the hard-core version of the EBHM and show the calculations of the scaling laws with respect to Q t in figure A2. There, t x (ˆ) and N t v (ˆ) fluctuate rather strongly. This behavior comes from the fact that fluctuations of the particle number at each site is smaller compared with the soft-core case, and as a result, the stability of the phase degrees of freedom of the SF order parameter is weakened.
We terminate the linear quench at t 300 . After t f , the system approaches to an equilibrium as the results in figures 3 and 4 indicate. It is interesting to see how the correlation length of the SF develops. As the results in figure 7 show, the correlation length increases after passing the critical point as it is expected. However, its increase gets weak at t t eq , and it saturates at t∼500 and keeps a finite value. To study the resultant phase, we measured N v and found that there exist no vortices at t>500. One may expect that the system settles in a finite-temperature (T) SF phase for sufficiently large t with an effective T,T eff . The finite-T SF in 2D has a quasi-long-range order and the correlation length diverges, i.e., the Kosterlitz-Thouless (KT) phase. The above result seems to indicate that some other state is realized in the final stage of the present process. However, the system behavior may strongly depend on the average particle density ρ. Further study is needed to clarify this interesting problem. In fact, we studied this problem in the case of the mean particle density ρ≈1 and V/U=0.375 [40]. In the quench process such as the DW→SS→SF, the correlation length continues to increase even for large t. This result seems to indicate that a KT phase of the SF is realized there.

Consideration by the GL theory
In the previous section, we showed that the results obtained by the GW methods indicate the scaling laws of t, t eq and the correlation length with respect to the quench time Q t . It is interesting and also important to study the origin of these observations from more universal and intuitive point of view. To this end, the GL theory is quite useful. In fact very recently, it was pointed out that the GL theory can drive the scaling laws for the second-order phase transition by analytical transformation of the associated equations of motion [50]. In this section, we first review the above derivation of the scaling laws for the ordinary second-order phase transition, and then give an intuitive picture of the scaling laws by using a classical solution representing decay of the false vacuum. Then, we extend the methods to the present case involving the SF and DW order parameters. This consideration also gives an insight about the physical meaning and limitation of the GW methods.

Second-order phase transition
Let us start with the stochastic GL equation for a complex order parameter (condensate T is the temperature of particles ensemble not participating the Bose-Einstein condensate. As in [50], we consider the critical parameter ò(t) such as where λ is a parameter for the quench protocol. Then, let us change variables as follows In equation (11), the Q t -dependence in equation (9) disappears except the last white-noise term. From the above fact, it is concluded in [50] that the Q t -dependence of t and t x (ˆ) are expected to follow the transformation in equation (10), and they are given as follows for sufficiently low T t t , .
. The above estimations agree with those of the KZ scaling with the mean-field exponents such as 1 2 n = and z=2. As we show, the above scaling transformation gives an intuitive picture that derives the KZ scaling law. To this end, we put r t , 0 Q  = ( ) in equation (8) and consider a static potential such as ò(t)=−ò 0 <0. In this case, the static ground state is given as To study the sudden quench dynamics, we consider the decay of the false vacuum f=0 to the true ground state In 1D case, a classical solution representing the decay is obtained as follows [38] t x x v t , 1 e x p 2 , . The solution equation (13) obviously represents the situation in which the true vacuum 0  f = born in the false vacuum expands with the speed v 0 . Let us consider the 'slow' quench dynamics and study bubble-nucleation-evolution process in the SF formation. We expect that this process corresponds to the numerical studies in the previous sections. We have to find the solution to equation (8) that describes a single SF-bubble evolution in the false vacuum f=0, but we cannot find an exact solution. However, the above solution in equation (13) suggests that a spherically symmetric solution in higher dimensions and also for the time-dependent ò(t) has the following form for ò(t)<0 3 where v C t t 0  = | ( )| with a certain constant C 0 , and F(x) is a decreasing function such as F 1 -¥ = ( ) and F 0 ¥ = ( ) . In fact, we can show that the function r t , s f  ( )in equation (14) satisfies the scaling transformation in equation (10) for the time-dependent ò(t) in equation (9), i.e.
does not depend on Q t . As far as the above picture holds in the time evolution of the system, equation (15) implies that typical events and phenomena are observed similarly in systems with various Q t ʼs, and corresponding times have Q t -dependence such as Q 1 t l l+ ( ) . For example, we numerically obtained t and t eq for various Q t ʼs in section 3 by starting with qualitatively the same initial states. These values are related to Q t -independent ĥ and eq h that are obtained by the rescaled picture from equation (15), i.e., t and t eq in the Q t -system are given by t . Furthermore, a typical linear size of the bubble at t, i.e., the correlation length at t, ξ(t), is given as and therefore, t After t eq , the merging and coarsening process of SF bubbles takes place [26], and therefore the above picture and also the resultant scaling laws do not hold anymore.

GL theory, GW methods and quantum MC simulation
Here, it is suitable to comment on the GW approximation. The GL theory and also the Gross-Pitaevskii (GP) equation consider only the mean field and totally ignore fluctuations around it. On the other hand in the GW approximation, we focus on a wave function of site factorization, and wave function at each site is obtained by solving the site-factorized Hamiltonian in which the NN operators are replaced with their expectation values [26]. The uncertainty relation between the particle number and phase at each site is faithfully taken into account although an equation of motion similar to the GL (GP) equation is derived by the GW methods. This is an advantage of the GW approximation over the GL and GP theories.
As more reliable methods, let us consider the quantum MC simulations of the coherent-state path-integral in the imaginary-time formalism. In this MC simulations, quantum operators are reduced into classical variables and the quantum superpositions are treated by the fluctuations in the imaginary-time direction. Large number of configurations are generated by the MC updates and physical quantities are calculated by averaging them over generated configurations. In the Metropolis MC algorithm, the local updates are applied to variables at each site by calculation a local energy around that site. In the vicinity of a phase transition point, a large number of configurations contribute equally, and calculations by large CPU times are required in order to take into account all relevant configurations. On the other hand away from the critical point, the number of important configurations is not so large. From the viewpoint of the MC simulation with the local update, we can get an interesting insight into the GW approximation. That is, let us imagine that we perform a GW calculation for a system with size 10 4 ×10 4 . When we calculate expectation values, we divide the 10 4 ×10 4 system into 10 4 number of 10 2 ×10 2 subsystems. We obtain the expectation values by averaging values calculated in each subsystem. Compared with the path-integral MC simulation, this method is more reliable as the uncertainty relation is faithfully respected. (In the path-integral MC simulation, this relates to the problem how accurately effects of the Berry phase are taken into account. See for example, [51].) However in the vicinity of the phase transition, 10 4 configurations are not sufficient to obtain physical quantities closely related to the singularities of the phase transition. The above consideration suggests that the GW methods are a fairly good approximation for calculating physical quantities that are finite even for the critical regime, e.g., finite order parameters. In other words, the estimation of the critical exponents by the GW methods is not reliable even for using very large systems.
The above consideration may over estimate the reliability and applicability of the GW methods, but it explains why the GW methods often succeed in obtaining correct results such as the phase diagrams, etc. We 3 Solution in equation (14) might be regarded a solution in the slow quench limit, in which the time-derivative of ò(t) is small. However, it also satisfies the scaling transformation with ò(t) in equation (9). See the discussion below. 4 Rough estimation of ĥ and η eq are the followings. As t is determined by the condition such as t constant for the 2D case. On the other hand, as t eq is the time at which the overlap of SF bubbles starts [26], v eq eq h h = ( ) constant. Simulation for various λ's is a future work. expect that the GW methods also works for the correlation functions as far as the correlation length is finite as the quantum MC simulations do, although at present there are no ways to verify it in the quench dynamics.

First-order phase transition in vicinity of triple point
As the phase diagram in figure 1 shows, the present first-order phase transition is located in the vicinity of the triple point of the DW, SF and SS. The GL theory for the quench dynamics in section 4.1 can be applied to this case with some modification. Besides the SF order parameter, we introduce a coarse-grained real DW order parameter, D r t n , . GL equations are given as where the positive parameters g 1 , g 2 and g 3 are phenomenological ones, which are to be determined by the parameters U and V. The positivity of g 3 comes from the fact that the SF and DW are competing orders in the original EBHM. On the other hand, ò(t) and m(t) are parameters that are determined by J(t), U and V. In the quench from the DW to SF, both ò(t) and m(t) are decreasing functions of t. Let us consider a slow quench, and denote the phase transition time from the DW to SF by t c . At ), the system is in the DW and then, t g D t t g has a phase coherence in the SF, which induces amplitude fluctuations, as the amplitude and phase of the SF order parameter are quantum conjugate variables with each other. On the other hand in the would-be DW region for t>0, the phase coherence is lost, and then the SF amplitude is stable. The DW order parameter Δ DW does not have a stable finite value even after passing through the critical point at t=0. These results indicate that some kind of domain structure forms there, i.e., small DW domains may coexist with local SF regions.
Calculations of the amplitude of Ψ i and the particle density at t Q t = are shown in figure 9. As expected above, DW domains and regions with finite SF amplitude coexist without overlapping with each other.
In Case A, the quench stops with J 0 Q t = ( ) , and therefore no movement of particles occurs after the quench, and the particle density snapshot in figure 9 continues to describe the states for t Q t > . Similarly, we expect that the coherence of the phase of Ψ i is destroyed at t Q t = because J 0 Q t = ( ) and also 300 Q t = is a slow quench. See figure 9. In order to verify the expected behavior of Ψ i , we measured the vortex density as a function of time. At t Q t = , N v ∼300 is sufficiently large. In summary, in Case A with 300 Q t = , an inhomogeneous state with local DW and SF domains forms after quench. SF order parameter gradually loses its phase coherence during the slow quench.
On the other hand for cases of smaller 100 Q t = and 50, the SF order parameter Ψ i is finite even at t Q t = , and it varies after t Q t = . The phase of Ψ i gradually loses its long-range coherence by the existence of the repulsive interactions for t Q t > .  (see figure 10). We also study how the system evolves after t f . Observed quantities are shown in figure 10 for 50 Q t = . The DW order parameter Δ DW develops but its value fluctuates in rather long period after passing J c as in Case A. The total energy slightly decreases until t f , and the kinetic and on-site energies exhibit fluctuating behavior for t<t f although the NN interaction energy is rather stable. This behavior mostly originates from the local density fluctuations, and the stability of the NN interaction comes from the cancellation mechanism between NN sites j iNN Î . After passing the critical point at t=0, the Ψ i keeps a coherent SF order for some period as the calculation of the vortex number N v indicates. At t≈100, it starts to lose the coherence and the SF is destroyed as the increase in N v indicates. The state at t∼t f is a supercooled state, and a coexisting phase of local domains of the DW and SF is realized there. The observed phenomenon after t>t f , therefore, has very similar nature to the glass transition, in which the phase coherence and superfluidity are getting lost as the supercooled state evolves after the quench. We call it quantum glass transition (QGT) as the hopping amplitude J, instead of temperature, is the controlled physical quantity and the relevant transition is quantum mechanical one instead of thermal one. We have verified that similar phenomenon is observed for other values of Q t , e.g., 20 Q t = and 200.
In both Case A and Case B, the above mentioned QGT is observed dynamically as a nonequilibrium phenomenon, i.e., the QGT point is passed through as the system evolves. Therefore as the next problem, it is interesting to see whether there exits a genuine glass transition point, J J c g < ( ). Below J g , the supercooled state is meta-stable or at least has a long life time, and the SF survives without losing its phase coherence. For Cases A and B, J<J g . Then as Case C, we studied the quench whose finial point is J(t f )=0.02, i.e., very close to the equilibrium critical point. Obtained order parameter Y | |and vortex number N v are shown in figure 11 for 50 Q t = , and time evolution of the particle density, amplitude and phase of Ψ i are shown in figure 12. After passing the critical point J=J c at t=0, the domain formation of the DW starts as shown by the particle density snapshot in figure 12, whereas the long-range coherence of the SF order parameter Ψ i exists there. Compared with the cases of J(t f )=0 and J(t f )=0.01, the destruction of SF and formation of the DW region are slow, but , the total energy of the system keeps a constant value as the system is an isolated one. Figure 11. Transition from SF to DW with J(t f )=0.02, Case C. Increase of N v is slow compared to the cases J(t f )=0 and J(t f )=0.01. SF amplitude Y | | also keeps a finite value even for t→large. However, N v increases smoothly, and therefore, the supercooled state formed in the quench is not a meta-stable state. after t>450, the quantum glass state forms. Local DW domains develop but also empty regions (voids) form. SF order loses a long-range coherence. This result indicates that J g cannot be observed. Similar results are obtained for the case of 20 Q t = and 200 Q t = .

Evolution from SF state with small phase fluctuations
In section 5.1, we studied dynamical evolution of the system from the SF to DW. In that study, the initial state is set to the ground state obtained by the equilibrium GW methods. It is interesting to see how the dynamical phenomena depend on the initial state as we are considering the first-order phase transition. In order to study this problem, we consider a SF state that is uniform and has almost perfect phase coherence with very small random fluctuations. For the practical calculation, we employ an initial state GW wave function in equation (3) corresponding to e The other condition is the same with the Case A, (please refer to the left panel in figure 8). We call the present study Case D. We investigated the time evolution of the system by the tGW methods, and obtained results are shown in figure 13. Interestingly enough, the system behavior after passing across the critical point J c is substantially different from that in Cases A. The SF order parameter Y | |decreases a finite amount at t∼100, and the density difference at even-odd sublattice increases there. On the other hand, the vortex number starts to increase rapidly at t∼150.
Snapshots of the particle density, SF amplitude and SF phase are shown in figure 13. Contrary to Case A, the DW pattern starts to form at t∼115 and it develops to the whole system at t∼300, even though there exist domain walls. It should be noticed that a similar behavior was observed for the classical first-order phase transition in [36]. On the other hand, the SF phase coherence exists at t<115, whereas it is destroyed at t∼300.
The initial state of Case D has higher energy than that of Case A. The above numerical result indicates that there exists an energy barrier between the supercooled SF state and the genuine DW, and some amount of energy is need to overcome the barrier. Furthermore, the above result also indicates that the existence of the SF phase coherence in large spatial regions prevents the formation of large size DW domains. In other words, local fluctuations of the superfluidity coherence substantially develops under a quench even if they are initially tiny, and the DW is preferred as a result.
We expect that the above interesting phenomenon is observed by experiments on ultra-cold atomic gases in the near future. In this appendix, we shall give numerical calculations of the physical quantities concerning to the static properties of the model in figure A1, and also the Q t -dependence of t, etc in the quench dynamics in figure A2.
The results in figure A1 obviously show that there is a first-order phase transition from the DW to SF for increasing J/V as the quantum MC simulations in [53,54] proved. On the other hand, the correlation length ξ and vortex density N v fluctuate rather strongly compared to the soft-core cases. This result comes from the fact that the HCEBHM has a small fluctuations in the particle number at each site, and as a result, the phase of the SF order parameter fluctuates rather randomly.  The correlation length ξ and vortex density N v fluctuate rather strongly compared to the soft-core cases. This result comes from the fact that particle-number fluctuation at each site is restricted by the hard-core constraint, and as a result, fluctuation in the phase of the SF order parameter Ψ i is getting large.