Universal front propagation in the quantum Ising chain with domain-wall initial states

We study the melting of domain walls in the ferromagnetic phase of the transverse Ising chain, created by flipping the order-parameter spins along one-half of the chain. If the initial state is excited by a local operator in terms of Jordan-Wigner fermions, the resulting longitudinal magnetization profiles have a universal character. Namely, after proper rescalings, the profiles in the noncritical Ising chain become identical to those obtained for a critical free-fermion chain starting from a step-like initial state. The relation holds exactly in the entire ferromagnetic phase of the Ising chain and can even be extended to the zero-field XY model by a duality argument. In contrast, for domain-wall excitations that are highly non-local in the fermionic variables, the universality of the magnetization profiles is lost. Nevertheless, for both cases we observe that the entanglement entropy asymptotically saturates at the ground-state value, suggesting a simple form of the steady state.


I. INTRODUCTION
The nonequilibrium dynamics of isolated many-body systems is at the forefront of developments within quantum statistical physics research, see Refs. [1][2][3] for recent reviews. Particularly interesting is the case of integrable models in one dimension, where the dynamics is constrained by a large set of conserved charges [4], leading to peculiar features in the transport properties [5] or the relaxation towards a stationary state [6]. A key paradigm of closedsystem dynamics is the quantum quench, where one typically prepares the pure ground state of a Hamiltonian which is then abruptly changed into a new one and the subsequent unitary time evolution is monitored [7]. The most common setting, well suited to study relaxation properties, is a global quench where the Hamiltonian is translational invariant both before and after the quench. However, if one is interested in transport properties far from equilibrium, some macroscopic inhomogeneity has to be present in the initial state.
In the context of spin chains, the simplest realization of such an inhomogeneity is a domain wall. In particular, for the XX chain a domain wall can be created by preparing the two halves of the system in their respective ground states with opposite magnetizations [8]. In the equivalent free-fermion representation, the simplest case of the maximally magnetized domain wall corresponds to a step-like initial condition for the occupation numbers. Under time evolution the initial inhomogeneity spreads ballistically, creating a front region which grows linearly in time. While the overall shape of the front is simple to obtain from a hydrodynamic (semi-classical) picture in terms of the fermionic excitations [9], the fine structure is more involved and shows universal features around the edge of the front [10,11] The melting of domain walls has been considered in various different lattice models, such as the transverse Ising [12,13], the XY [14] and XXZ chains [15][16][17][18], hard-core bosons [19][20][21], as well as in the continuum for a Luttinger model [22], the Lieb-Liniger gas [23] or within conformal field theory [24,25]. Instead of a sharp domain wall, the melting of inhomogeneous interfaces can also be studied by applying a magnetic field gradient, which is then suddenly quenched to zero [26][27][28]. Mappings from the time-evolved state of an initial domain wall to the ground state of a specific Hamiltonian have also been established [26,29]. Very recently, domain-wall melting in disordered XXZ chains has been studied as a probe of many-body localization [30].
Here we consider another realization of a domain wall which is created in the ordered ferromagnetic phase of the transverse Ising chain. Starting from one of the symmetry-broken ground states of the model, the order-parameter magnetization can be reversed along half of the chain. Due to the asymptotic degeneracy of the ordered states, this is still an eigenstate of the Ising chain locally, except for the neighbourhood of the kink in the magnetization where the domain-wall melting ensues.
The above setting has recently been studied numerically on infinite chains [31], using a matrix product state [32] related method, with two slightly different realization of the domain wall. For the excitation that is local in terms of Jordan-Wigner fermions, a very interesting observation on the magnetization profiles was made. Namely, it was pointed out that, after normalizing with the equilibrium value of the magnetization, the resulting snapshots of the profiles taken at times ht (i.e. rescaled by the value of the transverse field h) all collapse onto each other to almost machine precision [31]. Furthermore, the universal profile was conjectured to be identical to the one [8] obtained for the free-fermion chain with the step-like initial state.
In this paper we revisit this problem and show that these features can be understood analytically. First, we show that a very simple semi-classical interpretation of the front profiles in the hydrodynamic scaling regime can be found.
Moreover, in the limit of an infinite chain, even the fine structure of the profiles can be recovered by using a form-factor approach, providing an analytical support for the universality. For all of these results it turns out to be crucial that the domain-wall excitation is created by acting with a local fermion operator. Indeed, for a non-local realization of the same initial profile, the universality of the time-evolved front is lost and even the semi-classical picture breaks down.
The exact relation between the front profiles of the Ising and free-fermion domain-wall problems is quite remarkable. Indeed, in the latter case the time evolution is governed by a critical Hamiltonian whereas for the Ising chain we are always in the non-critical ferromagnetic regime. Despite the universal form of the magnetization profiles, one expects that this difference should clearly be reflected on the level of the time-evolved states. In fact, we will show that the entanglement entropy in the Ising chain always saturates for large times, in sharp contrast to the free-fermion case where it has a logarithmic growth in time [26,33,34]. Therefore, entanglement perfectly witnesses the non-criticality of the underlying Hamiltonian. Moreover, our results also indicate that the entropy in the non-equilibrium steady state of the Ising chain is equal to its ground-state value, suggesting that a non-trivial unitary transformation between these two states should exist.
The structure of the paper is as follows. In the next section we introduce the model and set up the basic formalism. The magnetization profiles for the Jordan-Wigner excitation are calculated in Sec. III using a number of different approaches. The results are then contrasted to those obtained for a non-local fermionic realization of the domain wall in Sec. IV. The time evolution of the entanglement entropy is discussed in Sec. V for both kinds of initial states. In Sec. VI we show that some of the above results can naturally be carried over to the XY chain by duality. We conclude with a discussion of our results and their possible extensions in Sec. VII. The manuscript is supplemented by three appendices with various details of the analytical calculations.

