Transport in out-of-equilibrium XXZ chains: non-ballistic behavior and correlation functions
Lorenzo Piroli, Jacopo De Nardis, Mario Collura, Bruno Bertini, Maurizio Fagotti
TTransport in out-of-equilibrium XXZ chains:non-ballistic behavior and correlation functions
Lorenzo Piroli, Jacopo De Nardis, Mario Collura, Bruno Bertini, and Maurizio Fagotti SISSA and INFN, via Bonomea 265, 34136, Trieste, Italy D´epartement de Physique, ´Ecole Normale Sup´erieure / PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France The Rudolf Peierls Centre for Theoretical Physics,Oxford University, Oxford, OX1 3NP, United Kingdom
We consider the nonequilibrium protocol where two semi-infinite gapped XXZ chains, initiallyprepared in different equilibrium states, are suddenly joined together. At large times, a generalizedhydrodynamic description applies, according to which the system can locally be represented byspace- and time- dependent stationary states. The magnetization displays an unusual behavior:depending on the initial state, its profile may exhibit abrupt jumps that can not be predicteddirectly from the standard hydrodynamic equations and which signal non-ballistic spin transport.We ascribe this phenomenon to the structure of the local conservation laws and make a prediction forthe exact location of the jumps. We find that the jumps propagate at the velocities of the heaviestquasiparticles. By means of time-dependent density matrix renormalization group simulations weshow that our theory yields a complete description of the long-time steady profiles of conservedcharges, currents, and local correlations.
The role of integrability in modern many-body quan-tum physics could hardly be overestimated, havingstrengthened our understanding of ground-state andthermal physics for the better part of the last century[1–4]. In the past two decades there has been an un-precedented and increasing interest in the nonequilibriumdynamics of isolated systems, and integrable models havebeen the ideal environment where to investigate out-of-equilibrium physics (see [5] for a volume of recent ped-agogic reviews on the subject). This is intimately con-nected with the development of new experimental tech-niques in cold atoms, which now allow us to probe almostunitary quantum nonequilibrium dynamics [6–16].One problem that has catalyzed an impressive amountof theoretical efforts is the characterization of the long-time steady state reached in an integrable system pre-pared in an out-of-equilibrium state. In the simplest in-stance, the system is initially prepared in a homogeneousequilibrium state, and perturbed by a sudden change inone of the Hamiltonian parameters (a so-called quantumquench [17]). It is now established that local observablesrelax and can be quantitatively described by a General-ized Gibbs Ensemble (GGE), a statistical ensemble builttaking into account all the local and quasi-local conservedoperators [18–20].More recently, increasing attention has been devoted tothe more general and complex case where the system isinitially prepared in a inhomogeneous state, for exampleby joining together two macroscopically different homo-geneous states. Initially, analytical understanding wasrestricted to either free systems [21–38] or conformallyinvariant models and Luttinger liquids [39–50], while gen-uinely interacting systems were mainly investigated nu-merically, or relying on ad hoc conjectures [51–58]. Abreakthrough came in the past year, with the introduc-tion of the so-called generalized hydrodynamics [59, 60].A significant number of studies have already been de-voted to further investigate some of its most interesting t ⇣ +1 ⇣ ⇣ +2 ⇣ . . .. . . ⇣ e L H +( h ) L S z e R H +( h ) R S z j FIG. 1. Pictorial representation of the quench protocol. Afterjoining together two semi-infinite XXZ chains, quasiparticleexcitations are created. Different types of quasiparticles movewith different maximum and minimum velocities ζ ± n . Theheaviest quasiparticles move with velocity ζ ∞ ( cf. Sec. II B 1). aspects [61–68]. Indeed, such approach is extremely flex-ible, allowing one to study, for instance, nonequilibriumdynamics in the presence of localized defects [69, 70], orof confining traps [61]. Furthermore, it has been shownthat the hydrodynamic picture not only gives an ex-act description at infinite length- and time- scales butcan also be a surprisingly good approximation at finitescales [65–67] or in the presence of small integrabilitybreaking terms [61]. These developments also boostedthe study of linear and non-linear transport in integrablesystems, which is a topic of long-standing interest andwith a long history [71–86]. Hydrodynamic approachesled to important results on several open questions, suchas the nature of spin and charge Drude weights [87, 88].In addition, related ideas were recently employed for thecomputation of the time evolution of the entanglemententropy after a global quench [89].According to the hydrodynamic picture, at large times t local subsystems at distance x from the junction reachdifferent stationary states depending on the “ray” x/t ,see Fig. 1. Stationary states describing observableson a fixed ray are characterized by appropriate GGEsor, equivalently, by appropriate distributions of quasi- a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p momenta (or rapidities ) of the elementary excitations.The derivation of the hydrodynamic equations pro-posed in [59] and [60] is based on the existence of a com-plete set of conservation laws. Remarkably, the conser-vation laws in the XXZ chain have a different structurein the gapped and gapless regimes [82, 90–95]. Only thelatter case, however, has been analyzed theoretically [60],so it is natural to wonder whether and why qualitativedifferences emerge in the space-time profiles of local ob-servables.In this work we settle this issue. We provide a de-tailed analysis of the space-time profiles of local con-served charges, currents and local correlations in thegapped regime. The most remarkable phenomena areobserved when the sign of the magnetization in the ini-tial state is not the same on the two sides of the junction.In that case, observables display a different behavior de-pending on how they transform under spin flip. In partic-ular, the magnetization exhibits abrupt jumps, which cannot be predicted directly from the hydrodynamic equa-tions obtained in [60]; the jumps are the signature of non-ballistic transport. We derive an equation that describesthe location of the jumps, and relate them to the velocityof the heaviest quasiparticles. This information comple-ments the structure unveiled in [60], so as to provide acomplete description of the long-time steady profiles ofall local observables in the gapped regime. Moreover, wediscuss the emergence of non-analyticities in the profilesof observables, revealing their connection with the quasi-particle content of the theory.The manuscript is laid out as follows. In Section I weintroduce the XXZ model and its thermodynamic Betheansatz (TBA) solution. In Section II we describe thequench protocol and review the hydrodynamic approach.Section III shows how to determine the jumps in observ-ables which are odd under spin flip. In Section IV weanalyze profiles of densities of charges, currents and cor-relation functions, and discuss the implications of thejumps for the spin transport. Section V contains ourconclusions. I. THE MODEL
We consider the XXZ spin-1 / H = L − (cid:88) j = − L (cid:2) s xj s xj +1 + s yj s yj +1 + ∆ s zj s zj +1 (cid:3) , (1)where s αj = σ αj ( α = x, y, z ) and σ αj are Paulimatrices. We assume periodic boundary conditions,namely s α L = s α − L . We focus on the case ∆ >
1, wherethe Hamiltonian (1) is gapped, and make use of theparametrization∆ = cosh η , η > . (2) We denote the energy density operator by e j , i.e. H = (cid:80) j e j . The eigenvalues of (1) can be constructedexactly by means of the Bethe ansatz method [2, 96]. Inthat framework, eigenstates are characterized in termsof rapidities λ ∈ [ − π/ , π/
2] of the so-called magnons,fulfilling some appropriate quantization conditions. Asthe theory is fully interacting, these particles can createbound states of any size n [2] ( n = 1 corresponds to un-bound magnons). In the thermodynamic limit L → ∞ ,the number of magnons diverges and the stationary stateswith nonnegative magnetization (we will come back tothis point later) are characterized by the densities of theirquasi momenta. For each bound state type there is a den-sity of rapidities (or root density) ρ n ( λ ). These densitiesare such that, in a large finite volume L , Lρ n ( λ )d λ givesthe number of magnonic bound states of length n withrapidities in the interval [ λ, λ + dλ ). The root densitiescan be thought of as the generalization to interactingmodels of the occupation numbers in free systems. Inthe XXZ chain with ∆ >
1, they characterize the ex-pectation value of any local operator which is invariantunder spin flip (see the discussion in Section I A).Along with the rapidity distributions ρ n ( λ ), it iscustomary to introduce the hole distribution functions ρ hn ( λ ); they describe allowed values of the rapiditieswhich are not occupied in the state. Distributions ofparticles and holes are related by the TBA equations ρ n ( λ ) + ρ hn ( λ ) = a n ( λ ) − ∞ (cid:88) m =1 ( T nm ∗ ρ m )( λ ) . (3)Here we introduced the so called scattering kernel T nm ( λ ), which reads as T nm ( λ ) = (1 − δ nm ) a | n − m | ( λ ) + 2 a | n − m | +2 ( λ )+ . . . + 2 a n + m − ( λ ) + a n + m ( λ ) , (4)where the functions a n ( λ ) are given by a n ( λ ) = 1 π sinh ( nη )cosh( nη ) − cos(2 λ ) . (5)The convolution between two functions is defined as fol-lows ( f ∗ g ) ( λ ) = (cid:90) π/ − π/ d µf ( λ − µ ) g ( µ ) . (6)Given the distributions ρ n ( λ ) and ρ hn ( λ ), the followingcombinations ρ tn ( λ ) = ρ n ( λ ) + ρ hn ( λ ) , ϑ n ( λ ) = ρ n ( λ ) ρ tn ( λ ) , (7)are usually called total root densities and filling func-tions, respectively. A. Conserved quantities and rapidity distributions
Being integrable, the model described by the Hamil-tonian (1) admits a macroscopically large set of localand quasi-local conserved charges { S z , Q ( n ) j } [90]. Here S z = (cid:80) (cid:96) s z(cid:96) indicates the total magnetization in the z direction and Q (1)1 + ∆ L is the Hamiltonian (1). Theexpectation values of these charges can be taken as thequantum numbers used to classify the eigenstates of theHamiltonian in the thermodynamic limit. As shown in[91, 92], there is a one-to-one correspondence between thedistribution of rapidities ρ n ( λ ) and the conserved quanti-ties. In the gapped regime, given the expectation valuesof the charges, the distributions ρ n ( λ ) read as ρ n ( λ ) = X + n ( λ ) + X − n ( λ ) − X n − ( λ ) − X n +1 ( λ ) . (8)Here n >
0, the quantities X n ( λ ) are the generatingfunctions of the expectation values of { Q ( n ) j } , and X [ ± ] n ( λ ) = X n ( λ ± iη/ { ρ n } , the expecta-tion values of the density q of a generic conserved charge Q ∈ { Q ( n ) j } is given by (cid:104) q (cid:105) = ∞ (cid:88) n =1 (cid:90) d λ q n ( λ ) ρ n ( λ ) , (9)where the functions q n ( λ ) are called single-particle eigen-values or bare charges. For example, the single-particleeigenvalues of the energy density (shifted by ∆ /
4) are q n ( λ ) = − π sinh( η ) a n ( λ ).For ∆ >
1, the conserved charges generated by X n ( λ )are invariant under spin flip O → Π O Π, withΠ = (cid:81) i σ xi , so the functions X n ( λ ) do not change if theyare computed in the state where all the spins are flipped.As a result, the stationary states with magnetization (cid:104) S z (cid:105) and −(cid:104) S z (cid:105) are described by the same set of rapidity dis-tributions ρ n ( λ ) [92, 93]. This is a crucial difference withrespect to the regime | ∆ | <
1, where also odd conservedcharges are generated by some X n ( λ ), and states withopposite magnetization are described by different distri-butions ρ n ( λ ) [87, 95]. For ∆ > ρ n ( λ ) are sufficientto characterize only the expectation values of the ob-servables which are even under spin flip, including theabsolute value of the magnetization |(cid:104) s z (cid:105)| = 12 − ∞ (cid:88) n =1 (cid:90) d λ n ρ n ( λ )= lim n →∞ (cid:90) dλ ρ hn ( λ ) ≥ , (10)where in the second step we used the TBA equations (3).Since, in the present understanding, S z is the only oddconserved charge, only an additional “bit” of informationis required to fully characterize the state. Specifically, itis widely accepted that it is sufficient to supplement theset of ρ n ( λ ) with a binary variable f = ± , which bearsinformation about the sign of the magnetization. Weindicate by | ρ, f (cid:105) a state with sign of the magnetizationequal to f and rapidity distributions given by ρ n ( λ ). Ex-pectation values in the state | ρ, f (cid:105) are denoted by (cid:104)·(cid:105) f . For operators O e , even under spin flip, we have (cid:104) O e (cid:105) f = (cid:104) O e (cid:105) + , (11)while for odd operators O o we have (cid:104) O o (cid:105) f = f (cid:104) O o (cid:105) + . (12)Note that (cid:104) s z (cid:105) + is that given in Eq. (10). B. Universal dressing equations and velocities
Let us consider the system in a large, finite volume L .In the thermodynamic limit, the state is described by theroot densities ρ n ( λ ), or, equivalently, by the filling func-tions ϑ n ( λ ). Excitations on this state can be constructedby injecting an extra string of size n with rapidity λ . Thisoperation induces a change in the expectation values ofthe conserved charges (cid:104) Q (cid:105) f → (cid:104) Q (cid:105) f + q dn ( λ ) . (13)Here q dn ( λ ) is called “dressed charge” and is an O ( L )deformation of the charge due to the presence of the newparticle of type n with rapidity λ . Its derivative withrespect to λ , q d (cid:48) n ( λ ), can be expressed as a linear integralequation which takes the following universal form q (cid:48) dn ( λ ) = q (cid:48) n ( λ ) − (cid:104) ∞ (cid:88) m =1 T nm ∗ ( q (cid:48) dm ϑ m ) (cid:105) ( λ ) . (14)Here q n ( λ ) is the bare charge ( cf . Eq (9)), i.e. , the chargecomputed with respect to the reference state with allspins up. We note that the bare charge p n ( λ ) for the mo-mentum is such that p (cid:48) n ( λ ) = 2 πa n ( λ ), so from Eq. (14)and Eq. (3) it follows p (cid:48) dn ( λ ) = 2 πρ tn ( λ ) . (15)We indicate with ε n ( λ ) the dressed energy of the particleexcitations. From the momentum and the energy we canextract the group velocity of a particle excitation of type n and rapidity λ [97] v n ( λ ) = ∂ε n ( λ ) ∂p dn ( λ ) = ε (cid:48) n ( λ )2 πρ tn ( λ ) . (16)We stress that the dressing equations (14) are valid forany integrable model (with diagonal scattering), providedthat its scattering kernel T nm ( λ ) is known. C. Currents
A current J q,(cid:96) is defined in terms of the density ofcharge q (cid:96) via the following continuity equation J q,(cid:96) +1 − J q,(cid:96) = i [ q (cid:96) , H ] . (17)Requiring J q,(cid:96) to vanish in the reference state, the opera-tor J q,(cid:96) is determined up to operators with zero expecta-tion value in any translationally invariant state. Impor-tantly, currents are generically not conserved and, aftera quantum quench, their expectation values undergo anon trivial time evolution.An important result of Refs [59, 60] was to suggesthow to compute the expectation value of a current in a“generic” stationary state in generic TBA solvable sys-tems. The result takes the simple form (cid:104) J q (cid:105) = ∞ (cid:88) n =1 (cid:90) d λ q n ( λ ) v n ( λ ) ρ n ( λ ) . (18)For ∆ >
1, this expression applies to the current of everycharge but S z . In particular, in the case of the energycurrent this expression can be rewritten as [60] (cid:104) J e (cid:105) f = ∞ (cid:88) n =1 (cid:90) d λ e n ( λ ) ρ n ( λ ) v n ( λ )= ∞ (cid:88) n =1 (cid:90) d λ q (1) n, ( λ ) ρ n ( λ ) = (cid:104) q (1)2 (cid:105) f , (19)where q (1)2 ,n ( λ ) is the bare charge corresponding to thesecond local conserved charge Q (1)2 . This is in accordancewith the well-known relation (cid:80) (cid:96) J e,(cid:96) = Q (1)2 . The spincurrent has to be supplemented with the information onthe sign of the magnetization ( cf . Sec. I A) and reads as (cid:104) J s (cid:105) f = f ∞ (cid:88) n =1 (cid:90) d λ n ρ n ( λ ) v n ( λ )= f n →∞ (cid:90) dλρ hn ( λ ) v n ( λ ) ; (20)in the second step, we used the equations (14). II. THE QUENCH PROTOCOL AND THEHYDRODYNAMIC EQUATIONSA. The Initial State
As discussed in the introduction, in this work we con-sider the nonequilibrium dynamics resulting from joiningtwo semi-infinite chains with different macroscopic prop-erties. In particular, we focus on the case where twochains are at thermal equilibrium with different values oftemperature and magnetic field. The initial state is thengiven by ρ = e − β L H L +( βh ) L S zL Z L (cid:79) e − β R H R +( βh ) R S zR Z R , (21)where the operators with the subscript L ( R ) are definedby restricting the sums of their density to the negative (positive) sites, while Z L and Z R are appropriate con-stants that ensure normalization.Starting from ρ , the region where local observablesare thermal remains macroscopically large: at any time t , as a consequence of the Lieb-Robinson bounds [98], faraway from the junction local observables are always de-scribed by thermal states. In integrable models, however,there is a region of width ∼ t around the origin where ob-servables are described by a family of non-thermal sta-tionary states, as pictorially represented in Fig. 1. Thecharacterization of this family, which was called LocallyQuasi-Stationary State in [69], is the subject of the nextsubsection. B. The Hydrodynamic Equations
In integrable models, like the XXZ spin-1/2 chain,there are stable quasiparticle excitations which propa-gate at different velocities and scatter elastically withone another. They are responsible for the propagation ofinformation throughout the system [97]. In many cases ofinterest, at large time- and length- scales, the quasipar-ticle excitations behave like free classical particles, andthe effects of the interactions can be taken into accountby letting the velocity of the quasiparticles to depend onthe state. In particular, if the initial state is the junc-tion of two homogeneous states, one can infer that, atlarge times, local observables moving on a certain “ray” ζ = j/t are characterized by a ζ -dependent steady state ρ s ( ζ ), c.f. Fig. 1. Indeed, different rays receive infor-mation from different quasiparticles. The hydrodynamicequations are derived under this assumption.We note that the state becomes equivalent to ρ s ( ζ )only when both the time and the distance from the junc-tion approach infinity at fixed ratio. By fixing the po-sition and increasing time, the observables explore theentire family of stationary states, eventually ending upin the state ρ s (0). The state ρ s (0) is known as nonequi-librium steady state (NESS).The state ρ s ( ζ ) has been characterized in Refs. [59,60], using that the expectation values of the local andquasi-local charges determine the rapidity distributions.Specifically, the root densities ρ n,ζ ( λ ) of the state ρ s ( ζ )have been shown to satisfy the following continuity equa-tion ζ∂ ζ ρ n,ζ ( λ ) = ∂ ζ (cid:104) v n,ζ ( λ ) ρ n,ζ ( λ ) (cid:105) , (22)where the velocity v n,ζ ( λ ) is given in Eq. (16). Here weare working in the limit of infinite times and distances atfixed ray ζ , where Eq. (22) is exact. One could also tryto extend this equation to describe finite-time dynamics.However, further terms would appear. In particular, wecan easily identify two kinds of finite-time corrections tothe naive finite-time version of (22) ∂ t ρ n,x,t + ∂ x (cid:104) v n [ ρ n,x,t ] ρ n,x,t (cid:105) = 0 . (23)The first type of corrections is related to the introductionof finite length scales, which make the thermodynamicdescription only approximate. While such correctionscould in principle be written in terms of root densities, itis difficult to estimate them in practice. The second kindof corrections are due to the fact that currents are notgenerically conserved. Indeed, as discussed in [68] for thespecific case of a noninteracting model, the expectationvalues of the currents take the form (18) only if the stateis stationary. These corrective terms can not be generi-cally written in terms of root densities. Nevertheless, forparticular classes of initial states, such corrections mightbe very small, leading to accurate quantitative predic-tions [65–67].Assuming that, for any λ and n , the equation ζ = v n,ζ ( λ ) has a unique solution (no exceptions areknown), the solution to Eq. (22) is most easily written interms of the filling functions ϑ n ( λ ) as follows ϑ n,ζ ( λ ) = ϑ n,R ( λ ) θ ζ − v n,ζ ( λ ) + ϑ n,L ( λ ) θ v n,ζ ( λ ) − ζ . (24)Here θ x is the Heaviside theta function, which is nonzeroand equal to 1 only if x >
0, while the “left” and “right”filling functions ϑ n,L ( λ ) and ϑ n,R ( λ ) are those charac-terizing the state at infinite distance from the junctionon the right and on the left hand side, respectively. Inour case, they correspond to thermal states with inversetemperatures β L and β R , and read as ϑ thn,L/R ( λ ) = 11 + e β L/R ε thn,L/R ( λ ) , (25)where ε thn,L/R ( λ ) is the thermal dressed energy, which sat-isfies ε thn,L/R ( λ ) = e n ( λ ) + h L/R + β − (cid:34) ∞ (cid:88) m =1 T nm ∗ ln (cid:104) e − βε thm,L/R (cid:105)(cid:35) ( λ ) . (26)We stress that the solution (24) is implicit: it dependson v n,ζ ( λ ), which in turn depends on the state. Theseequations can be generally solved numerically by simpleiterative schemes.Ref. [60] focused on the XXZ chain for | ∆ | <
1. Thederivation proposed is very general and can be appliedalso to the XXZ chain for ∆ >
1; (22) continues to holdalso there. To completely characterize the states, how-ever, there is a missing ingredient: we need to understandthe behavior of the sign f ζ ( cf . Sec. I A). Only once thebehavior of f ζ is known, the hydrodynamic descriptionbecomes complete. This problem will be addressed inthe next section.
1. Light Cones
The structure of the solution (24) allows us to infersome general properties of the profiles of the local observ-ables as a function of the ray ζ . To that aim, let us con-sider a ray ζ < min λ [ v n,ζ ( λ )], that is to say ζ < v n,ζ ( λ ) for any λ . From (24) it follows that the state at thatray has no information about the bound states of type n coming from the right hand side. Since we assumedthat the equation ζ = v n,ζ ( λ ) has a unique solution andlim ζ →−∞ v n,ζ ( λ ) is finite, we have ζ < v n,ζ ( λ ) ⇔ ζ < ¯ ζ n ( λ ) ∀ λ (27)where ¯ ζ n ( λ ) is the solution to the equation¯ ζ n ( λ ) = v n, ¯ ζ n ( λ ) ( λ ) ∀ λ . (28)Using (27), one can easily prove ζ < min λ [ v n,ζ ( λ )] ⇔ ζ < min λ ¯ ζ n ( λ ) (29)and min λ ¯ ζ n ( λ ) = ζ − n , (30)where ζ − n is the solution to the equation ζ − n = min λ [ v n,ζ − n ( λ )] . (31)We call ζ − n the “negative n -th light cone”. The ray ζ − n isthe one where the first particles of type n coming from theright become visible. Analogously, we define the “posi-tive n -th light cone” ζ + n as the solution to the equation ζ + n = max λ [ v n,ζ + n ( λ )] . (32)For ζ > ζ + n , there is no bound state of type n comingfrom the left hand side.When ζ is close to ζ ± n , the profiles of the local ob-servables have a typical square root behavior (cid:104) O (cid:105) ∼ θ ∓ ζ ± ζ ± n (cid:112) ∓ ζ ± ζ ± n . These non-analytic points are visiblein the numerical solutions of (22), see, e.g. , Figs 2 and3. This is essentially the same behavior seen in noninter-acting models, where, close the light cones, observablesdisplay universal properties that belong to the KPZ uni-versality class [26, 34].Finally, we note that, quite generally, the images ofthe velocities shrink in the limit of large n , and the ve-locities converge to a constant lim n →∞ v n,ζ ( λ ) = v ∞ ,ζ independent of λ . As a consequence, the following limitsexist lim n →∞ ζ + n = lim n →∞ ζ − n = ζ ∞ . (33)In the following, we analyze in detail the behavior ofspace-time profiles of local observables in correspondenceof this ray. In particular, we show that odd operatorsmight exhibit a discontinuous behavior depending on theinitial state, signalling the presence of non-ballistic trans-port. Further details on the fine structure of the profilesnear ζ ∞ will be presented elsewhere. III. BALLISTIC AND NON-BALLISTICTRANSPORT FOR ∆ ≥ The XXZ model (1) is integrable and, like any other in-tegrable model, is characterized by excitations that prop-agate ballistically. This allowed us to develop the hydro-dynamic theory presented in Section II B as a kinetic the-ory of particle excitations moving throughout the system.However, in some cases, symmetries may prohibit ballis-tic transport of certain quantities, leading to sub-ballistic(such as diffusive) behavior. This happens in the gappedregime | ∆ | ≥
1, where the spin-flip invariance of the rootdensities (8) results in a non-ballistic propagation of thespin degrees of freedom. Specifically, there is a region D , of size |D| ∼ t α with α <
1, where the magnetizationexperiences finite variations. Clearly, the magnetizationprofile as a function of the ray ζ = x/t becomes discon-tinuous at the ray ¯ ζ corresponding to the region D . Thedescription of the sublinear scaling region D goes beyondhydrodynamics, and most of the past investigations havebeen numerical [79, 99–101]. Here we show that the hy-drodynamic picture can provide useful information evenin such cases. In particular, we point out that a sub-ballistic region generically emerges by joining states withopposite total magnetization. Moreover, we demonstratethat such sub-ballistic behavior does not correspond al-ways to the NESS region ζ = 0, but it can be developedat finite rays ¯ ζ . A. The sign of the odd operators
Let us focus on the case where the two halves of the ini-tial state have magnetizations of opposite signs, f L f R < cf . I A). For our initial state (21), this situation is real-ized when h L h R <
0. Here we show that the profilesof all local operators O that are odd under spin-flip de-velop a discontinuity at a given ray ¯ ζ , as clearly visiblein Fig. 5. More precisely, we prove that the sign f ζ hasa single discontinuity at a ray ¯ ζ , whose position is fixedby the rapidity distributions ρ n ( λ ) of the left and rightstates.We start by considering the continuity equation for themagnetization ζ∂ ζ (cid:0) f ζ (cid:104) s z (cid:105) + ζ (cid:1) = ∂ ζ (cid:0) f ζ (cid:104) J s (cid:105) + ζ (cid:1) , (34)where (cid:104) s z (cid:105) + ζ and (cid:104) J s (cid:105) + ζ are the expectation values ina state with positive magnetization; they are given byEqs (10) and (20). From the continuity equation (22) forthe root densities it follows ζ∂ ζ (cid:104) s z (cid:105) + ζ = ∂ ζ (cid:104) J s (cid:105) + ζ . (35)Using this in (34) we find( (cid:104) J s (cid:105) + ζ − ζ (cid:104) s z (cid:105) + ζ ) ∂ ζ f ζ = 0 . (36)The solution to this equation is a piecewise constant func-tion of ζ equal to ± f ζ is a sign), which can be written as f ζ = f R θ ζ − v zζ + f L θ v zζ − ζ . (37)Here we have used that the equation ζ = v zζ ≡ (cid:104) J s (cid:105) + ζ (cid:104) s z (cid:105) + ζ , (38)has a unique solution. This can be proved by integratingthe continuity equation (35), which gives ζ − v zζ = (cid:82) ζ ¯ ζ (cid:104) s z (cid:105) + ζ (cid:104) s z (cid:105) + ζ , (39)where we called ¯ ζ a zero of ζ − v zζ . Since, by definition, (cid:104) s z (cid:105) + ≥
0, the right hand side is equal to zero only for ζ = ¯ ζ , that is to say, the solution to (38) is unique.The solution ¯ ζ has a nice interpretation in terms oflight cones ( cf . Sec. II B 1). Considering the velocity v zζ and using the identities (10) and (20), we have v zζ = lim n →∞ (cid:82) d λ v n,ζ ( λ ) ρ hn,ζ ( λ ) (cid:82) d λ ρ hn,ζ ( λ ) . (40)As observed in Section II B 1, the images of the velocitiesshrink in the limit of large n ( cf . Sec II B 1), thus we find v zζ = v ∞ ,ζ . (41)The solution ¯ ζ to (38) is then identified with the accu-mulation point ζ ∞ for the light cones.Eqs (37) and (41), together with (24), fully character-ize the state in the hydrodynamic limit. Despite the no-tion of f is outside the TBA description, Eq. (41) suggeststhat the information about the sign of the odd operatorsis carried by the heaviest bound state. This result is notsurprising if one looks at the behavior of the spin den-sity and related current in the gapless regime ( | ∆ | < π/n ), with n an integernumber. In that case there are n species of excitationsand the information about the sign of the magnetizationis encoded in the last two species [95]. This can be seenin Fig. 6, in the cases ∆ = cos( π/
3) and ∆ = cos( π/ → − , n approaches infinity, and the last two species are sentto infinity. If this property does not break down in thegapped regime, the sign of the odd operators should notchange before the light cones of the heaviest bound states.Since, for ∆ >
1, the corresponding velocities approacha constant, ζ ∞ has to be exactly the ray where the signchanges.Remarkably, the sign of the front’s velocity can giveglobal information about the magnetization profile. Forexample, if the front moves towards the side with largermagnetization (in modulus), the absolute value of themagnetization can not be monotonous inside the lightcone. This can be proved by reductio ad absurdum . Letus assume that the absolute value of the magnetizationis smaller on the right hand side, so the front is propa-gating to the left, i.e. it has negative velocity. If (cid:104) s z (cid:105) + ζ ismonotonous, using the continuity equation (35), we have0 ≥ | ζ | ∂ ζ (cid:104) s z (cid:105) + ζ = sgn ζ ∂ ζ (cid:104) J s (cid:105) + ζ . (42)Integrating this equation from −∞ to the accumulationpoint ¯ ζ gives (cid:104) J s (cid:105) +¯ ζ ≥ , (43)where we used that the current outside the light cone iszero. The inequality in (43) can not be satisfied because (cid:104) J s (cid:105) +¯ ζ has the sign of ¯ ζ = v z ¯ ζ , which was negative byassumption ( cf . Eq. (38)). IV. RESULTS
In this section we elaborate on our predictions for theprofiles of local observables as a function of the ray ζ andshow a comparison with time-dependent density matrixrenormalization group (tDMRG) simulations. The pre-dictions are obtained by taking the expectation value oflocal observables in the state ρ s ( ζ ), which we representmicrocanonically by | ρ ζ , f ζ (cid:105) , where ρ ζ and f ζ are com-puted by first solving (22) and then (37).The tDMRG simulations are obtained for finite latticesof L sites, with L ∈ [80 , (i) We prepare each half-chain in the mixed productstate (cid:81) j e ( βh ) L/R s zj . In terms of locally purified matrixproduct states (MPS), such a state only needs a two-dimensional ancilla and an auxiliary bond dimension χ = 1. (ii) We implement imaginary time evolution us-ing second-order Trotter decomposition of the operator ∝ e − β L/R H , with imaginary time-step dβ = 10 − . (iii) We evolve both the initial left and right mixed productstates up to the desired temperatures.After joining together the two open chains, the sys-tem is unitarily evolved using second-order Trotter de-composition of the evolution operator, with time-step dt = 10 − . During the time evolution, the bond di-mension of the MPS is dynamically updated, up to amaximum value χ max = 200. For this reason, the max-imum time we can reach keeping the accumulated errorreasonably small is t max (cid:39) ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●▲ - � - � � � ����������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ������������ - � - � � � � - ���� - ����� - � - � � � � - ���� - ����� FIG. 2. Space-time profiles of densities and currents of spinand energy. Solid black lines display the theoretical predic-tions, while points correspond to the exact time evolutioncomputed by tDRMG simulations up to times t = 20. Verti-cal dashed lines represent positive and negative light cones ofthe different n -quasiparticle bound states, see Sec. II B 1. Thecorresponding rays ζ = ζ ± n , with n = 1 , ,
3, are displayed asgray dashed lines, while the black dashed line corresponds tothe largest string ζ = ζ ∞ . A. Homogeneous magnetization signs: light cones
Let us start by considering the case where the sign ofthe magnetization is homogeneous in the initial state and f ζ is constant throughout the light cone. In this settingthe qualitative behavior of the space-time profiles doesnot differ much from the one in the gapless regime.In Fig. 2 we report the space-time profiles of localobservables after the sudden junction of two infinite-temperature states with different (but positive) magne-tizations. One immediately sees that the profiles are notsmooth, presenting a number of cusps. These are thenon-analytic points discussed in Sec. II B 1, and their po-sitions { ζ ± n } can be predicted by solving Eqs (30) and(32). As discussed in Sec. II B 1, these points have a nat-ural interpretation in terms of moving quasiparticles: ζ + n and ζ − n correspond to the rays where the quasiparticlesof species n , coming respectively from the left and fromthe right, become visible.The first light cone is where the profiles begin to devi-ate from a constant function. This ray corresponds to thevelocity of the fastest particles (the unbound magnons inour case). Note that, since the system is interacting, themaximal velocities on the two sides are generically differ-ent from one another. This is the case for the data re-ported in Fig. 2. As the dispersion law of quasi-particlesis smooth, the profiles are expected to remain smoothbetween two consecutive light cones.Cusps are also present in the gapless regime studied in[60]; depending on the initial state, they can be more orless marked.As the rapidity distributions ρ n ( λ ) completely charac-terize the state, the solution to the hydrodynamic equa-tion (22) allows us to investigate further light-cone prop-erties, going beyond the analysis of conserved chargesand currents. To that aim, we use some recently de-veloped formulae [102] for the expectation values of lo-cal observables in generic eigenstates of the gapped XXZHamiltonian. In particular, we have computed nearestand next-to-nearest neighbor correlations inside the light-cone. Our results are reported in Fig. 3. Once again,cusps are clearly visible. We also observe an interesting,non-monotonic behavior of transverse correlators. Fig. 3also displays data from tDMRG simulations, which arefound to be in very good agreement with our predictions,further corroborating the validity of our results. ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●▲ - � - � � � � - �������������������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � � - ���� - ������������ ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ����������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ����������� FIG. 3. Space-time profiles of local correlators, same nota-tions as in Fig. 2. Note that the absolute value of correlatorsalong the x -direction is two orders of magnitude smaller thanthat along the z -direction. In the former case the visible smallripples on the theoretical curves are numerical artifacts. Finally, before turning to the next section, we pro-vide a dedicated analysis of the celebrated NESS en-ergy current, corresponding to (cid:104) J e (cid:105) ζ =0 . Fig. 4 showsits value as a function of the anisotropy ∆, in the casewhere the two semi-infinite chains are initially preparedat different temperatures and with vanishing magneticfields. The energy current has a non-monotonic behaviorin ∆, reaching a peak when ∆ ∼ min (cid:0) β − R , β − L (cid:1) . Fur- ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲●▲ � � �� �� ���� - � �� - � ��������������� FIG. 4. NESS energy current (cid:104) J e (cid:105) ζ =0 as a function of theanisotropy ∆. The initial state is prepared by joining to-gether two semi-infinite chains with vanishing magnetic fieldand different temperatures. thermore, for the chosen values of the initial parameters,the maximum is reached for ∆ >
1. The current is al-ways seen to vanish exponentially for ∆ → ∞ , as onecan clearly see from the logarithmic plot in Fig. 4. As afunction of the temperatures, it approximately behavesas ∼ exp (cid:2) − ∆ min (cid:0) β R , β L (cid:1) / (cid:3) . B. Heterogeneous magnetization signs: spin-jumps
We now turn to presenting our results for the casewhere the semi-infinite spin chains are initially preparedin equilibrium states with different magnetization signs f R = − f L ≡ − f .In light of the discussion in Sec. III A, we expect theobservables that are odd under spin-flip to abruptly fliptheir sign at ζ = ζ ∞ . This is nicely displayed in Fig. 5,where we reported our theoretical predictions and nu-merical data from tDMRG simulations.In order to test our predictions of the jumps againstnumerics, we proceed as follows. We fix a local observable O i and compute its profiles starting from two differentinitial states ρ (1)0 and ρ (2)0 . These states are chosen todiffer only in the sign of the magnetic field on the right.For the first state we have f L = f and f R = f , while, forthe second one, f L = f and f R = − f . We then define theratio R O f ,ζ ( t ) ≡ tr( O ζt ( t ) ρ (2)0 )tr( O ζt ( t ) ρ (1)0 ) . (44)This ratio is such thatlim t →∞ R O f ,ζ ( t ) = (cid:40) sgn( ζ ∞ − ζ ) Π O i Π = − O i O i Π = O i . (45)We remark that the analytic calculation of the profile of O i is not required to test this prediction; this allows usto consider also observables for which we are not able tocalculate the profile. For example, in Fig. 5 we also report R O f ,ζ ( t ) for O i = σ zi σ zi +1 σ zi +2 . We see from the figurethat we are able to successfully test (45) against tDMRGdata, even though no formula involving the root densitiesis currently available for computing the expectation valueof this operator.In all the cases considered, the tDMRG simulationsare compatible with our predictions, but the correctionsare not always small. In particular, a slow, sub-ballisticbehavior is expected at the discontinuity of the profiles,which contributes to the presence of large finite-time ef-fects. As a result, the tDMRG simulations can not reachsufficiently long times to observe an actual discontinuousbehavior. We ascribe the differences between our pre-dictions and the tDMRG data to such numerical prob-lems; our analysis of how the tDMRG data approachtheir asymptotic values supports that conclusion. ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●▲ - � - � � � ����������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � � - ��� - ������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � � - �� - ��������� FIG. 5. Space-time profiles of energy and magnetization(top) and spin-spin correlation functions (bottom), same no-tations as in Fig. 2. The function R σ z σ z is computed asthe ratio between two different profiles: the first is the profileof the correlator (cid:104) σ zi σ zi +1 (cid:105) obtained by joining thermal stateswith µ L = 1, β L = 0 and µ R = − β R = 0; the second isthe profile for µ L = 1, β L = 0 and µ R = +2, β R = 0. Theplot for R σ z σ z σ z is obtained analogously from the correlator (cid:104) σ zi σ zi +1 σ zi +2 (cid:105) . Note that the odd operators show a genuinediscontinuity at ζ ∞ ∼ − .
086 (black vertical dashed line).
The abrupt jumps in the profiles of odd observablesdisplayed in Fig. 5 find no correspondence in the gap-less regime. It is then important to understand how suchdiscontinuities arise as the value of the anisotropy is con-tinuously varied from ∆ < >
1. Some data arereported in Fig. 6. We see that, while the profiles remainalways continuous for ∆ <
1, they become increasinglysharp as ∆ is increased, finally developing a discontinuity - � � �������� - � � � - ��� - ������� - � � � - ���� - ���� - ����� - � � ������������� FIG. 6. Space-time profiles of densities and currents of spinand energy. Different plots correspond to different values of∆, in the gapless regime ∆ = cos π/(cid:96) , with (cid:96) = 3 , .
2. The small circles on top ofthe profiles indicate the positions of the light cones ζ ± n foreach values of ∆. Note that the number of light cones in thegapless regime is finite as the number of species is also finite.In the gapped regime instead there is an infinite number oflight cones converging to the ray ζ ∞ , where the magnetizationdensity and the spin current change sign. at ∆ = 1. C. Zero to finite magnetization: sharp front
In this section we finally consider the situation whereone of the two semi-infinite chains (say, the left one) isinitially prepared in a thermal state with vanishing mag-netic field, while the other (the right one) has a non-zero magnetic field. This is a limiting case of the onespresented in the previous subsections. For this quenchprotocol, the long-time magnetization profiles have def-inite sign as a function of the ray ζ . Accordingly, theprofiles of all local observables are simply obtained fromthe solution to the hydrodynamic equation (22), in anal-ogy to the situation discussed in Sec. IV A. In this case,however, the solution displays some interesting proper-ties which are worth discussing in a detailed fashion.The first example is a problem of release into the vac-uum. The right part of the system is prepared in thestate with all spins up (the vacuum), while the left partis in an infinite temperature state with vanishing mag-netic field. The numerical solution to the hydrodynamicequations (22) is displayed in Fig. 7. We clearly see that0 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●▲ - � - � � � ������� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ��������� - � - � � � � - ��� - ��� - ���� - � - � � � � - ���� - ���� - ���� FIG. 7. Space-time profiles of densities and currents of spinand energy, same notations as in Fig. 2. Remarkably, we seethat the leftmost light cones of the magnetization and of theenergy profiles do not coincide. This is due to the specialproperties of the initial state, as explained in detail in themain text. the leftmost light cones of the magnetization and energyprofiles do not coincide. This remarkable property canbe seen as a corollary of our theory on the sign of theodd operators.In order to show this, we consider the two situa-tions where tiny magnetic fields, respectively positive( h L = h (cid:15) ) and negative ( h L = − h (cid:15) ), are initially turnedon in the left semi-infinite chain. On the left hand sideof the first light cone ζ − the magnetization will be non-vanishing, (cid:104) s z (cid:105) ±−∞ = ± (cid:15) . By integrating the continuityequation (42) for the magnetization from ζ − to ζ +1 , wefind ± (cid:90) ζ ∞ ζ − d z (cid:104) s z (cid:105) + z + (cid:90) ζ +1 ζ ∞ d z (cid:104) s z (cid:105) + z = ζ +1 (cid:104) s z (cid:105) + ζ +1 ∓ ζ − (cid:15) , (46)where we used (37) and that the spin current is zerooutside the light cone. Taking the difference between thetwo cases gives lim (cid:15) → (cid:90) ζ ∞ ζ − d z (cid:104) s z (cid:105) + z = 0 . (47)Since (cid:104) s z (cid:105) + ζ is nonnegative, (47) implieslim h (cid:15) → (cid:104) s z (cid:105) + ζ = lim (cid:15) → (cid:104) s z (cid:105) + ζ = 0 ∀ ζ < ζ ∞ . (48)It is now reasonable to assume that the magnetizationprofile for h L = 0 can be obtained as the limit h (cid:15) → ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ ●▲ - � - � � � � - ��� - ���� ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲▲ ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ - � - � � � ���������������� - � - � � � ����������������� - � - � � � ������������ FIG. 8. Space-time profiles of densities and currents of spinand energy, same notations as in Fig. 2. The magnetizationexhibits a non-monotonic behavior, which is naturally inter-preted as a thermoelectric effect. of the profile where the left magnetic field is positive butsmall. In fact, this is actually implicit in the numericalsolution to the hydrodynamic equation (22). This simpleargument shows that the magnetization profile for h L = 0is vanishing for all the rays ζ smaller than ζ ∞ .This is a general property, and is observed every timethe initial state has vanishing magnetization on one of itstwo halves. For example, this is also the case displayedin Fig. 8, where the right magnetic field is finite.In Fig. 8, the magnetization profile exhibits an inter-esting non-monotonic behavior. The latter is naturallyinterpreted as a thermoelectric effect and is observed alsowhen the initial halves of the chain have the same non-vanishing magnetization but one part is much colder thanthe other.Finally, we point out that the magnetization profile inFig. 8 seems to develop a discontinuity at the accumula-tion point ζ ∞ . Our numerical analysis of the profiles withincreasing right magnetic fields seems to suggest that thefront could in fact be continuous, albeit extremely sharp.Near the accumulation point ζ ∞ , the profile varies veryquickly over a region that approaches zero in the limitwhere the right magnetic field is sent to zero. This couldbe at the basis of the apparent discontinuity displayedin Fig. 8.1 V. CONCLUSIONS
We have considered the time evolution of local ob-servables after the junction of two semi-infinite XXZchains initially prepared in thermal equilibrium withdifferent temperatures and magnetic fields. We focusedon the gapped regime | ∆ | >
1, where interactions arelarger in the direction of the anisotropy. Our analysisis complementary to the one of [60], where the gaplesscase | ∆ | < ζ ∞ , can not be predicted directly from thehydrodynamic equations satisfied by the root densities.Our main result is to derive an equation that describesthe ray ζ ∞ , completing the characterization of thelong-time steady profiles of all the local observables.It would be interesting to analyze in greater detail thesub-ballistic behavior around ζ ∞ , as done in [100] for ζ ∞ = 0. This is a challenging numerical problem, whichwe leave to future investigations.The form of the hydrodynamic equations employed inour work is extremely general and is expected to hold inall Bethe ansatz integrable models or in models wherestable particle excitations can be constructed [103]. Inparticular, we expect discontinuities in the space-timeprofiles of odd operators in every model with root den-sities invariant under some discrete symmetry. This is, for instance, the case of the Hubbard model [104], wherethe root densities are in one-to-one correspondencewith spin-flip and charge-flip invariant commuting fusedtransfer matrices [87, 105].Finally, we note that inhomogeneous quantumquenches prove to be an excellent setting where tostudy the particle content of the model: the space-timeprofiles of local observables can be used as effective“spectroscopes” of collective excitations, as the ap-pearance of singularities in the profiles of the localobservables is connected with the presence of morespecies of quasiparticles. This method is complementaryto others suggested for homogeneous quenches, based,for instance, on the computation of local correlations[106] or entanglement entropy and mutual information[89, 107]. ACKNOWLEDGMENTS
We are grateful to J.-S. Caux for drawing our atten-tion to the effects of the bound states in the profiles oflocal observables. We thank Vincenzo Alba, BenjaminDoyon, Enej Ilievski, and Herbert Spohn for stimulatingdiscussions.B.B. acknowledges the financial support by the ERCunder Starting Grant 279391 EDEQS. J.D.N. andM.F. acknowledge support by LabEx ENS-ICFP:ANR-10-LABX-0010/ANR-10-IDEX-0001-02 PSL*. M.C.acknowledges support by the Marie Sklodowska-CurieGrant No. 701221 NET4IQ.
Author contributions.
L. Piroli, J. De Nardis, B.Bertini, and M. Fagotti contributed equally to the de-velopment of the theory; M. Collura took care of thetDMRG simulations. In addition, L. Piroli and J. DeNardis took charge of the numerical solutions to the hy-drodynamic equations; B. Bertini and M. Fagotti playeda major role in the organization of the work and in therefinement of the theory. [1] B. M. McCoy and T. T. Wu,
The Two-DimensionalIsing Model , Harvard University Press (1973);R. J. Baxter,
Exactly Solvable Models in Statistical Me-chanics , Academic Press (1982);B. Sutherland,
Beautiful Models
World Scientific (2004).[2] M. Takahashi,
Thermodynamics of one-dimensionalsolvable models , Cambridge University Press (1999).[3] M. Gaudin,
La fonction d’onde de Bethe , Masson(1983);M. Gaudin (translated by J.-S. Caux),
The Bethe wavefunction
Cambridge University Press (2014).[4] V.E. Korepin, N.M. Bogoliubov and A.G. Izergin,
Quantum inverse scattering method and correlationfunctions , Cambridge University Press (1993). [5] P. Calabrese, F. H. L. Essler, and G. Mussardo, J. Stat.Mech. (2016) 64001.[6] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[7] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac,and M. Rigol, Rev. Mod. Phys. , 1405 (2011).[8] A. Polkovnikov, K. Sengupta, A. Silva, and M. Ven-galattore, Rev. Mod. Phys. , 863 (2011).[9] T. Langen, T. Gasenzer, and J. Schmiedmayer, J. Stat.Mech. (2016) 64009.[10] T. Kinoshita, T. Wenger, and D. S. Weiss, Nature ,900 (2006).[11] S. Hofferberth, I. Lesanovsky, B. Fischer, T. Schumm,and J. Schmiedmayer, Nature , 324 (2007). [12] M. Gring, M. Kuhnert, T. Langen, T. Kitagawa, B.Rauer, M. Schreitl, I. Mazets, D. A. Smith, E. Demler,and J. Schmiedmayer, Science , 1318 (2012).[13] T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Che-neau, I. Bloch, and C. Gross, Nature , 76 (2013).[14] T. Langen, R. Geiger, M. Kuhnert, B. Rauer, and J.Schmiedmayer, Nature Phys.
640 (2013).[15] R. Geiger, T. Langen, I. E. Mazets, and J. Schmied-mayer, New J. Phys. , 053034 (2014).[16] T. Langen, S. Erne, R. Geiger, B. Rauer, T. Schweigler,M. Kuhnert, W. Rohringer, I. E. Mazets, T. Gasenzer,and J. Schmiedmayer, Science , 207 (2015).[17] P. Calabrese and J. Cardy, Phys. Rev. Lett. , 136801(2006);P. Calabrese and J. Cardy, J. Stat. Mech. (2007)P06008.[18] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii,Phys. Rev. Lett. , 050405 (2007).[19] L. Vidmar and M. Rigol, J. Stat. Mech. (2016) 64007.[20] F. H. L. Essler and M. Fagotti, J. Stat. Mech. (2016)64002.[21] T. Antal, Z. R´acz, A. R´akos, and G. M. Sch¨utz, Phys.Rev. E , 4912 (1999).[22] W. H. Aschbacher and C.-A. Pillet, J. Stat. Phys. ,1153 (2003).[23] W. H. Aschbacher and J.-M. Barbaroux, Lett. Math.Phys. , 11 (2006).[24] T. Platini and D. Karevski, J. Phys. A: Math. Theor. , 1711 (2007).[25] J. Lancaster and A. Mitra, Phys. Rev. E , 61134(2010).[26] V. Eisler and Z. R´acz, Phys. Rev. Lett. , 60602(2013).[27] A. De Luca, J. Viti, D. Bernard, and B. Doyon, Phys.Rev. B , 134301 (2013).[28] M. Collura and D. Karevski, Phys. Rev. B , 214308(2014).[29] V. Eisler and Z. Zimbor´as, New J. Phys. , 123020(2014).[30] M. Collura and G. Martelloni, J. Stat. Mech. (2014)P08006.[31] A. De Luca, G. Martelloni, and J. Viti, Phys. Rev. A , 21603 (2015).[32] B. Doyon, A. Lucas, K. Schalm, and M. J. Bhaseen, J.Phys. A: Math. Theor. , 95002 (2015).[33] J. Viti, J.-M. St´ephan, J. Dubail, and M. Haque, EPL , 40011 (2016).[34] N. Allegra, J. Dubail, J.-M. St´ephan, and J. Viti, J.Stat. Mech. (2016) 053108.[35] M. Kormos and Z. Zimbor´as, J. Phys. A: Math. Theor. , 264005 (2017).[36] B. Bertini, Phys. Rev. B , 75153 (2017).[37] M. Kormos, arXiv:1704.03744 (2017).[38] G. Perfetto and A. Gambassi, Phys. Rev. E , 012138(2017).[39] S. Sotiriadis and J. Cardy, J. Stat. Mech. (2008) P11003.[40] P. Calabrese, C. Hagendorf, and P. L. Doussal, J. Stat.Mech. (2008) P07013.[41] M. Mintchev, J. Phys. A: Math. Theor. , 415201(2011).[42] M. Mintchev and P. Sorba, J. Phys. A: Math. Theor. , 95006 (2013).[43] B. Doyon, M. Hoogeveen, and D. Bernard, J. Stat.Mech. (2014) P03002. [44] D. Bernard and B. Doyon, J. Phys. A: Math. Theor. ,362001 (2012).[45] D. Bernard and B. Doyon, Ann. Henri Poincar´e , 113(2015).[46] D. Bernard and B. Doyon, J. Stat. Mech. (2016) 33104.[47] D. Bernard and B. Doyon, J. Stat. Mech. (2016) 64005.[48] E. Langmann, J. L. Lebowitz, V. Mastropietro, and P.Moosavi, Commun. Math. Phys. , 551 (2017).[49] J.-M. St´ephan, J. Dubail, P. Calabrese, and J. Viti, Sci-Post Physics , 2 (2017).[50] J. Dubail, J.-M. St´ephan, and P. Calabrese,arXiv:1705.00679 (2017).[51] D. Gobert, C. Kollath, U. Schollw¨ock, and G. Sch¨utz,Phys. Rev. E , 36102 (2005).[52] T. Sabetta and G. Misguich, Phys. Rev. B , 245114(2013).[53] A. De Luca, J. Viti, L. Mazza, and D. Rossini, Phys.Rev. B , 161101 (2014).[54] V. Alba and F. Heidrich-Meisner, Phys. Rev. B ,75144 (2014).[55] O. Castro-Alvaredo, Y. Chen, B. Doyon, and M.Hoogeveen, J. Stat. Mech. (2014) P03011.[56] A. Biella, A. De Luca, J. Viti, D. Rossini, L. Mazza,and R. Fazio, Phys. Rev. B , 205121 (2016).[57] X. Zotos, arXiv:1604.08434 (2016).[58] L. Vidmar, D. Iyer, and M. Rigol, Phys. Rev. X , 21012(2017).[59] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura,Phys. Rev. X , 41065 (2016).[60] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti,Phys. Rev. Lett. , 207201 (2016).[61] B. Doyon and T. Yoshimura, SciPost Physics , 14(2017).[62] B. Doyon, H. Spohn, and T. Yoshimura,arXiv:1704.04409 (2017).[63] B. Doyon and H. Spohn, arXiv:1703.05971 (2017).[64] B. Doyon, T. Yoshimura, and J.-S. Caux,arXiv:1704.05482 (2017).[65] B. Doyon, J. Dubail, R. Konik, and T. Yoshimura,arXiv:1704.04151 (2017).[66] V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E.Moore, arXiv:1704.03466 (2017).[67] V. B. Bulchandani, R. Vasseur, C. Karrasch, and J. E.Moore, arXiv:1702.06146 (2017).[68] M. Fagotti, J. Phys. A: Math. Theor. , 130402(2016).[70] A. Bastianello and A. De Luca, arXiv:1705.09270(2017).[71] A. L. de Paula, H. Bragana, R. G. Pereira, R. C. Dru-mond, and M. C. O. Aguiar, Phys. Rev. B , 45125(2017).[72] R. Vasseur, C. Karrasch, and J. E. Moore, Phys. Rev.Lett. , 267201 (2015).[73] H. Castella, X. Zotos, and P. Prelovˇsek, Phys. Rev. Lett. , 972 (1995);X. Zotos, F. Naef, and P. Prelovˇsek, Phys. Rev. B ,11029 (1997).[74] F. Heidrich-Meisner, A. Honecker, D. C. Cabra, and W.Brenig, Phys. Rev. B , 134436 (2003).[75] T. Prosen and M. ˇZnidariˇc, J. Stat. Mech. (2009)P02035.[76] J. Sirker, R. G. Pereira, and I. Affleck, Phys. Rev. Lett. , 216602 (2009). [77] T. Prosen, Phys. Rev. Lett. , 217206 (2011).[78] S. Langer, M. Heyl, I. P. McCulloch, and F. Heidrich-Meisner, Phys. Rev. B , 205115 (2011).[79] M. ˇZnidariˇc, Phys. Rev. Lett. , 220601 (2011).[80] R. Steinigeweg and W. Brenig, Phys. Rev. Lett. ,250602 (2011).[81] C. Karrasch, J. H. Bardarson, and J. E. Moore, Phys.Rev. Lett. , 227206 (2012).[82] T. Prosen and E. Ilievski, Phys. Rev. Lett. , 57203(2013).[83] C. Karrasch, J. E. Moore, and F. Heidrich-Meisner,Phys. Rev. B , 75139 (2014).[84] B. Doyon, Nuclear Physics B , 190 (2015).[85] R. Steinigeweg, J. Herbrych, X. Zotos, and W. Brenig,Phys. Rev. Lett. , 17202 (2016).[86] C. Karrasch, Phys. Rev. B , 115148 (2017).[87] E. Ilievski and J. De Nardis, Phys. Rev. Lett. ,020602 (2017).[88] B. Doyon and H. Spohn, arXiv:1705.08141 (2017).[89] V. Alba and P. Calabrese, PNAS , 7947 (2017).[90] E. Ilievski, M. Medenjak, T. Prosen and L. Zadnik, J.Stat. Mech. (2016) 064008.[91] E. Ilievski, J. De Nardis, B. Wouters, J.-S. Caux, F. H.L. Essler, and T. Prosen, Phys. Rev. Lett. , 157201(2015).[92] E. Ilievski, E. Quinn, J. De Nardis, M. Brockmann, J.Stat. Mech. (2016) 063101.[93] L. Piroli, E. Vernier, and P. Calabrese, Phys. Rev. B , 54313 (2016). [94] L. Zadnik, M. Medenjak, and T. Prosen, NuclearPhysics B , 339 (2016).[95] A. De Luca, M. Collura, and J. De Nardis, Phys. Rev.B , 020403 (2017).[96] H. Bethe, Z. Phys. , 205 (1931); R. Orbach, Phys.Rev. , 309 (1958).[97] L. Bonnes, F.H.L. Essler and A. M. L¨auchli, Phys. Rev.Lett. , 187203 (2014).[98] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. , 251 (1972).[99] K. Fabricius and B. M. McCoy, Phys. Rev. B , 16117 (2017).[101] M. Medenjak, C. Karrasch, T. Prosen, arXiv:1702.04677(2017).[102] M. Mesty´an and B. Pozsgay, J. Stat. Mech. (2014)P09020;B. Pozsgay, J. Phys. A: Math. Theor. , 74006 (2017).[103] L. Vanderstraeten, J. Haegeman, T. J. Osborne, and F.Verstraete, Phys. Rev. Lett. , 257202, (2014).[104] F. H. L. Essler, H. Frahm, F. G¨ohmann, A. Kl¨umperand V. E. Korepin, The One-Dimensional HubbardModel , Cambridge University Press (CUP), (2005).[105] A. Cavagli`a, M. Cornagliotto, M. Mattelliano and R.Tateo, J. High. Energy. Phys. 6, (2015).[106] V. Gritsev, E. Demler, M. Lukin, and A. Polkovnikov,Phys. Rev. Lett.99