Quantum XX-model with competing short- and long-range interactions:
Phases and phase transitions in and out of equilibrium
We consider the quantum XX-model in the presence of competing nearest-neighbour and global-range interactions, which is equivalent to a Bose-Hubbard 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 XY-phase having quasi-long-range 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 AF-phase beyond a static metastability line AF order arises on top of remanent XY quasi-long-range order, which corresponds to dynamically generated supersolid state in the equivalent Bose-Hubbard model with hard-core bosons.
Recently there is an increased interest to study the phase diagram and non-equilibrium dynamics of quantum many-body systems with competing short- and long-range interactions. Experimentally, such systems have been realised with ultracold atoms in optical lattices inside a high-finesse optical cavityScience-Esslinger ; Klinder2015 ; Landig2016 . The strength of the short-range (on-site) 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 long-range interactions may result in a rich phase-diagram with exotic phases and interesting non-equilibrium dynamics.
Classical many-body systems with competing short- and ferromagnetic global-range interactions have been studied earlier and in the thermodynamic limit the order-parameter is obtained through a self-consistent treatment, like in mean-field models. Kardar1983 ; Mukamel2005 . Theoretical results for the phase diagrams of quantum many body systems with competing short-range and global-range interactions are rare and up to now confined to the aforementioned bosonic systemChen2016 ; Dogra2016 ; Niederle2016 ; Sundar2016 ; Panas2017 ; Flottat2017 .
The non-equilibrium 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 time-evolution 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 non-integrable quantum systems. Non-integrable 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 so-called generalized Gibbs ensemble, which includes all the conserved quantities of the system. Recently it has been observed that in certain models both local and quasi-local 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 BH-dynamics-exp and theoretically BH-dymamics-theory . The underlying Hamiltonian, the Bose-Hubbard model with nearest-neighbour interactions, is known to be non-integrable, and thus the dynamics expected to thermalize, but numerical studies, comprising DMRG in one dimension DMRG , t-VMC in higher dimensions tVMC or numerical dynamical MFT dMFT indicate non-thermal behaviour for strong quenches. Non-equilibruim dynamics of the bosonic system in the cavity-setup 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 Esslinger-Metastab ; QMC-trap .
In this paper we study theoretically quantum many-body systems with competing short- and long- (infinite) range interactions. The system we consider is the quantum XX-model, which is closely related to the Bose-Hubbard 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 XX-model is put on a one dimensional lattice and we solve its ground state and quench dynamics exactly by free-fermionic techniques. A brief account of our results translated into the bosonic language has been given by us recently Blass-PRL .
The paper is organised as follows. In Sec.II we present the model and show, how the term with global-range interaction is transformed to an effective field, the strength of which is calculated self-consistently in the thermodynamic limit. In Sec.III we solve the one-dimensional model exactly and in Sec.IV we calculate the quantum phases and phase transitions in the ground state. In Sec.V we study non-equilibrium 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 XX-model with global range interactions
Let us consider the quantum XX-model with global range interactions in a one-dimensional lattice, defined by the Hamiltonian:
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 anti-parallel -orientation of the spins on even and odd lattice sites. It is expressed as the square of the staggered magnetization operator:
We note that the Hamiltonian in Eq.(1) is equivalent to the Bose-Hubbard 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 ; BH-Cavity1 ; BH-Cavity2 and therefore has immediate experimental relevanceScience-Esslinger ; Klinder2015 ; Landig2016 . Consequently our results on the XX-model can be translated to hard core lattice bosons in 1d with cavity mediated long range interactions Blass-PRL .
In the next step we transform the Hamiltonian in Eq.(1) in an equivalent form, by linearizing the global-range interaction term in the thermodynamic limit. Following the steps of the derivation in the Appendix we arrive to the Hamiltonian:
where has to be determined self-consistently via
is the average in the ground state of the system defined by . This condition is equivalent to the requirement that the ground-state energy of is minimal with respect of , which follows from the Hellmann-Feynmann theorem:
At this point we make two comments about possible generalisation of the treatment of the cavity-induced global-range 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 global-range interaction term can be linearised in the thermodynamic limit for other models, too, such as for the extended Bose-Hubbard model (with soft- or hard-core bosons), as defined in Eq.(48) the Appendix A, and also to fermionic models as for instance free fermions on bipartite lattices.
Iii Free-fermion solution
Using the Jordan-Wigner transformation the Hamiltonian in Eq.(3) is transformed in terms of the fermion creation () and annihilation () operators in quadratic form:
with and . In the next step we introduce the Fourier representation:
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:
which is separated into independent sectors.
The operators are diagonalized by the canonical transformation:
and . Then:
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:
Iv Ground-state 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:
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: .
where is the elliptic integral of the first kind. In the ground state . We note that using the representation of the staggered magnetization:
the same self-consistency criterion can be obtained through Eq.(4).
The stability of the self-consistent solution depends on the sign of the second-derivative:
The trivial solution, represents a (local) minimum, if
which is satisfied for , with
One can show similarly, that for the non-trivial self-consistent 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 self-consistent solution, which is the trivial one, , in the first regime, , and it is the non-trivial 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 self-consistently. 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:
The value of the first derivative at the reference points: and are given by:
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:
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 spin-spin correlation function is zero.
Anti-ferromagnetic (AF) ground state with and . Here the energy gap is a finite, , and the spin-spin correlation function decays exponentially. This can be shown exactly for the end-to-end 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 spin-spin correlation function decays algebraically. The end-to-end 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 XY-phase is given by: , which corresponds according to Eq.(22) .
The AF-phase 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 Spin-spin correlation function
We have calculated the spin-spin correlation function, defined as
which in the free-fermionic description is given by a determinant of order . The simplest form is given for end-to-end 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:
Evidently in the ferromagnetic phase with we have , which follows also from symmetry. In the XY-phase with and the sum in Eq.(25) with terms of alternating signs will result in an uncompensated term of , thus the end-to-end correlations decay algebraically as
Thus in the XY-phase there is quasi-long-range spin order.
Finally in the AF-phase 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 end-to-end correlations decay exponentially:
We have demonstrated this by evaluating Eq.(25) numerically. The correlation length is a monotonously increasing function of . In the limit it goes like:
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 Kosterlitz-Thouless 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:
Here using the fact, that:
wich defines a correlation length
iv.2 Phase-diagram of the non-linearized Hamiltonian
We have checked the role of finite-size 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 ground-state 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 finite-size phase-diagrams 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 finite-size corrections are present for larger values of close to the phase-boundary between the and the phases.
V Non-equilibrium dynamics after a quench
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 time-evolution is governed by the after-quench Hamiltonian, so that the state of the system at time is given by:
In our calculation we have used an equivalent Hamiltonian in Eq.(3) in which the global-range 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 self-consistently 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 time-step, thus the effective linearized Hamiltonian governing the dynamics assumes the same form as in Eq.(3), however the parameter is replaced by a time-dependent function , which satisfies the self-consistency criterion:
One can easily show, that the total energy is conserved under the process. Indeed using the Hellmann-Feynmann theorem one obtains:
where in the last step the self-consistency equation in Eq.(34) is used.
Calculation of time-dependent 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 free-fermion operators are created, and , which are related to the original ones by a rotation, see in Eq.(76). The time-dependent fermion operators, and are expressed with and through time-dependent 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 fourth-order Runge-Kutta method and the step-size was chosen appropriately to obtain stable results. We have checked by comparing results of non-equilibrium relaxation with sizes and , that finite-size effects are negligible until time . In the vicinity of non-equilibrium critical points where the time-scale 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 global-range 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 time-dependence 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 non-equilibrium dynamical phase transition. The variation of with is shown in Fig.4: it vanishes linearly at the phase-transition point:
We have observed similar behaviour of for different initial AF states, the non-equilibrium dynamical phase-transition is found numerically to satisfy the relation:
We have calculated the non-equilibrium spin-spin correlation function: after the quench at time . In the AF phase using periodic chains the equilibrium spin-spin correlation function: has an exponential -dependence, similarly to the end-to-end 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 non-equilibrium dynamical phase transition. The numerical results for vs. in log-log scale are collected in Fig.6.
At the critical point the staggered magnetization for large time decays algebraically in an oscillatory fashion:
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:
see in Fig.7.
After passing the minimum the dynamical staggered magnetization grows to a stationary value and start to oscillate in the form:
in which the amplitude of the oscillations goes to zero as , which is illustrated in the inset of Fig.7.
The time-scale in the stationary region, is in the same order of magnitude as , thus there is just one time-scale in the problem, which is divergent at the dynamical phase-transition point:
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:
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 Bogoliubov-parameters 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 XY-phase 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 order-parameter has been observed also in Ref.Schuetz_Morigi .
This process represents a dynamical phase-transition separating a region in which the solution is stable from a region, in which the time-average value of the dynamically generated staggered magnetization is finite. The dynamical phase-transition coincides with the metastability line in the AF phase, thus dynamically generated staggered magnetization takes place only in such regions, in which no metastable XY-type solution exist.
Denoting the reduced control-parameter by the maximum value of the dynamically generated staggered magnetization vanishes as a power of , but at the same time the time-scale, , is divergent. Our numerical results in Fig.9 are consistent with the asymptotic relations:
Close to the dynamical phase-transition point the dynamical staggered magnetization follows the scaling form:
which is illustrated in the left inset of Fig.8. Consequently the time-average value of the staggered magnetization for small behaves as:
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 time-scales are shifted by an amount of
Consequently any non-zero perturbation causes a measurable increase of the dynamical staggered magnetization if the quench is performed to .
For periodic boundary conditions we have measured the non-equilibrium spin-spin correlation function after a quench protocol from to different values of . The equilibrium spin-spin correlations in the XY phase are found to show an algebraic decay, similarly to the end-to 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 phase-transition point where the non-equilibrium spin-spin 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 order-parameter, , the larger the difference between the ratios.
We can thus conclude that after a quench from the XY-phase to the AF-phase above the metastability line such a state is created, in which AF order and XY quasi-long-range order coexist. The analogous state in the Bose-Hubbard model is the supersolid phase.