II. MODEL AND SETTING
We consider a finite transverse Ising (TI) chain of length N with open boundaries, defined by the Hamiltonian The diagonalization of H T I follows standard practice by mapping the spins to fermions via a Jordan-Wigner transformation [35]. For the open chain it will be most convenient to work with Majorana fermions defined by and satisfying anticommutation relations {a m , a n } = 2δ m,n . The Jordan-Wigner transformation brings the Hamiltonian into a quadratic form in terms of the Majorana operators which can be further diagonalized via The η k are standard fermionic operators satisfying {η k , η † l } = δ k,l and bring the Hamiltonian into the diagonal form The spectrum ǫ k in Eq. (4) and the vectors φ k and ψ k in Eq. (3) follow from the eigenvalue equations that are solved numerically with the matrices We will now consider the ordered phase (h < 1) of the TI model. It is well known that one has an exponentially vanishing gap in the system size N , and the ground and first excited states become degenerate in the thermodynamic limit N → ∞. For finite sizes, however, one has where | ⇑ and | ⇓ denote the macroscopically ordered states with a finite magnetization pointing in the ±x direction. Note, that for both |0 and |1 the magnetization vanishes since they respect the spin-flip symmetry of the Hamiltonian. The magnetization in the symmetry-broken ground states can thus be computed as Starting from the symmetry broken ground-state, we introduce two different types of initial states with a domainwall magnetization profile Here JW stands for Jordan-Wigner, since the excitation is created by applying a single Majorana fermion, see (2). In contrast, DW is a simple domain-wall excitation which is, however, non-local in terms of the Majorana operators. It is easy to check that both of the above excitations simply flip the magnetization in the x-direction for all spins up to site n 0 − 1 Additionally, the JW excitation creates a spin-flip in the z-direction at site n 0 . To simplify the setting, we will consider a symmetric domain wall, n 0 = N/2 + 1 with N even, in all of our numerical calculations. It should also be stressed that the domain wall is now in the longitudinal direction, in contrast to previous studies where a domain wall of the transverse magnetization was prepared. The equilibrium magnetization can be computed by evaluating the matrix element in (9). Rewriting σ x n with Majorana operators one has where we used |1 = η † 1 |0 , corresponding to the the lowest-lying excitation with ǫ 1 → 0 for N ≫ 1. Note that the vectors φ 1 (m) and ψ 1 (m) defining the mode η † 1 are localized around the left/right boundary of the chain, with elements decaying exponentially on a characteristic boundary length scale ξ b ∝ | ln h| −1 [36,37].
We thus have to evaluate the expectation value of a string of Majorana operators in the ground state, which can be factorized, according to Wick's theorem, into products of two-point functions. The latter can be calculated as 0|a j a l |0 = δ j,l + iΓ j,l , where the antisymmetric covariance matrix has a 2 × 2 block structure with matrix elements given by One further needs the matrix elements with the edge mode Finally, the magnetization at site n can be written as a Pfaffian of a 2n × 2n antisymmetric matrix [38,39] ⇑ |σ x n | ⇑ = Pf(M 0 ), Here Γ denotes the (2n − 1) × (2n − 1) reduced covariance matrix, whereas H (resp. its transpose) is a column (row) vector of length 2n − 1. The expression in (16) turns out to be real. Indeed, due to the simple checkerboard structure (14) of Γ, with nonvanishing elements only in the offdiagonals of the 2 × 2 blocks, the evaluation of the Pfaffian actually reduces to the calculation of the following n × n determinant

III. EVOLUTION OF MAGNETIZATION AFTER JORDAN-WIGNER EXCITATION
After having set up the basic formalism, we are now ready to consider the time evolution. First we deal with the JW excitation, where the time-evolved state reads The most important observable we are interested in is the order parameter magnetization σ x n , for which results can be obtained using a number of different approaches. First, we follow along the lines of the previous section and derive analogous Pfaffian formulas for the time-evolved magnetization which are exact for open chains of finite size. The scaling behaviour of the results suggests that a simple interpretation within a semi-classical approach should exist, which is presented in the second subsection. To study the fine structure of the profile directly in the thermodynamic limit, N → ∞, one has to follow a different route using the form-factor approach. At the end of the section we shortly discuss also the time evolution of the transverse magnetization σ z n .

