Front dynamics in the XY chain after local excitations

We study the time evolution of magnetization and entanglement for initial states with local excitations, created upon the ferromagnetic ground state of the XY chain. For excitations corresponding to a single or two well separated domain walls, the magnetization profile has a simple hydrodynamic limit, which has a standard interpretation in terms of quasiparticles. In contrast, for a spin-flip we obtain an interference term, which has to do with the nonlocality of the excitation in the fermionic basis. Surprisingly, for the single domain wall the hydrodynamic limit of the entropy and magnetization profiles are found to be directly related. Furthermore, the entropy profile is additive for the double domain wall, whereas in case of the spin-flip excitation one has a nontrivial behaviour.


Introduction
The nonequilibrium dynamics of integrable quantum many-body systems has been the focus of intensive research [1]. The interest in these peculiar models, characterized by the existence of a large set of conservation laws, comes from two main perspectives. On one hand, they show relaxation towards generalized stationary ensembles that are not described by conventional statistical mechanics [2]. On the other hand, owing to the presence of stable quasiparticle excitations, integrable models have anomalous transport properties [3]. A recent milestone in understanding the transport driven by an initial inhomogeneity has been the formulation of generalized hydrodynamics (GHD) [4,5], which gives accurate predictions for the profiles of conserved densities in an appropriate spacetime scaling limit.
The simplest paradigm of an inhomogeneous initial state is a domain wall, separating domains of spins with different magnetizations. Letting the system evolve, the domain wall starts to melt, giving rise to an expanding front region characterized by a nonzero spin current. The resulting magnetization profiles were studied in various integrable spin models such as the XX chain [6][7][8], the transverse Ising (TI) [9][10][11], the XY [12] as well as the XXZ chains [4,[13][14][15][16][17]. Rather generically one finds ballistic transport, with the exception of the isotropic Heisenberg chain where a diffusive behaviour is observed instead [18][19][20][21][22][23]. The common feature in all of the examples above is that the domain wall is oriented along the z-axis, and thus the magnetization is a local operator in the fermionic representation of the corresponding spin chain. In particular, for models with fermion-number conservation, the transverse magnetization itself corresponds to a locally conserved density, which makes the problem directly amenable to GHD techniques.
Recently, however, domain walls created upon the symmetry-broken ferromagnetic ground states of TI or XY chains have been considered [24][25][26]. The ordering in these chains occurs in the longitudinal component of the magnetization, which is a highly nonlocal string operator in the fermionic picture, being nontrivially related to the local conserved densities. Hence, even though one has a free-fermion model at hand, it is a priori unclear whether a hydrodynamic description still holds for this observable. Nevertheless, in [25,26] it has been shown that, for domain walls excited by a single local fermion operator, the longitudinal magnetization profile has the usual hydrodynamic scaling limit one would naively expect. Namely, the profile is determined by noninteracting quasiparticles carrying the fraction of a spin-flip and traveling at the corresponding group velocity.
In the present work we extend these studies to excitations that can be written as the product of two local fermion operators. In the spin language they describe a double domain wall, and if the distance between them is sufficiently large, we find that the magnetization profile factorizes in the hydrodynamic scaling limit. In other words, the quasiparticle excitations created at the two domain walls are completely independent. In contrast, the situation becomes nontrivial if the fermionic excitations act on neighbouring sites, even though the product of two adjacent domain walls is just a spin-flip and thus perfectly local in the spin-representation. Indeed, it turns out that this composite fermionic excitation leads to interference effects between the quasiparticle modes, encoded in the form factors of the spin operator. This interference term yields a significant contribution to the hydrodynamic profile, which can be found analytically via stationary phase analysis.
We also study in detail the correlation functions and the entanglement entropy for the single domain wall excitation. Interestingly, both of them can be directly related to the magnetization. For the correlations we derive a relation which holds also for finite times if the separation of the spins is much larger than the correlation length. On the other hand, for the entropy we propose an ansatz that is motivated by recent results for single-mode quasiparticle excitations in a free massive quantum field theory (QFT) [27,28]. Our ansatz works perfectly in the hydrodynamic regime, thereby creating an exact relation between the magnetization and entanglement profiles. Furthermore, we observe that the entropy becomes additive for the double domain wall excitation, whereas for the spin-flip one has again a nontrivial behaviour due to the above mentioned interference terms.
The paper is structured as follows. We start by introducing the model in Sec. 2. The magnetization dynamics is studied in Sec. 3 for three different local excitations as well as for a local quench. The correlation functions are investigated in Sec. 4, followed by the study of the entropy profiles in Sec. 5. We discuss our findings in Sec. 6, and the technical details of the calculations are reported in three Appendices.

Model
We consider an XY spin chain of length N described by the Hamiltonian where σ α n are Pauli matrices located at site n, h and γ denote the transverse magnetic field and the XY anisotropy, respectively. We restrict ourselves to the parameter regime 0 < h < 1 and 0 < γ ≤ 1 where the chain is in a gapped ferromagnetic phase, with γ = 1 corresponding to the TI chain.
The Hamiltonian (1) is diagonalized through a standard procedure [29], by first introducing Majorana fermions via a Jordan-Wigner transformation satisfying anticommutation relations {a k , a l } = 2δ k,l . While (1) describes an open chain which is most suitable for our numerical calculations, the analytical treatment of the problem requires to consider either periodic (s = +) or antiperiodic (s = −) boundary conditions, σ x N +1 = sσ x 1 and σ y N +1 = sσ y 1 . Due to the global spin-flip symmetry of the model, the corresponding Hamiltonians can then be split into two parts In terms of the Majorana fermions, the corresponding symmetry sectors are described by the Hamiltonians which differ in the boundary conditions a 2N +1 = ±a 1 and a 2N +2 = ±a 2 being periodic for the Ramond (R) and antiperiodic for the Neveu-Schwarz (NS) sectors. In order to diagonalize (4), one performs a Fourier transformation followed by a Bogoliubov rotation where the Bogoliubov angle and the dispersion are given by Note that the above definition ensures that the function θ q is continuous within the Brillouin zone q ∈ [−π, π]. To satisfy the proper boundary conditions, the allowed values of the momenta are q k = 2π N k for R and q k = 2π N (k + 1/2) for NS, respectively, with k = −N/2, . . . , N/2 − 1 and N even. The diagonalized Hamiltonian and its K-particle eigenstates are then given by It should be stressed that the eigenstates with K even belong to the spin-periodic Hamiltonian H + , whereas the eigenstates of the spin-antiperiodic H − have odd K. In the thermodynamic limit N → ∞, the periodic chain H + has a doubly degenerate ground state with ferromagnetic ordering along the x-axis, denoted by |⇑ and |⇓ , respectively. Note however, that for finite N the actual ground states in both symmetry sectors are given by which are separated by an exponentially small gap and both have vanishing magnetizations.

Magnetization dynamics
We are interested in the dynamics of the magnetization of various initial states, excited locally from the ferromagnetic ground state |⇑ and time-evolved under the Hamiltonian H in (1). The locality of the excitation is understood in terms of the Majorana basis, which implies that these excitations may become highly non-local in the spin-basis representation.
In fact, the latter will correspond to domain-wall excitations and one is interested in how the inhomogeneity spreads out under unitary time evolution. On the other hand, since the orderparameter magnetization is not conserved, even a single spin-flip excitation (which is local in terms of the spins) will lead to nontrivial dynamics. For the study of domain-wall melting, we will also consider for comparison a local quench setup where two separate chains are initially prepared in oppositely magnetized ground states, and subsequently joined together. The time-evolved magnetization can be extracted in a number of different ways. On the numerical side, we apply matrix product state (MPS) calculations [30,31] in an open-chain geometry. To ensure that we obtain the proper ferromagnetic (symmetry-broken) ground state |⇑ , we introduced a small longitudinal field h x > 0 in the Hamiltonian H − h x i σ x i for the first few sweeps and set h x = 0 afterwards, until convergence is reached. The excitations are then created by acting with the matrix product operator representation of the corresponding spin-excitation. Finally, the time evolution was implemented with the finite two-site timedependent variational principle (TDVP) algorithm [32].
On the other hand, we also employed Pfaffian techniques for the numerical evaluation of the magnetization. For the simple domain-wall excitation these were described in Ref. [25], but the calculations can easily be generalized for the other local excitations we deal with. In all of the examples we observed a perfect agreement with the results of MPS calculations.
Finally, we also present analytical results based on form-factor calculations. To this end, one has to first express the excited initial state |ψ 0 = (|ψ 0 R +|ψ 0 NS )/ √ 2 in the fermion basis, which is then time-evolved with the corresponding Hamiltonian in both symmetry sectors as Once |ψ 0 R/NS is written as a linear combination of the K-particle eigenstates (7), the time evolution is trivial e −itH R/NS |q 1 , q 2 , . . . , q K R/NS = e −it K k=1 ǫq k |q 1 , q 2 , . . . , q K R/NS , since the Hamiltonian H R/NS is diagonal in this basis. It is useful to introduce the normalized magnetization which can be evaluated as Note that, since the operator σ x n changes the parity of the state, the only non-vanishing contribution to the expectation value is between different parity sectors. In turn, the calculation of M n (t) boils down to evaluating multiple sums over the momenta with the form factors R p 1 , . . . , p L |σ x n |q 1 , . . . , q K NS , which are known explicitly from previous studies [33][34][35]. In the following we always consider the thermodynamic limit N → ∞, where the sums over momenta can be turned into integrals and the expressions for the form factors are summarized in Appendix A.

Single domain wall
Our first example is a single domain wall, which has already been considered for the TI [25] as well as for the XY chains [26]. For completeness, we revisit here the results obtained previously for the normalized magnetization. The single domain wall is an excitation |ψ 0 = D n 1 |⇑ created by the operator As remarked before, D n 1 is strictly local in terms of the fermions, whereas in the spin representation it creates spin-flips all over the sites j < n 1 . In the eigenbasis of the Hamiltonian it corresponds to a linear combination of one-particle states where we have suppressed the subscripts R/NS of the symmetry sector for notational simplicity. One thus only needs the form factors between one-particle states, which has a relatively simple form (51) given in Appendix A. Performing the time evolution (9) via (10) and inserting the result into (11), one arrives at The above expression simplifies considerably in appropriate scaling limits. Indeed, noting that the integral receives the dominant contribution due to a pole at q = p in the integrand of (14), one can change variables as Q = q − p and P = (q + p)/2, and perform a stationary phase analysis as described in Appendix B. In turn, one obtains which is the so-called hydrodynamic scaling limit. Here Θ(x) is the Heaviside step function, v P = dǫ P dP is the group velocity of the single-particle excitations and ν is the ray variable, with the distance measured from the initial location n 1 − 1/2 of the domain wall. The result (15) has a simple semiclassical interpretation, which has been applied many times to understand front dynamics in quantum chains [36][37][38][39]. Namely, the magnetization is transported by single-particle excitations, each carrying an elementary spin-flip, which contribute to the hydrodynamic profile at a given ray only if their velocity v P > ν.
Another interesting scaling regime emerges around the edge of the front ν ≈ v max , given by the maximum speed of excitations. In order to understand the fine structure of the edge, a higher order stationary phase analysis has to be performed around the momentum q * which yields the maximum velocity v q * = v max . As shown in Appendix B, this leads to the following result In other words, with the proper choice of the scaling variable X measuring the distance from the edge, and after appropriate rescaling, the fine structure of the magnetization front is given via the function Note that ρ(X) is nothing else but the diagonal part of the Airy-kernel K Ai (X, Y ) [40], which appears in a number of front evolution problems related to free-fermion edge universality [8,11,[41][42][43][44][45][46].
The results (15) and (16) have already been tested against numerical calculations for various parameters of the XY chain, where the notable feature of a hydrodynamic phase transition at h c = 1 − γ 2 was observed [26]. Indeed, this phase transition can be understood by the appearance of a second local maximum in the group velocities v q for h < h c , which in turn leads to kinks in the bulk of the hydrodynamic magnetization profile [26].
Finally, it should be noted that the analytical result was obtained by following the time evolution of one-particle states building up the domain wall. Strictly speaking, these states are eigenstates of H − only, i.e. the time evolution has to be performed with antiperiodic boundary conditions on the spin chain. However, since the form factor calculations are carried out directly in the thermodynamic limit, the boundaries actually do not play any role.

Double domain wall
We now move on to consider more complicated excitations, that are created by acting with the operator where n 2 > n 1 is assumed. In terms of fermions this is a two-local operator, i.e. supported on two sites only. In contrast, D n 1 ,n 2 is again nonlocal in the spin representation, and it is easy to see that it describes a double domain wall, located at sites n 1 and n 2 , respectively. Using (5), the excited initial state can be written as We shall restrict ourselves to the case n 2 − n 1 ≫ 1, i.e. when the two domain walls are spatially well separated, such that the sum in the first term of (19) becomes highly oscillatory and can be neglected. The initial state then involves only two-particle excitations and the time evolved state can be written as The magnetization M n (t) can thus be expressed as a quadruple integral via two-particle form factors R p 1 , p 2 |σ x n |q 1 , q 2 NS , that are reported in (53) in Appendix A. The result can be simplified, similarly to the single domain wall case, by analyzing the pole-structure of the form factors combined with a stationary phase approximation. The poles appear for momenta satisfying q 1 = p 1 and q 2 = p 2 or q 1 = p 2 and q 2 = p 1 . For the first pole one obtains two independent stationary phase conditions where P i = (q i + p i )/2 for i = 1, 2. Note that this pole corresponds to a process where the incoming momenta are matched with the outgoing ones at each domain wall separately. In contrast, at the second pole an incoming momentum of the first domain wall must match with an outgoing momentum of the second domain wall. However, as shown in Appendix B, after the exchange of the outgoing momenta and under the assumption n 2 − n 1 ≫ 1, the stationary phase condition cannot be satisfied. Thus only the first pole gives a contribution to the integral and leads to the result The hydrodynamic scaling limit of the profile in (22) has thus a factorized form with again a very simple physical interpretation. The ray variables ν i now measure the distances from the corresponding initial domain wall locations n i − 1/2, where quasiparticles with velocity v P i are emitted, each carrying a spin-flip. If, for a given pair of particles, one has v P 1 > ν 1 and v P 2 > ν 2 then both of the particles have reached site n at time t, hence the spin is flipped twice and one has a positive contribution. If, on the other hand, v P 1 < ν 1 and v P 2 > ν 2 , then only one particle has arrived and the contribution is negative. The profile is then obtained by summing the contributions over all pairs. In Fig. 1 we show the results of our MPS simulations together with the result (22). One can see a perfect agreement, even after the two fronts propagating from different locations overlap in the middle. In particular, one observes the emergence of two cusps at the ends of the overlap region, which follows from the factorized form of (22), i.e. one multiplies two single domain wall front profiles, each having square-root singularities at their edges. Moreover, this also implies that the outer edge of the front is still described by the same scaling (16) as for the single domain wall. On the right of Fig. 1 there are extra kinks to be seen, which is due to the fact that one has h < h c there, i.e. one is beyond the hydrodynamical phase transition point.

Single spin-flip
After having discussed the evolution of domain walls, we now study a very simple excitation, in the form of a single flipped spin. Naively, one would think that this excitation has a trivial hydrodynamic limit, and the flipped spin just disperses. However, since the magnetization is not conserved under the XY dynamics, it turns out that the profile is far from being trivial. In fact, the operator that creates a spin-flip at site n 1 is just σ z n 1 = −ia 2n 1 −1 a 2n 1 , which is strictly local in the spin representation, but is again two-local, i.e. a product of two adjacent Majoranas in the fermionic picture. Hence, this form is more reminiscent of a double domain wall excitation, with the exception that they are now created at neighbouring sites. Rewriting the excitation in the fermionic basis one has where the ground-state contribution is now proportional to the transverse magnetization and thus cannot be neglected. The calculation of M n (t) follows the same steps as in the previous cases. Note, in particular, that the two-particle contribution in (23) has almost the same form as (19) for the double domain wall with n 2 = n 1 + 1, except for the sign of the Bogoliubov phases. After time evolving and taking the expectation value with |ψ t , one has now cross terms where the form factors R 0|σ x n |q 1 , q 2 NS appear, see (52). However, since they have no poles, it is easy to see that their contribution is negligible in the scaling limit we are interested in. On the other hand, the two-particle form factors now yield a contribution from both of the poles. Indeed, the stationarity condition is, up to the sign of the θ ′ P i term, is the same as (21) for the double domain wall with n 2 = n 1 + 1. However, in the limit of t ≫ 1 and |n − n 1 | ≫ 1, the two equations are essentially the same. Hence, the process in which an incoming momentum of the first domain wall scatters into an outgoing momentum of the neighbouring one is equally well permitted and yields a sizable contribution.
Carrying out the stationary phase analysis in detail (see Appendix B), one arrives at the following result in the hydrodynamic limit where the ray variableν = n−n 1 t is slightly changed compared to (15), since the distance is now measured from the location n 1 of the spin-flip. The profile can be written as the sum of three terms, where the first one is simply the ground-state contribution. The second one corresponds to the factorized result for the double domain wall and the third one describes a kind of interference term, where the momenta of the excitations building up the two domain walls are exchanged. There is no simple semiclassical interpretation of this interference term, since the quasiparticles contribute with a phase factor. The result (25) is compared against our numerical calculations in Fig. 2 with an excellent agreement.
It is also interesting to have a look at the edge behaviour of the profile. Performing the higher order stationary phase analysis (see Appendix B), one is led to the following result where the scaling function is given bỹ  The result is thus very similar to the one for the domain wall in (16), however the scaling functionρ(X) acquires a nontrivial prefactor, which depends explicitly on the transverse magnetization m z , and even on the Bogoliubov phase evaluated at q * where the quasiparticle velocity has its maximum. In particular, this phase factor vanishes for the TI chain and one has a factor of 2 difference with respect to ρ(X). This explains the numerical findings of Ref. [24] where the very same setup was studied. We checked the validity of the edge scaling (26) in Fig. 3 for various parameter values and found a very good agreement, there are however some differences in the convergence towards the scaling functionρ(X).

Local quench
As a final example, we show here the results for the magnetization profile resulting from a local quench. That is, instead of applying a local excitation to the symmetry-broken ferromagnetic state, we rather prepare the two halves of our chain in oppositely magnetized ground states and join them together. Our goal is to check whether this protocol yields a similar result for the hydrodynamic profile as the one found for the single domain wall excitation.
The initial and time-evolved states are now given by Since our initial state is not prepared as an excitation upon the bulk vacuum state, it is a nontrivial question how |ψ 0 can be written in the basis of the full Hamiltonian H. Thus we shall only perform numerical (MPS and Pfaffian based) calculations for the quench. The results, shown in Fig. 4, turn out to be rather surprising. Namely, we find that in the TI limit (γ = 1) the profiles after the local quench (full symbols) almost exactly coincide with the ones for the domain wall excitation (empty symbols). The only deviations visible at the scale of the figure are around the front edges. In sharp contrast, for γ = 0.5 one has a huge deviation between the profiles for all the values of h we considered. This signals that in the latter case the factorized initial state is not well approximated by a single-particle excitation in the fermionic basis. We observe that the mismatch between the profiles gradually increases as one moves away from the TI limit. However, we have no clear explanation of this phenomenon which needs further studies.

Correlation functions
The form-factor approach is not restricted to the study of the magnetization profile. The next simplest physically interesting observable is the correlation function between the spins. Here we shall restrict ourselves to equal-time correlations between the x-components of the spin, which have already been addressed briefly in [26]. It is useful to work with the normalized correlation functions C m,n (t) = NS ψ t |M mMn |ψ t NS , where the expectation value is now taken between the NS components only, since the operator σ x m σ x n does not change the parity. Note that we use here that the corresponding expectation value between the R components is equal to (29) in the thermodynamic limit.
In order to get a form-factor expansion of (29), we shall insert the resolution of the identity Note that the resolution must be taken within the R sector, but we omit here the subscripts for notational simplicity. The form-factor expansion can be obtained by inserting the expression of |ψ t NS in terms of the fermionic basis. We focus here on the case of a single domain wall, since the calculations become rather cumbersome for more complicated excitations. In this case |ψ t NS is a superposition of single-particle states only and it is reasonable to assume that, for distances much larger than the correlation length |n − m| ≫ ξ, the dominant contribution to the correlations comes from the single-particle terms in (30) as well. To lowest order in the form-factor expansion we thus arrive at the result The hydrodynamic limit of (31) can be obtained in a similar fashion as was done for the magnetization profile. Expanding around the poles of the integrand and using the properties of the Θ function (see Appendix C for details) one obtains where the ray variables are measured from the initial domain wall location and the expression has a very simple interpretation. Let us assume ν > µ and consider the contribution of a single quasi-particle traveling at speed v P . Now, for short times v P < µ the excitation has not yet reached the first spin and thus the correlations are ferromagnetic. Once µ < v P < ν, the first spin has been flipped while the second one is still untouched, hence the correlation is antiferromagnetic. Finally, after the excitation has traveled through, v P > ν, the second spin is also flipped and the correlation becomes ferromagnetic again. It turns out that, instead of approximating the integrals in (31), there is a way to directly relate C m,n (t) to the profile M n (t). Indeed, by turning the integral over p into a contour integral and applying the residue theorem, one obtains the formula (80) reported in Appendix C, which is an exact relation at the level of one-particle form factors. However it is easy to see that, similarly to the hydrodynamic approximation in (32), it yields perfect ferromagnetic correlations C m,n (t) ≃ 1 when both spins are outside the front region. Indeed, it can be shown that the many-particle form factors are the ones responsible for the exponentially decaying correlations C 0 m,n in the ground state [34]. One can thus reincorporate these correlations into the approximation as The relation in (34) is tested against exact numerical calculations for the TI chain in Fig. 5. We have calculated the correlations along the front region while keeping the distance d between the spins fixed. One can see that, for d = 1, there is still a slight deviation from (34) which, however, decreases with increasing d. For d = 10 one has already an excellent agreement with no visible deviations. In fact, for |n − m| ≫ ξ one has C 0 m,n → 1, and one recovers the one-particle result (80) which should become exact. Note, however, that calculating the corrections to (34) would require to evaluate multiple integrals with higher-order offdiagonal form factors and is thus a difficult task. Nevertheless, a closer investigation of the form-factor structure in (48) confirms, that the dominant pole contribution is suppressed and thus one indeed obtains subleading terms.

Entanglement dynamics
So far we have studied the simplest observables. One can, however, gather important information about the time-evolved state by looking at the entanglement dynamics. In particular, we are interested in the entanglement profiles along the front region, considering a bipartition into two disjoint segments A = [1, N/2 + r] and its complement B, and calculating the resulting von Neumann entropy. Entanglement profiles for domain-wall type initial conditions have been studied extensively for time evolution under critical Hamiltonians [13,15,17,22,24,[47][48][49][50], and even a description in terms of CFT has been given [51,52]. However, much less is known about the non-critical case, such as the one at hand. The calculation of the entanglement profile is straightforward in MPS calculations, however, extracting the entropy via covariance-matrix techniques for Gaussian states [53,54] requires some extra considerations. Indeed, the problem lies in the nature of the initial state, since the excitations are created upon the symmetry-broken ground state, which is inherently non-Gaussian [55]. Nevertheless, this difficulty can be overcome by relating the problem to the one where the very same excitations are created upon the Gaussian, non-magnetized ground states in (8). The method has already been outlined in [26] but we expand here the arguments for completeness.
Let us consider initial states corresponding to the two symmetry-broken ground states of the system. Using (8), the density matrices are given by where the even and odd parity components, satisfying [P, ρ e ] = 0 and {P, ρ o } = 0, respectively, are defined as Clearly, the problem is with the odd component ρ o , since a Gaussian density operator is by definition even. One can, however, eliminate ρ o by considering an equal-weight convex combination of the density matrices in (35). The resulting density matrix ρ e is itself still a convex combination of two Gaussian states from the NS and R sectors. However, working in the thermodynamic limit, these two states become indistinguishable [55], and one concludes that ρ e is equivalent to a proper Gaussian state. Furthermore, as shown in Ref. [56], excitations that can be written as a product of Majorana fermions where J is an arbitrary index set, preserve Gaussianity. So does unitary time evolution governed by a quadratic Hamiltonian. Hence, introducing the notation for the reduced density matrices of a given bipartition, after exciting and time evolving the initial states in (35), we finally come to the conclusion that is a well-defined Gaussian state living on the Hilbert space of segment A. Our goal is now to relate the entropy S(ρ ⇑ A ) = −Tr ρ ⇑ A ln ρ ⇑ A of our target state to that S(ρ A ) of the Gaussian state in (39). To this end we use the inequality for convex combinations of density matrices [57,58] First, we note that from trivial symmetry arguments one has S(ρ ⇓ A ) = S(ρ ⇑ A ). Furthermore, it is also known [58] that the inequality (40) is saturated if the ranges of ρ i are pairwise orthogonal, which is again clearly satisfied in our case due to ⇑ | ⇓ = 0. Hence one finds Finally, it remains to calculate the covariance matrix Γ A corresponding to ρ A , from which the calculation of the entropy S(ρ A ) follows standard procedure [53,54]. Since ρ A is the reduced density matrix of the time-evolved and excited ground state |ψ t NS , Γ A is just the reduction of the full covariance matrix with elements Γ k,l = NS ψ t | [a k , a l ] |ψ t NS /2. This can be obtained by working in the Heisenberg picture. Since D J is unitary, D J D † J = 1, the effect of the excitation can be absorbed by a change of the Majorana basis [56] The orthogonal transformation Q has a simple diagonal matrix form with entries ±1, depending on whether the corresponding column is part of the index set J or not. In complete analogy, the unitary time evolution corresponds to the basis rotation where the explicit form of the orthogonal matrix R was reported in Ref. [25]. Putting everything together, one finds that where Γ 0 is the ground-state covariance matrix with elements (Γ 0 ) k,l = NS 0| [a k , a l ] |0 NS /2.
We are now ready to discuss the entanglement dynamics for the simple excitations introduced in Sec. 3. In each case we have verified that the entropy obtained by the procedure outlined above agrees perfectly with the results of our MPS calculations.

Single domain wall
The entropy profiles for the single domain wall, located initially in the center (r = 0) of the chain, have already been considered in [26] and are shown in the left of Fig. 6 for γ = 0.5 and several values of h. The profile ∆S(r) = S(ρ ⇑ A ) − S 0 is always measured from the initial entropy S 0 of the bulk ferromagnetic state, and is plotted against the rescaled distance ζ = r/t from the center of the chain. The main feature to be seen is the emergence of a kink in the profile for h < h c , at the value ζ * that equals the local maximum of the quasiparticle velocity, in complete analogy to the case of the magnetization.
Due to the similar features observed in the entropy and magnetization profiles, one is naturally led to the question whether there is a simple relation between the two of them. We are also motivated by recent results of Refs. [27,28], where the entanglement content of particle excitations in 1 + 1-dimensional massive quantum field theories was studied, with a surprisingly simple result. Namely, it has been found that the entropy difference (relative to the ground state) of a single-mode excitation is independent of the wavenumber and given by the binary entropy formula involving the ratio of the subsystem and full system lengths [27,28]. This ratio is just the density fraction of the single-mode excitation that is contained within the subsystem.
Inspired by these findings, we put forward the following ansatz In other words, we assume that the static results of [27,28] would generalize to our dynamical scenario, and the entropy difference for bipartitions along the ray ζ is just given by the same binary formula, with the density ratio N (ζ) being the fraction of the quasiparticles that have reached the entangling point. Surprisingly, we find that the simple-minded ansatz (46), shown by the red solid lines in the left of Fig. 6, gives a very good description of the entropy profiles. Via the density fraction N (ζ), the entropy profiles are thus directly related to those of the magnetization (15). In case h < h c , one observes some deviations from the ansatz (46), which are only visible in the regime ζ < ζ * and are assumed to be finite-time effects. In order to better understand the convergence, on the right of Fig. 6 we also studied the time evolution of the half-chain entropy ∆S(0), for the same parameter values. Although each of them can be seen to converge towards the asymptotic value ln 2, their approach is rather different. For h > h c the convergence is fast and steady, with rapid oscillations only, whereas for h < h c there is a smaller frequency appearing with a larger amplitude, and the curve bounces back from its asymptotical value repeatedly. Interestingly, at the critical point h = h c = 0.75 one can see a slowing down in the convergence, which becomes most evident on a logarithmic scale as shown on the inset of the figure. Indeed, the approach seems to be a power law t −1/2 , as opposed to t −1 in the h = h c case. This critical slowing down is responsible for the dip around ζ = 0 in the profile for h = h c on the left of Fig. 6.
One should stress the marked difference of the entropy profiles as compared to domainwall evolution in critical systems, such as the XX chain. Indeed, in the latter case the entropy was found to grow logarithmically in time in the entire front region [48,52], whereas here the profiles converge to the scaling function (46) when plotted against ζ = r/t. In particular, the result ∆S(0) = ln 2 for ζ = 0 implies that the entropy converges towards the value attained in the ground state |0 NS , which has been studied in [59,60]. Indeed, applying the relation (41) at t = 0, one finds that the entropy S 0 in the initial symmetry-broken ground state is exactly ln 2 less than that of the NS ground state. This strongly suggests that the steady state is nothing but the ground state with its symmetry restored.

Double domain wall
The profiles for the double domain wall are shown in Fig. 8 for various times and two different model parameters. In both cases, the profiles resemble those of two separate single domain walls for short times, while for large times the main feature is the emergence of an additional plateau in the overlap region. This strongly suggests the relation ∆S n 1 ,n 2 (r) = ∆S n 1 (r) + ∆S n 2 (r) , where ∆S n 1 ,n 2 (r) and ∆S n i (r) denote the entropy differences for double and single domain walls, respectively, with the indices referring to the initial locations of the excitations. In other words, one expects the entropy differences to be additive, which is indeed perfectly confirmed by the numerics.

Single spin-flip
Finally, we consider the entropy profiles for the spin-flip excitation, with the results shown in Fig. 8, for the same choice of parameters as for the magnetization profiles in Fig. 2. When plotted against the scaling variable ζ, the profiles show a different behaviour as compared to those of the single domain wall excitation in Fig. 6. In particular, the additivity (47) is not satisfied, analogously to the corresponding result (25) for the magnetization, which does not have a factorized form. Indeed, as explained under Sec. 3.3, this has to do with an interference effect in the dynamics, where an incoming momentum of the first excitation can travel forward as an outgoing momentum of the second one. Clearly, such a process creates entanglement between the quasiparticles building up the two domain-wall excitations, which spoils the additivity and reduces the overall entropy of the state. Unfortunately, despite the qualitative understanding of the origin of the nontrivial entropy behaviour, we have not been able to find an ansatz analogous to (46) that captures the profiles quantitatively.

Discussion
We studied the time evolution of the magnetization and entanglement profiles in the XY chain for simple initial states that can be written as a product of one or two local fermionic excitations. The former corresponds to a single domain wall in the spin-picture and the magnetization profile has a simple hydrodynamic limit (15), corresponding to the motion of independent quasiparticles. Furthermore, in the very same limit we find that the entropy is given by the simple ansatz (46) and is thus directly related to the magnetization profile. The correlation function is also found to be related via (29) to the magnetization, which gives a very good approximation even for finite times and distances. For double domain walls, excited by the product of two fermions separated by a large distance, we simply find the factorized form (22) for the magnetization, as well as the additivity (47) of the entropy differences. For a single spin-flip, however, the fermions are located on neighbouring sites and the excitation cannot be considered strictly local any more. As a consequence, we find an interference term in the magnetization profile (25). Furthermore, the additivity of the entropy is lost, and we find convergence towards a nontrivial profile We have also compared the profiles for the single domain wall to the ones obtained via a local cut and glue quench, where the two ferromagnetic ground states are prepared on halfchains and joined together. Rather surprisingly we found that, while for the TI chain the respective profiles almost coincide, for the generic XY case they become completely different (see Fig. 4), with the discrepancy growing with the distance from the TI limit. Apparently the local quench is well approximated by a single fermionic excitation for the TI but not any more for the XY case. The precise origin of this phenomenon is unclear to us and requires further studies.
It would be also interesting to see if a QFT treatment of the entropy increase could be given. Even though our ansatz (46) was inspired by QFT results [27,28] on the entanglement content of particle excitations, those particles are single wave modes and there is no dynamics involved. On the other hand, for the case of critical Hamiltonians there exists a CFT framework for calculating the time evolution of entropy after spatially local excitations [61]. Whether this approach could be generalized to a massive QFT to predict the asymptotic entropy increase after the excitations considered in this paper is a puzzling question to be addressed.
One could also think about extending the studies to excitations composed of a product of more than two fermions. While being a straightforward generalization, the form-factor calculations are likely to be very cumbersome, due to the increasing number of the pole contributions one has to account for. Finally, it is natural to ask how the setup could be extended to interacting integrable systems, and if the treatment of such composite but still essentially local excitations could be incorporated into the theory of GHD.

A Form factors for the TI and XY chains
Here we present the form factors used in the calculations of the main text. Although for our simple excitations we required only few-particle form factors, the general expression is reported for completeness. The formula is rather involved even after taking the thermodynamic limit N → ∞, and for the TI chain (γ = 1) it reads [33] R p 1 , . . . , p L |σ x n |q 1 , . . . , q K NS R 0|σ x n |0 NS We have assumed here that the number of momenta K and L on the right and left hand side have the same parity, otherwise the form factor vanishes. Note that we have normalized with the vacuum form factor, i.e. with the expectation value of the ground-state magnetization. For K = L the form factors (48) depend only on the dispersion relation ǫ q , given in Eq. (6), and the values of the momenta. For the generic case of the XY chain, the expressions become even more complicated. In the limit N → ∞ they can be written as [34,35] R p 1 , . . . , p L |σ x n |q 1 , . . . , q K NS R 0|σ x n |0 NS where we have defined The above definition is valid in the parameter regime 1 − γ 2 < h < 1, i.e. in the nonoscillatory ferromagnetic phase. In the oscillatory phase 0 < h < 1 − γ 2 the corresponding expressions can be obtained by analytic continuation [34]. One can also check that, in the singular TI limit γ → 1, the expression (49) goes over to the one in (48). While in general they differ in the details, these will turn out to be irrelevant for the hydrodynamic limit, since their pole structure is exactly the same. We now discuss the form factors needed in the main text. The simplest is the one-particle form factor (K = L = 1), where using some hyperbolic identities in (49), one can show that the TI and XY cases yield the same expression Thus the formula (14) for the single domain wall excitation is valid for arbitrary parameter values of the XY chain. In general, no such simplification occurs and in the following we restrict ourselves to the TI case for the sake of simplicity. For the spin-flip excitation one needs the off-diagonal form factor with K = 2 and L = 0 which reads One can see immediately, that this form factor does not have any poles which implies that it will only give a subleading contribution. The diagonal two-particle form factor (K = L = 2), on the other hand, has the form with two possible poles for q 1 = p 1 and q 2 = p 2 , or with an exchange of momenta for q 1 = p 2 and q 2 = p 1 . It should be noted that, for the generic diagonal K-particle form factors in (48), an arbitrary permutation between the incoming and outgoing momenta yields a pole, which makes the analysis of the contributions increasingly complicated.

B Stationary phase calculations for the profile
In this appendix we summarize the calculations leading to the approximations of the magnetization profile in the hydrodynamic regime. The simplest case is the single domain wall, where M n (t) is given by a double integral (14). The integrand has a pole due to the form factor, which can be regularized as M n (t) = 1 + π −π dp 2π π −π dq 2π by introducing the infinitesimal shift ε > 0. The integrand of (54) is highly oscillatory for |n − n 1 | ≫ 1 and t ≫ 1, and the location of the pole at q = p suggests the change of variables Q = q − p and P = (q + p)/2. The phase factors become stationary at Q = 0, thus the integrand should be expanded around this value. Keeping only the most singular term one has where we have extended the integration in the relative momentum up to infinity. Thanks to the definition (6), the function θ ′ P varies smoothly and one can neglect it in the hydrodynamic regime. Then using the integral representation of the Heaviside theta function and introducing the ray variable ν = (n − n 1 + 1/2)/t brings us to the result (15) in the main text.
The bulk hydrodynamic profile is thus recovered by solving the equation v q = ν. Special attention is needed around the maximum v q * = v max of the velocities, where the solutions coalesce at momentum q * . To get the fine structure of the front edge, one has to expand the dispersion around q * as Furthermore, one can introduce the following rescaled variables Substituting (57) and (58) into (54), one arrives at the following integral Using the integral representation of the Airy kernel [40] K Ai (X, Y ) = lim ε→0 dP 2π dQ 2π one recovers (16) of the main text, with ρ(X) = lim Y →X K Ai (X, Y ) given by the diagonal terms of the Airy kernel.
The hydrodynamic limit (22) for the double domain wall can be obtained in a similar fashion, however, one has now a quadruple integral to start with. The poles are contained in the two-particle form factor (53). First, we consider the pole with q 1 = p 1 and q 2 = p 2 . Changing variables as and expanding the phases around the stationary points Q i = 0, one has where we defined The function f in (62) describes the slowly varying part of the form factor in (53). It is easy to see, that the terms containing the dispersion ǫ q i and ǫ p i can be approximated by 1 to leading order. It remains to analyze the contribution of the trigonometric factors that do not contain the poles, which can be rewritten as Thus, again to leading order around Q i = 0, one has f (P 1 , P 2 , Q 1 , Q 2 ) ≈ −1 + O(Q 1 Q 2 ), meaning that the first correction would already remove the singularity in the integral (62), and can be neglected. Setting f = −1, one recovers immediately the factorized result (22). The second pole of the form factor (53) is given by q 1 = p 2 and q 2 = p 1 and corresponds to an exchange of the outgoing momenta. The form factor itself transforms trivially under this exchange, acquiring only a sign. The time-evolved state (20), however, has phase factors attached to the locations of the domain walls and thus transforms nontrivially under exchange of the momenta. Indeed, introducing the variables this phase factor can now be rewritten as The second term contains the center of mass momenta and becomes highly oscillatory for |n 2 − n 1 | ≫ 1. This phase, however, cannot be made stationary, since the time-dependent part of the phase in (20) is symmetric under the exchange of the momenta. One thus concludes that, for large separations of the domain walls, the second pole gives a negligible contribution. The situation for the spin-flip excitation is different. As discussed in the main text, except for a sign change of the Bogoliubov angles, the state (23) is a double domain wall with n 2 = n 1 + 1. The first pole thus yields the very same factorized result as in (62), with the corresponding changes in x i . In the hydrodynamic limit, however, it is more natural to measure distances from the spin-flip location n 1 (instead of n 1 ± 1/2) and use the ray variableν = n−n 1 t , which gives the second term in (25). The second pole, however, has also a significant contribution, since n 2 − n 1 = 1 and the phase factor in (66) now varies slowly. Expanding around Q ′ i = 0, one finds where x ′ i = v P ′ i t − (n − n 1 ) and the sign change in the form factor has been taken into account. It is easy to see that Regularizing the pole via the identity Q ′−1 = iπδ(Q ′ ) + lim ε→0 (Q ′ + iε) −1 , using (56) and the expression of the transverse magnetization in (24), the third term of (25) follows. It remains to investigate the edge scaling regime for the spin-flip excitation. The second term of (25) is simply the square of the profile for a single domain wall, where the edge scaling is given by (16). To leading order, this just yields a factor 2. The situation is similar for the third term in (25) where, additionally, the phase factors in the integral must be evaluated at the momentum q * where the velocity has its maximum, v q * = v max . Collecting the terms, one obtains the prefactor in (27).
Finally it should be noted that, although the calculation above has been carried out using the form factors for the TI chain, the result generalizes to the XY case. Indeed, the pole structure of the form factors is exactly the same, whereas the differences in the slowly varying part are irrelevant in the hydrodynamic limit, since they have the same trivial limit after expanding around the pole.
The stationary phase approximation of this integral is very similar to that of the magnetization profile. Introducing the new set of variables and expanding around the poles Q 1 = 0 and Q 2 = 0, one obtains C m,n (t) ≃ 4 dP 2π dQ 1 2π Applying (56) in both the Q 1 and Q 2 integrals, the result can again be written with the help of step functions where the scaling variable µ = (m − n 1 + 1/2)/t is introduced analogously to ν. Assuming µ < ν and using the identities for the step function the result (32) of the main text follows immediately.
Instead of applying the stationary phase argument, one can also do a more precise analysis. Indeed, it turns out that the integral over p in (69) can be carried out explicitly. We first regularize the factor containing the pole as (74) Multiplying out this expression, the terms containing the delta functions can be plugged back into (69) and integrated over. Comparing to (54), one can identify the resulting double integrals as M m (t) − 1 and M n (t) − 1, respectively, while the product of the delta functions trivially yields one. The remaining factor from (74) can be rewritten as the integral over p is transformed into the contour integral where f (z) is the slowly varying regular part of the integrand in (69), and the contour is the unit circle. Now the two poles are located at z = z 0 and z = z −1 0 . However, for ε > 0, only z = z 0 lies inside the contour and contributes to the integral. We have thus to obtain the residue around this pole. Rewriting and the two poles correspond to p = q 1 and p = q 2 , respectively. Hence the result of the contour integral is Finally, noting that I enters with a minus sign (see (74)), and inserting the result back into (69), one can easily identify the contribution as −2(M n (t) − 1). Collecting all the terms, one arrives at the result C m,n (t) ≃ 1 + M m (t) − M n (t) .
As a closing remark, we give a simple argument why the many-particle contributions in the form-factor expansion of the correlation functions can be neglected. In the one-particle expression (69), the dominant contribution is obtained from momenta satisfying q 1 = p = q 2 , where the stationary phase conditions match the poles of the integrand. The next nonvanishing term in the expansion involves three intermediate particles, where the phase factor could be made stationary for q 1 = p 1 = q 2 and p 2 = −p 3 . However, from (48) one can see that there is no pole in the form factor at p 2 = −p 3 , and thus the contribution is subleading.