Quantum XXmodel with competing short and longrange interactions:
Phases and phase transitions in and out of equilibrium
Abstract
We consider the quantum XXmodel in the presence of competing nearestneighbour and globalrange interactions, which is equivalent to a BoseHubbard model with cavity mediated global range interactions in the hard core boson limit. Using fermionic techniques the problem is solved exactly in one dimension in the thermodynamic limit. The ground state phase diagram consists of two ordered phases: ferromagnetic (F) and antiferromagnetic (AF), as well as an XYphase having quasilongrange order. We have also studied quantum relaxation after sudden quenches. Quenching from the AF phase to the XY region remanent AF order is observed below a dynamical transition line. In the opposite quench, from the XY region to the AFphase beyond a static metastability line AF order arises on top of remanent XY quasilongrange order, which corresponds to dynamically generated supersolid state in the equivalent BoseHubbard model with hardcore bosons.
pacs:
I Introduction
Recently there is an increased interest to study the phase diagram and nonequilibrium dynamics of quantum manybody systems with competing short and longrange interactions. Experimentally, such systems have been realised with ultracold atoms in optical lattices inside a highfinesse optical cavityScienceEsslinger ; Klinder2015 ; Landig2016 . The strength of the shortrange (onsite) interaction is related to the depth of the optical lattice, while long (infinite) range interactions are controlled by a vacuum mode of the cavity. The interplay of short and longrange interactions may result in a rich phasediagram with exotic phases and interesting nonequilibrium dynamics.
Classical manybody systems with competing short and ferromagnetic globalrange interactions have been studied earlier and in the thermodynamic limit the orderparameter is obtained through a selfconsistent treatment, like in meanfield models. Kardar1983 ; Mukamel2005 . Theoretical results for the phase diagrams of quantum many body systems with competing shortrange and globalrange interactions are rare and up to now confined to the aforementioned bosonic systemChen2016 ; Dogra2016 ; Niederle2016 ; Sundar2016 ; Panas2017 ; Flottat2017 .
The nonequilibrium dynamics of closed quantum many body systems after sudden quenches has attracted a lot of attention in the last decade. Here one is interested in the timeevolution of different observables, such as the order parameter or some correlation function, after the quench. Fundamental questions concerning quantum quenches include i) the functional form of the relaxation process in early times, and ii) the properties of the stationary state of the system after sufficiently long time. The latter problem is related to the question of thermalization, which is expected to be different for integrable and nonintegrable quantum systems. Nonintegrable systems are expected to evolve into a thermalized state, in which local observables can be characterized by thermal expectation values Rigol_07 ; Calabrese_06 ; Calabrese_07 ; Cazalilla_06 ; Manmana_07 ; Cramer_08 ; Barthel_08 ; Kollar_08 ; Sotiriadis_09 ; Roux_09 ; Sotiriadis_11 , but there are some counterexamples larson ; hamazaki ; blass_rieger . On the other hand, integrable systems in the stationary state are generally described by a socalled generalized Gibbs ensemble, which includes all the conserved quantities of the system. Recently it has been observed that in certain models both local and quasilocal conserved quantities have to be taken into account to construct an appropriate generalized Gibbs ensemble wouters ; pozsgay ; goldstein ; pozsgay1 ; pozsgay2 ; essler_mussardo_panfil ; ilievski1 ; ilievski2 ; doyon ; vidmar .
Sudden quenches have been studied for bosons in optical lattices experimentally BHdynamicsexp and theoretically BHdymamicstheory . The underlying Hamiltonian, the BoseHubbard model with nearestneighbour interactions, is known to be nonintegrable, and thus the dynamics expected to thermalize, but numerical studies, comprising DMRG in one dimension DMRG , tVMC in higher dimensions tVMC or numerical dynamical MFT dMFT indicate nonthermal behaviour for strong quenches. Nonequilibruim dynamics of the bosonic system in the cavitysetup has been studied recently and after very long times the dynamics is observed to become incoherent, which is explained that dissipation due to photon loss can take place EsslingerMetastab ; QMCtrap .
In this paper we study theoretically quantum manybody systems with competing short and long (infinite) range interactions. The system we consider is the quantum XXmodel, which is closely related to the BoseHubbard model. For hard core bosons, when the lattice sites have only single occupation the boson operators can be represented by spin operators. In the actual calculation the quantum XXmodel is put on a one dimensional lattice and we solve its ground state and quench dynamics exactly by freefermionic techniques. A brief account of our results translated into the bosonic language has been given by us recently BlassPRL .
The paper is organised as follows. In Sec.II we present the model and show, how the term with globalrange interaction is transformed to an effective field, the strength of which is calculated selfconsistently in the thermodynamic limit. In Sec.III we solve the onedimensional model exactly and in Sec.IV we calculate the quantum phases and phase transitions in the ground state. In Sec.V we study nonequilibrium dynamics of the model after a quench and properties of dynamical phase transitions are calculated. Our results are discussed in Sec.VI and detailed derivation of different results are put in the Appendices.
Ii The quantum XXmodel with global range interactions
Let us consider the quantum XXmodel with global range interactions in a onedimensional lattice, defined by the Hamiltonian:
(1) 
Here the are Pauli matrices at site , the nearest neighbour coupling constant and the strength of the transverse field are denoted by and , respectively. The last term of the r.h.s. represents a global range antiferromagnetic interaction of strength , favoring antiparallel orientation of the spins on even and odd lattice sites. It is expressed as the square of the staggered magnetization operator:
(2) 
We note that the Hamiltonian in Eq.(1) is equivalent to the BoseHubbard model for hard core bosons, which is explained in Appendix A. Here the global range interaction represents a cavity mediated long range interaction Cavity_review ; BHCavity1 ; BHCavity2 and therefore has immediate experimental relevanceScienceEsslinger ; Klinder2015 ; Landig2016 . Consequently our results on the XXmodel can be translated to hard core lattice bosons in 1d with cavity mediated long range interactions BlassPRL .
In the next step we transform the Hamiltonian in Eq.(1) in an equivalent form, by linearizing the globalrange interaction term in the thermodynamic limit. Following the steps of the derivation in the Appendix we arrive to the Hamiltonian:
(3) 
where has to be determined selfconsistently via
(4) 
is the average in the ground state of the system defined by . This condition is equivalent to the requirement that the groundstate energy of is minimal with respect of , which follows from the HellmannFeynmann theorem:
(5) 
At this point we make two comments about possible generalisation of the treatment of the cavityinduced globalrange interaction term. First, the equivalence of the two Hamiltonians and with (4) is valid for bipartite lattices in arbitrary dimensions in the thermodynamic limit . Second, using the technique of Appendix B a cavity induced globalrange interaction term can be linearised in the thermodynamic limit for other models, too, such as for the extended BoseHubbard model (with soft or hardcore bosons), as defined in Eq.(48) the Appendix A, and also to fermionic models as for instance free fermions on bipartite lattices.
Iii Freefermion solution
Using the JordanWigner transformation the Hamiltonian in Eq.(3) is transformed in terms of the fermion creation () and annihilation () operators in quadratic form:
(6) 
with and . In the next step we introduce the Fourier representation:
(7) 
where the values are in the range: . (For these are , , while for these are , with ) In terms of the Fourier operators the Hamiltonian assumes the form:
(8) 
which is separated into independent sectors.
The operators are diagonalized by the canonical transformation:
(9) 
with
(10) 
and . Then:
(11) 
with:
(12) 
Note that the energy of modes is symmetric to , thus we can restrict ourselves to the range: , however with . The inverse of Eq.(9) is given by:
(13) 
Iv Groundstate phase diagram
The energy of modes of the diagonalised Hamiltonian are for all , but the are positive, if . In the following we characterise a state by a wavenumber , so that for all and for and for . The energy per site is given by:
(14) 
Here the second term of the r.h.s. of the last equation can be expressed as , in terms of the elliptic integral of the second kind: .
The selfconsistency criterion is obtained from Eq.(14) through Eq.(5):
(15) 
where is the elliptic integral of the first kind. In the ground state . We note that using the representation of the staggered magnetization:
(16) 
the same selfconsistency criterion can be obtained through Eq.(4).
The stability of the selfconsistent solution depends on the sign of the secondderivative:
(17) 
The trivial solution, represents a (local) minimum, if
(18) 
which is satisfied for , with
(19) 
One can show similarly, that for the nontrivial selfconsistent solution, , is also a (local) minimum. This follows from the fact, that for the first and third terms at the r.h.s. of Eq.(17) cancel and the remaining second term is positive. We conclude that at a fixed value of there is always one stable selfconsistent solution, which is the trivial one, , in the first regime, , and it is the nontrivial one, , in the second regime, .
For fixed values of the parameters, , and , the ground state has the lowest energy, thus it is selected by the condition: . The ground state is characterised by its filling value, and the staggered magnetization, , which is calculated selfconsistently. In the following we calculate the minimal values of in the two regimes separately, and then comparing those we select . To get information about the behavior of we calculate its first two derivatives:
(20)  
(21) 
The value of the first derivative at the reference points: and are given by:
(22) 
In the first regime with , is a concave function, since . The minimal value of is at , if , thus . If , but at the same time then the minimal value is located in the interior of the first regime. Finally, if the minimum of in the first regime is at .
In the second regime with , is a convex function, since . This can be shown by differentiating the two sides of Eq.(15) which leads to the relations:
(23) 
Putting this into Eq.(21) we obtain the announced relation. In the second regime the minimum of can only be at the boundaries, either at or at . The dependence of the energy is shown in Fig.1 at different values of . Varying the parameter we explore the different phases and phase transitions.
The absolute minimum of can not be at , since it is an inflexion point, thus the possible ground states are of three types.