A. Pfaffian approach
Instead of using the time-evolved state of Eq. (18), it is easier to work in a Heisenberg picture where the operators evolve as σ x n (t) = e iHT I t σ x n e −iHT I t . The time-dependent magnetization can then be obtained by taking expectation values in the initial state and can be written as The above formula is analogous to that of Eq. (12), one has, however, Heisenberg operators in the product surrounded by two extra a 2n0−1 and thus one has to evaluate a string of 2n + 2 operators. In order to apply Wick's theorem, one first needs the time evolution of the Majorana operators where the matrix elements of the propagator R are given by It is easy to show that the two-point functions of the Heisenberg operators do not change in time, 0|a j (t)a l (t)|0 = 0|a j a l |0 . Indeed, since the Hamiltonian is unchanged in our protocol (i.e. there is no quench involved), the exponential factors in the Heisenberg operators act trivially on the ground state. However, the expectation values of products of operators at different times becomes nontrivial and, using (20) and (13), can be evaluated as The remaining two-point functions are given by 0|a 2n0−1 a 2n0−1 |0 = 1 and 0|a 2n0−1 η † 1 |0 = φ 1 (n 0 ). With all the ingredients at hand, one can again arrange the two-point functions in an antisymmetric matrix M of size (2n + 2) × (2n + 2) and calculate its Pfaffian. However, the calculation can be simplified using the special properties of Pfaffians, as shown in detail in Appendix A. In turn, the result can be written in a form analogous to the equilibrium case whereM is a matrix of size 2n × 2n and its entries are given bỹ HereΓ andH are again a matrix of size (2n − 1) × (2n − 1) and a vector of length 2n − 1, respectively. However, due to the extra terms in Eq. (24), the matrixM does not have the simple checkerboard structure as in equilibrium, and thus the result cannot be rewritten as a n × n determinant. Nevertheless, it is easy to show (see Appendix A) that all the elements ofΓ are real, and the only imaginary entries inH are due to H 2m , see Eq. (15). Therefore, taking the real part in (23) is equivalent to setting H 2m = 0 and calculating a real-valued Pfaffian, which can be performed by efficient numerical algorithms [40].
The results for the magnetization are shown in Fig. 1 for a chain of length N = 200. One can observe a number of features from the unscaled profiles at fixed t = 50 (shown on the left). In particular, it is easy to see that the edges of the expanding front are located at a distance ≈ ht measured from the initial location n 0 − 1/2 = (N + 1)/2 of the domain wall. From this it is easy to infer that the maximum speed of propagation is given by v = h which will be verified by the semi-classical approach of the next subsection. Beyond the edge of the front one recovers, up to exponential tails, the equilibrium profile, which shows well-known boundary effects [36,37] on a length scale ξ b close to the ends of the chain.
To better understand the behaviour of the front, one should compare snapshots of the magnetization, normalized by the equilibrium value, for various fields h but keeping the scaling variable ht = 50 fixed, as depicted on the right of Fig. 1. Remarkably, as already noted in Ref. [31] for infinite chains, the data sets collapse close to machine precision on a single curve, which turns out to be identical to the one [8] for the free-fermion chain with a step-like initial state.

B. Semi-classical approach
To interpret the above results, we now present a very simple semi-classical argument which yields the correct magnetization profile for the JW excitation in the scaling regime, i.e. |n − n 0 | → ∞ and ht → ∞ with (n − n 0 )/ht kept fixed. To simplify the discussion, here we work directly with an infinite chain, with no boundary conditions imposed. Due to perfect translational invariance, the eigemodes created by η † q are now propagating waves with continuous momenta chosen from the interval q ∈ [−π, π].
In the context of the TI chain, the semi-classical reasoning was originally presented by Sachdev and Young [41], and has since been used to obtain the magnetization for various (global or local) quench protocols [42,43]. The argument is as follows: the initial Majorana operator a 2n0−1 which excites the domain wall is, in fact, a mixture of the various eigenmodes excited by η † q . One could think of each of these modes as an elementary domain-wall excitation which propagates at a given speed with ǫ q the dispersion of the TI chain. To get the magnetization at a given site n, one simply has to determine the density N of excitations that have sufficient velocities to arrive from the initial location n 0 to the point of observation in time t. Introducing the scaling variable v = (n − n 0 )/t and focusing on v > 0, one has In other words, N (v) is just the fraction of domain-wall excitations that have a velocity larger than v. This can be obtained via (25), by solving the equation v q± = v for q ± (v), which is represented graphically on Fig. 2. The analytical solution can be found by introducing the variable z = sin q, which leads to the following quadratic equation with the roots given by Finally, the solution for the wavenumbers reads where v 0 is the solution of the equation z + (v 0 ) = 1.
Using the identity arcsin(z) + arccos(z) = π/2 and the addition formulas for the arccos function, the difference of the wavenumbers can be written as From the solutions (28) one finds It is easy to see that the right hand side term within the absolute value changes sign exactly at v = v 0 , hence from (30) one finds that the minus sign applies for all values of v. Finally, substituting both terms one arrives at which is exactly the free-fermion result [8] with the velocity rescaled by h.

