Inhomogeneous quenches in the transverse field Ising chain: scaling and front dynamics

We investigate the non-equilibrium dynamics of the transverse field quantum Ising chain evolving from an inhomogeneous initial state given by joining two macroscopically different semi-infinite chains. We obtain integral expressions for all two-point correlation functions of the Jordan-Wigner Majorana fermions at any time and for any value of the transverse field. Using this result, we compute analytically the profiles of various physical observables in the space-time scaling limit and show that they can be obtained from a hydrodynamic picture based on ballistically propagating quasiparticles. Going beyond the hydrodynamic limit, we analyze the approach to the non-equilibrium steady state and find that the leading late time corrections display a lattice effect. We also study the fine structure of the propagating fronts which are found to be described by the Airy kernel and its derivatives. Near the front we observe the phenomenon of energy back-flow where the energy locally flows from the colder to the hotter region.


Introduction
The last decade witnessed an ever increasing interest in the out of equilibrium dynamics of quantum many-body systems [1,2]. To a large extent this is due to the spectacular advances in experiments performed with cold atoms which makes it possible to study the coherent time evolution of quantum systems that can be, to a good precision, considered isolated from their environment. Lattice and continuum systems can be studied and their parameters such as their coupling strengths, the confining potential etc. can be changed in real time, opening the way to study a large variety of dynamical phenomena [3,4,5].
The mechanism of equilibration and thermalization has been in the focus of attention of many theoretical works. The main paradigm has been the so-called quantum quench [6], i.e. a sudden change of a system parameter and the subsequent out of equilibrium time evolution. It was soon realized that integrable systems are peculiar [7] because in general they cannot reach thermal equilibrium, even locally, due to the presence of an extensive number of local conserved quantities [8].
In most of these studies, the initial state was taken to be the ground state of a locally interacting Hamiltonian. To be able to study transport properties, however, one needs to consider either an open system with driving at its boundaries, or in the case of isolated systems, an inhomogeneous initial state. The simplest inhomogeneous setup is arguably provided by the so-called partitioning protocol or tensor product initial state. This means that two systems with macroscopically different properties (temperature, magnetization, etc.) are suddenly connected to each other 1 .
Integrable systems behave in a special way also in this setup. Fourier's law does not hold in them: they can support currents even in the absence of a gradient and they can feature exotic transport properties. In many integrable systems the transport can be ballistic, however, this is not necessary [9,10,11]. In particular, the non-equilibrium steady state (NESS) forming after a long time in the partitioning protocol is typically a translationally invariant state which however supports a finite current. Many works focused on the non-equilibrium steady state in the two-temperature scenario, and in particular on the computation of the expectation value of the energy current in the NESS in the XX [12,13], XY [14,15], Ising [16], XXZ [17,18,19] spin chains, in conformal [20,21,22,23,24] and integrable [25,26] field theories. Fluctuation relations for the currents have also been derived [22,27], and various bounds for the ballistic currents were discussed [28]. Subsystem correlations [29,30,31] and higher dimensional models with similar setups have also been investigated [21,32,33].
Apart from the NESS, the details of the non-equilibrium time evolution and the spacetime dependent profile of physical quantities is also of great interest. When the transport is ballistic, at large times these profiles are given by functions of the scaling variable x/t, i.e. they are constant along a "ray" defined by this ratio. Profiles of the magnetization, particle and energy density were studied in various systems, including the XX spin chain (free hopping fermions) [34,35,36,37,38,39], the XY and Ising spin chain [36,40,41,42,43], continuum free fermions [44], and conformal field theory [45,46]. Inhomogeneous initial states were studied extensively in the XXZ spin chain [9,47,48,38,49,50,51,52,53,54,10]. The fact that the profiles are functions of the scaling variable x/t raises the possibility of a hydrodynamic or semiclassical description in which the space-time dependence of various quantities can be understood in terms of ballistically propagating quasiparticles [55]. Recently this description has been improved further into a kind of generalized hydrodynamics [53,56,57,58] suitable to describe ballistic transport in Bethe Ansatz integrable systems. The notion of emerging eigenstate solutions [59,60] constitutes another interesting direction.
In this work we present a thorough analysis of the partitioning protocol for the transverse field Ising model, a paradigmatic integrable quantum spin chain. This system can be mapped to free fermions, which makes it possible to give a microscopic derivation of the space-time profiles in inhomogeneous quenches. First principles derivations of this kind provide a valuable testing ground for more general frameworks as the hydrodynamic approach. Somewhat surprisingly, apart from the energy current and its fluctuations in the NESS [16], not much is known about the time evolution and the approach to the NESS. Only the (transverse) magnetization profiles have been computed for domain wall initial states [36,43] and for a special two-temperature initial state at the quantum critical point [40]. The work [42] went beyond the NESS and studied the space and time dependent correlation matrix of Majorana fermions in a setup slightly different form ours.
In this paper we fill this gap by providing analytical expressions for the profiles of various quantities at or away from the critical point in both the paramagnetic and ferromagnetic phase, and for a broad class of initial states defined by the initial fermionic occupations. As a byproduct, this provides a derivation of the NESS using the standard toolkit of manybody physics. Compared to the currents and other expectation values, two-point correlation functions are much less studied (notable exceptions are [15,46,38,61,50,39]), here we discuss these as well. Note that due to the non-locality of the mapping to free fermions, the interacting nature of the model shows up in correlation functions.
We also go beyond the hydrodynamic scaling limit by computing the leading corrections to the NESS corresponding to the late time approach to the steady state. These corrections are very difficult to determine in general, and the Ising model gives us the rare opportunity to derive them analytically. This can provide insights to the applicability of (generalized) hydrodynamic approaches based on the local density approximation to describe corrections to the scaling limit.
We also study the fine structure of the propagating front which for free fermions shows universal features and a characteristic subdiffusive scaling related to the Airy kernel [62,63,39]. It is an interesting question how much of these features appear in the more complicated Ising model.
The paper is organized as follows. In Sec. 2 we review the diagonalization of the Ising chain with open boundary conditions. In Secs. 3-5 we analyze the non-equilibrium dynamics of the fermionic correlators, the building blocks of all spin correlation functions. In Sec. 3 we discuss the computation of the time evolution of fermionic correlation functions. The resulting integral expressions are analyzed in the semiclassical limit in Sec. 4 yielding the profiles of the fermionic correlators. In Sec. 5 we focus on the late time approach of the NESS and on the fine structure of the front. In Sec. 6 we use the results for the fermionic correlation functions obtained in Secs. [3][4][5] to derive the profiles of various physical quantities, including the energy density and current, the transverse magnetization and spin-spin correlation functions. The domain wall initial state is discussed in Sec. 7. We give our summary and outlook in Sec. 8. Details of the calculations and plots in support of the analytical results are collected in the appendices.