Ferromagnetic (F) ground state with and . Here there is a finite energy gap, , and the spinspin correlation function is zero.

Antiferromagnetic (AF) ground state with and . Here the energy gap is a finite, , and the spinspin correlation function decays exponentially. This can be shown exactly for the endtoend correlation function using the expression in Eq.(25).

XY ground state with and . This is a critical ground state, by changing the parameters the filling parameter, is continuously changing. The energy gap is vanishing, , and the spinspin correlation function decays algebraically. The endtoend correlations in Eq.(25) decays as .
We note that in our model no ground state with simultaneous XY and AFM order – corresponding to supersolid order in the equivalent hard core BH model – is realised, which would have and .
The phase transition between the AF state and the F or the XY state is of first order, there is a jump in the value of at the transition point. On the other hand transition between the XY and the F states is continuous. The phase diagram is shown in Fig.2. The metastability limit of the XYphase is given by: , which corresponds according to Eq.(22) .
The AFphase for small is extremly narrow in as can be seen in the upper panel of Fig.2. Its extension can be estimated from the analytical form of the metastability limit as . An estimate for the staggered magnetization follows by requiring , which from the last equation of (22) gives . This means that the jump of the staggered magnetization at the AF XY transition is extremely small for small and it goes to zero in a special, exponential form.
iv.1 Spinspin correlation function
We have calculated the spinspin correlation function, defined as
(24) 
which in the freefermionic description is given by a determinant of order . The simplest form is given for endtoend correlations, i.e. with and , and for free chains. In this case with sites the calculation is performed in Appendix C, which leads to the expression:
(25) 
Evidently in the ferromagnetic phase with we have , which follows also from symmetry. In the XYphase with and the sum in Eq.(25) with terms of alternating signs will result in an uncompensated term of , thus the endtoend correlations decay algebraically as
(26) 
Thus in the XYphase there is quasilongrange spin order.
Finally in the AFphase with and the expression in Eq.(25) is similar to the form of sine square deformationsine_square and this leads to an exponential correction in , thus the endtoend correlations decay exponentially:
(27) 
We have demonstrated this by evaluating Eq.(25) numerically. The correlation length is a monotonously increasing function of . In the limit it goes like:
(28) 
Here we can use the estimate of at the end of the previous section, which leads to , which grows exponentially, in somewhat similar way as at the KosterlitzThouless transition.
The correlation length can be calculated analytically in the opposite limit , when in the r.h.s. of Eq.(25) we perform a Taylor expansion:
(29) 
Here using the fact, that:
(30) 
we obtain
(31) 
wich defines a correlation length
(32) 
iv.2 Phasediagram of the nonlinearized Hamiltonian
We have checked the role of finitesize effects by solving the problem in the original form, given by the Hamiltonian in Eq.(1). In particular we have calculated the component of the magnetization defined by: , which is shown in Fig.3 for finite chains with and . Since is a conserved quantity, it could have only integer values in the groundstate of the system, therefore in Fig.3 there are possible discrete values of Also in finite systems the F, XY and AF phases are identified with , and , respectively. The finitesize phasediagrams and the values of the are qualitatively similar for finite values of , as well as in the thermodynamic limit, see in the lower panel of Fig.2. Somewhat larger finitesize corrections are present for larger values of close to the phaseboundary between the and the phases.
V Nonequilibrium dynamics after a quench
v.1 Preliminaries
Let us consider our model defined in Eq.(1) in which the parameters are suddenly changed at time : from to . The actual state of the system at is , the ground state of the initial Hamiltonian but its timeevolution is governed by the afterquench Hamiltonian, so that the state of the system at time is given by:
(33) 
In our calculation we have used an equivalent Hamiltonian in Eq.(3) in which the globalrange interaction term is linearised in the thermodynamic limit. This linearised Hamiltonian depends formally on the value of the staggered magnetization, which has to be calculated selfconsistently at . For infinitesimal times the time evolution operator can formally be treated in the thermodynamic limit in the same way as the partition function, as detailed at the end of Appendix B. Then the saddle point equation in Eq.(63) holds at each timestep, thus the effective linearized Hamiltonian governing the dynamics assumes the same form as in Eq.(3), however the parameter is replaced by a timedependent function , which satisfies the selfconsistency criterion:
(34) 
with
(35) 
One can easily show, that the total energy is conserved under the process. Indeed using the HellmannFeynmann theorem one obtains:
(36) 
where in the last step the selfconsistency equation in Eq.(34) is used.
Calculation of timedependent quantities in the fermionic basis is shown in Appendix D, here we shortly recapitulate the main steps of the derivation. During the quench at new set of freefermion operators are created, and , which are related to the original ones by a rotation, see in Eq.(76). The timedependent fermion operators, and are expressed with and through timedependent Bogoliubov parameters in Eq.(77). These generally complex parameters satisfy a set of differential equations in Eq.(80), which contain in Eq.(81) and has to be integrated with the known initial conditions at .
v.2 Numerical results
The calculation of the state of the system after time from the quench necessitates the integration of a set of linear differential equations of complex variables. This integration has been performed numerically by the fourthorder RungeKutta method and the stepsize was chosen appropriately to obtain stable results. We have checked by comparing results of nonequilibrium relaxation with sizes and , that finitesize effects are negligible until time . In the vicinity of nonequilibrium critical points where the timescale is divergent (see in Eqs.(40) and (42)) we went up to .
We note, that the term with the transverse field: commutes with the Hamiltonian: , therefore the wavefunction of a given state of does not depend on (but its energy naturally does). Consequently a quench from to with fixed does not modify the stationary properties of the system. In the following we concentrate on those quenches in which both and are kept fixed and only the parameter of the globalrange interaction term changes from to .
v.3 Quench from the AF phase
In this subsection the initial state we consider belongs to the ordered AF phase, thus and the staggered magnetization at is given by . First we consider the case when , thus the initial state is not fully antiferromagnetic, . As an example we choose , fix and vary for the final state. The timedependence of the staggered magnetization for different values of the reduced control parameter with are shown in the inset of Fig.4.
If the strength of the global range interaction term is reduced, , the staggered magnetization shows a fast decay, and after this initial period, characterised by the time and value of the absolute minima, and , respectively, has an oscillatory behaviour and after sufficiently long time it attains a stationary value . As a general trend is monotonously decreasing with , and for too strong quenches, , thus , vanishes, thus the system exhibits a nonequilibrium dynamical phase transition. The variation of with is shown in Fig.4: it vanishes linearly at the phasetransition point:
(37) 
We have observed similar behaviour of for different initial AF states, the nonequilibrium dynamical phasetransition is found numerically to satisfy the relation:
(38) 
We have calculated the nonequilibrium spinspin correlation function: after the quench at time . In the AF phase using periodic chains the equilibrium spinspin correlation function: has an exponential dependence, similarly to the endtoend correlation function in Sec.IV.1. This is illustrated in Fig.5. If the quench is performed to the regime with no dynamically generated AF order, i.e. with , then for sufficiently long time approaches a stationary behavior with an exponential decay, see in Fig.5. The correlation length increases with .
The behavior of changes, if the quench is performed to the region of dynamically generated AF order, i.e. for . In this case, as illustrated in the inset of Fig.5 changes sign and has an oscillatory dependence.
In the following we keep and concentrate on the behaviour of in the vicinity of the nonequilibrium dynamical phase transition. The numerical results for vs. in loglog scale are collected in Fig.6.
At the critical point the staggered magnetization for large time decays algebraically in an oscillatory fashion:
(39) 
Our numerical results indicate and , see in Fig.6.
In the AF ordered regime with the curves start to deviate from the critical one and have a minimum at a characteristic time , having a value . Close to the transition point these scale as:
(40) 
see in Fig.7.
After passing the minimum the dynamical staggered magnetization grows to a stationary value and start to oscillate in the form:
(41) 
in which the amplitude of the oscillations goes to zero as , which is illustrated in the inset of Fig.7.
The timescale in the stationary region, is in the same order of magnitude as , thus there is just one timescale in the problem, which is divergent at the dynamical phasetransition point:
(42) 
Since the prefactor, has only a weak dependence close to the critical point the curves in the main Fig.6 could be scaled together, which is shown in the inset.
We note that the staggered magnetization after passing the absolute minima has a strong revival, since the ratio of its values at the final (stationary) period and the initial (minimum) period is given by : , which is divergent as the transition point is approached.
We have checked, that the values of the scaling exponents are universal for . If, however the starting state is fully antiferromagnetic, i.e. and thus , than the appropriate scaling combinations in the stationary regime are the following:
(43) 
which can be illustrated by an appropriate scaling plot of the curves (not shown here). This difference is due to the fact, that for the Bogoliubovparameters are symmetric: and , which is not the case for .
v.4 Quench from the XY phase
In this subsection the initial state belongs to the XYphase thus we have and . If is exactly , then according to Eqs.(80) the Bogoliubov parameters decouple from each other and stays zero for . In the following we test the stability of this solution by adding a small perturbation, , to the staggered magnetization. Having fixed and we have quenched the system to various values of . For smaller values of the resulting staggered magnetization oscillates around the mean value of and its amplitude stays in the order of . If, however, we quench the system to , then grows within a time to a maximum value of and then oscillates between and with a period . This is illustrated in the right inset of Fig.8. We note, that similar type of macroscopic revival of an orderparameter has been observed also in Ref.Schuetz_Morigi .
This process represents a dynamical phasetransition separating a region in which the solution is stable from a region, in which the timeaverage value of the dynamically generated staggered magnetization is finite. The dynamical phasetransition coincides with the metastability line in the AF phase, thus dynamically generated staggered magnetization takes place only in such regions, in which no metastable XYtype solution exist.
Denoting the reduced controlparameter by the maximum value of the dynamically generated staggered magnetization vanishes as a power of , but at the same time the timescale, , is divergent. Our numerical results in Fig.9 are consistent with the asymptotic relations:
(44) 
Close to the dynamical phasetransition point the dynamical staggered magnetization follows the scaling form:
(45) 
which is illustrated in the left inset of Fig.8. Consequently the timeaverage value of the staggered magnetization for small behaves as:
(46)  
with .
For quenches more deep into the AF phase goes over a maximum and then decays for large as , but we still have , as shown in Fig. 9.
We have also checked the effect of the strength of the perturbation on the relaxation process. As shown in Fig.10 with decreasing the curves are simply shifted in time, thus stays the same but the timescales are shifted by an amount of
(47) 
Consequently any nonzero perturbation causes a measurable increase of the dynamical staggered magnetization if the quench is performed to .
For periodic boundary conditions we have measured the nonequilibrium spinspin correlation function after a quench protocol from to different values of . The equilibrium spinspin correlations in the XY phase are found to show an algebraic decay, similarly to the endto end correlations in Sec.IV.1, and the value of the decay exponent is consistent with the exact asymptotic relation: , see in Fig.11. If the quench is performed below the dynamical phasetransition point where the nonequilibrium spinspin correlation function is identical with . If, however after the quench there is a dynamically generated AF order, i.e. and , then the shape of is different from . The algebraic decay is preserved, but for the prefactor is different for even and odd distances between the spins. This is shown in Fig.11. We have also calculated the ratio at different times after the quench, which is shown in the inset of Fig.11. This ratio is different for even and odd distances between the spins but practically independent of the value of of the given parity. The larger the orderparameter, , the larger the difference between the ratios.
We can thus conclude that after a quench from the XYphase to the AFphase above the metastability line such a state is created, in which AF order and XY quasilongrange order coexist. The analogous state in the BoseHubbard model is the supersolid phase.