C. Form-factor approach
The semi-classical approach yields a very simple physical explanation for the magnetization profile in the scaling limit |n − n 0 | → ∞ and t → ∞ with the ratio v = |n − n 0 |/t kept fixed. However, it does not account for the perfect collapse of the normalized magnetization curves JW|σ x n (t)|JW / ⇑ |σ x n | ⇑ even at finite times for a fixed value of ht, see Fig. 1. To capture the fine-structure of the profile, one can follow a form-factor approach which was used successfully to obtain results for the magnetization in case of a global quench [44]. In contrast to the Pfaffian approach, which is well-suited for open chains of finite size, the form-factor approach works most efficiently in the thermodynamic limit.
Our starting assumption for the semi-classical approach was that the JW excitation is a mixture of the various single-particle eigenmodes. It turns out that, to make this statement rigorous, one has to consider a TI chain with antiperiodic boundary conditions σ x N +1 = −σ x 1 . The Hamiltonian H for the antiperiodic chain can be diagonalized by the very same procedure as the periodic one, which is summarized in Appendix B. The main feature of both geometries is that the Hilbert space splits up into the Neveu-Schwarz (NS) and Ramond (R) sectors, which differ by their symmetry properties with respect to a global spin-flip transformation. In particular, the vacua of the two sectors, |0 NS and |0 R , are analogous to those |0 and |1 of the open chain in Eq. (8), and the symmetry-broken ground states are obtained as their linear combinations. In turn, the time-evolved magnetization after the JW excitation is given by The Jordan-Wigner excitation a 2n0−1 can now be rewritten in the basis that diagonalizes the Hamiltonian, as shown in (B13) of Appendix B. It will populate the vacua as where the θ q and θ p are Bogoliubov angles defined in (B10). The momenta q of the NS sector (respectively p of the R sector) are quantized differently: they are half-integer (integer) multiples of 2π/N . Most importantly, the singleparticle states |q NS and R p| are exact eigenvectors of the antiperiodic Hamiltonian. Hence, their time evolution becomes trivial with the dispersion relation defined in (25). The role of the antiperiodic boundary conditions should be stressed at this point, since the eigenvectors of the periodic TI chain always have an even number of single-particle excitations. Clearly, thanks to the simple time evolution in (35), the only remaining ingredients we need are the form factors R p|σ x n |q NS between the single-particle states. Fortunately, for the particular fermionic basis at hand, the form factors are known exactly and in the limit N → ∞ are given by [45] The vacuum matrix element in the denominator of the left hand side is simply the equilibrium magnetization. In fact, the form factors for finite N are also known exactly [45], but we are only interested in the thermodynamic limit.
Combining the results (34)- (36) and turning the sums into integrals, one finally arrives at the result for the normalized magnetization To show the identity with the free-fermion result, one has to evaluate the above double integral. First, one notices that the Dirichlet kernel appears in the integrand of (37) which can be rewritten as leaving us with a sum of integrals with simpler integrands to evaluate. Assuming that the result depends on the scaling variable ht only (see Fig. 1), one can show after a rather tedious exercise, the details of which are given in Appendix C, that each of these integrals reproduce the square of a Bessel function π −π dp 2π π −π dq 2π Consequently, the normalized magnetization profile is obtained in the simple form which is indeed the free-fermion result of Ref. [8].
Finally, it should be pointed out that the semi-classical result (26) in the scaling limit, with the scaling function (32), could also be obtained from a saddle-point approximation of the double integral (37) along the lines of Ref. [46].

D. Transverse magnetization
To conclude this section, we shortly discuss the time-evolution of the transverse magnetization σ z n for the open-chain geometry. As remarked earlier, the JW excitation flips the transverse spin at site n 0 and the disturbance spreads out in time. Since σ z n is an even operator in terms of the fermions, it has only diagonal matrix elements w.r.t. the ground states |0 or |1 , and these will coincide for N ≫ 1. In turn one has where the matrix elementΓ 2n−1,2n can be calculated according to Eq. (A4) of Appendix A. In the limit N → ∞, Eq. (41) can be shown to coincide with the corresponding result of Ref. [31]. The normalized profiles of the transverse magnetization are shown in Fig. 3. The spreading of the initially flipped spin can be seen on the left for h = 0.5. Interestingly, the total transverse magnetization seems to be conserved to a very good precision, at least until the front reaches the boundary region. On the right we show the profiles for various h but with a fixed value of ht = 25. Clearly, in contrast to the order-parameter, the normalized transverse magnetization is not a function of ht only.

IV. MAGNETIZATION PROFILES FOR DOMAIN-WALL EXCITATION
In the previous section we have shown that the profiles for the JW excitation can be obtained using various approaches and the underlying physics can be understood by a simple semi-classical argument. We now turn our attention towards the simple domain-wall excitation [31], defined on the right of Eq. (10). Although the difference from the JW excitation seems innocuous in the spin-representation, due to the non-locality of the Jordan-Wigner transformation, the DW excitation becomes a string of Majorana operators. Analogously to Eq. (19), the magnetization can now be written as The expectation value in (42) can still be written as a Pfaffian, albeit with a matrix of much larger size. To this end we define the rectangular matrices C and D of unequal-time two-point functions with elements where the capitalized index J runs over the set J = 1, . . . , 2n 0 −2, whereas j = 1, . . . , 2n−1 as before. Similarly, one can introduce the reduced covariance matrix Γ 0 (with elements Γ J,J ′ ) and the column vector H 0 (with elements H J ) where again J, J ′ = 1, . . . , 2n 0 −2. Using these definitions, the magnetization can be written as a (4n 0 −4+2n)×(4n 0 −4+2n) Pfaffian, given explicitly in (A6). Furthermore, performing manipulations similar to the JW case (see Appendix A), the expression can again be reduced to a Pfaffian of size 2n × 2n given by The Pfaffian in (44) can be evaluated numerically with the results shown in Fig. 4. The normalized magnetization profiles are plotted against the distance from the initial location of the domain-wall, rescaled by ht. On the left of Fig. 4, we show the profiles at fixed ht = 50 and for various values of h. From the figure it becomes evident that the universality is lost as one finds no data collapse. Moreover, even the semi-classical result valid for the JW case and shown by the solid line, breaks down for the DW excitation: while the agreement for h = 0.5 still seems to be fairly good, the deviations increase dramatically when approaching the critical value of the field h → 1.
The breakdown of the semi-classical picture does not come entirely unexpected. In fact, the DW initial state consists of a string of Majorana excitations extending over the left half-chain, which cannot any more be considered as a mixture of single-particle excitations in the momentum space. This becomes even more obvious in terms of the form-factor approach of the previous section. Indeed, in the DW case one has to consider many-particle form factors instead of the single-particle matrix elements of Eq. (36). Since these form factors become highly involved with increasing particle number [45], such a calculation is beyond our reach. Nevertheless, for a fixed value of h, one still has a ballistic expansion with the maximal signal velocity given by h, as demonstrated by the rescaled data on the right of Fig. 4.
Finally, one could also have a look at the transverse magnetization. Although, in contrast to the JW case, the initial state does not have any flipped spin in the z-direction, the profile will not remain constant. Indeed, one observes a signal front (albeit much weaker than in the JW case) propagating outwards from the location of the domain wall with the same speed v = h. In complete analogy with Eq. (41), the transverse magnetization for the DW case is given by DW|σ z n (t)|DW =Γ 2n−1,2n , with the corresponding matrix element defined in (A9).