Diagonalizing the Ising spin chain
We start by giving a brief summary of the diagonalization of the transverse field Ising chain with open boundary conditions through a non-local mapping to free fermions [64,65]. Further details are given in Appendix A.
The Hamiltonian of the transverse field Ising spin chain is where σ α j are the Pauli matrices, and we assume open boundary conditions. J is the strength of the Ising coupling, and h sets the transverse field. Energies and times are measured in J and J −1 , respectively. We set J = 1/2 and work with dimensionless quantities throughout the paper.
The Hamiltonian (1) can be rewritten in terms of canonical spinless fermions c j , c † j after a Jordan-Wigner transformation: where σ ± j = (σ x j ± iσ y j )/2 are the usual ladder operators, yielding It turns out to be useful to introduce the combinations 2 in terms of which The Hamiltonian (3) can be diagonalized by a linear transformation leading to the modes {η k , η † k } : The (real) mode functions are where A k is a normalization constant such that L j=1 φ k (j) 2 = L j=1 ψ k (j) 2 = 1, and θ k is the Bogoliubov angle where we select the branch such that θ k ∈ [0, π). The variable k is quantized as In terms of the modes η k the Hamiltonian reads where the energy eigenvalues are 3 2 In the literature it is common to include a factor of i in the definition of Bj so that both Aj and Bj are Majorana operators. Here we use the original definitions of Ref. [64]. 3 Our expressions for θ k and ε k differ in minus signs from many works on the subject. This comes from a difference in the mapping (2), in particular, dropping the minus signs in the Jordan-Wigner string amounts to a π-shift in the definition of the momentum k.
We note that for h < L/(L + 1) the quantization condition (9) has a complex root corresponding to a localized edge mode which however does not affect the results below obtained in the L → ∞ limit 4 . The group velocity of these excitations is It follows that the maximal velocity is h for h < 1 and 1 for h ≥ 1, 3 Time evolution of the fermionic correlation functions The inhomogeneous initial state consists of two independent, disjoint chains of length L having macroscopically different properties, e.g. thermalized at different temperatures T L and T R . So the initial density matrix is a tensor product, At time t = 0 the two chains are joined and let evolve by the Hamiltonian H of the chain of length 2L. In other words, we turn on the coupling between site 0 and site 1, where H L is the Hamiltonian of the left chain of sites j = −L+1, . . . , 0, H R is the Hamiltonian of the right chain of sites j = 1, . . . , L. This setup is sometimes referred to as the "two-reservoir initial state" or the "partitioning protocol". In Eq. (15) we introduced the fermionic mode operators γ k that diagonalize H on the chain of length 2L. Let us denote the mode functions of H R by φ R q and ψ R q , then due to swapping the boundary conditions, the mode functions of H L are given by The mode functions of the final Hamiltonian on the full chain H are ϕ k (j) = φ k (j + L) and χ k (j) = ψ k (j + L), where the momenta {k n } are quantized in volume 2L, which amounts to the substitution L → 2L in Eq. (9). Note that the functional form of ε k and θ k are unchanged as we do not perform a quench in the Ising interaction or the transverse field. The fermion operators A j , B j on the full chain are related to the modes by In order to compute the time evolution of spin correlation functions, we first need to compute the building blocks given by the correlations The time evolved operators in the Heisenberg picture are obtained by writing them in terms of γ k using (17), exploiting the simple time dependence of the mode operators and then rewriting everything in terms of A j , B j . After these steps we find [36] In the infinite volume limit, L → ∞, the sum over k can be written as an integral. After dropping highly oscillating terms in the integrands 5 we find where we defined the infinite volume mode functions Using Eqs. (18), the time-evolved two-point functions can be written as where X and Y stand for A or B. The initial correlations will be non-zero only if j, l ≥ 1 or j, l ≤ 0, so the correlation function (22) can be split into two parts corresponding to the contributions of the left and right chains. We shall denote these contributions by X n (t)Y m (t) L/R . The initial correlations are calculated by rewriting the fermion operators in terms of the mode operators η q that diagonalize the left and right half chains. Here we assume that for each half chain the anomalous correlations are zero, so For the two-temperature setup is the thermal Fermi-Dirac distribution function. Exploiting the completeness of the mode functions, we arrive at Due to the orthonormality of the mode functions, the terms containing the correlations (25a) yield a Kronecker δ n,m in Eq. (22), resulting in Since the two half-chains initially differ solely in the initial fermionic occupations, the left/right parts of the fermionic correlation functions are related to each other by (for a proof, see Appendix B) Expression (26) can be used to numerically evaluate all equal-time correlation functions on a finite chain by using the finite volume mode functions (7) as well as the edge modes in the initial correlations (25) and in the coefficients (19). All our numerical results presented in the forthcoming sections below were obtained in this way.
Equation (26) is also the starting point of analytic manipulations in the infinite system size limit. Then the coefficients in Eqs. (20) are used and the sums in Eqs. (25) are straightforwardly turned into integrals. The remaining sums over lattice sites run from −∞ to ∞, eliminating all explicit dependence on L. The correlation functions are thus given by a double sum of triple integrals. The sums are essentially geometric series and can be performed analytically. The integral over q coming from expressions (25) can be manipulated using contour integral techniques. For details of the derivation we refer the reader to Appendix B. The final result is a set of double integral expression that give all the equal-time correlations exactly at any time after the quench: Here (ν = L,R) where I ν (k) are even functions, expressed as contour integrals in Eq. (158) of Appendix B. The B n B m correlations can be obtained using relations (28). These expressions were not known and thus are new results of the paper. They can be used to obtain further analytic results in various limits which is the subject of the next section.

