Ballistic front dynamics after joining two semi-infinite quantum Ising chains
BBallistic front dynamics after joining two semi-infinite quantum Ising chains G abriele P erfetto , A ndrea G ambassi DISAT, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy SISSA - International School for Advanced Studies, via Bonomea 265, 34136 Trieste, Italy INFN, Sezione di Trieste, Trieste, ItalyE-mail: [email protected], [email protected]
We consider two semi-infinite quantum Ising chains initially at thermal equilibrium at two dif-ferent temperatures and subsequently joined by an interaction between their end points. Transportproperties such as the heat current are determined by the dynamics of the left- and right-movingfermionic quasi-particles which characterize the ensuing unitary dynamics. Within the so-calledsemi-classical space-time scaling limit we extend known results by determining the full space andtime dependence of the density and current of energy and of fermionic quasi-particles. Upon ap-proaching the edge of the propagating front, these quantities as well as the two-point correlationfunction display qualitatively different behaviors depending on the transverse field of the chain be-ing critical or not. While in the latter case corrections to the leading behavior are described, asexpected, by the Airy kernel, in the former a novel scaling form emerges with universal features.
I. I ntroduction
The non-equilibrium dynamics of isolated quan-tum statistical systems has recently captured a lotof attention from both the theoretical and the ex-perimental point of view. In fact, significant exper-imental advances have made it possible to finely con-trol trapped ultra-cold atomic gases (see, for instance,Refs. [1–5]). These systems are so weakly coupledwith the surrounding environment that they allowthe observation of a unitary non-equilibrium timeevolution, with the consequent remarkable phenom-ena which cannot be observed in standard condensedmatter systems due to decoherence and dissipativetransport.Within this context, one-dimensional systemsreaching a non-equilibrium steady state and support-ing a current, e.g., of energy, particles and charge,are now topical. In fact, on the one hand, they pro-vide an approximation of actual three-dimensionalsystems with strong anisotropy and on the other,they display an anomalous heat conduction whichviolates Fourier’s law because of the ballistic trans-port of energy. The role of spatial dimensionalityin determining the features of the quantum dynam-ics of many-body systems out of equilibrium hasbeen also experimentally demonstrated: for example,one-dimensional systems may relax towards a non-canonical distribution due to the possible presence ofadditional conservation laws which make them inte- grable compared to systems in higher dimensional-ity, see, e.g., Refs. [6–8]. In this work, we investi-gate within this framework the non-equilibrium dy-namics and transport properties of perhaps the bestknown integrable lattice model, i.e., the transversefield Ising chain (TFIC). In order to realize a non-equilibrium steady state (NESS) we adopt the pro-tocol which involves two Hamiltonian reservoir, see,e.g., Refs. [9, 10]: the system consists of two adja-cent TFICs, referred to as the left and the right chain,of finite length and initially disconnected and ther-malized at two different temperatures β − l and β − r ,respectively. Apart from this different initial condi-tions, the two chains are otherwise identical. The ini-tial probability distribution is encoded by the densityoperator ρ = e − β r H r ⊗ e − β l H l / Z , (1)where H r and H l are the Hamiltonians of the rightand left chain, respectively, while Z is a normaliza-tion constant. At time t = a r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r cal quantum quench, see, e.g., Refs. [11–13].The protocol described above has been extensivelystudied in various contexts: in particular, it hasbeen shown that for critical one-dimensional quan-tum systems a NESS with a factorized density ma-trix emerges in terms of right- and left-moving exci-tations [9, 14–16]. Similarly, the persistence of a non-equilibrium current at long times has been explainedin terms of the rightward and leftwards ballistic prop-agation of the excitations of the initial state. Theseexcitations enter into the adjacent chain from the con-tact point and they establish a non-equilibrium statewithin a spatial region, the extent of which is de-termined by the velocity of the propagating front[15, 16]. The resulting energy current and the cu-mulant generating function of the energy transferredalong the chain has been also determined in the long-time NESS [9, 14].For models of free fermions, instead, the completedynamics of the two-point correlation function and ofsome transport properties such as the energy density,the fermion concentration and the transverse mag-netization has been obtained also in the transientregime preceding the NESS, both on a lattice witheither an initial domain-wall state [17, 18] or a moregeneral factorized Fermi sea state [19] and on the con-tinuum [20, 21]. For the specific case of the TFIC,the dynamics of the magnetization has been derivedin Ref. [22] starting from a domain-wall initial state.In Ref. [23], instead, it has been shown that, startingfrom ρ in Eq. (1), a NESS with a factorized densitymatrix develops also if the TFIC is not critical and thestatistical properties of the stationary current (includ-ing the large-deviaton function) have bee n calculatedanalytically.Here we consider the TFIC, not necessarily at itscritical point, with an initial factorized thermal statedescribed by ρ in Eq. (1) and we extend the afore-mentioned results by describing the complete space-time evolution of various relevant quantities such asthe heat current, studying in detail how the eventualand known NESS is reached. Concerning the energycurrent, for example, we show that its non-analyticapproach to the propagating front of the excitationsdepends qualitatively on whether the transverse field h is critical ( h = h c =
1) or not. For the two-point cor-relation function we further investigate the behaviourupon approaching the edge of the front, showingthat, due to the initial finite temperatures of the ini- tially separate chains, these correlations acquire a cor-rection compared to the case at zero temperature,known to be described by the Airy kernel [24]. Forthe energy current, as long as h (cid:54) = h c , the edge be-haviour turns out to be well described by the Airykernel which determines a staircase structure of theprofile beyond the semi-classical approach describ-ing the dynamics far from the propagating front. For h = h c , instead, such a profile changes qualitativelyand in particular the staircase structure which char-acterizes the aforementioned Airy kernel is smoothedout and lost.The rest of the presentation is organized as follows:In Sec. II and Appendix A we briefly recall the ex-act solution of the TFIC, following Ref. [23], in orderto set the stage for studying the non-equilibrium dy-namics. In Sec. III, we determine, within the space-time scaling limit, the time evolution of the relevanttransport quantities, such as the energy current andthe related energy density. The details of the corre-sponding calculations are reported in Appendix B.In Sec. IV, we study the two-point correlation func-tion and the energy current close to the edge of thepropagating front, while the details of this analysisare collected in Appendix C. Section V summarizesour results and presents our conclusions. Part of thiswork is based on the unpublished results of Ref. [25]. II. T he I sing chain in a transversefield : exact solution In order to study the protocol discussed in the pre-vious section, we assume that the two TFIC of length N are originally disconnected, with the right one (de-noted by the subscript r ) occupying the lattice sites {
1, 2, . . . , N } along a line, while the left one ( l ) thesites {− N + − N +
2, . . . , 0 } . Accordingly, the cor-responding Hamiltonians are, respectively, H r = − J (cid:34) N − ∑ q = σ xq σ xq + + h N ∑ q = σ zq (cid:35) , (2a) H l = − J (cid:34) N − ∑ q = σ x − q σ x − q + + h N − ∑ q = σ z − q (cid:35) , (2b)where σ x , y , zq are the usual Pauli matrices, J is thecoupling strength and h the transverse field. Openboundary conditions are assumed for both chains.The pre-quench Hamiltonian H = H l + H r consistsof the two disconnected and independent chains.It is well known (see, e.g., Ref. [26]) that the firststep in order to diagonalize this model consists of theJordan-Wigner transformation: c q = (cid:18) e i π ∑ q − l = σ + l σ − l (cid:19) σ + q = (cid:32) q − ∏ l = σ zl (cid:33) σ + q , (3)where we have introduced the usual spin raising andlowering operators σ ± q = ( σ xq ± i σ yq ) /2. In terms ofthese new fermionic operators c q , with { c q , c † q (cid:48) } = δ q , q (cid:48) , the Hamiltonians H r , l acquire the bilinear form H r = − J N − ∑ q = (cid:16) c † q c † q + + c † q c q + + h . c . (cid:17) + Jh N ∑ q = c † q c q ,(4)where h . c . denotes the Hermitian conjugate of thepreceding expression and an analogous form holdsfor H l . As discussed in Ref. [27] and detailed in Ap-pendix A, Hamiltonians of this type can be diagonal-ized via a Bogoliubov transformation which suitablyintroduces two fields φ r , l ( k ) as φ r ( k ) = N ∑ q = (cid:104) ω qr ( k ) c q + ξ qr ( k ) c † q (cid:105) , (5)with an analogous expression for φ l ( k ) , but with co-efficients { ω ql ( k ) , ξ ql ( k ) } and q = − N +
1, . . . , 0. Interms of these fields one finds H r , l = ∑ k ε ( k ) φ † r , l ( k ) φ r , l ( k ) , (6)where the single-particle energy spectrum is given by ε ( k ) = J (cid:112) h − h cos k +
1. (7)Due to the finite length N of both chains, the set ofallowed values of k in the sum of Eq. (6) is discreteand, as a consequence of the open boundary condi-tions, determined by the implicit condition k n = n π N + + f ( k n ) N + n =
0, 1, ... N , (8)where f ( k ) ≡ arctan (cid:18) sin k cos k − h (cid:19) . (9) In the thermodynamic limit N → ∞ , both chains be-come semi-infinite, either to their right or to theirleft. Correspondingly, the set of allowed values k n becomes continuous within the interval [ π ] and,upon redefining Φ r , l ( k ) = lim N → ∞ ( N / π ) φ r , l ( k ) ,the Hamiltonians take the diagonal form H r , l = (cid:90) π dk ε ( k ) Φ † r , l ( k ) Φ r , l ( k ) , (10)where the single-particle energy spectrum ε ( k ) is thesame as in Eq. (7).At time t = H = H + δ H = H − J σ x σ x (11)where δ H represents the energy cost of connectingthe two chains through their closest end points at q = q =
1. Note that this operator is local, as ithas support only across the connection between thesetwo points. After the quench, the chain becomestranslationally invariant in the thermodynamic limitand thus [ H , P tr ] =
0, where P tr is the translation op-erator along the chain defined by the action: σ α q − = P † tr σ α q P tr , with α = x , y , z . (12)Since H is also invariant under the spatial inversion P , i.e., [ H , P ] =
0, one realizes that for each value ofthe wavevector k a two-fold degeneracy of the energyspectrum arises. Accordingly, one can introduce twofermionic operators Ψ R , L ( k ) (see Appendix A) whichare obtained via a suitable linear combinations of thepre-quench operators φ r ( k ) and φ l ( k ) (equivalently, Φ r ( k ) and Φ l ( k ) ) and which acquire opposite phasesunder the action of the translation operator: P † tr Ψ R , L ( k ) P tr = e ∓ ik Ψ R , L ( k ) ; (13)in terms of these Ψ R , L , the post-quench Hamiltonianbecomes: H = (cid:90) π dk ε ( k ) (cid:104) Ψ † R ( k ) Ψ R ( k ) + Ψ † L ( k ) Ψ L ( k ) (cid:105) ≡ H R + H L . (14) III. D ynamics in the semi - classicallimit As described in Sec. I, the NESS is obtained by join-ing at time t = ρ reported in Eq. (1). In order to be able to access thestationary state, both the time t and the system size N must be large, the latter being always larger than themaximal distance v max t travelled by the fermionic ex-citations at time t , where v max is their maximal veloc-ity (this quantity will be discussed further below inEq. (31)). In this case, the steady state density matrix ρ stat is formally defined by requiring that [14, 28, 29] (cid:104) O (cid:105) stat ≡ lim t → ∞ lim N → ∞ Tr [ O ρ ( t )] = Tr [ O ρ stat ] , (15)for any generic local observable O , as discussed fur-ther below. Accordingly, ρ stat can be formally ex-pressed as ρ stat = S ρ S † , (16)where the operator S = lim t → ∞ lim N → ∞ e − iHt e iH t (17)evolves states to time t = − ∞ according to the dy-namics of H and then brings them back to t = H . In order to ob-serve the stationary behaviour, measurements have tobe performed within the spatial region which has al-ready reached a stationary state, the typical extensionof which is given by v max t , since excitations propa-gate ballistically, as we shall see below. As a con-sequence, the spatial support of the observable O should include points at a maximal distance (cid:96) fromthe junction between the two chains which is muchsmaller than the distance v max t within which thesteady state is established at time t , i.e., (cid:96) (cid:28) vt (cid:28) N .Under these conditions, Eq. (16) defines the steadystate ρ stat which describes the steady average of anyoperator O with a finite support.The protocol described above for realising a sta-tionary state after joining two thermalised chainswhich act as asymptotic thermal baths, i.e., as Hamil-tonian reservoirs, is usually referred to as partition-ing protocol and it has been extensively studied[9, 16, 30]: in particular, for this type of quench, itis well-known that the stationary density matrix ρ stat eventually takes the form ρ stat = e − β r H L ⊗ e − β l H R / Z , (18)with H L , R given in Eq. (14). Essentially, this tellsus that the right- (viz. left-)moving excitations of the Hamiltonian are characterised by the initial tempera-ture of the left (viz. right) chain. A. Energy current and density
In the present work we will mainly focus on thedynamics of the energy current which emerges afterjoining the two chains. In particular, in order to de-fine the energy current flowing from the left to theright of the contact point x = H r of the right chain j rE ( t ) = dH r dt ; (19)however, one can define also the energy current as theopposite of the energy variation of the left chain, i.e.,as j lE ( t ) = − dH l / dt . Here we consider the sym-metric combination of these two equivalent contribu-tions, that as well quantifies the total energy trans-ferred from the left to the right chain across theirjunction link at time tj E ( t ) = j rE ( t ) + j lE ( t ) = dE ( t ) dt , (20)where we have introduced the energy difference E ( t ) ≡ ( H r − H l ) /2 between the left and right chain.This definition can be naturally extended to a genericpoint x of the chain with j E ( x , t ) = dE x ( t ) / dt beingthe energy current across the link between sites atposition x and x + j E takes the form j E ( x , t ) = ihJ e iHt (cid:16) c † x + c x − c † x c x + (cid:17) e − iHt . (21)Accordingly, the average energy current J E ( x , t ) attime t and point x along the chain is given by J E ( x , t ) = Tr [ j E ( x , t ) ρ ] . (22)The stationary and space-independent average value J NESSE of the current operator J E ( x , t ) in the NESSspecified by ρ stat in Eq. (18) was calculated in Ref. [23]according to the prescription of Eq. (15) and it turnsout to depend on the initial inverse temperatures as[31], J NESSE = g ( β l ) − g ( β r ) , (23)where the function g will be discussed further below,after Eq. (37). However, the approach of J E ( x , t ) to J NESSE was not previously investigated and here wefill in this gap. Equation (23) can be specialised tothe case of chains with a critical field h = h c = β r , l J (cid:29)
1, it turns outto agree with the general prediction J NESSE = π ( β − l − β − r ) /24 of conformal field theory [9, 14, 32–34].In order to determine the dynamics of the currentoperator in Eq. (21) and other similar observables,the approach described above — based on Eq. (18)— is not viable, as it provides information only onthe NESS. Accordingly, we directly calculate the av-erage of space-time dependent observables from theinitial density operator: in particular, the energy cur-rent J E ( x , t ) is determined from the trace in Eq. (22).In order to do this, one first expresses the operator j E in Eq. (21) in terms of right- and left-moving fermions Ψ R , L appearing in Eq. (14), by inverting the trans-formations reported in Sec. II which relate them tothe original fermionic operators c x , as detailed in Ap-pendix A. This is done via the following Bogoliubovtransformation: c x = (cid:90) π − π dk (cid:104) Ψ R ( k ) ( ω xR ( k )) ∗ + Ψ † R ( k ) ξ xR ( k ) (cid:105) , (24)where ω xR ( k ) and ξ xR ( k ) are given in Eq. (A.19). Re-membering that the dynamics of Ψ R ( k ) under thepost-quench Hamiltonian H in Eq. (14) is trivial, i.e., e − iHt Ψ R ( k ) e iHt = e i ε ( k ) t Ψ R ( k ) , the time evolution of j E in Eq. (21) can be explicitly determined (see Ap-pendix B). The remaining average over the initial den-sity matrix turns out to be J E ( x , t ) = (cid:90) π − π dkdk (cid:48) π ihJ (cid:16) e − ik − e ik (cid:48) (cid:17) I ( k , k (cid:48) ) e i ϕ x , t ( k , k (cid:48) ) (25)where ϕ x , t ( k , k (cid:48) ) = [ ε ( k ) − ε ( k (cid:48) )] t + x ( k (cid:48) − k ) , (26)with ε ( k ) defined in Eq. (7), while I ( k , k (cid:48) ) = Tr (cid:2) ρ Ψ † R ( k ) Ψ R ( k (cid:48) ) (cid:3) (see Eq. (B.9)) encodes the infor-mation about the initial state and is given by Eq. (B.4).In particular, we consider the so-called space-timescaling limit, also referred to as semi-classical or hy-drodynamic approach [19–22, 33, 35], in which boththe time t and the coordinate x are assumed tobe much larger than the corresponding microscopicscales, respectively set by J − — the inverse of the energy of a single link of the chain — and v max J − — the typical velocity of the excitations, introducedfurther below in Eq. (31) — but such that the ratio x / t takes arbitrary finite values. In this limit, J E ( x , t ) in Eq. (25) is determined by the values of k and k (cid:48) within the integration domains at which the phase ϕ x , t ( k , k (cid:48) ) in the exponential is stationary and by thepossible singularities of I ( k , k (cid:48) ) . Since ϕ x , t ( k , k (cid:48) ) turnsout to be stationary for k = k (cid:48) , the integral is then de-termined by the behaviour of I ( k , k (cid:48) ) for k (cid:39) k (cid:48) . Ac-cordingly, it is convenient to introduce the variables Q = k − k (cid:48) and K = ( k + k (cid:48) ) /2 and to consider the in-tegrand in Eq. (25) for Q (cid:39) ϕ x , t upto first order in Q and K one eventually gets: J E ( x , t ) = (cid:90) π − π dk π ε ( k ) v g ( k ) × (cid:2) f β l ( k ) Θ ( v g ( k ) t − x ) − f β r ( k ) Θ ( v g ( k ) t + x ) (cid:3) ,(27)where v g ( k ) = d ε ( k ) / dk is the group velocity ofthe relevant excitations, with ε ( k ) in Eq. (7), and f β ( k ) = [ + e βε ( k ) ] is the usual Fermi-Dirac dis-tribution at inverse temperature β which encodes thedistributions of quasi-particles in the chains beforethey are joined. In Eq. (27), Θ ( v ) indicates the stepfunction with Θ ( v ≥ ) = Θ ( v < ) =
0. Thisequation gives the exact profile of the average energycurrent at a certain time t and point x along the chainwithin the space-time scaling regime. Exploiting thecontinuity equation ∂∂ t u ( x , t ) = − ∂∂ x J E ( x , t ) , (28)expected to hold for the heat current, we can calcu-late the energy density u ( x , t ) along the chain, witha temporal integration of the r.h.s. of this equation.The initial condition u ( x , t = ) required in the inte-gration derives from the fact that at t = ε ( k ) are distributed accordingto the distribution f β l ( k ) on the left chain x ≤ f β r ( k ) on the right one x >
0, i.e., u ( x , 0 ) = (cid:90) π − π dk π ε ( k )[ f β l ( k ) Θ ( − x ) + f β r ( k ) Θ ( x )] .(29)After integration in time of Eq. (28) and the restrictionof the domain of integration, one eventually finds u ( x , t ) = (cid:90) π dk π ε ( k ) (cid:8) f β l ( k ) + f β r ( k )+ (cid:2) f β l ( k ) − f β r ( k ) (cid:3) (cid:2) Θ ( v g ( k ) t − x ) − Θ ( v g ( k ) t + x ) (cid:3) (cid:9) .(30)Expressions similar to Eqs. (27) and (30) have beenobtained with different approaches for other spinmodels in one spatial dimension: specifically, for theXX model (equivalent to a model of free fermions)it has been found [18] that the magnetization andmagnetization current evolve in the space-time scal-ing limit according to a scaling function which, apartfrom constants and prefactors, is similar to Eq. (39)and (35), respectively, discussed further below. Forthe TFIC in an initial thermal tensor state (of theform in Eq. (1)) for the two halves of the chain, in-stead, the space-time dependence of the transversemagnetization m ( q , t ) = (cid:104) σ zq ( t ) (cid:105) has been numericallystudied [36] for h = h c =
1, with one vanishing andone infinite initial temperature. The ensuing wave-front exhibits a light cone analogous to the one an-alyzed in this work, while the interpolation betweenthe asymptotic values of m ( q , t ) for | q | > v max t oc-curs linearly as a function of q / t . Since the space andtime dependence of the current in Eq. (27) is fully en-coded by the two functions Θ , their arguments canbe rescaled by a factor t > J E ( x , t ) turns out to be a function solely of the scaling vari-able v ≡ x / t , with v → ± ∞ corresponding to bothshort times and large distances and v → J E ( x , t ) ≡ J E ( v ) and the associated energy density u ( x , t ) ≡ u ( v ) for two different values of the inversetemperature β r with fixed β l = J − ). The profile of the heat current turnsout to be an even function of v , i.e., J E ( v ) = J E ( − v ) ,as it is clearly shown by the figure and by a careful in-spection of Eq. (27), see the discussion after Eq. (B.16).In addition, Fig. 1 shows that, because of causalityand the finite maximum value v max of the propaga-tion velocity v g ( k ) of the excitations, there is alwaysa region in space within which the initial state is notperturbed and, correspondingly, J E vanishes. Theseregions act as unperturbed "thermal reservoirs" forthe central part of the chain and, with the dispersion β r - - - J E ( v ) (a) β r - - - u ( v ) (b)FIG. 1. Dependence of (a) the heat current J E and (b) theenergy density u at time t and point x along the Ising chainon the scaling variable v ≡ x / t , within the space-time scal-ing limit. The parameters of the chain are h = J =
1, with the left part v < β l = v > β r = relation in Eq. (7), one finds (see also Ref. [12]) v max = max k ∈ [ − π , π ] | v g ( k ) | = J min ( h , 1 ) . (31)Accordingly, the "front" of the perturbation due tojoining the two chains propagates with velocity v max and the central perturbed region expands with veloc-ity 2 v max . Note that v max is at most J and with J = h > v max =
1, as it is clearlyshown by the plots in Fig. 1.Close to the joining point x = J E ( v = ) corresponding to the one J NESSE ( β l , β r ) eventually attained in the stationary state reacheduniformly over the entire chain; indeed, the value J E ( v = ) resulting from Eq. (27) coincides with J NESSE ( β l , β r ) reported in Eq. (23) and determined inRef. [23]. Equation (27) extends this known result onthe steady-state value of the current as to describe itstransient behaviour in both space and time, i.e., it pre-dicts how the far region of the thermal reservoirs andthe steady state value are asymptotically approachedby the dynamics.The physical interpretation of Eq. (27) is straight-forward in terms of the propagation of quasi-particleswith momentum k and velocity ± v g ( k ) which areproduced in the initial state with a statistics f β r ( k ) and f β l ( k ) for x > x <
0, respectively, andwhich contribute with ε ( k ) v g ( k ) dk to the flow of en-ergy. In particular, this interpretation has fist beenproposed for the TFIC in Ref. [37] and then used inRef. [22] for the same model and in Ref. [38] for theXXZ chain. In fact, since the post-quench Hamilto-nian H in Eq. (14) is diagonal in terms of the oper-ator Ψ R ( k ) , the states | k (cid:105) R = Ψ † R ( k ) | (cid:105) of the Fockspace (where | (cid:105) indicates the ground state of thechain) have an infinite lifetime and therefore prop-agate freely, with no scattering. Based on this pic-ture, Eq. (27) (as well as all the analogous equa-tions which are presented further below) could havebeen derived without the explicit calculations re-ported above. In fact, consider the space-time dia-gram in Fig. 2: the excitations with wavevector k > t = ± v g ( k ) for t > f β r ( k ) (bluerays in Fig. 2) [viz. f β l ( k ) (red rays)] originating from x > x <
0] also propagate into the comple-mentary part of the chain. As a result, the flux ofenergy (i.e., the energy current) produced by eachof these modes at a point with coordinate x (e.g., x = | v g ( k ) | t < | x | because theflux of energy carried by the particles with wavevec-tor + k cancels out the one of particles with wavevec-tor − k moving in the opposite direction and havingthe same statistics. This cancellation no longer oc-curs for | v g ( k ) | t > | x | because, for x >
0, the statis-tics of the excitations with veloc ity − v g ( k ) crossingthe world line of the point x (in green in Fig. 2) isgiven by f β r ( k ) while that of the excitations with ve-locity + v g ( k ) by f β l ( k ) , as they were originally gen-erated in the left part of the chain, see the sketchin Fig. 2. As a consequence, for each value of k ∈ - - x t FIG. 2. Space-time diagram in which the coordinate x alongthe chain and the time t > t = ± v g ( k ) , with k >
0, energy ε ( k ) and statistics f β l ( k ) for x < f β r ( k ) for x > t > x = [ π ] , the contribution to the energy flux is given by ε ( k ) v g ( k ) dk × [ f β l ( k ) − f β r ( k )] Θ ( v g ( k ) t − x ) for x > ε ( k ) v g ( k ) dk × [ f β l ( k ) − f β r ( k )] Θ ( v g ( k ) t + x ) for x <
0, which is equivalent to the integrand of Eq. (27).Analogous interpretation can be given to the expres-sion for the energy density u ( x , t ) reported in Eq. (30),which can be actually derived without the analysispresented above in this and in the previous section.The integral over k in Eq. (27) can be calculated inanalytic form, as detailed in Appendix B, and there-fore J E ( v ) can be written in the form J E ( v ) = Θ ( v max − | v | ) [ J ( β l , v ) − J ( β r , v )] (32)where J ( β , v ) = πβ { G ( β [ ε > ( v ) − ε < ( v )])+ − G ( β [ ε > ( v ) + ε < ( v )]) } , (33)with G ( x ) = − Li ( − e − x ) + x log ( + e − x ) , (34) ε > ( v ) = (cid:112) [ J max ( h )] − v , and ε < ( v ) has the sameexpression as ε > but with max replaced by min, suchthat (see Eq. (31)) ε < ( v ) = (cid:112) v max − v . Since ε ≷ ( v ) and therefore J E ( v ) depends only on v , the energycurrent in Eq. (32) is confirmed to be an even functionof v , as anticipated above.Figure 1 clearly shows that J E ( v ) , upon approach-ing the values ± v max of the variable v which corre-spond to the edge of the propagating front, displaysa non-analytic behaviour, which can be determinedon the basis of Eqs. (32) and (33). In particular, for v → ± v ∓ max one finds, at the leading order, J E ( v ) = C (cid:16) v max − v (cid:17) + O (( v max − | v | ) ) ,(35)where C is given in Eq. (B.27) [see also Eq. (B.29)]and depends on h and β r , l . Note that J E ( v ) vanishesat the edge according to a semi-circular law, as shownin Fig. 3(a), and consistently with what is observedin Ref. [18] for the XX chain evolving from a domain-wall initial state and in Ref. [22] for the TFIC in aninitial domain-wall state created by the action of a lo-cal Jordan-Wigner fermion operator. However, whenthe transverse field h of the Ising chain is poised atits critical value h c =
1, the constant C in Eq. (35)vanishes and the approach of J E ( v ) to the edge turnsout to change qualitatively (see the discussion afterEq. (B.27)), with J E ( v ) = β r − β l π (cid:16) v max − v (cid:17) + O (( v max − | v | ) ) .(36)This is clearly shown in Fig. 3(b), where we plot thebehaviour of J E ( v → v − max ) for the same conditions asin panel (a), but with h = J NESSE of the current within the space-time scal-ing limit corresponds to J E ( v = ) : from Eqs. (32)and (33) it takes the form J NESSE = G ( β l J | h − | ) − G ( β l J ( h + )) πβ l − G ( β r J | h − | ) − G ( β r J ( h + )) πβ r , (37)which reproduces the known expression reported inEq. (24) of Ref. [23] (which assumes h > G ( x ) here equals − j ( x ) therein.In order to highlight the qualitative differences inthe normalized profile J E ( v ) / J NESSE (with J NESSE = J E ( ) ) as a function of v / v max upon varying h , we plot vv max J E (a) vv max J E (b)FIG. 3. Energy current J E ( v ) as a function of the scalingvariable v / v max for a TFIC with (a) h = h = h c =
1, prepared in an initial state with β r = β l =
2, where J = v → v − max = Jh in Eqs. (35) and (36) for panels (a) and (b), respectively. it in Fig. 4 for three different values of the magneticfiled h . Within the central part of the interval, the sta-tionary state has already been reached and in fact thecurve approaches one. Near the edges, instead, thebehaviour of the normalized current changes signif-icantly at the critical point h = h c =
1, as the curveapproaches these edges with a vanishing slope, dif-ferently from the non-critical case where the semicir-cular behaviour of Eq. (35) causes its divergence.Exploiting the continuity equation (28), rewrittenin terms of the scaling variables ddv u ( v ) = v ddv J E ( v ) , (38)Eq. (35) can be used in order to derive the edge be- h0.811.3 - - vv max J E ( v ) J E ( ) FIG. 4. Dependence of J E ( v ) / J E ( ) on v / v max for h = h = h = h c =
1, from top to bottom, for the TFICwith J = β l =
1, and β r = haviour of u ( v ) ; for h (cid:54) =
1, it turns out to be u ( v → ± v ∓ max ) = u ( ± ∞ ) ± C arccos (cid:18) | v | v max (cid:19) + O (( v max − | v | ) ) ,(39)where u ( ± ∞ ) is the value of the spatially constantinitial energy density of the chain for x > x <
0, determined by the initial temperatures β − r and β − l , respectively [the explicit expression of u ( x , t = ) is reported in Eq. (29)]. Equation (39) is com-pared in Fig. 5(a) with the actual energy density pro-file which can be obtained from Eq. (30) and it dis-plays the same qualitative features as the recent re-sult of Ref. [22] concerning the magnetization (cid:104) σ xn ( t ) (cid:105) of the TFIC with an initial domain-wall state withinthe same semi-classical approach.When the field h is tuned to its critical value h c = u ( v → ± v ∓ max ) = u ( ± ∞ ) ± β r − β l π v max (cid:34) arccos (cid:18) | v | v max (cid:19) − | v | v max (cid:115) − v v max (cid:35) ++ O (( v max − | v | ) ) (40)instead of Eq. (39). Equation (40) is comparedin Fig. 5(b) with the actual energy obtained from vv max u (a) vv max u (b)FIG. 5. Energy density u ( v ) as a function of the scalingvariable v / v max for the same values of the parameters asthose in the corresponding panels of Fig. 3. In both panels,the lower (blue) curves correspond to Eq. (30), while theupper (red) curves to their approximation close to the edge v → v − max = Jh in Eqs. (39) and (40) for panels (a) and (b),respectively. Eq. (30). Their qualitative behaviour upon approach-ing the edge v (cid:39) v max is markedly different from theone reported in panel (a) of the same figure for thenon-critical chain. B. Correlation function and density of quasi-particles
The procedure outlined above for predicting theevolution of the energy current and density in thespace-time scaling limit can be extended to other rel-evant quantities. Here we focus on the quasi-particleexcitations and we study their dynamics along thechain by determining the corresponding spatial den-sity and current. In particular, we consider the0Fourier transform of the operator Ψ R ( k ) , i.e., γ x = (cid:90) π − π dk √ π e ikx Ψ R ( k ) . (41)In terms of γ x , the total number operator ˆ N of quasi-particles excitations along the chain can be easily ex-pressed asˆ N = (cid:90) π − π dk Ψ † R ( k ) Ψ R ( k ) = + ∞ ∑ x = − ∞ γ † x γ x , (42)which is conserved by the post-quench dynamics dic-tated by Eq. (14) since the quasi-particles γ x prop-agate freely along the chain, without experiencingscattering. According to Eq. (42), γ † x γ x can be inter-preted as the density of quasi-particles on the lattice.In order to calculate the evolution of its expectationvalue, we first consider the two-point, equal-time cor-relation function (cid:104) γ † x ( t ) γ y ( t ) (cid:105) = Tr [ ρ γ † x ( t ) γ y ( t )] ofthe fermionic operators γ x and γ y at two distinct po-sition x and y along the chain. The calculation of thisquantity proceeds exactly as described above for theenergy current and density in the space-time scalinglimit and, for brevity, we do not reproduce it here butwe report only the final expression: (cid:104) γ † x ( t ) γ y ( t ) (cid:105) = (cid:90) π − π dk π e ik ( x − y ) Θ ( v g ( k ) t − x + y ) f β l ( k )+ (cid:90) π − π dk π e ik ( x − y ) Θ ( x + y − v g ( k ) t ) f β r ( k ) .(43)The density n ( x , t ) = (cid:104) γ † x ( t ) γ x ( t ) (cid:105) of these quasi-particle excitations can be readily obtained by setting y = x in the previous expression: n ( x , t ) = (cid:90) π − π dk π (cid:2) Θ ( v g ( k ) t − x ) f β l ( k )+ Θ ( x − v g ( k ) t ) f β r ( k ) (cid:3) , (44)which, in the stationary limit t → ∞ , agrees withEq. (5) of Ref. [39] for the mode occupation num-bers in the NESS of the TFIC generated as discussedhere. Equation (43) for t → ∞ , instead, has a sim-pler structure compared to the analogous expressionwhich was derived in Ref. [39] (see Eq. (7) therein)for the two-point correlation function. This is due tothe fact that Eqs. (43) and (44) refer to the fermionicoperator γ x in Eq. (41), while Ref. [39] provides thecorresponding expression for the correlation function of the operators c x , y ( t ) introduced with the Jordan-Wigner transformation in Eq. (3), which the γ x ’s arelinearly related to. Note that also the expression inEq. (44) for n ( x , t ) can be given a semi-classical inter-pretation, analogous to the one explained in Fig. 2.In fact, at t = f β r ( k ) to be uniformlydistributed in space on the right part of the chain,i.e., at x >
0, while those with distribution f β l ( k ) tobe on the left one at x <
0. (This initial condition iscorrectly reproduced by Eq. (44) upon setting t = t >
0, the particles with a certain k propagate bal-listically and independently with their characteristicvelocity v g ( k ) and therefore the ensemble of particleswith initial statistics f β l ( k ) reach all the points with x < v g ( k ) t , while those with statistics f β r ( k ) , all thepoints with x > v g ( k ) t . Translated into equations,this picture yields directly Eq. (44).Proceeding as done above for the energy density, aparticle current J N ( x , t ) can be associated with n ( x , t ) on the basis of a continuity equation in which thetransported quantity is now the number of quasi-particles. Taking into account the boundary condi-tion J N ( ± ∞ , t ) =
0, one eventually finds J N ( x , t ) = − (cid:90) x − ∞ dx (cid:48) ∂ n ( x (cid:48) , t ) ∂ t (45)and, by using Eq. (44), J N ( x , t ) = (cid:90) π − π dk π v g ( k ) (cid:2) f β l ( k ) − f β r ( k ) (cid:3) Θ ( v g ( k ) t − x ) .(46)This expression has exactly the same interpretationas Eq. (27) and, after having restricted the integra-tion domain to 0 ≤ k ≤ π in both expressions, it canbe shown to coincide with the latter [see Eq. (B.17)]upon replacing ε ( k ) (i.e., the energy) with 1, as ap-propriate for a quasi-particle density. The plots ofthe scaling functions J N ( x , t ) ≡ J N ( v = x / t ) and n ( x , t ) ≡ n ( v = x / t ) in Eq. (44) display the samequalitative features as the energy current and energydensity in Fig. 1, respectively, with a marked prop-agating front which moves ballistically. In fact, bylooking specifically at Fig. 1, one realizes that as timeelapses, the gradient of the particle number n ( x , t ) around the junction point x = t → + ∞ . Since the particlecurrent J N keeps anyhow a non-zero value, it cannotbe proportional to the gradient of n ( x , t ) , as requiredby diffusion, but to n ( x , t ) , as expected by ballistic1transport. Accordingly, every possible diffusive com-ponent of the current J N gets suppressed, making thesteady-state transport of fermions purely ballistic.As argued above and discussed in more detail inAppendix B, J E ( v ) and J N ( v ) have analogous formsand actually they can be obtained from the function J a ( v ) introduced in Eq. (B.18) by setting a = J N ( v ) has the same expres-sion as Eqs. (32) and (33) with J replaced by J J ( β , v ) = πβ { G ( β [ ε > ( v ) − ε < ( v )])+ − G ( β [ ε > ( v ) + ε < ( v )]) } , (47)and G by G ( x ) = ln (cid:0) + e − x (cid:1) . (48)In particular, also J N ( v ) shows a non-analytic be-haviour as v approaches ± v max , which has the sameform as Eq. (35), with C replaced by the constant C reported in Eq. (B.27). Moreover, similarly to J E ( v ) ,the qualitative features of the approach of J N to theedge do change if the transverse field h is tuned to itscritical value h c =
1, with J N ( v ) = ( β r − β l ) π ( v max − v ) + O (( v max − | v | ) ) .(49)In Fig. 6 we compare the actual approach of J N ( v ) to the edge v max with the approximation providedby the semicircular law for the non-critical case inpanel (a) and by Eq. (49) for the critical case in panel(b). From the continuity equation relating the par-ticle current J N ( v ) with the fermionic density n ( v ) ,exploiting the edge behaviour of Eq. (35) with C re-placed by C , one can determine the behaviour of n ( v ) upon approaching the edge, based on the pre-vious results. For h (cid:54) =
1, it turns out to be n ( v → ± v ∓ max ) = n ( ± ∞ ) ± C arccos (cid:18) | v | v max (cid:19) + O (( v max − | v | )) , (50)where n ( ± ∞ ) has, for the fermionic density, the samemeaning as u ( ± ∞ ) in Eq. (39) for the energy density.When h is tuned to its critical value h c =
1, the pre-vious equation must be corrected taking into accountEq. (49): n ( v → ± v ∓ max ) = n ( ± ∞ ) ∓ β r − β l π ( | v | − v max )+ O (( v max − | v | ) ) . (51) / v max J N (a) / v max J N (b)FIG. 6. Fermion current J N ( v ) as a function of the scalingvariable v / v max for a TFIC with (a) h = h = h c =
1, prepared in an initial state with β l = β r = J = v → v − max = Jh . In Fig. 7 the behavior of n ( v ) as v approaches theedge v max is compared with that of the expansionprovided by Eq. (50) for the non-critical case in panel(a), and by Eq. (51) for the critical one in panel (b).Analogously to the case of the energy current dis-cussed above, the value J NESSN of the particle current J N in the steady state can be exactly calculated by set-ting v = G (cid:55)→ G ,see Eq. (48); this result into Eq. (37) with G (cid:55)→ G ,i.e., J NESSN ( β l , β r ) = J ( β l , 0 ) − J ( β r , 0 ) , (52)where J ( β , 0 ) = πβ ln (cid:32) + e − β J | h − | + e − β J ( h + ) (cid:33) ; (53)2 / v max n (a) / v max n (b)FIG. 7. Fermion concentration n ( v ) as a function of thescaling variable v / v max for the same TFIC as in Fig. 6. Thelower (blue) curves correspond to the analytic expressiondiscussed in the main text for (a) a non-critical and (b) thecritical value of the transverse field h , while the upper (red)curves to the corresponding approximations close to theedge v → v − max = Jh . this expression, as well as the one of the energy cur-rent in Eq. (37), takes the general form of Eq. (23),which often appears in studies of transport in non-equilibrium quantum stationary states, see, e.g.,Refs. [20, 21, 23]. In the case of a chain with criti-cal transverse field h = h c =
1, in the limit β r , l J (cid:29) J NESSN = ln22 π ( β − l − β − r ) , (54)which coincides with the result obtained in Refs. [20,21] for a local quench of non-interacting Fermi gasesin one spatial dimension. IV. C orrelations and energy currentnear the edge
In this section we investigate in more detail the dy-namics of the propagating front near the edge, i.e., for x (cid:39) ± v max t : in fact, experience with analogous cases[17, 40] suggests that both the correlations and theenergy current acquire non-trivial corrections withina distance ∆ x ∝ t from the edge and we shall seebelow how they emerge in the present setting. Asthe relative width ∆ x / x ∼ t − of the spatial re-gion interested by these corrections vanishes in thespace-time scaling limit, these features are not cap-tured by the previous analysis and therefore they re-quire a separate treatment. A. Correlation functions
Consider the two-point correlation function (cid:104) γ † x ( t ) γ y ( t ) (cid:105) = (cid:90) π − π (cid:90) π − π dkdk (cid:48) π e i [ φ x , t ( k ) − φ y , t ( k (cid:48) )] I ( k , k (cid:48) )= (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l + (cid:104) γ † x ( t ) γ y ( t ) (cid:105) r , (55)where φ x , t ( k ) ≡ ε ( k ) t − kx and I ( k , k (cid:48) ) is definedin Eq. (B.9). Note that the r.h.s. of this equationnaturally decomposes as the linear superposition oftwo distinct contributions (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l , r correspond-ing to the effect of considering separately one of thetwo half chains initially populated according to thecorresponding thermal distribution, while the otherunoccupied (i.e., with f β r , β l = Q = k − k (cid:48) and K = ( k + k (cid:48) ) /2 and after expanding around Q = Q : this renders Eq. (43) discussedin the previous section. Here, instead, we are inter-ested in the behaviour of this quantity near the edgeof the propagating front, corresponding to having | x | , | y | (cid:39) v max t : in this case, higher-order correctionsin the expansion of the phases φ x , t ( k ) and φ y , t ( k (cid:48) ) around the respective stationary points become im-portant and therefore they have to be accounted for.Namely, as v ≡ x / t approaches ± v max , the two solu-tions k + s ( v ) and k − s ( v ) of the stationary phase equa-3 h = = π π kv g ( k ) v max k s k s + k s - FIG. 8. Graphical representation of the solution of Eq. (56).The dashed and solid curves represent the group velocity v g ( k ) as a function of k for h = (cid:39) h c and h = J =
1. The lowest (green) horizontaldashed line indicates a certain assigned value of | v = x / t | and, correspondingly, the values k ± s ( v ) of k at which itcrosses the previous curves are the solution of the equation.As | v | approaches v max , indicated by the upper dashedhorizontal line, the two solutions k + s and k − s merge intoa unique value k s . tion (see Appendix B for details) v g ( k ± s ( v )) = | v | (56)merge into a unique stationary point k s = k ± s ( v max ) obtained by taking v = v max into Eq. (B.22), atwhich the second derivative of the phases φ x , t ( k ) and φ y , t ( k (cid:48) ) vanishes, as shown in Fig. 8. Accordingly, oneexpects non-trivial corrections due to higher-orderterms in the expansion of φ x , t ( k ) and φ y , t ( k (cid:48) ) aroundthe stationary point k s . In particular, the third-ordercorrection is expected to provide the dominant con-tribution and, accordingly, the phase is approximatedas φ x , t ( k ) = ε ( k s ) t − k s x + ( k − k s )( v max t − x ) − v max ( k − k s ) t + O (( k − k s ) ) , (57)where we used the fact that ε (cid:48)(cid:48)(cid:48) ( k s ) = − v max and ε (cid:48)(cid:48) ( k s ) =
0, see Eq. (C.2). In order to evaluate the inte-gral in Eq. (55) within this approximation, it is conve-nient to introduce the variables K = ( v max t /2 ) ( k − k s ) and Q = ( v max t /2 ) ( k (cid:48) − k s ) instead of k and k (cid:48) : by expanding the Fermi-Dirac distributions f β l , r which, via I ( k , k (cid:48) ) [see Eq. (B.9)], appears on ther.h.s. of Eq. (55), one obtains an expansion of the "left" contribution (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l to the correlation functionin Eq. (55) (an analogous result holds for the "right"contribution (cid:104) γ † x ( t ) γ y ( t ) (cid:105) r ). Up to order t − , onefinds (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l e − ik s ( x − y ) f β l ( k s ) = (cid:18) v max t (cid:19) K A ( X, Y )+ (cid:18) v max t (cid:19) K β l ( X, Y ) + O (cid:18) v max t (cid:19) (58)where the relevant scaling variables are X ≡ x − v max t ( v max t /2 ) , and Y ≡ y − v max t ( v max t /2 ) . (59)The analogous expression for (cid:104) γ † x ( t ) γ y ( t ) (cid:105) r is ob-tained from the previous one upon replacing β l with β r . The leading contribution on the r.h.s. of Eq. (58)is the so-called Airy kernel, which is given by K A ( X , Y ) = Ai ( X ) Ai’ ( Y ) − Ai’ ( X ) Ai ( Y ) X − Y , (60)in terms of the Airy function Ai ( X ) (see Eq. (C.9)for its integral representation); as Eq. (58) dependson x , y , and t via the scaling variables in Eq. (59),it expresses the scaling behaviour of the correlationfunction, which emerges within a spatial region ofthickness ∝ t from the location ∝ t of the edge andtherefore it cannot be captured by the semi-classical,hydrodynamic, approach discussed in the previousSection, as anticipated above. In particular the Airykernel results from the fact that the two solutions ofEq. (56) coincide near the edge of the light-cone and,in fact, this kernel emerges rather generically in theliterature concerning free spinless fermionic chains[17, 41], where it has been reported for the case of aninitial state consisting of a fully occupied half chainand for a more general initial factorized Fermi seastate [19]. Equation (58) shows that the leading ef-fect of an initial state with two different (finite) tem-peratures β − l and β − r is the presence of the cor-responding distributions f β l ( k s ) [and f β r ( k s ) in theanalogous expression for (cid:104) γ † x ( t ) γ y ( t ) (cid:105) r which we donot report here] as a multiplicative factor of the Airyscaling function which emerges also for the differentinitial conditions mentioned above, corresponding to f β l ( k ) = f β r ( k ) = f β ( k ) = Θ ( k F − k ) for the Fermi sea state4where k F denotes the Fermi momentum. In passingwe mention that the Airy kernel and a generalizationof it emerge at the spatial edge of a system at zeroand finite temperature, respectively, also in the caseof a one-dimensional gas of free fermions confinedby an harmonic potential [42–44]. However, in thiscase, the edge does not expand in time but is ratherfixed by the presence of the harmonic potential whichmakes the fermion density vanish beyond a certaindistance from the center of the trap. Close to thatedge, the correlation function is expressed as a deter-minantal process whose kernel can be interpreted asan extension of the Airy one. In the present case oftwo Ising chains with an initial thermal distributionbut with two different temperatures β − l and β − r , theleading-order behaviour at the edge is modified in asomehow expected way, i.e., it has the same form asin the aforementioned cases with domain wall andfactorized Fermi sea except for the multiplicative fac-tor determined by the corresponding thermal distri-bution f β r , l ( k s ) evaluated at the saddle point k s . Inorder to highlight the effects of the first non-trivialcontributions due to the finite initial temperatures,we report in Eq. (58) also the rescaled first-order cor-rection ( v max t ) K β l [see its explicit expression inEq. (C.7)], which is compared with the leading order K A in Fig. 9. One can notice that, close to X =
0, i.e.,for | x | (cid:39) v max t , the first-order correction K β l turns outto be one order of magnitude smaller than the AiryKernel K A , for both Y = Y = t = Y increasesone gets progressively away from the edge region forthe Y variable and thus correlations are expected tobe captured in a less precise way by the expansion inEq. (58). B. Energy current
The procedure described above for studying thecorrelation functions close to the edge, can be also ap-plied to the energy current. In this case, at the leadingorder for x (cid:39) v max t (see Appendix C for details, withan analogous expression holding for x (cid:39) − v max t dueto the symmetry J E ( x , t ) = J E ( − x , t ) ), one finds J E ( X , t ) = ε ( k s ) v max [ n l ( X , t ) − n r ( X , t )] , (61) K A ( X,0.5 ) K l ( X,0.5 ) - - - - - K A ( X,1 ) K l ( X,1 ) - - - - - K A ( X,2 ) K l ( X,2 ) - - - - - FIG. 9. Airy Kernel K A ( X , Y ) (blue) and first-order cor-rection ( v max t ) K β l ( X , Y ) (red) as functions of X for Y = h = J =
1, and t = where the scaling variable X is given in Eq. (59), n l , r ( X , t ) = (cid:104) γ † x ( t ) γ x ( t ) (cid:105) l , r = ( v max t /2 ) − f β l , r ( k s ) K A ( X , X ) , (62)(see Eq. (55)) and the kernel K A ( X , X ) = (cid:2) Ai (cid:48) ( X ) (cid:3) − X [ Ai ( X )] (63)is obtained as the limit Y → X of Eq. (60). As itis clear from the structure of Eq. (61), the current J E at the leading order can be given a simple inter-pretation as resulting from the superposition of the5 K A ( X,X ) π - X - - - - - FIG. 10. Dependence of the energy current J E ( X , t ) on therescaled coordinate X [see Eq. (59)] near the edge X = t = h = J = β l =
2, and β r = energy current ε ( k s ) v max n l ( X , t ) close to the edge,due to the fastest excitations (at temperature β − l )propagating rightward and originally produced onthe left part of the chain and the one with oppo-site sign − ε ( k s ) v max n r ( X , t ) due to those originat-ing on the right part of the chain and propagatingleftward. In particular, the semi-classical limit dis-cussed in Sec. III A — in which v ≡ x / t ≤ v max is kept constant as t → ∞ — corresponds to hav-ing X ∝ ( v − v max ) t → − ∞ here and, in fact, inthis limit, Eq. (61) renders the behavior of J E ( v ) closeto the edge, reported in Eq. (35). This can be eas-ily seen by using Eq. (B.29) in Eqs. (61) and (62) andby taking into account the asymptotic behaviour ofthe Airy Kernel K A ( X , X ) → √− X / π , which followsfrom Eq. (C.11).Figure 10 presents a plot of the current J E ( X , t ) × ( v max t /2 ) as a function of the scaling variable X ,which is compared with Eq. (61), while the dashedline indicates the asymptotic behavior for X → − ∞ .The solid line features the typical staircase behaviourcaused by the cubic term in the expansion of thephase in Eq. (57) which is therefore not capturedby the semi-classical limit discussed in Sec. III A,which actually corresponds to the square-root enve-lope (dashed line) of the boundary scaling regime.These oscillations are similar to those obtained inRef. [17] for a free-fermionic chain starting from adomain-wall initial state, in which case the subse-quent steps in the staircase have been explained onthe basis of the correspondence existing between thecounting statistics of free fermions at the edge of apropagating front and that one of the eigenvalues of a random matrix. As we noted in Sec. III A (see, inparticular, Fig. 4), the behavior of the current J E in thespace-time scaling limit changes qualitatively whenthe transverse field h takes its critical value h c = k s → h → ε ( k s ) → J E ( X , t ) in Eq. (61) vanishidentically. In this case, within the stationary-phaseapproximation adopted here, one has to keep termsup to the first non-vanishing order ∝ t in the ex-pansion in k and k (cid:48) around k s . Proceeding in this way(see Appendix C for details), one finds J E ( X , t ) = t v max ( β r − β l ) K c ( X ) (64)where X in Eq. (61) (see Eq. (59)) is replaced by X ≡ x − v max t ( v max t /8 ) (65)and K c ( X ) = (cid:8) X [ Ai ( X )] −
12 Ai ( X ) Ai’ ( X ) − X [ Ai’ ( X )] (cid:9) . (66)By using the property of the Airy functionAi (cid:48)(cid:48) ( X ) = X Ai ( X ) , (67)one can easily show by differentiating the previousequation that K c ( X ) is related to K A ( X , X ) in Eq. (63)by (see also Appendix C) K c ( X ) = (cid:90) + ∞ X dY K A ( Y , Y ) . (68)As in the case h (cid:54) = h c discussed above, Eq. (64) ren-ders Eq. (36) after taking into account the asymptoticbehaviour K c ( X → − ∞ ) (cid:39) √− X / ( π ) , which canbe again obtained from Eq. (C.11). Figure 11 presentsa plot of the current J E ( X , t ) × t / [ Jv max ( β r − β l )] ,i.e., of K c ( X ) (solid line) as a function of the scalingvariable X , which is compared with the asymptoticbehavior for X → − ∞ (dashed line).One immediately notes that for X > K c ( X ) is qualitatively similar to the noncritical one K A ( X , X ) which determines, up to con-stants, J E in Eq. (61) and they both decay exponen-tially upon increasing X , as can be readily checked6 K C ( X ) π - X - - - - - FIG. 11. Dependence of the energy current J E ( X , t ) on therescaled coordinate X [see Eq. (65)] near the edge X = h = h c = from Eq. (C.10). On the contrary, for X <
0, the typ-ical staircase structure of the Airy kernel K A ( X , X ) shown in Fig. 10 is absent in the critical case K c ( X ) reported in Fig. 11, due to the fact that the integra-tion in Eq. (68) smoothens it to the extent that it is nolonger visible. Accordingly, when h is set to its crit-ical value h c =
1, although the relevant scaling vari-able is ∝ ( x − v max t ) / t as in the non-critical caseof Eq. (59) (this is due to the fact that Eq. (56) stilladmits a unique solution k s = x → v max t ), the qualitative features of J E ( X , t ) changesignificantly as a consequence of the fact that the en-ergy gap vanishes and a novel kernel K c ( X ) (relatedto the Airy kernel by Eq. (68)) emerges. Note that K c ( X ) is as "universal" as the Airy kernel, as it doesnot depend on the specific properties of the systemunder investigation, i.e., J , v max , etc. and, togetherwith the scaling function which involves it, is essen-tially determined by ε (cid:48) ( k s ) and ε ( ) ( k s ) . V. S ummary and perspectives
In this work we have investigated the non-equilibrium dynamics induced by joining via a localinteraction two independent transverse field quan-tum Ising chains, initially thermalized at two differ-ent temperatures. In Sec. II the model and the quenchprotocol have been introduced: in particular it hasbeen shown, following Ref. [23], how the restoredtranslational invariance of the final Hamiltonian nat- urally determines a two-fold degeneracy of the singleparticle spectrum ε ( k ) , which allows the descriptionof the excitations of the full chain in terms of right-and left-moving quasi-particles.In terms of these excitations, in Sec. III, the cur-rent of energy J E ( x , t ) , of fermionic quasi-particles J N ( x , t ) and the corresponding densities u ( x , t ) and n ( x , t ) , respectively, have been exactly calculated inthe space-time scaling limit in which, formally, both x and t are assumed to be large on the correspond-ing microscopic scale, with a fixed ratio v = x / t . Inparticular, one of the main result of the present workconcerns the form of the profile of the energy andthe quasi-particles currents and densities in Fig. 1 asfunctions of v and the fact that their qualitative be-havior depends on the transverse field h being criticalor not, as shown in Figs. 3, 4, 5, 6, and 7. All these fig-ures show how these quantities propagate along thechain in the form of a front travelling with the char-acteristic velocity given in Eq. (31). At any finite t ime t , the sites of the chain which are further away fromthe origin than v max t are unperturbed and still retaintheir initial features. These asymptotically far regionsplay the role of thermal baths which allow the devel-opment, and eventually the persistence in the steadystate, of a non-equilibrium dynamics. In this con-text, our results in Eqs. (27), — explicitly calculatedin Eqs. (32), (33), and (34) — (30), (44), and (46), for J E , u , J N , and n , respectively, generalize the well-knownpicture of current-carrying steady states [9, 14, 15] bydescribing the whole dynamics, both in space andtime, and also how the stationary state is actuallyapproached. In Sec. IV we investigate the behav-ior of the two-point correlation function in Eq. (55)and of the energy current J E ( x , t ) in Eq. (25), closeto the edge of the front, i.e., for x (cid:39) v max t and be-yond the space-time scaling limit discussed in Sec. IIIsummarized above. In particular we show that, asit occurs in similar cases investigated in the litera-ture [17, 19, 40, 41], these quantities acquire a "uni-versal" behavior within a region of width ∆ x ∝ t − around the edge at x (cid:39) v max t , conveniently expressedin terms of the scaling variable X ∝ ( x − v max t ) / t .Since the two Ising chains are initially in a thermalstate, the Fermi-Dirac statistics enters in the expres-sion of these two physical observables, respectivelyin Eq. (58) and Eq. (61). At the leading order, for longtimes on the microscopic scale, they coincide (up tomultiplicative constants) with the kernel which char-7acterizes the behaviour of the edge of either a systemof free fermions initially prepared either in a factor-ized Fermi-sea ground state with different fillings onthe two parts of the chain [19] or in a domain-wallinitial state [17, 41]. In order to investigate the ef-fects of temperature on correlations beyond this sim-ple rescaling, we also determined the first-order non-trivial correction for the two-point correlation func-tion, which can be expressed in terms of the Airykernel. As observed in Sec. III for the profiles of var-ious quantities as functions of v = x / t in the space-time scaling limit, a change in the qualitative behav-ior of the current J E occurs as a function of X as thetransverse field takes its critical value h c =
1; in fact,the leading term in Eq. (61) vanishes and thereforehigher-order terms must be considered. Specifically,a scaling form emerges which involves a novel ker-nel K c ( X ) , see Eq. (66), which is actually the integralof the Airy kernel and in which the staircase struc-ture characterising the latter basically disappears, seeFigs. 10 and 11.Among the possible extensions of the present workwe plan to carry out the analysis of the scaled cumu- lant generating function [9, 45] of the energy in thespace-time scaling limit, in order to study how thewhole statistics of this quantity and therefore its fluc-tuations changes upon approaching the edge of thepropagating front and upon tuning h to its criticalvalue. So far, only fluctuations in the steady statehave been studied [23]. Moreover, we will also ex-plore the possibility to extend the general diagonal-ization procedure explained in Sec. II to other one-dimensional quantum chains models, even not trans-lationally invariant like that one studied in Ref. [46].A description of these models in terms of right/leftmoving quasi-particles is still lacking and it is in factnot yet clear which form these excitations may takefor such systems and how this could influence trans-por t properties and their statistics. A cknowledgments The authors would like to thank V. Eisler, M. Kor-mos, A. De Luca and J. Viti for useful discussions.
Appendix A. Exact solution of the transverse field Ising chain
The Hamiltonian in Eq. (4) can be written as a quadratic form in the fermionic operators c i and c † i as H r = ∑ i , j (cid:20) c † i M ij c j + (cid:16) c † i N ij c † j + h . c . (cid:17)(cid:21) , (A.1)where M ij = − J ( δ i + j + δ i , j + ) − Jh δ ij and N ij = − J ( δ i , j + − δ i + j ) . (A.2)An analogous expression holds for H l . In order to diagonalize this type of Hamiltonians we perform the canonicaltransformation reported in Eq. (5) and then impose that the resulting Hamiltonian, once expressed in terms of φ r ( k ) , takes the diagonal form reported in Eq. (6). This amounts at requiring that [ φ r , H r ] − ε ( k ) φ r ( k ) =
0, (A.3)which, implies the following conditions for the coefficients ω qr ( k ) and ξ qr ( k ) in Eq. (5): ε ( k ) ω qr ( k ) = ∑ j (cid:16) ω jr ( k ) M jq − ξ jr ( k ) N jq (cid:17) , (A.4) ε ( k ) ξ qr ( k ) = ∑ j (cid:16) ω jr ( k ) N jq − ξ jr ( k ) M jq (cid:17) . (A.5)8At this point it is convenient to express ω qr ( k ) and ξ qr ( k ) in terms of their sum and difference, i.e., as ω qr ( k ) = A qr ( k ) + B qr ( k ) ξ qr ( k ) = A qr ( k ) − B qr ( k ) A r ( k )( M − N ) = ε ( k ) B r ( k ) , (A.7a) B r ( k )( M + N ) = ε ( k ) A r ( k ) , (A.7b)where A r = (cid:0) A r , A r , A r ... A Nr (cid:1) is the vector of the coefficients { A qr } q = N , and analogous for B r . Written interms of components, these equations become: (cid:40) hJ A qr ( k ) − J A q + r ( k ) = ε ( k ) B qr ( k ) , − JB qr ( k ) + hJB q + r ( k ) = ε ( k ) A qr ( k ) , (A.8)with q =
1, 2, . . . , N and boundary conditions A N + r ( k ) = B r ( k ) =
0. The solution of this system ofequation is [27] A qr ( k ) = a k sin ( qk − f ( k )) and B qr ( k ) = b k sin ( qk ) , (A.9)where a k and b k are (ortho-)normalization constants which are determined by requiring the normalization ofthe vectors, i.e., ∑ Nq = A qr ( k ) A qr ( k (cid:48) ) = ∑ Nq = B qr ( k ) B qr ( k (cid:48) ) = δ k , k (cid:48) , while k and f ( k ) are given by Eqs. (8) and (9),respectively.In the thermodynamic limit N → ∞ the set of allowed values of k becomes continuous and covers the interval [ π ] , while the functions in Eq. (A.9) become A qr ( k ) = (cid:114) π sin ( qk − f ( k )) and B qr ( k ) = (cid:114) π sin ( qk ) . (A.10)The procedure outlined above carries over to the left chain in Eq. (2); in particular one introduces φ l ( k ) = ∑ q = − N + (cid:16) ω ql ( k ) c q + ξ ql ( k ) c † q (cid:17) (A.11)and the matrices M ij and N ij analogous to Eq. (A.2), but appropriate for the right chain, as well as the vectors A l and B l as discussed above for the right chain. Due to the relationship between the expressions of H r and H l inEqs. (2), the functions { A ql , B ql } turn out to be related to those of the right chain as A ql ( k ) = B − qr ( k ) and B ql ( k ) = A − qr ( k ) , (A.12)which follows from the fact that the left chain has boundary conditions dual to those of the right chain, i.e., A l ( k ) = B − Nl ( k ) = H = − J (cid:32) N − ∑ q = − N + σ xq σ xq + + h N ∑ q = − N + σ zq (cid:33) , (A.13)which has the same expression as H r in Eq. (2) but with N → N and q → q + N and therefore it can bediagonalized as explained above. More precisely, by expanding Eq. (8) for k n in the thermodynamic limit N → + ∞ , after the replacement N → N , one finds k n = π N − π n N + N arctan (cid:18) sin ( π n N ) cos ( π n N ) − h (cid:19) + O (cid:18) N (cid:19) . (A.14)9Accordingly, k n = π n /2 N → k ∈ ( π ) at the lowest order in 1/ N . By using Eq. (A.14) into Eq. (A.9) and aftershifting the lattice index q → q + N , one readily observes that the analogous of Eqs. (A.7) for H are satisfied bytwo sets of functions A q ( k ) , B q ( k ) and A q ( k ) , B q ( k ) according to the parity of the index n labelling the discretemomenta before taking the thermodynamic limit: A q ( k ) = (cid:114) π sin (cid:18) qk − f ( k ) + k (cid:19) , B q ( k ) = (cid:114) π sin (cid:18) qk + f ( k ) − k (cid:19) , (A.15) A q ( k ) = (cid:114) π cos (cid:18) qk − f ( k ) + k (cid:19) , B q ( k ) = (cid:114) π cos (cid:18) qk + f ( k ) − k (cid:19) . (A.16)As a consequence, we can construct, via a Bogoliubov transformation similar to the one introduced in Eqs. (5)and (A.11), the two following operators: Ψ ( k ) = + ∞ ∑ q = − ∞ (cid:104) ω q ( k ) c q + ξ q ( k ) c † q (cid:105) , Ψ ( k ) = + ∞ ∑ q = − ∞ (cid:104) ω q ( k ) c q + ξ q ( k ) c † q (cid:105) , (A.17)in terms of which H is diagonal. The functions ω q ( k ) and ξ q ( k ) are determined in the same way as in Eq. (A.6).Due to translational symmetry of the total chain for N → ∞ , it is convenient the look for a linear combination Ψ R , L ( k ) of the operators Ψ ( k ) , which transforms according to Eq. (13) under spatial translations, which turnsout to be given by (cid:18) Ψ R ( k ) Ψ L ( k ) (cid:19) = e − i f ( k ) − k √ (cid:18) − i i (cid:19) (cid:18) Ψ ( k ) Ψ ( k ) (cid:19) . (A.18)These two operators can be interpreted as fermionic quasi-particles excitations delocalized along the chain. Inparticular, from the explicit expression of the functions in Eq. (A.16), one eventually finds Ψ R , L ( k ) = + ∞ ∑ q = − ∞ (cid:104) c q ω qR , L ( k ) + c † q ξ qR , L ( k ) (cid:105) , A qR ( k ) = (cid:114) π e i ( − qk + k ) , B qR ( k ) = (cid:114) π e i [ − qk + k − f ( k )] , A qL ( k ) = (cid:114) π e i [ qk − f ( k )] , B qL ( k ) = (cid:114) π e iqk . (A.19)For future reference, matrix elements between the pre- and post-quench single-particle fermionic states arecalculated here, since they play a fundamental role in determining the dynamics of transport properties discussedfurther below. In order to do this, one can invert Eq. (5) and express the local fermionic operators c q with q =
1, . . . , N in terms of φ r ( k ) , i.e., of the pre-quench operators: c q = (cid:90) π dk (cid:104) φ r ( k ) ω qr ( k ) + φ † r ( k ) ξ qr ( k ) (cid:105) , (A.20)where the coefficients of the linear combination of φ r ( k ) and φ † r ( k ) are determined such that to preserve thecanonical fermionic anti-commutation relation { c q , c † q (cid:48) } = δ q , q (cid:48) ; in particular, this implies the following complete-ness relation (cid:90) π dk (cid:104) ω qr ( k ) ω q (cid:48) r ( k ) + ξ qr ( k ) ξ q (cid:48) r ( k ) (cid:105) = δ q , q (cid:48) , (A.21)0for ω qr ( k ) and ξ qr ( k ) . Substituting Eq. (A.20) into Eqs. (A.17) and doing the same with r → l one can express theright- and left-moving fermions in terms of the pre-quench operators: (cid:18) Ψ R ( k ) Ψ L ( k ) (cid:19) = (cid:90) π dk (cid:48) (cid:20) m ( k , k (cid:48) ) (cid:18) φ r ( k (cid:48) ) φ l ( k (cid:48) ) (cid:19) + m ( k , k (cid:48) ) (cid:18) φ † r ( k (cid:48) ) φ † l ( k (cid:48) ) (cid:19)(cid:21) , (A.22)where m ( k , k (cid:48) ) and m ( k , k (cid:48) ) are 2 × m i ( k , k (cid:48) ) = (cid:18) m i , Rr ( k , k (cid:48) ) m i , Rl ( k , k (cid:48) ) m i , Lr ( k , k (cid:48) ) m i , Ll ( k , k (cid:48) ) (cid:19) with i = {
1, 2 } . (A.23)Each matrix element in m ( k , k (cid:48) ) can be expressed as a series involving the functions in Eqs. (A.16) and (A.6).In particular m Rr and m Rl turn out to be the relevant ones for the calculation of the heat current presented inSec. III and they are given by m Rr ( k , k (cid:48) ) = + ∞ ∑ q = (cid:104) ω qR ( k ) ω qr ( k (cid:48) ) + ξ qR ( k ) ξ qr ( k (cid:48) ) (cid:105) , (A.24a) m Rl ( k , k (cid:48) ) = − ∞ ∑ q = (cid:104) ω qR ( k ) ω ql ( k (cid:48) ) + ξ qR ( k ) ξ ql ( k (cid:48) ) (cid:105) . (A.24b)As anticipated, these coefficients can indeed be interpreted as matrix elements: in fact, given the pre-quenchsingle particle state | k (cid:48) (cid:105) r , l ≡ φ † r , l ( k (cid:48) ) | (cid:105) and the post-quench one | k (cid:105) R = Ψ † R ( k ) | (cid:105) , where | (cid:105) is the ground state ofthe complete chain, from Eq. (A.22) the various coefficients are recognized to correspond to the following scalarproduct, where α ∈ { r , l } , R (cid:104) k | k (cid:48) (cid:105) α = (cid:104) | Ψ R ( k ) φ † α ( k (cid:48) ) | (cid:105) = (cid:90) π dk (cid:48)(cid:48) (cid:104) | m R α ( k , k (cid:48)(cid:48) ) φ α ( k (cid:48)(cid:48) ) φ † α ( k (cid:48) ) | (cid:105) = m R α ( k , k (cid:48) ) . (A.25) Appendix B. Space-time scaling limit
In this appendix we report the details of the calculations leading to the expression of the energy current J E ( x , t ) in Eq. (27). One starts by substituting Eq. (24) into Eq. (21) and by evolving in time the corresponding operatorunder the unitary dynamics prescribed by H : j E ( x , t ) = ihJ (cid:16) c † x + c x − c † x c x + (cid:17) = ihJ (cid:90) π − π dk (cid:90) π − π dk (cid:48) (cid:16) e − ik − e ik (cid:48) (cid:17) e i [ ε ( k ) − ε ( k (cid:48) )] t Ψ † R ( k ) Ψ R ( k (cid:48) ) (cid:110) ω xR ( k ) (cid:2) ω xR ( k (cid:48) ) (cid:3) ∗ + ξ xR ( k ) (cid:2) ξ xR ( k (cid:48) ) (cid:3) ∗ (cid:111) + ihJ (cid:90) π − π dk (cid:90) π − π dk (cid:48) (cid:16) e − ik − e − ik (cid:48) (cid:17) e i [ ε ( k )+ ε ( k (cid:48) )] t Ψ † R ( k ) Ψ † R ( k (cid:48) ) ω xR ( k ) ξ xR ( k (cid:48) )+ ihJ (cid:90) π − π dk (cid:90) π − π dk (cid:48) (cid:16) e ik − e ik (cid:48) (cid:17) e − i [ ε ( k )+ ε ( k (cid:48) )] t Ψ R ( k ) Ψ R ( k (cid:48) ) (cid:2) ω xR ( k (cid:48) ) (cid:3) ∗ ξ xR ( k ) , (B.1)where we have used the fermionic anticommutation relation { Ψ R ( k ) , Ψ † R ( k (cid:48) ) } = δ ( k − k (cid:48) ) . As we shall explainlater in this appendix the last two terms of Eq. (B.1) are vanishing in the limit k = k (cid:48) and therefore they will notcontribute in the space-time scaling limit that will be introduced below. As a consequence, we will henceforthmainly consider the first term of the sum in Eq. (B.1). The integral is calculated over the square domain [ − π , π ] × [ − π , π ] , which for brevity, will not be indicated further below.The trace in Eq. (22) can be calculated based on the knowledge of Tr [ Ψ † R ( k ) Ψ R ( k (cid:48) ) ρ ] ,Tr [ Ψ † R ( k ) Ψ † R ( k (cid:48) ) ρ ] andTr [ Ψ R ( k ) Ψ R ( k (cid:48) ) ρ ] . In turn, considering the first contribution, one exploits the change of basis in Eq. (A.22) bywriting: Ψ † R ( k ) Ψ R ( k (cid:48) ) = (cid:90) π dk dk (cid:104) m ∗ Rr ( k , k ) m Rr ( k (cid:48) , k ) Φ † r ( k ) Φ r ( k ) + m ∗ Rl ( k , k ) m Rl ( k (cid:48) , k ) Φ † l ( k ) Φ l ( k ) + ... (cid:105) ,(B.2)where ”...” denotes terms which are regular within the integration domain and therefore, as explained furtherbelow, do not contribute to the integral. Taking the trace of this expression as in Eq. (22) and remembering that(since from Eq. (10) one has that the operators Φ r , l represent the excitations of a free fermion system)Tr [ ρ Φ † α ( k ) Φ γ ( k (cid:48) )] = δ αγ δ ( k − k (cid:48) ) f α ( k ) = δ αγ δ ( k − k (cid:48) ) e β α (cid:101) ( k ) + α ∈ { r , l } (B.3)one gets Tr [ ρ Ψ † R ( k ) Ψ R ( k (cid:48) )] = I ( k , k (cid:48) ) = I r ( k , k (cid:48) ) + I l ( k , k (cid:48) ) , (B.4)where we defined I α ( k , k (cid:48) ) = (cid:90) π dk m ∗ R α ( k , k ) m R α ( k (cid:48) , k ) f α ( k ) . (B.5)Within the space-time scaling limit introduced in Sec. III to calculate Eq. (25) we are led to consider the pointswhere the phase ϕ x , t ( k , k (cid:48) ) = [ ε ( k ) − ε ( k (cid:48) )] t + x ( k (cid:48) − k ) , appearing in the expression of J E ( x , t ) which follows fromcalculating the trace of Eq. (B.1) according to Eq. (22), is stationary, taking into account the explicit expression of ω xR ( k ) and ξ xR ( k ) from Eq. (A.19), one obtains: v g ( k ) t − x = − v g ( k (cid:48) ) t + x =
0, (B.6)where we introduced the group velocity v g ( k ) of the excitations as already done in Eq. (27). Depending on thevalue of the ratio v ≡ x / t , each of the two equations in Eq. (B.6) has either no or two solutions k + s ( v ) and k − s ( v ) (as shown in the main text in Fig. 8). In the latter case we have four possible stationary points for the pair ( k , k (cid:48) ) ,i.e., ( k + s , k + s ), ( k − s , k − s ), ( k + s , k − s ), ( k − s , k + s ). Focusing now on Eq. (B.5), one needs to specify further the structure ofthe matrix elements m , Rr ( k , k (cid:48) ) and m , Rl ( k , k (cid:48) ) involved. From Eq. (A.24) calculating the series one finds m Rr ( k , k (cid:48) ) = π i (cid:34) e i [ k − f ( k (cid:48) )] − e i ( k (cid:48) − k + i δ ) − e i [ k + f ( k (cid:48) )] − e − i ( k + k (cid:48) − i δ ) + e i [ k − f ( k )] − e i ( k (cid:48) − k + i δ ) − e i [ k − f ( k )] − e − i ( k + k (cid:48) − i δ ) (cid:35) , m Rl ( k , k (cid:48) ) = π i (cid:34) e − i [ f ( k )+ f ( k (cid:48) )] − e i ( k (cid:48) + k + i δ ) − e i [ f ( k (cid:48) ) − f ( k )] − e i ( k − k (cid:48) + i δ ) + − e i ( k (cid:48) + k + i δ ) − − e i ( k − k (cid:48) + i δ ) (cid:35) , (B.7)here δ > m R α ( k , k (cid:48) ) = − m R α ( k , − k (cid:48) ) and f α ( k ) = f α ( − k ) , with α ∈ { l , r } , one finds fromEq. (B.5) that: I α ( k , k (cid:48) ) = = (cid:90) π − π dk m ∗ R α ( k , k ) m R α ( k (cid:48) , k ) f α ( k ) == (cid:73) C dz m ∗ R α ( k , − i ln ( z )) m R α ( k (cid:48) , − i ln ( z )) f α ( − i ln z ) iz , (B.8)2where the original integral has been extended to the complex plane, along the circumference C centred in theorigin and with unitary radius, via the change of variable z = e ik . Since the integration path is a closed contour, I ( k , k (cid:48) ) is determined by the singularities of the integrand which are located inside C . Applying the residuetheorem to Eq. (B.8) with the expressions in Eq. (B.7), the final form of I ( k , k (cid:48) ) from Eq. (B.4) is: I ( k , k (cid:48) ) = π i (cid:32) f β l ( k ) + f β l ( k (cid:48) ) k − k (cid:48) − i δ − f β r ( k ) + f β r ( k (cid:48) ) k − k (cid:48) + i δ (cid:33) + regular terms as k → k (cid:48) . (B.9)In the space-time scaling limit introduced to calculate Eq. (B.1) one expects the double integral in k and k (cid:48) to bedominated by the pairs ( k , k (cid:48) ) solution of Eq. (B.6) where the phase ϕ ( k , k (cid:48) ) is stationary as well as possible singu-larities in I ( k , k (cid:48) ) . By inspection of Eq. (B.9) one notices that the only saddle points where I ( k , k (cid:48) ) is stationary arethe couples ( k + s , k + s ), ( k − s , k − s ), therefore one concludes that the double integral in k and k (cid:48) of Eq. (B.1) is dominatedby the region where k (cid:39) k (cid:48) . This is also the reason why we omitted in Eq. (B.9) terms which are not singular as k → k (cid:48) since they do not contribute in the space-time scaling limit, as just explained. For the same reason, in thespace-time scaling limit, we can neglect the second and the third term of Eq. (B.1) since they vanish in the limit k → k (cid:48) as anticipated at the beginning of th e section.More specifically, considering the term I ( k , k (cid:48) ) ≡ Tr [ ρ Ψ † R ( k ) Ψ † R ( k (cid:48) )] , it can be evaluated following the samesteps considered for I ( k , k (cid:48) ) Ψ † R ( k ) Ψ † R ( k (cid:48) ) = (cid:90) π dk dk (cid:104) m ∗ Rr ( k , k ) m ∗ Rr ( k (cid:48) , k ) Φ † r ( k ) Φ r ( k ) + m ∗ Rl ( k , k ) m ∗ Rl ( k (cid:48) , k ) Φ † l ( k ) Φ l ( k ) + ... (cid:105) (B.10)where the new matrix elements m Rr ( k , k (cid:48) ) , m Rl ( k , k (cid:48) ) are given by the series m Rr ( k , k (cid:48) ) = + ∞ ∑ q = (cid:104) ω qR ( k ) ξ qr ( k (cid:48) ) + ξ qR ( k ) ω qr ( k (cid:48) ) (cid:105) , m Rl ( k , k (cid:48) ) = − ∞ ∑ q = (cid:104) ω qR ( k ) ξ ql ( k (cid:48) ) + ξ qR ( k ) ω ql ( k (cid:48) ) (cid:105) , (B.11)that turns out to sum to m Rr ( k , k (cid:48) ) = π i (cid:34) e i [ k − f ( k (cid:48) )] − e i ( k (cid:48) − k + i δ ) − e i [ k + f ( k (cid:48) )] − e − i ( k + k (cid:48) − i δ ) − e i [ k − f ( k )] − e i ( k (cid:48) − k + i δ ) + e i [ k − f ( k )] − e − i ( k + k (cid:48) − i δ ) (cid:35) , m Rl ( k , k (cid:48) ) = π i (cid:34) e − i [ f ( k )+ f ( k (cid:48) )] − e i ( k (cid:48) + k + i δ ) − e i [ f ( k (cid:48) ) − f ( k )] − e i ( k − k (cid:48) + i δ ) − − e i ( k (cid:48) + k + i δ ) + − e i ( k − k (cid:48) + i δ ) (cid:35) . (B.12)From Eq. (B.10) I ( k , k (cid:48) ) ≡ Tr [ ρ Ψ † R ( k ) Ψ † R ( k (cid:48) )] = I l ( k , k (cid:48) ) + I r ( k , k (cid:48) ) can be eventually calculated: I r ( k , k (cid:48) ) = π (cid:20) e − i [ k − f ( k )] (cid:18) e − i [ k (cid:48) − f ( k )] f β r ( k ) i ( k − k (cid:48) ) − e − i [ k (cid:48) − f ( k (cid:48) )] f β r ( k ) i ( k − k (cid:48) ) (cid:19)(cid:21) , I l ( k , k (cid:48) ) = π (cid:20) e i [ f ( k (cid:48) ) − f ( k )] f β l ( k ) i ( k − k (cid:48) ) − f β l ( k ) i ( k − k (cid:48) ) (cid:21) . (B.13)Once the trace over ρ in Eq. (B.3) has been taken, it is useful to define Q = k − k (cid:48) and K = ( k + k (cid:48) ) /2 andexploit the fact that Q is infinitesimal in order to Taylor expand the dependence of the integral around Q = I r ( k , k (cid:48) ) and I l ( k , k (cid:48) ) , once they are multiplied by the factor e − ik − e − ik (cid:48) in Eq. (B.1), vanish in the space-time scaling limit.The very same reasoning, that here is not reported again for brevity, can be also applied to the term determinedby Tr [ ρ Ψ R ( k ) Ψ R ( k (cid:48) )] .3Concentrating as a consequence on the first term of Eq. (B.1), in particular, with reference to the phase ε ( k ) − ε ( k (cid:48) ) ,one has: ε ( k ) − ε ( k (cid:48) ) = ε (cid:18) K + Q (cid:19) − ε (cid:18) K − Q (cid:19) = Qv g ( K ) + O ( Q ) . (B.14)In terms of the new variables, the heat current takes the form (see Eq. (B.1) and (22)): J E ( x , t ) = hJ π (cid:90) π − π dK (cid:90) + ∞ − ∞ dQ ( sin ( K )) e iQ [ v g ( K ) t − x ] π i (cid:18) f β l ( K ) Q − i δ − f β r ( K ) Q + i δ (cid:19) . (B.15)Recalling the integral representation of the Heaviside step function Θ ( x ) = lim δ → + (cid:82) ∞ − ∞ dy π i e ixy y − i δ , we eventuallyfind: J E ( x , t ) = (cid:90) π − π dk π ε ( k ) v g ( k ) (cid:2) f β l ( k ) Θ ( v g ( k ) t − x ) − f β r ( k ) Θ ( x + v g ( k ) t ) (cid:3) , (B.16)which represents the exact profile of the energy current in time and space within the space-time scaling regime.The integral over k above can be restricted to the domain [ π ] by using the fact that ε ( − k ) = ε ( k ) and therefore v g ( − k ) = − v g ( k ) . The resulting expression clearly shows that J E ( − x , t ) = J E ( x , t ) and therefore x in the r.h.s. ofEq. (B.16) can be replaced by | x | . After introducing the convenient scaling variable v = x / t and using the factthat Θ ( v − v g ( k )) = − Θ ( v g ( k ) − v ) , one has J E ( v ) = (cid:90) π dk π ε ( k ) v g ( k ) (cid:2) f β l ( k ) − f β r ( k ) (cid:3) Θ ( v g ( k ) − | v | ) = Θ ( v max − | v | ) (cid:90) k + s ( v ) k − s ( v ) dk π ε ( k ) v g ( k ) (cid:2) f β l ( k ) − f β r ( k ) (cid:3) (B.17)where k ± s ( v ) are the solutions of Eq. (56), as shown in Fig. 8. For later convenience, let us consider a generalization— henceforth denoted by J a ( v ) — of the previous expression in which ε ( k ) in the integrand is replaced by ε a ( k ) ,with a = J = J E , it turns out that J = J N , i.e., it corresponds to the particle currentdiscussed in Sec. III B, see Eq. (46). Since v g ( k ) = ε (cid:48) ( k ) one can perform the change of variable k (cid:55)→ ε ( k ) endingup with J a ( v ) = Θ ( v max − | v | ) (cid:90) ε + ( v ) ε − ( v ) d ε π ε a (cid:2) f β l ( ε ) − f β r ( ε ) (cid:3) , with ε ± ( v ) ≡ ε ( k ± s ( v )) , (B.18)which after some rescaling can be written as J a ( v ) = Θ ( v max − | v | ) [ J a ( β l , v ) − J a ( β r , v )] , (B.19)with J a ( β , v ) = G a ( βε − ( v )) − G a ( βε + ( v )) πβ a + , (B.20)where we introduced G a ( x ) ≡ (cid:90) + ∞ x d ξ ξ a e ξ + a = a =
0. In order to make this expressionmore explicit, one still need to determine ε ± ( v ) defined in Eq. (B.18), which is readily done by observing that thesquare of Eq. (56) is a second order equation for cos ( k ± s ( v )) , the solutions of which arecos ( k ± s ( v )) = v J h ∓ (cid:115)(cid:18) − v J h (cid:19) (cid:18) − v J (cid:19) . (B.22)4Replacing this expression into the dispersion relation Eq. (7), one obtains ε ± ( v ) = ε > ( v ) ± ε < ( v ) , (B.23)where (see also Eq. (31)) ε > ( v ) = (cid:113) [ J max ( h )] − v and ε < ( v ) = (cid:113) [ J min ( h )] − v = (cid:113) v max − v . (B.24)Inserting Eq. (B.23) into Eq. (B.19) (see also Eq. (B.20)) for a = J a ( v ) as v approaches the edges at ± v max , note that,correspondingly, ε < ( v ) approaches zero ∝ (cid:112) | v − v max | while ε > ( v ) = ε > ( v max ) + O ( | v − v max | ) . Accordingly,the leading term of the current J a in Eq. (B.19), upon approaching the edge | v | (cid:46) v max , can be obtained from theexpansion G a ( β [ ε > ( v ) ± ε < ( v )]) = G a ( βε > ( v max )) ± βε < ( v ) G (cid:48) a ( βε > ( v max )) + O (( v max − | v | ) ) , (B.25)with G (cid:48) a ( x ) = − x a / ( e x + ) [see Eq. (B.21)] and therefore J a , according to Eqs. (B.19) and (B.20), can be expressedas J a ( v ) = C a (cid:16) v max − v (cid:17) + O (( v max − | v | ) ) , (B.26)where C a = G (cid:48) a ( β r ε > ( v max )) πβ ar − G (cid:48) a ( β l ε > ( v max )) πβ al = J a | − h | a /2 π (cid:34) e β l J √ | − h | + − e β r J √ | − h | + (cid:35) . (B.27)Note that the constant C a reported above vanishes if the field is turned to its critical value, i.e., if h =
1. In thiscase, one has to consider in Eq. (33) with ε > ( v ) = ε < ( v ) the next order in the expansion of G a ( βε < ( v )) for ε < ( v ) →
0. Taking into account that G ( x ) = π /12 − x /4 + x /12 + O ( x ) [see Eq. (B.21)], one finds J ( β , v ) = ε < ( v ) π − β π ε < ( v ) + O ( ε < ( v )) , (B.28)which, inserted into Eq. (33), results into the different form of the edge behaviour reported in Eq. (36). We alsonote that C can be equivalently rewritten in terms of the single-particle energy spectrum and the Fermi-Diracdistributions as: C = ε ( k s ) π ( f β l ( k s ) − f β r ( k s )) ; (B.29)the vanishing of C — which implies different behaviour of the energy current J E ( v ) as v approaches the edge | v | → v − max for h = ε ( k s ) vanishes in this case.Proceeding in exactly the same way, since G ( x ) = ln2 − x /2 + x /8 − x /192 + O ( x ) , for the particle currentone has J ( β , v ) = ε < ( v ) π − β π ε < ( v ) + O ( ε < ( v )) , (B.30)which gives the critical edge behavior of Eq. (49), once inserted into Eq. (47).5 Appendix C. Stationary phase calculation at the front edge
A Correlation function
In order to determine the two-point correlation function near the right edge of the front with x , y (cid:39) v max t (forthe front to the left with x , y (cid:39) − v max t the same procedure applies) we use a stationary-phase approximatiomstarting from the expression (cid:104) γ † x ( t ) γ y ( t ) (cid:105) = (cid:90) π − π (cid:90) π − π dkdk (cid:48) π e i [ φ x , t ( k ) − φ y , t ( k (cid:48) )] I ( k , k (cid:48) )= (cid:90) π − π (cid:90) π − π dkdk (cid:48) π e i [ φ x , t ( k ) − φ y , t ( k (cid:48) )] π i (cid:34) f β l ( k ) + f β l ( k (cid:48) ) k − k (cid:48) − i δ + f β r ( k ) + f β r ( k (cid:48) ) k (cid:48) − k − i δ (cid:35) = (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l + (cid:104) γ † x ( t ) γ y ( t ) (cid:105) r (C.1)which can be obtained following exactly the same steps which led to Eq. (25), with I ( k , k (cid:48) ) given in Eq. (B.9) and φ x , t ( k ) = ε ( k ) t − kx with analogous definition for φ y , t ( k (cid:48) ) . Since, as explained in Sec. IV, the stationary phaseequation (56) admits only one solution k s for x , y (cid:39) v max t , such that ε (cid:48) ( k s ) = v g ( k s ) = v max and ε (cid:48)(cid:48) ( k s ) = v (cid:48) g ( k s ) = φ (cid:48)(cid:48) x , t ( k ) and φ (cid:48)(cid:48) y , t ( k ) vanish, we need to expand them beyond the second order; actuallywe expand below up to the fourth order in order to account also the first correction beyond the leading behaviour,finding, e.g., φ x , t ( k ) = ε ( k s ) t − k s x + ( k − k s )( v max t − x ) − v max ( k − k s ) t + ε ( ) ( k s ) ( k − k s ) t + O (( k − k s ) ) , (C.2)and analogous for φ y , t ( k (cid:48) ) , where we used the fact that ε ( ) ( k ) = − ε (cid:48) ( k )[ + ε (cid:48)(cid:48) ( k ) / ε ( k )] , which implies (for h (cid:54) = h = ε ( ) ( k s ) = − v max and ε ( ) ( k s ) = v max / ε ( k s ) = v max / ( J (cid:112) | h − | ) (see Eq. (B.23)). Accordingly, the first term in Eq. (C.1) can be written as (cid:104) γ † x ( t ) γ y ( t ) (cid:105) l = e − ik s ( x − y ) (cid:90) + ∞ − ∞ d (cid:101) k π (cid:90) + ∞ − ∞ d (cid:101) q π f β l ( (cid:101) k ) + f β l ( (cid:101) q ) i ( (cid:101) k − (cid:101) q − i δ ) × exp (cid:40) − i (cid:34)(cid:101) k ( x − v max t ) − (cid:101) q ( y − v max t ) + (cid:101) k v max t − (cid:101) k ε ( ) ( k s ) t − (cid:101) q v max t + (cid:101) q ε ( ) ( k s ) t (cid:35)(cid:41) ,(C.3)with (cid:101) k = k − k s and (cid:101) q = k (cid:48) − k s . Similarly, the Fermi-Dirac functions have also to be expanded around thestationary point f β l ( (cid:101) k ) = f β l ( k s ) + f (cid:48) β l ( k s ) (cid:101) k + f (cid:48)(cid:48) β l ( k s ) (cid:101) k + O ( (cid:101) k ) . (C.4)Since one expects the term proportional to (cid:101) k in the exponential to be the leading order, it is convenient to rescalethe variables as K = ( v max t /2 ) (cid:101) k , Q = ( v max t /2 ) (cid:101) q and introduce the scaled coordinates X and Y as in6Eq. (59), such that the exponentials in Eq. (C.3) can be written as:exp (cid:32) − iKX − iK + i ε ( ) ( k s ) K v max t (cid:33) = exp (cid:18) − iKX − iK (cid:19) (cid:32) + i ε ( ) ( k s ) K v max t (cid:33) + O (cid:16) ( v max t ) − (cid:17) ,(C.5)exp (cid:32) + iQY + iQ − i ε ( ) ( k s ) Q v max t (cid:33) = exp (cid:18) + iQY + iQ (cid:19) (cid:32) − i ε ( ) ( k s ) Q v max t (cid:33) + O (cid:16) ( v max t ) − (cid:17) .(C.6)Inserting Eqs. (C.4), (C.5), and (C.6) into Eq. (C.3) and keeping terms up to order t − we end up with the resultEq. (58) of the main text with K l ( X , Y ) = i (cid:34) f (cid:48) β l ( k s ) f β r ( k s ) (cid:18) ∂ K A ( X , Y ) ∂ X − ∂ K A ( X , Y ) ∂ Y (cid:19) + ε ( ) ( k s ) (cid:18) ∂ K A ( X , Y ) ∂ X − ∂ K A ( X , Y ) ∂ Y (cid:19)(cid:35) , (C.7)where we have used the integral representation of the Airy kernel K A ( X , Y ) = (cid:90) + ∞ − ∞ dK π (cid:90) + ∞ − ∞ dQ π e − iKX − iK /3 + iQY + iQ /3 i ( K − Q − i δ ) , (C.8)which coincides with Eq. (60) once we observe that − ( ∂ X + ∂ Y ) K A ( X , Y ) = Ai ( X ) Ai ( Y ) with the usual integralrepresentation of the Airy function: Ai ( X ) = (cid:90) + ∞ − ∞ dK π e iKX + iK /3 . (C.9)For completeness we report here the well-known asymptotic behaviours of the Airy function that have been usedin the main text: Ai ( X → + ∞ ) (cid:39) √ π X exp (cid:26) − X (cid:27) , (C.10)Ai ( X → − ∞ ) (cid:39) √ π | X | cos (cid:18) − | X | + π (cid:19) . (C.11) B Energy current
In order to study the behaviour of the energy current J E ( x , t ) near the edge x (cid:39) v max t , the procedure iscompletely analogous to the one presented above; in particular, starting from Eq. (25), with I ( k , k (cid:48) ) given byEq. (B.9) and ϕ x , t ( k , k (cid:48) ) = φ x , t ( k ) − φ x , t ( k (cid:48) ) , we can expand (as done in Appendix C A, see Eq. (C.2)) the phases φ x , t ( k ) and φ x , t ( k (cid:48) ) around the stationary point k s , as we are interested in the leading-order correction to thespace-time scaling limit. Following the same steps as above, one finds Eq. (61), where K A is given by Eq. (63),obtained by taking the limit X → Y of Eq. (60).As outlined in the main text, on the other hand, the expression in Eq. (61) vanishes when h is set to its criticalvalue h c =
1. Accordingly, one has to expand Eq. (25) up to the first non-vanishing order. In particular, for h = h c , the stationary point k s of the phase φ x , t ( k ) turns out to be k s =
0, approached from above for x (cid:39) v max t and from below for x (cid:39) − v max t (with v max = J ). Correspondingly, the odd derivatives of the dispersion relation7 ε ( k ) become discontinuous at k = ε ( k ) = v max | sin ( k /2 ) | , and therefore one has to consider the properlimits, i.e., lim k → ± ε (cid:48) ( k ) = v g ( ± ) = ± v max and lim k → ± ε ( ) ( k ) = ∓ v max /4, (C.12)while all the even derivatives vanish. As a consequence, by expanding up to third order in the phase φ x , t ( k ) , for x (cid:39) v max t , one finds (instead of Eq. (C.2) with k s = φ x , t ( k ) = k ( v max t − x ) − v max k t + O ( k )= − K − XK + O (cid:16) t − (cid:17) , (C.13)where we defined k = ( v max t /8 ) − K and the scaling variable X = ( x − v max t ) / ( v max t /8 ) which is analogousto Eq. (59), except for a numerical factor due to the fact that ε ( ) ( k s ) at the critical point is no longer − v max asin the non-critical case, but it is given by Eq. (C.12). Analogous expansion is done for φ x , t ( k (cid:48) ) = − Q /3 − XQ + O ( t − ) , where k (cid:48) = ( v max t /8 ) − Q . Keeping into account that the factor e − ik (cid:48) − e ik in Eq. (25) must beexpanded up to first order in K and Q , since it vanishes identically at the lowest order, e − ik (cid:48) − e ik = − i ( K + Q )( v max t /8 ) − + O (cid:16) t − (cid:17) , (C.14)we get the following expression for Eq. (25) J E ( X , t ) = (cid:18) v max t (cid:19) J (cid:90) + ∞ − ∞ dK π (cid:90) + ∞ − ∞ dQ π ( Q + K ) e − iKX − iK /3 e iQX + iQ /3 × (cid:34) f β l ( K ( v max t /8 ) − ) + f β l ( Q ( v max t /8 ) − ) i ( K − Q − i δ ) − f β r ( K ( v max t /8 ) − ) + f β r ( Q ( v max t /8 ) − ) i ( K − Q − i δ ) (cid:35) .(C.15)In order to determine the first non-vanishing order, one needs to expand the Fermi-Dirac distributions f β l ( K ( v max t /8 ) − ) + f β l ( Q ( v max t /8 ) − ) − f β r ( K ( v max t /8 ) − ) − f β r ( Q ( v max t /8 ) − )) == [ f (cid:48) β l ( + ) − f (cid:48) β r ( + )]( K + Q )( v max t /8 ) − + O (cid:16) t − (cid:17) (C.16)with f β ( + ) = − β v max /4. By combining Eqs. (C.15) and (C.16), one eventually finds J E ( X , t ) = t v max ( β r − β l ) K c ( X ) , (C.17)where K c ( X ) = (cid:90) + ∞ − ∞ dK π (cid:90) + ∞ − ∞ dQ π ( Q + K ) e − iKX − iK /3 e iQX + iQ /3 i ( K − Q − i δ ) . (C.18)One can make the expression of K c more explicit by taking a derivative with respect to X and then by usingEqs. (C.9), (67), and (63) in order to show that ∂ K c ( X ) ∂ X = − K A ( X , X ) , (C.19)which renders Eqs. (68) and, by integration, Eq. (66). [1] Pasquale Calabrese, Fabian H L Essler, and GiuseppeMussardo. Introduction to quantum integrability in out of equilibrium systems. Journal of Statistical Me- chanics: Theory and Experiment , 2016(6):064001, 2016.[2] Jens Eisert, Mathis Friesdorf, and Christian Gogolin.Quantum many-body systems out of equilibrium. Na-ture Physics , 11(2):124–130, 2015.[3] Anatoli Polkovnikov, Krishnendu Sengupta, Alessan-dro Silva, and Mukund Vengalattore. Nonequilibriumdynamics of closed interacting quantum systems.
Re-views of Modern Physics , 83(3):863, 2011.[4] Ehud Altman. Non equilibrium quantum dynamics inultra-cold quantum gases. preprint arXiv:1512.00870 ,2015.[5] Pasquale Calabrese. Non-equilibrium dynamics of iso-lated quantum systems. In
EPJ Web of Conferences , vol-ume 90, page 08001. EDP Sciences, 2015.[6] Lev Vidmar and Marcos Rigol. Generalized gibbs en-semble in integrable lattice models.
Journal of StatisticalMechanics: Theory and Experiment , 2016(6):064007, 2016.[7] Fabian HL Essler and Maurizio Fagotti. Quench dy-namics and relaxation in isolated integrable quantumspin chains.
Journal of Statistical Mechanics: Theory andExperiment , 2016(6):064002, 2016.[8] Enej Ilievski, Marko Medenjak, Tomaž Prosen, andLenart Zadnik. Quasilocal charges in integrable lat-tice systems.
Journal of Statistical Mechanics: Theory andExperiment , 2016(6):064008, 2016.[9] Denis Bernard and Benjamin Doyon. Energy flowin non-equilibrium conformal field theory.
Journal ofPhysics A: Mathematical and Theoretical , 45(36):362001,2012.[10] Herbert Spohn and Joel L Lebowitz. Stationary non-equilibrium states of infinite harmonic systems.
Com-munications in Mathematical Physics , 54(2):97–120, 1977.[11] Pasquale Calabrese, Fabian HL Essler, and MaurizioFagotti. Quantum quench in the transverse-field Isingchain.
Physical Review Letters , 106(22):227203, 2011.[12] Pasquale Calabrese, Fabian HL Essler, and MaurizioFagotti. Quantum quench in the transverse field Isingchain: I. time evolution of order parameter correlators.
Journal of Statistical Mechanics: Theory and Experiment ,2012(07):P07016, 2012.[13] Mario Collura, Spyros Sotiriadis, and Pasquale Cal-abrese. Quench dynamics of a Tonks-Girardeau gasreleased from a harmonic trap.
Journal of Statistical Me-chanics: Theory and Experiment , 2013(09):P09025, 2013.[14] Denis Bernard and Benjamin Doyon. Non-equilibriumsteady states in conformal field theory.
Annales HenriPoincaré , 16(1):113–161, 2015.[15] Benjamin Doyon. Nonequilibrium density matrix forthermal transport in quantum field theory. preprintarXiv:1212.1077 , 2012.[16] Denis Bernard and Benjamin Doyon. Conformalfield theory out of equilibrium: a review. preprintarXiv:1603.07765 , 2016.[17] Viktor Eisler and Zoltán Rácz. Full counting statistics in a propagating quantum front and random matrixspectra.
Physical Review Letters , 110(6):060602, 2013.[18] T Antal, Z Rácz, A Rákos, and GM Schütz. Transportin the XX chain at zero temperature: Emergence of flatmagnetization profiles.
Physical Review E , 59(5):4912,1999.[19] Jacopo Viti, Jean-Marie Stéphan, Jérôme Dubail,and Masudul Haque. Inhomogeneous quenchesin a fermionic chain: exact results. preprintarXiv:1507.08132 , 2015.[20] Mario Collura and Dragi Karevski. Quantum quenchfrom a thermal tensor state: boundary effects andgeneralized gibbs ensemble.
Physical Review B ,89(21):214308, 2014.[21] Mario Collura and Gabriele Martelloni. Non-equilibrium transport in d-dimensional non-interacting fermi gases.
Journal of Statistical Mechanics:Theory and Experiment , 2014(8):P08006, 2014.[22] Viktor Eisler, Florian Maislinger, and Hans Gerd Ev-ertz. Universal front propagation in the quantumIsing chain with domain-wall initial states. preprintarXiv:1610.01540 , 2016.[23] Andrea De Luca, Jacopo Viti, Denis Bernard, and Ben-jamin Doyon. Nonequilibrium thermal transport in thequantum Ising chain.
Physical Review B , 88(13):134301,2013.[24] Craig A Tracy and Harold Widom. Level-spacing dis-tributions and the Airy kernel.
Communications inMathematical Physics , 159(1):151–174, 1994.[25] Gabriele Perfetto. Statistics and ballistic transport froma quench of two ising chains.
Master Thesis, Polytechnicof Turin , July 2016.[26] Subir Sachdev.
Quantum phase transitions . Wiley OnlineLibrary, 2007.[27] Elliott Lieb, Theodore Schultz, and Daniel Mattis. Twosoluble models of an antiferromagnetic chain. In
Con-densed Matter Physics and Exactly Soluble Models , pages543–601. Springer, 2004.[28] Olalla Castro-Alvaredo, Yixiong Chen, BenjaminDoyon, and Marianne Hoogeveen. Thermodynamicbethe ansatz for non-equilibrium steady states: ex-act energy current and fluctuations in integrable qft.
Journal of Statistical Mechanics: Theory and Experiment ,2014(3):P03011, 2014.[29] Denis Bernard and Benjamin Doyon. Time-reversal symmetry and fluctuation relations in non-equilibrium quantum steady states.
Journal of PhysicsA: Mathematical and Theoretical , 46(37):372001, 2013.[30] Benjamin Doyon, Andrew Lucas, Koenraad Schalm,and MJ Bhaseen. Non-equilibrium steady states in theKlein–Gordon theory.
Journal of Physics A: Mathematicaland Theoretical , 48(9):095002, 2015.[31] Christoph Karrasch, Roni Ilan, and JE Moore.Nonequilibrium thermal transport and its relation to linear response. Physical Review B , 88(19):195129, 2013.[32] Alberto Biella, Andrea De Luca, Jacopo Viti, DavideRossini, Leonardo Mazza, and Rosario Fazio. Energytransport between two integrable spin chains.
PhysicalReview B , 93(20):205121, 2016.[33] Olalla A Castro-Alvaredo, Benjamin Doyon, andTakato Yoshimura. Emergent hydrodynamics in inte-grable quantum systems out of equilibrium.
PhysicalReview X , 6(4):041065, 2016.[34] Andrea De Luca, Jacopo Viti, Leonardo Mazza, andDavide Rossini. Energy transport in Heisenberg chainsbeyond the Luttinger liquid paradigm.
Physical ReviewB , 90(16):161101, 2014.[35] Bruno Bertini, Mario Collura, Jacopo De Nardis, andMaurizio Fagotti. Transport in out-of-equilibrium XXZchains: Exact profiles of charges and currents.
PhysicalReview Letters , 117(20):207201, 2016.[36] T Platini and D Karevski. Out of equilibrium processin Ising quantum chains.
Journal of Physics: ConferenceSeries , 40(1):93, 2006.[37] Subir Sachdev and AP Young. Low temperature re-laxational dynamics of the Ising chain in a transversefield.
Physical Review Letters , 78(11):2220, 1997.[38] Tibor Antal, PL Krapivsky, and A Rákos. Logarithmiccurrent fluctuations in nonequilibrium quantum spinchains.
Physical Review E , 78(6):061115, 2008.[39] Felipe Barra Shigeru Ajisaka and Bojan Žunkoviˇc.Nonequilibrium quantum phase transitions in the XY model: comparison of unitary time evolution and re-duced density operator approaches.
New Journal ofPhysics , 16(3):033028, 2014.[40] V Hunyadi, Z Rácz, and L Sasvári. Dynamic scalingof fronts in the quantum XX chain.
Physical Review E ,69(6):066103, 2004.[41] Nicolas Allegra, Jérôme Dubail, Jean-Marie Stéphan,and Jacopo Viti. Inhomogeneous field theory insidethe arctic circle.
Journal of Statistical Mechanics: Theoryand Experiment , 2016(5):053108, 2016.[42] Viktor Eisler. Universality in the full counting statis-tics of trapped fermions.
Physical Review Letters ,111(8):080402, 2013.[43] David S Dean, Pierre Le Doussal, Satya N Majumdar,and Grégory Schehr. Finite-temperature free fermionsand the Kardar-Parisi-Zhang equation at finite time.
Physical Review Letters , 114(11):110402, 2015.[44] David S Dean, Pierre Le Doussal, Satya N Majumdar,and Gregory Schehr. Non-interacting fermions at finitetemperature in a d -dimensional trap: universal corre-lations. preprint arXiv:1609.04366 , 2016.[45] Hugo Touchette. The large deviation approach to sta-tistical mechanics. Physics Reports , 478(1):1–69, 2009.[46] Ming-Chiang Chung, Anibal Iucci, and MA Cazalilla.Thermalization in systems with bipartite eigenmodeentanglement.