V. ENTANGLEMENT EVOLUTION
Given the universal result (40) for the magnetization profile in the JW case, the question naturally emerges whether one has a deeper connection to the free-fermion domain-wall problem on the level of the time-evolved state. To answer this question, we shall now consider the entanglement entropy, which carries important information about the state itself. Entanglement evolution has been considered previously in Refs. [15,26,28,31,33,34,47] for various domainwall initial states.
The time evolution of the entanglement entropy is given by Even though there are several analytical approaches to obtain the entanglement entropy for Gaussian states of the TI chain (see e.g. Ref. [48] and references therein), the situation here is more subtle. Indeed, the initial state |JW is defined in terms of the symmetry-broken ground state which, however, is not itself Gaussian but rather the superposition of two Gaussian states. Thus we will determine the entropy via density matrix renormalization group (DMRG) related calculations [49].
Within the matrix product state (MPS) formalism of DMRG [32], the ground state can be approximated to a very high precision by the ansatz where |s denotes the spin basis states and A si are auxiliary matrices with a variable bond dimension. They can be obtained by minimizing ⇑ |H| ⇑ with respect to the product A si A si+1 for a given site index i, while keeping all the other matrices fixed. Repeating the procedure for every pair of neighbouring lattice sites, the MPS will converge to the ground state after several sweeps. To ensure that we end up in the symmetry-broken ground state | ⇑ , we introduced a small longitudinal field h x > 0 in the Hamiltonian H = H T I − h x i σ x i for the first few sweeps and set h x = 0 afterwards, until convergence is reached. The JW and DW excitations can then be created by acting with their (trivial) matrix product operator representations on the ground-state MPS. Finally, the time-evolution of the states was implemented with the finite two-site time-dependent variational principle (TDVP) algorithm [50].
We start by looking at the initial entropy profile at time t = 0. One should point out that the state |JW is created by acting on the symmetry-broken ground state by a product of strictly local (on-site) terms, which do not modify the entropy. One thus expects that the result for the real ground state |0 should be recovered, except for a ln 2 contribution coming from the degeneracy, which is now removed. In the limit of an infinite chain N → ∞, this is given via elliptic integrals by the analytical expression [51,52] with h ′ = √ 1 − h 2 , which we recover in the bulk (ξ b ≪ n ≪ N − ξ b ) of the chain. For cuts within about a distance ξ b from the boundaries, the profile becomes inhomogeneous with increasing entropies for smaller subsystems.
The time evolution is depicted on Fig. 5. On the left we plot the entropy of a half-chain (n = n 0 ), with S(0) subtracted, against the scaling variable ht. After a sudden increase, the curves show a slower, oscillatory approach towards an asymptotic value which seems to be given by ln 2. Interestingly, the approach takes place from below, with the curves never crossing the asymptotic value. The distance of the maxima is given by ht = π to a good precision. Note, however, that the collapse of the curves against the variable ht is good but not exact. On the right of Fig. 5 we show the full profiles for h = 0.5 and various times, again with S(0) subtracted and with the distance n − n 0 of the cut from the centre rescaled by ht. The rescaled profiles converge towards a scaling function for large times, which remains unchanged for other values of h (not shown on the figure) as well. To interpret the above results, it is useful first to compare them to the corresponding result for the free-fermion case. There the entropy profile of an infinite chain is found to be given by the function [33,34] with the rescaled distance v = (n − n 0 )/t and |v| < 1. The latter profile is not only a function of v, but one has a contribution which grows logarithmically in time. Obviously, this is not the case for the JW excitation of the TI chain. In fact, the difference in the results gives perfect account about the underlying Hamiltonians: while for the free fermion the time-evolution is governed by a critical Hamiltonian, for the TI one is always in the h < 1 non-critical regime and thus the entropy saturates. It is important to stress that this difference is not revealed by looking only at the magnetization profiles, which are identical after rescaling. Despite the difference in the entropy profiles, there is one important analogy which can be uncovered. We have observed (see left of Fig. 5) that the t → ∞ result for the half-chain entropy is given by S = S(0) + ln 2. However, as pointed out before, this is nothing else but the entropy of the real ground state |0 . Moreover, from the scaling behaviour (see right of Fig. 5) one infers, that the same is true for arbitrary cuts with |n − n 0 |/ht → 0, i.e. finite distances from the centre and infinite time. This is exactly the regime, where a translational invariant current-carrying steady state is formed. Hence, no matter where we cut the system within the steady-state regime, we always get an entropy that is equal to the ground state value. Since the entropy gets contributions only from a distance of order ξ of the correlation length measured from the cut, this suggests that the steady state, i.e. the reduced state of a finite segment of size L ≫ ξ in the limit t → ∞, is unitarily equivalent to the reduced density matrix of the ground state In fact, this is exactly the case for the free-fermion chain, where the steady state is simply given by a boosted Fermi sea [28,53]. Thus, taking a finite subsystem of length L on the right-hand side of site n 0 in the free-fermion chain, the asymptotic entropy for t ≫ L is given by the ground-state value S = 1/3 ln L + const., with the non-universal constant ≈ 0.726. In this sense, the two results are completely analogous.
Finally, we study the entropy evolution also for the DW case. The results for the half-chain entropy as well as for the profiles are shown in Fig. 6. Although for smaller values of h the half-chain entropy looks qualitatively similar to the JW case, for h = 0.9 there are noticeable differences. Namely, the increase for early times becomes slower, whereas for large times one has additional oscillations. Nevertheless, the asymptotical value of S(t) − S(0) still seems to be given by ln 2. Although the rescaled profiles collapse onto each other for fixed h and various times, the shape of the scaling curves slightly changes for different values of h in the DW case, as shown on the right of Fig. 6.