The semiclassical limit
The existence of a scaling limit for the magnetization profile evolving from a domain wall initial state in the XX spin chain was first observed in Ref. [34]. The same behavior was found for the two-temperature initial state [12] and for the Ising chain [36,42].
This semiclassical or hydrodynamic limit corresponds to the case when n, m, t → ∞ with lim n/t = lim m/t fixed implying (n − m)/t → 0. Physically, we focus on the large time behavior near a given ray defined by the ratio n/t ≈ m/t. The appearance of this scaling variable suggests a hydrodynamic description in terms of quasiparticles [55].
The behavior in this limit can be derived through a stationary phase analysis of the double integral expressions [12,39]. In each integral, the stationary points are given by As (n − m)/t → 0, the stationary points coalesce, k s → k s , and the integrand becomes singular. This singularity governs the leading order behavior in the semiclassical limit. We note that another singularity could arise as k s → −k s , however, the integrand vanishes due to G ν (k, −k) = 0 so this gives a subleading contribution. Let us start with the right contribution to the correlation (29) given by the first line of the right hand side. We first introduce the center of mass and relative coordinates along with the momentum difference and total momentum: Now we expand around Q = 0 to find Using the integral expression for the Heaviside theta function, In the semiclassical limit x, t → ∞, so the constant shift in the argument of the step functions can be neglected. The left contribution can be obtained from Eqs. (28) via the substitution The total result is then where u = x/t = (n + m)/(2t) specifies the ray. In a similar fashion, from (30) we obtain where we used that the integral of sin(Kr) vanishes.
For a fixed separation r, the resulting expressions (36), (37) are functions of u = x/t, i.e. they depend on the ray only, which is the hallmark of ballistic behavior. Indeed, a natural interpretation of these expressions can be given in terms of free quasiparticles [55] with dispersion relation ε k that are present on each side before the quench following the original left and right distributions. Around a given space-time point (x, t) only those left particles can contribute which have velocities larger than u = x/t. Similarly, contributing right particles cannot have velocities greater than u. This is perfectly reflected by the Heaviside theta functions of the integrands. The same observation was made in [42] for translationally invariant Gaussian initial states evolving under a Hamiltonian with a localized defect. As we will see below, this interpretation is even clearer for simple physical quantities such as the energy density and energy current.
If u > v max = min(1, h), or u < −v max , then one of the theta functions in the integrands is equal to 1 while the other vanishes for all K, so only the right or the left contribution survives. This is a horizon effect: the physics outside of a light cone set by the maximal group velocity is unaffected by the quench of joining the two sides. If the initial half chains were in equilibrium before the quench then all observables assume their initial expectation values and they change only when the ballistically propagating front arrives 6 .
For asymptotically large times, a non-equilibrium steady state (NESS) is formed around the junction in the middle of the chain. In accordance with the quasiparticle picture, the region in which the system is asymptotically close to the NESS grows ballisitically. The NESS is obtained in the limit t → ∞ and n, m fixed which amounts to u = 0, implying and As the NESS is translationally invariant, the correlations in the NESS depend only on the separation r = m − n. The NESS correlations can be interpreted in terms of quasiparticles with the left momentum distribution (e.g. thermalized with the left temperature) going to the right and by quasiparticles with the right momentum distribution (e.g. thermalized with the right temperature) going to the left. These expressions were first derived using C * -algebra methods in [14,66]. In our approach, the NESS appears as a single special member of a continuous family that gives the spacetime profile of correlations in the semiclassical limit. The same was achieved in [42] for a homogeneous initial state evolving under a non-homogeneous Hamiltonian.
Finally, let us write down the steady state correlations of the original fermion operators:

Beyond the semiclassical limit
In this section we discuss the leading finite time corrections to the NESS and the structure of the propagating front.

Approach to the NESS
Although the correlations in the NESS were known, the finite time corrections characterizing the late time approach to the NESS have not been analyzed before. The analytic expressions (29), (30) can be used to study the leading corrections in the limit t → ∞, n/t, m/t → 0. The details of the calculations are given in Appendix D, here we only quote the results.
Note that in the critical case (h = 1) the result is much simpler, where we used that θ 0 = 1/2. The correction to the left part A n B m L NESS can be obtained from Eq. (42) via the mapping (28). These expressions are checked by comparing them to numerical results in Fig. 7 in Appendix D.
The first correction in the second and third line is oscillating around the NESS value both in time and space. The wavelength of the spatial oscillation is 2 (lattice spacings) while the frequency of the temporal oscillation is either 2 or 2h. Some of the terms in the last two lines can be written as where for the derivatives of g(k) and θ k we used Eqs. (174),(187) and that θ(π) = 0 for h > 1 and θ(π) = π for h < 1. This means that the profile of A n B m for fixed time and separation m − n is linear near the origin. Since (m + n)/t = 2u, this is the first order correction in u to the NESS given by u = 0. Indeed, (44) can be obtained by a Taylor expansion of the semiclassical result (36) in u (see Appendix D). On top of this there is an additional correction in the last two lines of (42): This amounts to a finite shift at any large but finite t, in particular, in the origin the NESS is not reached immediately but only asymptotically as ∼ 1/t. All these different types of corrections can be observed in Figs. 1b,1d where we plotted the energy density, a combination of A n B m correlators (c.f. Eq. (63)). With time, the flat profile of the translationally invariant NESS develops: the oscillations get damped, and both the slope and the offset of the profile near the origin tends to zero.
In a similar fashion, we obtain The corresponding expressions for B n (t)B m (t) can be determined using the mapping (28). These corrections are compared with the numerical data in Fig. 7. In contrast to A n B m , the correlation functions A n A m and B n B m are constant around the origin and are not shifted with respect to their NESS value. The leading order corrections are given by oscillations whose amplitude decays as ∼ 1/t. However, note that if n − m is even, these terms vanish and the leading correction is of order O(t −2 ). This is also true in the critical case h = 1 for any n, m. We note that the leading order corrections in Eqs. (42), (46) show a parity effect for h = 1, i.e. oscillations with wavelength of two lattice sites. This may cast doubt on the ability of the local density approximation and a hydrodynamic description to capture the corrections to the ballistic semiclassical results. 7 Finally, we mention that a very similar derivation can be performed in principle at finite u = x/t values to obtain the leading late time corrections to the semiclassical profile. The essential difference is that the stationary points would be shifted. We have not carried out this calculation but performed numerical checks which show that the leading corrections are O(1/t) also in this case.

Fine structure of the front
It is an interesting question whether the fine structure of the propagating front shows some universal behavior. In Ref. [62] the XX spin chain evolving from a domain wall initial state was studied, and it was found that the shape of the front can be described by a scaling function. In particular, the width of the front scales sub-diffusively, as t 1/3 . In Ref. [63] the full distribution of the magnetization was determined and it was related to the eigenvalue statistics of random matrices. The rescaled correlations were shown to be described by the Airy kernel (see also [39]). In Ref. [40] a similar behavior was observed numerically for the two-temperature setup in the Ising spin chain for h > 1 and T L = 0, T R = ∞.
Here we analyze this edge behavior within our formalism for any h and any left/right temperatures, following the method of Ref. [39]. The fine structure of the front is given by a correction to the time-independent equilibrium value in the region that has not yet been reached by the front. This correction can be studied by taking the limit x, t → ∞, x/t, y/t ≈ v max , (x − y)/t → 0. In this limit the stationary points for x and y approach each other giving rise to a specific scaling form of the edge of the front.
We present the analysis of A n (t)B m (t) R near the right front. This is the most complicated case, as the right equilibrium value of this correlation is non-zero (The left contributions are zero outside of the light cone and the A n A m , B n B m correlations are zero for n = m.) Clearly, this value must be subtracted before we can focus on the correction. The starting point is the first double integral expression in the general formula (29). We show that the (right) equilibrium correlation is related to the residue of the integrand at the pole k = k .
which coincides with the semiclassical result for u = (n + m)/(2t) > v max that gives the equilibrium expectation value to the right of the front.
In the double integrals in (29) the infinitesimal imaginary shifts in the denominator mean that the integration contour of the k integral avoids the pole at k = k from below. By virtue of the residue theorem, the expression (47) is the result of performing the k integral along the contour encircling the point k on the real axis. Subtracting this value is equivalent to pulling the contour through the pole at k = k so that it runs above the pole. This amounts to changing the sign of δ in the double integral. Then can be recast as (note the sign of δ) The stationary phase for the term in the first bracket is given by ε k = n/t, ε q = −m/t or ε k = −n/t, ε q = m/t. As we are interested in the front, n and m are close to each other so k ≈ −q follows. But then both ε k −ε q ≈ 0 and G(k, q) ≈ 0, so this gives a subleading correction. The leading term is given by the second bracket and the stationary points ε k = n/t, ε q = m/t, and ε k = −n/t, ε q = −m/t. Let us introduce the wave number κ corresponding to the maximal group velocity, Then expanding around ±κ we have where we used that ε κ = 0 and that ε k is even in k, and analogous expressions hold for ε q t ± qm. We have to treat differently the critical and off-critical cases.

Off-critical case
For h = 1 the leading contribution can be written as where we definedk = k − κ,q = q − κ in the first term andk = k + κ,q = q + κ in the second term. As only the smallk,q region gives non-negligible contribution to the integral, the (k +q) cot κ term can be dropped in the coefficient of iδ. We introduce rescaled variables, and we extend the domain of integration to the real line since the large K, Q region does not contribute due to the fast oscillations. where is the celebrated Airy kernel [67].
Following the same steps we find The functional form of the left contributions for all correlation functions turns out to be simply opposite in sign (and obviously, f R κ is replaced by f L κ ). The appearance of the Airy kernel is quite generic in free fermionic systems near an edge originating either from a moving front as in our case, or from an external confinement [68,69,70]. The mathematical reason is the presence of a degenerate stationary point of the phase in the integrand of the fermionic propagator where the denominator also becomes singular [71]. Even though the integrands in question are slightly more complicated for the Ising chain, the necessary ingredients for the appearance of the Airy kernel are present. The finite temperature only enters in the prefactor of the Airy kernel.
However, not all correlations are described by the Airy kernel. Interestingly, in the ferromagnetic case (h < 1), cos θ κ = 0 as can be seen from Eqs. (49) and (8). This implies that the right hand side of (54) vanishes for n = m and the edge behavior of A n B n is not described by the Airy kernel. As this correlator is related to the density of Jordan-Wigner fermions, this behavior is markedly different from the case of free spinless fermions. This difference may be related to the fact that the fermion number is not conserved in the Ising spin chain.
We note that the measurable excess correlations with respect to the equilibrium values are not identically zero outside the classical light cone. Using the asymptotic behavior of the Airy kernel for large positive arguments, we see that their decay is faster than exponential, in accordance with Lieb-Robinson bounds [72].
Corrections to the scaling forms (54), (56) can be obtained by expanding the non-singular part of the integrand to first order ink andq. For example, in Eq. (48) we expand The extra terms linear ink andq can be written as derivatives with respect to n and m, and eventually lead to the appearance of the partial derivatives of the Airy kernel: The derivatives of F (k, q) are computed using Eqs. (31) and (173). We note that even these corrections vanish for A n B n if h < 1.
Both the scaling forms (54), (56) and the results including these corrections are compared with the numerical data for finite t in Fig. 8 of Appendix E. It is clear that the agreement with the numerical results is improved significantly by the leading order corrections. As we shall see below, they are also responsible for interesting physical phenomena.