VI. DUALITY WITH THE ZERO-FIELD XY CHAIN
It is natural to ask whether the results found for the magnetization and the entropy of the TI chain could exist for a broader universality class of spin models with ferromagnetic ground states. In the following we will show that the result naturally carries over to JW-type excitations of the anisotropic XY chain in zero magnetic field. Let us consider a chain of 2N sites defined by the Hamiltonian Applying the duality transformations [54][55][56][57] the XY Hamiltonian decomposes into the sum of two TI chains, defined in terms of the dual variables as Here the magnetic fields are defined as Thus, the ground state of the XY Hamiltonian corresponds to the direct product of TI ground states on the corresponding sublattices. Note that, for 0 < γ < 1, the Hamiltonian H T I,1 is in its ordered phase, whereas H T I,2 is in the disordered phase.
To calculate the magnetization for the XY chain, one has to rewrite the σ x operators in the dual variables Since both the ground states as well as the operators factorize on the two sublattices, one can write where we have used the fact that the lowest lying excitation of XY corresponds to exciting the ordered T I chain H T I,1 only. The result for σ x 2i is similar. Furthermore, one can also construct the Majorana operators using the representation of the string variables in the dual language In fact, this operator creates nothing else but a JW excitation of H T I,1 (up to an irrelevant sign factor) while the ground state of H T I,2 is left untouched Finally, since the two TI Hamiltonians commute [H T I,1 , H T I,2 ] = 0, the time evolution operator also factorizes and one arrives at the relation Hence, after proper rescaling, one indeed finds the universal free-fermion result (40) both on even and odd lattice sites. The relation (60) has also been checked against DMRG calculations with an excellent agreement. The same argument also applies to the entanglement entropies and yields S XY (t) = S T I,1 1 + γ 2 t + S T I,2 (0) .
Note that similar duality relations between entropies of XY and T I chains were found earlier for the ground state [58] as well as for local quenches [59].

VII. CONCLUSIONS
We have studied the domain-wall melting for particular initial states of the ferromagnetic TI chain. For the JW excitation that is local in terms of the fermion operators that diagonalize the Hamiltonian, the longitudinal magnetization profiles after proper rescaling are completely identical to the ones observed for a fermionic hopping chain with step initial condition. The result carries over to the anisotropic XY chain with h = 0. The entanglement entropy is, however, found to saturate during time evolution and signals the non-criticality of the underlying Hamiltonian.
The case of the non-local DW excitation is quite different. In particular, the semi-classical approach, that yields the correct JW profiles in the scaling limit, breaks down and thus we have not been able to find an analytical result for the DW profiles. It might be possible to derive some results via the form factor approach which, however, also becomes highly involved and we have thus left this question open for future studies.
There are also a number of natural extensions of this work. First of all, one should check if the universality of the JW magnetization profiles extends to the full ferromagnetic phase of the XY model. A further step would be to investigate more general spin chains, such as the XXZ chain, that cannot be transformed into free fermions. While we do not expect the full universality for the fine structure of the profile to hold in this case, some essential features might still be inherited. It would also be instructive to compare the results to a quench setting, where the | ⇑ and | ⇓ states are prepared as the symmetry-broken ground states of two half chains which are then joined together.
Finally, our results for the entropy lead us to the conjecture that the non-equilibrium steady state is locally (i.e. in the region where the front has already swept through) related to the symmetry unbroken ground state of the TI chain. It would be interesting to find further evidence by comparing more complicated observables, such as spin correlation functions, which could also be obtained from the Pfaffian formalism.
• Multiplication of a row and a column by a constant is equivalent to multiplication of the Pfaffian by the same constant.
• Simultaneous interchange of two different rows and corresponding columns changes the sign of the Pfaffian.
• A multiple of a row and corresponding column added to another row and corresponding column does not change the value of the Pfaffian.
• For a 2n × 2n antisymmetric matrix M and constant λ one has Pf(λM ) = λ n Pf(M ) • The Pfaffian of a 2n × 2n antisymmetric matrix M can be expanded into minors according to the reduction rule where M (1,j) is a (2n − 2) × (2n − 2) antisymmetric matrix obtained by removing the first and j-th rows and columns of M .
The above rules are very similar to the properties of determinants, except that one has to manipulate the rows and columns simultaneously.