Critical case
For h = 1 the stationary point is at κ = −π. The derivation has to be modified because of the appearance of formally singular terms and because the stationary point is now located at the edge of the Brillouin zone [75]. We shall need the values Due to the last equality G L/R = 0, so we have to Taylor expand it alongside with ε k . Moreover, in this case the the second line in Eq. (48) is of the same order as the third line and cannot be dropped. The reason is that both ε k ± ε q are zero at k = q = κ, and out of the otherwise rapidly oscillating factors in the second line e i(ε k +εq)t is simply absent for k = q = κ while e i(kn+qm) = e −iπ(n+m) = e −iπ(n−m) is the same as the analogous factor in the third line. Thus all the four terms in (48) corresponding to the four possibilities k = ±π, q = ±π must be taken into account. Because ±π are the edges of the Brillouin zone, after the rescaling (52) the integrals over K and Q are on the positive or negative axis only, again realizing the four possibilities. Luckily, the integrands turn out to be the same in the four terms and their sum can be rewritten as a double integral over the real axis, as in the off-critical case.
Using ∂ 1,2 G(k, q) = g (∓k) from Eq. (31) and g (π) = g (−π), the Taylor expansion in all the four terms gives where the + sign is for the (k, q) ≈ (∓π, ±π) terms and the − sign is for the (k, q) ≈ (∓π, ∓π) terms. From here the derivation follows the off-critical case. The surviving (k +q) factor translates into derivatives just like in the off-critical corrections, and using g (π) = −1/(2T ) we finally obtain where X = (n−t)( while the left contributions are opposite in sign. We found that in the critical case the front is not described by the Airy kernel but by its derivative, and the decay is ∼ t −2/3 instead of ∼ t −1/3 . Moreover, the front in this case is featureless and lacks the typical staircase structure of the Airy kernel (see Fig. 2b).

Physical observables in the semiclassical limit
In this section we use the general expressions derived in the previous section to calculate the space and time dependence of various physical quantities, including the energy density, energy current, and various correlation functions. The analytical expressions are compared with numerical results obtained by computing the sums in Eq. (26) for finite systems in the the two-temperature case 8 .

Energy density
The energy density in the bulk can be written as (c.f. Eq. (5)) Using the semiclassical result in Eq. (38) we find where we used the identity (c.f. Eq. (129a) in Appendix A) The constant term −1/2 π −π dk/(2π)ε k corresponds to the ground state (zero point) energy. The rest of the formula has an obvious semiclassical interpretation: in the vicinity of the space-time point (x, t) the energy density is given by the particles arriving at point x at time t coming either from the left or the right with distributions corresponding to the left or right reservoir. In the NESS, The semiclassical formula (64) is compared with the numerical results at different values of h and temperatures in Fig. 1. Expression (64), plotted in dashed line, is seen to capture accurately the actual energy density profiles (dots and solid lines). For lower temperatures, the deviations from the semiclassical result are relatively larger. These corrections are however captured by our analytic results given in Sec. 5.1. The approach to the NESS near the origin is well described by the analytic expression (42), see Fig. 1d. The oscillations, the finite slope and the finite shift are all clearly visible. The numerical data show a longer wavelength modulation as well which is a subleading correction beyond our result. Near the front the oscillations are enhanced, which is described by the results in Sec. 5.2. Using Eq. (54) in the expression (63) for the energy density and invoking relation (65), we find the scaling behavior for h = 1 near the right edge where h R is the bulk energy density on the right and we used that in the large time limit This is plotted in Fig. 2a (dashed line) together with numerical results for h = 0.2 (dots). For finite time, the corrections (59) involving the derivatives of the Airy kernel and taking into account that X = Y, shown in solid line in the plots, yield significantly better agreement. It turns out that these corrections are responsible for the sizeable oscillations in the low temperature case (see Fig. 2b).
In the critical case from Eq. (62) we obtain with X = (n − t)(8/t) 1/3 , Y = (n + 1 − t)(8/t) 1/3 . This is shown in Fig. 2b together with the numerical results. The front is featureless in the critical case and does not have the staircase structure of the Airy kernel.

Energy current
The energy current can be obtained by taking the commutator of the Hamiltonian with the energy density and writing it as a divergence. This yields Using Eqs. (37), we obtain the current in the semiclassical limit, where we used that v k = −h sin k/ε k . As v k ε k is the semiclassical contribution of a quasiparticle to the energy current, this result also follows from a semiclassical picture similarly to the energy density: at a given space-time point the energy current is a sum of energy currents carried by single particles that arrive at time t at point x. On both sides, the current is zero outside the light cone set by the largest velocity v max . It is stratightforward to show that the semiclassical expressions for the energy density and the energy current satisfy the continuity equation, The current in the stationary steady state is in agreement with the results of Ref. [16]. At the quantum critical point we recover the conformal field theory result [20].  Typical profiles for different values of h and left/right temperatures are shown in Fig. 3. As the right temperature is higher, the current is negative so we plot −J e n instead. In each case, the semiclassical expression (70) (dashed line) agrees well with the numerical results (solid line). In Fig. 3d, the leading order correction based on Eqs. (46) is shown to capture the oscillations near the junction. The amplitudes of the oscillations are relatively large when one of the temperatures is low (c.f. Fig. 3b), so the corrections to the semiclassical limit are more significant. The oscillations are especially pronounced near the edges. Here a surprising phenomenon can also be observed on the hotter side: right after the arrival of the front, the current is initially negative, i.e. it starts to flow in the "wrong" direction, from the colder system to the hot reservoir.
Similarly to the case of the energy density, the edge behavior of the current can be obtained from the results of Sec. 5.2. Combining Eqs. (69) with (56) we find where we used v max = −h sin κ/ε κ . This result is shown in dashed line in Fig. 4 where it is compared with the numerical data shown in dots. We also plot in solid line the improved result including the first corrections to the scaling form, similar to Eq. (59). We conclude that the corrections to the scaling form (73) are responsible for the large oscillations, and in particular, for the energy backflow phenomenon.
In the critical case, with X = (n − t)(8/t) 1/3 , Y = (n + 1 − t)(8/t) 1/3 , and the oscillations are absent. The edge behavior of the energy current was first studied in [75] but the result obtained there seems to be different from ours.

Transverse magnetization
The transverse magnetization operator is simply related to the density of the Jordan-Wigner fermions: σ z n = 2σ + n σ − n − 1 = 2c † n c n − 1 = −A n B n . Its expectation value in the semiclassical limit is given by (36) for n = m : The cos θ k factor accounts for the Bogoliubov rotation connecting the original fermions with the modes diagonalizing the Hamiltonian. Indeed, using Eqs. (6) and (7), it is straightforward to show that in the L → ∞ limit, in a translationally invariant state In the semiclassical picture this correlation is transported through the system by the quasiparticles, which immediately leads to Eq. (75). The magnetization profile is shown in Fig. 5a for h = 3, T L = 0.5, T R = 2, for other parameters the profile is similar. The agreement between the semiclassical result (75) and the numerical data is excellent. The front of the transverse magnetization is described by Eq. (54) and the analogous left contribution: In the light of this result we can revisit the case studied in Ref. [40]. If T L = 0, T R = ∞, then f L = 0 and f R = 1/2 constants, and we arrive at σ z n (t) ∼ cos θ κ (v max t/2) −1/3 K(X, X). Using the asymptotic expansion of the Airy kernel for large negative arguments, K(X, X) ≈ √ −X/π, we recover the square root dependence observed numerically in Ref. [40].
Finally we note that in this simple case (when T L = 0, T R = ∞ and f L k = 0, In the critical case ε k = 2 cos(k/2), v k = − sin(k/2), θ k = k/2, and with u = n/t, so the profile is linear within the light cone, for −1 ≤ u ≤ 1 [40]. In this special case the energy density and the current also take simple forms: We note that without the zero point energy the energy density is h n (t) sc = (u + 1)/π.

Transverse spin-spin correlations
The transverse spin-spin correlation function is also easy to compute due to the fact that the relation between σ z and the fermion creation/annihilation operator is local. Using Wick's theorem we obtain It is straightforward, albeit lengthy, to write down the semiclassical limit using Eqs. (36), (37). Comparison with numerical data for the nearest neighbor correlator is plotted in Fig. 5b-d.
Interestingly, we find that the profile can be very non-monotonic, and the correlator can take values that are very far from the initial equilibrium values on the two half-chains, and can even change sign.

Longitudinal spin-spin correlation function
In contrast to the transverse magnetization, the longitudinal one is non-local in terms of the fermionic degrees of freedom that diagonalize the Hamiltonian. Consequently, correlation functions involving σ x n are much harder to calculate. In the equal time two-point function, the Jordan-Wigner string survives only between the two operators so it can be written as Using Wick's theorem, the multi-point function on the right hand side can be written as the Pfaffian of the matrix where Γ j,l are two-by-two matrices given by In the semiclassical limit, translational invariance is restored around a given ray. Consequently, Γ becomes a block Toeplitz matrix with blocks where so the block symbol isΓ It is of interest to obtain the asymptotic behavior of the correlation function for large separation of the two operators in the semiclassical limit, i.e. when m−n 1 but (m−n)/t → 0. First we note that the Pfaffian of an antisymmetric matrix is, up to a sign, the square root of its determinant. The asymptotic behavior of the determinant of a × block Toeplitz matrix is given by a generalization of the Szegő lemma [73], provided that the matrix satisfies certain properties. Even in thermal equilibrium the Szegő lemma can be directly applied only in the ferromagnetic phase, h < 1 [74], so we consider this case only. Generously forgetting about these conditions, a naive application of the lemma gives This implies that for large separations the correlation function decays exponentially [15], where the correlation length is The argument of the logarithm is a long expression, but the integral can be simplified drastically (see Appendix F) and can be shown to be equal to By setting f L k = f R k we reproduce the space and time independent equilibrium correlation length ξ −1 = − π −π dk 2π ln(1 − 2f k ). The inverse correlation length is given by the sum of contributions from particles with left and right momentum distributions. The NESS value, obtained by setting u = 0, is simply the average of the left and right inverse correlation lengths.
In Fig. 6 we compare the result of the Szegő lemma, Eqs. (90) and (92), with the profile obtained numerically for two different times. The prefactor of the exponential decay in (90) is set to be the square of the ground state magnetization,σ 2 = (1 − h 2 ) 1/8 . The agreement is quite satisfactory despite the fact that (90) is an asymptotic result while we plot it for a separation of a few lattice spacings.

Domain wall initial state
Even though in the previous section we showed results for the two-temperature scenario, our analytical results apply to a broader class of initial states, namely, to states characterized by fermionic occupations with no anomalous correlations, c.f. Eq. (23). The explicit form of f (k) was only used in the evaluation of some terms appearing in the corrections to the NESS (see Eqs. (172),(177) in Appendix C), but even the formulas (42), (54), (56) remain valid.
As an example, we consider the domain wall or kink initial state where initially the two half chains are fully polarized in the transverse direction: σ z j = +1 on the left and σ z j = −1 on the right. An important difference with respect to the two-temperature case is that these initial states are not eigenstates of H L/R , so they have non-trivial time evolution even without joining the two systems. As the semiclassical limit corresponds to t → ∞, the regions outside the light cone have already relaxed by the time the front reaches them. In other words, the scaling functions for the profiles interpolate between the bulk equilibrium values and not between the initial ones.
Our first task is to compute the initial correlations in momentum space. Using Eqs. (6), (94) Since A j and B j flip spins in the z-basis, it is easy to see that where the upper sign is for the all spins up and the lower sign is for the all spins down state.
Using the explicit expressions for the mode functions (7), we obtain in the L → ∞ limit Plugging this in Eq. (94) we find so the anomalous correlations are zero and 1 − 2f k = ± cos θ k . Using this in Eq. (64) for the energy density we find Recalling that ε k cos θ k = h + cos k we arrive at given by Then v k < u implies −π ≤ k < k 1 or k 2 < k ≤ π if u > 0, and k 2 < k ≤ k 1 if u < 0. Then the energy density can be written in the closed form In the critical case ε k = 2 cos(k/2), v k = − sin(k/2), θ k = k/2, so k 1 = −sgn(u)π, k 2 = −sgn(u) arccos(1 − 2u 2 ) = −2 arcsin u, which gives Similarly, from Eq. (70) a result first obtained in [36]. The magnetization profile interpolates between the equilibrium values ±1/2 which are half of the initial magnetizations on the two sides.

Conclusion
In this paper we studied the out of equilibrium evolution of the transverse field Ising model in the partitioning protocol when the initial state is the tensor product of two states with different macroscopic properties. We focused on the case when initially the anomalous correlations vanish (c.f. Eqs. (23)). We presented several analytical results and performed numerical simulations. Our results are directly applicable to initial states with vanishing anomalous fermionic correlations.
We derived analytic double integral expressions (29), (30) for the fermionic correlators valid for any time and for arbitrary values of the transverse field. Based on these expressions, we obtained the profiles of the fermionic correlators in the space-time scaling limit and showed that the correlations are given by semiclassical expressions (36), (37). As a by-product, we presented a down-to-earth derivation of the correlations in the non-equilibrium steady state.
We also derived the finite time corrections to this current carrying steady state and found a ∼ 1/t asymptotic approach for most correlations, although some of them have a ∼ 1/t 2 approach (see Eqs. (42), (46)). Interestingly, we found that the leading finite time corrections display lattice effects in the off-critical case, which may imply that hydrodynamic approaches based on the local density approximation are not able to capture the leading corrections to the scaling limit, and in particular, they cannot account for sub-ballistic (e.g. diffusive) corrections.
We also investigated the properties of the ballistically propagating front and found that in most off-critical cases it can be described in terms of the Airy kernel. This provides evidence for the universality of the Airy kernel not only for edge phenomena in free fermion systems but also for spin systems that can be mapped to free fermions. It would be interesting to study interacting integrable models in this respect. Comparison with numerical simulations showed the necessity to determine subleading corrections to the universal Airy kernel behavior.
However, perhaps surprisingly and in contrast to the case of free spinless fermions, in the critical case none of the correlations are described the Airy kernel near the edges (see Sec. 5.2.2). The edge behavior is captured by a different kernel given by the derivative of the Airy kernel (see also Ref. [75]). In particular, it is featureless and does not display the typical staircase structure (c.f. Fig. 2b). Moreover, the fermionic density or transverse magnetization is not described by the Airy kernel in the ferromagnetic phase. These findings may be related to the fact that the fermion number is not conserved in the transverse Ising chain.
Based on these general results for the fermionic building blocks, in Section 6 we investigated various physical quantities of the spin chain, including the energy density, the energy current, the transverse magnetization, and two-point correlation functions of the longitudinal and transverse components of the spin. We found that the energy current can show a back-flow phenomenon which means that it can flow in the wrong direction (c.f. Figs. 3b and 4a). This effect is not accounted for by the universal scaling form of the front in terms of the Airy kernel but it is described by the corrections to it involving the derivatives of the kernel. Another interesting finding is the strongly non-monotonic behavior of the transverse spin-spin correlation function (c.f. Fig. 5c).
As for the longitudinal correlation function, we derived a non-rigorous expression (92) for the correlation length based on the naive application of the Szegő lemma. Comparing with numerics, we found very good agreement even for short separations (see Fig. 6).
In the future it would be interesting to study the evolution of dynamical (two-time) spin-spin correlation functions in this and similar setups. The methods of the paper can also be applied to more general inhomogeneous initial states and to situations that involve a global quench as well. As a simple generalization, initial states with non-zero anomalous correlations, η k η k = 0, could be considered. The techniques of this work can be naturally carried over to the XY spin chain. The space and time dependent fluctuations and the large deviation functions of different quantities can also be computed, and it would be interesting to check the corresponding extended fluctuation relations [27].
Note added Recently, the work [75] appeared in which some of the results in the space-time scaling limit presented in our paper are derived using a different approach.
Funding information The author acknowledges funding from the "Bolyai" Scholarship and a "Prémium" Postdoctoral Grant of the Hungarian Academy of Sciences. He was partially supported by NKFIH grant no. K 119204.

A Diagonalizing the Ising chain with open boundary conditions
In this Appendix we provide, in order to be self-contained, details of the diagonalization of the Hamiltonian (3) based on [64,65]. The Hamiltonian is a bilinear form with the real matrices A and B given by (108) It follows that the Hamiltonian can be diagonalized by a linear transformation leading to the modes Here we anticipated that g k (j), h k (j) are real. In terms of these modes the Hamiltonian reads If the η k are canonical fermionic operators (see below), then Substituting Eqs. (107),(109) in this relation, we arrive at Let us introduce the combinations in terms of which Eqs. (112) read implying The product matrices (A ± B)(A ∓ B) are real and positive semi-definite, so all the eigenvalues ε k are real and it is possible to choose φ k and ψ k to be real and separately orthogonal. If, say, φ k is orthonormal, then ψ k will be orthonormal since It is a simple exercise to show that due to the orthogonality of the eigenvectors, hold, so the transformation (109) is indeed canonical. Using the completeness of the eigenvectors, the inverse transformation can be written as The constant in (110) can be determined exploiting the fact that under a canonical transformation the trace of H is invariant. 9 It will turn out to be very useful to introduce the combinations in terms of which The Jordan-Wigner string is a product of −σ z j so σ x j , σ y j can also be expressed as a product of A j and B j operators. Their relation to the mode operators is and The matrix equations (114) imply and These can be translated to with boundary conditions Plugging in the Ansatz ψ k (j) = −e ikj , (127b) 9 The notation is hiding the fact that the problem originally involved 2L×2L matrices acting on vectors with entries cj and c † j . This is reflected by the fact that we have two sets of L basis vectors. In principle we need 2L vectors to diagonalize the Hamiltonian. But notice that Eqs. (114) have the symmetry ε k → −ε k , ψ k → −ψ k . This amounts to interchanging g k and h k or η k with η † k , thus it is a particle-hole transformation, which renders the state with all negative energy states filled the new vaccum.
we obtain or taking the real and imaginary parts, which gives for the Bogoliubov angle tan θ k = sin k h + cos k .
As we will see, k ∈ [0, π], so if we set θ k ∈ [0, π] then ε k ≥ 0. Moreover, sin θ k > 0 and tan θ k , cos θ k are negative when h + cos k < 0. For the energy eigenvalues we get Imposing the boundary condition for ψ k we find Imposing the boundary condition for φ k leads to the quantization condition or k = π L + 1 n + 1 π arctan sin k h + cos k .
Note that this means that The normalization constant is Let us collect some useful relations: so ε k = sin k/ sin θ k and the quantization condition can be written as For h > 1 there are always L real solutions, while for h < 1 it is possible to have a complex root corresponding to a localized mode. To see this let us consider the last quantum number, n = L. The value k = π − 0 is always a solution (remember that 0 < θ k < π,) but it is unphysical. Expanding Eq. (9) around k = π we find the solution This is always satisfied for h > 1, but for h < 1 it is true only if L < h/(1 − h). So for a large enough L at fixed h there will always be a complex solution of the form k 0 = π + iv with the imaginary part satisfying For L large, v quickly approaches − ln h. Expanding around it we find The energy eigenvalue is exponentially small: The corresponding eigenvectors are which shows that these modes are localized at the edges of the chain. The normalization factor behaves for large L as A k 0 ≈ 2 √ 1 − h 2 h L . At criticality h = 1, θ k = k/2, and the quantization condition becomes The eigenfunctions are with their energy being ε k = 2 cos k 2 .

B Exact time evolved fermionic correlations
In this Appendix we explain how to arrive at the exact double integral expressions for the fermionic correlation functions. First note that the contributions of the left and right half-chains can be related to each other. For example, using Eq. (16), where we changed variables toj = 1 − l,l = 1 − j, and used X 1−n |Y 1−m (t) = Y n |X m (t) and A j |A n (t) = B j |B n (t) which follow from Eqs. (20).
We now turn to the manipulation of the correlation functions. Our starting point is the exact expression (26) with (27). Substituting the coefficients from Eqs. (20) and the initial correlations from Eqs. (25) we obtain Here we introduced They can be interpreted as time evolved eigenfunctions with being the overlaps between the eigenmodes of the half systems and those of the full system. Explicitly, Let us list all the products of the time dependent eigenfunctions that appear in Eqs. (149): Interchanging the integral over q with the integrals over k, k in Eqs. (149) we see that for all correlation functions we need to evaluate the integral where ν stands for R,L and we introduced Changing variables to z = e iq , the integral becomes a contour integral on the unit circle going around the origin: .
Now we separate the contributions of the simple pole at z = e ik−δ : where is a contour integral around a circle C of radius min(h, h −1 ) < r < 1. We will not evaluate the integral, only note that it is an even function of k. Now we are in the position to write down fermionic correlations. Introducing the notation the contribution of the right half to the correlators are where we made the variable change k → −k , k → −k for the second exponential in the numerator. Similarly, and The final result is obtained by adding the left and right contributions: It is easy to check that relations (28) hold: C Some properties of the g(k) and I(k) functions Here we analyze the integral expressions (155) and (158) and their derivatives at some special points. The resulting expressions can be used in the computation of the corrections to the semiclassical results (see Appendix D). The first combination we consider is Since there are no poles any more at |z| = 1, we can rewrite this as (z = e iq ) We note that for f = 0 (T = 0) the integral can be evaluated in closed form: where K(m) is the complete elliptic function of the first kind. If f (ε) = 0, the analytic properties of the integrand in (167) depends on the function f (ε). For thermal distributions f (ε) = 1/(1 + e ε/T ), then 1 − 2f (ε) = tanh[ε/(2T )] which is an odd function of ε and this eliminates the branch cut coming from the denominator. In return the tanh function has poles on the real axis where ε(z)/(2T ) = ikπ/2 where k is odd. The relevant solution of this equation is Expanding the integrand around these points we find so using the residue theorem we have to add the contributions of these infinitely many poles: Let us now investigate the derivative of g(k). Differentiating the integrand of (155) we obtain This expression can be used to compute numerically g (k) for all k except k = 0, π. At these points the integrand is not well-behaved (the cosine in the denominator becomes real). We need to get rid of the singularities, i.e. we need to start from Eq. (157). Here the integrand of I (k) contains a sin k factor (without iδ) that vanishes at k = 0, π. We are left with the derivative of the first term in Eq. (157) which gives We will also need the second derivatives at k = 0 and k = π. The second derivative of the first term in (157) vanishes at these momenta due to ε (0) = ε (π) = 0, so we are left with the second derivatives of I(k). Differentiating the integrand in (158), Here the integrand still has a pole at z = ±1, so we cannot extend the contour to the unit circle in order to reintroduce q. For the Fermi-Dirac distribution f (ε) = 1/(1 + e ε/T ) we can again add the residues of the infinitely many simple poles of the tanh function to obtain

D Approach to the NESS
In this appendix we deal with the finite time corrections to the NESS, or in other words, with the late time approach of the steady state. As a warm-up, let us focus on the profiles of the fermionic correlators near the origin. The Heaviside theta functions in Eqs. (36), (37) restrict the integrals on finite intervals set by k 1 and k 2 defined in Eq. (102). The NESS is given by u = 0 when k 1 = −π, k 2 = 0. For small u Then where we used that θ 0 = 0 and θ π = π for h < 1 and θ π = 0 for h > 1. The leading order correction gives a linear behavior near u = 0. For the correlators A n (t)A m (t) and B n (t)B m (t) the integrand is sin(rK) which vanishes both at K = 0 and K = π thus the first order correction vanishes. Now we analyze the behavior of the correlators in the limit t → ∞, n/t, m/t → 0. Following the method of Ref. [39], we take the derivative of our expressions with respect to time which results in the cancellation of the singular denominator.
where osc./t 2 stands for terms of order t −2 which are oscillating. After integration with respect to t, these will remain in the same order giving subleading contribution 10 . However, nonoscillating O(t −2 ) terms need to be separated as they will become O(t −1 ) after integration.
Substituting the various integrands for the function F (k), expanding the products and collecting the terms we find, using θ(0) = 0, Unlike A n B m , in the t −2 order there are only oscillating terms that are subleading after integration. Now we can integrate the result to obtain Eq. (46a) of the main text. Formula (46b) for A n (t)A m (t) L is obtained in an analogous fashion. Using the relation Eq. (28) we readily obtain the corresponding expressions for B n (t)B m (t) . These results are also compared with numerics in Fig. 7.

E Edge behavior of the fermionic correlations
In Sec. 5.2 we derived the universal scaling form of the front of the fermionic correlations in terms of the Airy kernel. However, comparison with the numerical simulations at finite time shows large deviations between the scaling form and the numerical data indicating that the leading order corrections are important. The calculation of these corrections are explained in Sec. 5.2. For clarity, let us spell out in a bit more detail the derivatives of the function F (k, q) in Eq. (58). Using θ κ = θ −κ = 1, The derivative of g(k) can be computed from Eq. (173). In Fig. 8 we compare the numerical data with the analytic results with and without the corrections. In all the cases taking into account the corrections yields significantly better agreement.

F Correlation length via the Szegő lemma
We would like to simplify expression (91) for the correlation length. Using Θ(x) 2 = Θ(x), Θ(−x) = 1 − Θ(x) we find (194) Let us distinguish the cases u ≥ 0 and u < 0. When u ≥ 0, (195) Then we can split the domain of integation in (91) into three parts as (196) In each integral, the two step functions in the first term give equal contributions due to the symmetry of the integrand. In the second terms we use again identity Θ(−x) = 1 − Θ(x) for one of the step functions. These manipulations lead to identical integrands in the two lines, and we arrive at Eq. (92).