JW excitation
We first deal with the simpler JW excitation. According to (19), the magnetization is given by the expectation value of a string of 2n + 2 Majorana operators. Hence, it can be rewritten as the following Pfaffian where we used a block-notation with (2n − 1) × (2n − 1) matrix Γ and column-vectors H, C and D of length 2n − 1 defined in (15) and (22). Note that the transpose of the above vectors simply give the corresponding row-vectors. The remaining entries correspond to the expectation values 0|a 2n0−1 a 2n0−1 |0 = 1 and 0|a 2n0−1 η † 1 |0 = φ 1 (n 0 ).
We can now use the Pfaffian rules above to transform the matrix M into matrices of simpler structure M ′ and M ′′ given by In the first step, we have subtracted the third row/column of M from the first ones which yields M ′ on the left of Eq. (A3). Subsequently, one can subtract the third row (column) of M ′ multiplied by C + D (respectively C T + D T ) from the second row (column) which then leads to M ′′ on the right of Eq. (A3), withΓ andH given in (24) of the main text. Clearly there is now only one nonzero entry in the first row/column of M ′′ , it can thus be reduced to a smaller matrix of size 2n × 2n by removing the first and third rows and columns. Indeed, using Eq. (A1), the only nonvanishing contribution is with j = 2n + 1 which gives and extra sign for the reduced Pfaffian. Finally, the factor (−i) n−1 in (A2) can be absorbed by multiplying all the matrix elements by −i, except for the last row and column. This yields the final result in Eq. (23). It is instructive to write out explicitly the matrix elements ofΓ andH Note that all the entriesΓ i,j are real, and the only imaginary entries inH j appear for j = 2n due to H 2n = iψ(n). In particular, for t = 0 the propagator R i,j = δ i,j is given by the identity and one has Now, if n < n 0 , the extra terms in (A5) do not give any contribution such that theM = M 0 and so the magnetization −Re Pf(M ) is given by −1 times the equilibrium one. On the other hand, for n ≥ n 0 , the extra contributions simply reverse the sign of the 2n 0 − 1-th row and column of M 0 , giving an extra sign and reproducing the equilibrium magnetization.
We will again manipulate the matrix M and transform it to simpler forms M ′ and M ′′ given by In the first step, we do a row-by-row (resp. column-by-column) subtraction of the matrices in the third row (column) from the first ones in the block matrix M . This zeroes out the entries iΓ 0 and H 0 in the first row/column and transforms −C → −(C + D) (resp. C T → C T + D T ). The remaining ±1 1 can be used to cancel out the iΓ 0 matrix in the third diagonal entry of M , by subtracting iΓ 0 /2 (resp. its transpose) times the first row/column from the third ones. This yields M ′ of Eq. (A7) with a modified rectangular matrix defined as In the next step, we can cancel out the remaining entries C T + D T and its transpose from the first row and column by subtracting the respective multiple of the third column/row from the second ones, which leads to M ′′ in Eq. (A7) withΓ andĤ defined in (45). Now, we can continue with the reduction of the matrix. The ±1 1 in the first row/column shows that one can eliminate 2 × (2n 0 − 2) rows/columns consecutively, reducing again the matrix to a size of 2n × 2n. According to (A1), every second step in the reduction gives a sign, which amounts to a factor (−1) n0−1 and cancels out with the respective sign term in (A6). Finally, the (−i) n−1 can again be absorbed just like in case of the JW calculation, and leads to the result in Eq. (44) in the main text.
One can again have a look at the matrix elementsΓ i,j andĤ j . Evaluating the matrix products in (45), one is left with the following simple expression where the sum over J runs on the index set J = 1, . . . , 2n 0 − 2 whereas the sum overJ runs on the complement set J = 2n 0 − 1, . . . , 2N . It is easy to check how this again gives the correct result for t = 0, where R i,j = δ i,j . Indeed, setting n < n 0 , then since i, j ≤ 2n − 1 one has R i,J = 0 and R j,J = 0 for all i, j and thusΓ = Γ. However,Ĥ = −H and thus the last row/column of the Pfaffian is multiplied by −1 which changes its sign and thus the magnetization is given by −Pf(M 0 ). On the other hand, for n ≥ n 0 some of the matrix elements ofΓ will be changed. Indeed, one haŝ The above transformation simply amounts to multiplying all the columns/rows between 2n 0 −1 and 2n of the Pfaffian, each of which giving a sign. However, since there are an even number of rows and columns involved, in the end the value of the Pfaffian is unchanged and we get back the correct result Pf(M 0 ) for the magnetization.

Appendix B: Diagonalization of HT I with antiperiodic boundary conditions
The TI chain with antiperiodic boundary conditions is given by the same Hamiltonian as in Eq. (1), except that both sums run until m = N and we set σ x N +1 = −σ x 1 . To diagonalize it, we follow a slightly different route along the lines of Ref. [45]. Instead of working with Majorana fermions, we define creation/annihilation operators where σ ± n = (σ x n ± iσ y n )/2 and the commutation relations are given by c m , c † n = δ m,n . We also introduce the global spin-flip operator which commutes with the Hamiltonian [H, W ] = 0. In terms of the fermion operators it reads and the boundary condition for the fermions becomes c N +1 = W c 1 . Since W 2 = 1, the eigenstates of the Hamiltonian split up into two sectors: the Ramond (R) sector corresponding to eigenvalue W = 1 has periodic, whereas the Neveu-Schwarz (NS) sector with W = −1 has antiperiodic boundary conditions for the fermions. For our purposes it will be more convenient to work in a dual basis defined by The dual transformation interchanges the two terms in the Hamiltonian where the dual fermions satisfy {d m , d † n } = δ m,n and the same boundary condition d N +1 = W d 1 . One then introduces the Fourier modes where the momenta are quantized depending on which sector of the Hilbert space one chooses In terms of the Fourier modes, (B5) can be rewritten as where the summation goes over the momenta defined by (B7), but we omitted the k indices for notational simplicity. The above Hamiltonian can be diagonalized by a Bogoliubov transformation where the dual Bogoliubov angle is given by The diagonal form of the Hamiltonian and the one-particle spectrum read The many-particle eigenstates of the antiperiodic Hamiltonian can then be constructed as In fact, all the eigenstates have an odd number of excitations, as opposed to the periodic chain where the number of excitations is always even. Finally, it is useful to rewrite the Majorana fermions of section II in terms of the eigenmodes of the Hamiltonian which then leads directly to Eq. (34) in the main text.

Appendix C: Integral formulas
In this appendix we will evaluate the integral The factors involving the Bogoliubov angle can be written for q > 0 as cos First we will consider the simplest case k = 0. The integral then simplifies to I 0 = π 0 dp π π 0 dq π ǫ p ǫ q cos θ q 2 cos θ p 2 cos(ǫ q − ǫ p )t , where we made use of the symmetry under exchange of p and q and the fact that the similar integral with sin θq 2 sin θp 2 vanishes due its oddness under reflections θ −q = −θ q or θ −p = −θ p . To evaluate (C3) it is more convenient to introduce ǫ q = 1 + hǫ q (similarly for ǫ p ) and rewrite the integral in terms of theǫ variables. The change of the integration measure can be derived from In terms of the new variables the integral reads to arrive at the result I 0 = J 2 0 (ht). Unfortunately, the treatment of the general case k > 0 is much more cumbersome. On one hand, there are no simplifications due to symmetries of the integrand and thus one has many more terms appearing. On the other hand, even though the transformation to theǫ variables yields the natural scaling variable ht in the argument of the time-dependent cosine in (C1), it also transforms the term cos k(q − p) to a more complicated expression. Indeed, using trigonometric identities, the extra factors can be rewritten in terms of the Chebyshev polynomials cos kq = T k (cos q), sin kq = sin qU k−1 (cos q) (C7) where, however, the argument has to be reexpressed as and similarly for p. Applying trigonometric addition formulas in the other cosine terms as well, the integral splits into a number of terms I k = I 1,kÎ1,k + I 2,kÎ2,k + I 3,kÎ3,k + I 4,kÎ4,k , where we defined The integrals with the hat symbols are very similar to the ones defined above, but with an additional factor (1 + hǫ) in the integrand, analogously to (C5). Note that we used the shorthand notation z, defined in Eq. (C8), in the arguments of the Chebyshev polynomials to simplify formulas. The exact evaluation of the above integrals is a very cumbersome task, due to the fact that the variable z appears in the argument of the Chebyshev polynomials. Hence, the individual integrals I α,k (h, τ ) andÎ α,k (h, τ ) for α = 1, . . . , 4 depend on both variables h and τ = ht. Nevertheless, as it is clear from Fig. 1, the final result I k in (C9) depends only on the scaling variable τ = ht. To show this analytically, one has to use the explicit form of the Chebyshev polynomials and expand the powers of z, which then lead to integrals that can be evaluated via [60] for arbitrary integers m and n. In turn, each of the integrals I α,k (h, τ ) andÎ α,k (h, τ ) can be rewritten as a double sum of terms containing various powers of h and expressions of the form (C12). Due to the huge amount of terms appearing, we were able to verify the relation ∂ ∂h I k (h, τ ) = 0 only using Mathematica, for k < 20. Using this property, one can also obtain the final result by setting h = 0 with τ = ht fixed in all of the integrals. Then the argument of the Chebyshev polynomials simplifies to −ǫ, the integrals with the hat symbols are identical to the ones without, and both I 3,k and I 4,k in (C11) vanish explicitly. The remaining terms can be evaluated via the integral identities [60] 1 −1 dǫ π T 2l (ǫ) cos(ǫht) √ 1 −ǫ 2 = (−1) l J 2l (ht), leading to the final result I k = J 2 k (ht).