Entanglement spreading after local fermionic excitations in the XXZ chain

We study the spreading of entanglement produced by the time evolution of a local fermionic excitation created above the ground state of the XXZ chain. The resulting entropy profiles are investigated via density-matrix renormalization group calculations, and compared to a quasiparticle ansatz. In particular, we assume that the entanglement is dominantly carried by spinon excitations traveling at different velocities, and the entropy profile is reproduced by a probabilistic expression involving the density fraction of the spinons reaching the subsystem. The ansatz works well in the gapless phase for moderate values of the XXZ anisotropy, eventually deteriorating as other types of quasiparticle excitations gain spectral weight. Furthermore, if the initial state is excited by a local Majorana fermion, we observe a nontrivial rescaling of the entropy profiles. This effect is further investigated in a conformal field theory framework, carrying out calculations for the Luttinger liquid theory. Finally, we also consider excitations creating an antiferromagnetic domain wall in the gapped phase of the chain, and find again a modified quasiparticle ansatz with a multiplicative factor.


I. INTRODUCTION
The non-equilibrium dynamics of integrable models has developed into a vast field of research [1]. Among the numerous aspects, the understanding of local relaxation and equilibration in closed quantum systems has become a central topic of investigation [2,3]. In this respect, integrable systems show a rather peculiar behaviour, as the dynamics is characterized by the existence of stable quasiparticle excitations. This is intimately related to the extensive number of nontrivial conservation laws, which nevertheless allow for a local relaxation in a generalized sense [4].
Starting from the early studies of this topic, it was identified that the spreading of entanglement must play a key role in our understanding of integrable dynamics. Ground states of homogeneous, local Hamiltonians have a low amount of entanglement, typically satisfying an area law [5]. However, considering the time evolution with respect to a different Hamiltonian as in the context of a global quantum quench [6], the rapid linear growth of entanglement was attributed to the ballistic propagation of entangled quasiparticle pairs [7]. These quasiparticles transmit entanglement over large distances, contributing to the buildup of an extensive entropy within any given subsystem, which signals the onset of some local thermalization. Specifically, in one-dimensional integrable chains it has been verified that the entanglement entropy accumulated in a subsystem actually plays the role of the thermal entropy as described by the generalized Gibbs ensemble [8][9][10].
The global quench is the simplest representative of an initial state that has an extensive amount of energy above the ground state of the Hamiltonian governing the dynamics, thus acting as a reservoir of quasiparticle excitations. The interpretation, however, becomes more complicated if the initial state lies in the low-energy regime. A particular example is the local quench, where the final Hamiltonian is disturbed only locally with respect to the initial one, such as joining two initially separated quantum chains. At criticality, the entanglement spreading can be captured via conformal field theory (CFT) [11][12][13], predicting a slow logarithmic growth of the entropy, which was indeed observed in free-fermion chains [14]. However, despite signatures of the underlying quasiparticle dynamics, such as a light-cone spreading with the maximal group velocity, it is unclear how the individual quasiparticles contribute to the entropy.
Yet another situation that has been studied intensively within CFT is the so-called local operator excitation [15][16][17]. Here the low-energy initial state is excited from the vacuum of the CFT by the insertion of a local primary operator, while the Hamiltonian is left untouched. The disturbance has then a linear propagation, increasing the entanglement of a segment only while passing through it, with a constant excess entropy determined by the quantum dimension of the local primary [15][16][17]. The calculations have been extended in various directions, considering fermionic [18] or descendant fields [19,20], multiple excitations [21], as well as the effects of finite temperatures [22] or boundaries [23].
Despite this increased attention, there have been much less studies on entanglement spreading after local excitations in integrable quantum chains. The CFT predictions have been tested on the critical transverse Ising [24] and XX chains [25], for various local operators that are lattice analogs of primary or descendant fields. On the other hand, entanglement spreading has also been considered in the non-critical ordered phase of the Ising [26] and XY chains [27,28], starting from a domain-wall initial state excited by a local Majorana operator. Remarkably, the emerging profile of the excess entropy was shown to be captured by a simple probabilistic quasiparticle ansatz [28]. Indeed, taking into account the dispersive spreading of quasiparticles, only a certain fraction of the initially localized excitation will cross the subsystem boundary located at a certain distance. Interpreting this quasiparticle fraction as the probability of finding the excitation within the subsystem, the excess entropy is simply given by a binary expression [28].
Here we aim to extend the quasiparticle description of entanglement spreading to local fermionic excitations in the XXZ chain. Being a Bethe ansatz integrable interacting model [29,30], its quasiparticle content is much more complex than in the free-fermion systems considered so far. Nevertheless, since our local excitations probe the lowenergy physics, it seems reasonable that the dominant weight is carried by low-lying spinon excitations, which we shall assume to build our quasiparticle ansatz. Compared against the profiles of the excess entropy, as obtained from density-matrix renormalization group (DMRG) calculations [31][32][33], we observe a good agreement after a local fermion creation for moderate values of the interaction. For larger interactions in the gapless phase, one finds deviations that can be attributed to different types of quasiparticles with higher energy.
We also study the profiles after a local Majorana excitation, which seem to be given by a simple rescaling of the spinon ansatz. This result is supplemented by CFT calculations carried out for the Luttinger liquid theory, which describes the low-energy physics of the XXZ chain. We find that, due to the left-right mixing of the chiral bosonic modes, the asymptotic excess entropy is doubled for the Majorana excitation, although with a very slow convergence towards this value. Finally, in the gapped phase of the chain we study the excess entropy profile after a local Majorana operator that excites an antiferromagnetic domain wall. Here our numerical results suggest that the spinon ansatz is multiplied by a nontrivial factor, related to the ground-state entropy.
The rest of the manuscript is structured as follows. In section II we introduce the XXZ chain and discuss its low-lying excitations. Section III is devoted to the study of entanglement spreading after local excitations in the gapless phase: we first introduce a quasiparticle ansatz for the excess entanglement, followed by our numerical studies of a fermion creation as well as a Majorana excitation. Our results for the gapless regime are complemented by a calculation of the Rényi entropy within a CFT framework in section IV. Finally, in section V we consider entanglement and magnetization profiles after a domain-wall excitation in the gapped regime. Our closing remarks are given in section VI, followed by an appendix containing the details of the CFT calculations.

II. XXZ CHAIN AND LOW-ENERGY EXCITATIONS
We consider an XXZ chain of length L with open boundary conditions that is given by the Hamiltonian where S α j = σ α j /2 are spin-1/2 operators acting on site j, and ∆ is the anisotropy. The energy scale is set by the coupling J which we fix at J = 1. The XXZ Hamiltonian (1) conserves the total magnetization S z in z-direction and we will be interested in its ground state in the zero-magnetization sector S z = 0. Equivalently, the XXZ spin chain can be rewritten in terms of spinless fermions by performing a Jordan-Wigner transformation, which brings (1) into the form where c † j (c j ) are fermionic creation (annihilation) operators, satisfying anticommutation relations {c i , c † j } = δ ij . One then has a half-filled fermionic hopping chain with nearest-neighbour interactions of strength ∆. For |∆| ≤ 1 the system is in a critical phase with gapless excitations above the ground state, whereas a gap opens for |∆| > 1. The case ∆ = 1 corresponds to the isotropic Heisenberg antiferromagnet.
In the following we give a short and non-technical introduction to the construction of the ground state and lowlying excited states of the XXZ chain. To keep the discussion simple, we shall rather consider a periodic chain, and focus on the behaviour in the thermodynamic limit L → ∞. The exact eigenstates of the XXZ chain can be found from Bethe ansatz [29,30]. These are constructed as a superposition of plane waves, the so-called magnons, labeled by their rapidities λ i which provide a convenient parametrization of the quasimomenta. The allowed values of the rapidities follow from the Bethe equations, with real solutions corresponding to spin-wave like states. Complex solutions organize themselves into strings and correspond to bound states.
For |∆| < 1 the half-filled ground state is obtained by occupying all the allowed vacancies of the L/2 real rapidities, thus forming a tightly packed Fermi sea. Low-energy excitations in the S z = 1 sector are called spinons and are created by removing a rapidity. This creates two holes in the Fermi sea, with all the remaining rapidities moving slightly with respect to their ground-state values, and the energy difference can be calculated from this back-flow effect. In the thermodynamic limit, the result can be found analytically and written directly in terms of the quasimomenta q 1 and q 2 of the two spinons as [29] where the spinon dispersion relation in the gapless regime with ∆ = cos(γ) is given by Note that spinons are always excited in pairs, with the individual momenta confined to 0 ≤ q 1,2 ≤ π. The total momentum is then given by q 1 + q 2 , and due to the additivity of (3) one actually has a band of excitation energies. In particular, the lower edge of the two-spinon band is obtained by setting q 2 = 0 or q 2 = π, and thus simply corresponds to shifting the dispersion in (4) for q > π. The group velocity of the spinons can be directly obtained from the derivative of the dispersion Further low-energy excitations with S z = 1 can be created by removing a single rapidity from the real axis and placing it onto the Im λ = π axis. The energy of this particle-hole excitation can be obtained, similarly to the spinon case, from the back-flow equations of the rapidities and yields the dispersion [29] ε ph (q) = π sin(γ) However, in contrast to spinons, particle-hole excitations are not composite objects and their momentum range is thus 0 ≤ q < 2π. Note that these spin-wave like excitations are only physical for −1 < ∆ < 0, i.e. in case of attractive interactions. For low momenta q → 0, the dispersion relation Eq. (6) approaches the one for spinons in Eq. (4). The group velocities of particle-hole excitations are obtained by taking the derivative of ε ph (q). Interestingly, it was found that the maximum particle-hole velocity can exceed the maximum spinon velocity only if the anisotropy satisfies ∆ < ∆ * ≈ −0.3, which was demonstrated in a particular quench protocol [34]. Finally, we consider the gapped phase where we focus exclusively on the antiferromagnetic regime ∆ > 1, with the standard parametrization ∆ = cosh φ. For even L the ground state has S z = 0 and is again given by L/2 magnons with real rapidities. However, the allowed number of vacancies is now L/2 + 1, which allows to construct a slightly shifted Fermi sea. In the Ising limit ∆ → ∞, this yields an exact twofold degenerate ground state, given by the linear combinations of the two Néel states For finite ∆, the two states |ψ ± constructed this way are only quasi-degenerate, with an energy difference decaying exponentially in the system size L. Considering the thermodynamic limit one can write where |ψ ↑ and |ψ ↓ correspond to ground states with spontaneously broken symmetry, displaying antiferromagnetic ordering. In fact, the bulk expectation value of the staggered magnetization can be calculated analytically as [35,36] The low-lying excitations in the gapped phase are given again by spinons, by creating two holes in the Fermi sea. The excitation energy is still given by Eq. (3), with the dispersion in the gapped phase obtained as [29] where the complete elliptic integral of the first kind reads K(u) = π/2 0 dp 1 − u 2 sin 2 (p) (11) and the elliptic modulus u satisfies The spinon velocity is obtained from the derivative of (10) and reads v s (q) = sinh(φ) π K(u) u 2 sin(q) cos(q)

III. ENTANGLEMENT DYNAMICS IN THE GAPLESS PHASE
The goal of this section is to study the entanglement dynamics after a particular class of excitations. Namely, we first initialize the chain in its gapless ground state |ψ 0 , which is then excited by an operator that is strictly local in terms of the creation/annihilation operators c † j and c j appearing in the fermionic representation (2) of the XXZ chain. The system is then let evolve freely and we are interested in the emerging entanglement pattern in the time-evolved state |ψ(t) . For a bipartition into a subsystem A and the rest of the chain B, this is characterized by the von Neumann entropy with the reduced density matrix ρ A (t) = Tr B ρ(t) and ρ(t) = |ψ(t) ψ(t)|. In particular, we consider the bipartition A = [−L/2 + 1, r] and B = [r + 1, L/2] and study the entropy profiles along the chain by varying r, where r = 0 corresponds to the half-chain. Note that by subtracting the ground-state entropy S(0), we aim to extract information about the excess entanglement created by a local excitation.
In the following subsections we first introduce an intuitive picture for the description of the entanglement spreading in terms of the low-lying quasiparticle excitations introduced in Sec. II. We then proceed to the numerical study of the entanglement profiles after exciting the ground state with a fermionic creation operator, and compare the results to our quasiparticle ansatz. In the last part we consider an excitation created by a local Majorana fermion operator.
A. Entanglement spreading in the quasiparticle picture Let us consider an excitation above the ground state of the XXZ chain by acting with a fermion creation operator c † j . To capture the dynamics, one would have to first decompose the initial local excited state in the eigenbasis of the Hamiltonian. As discussed in the previous section, these eigenstates are described by quasiparticles parametrized by their rapidities or quasimomenta. The entanglement properties of various eigenstates in the XXZ chain were studied before in [37,38], whereas a systematic CFT treatment of low-energy excitations was introduced in [39,40].
In the framework of free quantum field theory, a surprisingly simple result on quasiparticle excitations was recently found in [41,42]. Namely, the excess entanglement measured from the ground state was found to be completely independent of the quasiparticle momenta, depending only on the ratio p of the subsystem and full chain lengths. Moreover, for quasiparticles described by a single momentum, the excess entropy is given by a binary formula ∆S = −p ln p − (1 − p) ln(1 − p), which allows for a simple probabilistic interpretation. Indeed, the ratio p is just the probability of finding the quasiparticle within the subsystem. Motivated by these results, we now put forward a simple ansatz for the spreading of entanglement after the local excitation. Under time evolution, the quasiparticles involved in the decomposition of the initial state spread out with their corresponding group velocities. However, our main assumption is that their contribution to entanglement is still independent of the momentum. Furthermore, we shall also assume that the dominant part of the entanglement is carried by the lowest-lying spinon modes, and that a spatially localized excitation translates to a homogeneous distribution of the momenta in the initial state. Under these assumptions we expect that the entanglement profile at time t 1 and distance r 1 from the excitation, in the space-time scaling limit ζ = r/t fixed, is determined exclusively via where Θ(x) is the Heaviside step function and v s (q) is the spinon velocity. In fact, this is nothing else but the fraction of the spinon modes with sufficient velocity to arrive at the subsystem. The simple probabilistic interpretation of the entanglement then leads to the binary entropy formula for the profile In particular, for the gapless case considered here, inserting the expression (5) of the spinon velocity into (16), the spinon fraction can immediately be found as where v = v s (0) denotes the maximal spinon velocity. In summary, our simplistic ansatz (17) provides an interpretation of the excess entropy based on the dispersive dynamics of the quasiparticle modes, where N is the fraction of the initially localized excitation that arrives at the subsystem. In fact, the very same ansatz has recently been suggested for the description of entanglement spreading after local fermionic excitations in the XY chain, finding an excellent agreement with numerics [28]. Note, however, that the XY chain is equivalent to a free-fermion model and thus all the single-particle modes can exactly be included in N . In contrast, for the interacting XXZ chain, restricting ourselves to the spinon modes should necessarily introduce some limitations to the quasiparticle ansatz, as demonstrated in the following subsection.

B. Local fermionic excitation
We continue with the numerical study of the excitation produced by the fermionic creation operator c † j . The fermion operators are related to the spin variables via the Jordan-Wigner transformation where σ α j are the Pauli matrices and σ ± j = σ x j ± iσ y j /2. For simplicity, we shall only consider the case where the excitation is created by c † 1 in the middle of the chain. The time-evolved state after the excitaiton is then given by where |ψ 0 is the ground state and the normalization is given by as the ground state is half filled. The time evolution is actually implemented via time-dependent DMRG (tDMRG) [43,44] in the spin-representation of the XXZ chain, by first carrying out the ground-state search and applying the string operator (19) onto the MPS representation of |ψ 0 . The calculations were performed using the ITensor C++ library [45] and a truncated weight of 10 −9 .
The results of our simulations are shown in Fig. 1 for various interaction strengths ∆. The different symbols correspond to snapshots of the entropy profile ∆S at different times, plotted against the scaled distance ζ = r/t. The quasiparticle ansatz (17) computed using the spinon fraction (18) is shown by the red solid lines. For moderate values of |∆|, one observes a very good agreement with the numerical profiles, except for a peak around ζ = 0. Note that this peak rises above the maximum value ln(2) that can be obtained from the spinon ansatz. A closer inspection for r = 0 indicates that the entropy peak also converges to a finite value for large times, with a nontrivial dependence on ∆. Moreover, one can also observe a slight broadening of the peak for larger ∆. However, the precise origin of the peak cannot be understood within our simple quasiparticle ansatz.
Systematic deviations from (17) also occur for larger ∆, especially in the attractive regime. Indeed, for ∆ = −0.5 one already observes that the edges of the profile obtained from numerics fall slightly outside of the spinon edge, whereas the bulk profile still shows a good agreement. For ∆ = −0.8 the mismatch becomes more drastic both in the bulk and around the edges, signaling the breakdown of the naive spinon ansatz. Clearly, for strong attractive interactions the local excited state should have significant overlaps with other quasiparticle excitations of the XXZ chain. In fact, as discussed in Sec. II, in this regime the maximum velocity of particle-hole excitations exceeds the spinon velocity and matches perfectly the edges of the profile, as indicated by the black dashed lines in Fig. 1. Hence, the entropy spreading should be determined by the coexistence of the spinon and particle-hole excitations, allowing to reach values beyond ln (2). Presumably, improving the ansatz (17) would require the knowledge of the overlaps with the different families of quasiparticles. Finally, it should be noted that, even though the edge locations of the profile seem to be captured, significant deviations in the bulk also occur for large repulsive interactions (see ∆ = 0.8 in Fig. 1), which might be due to bound-state contributions.

C. Local Majorana excitation
As a second example, we are going to consider local Majorana excitations, given in terms of the spin variables via and satisfying the anticommutation relations {m k , m l } = 2δ kl . Majorana operators are Hermitian and related to the fermion creation/annihilation operators as m 2j−1 = c j + c † j and m 2j = i c j − c † j . Focusing again on an excitation m 1 in the middle of the chain, the time-evolved stated is now given by The entanglement profiles ∆S obtained from tDMRG simulations of (23) are depicted in Fig. 2 for four different values of ∆. To visualize the spreading of the profile, we now plot the unscaled data against the location of the subsystem boundary. For ∆ = 0, the profile looks similar to that of the corresponding c † 1 excitation and is indeed perfectly reproduced by the quasiparticle ansatz (17). However, in the interacting case ∆ = 0, one observes a marked difference when compared to the corresponding panels in Fig. 1. Namely, the profiles in Fig. 2 clearly exceed the value ln(2), indicated by the dashed horizontal lines, which is the maximum of the ansatz (17). Nevertheless, we observe that the profiles after the m † 1 excitations can be well described by a simple rescaling of the spinon ansatz (17), as shown by the solid lines in Fig. 2. The constant factor multiplying the ansatz is chosen such that the maxima of the profiles at r = 0 are correctly reproduced. Note also that the central peak observed for the c † 1 excitation in Fig.  1 is missing for the Majorana excitation. To better understand the behaviour of the maxima, on the left of Fig. 3 we plot the time evolution of the excess entropy ∆S in the middle of the chain (r = 0) with L = 200 and for various ∆. One observes that the asymptotic value of the excess entropy grows with increasing |∆|, approaching its maximum very slowly in time. In fact, for even larger times the entropy starts to decrease again as one approaches vt ≈ L, when the fastest spinons leave the subsystem after a reflection from the chain end. This is demonstrated on the right of Fig. 3 by repeating the calculations for a smaller chain with L = 50. The emergence of a plateau is clearly visible, which then immediately repeats itself for vt > L due to the symmetry of the geometry, with the spinons reflected from the other end of the chain entering the subsystem again. However, the question why the height of the plateau depends on the interaction strength ∆ can only be answered via a more involved CFT analysis of the problem, which is presented in the next section.

IV. ENTANGLEMENT AFTER LOCAL EXCITATIONS IN CFT
The low-energy physics of the gapless XXZ chain can be captured within quantum field theory via the bosonization procedure [46]. Using the fermionic representation (2) of the chain, one introduces the Heisenberg operators c(x, τ ) = e τ H c x e −τ H , where x is the spatial coordinate along the chain and we introduced the imaginary time τ = it. Linearizing the dispersion around the Fermi points, one can approximate where ψ(x, τ ) andψ(x, τ ) are the right and left-moving components of a fermion field. The phase factors with the Fermi momentum, where k F = π/2 for a half-filled chain, are included to ensure that the chiral fermions are described by slowly varying fields. Introducing the complex coordinates w = vτ − ix andw = vτ + ix, where v denotes the Fermi velocity, they can be written in a bosonized form [46] where ϕ(w) andφ(w) are the chiral boson fields. In terms of the new bosonic variables Apart from the velocity v, the Hamiltonian (27) is characterized by the Luttinger parameter K. Both of them can be fixed from the exact Bethe ansatz solution as v = π 2 sin(γ) γ , with the usual parametrization ∆ = cos(γ). Note that v = v s (0) is just the maximum of the spinon velocity (5). In CFT language, the Luttinger liquid corresponds to a free compact boson field theory. In order to study entanglement evolution after local operator excitations, we shall thus use the framework developed for a generic CFT [15,16]. In the following we summarize the main steps of the procedure. Let us consider the state excited from the CFT vacuum |0 by insertion of the local operator O(−d), where N accounts for the normalization of the state. For the sake of generality, we consider the situation where the excitation is inserted at a distance d measured from the center of the chain. After time evolution, the density matrix reads where is a UV regularization that is required for the state to be normalizable. Working in a Heisenberg picture, the time evolution can be absorbed into the operators, and the state can be represented as where the complex coordinates of the operator insertions are given by It should be stressed that thew j coordinates are actually not the complex conjugates of w j , as we are assuming τ = it to be real, such that we can work with Euclidean spacetime. With the expression (31) at hand, one can proceed to construct the path-integral representation of the reduced density matrix, by opening a cut at τ = 0 along the spatial coordinates of the subsystem A. The Rényi entropy for integer n can then be obtained by applying the replica trick [48], i.e. sewing together n copies of the path integrals cyclically along the cuts. In turn, one can express the excess Rényi entropy ∆S n = S n (t) − S n (0) via correlation functions of the local operator as [15,16] where Σ n denotes the n-sheeted Riemann surface, with w 1 , . . . , w 2n andw 1 , . . . ,w 2n being the replica coordinates of the insertion points (32). Although the expression (34) for the excess Rényi entropy is very general, the calculation of 2n-point functions on the complicated Riemann surface Σ n may become rather involved. However, if the subsystem A is given by a single interval 0 ≤ x ≤ in an infinite chain, the geometry can be simplified by the conformal transformation which maps the n-sheeted surface onto a single Riemann sheet. This transformation leads to the holomorphic coordinates of the operator insertions while the anti-holomorphic ones are given bȳ Furthermore, if the local operators are primary fields of the CFT with respective conformal dimensions h O andh O , the 2n-point function transforms as In the end, one is left with a problem of calculating 2n-point functions on the complex plane. For the sake of simplicity, in the following we shall only consider the case n = 2, and apply the procedure outlined above to the Luttinger liquid theory, with the local excitations considered in section III.

A. Fermionic excitation
We start with the fermion creation operator, which after bosonization (25) corresponds to the field insertion where we omitted normalization factors that cancel in the expression (34). Clearly, O f (w,w) is not itself a primary operator but rather a linear combination of two. Hence, the calculation of the four-point function that appears in ∆S 2 involves a number of terms with primaries, each of which can be mapped from Σ 2 to the complex plane using the transformation rule (38). The calculation of these correlation functions can be facilitated by first performing a canonical transformation which absorbs the Luttinger parameter K in the Hamiltonian (27). However, since the variables θ and φ are actually linear combinations (26) of the chiral bosons, the change of variables corresponds to the Bogoliubov transformation where K = e 2ξ . Thus, the transformation of the Luttinger liquid Hamiltonian induces a left-right mixing of the chiral bosonic modes. In the following we shall use the shorthand notations c = cosh(ξ) and s = sinh(ξ). Clearly, our task now boils down to evaluate correlation functions of vertex operators on the complex plane with respect to the Luttinger liquid theory scaled to the free-fermion point. The n-point function of vertex operators is then well known and given by [49] n j=1 where the neutrality conditions must be satisfied, otherwise the correlator vanishes. In particular, considering the two-point function one immediately sees that the vertex operator (42) is a primary with scaling dimensions h = α 2 /2 andh = β 2 /2. With all the ingredients at hand, performing the calculation for ∆S 2 is a straightforward but cumbersome exercise, and we refer to Appendix A for the main details. It turns out that the result depends only on the cross-ratios of the holomorphic and anti-holomorhic coordinates (36) and (37), where z ij = z i − z j andz ij =z i −z j , respectively. In terms of the cross-ratios, the final result reads It is important to stress that the notation |η| should be understood as (ηη) 1/2 , since the two cross ratios are not conjugate variables. In particular, in the limit → 0 of the regularization, one has the behaviour [15,16] This yields the following limit for the Rényi entropy The result has a very simple interpretation. Namely, our excitation is an equal superposition of a left-and rightmoving fermion, and the entanglement is changed by ln(2) only when the right-moving excitation is located within the interval. In fact, this is exactly the same picture that lies behind the quasiparticle ansatz (17), without the dispersion of the wavefront. Interestingly, apart from the presence of the spinon velocity v, the limiting result (48) is independent of the anisotropy ∆. The only effect of the left-right boson mixing appears in the exponents of the cross-ratios in (46), which simply determines how the sharp step-function for ∆S 2 is rounded off for finite UV regularizations. In fact, this result is very similar to the one obtained for a non-chiral EPR-primary excitation in Ref. [16,19]. Moreover, this is also a simple generalization of the result in Ref. [25], where the superposition of purely holomorphic and anti-holomorphic primaries was considered.

B. Majorana excitation
We move on to consider the Majorana excitation The calculation of ∆S 2 follows the exact same procedure as for O f (w,w), however, one has now an even larger number of terms to consider. The main steps are again outlined in Appendix A, which lead to the result where the terms in the logarithm are given by and a new variable is introduced as .
The result is thus rather involved and cannot be written as a function of the cross-ratios alone. However, in the limit → 0, the factors in A, B, and C can trivially be evaluated using (47), as well as using Z → 1 andZ → 1. For the case ∆ = 0, this leads to the following simple result In sharp contrast, for ∆ = 0, where c = 1 and s = 0, one recovers the result (48). Hence, one arrives at the rather surprising result that the excess entropy is doubled in case of interactions, which must be a consequence of the left-right boson mixing.
Obviously, for finite values of the regularization , this transition should take place continuously, rather than giving an abrupt jump. The behaviour of ∆S 2 for = 0.1 is shown in Fig. 4 for an interval of length = 20 at a distance d = 10 from the excitation. One can clearly see the development of a plateau for times d < vt < d + , the height of which increases monotonously with ∆. Nevertheless, even for the largest value ∆ = 0.8, the expected maximum of 2 ln(2) is by far not reached. The very slow convergence towards the → 0 (or, equivalently, t → ∞) limit can be understood by looking at the structure of the terms appearing in (50). In fact, for smaller values of |∆|, the slowest converging pieces are given by η 2c 2η 2s 2 as well as (1 − η) in the expression (52) of B, due to the large-time behaviourη ≈ 1 − η ≈ ( /2vt) 2 for d vt + d. Hence, the apparent nontrivial values of the plateau in Fig. 4 is a consequence of the very slow decay ( /vt) 4s 2 , where the exponent for e.g. ∆ = 0.5 is given by 4s 2 ≈ 0.08. Clearly, observing convergence towards ∆S 2 → 2 ln(2) would require enormous time scales as well as interval lengths. Despite the different geometry considered for the CFT calculations, we expect that the result (50) should also give quantitative predictions for the finite XXZ chain in a certain regime. First of all, for the half-chain bipartition where the excitation is applied directly at the boundary, the role of the dispersion should not play an important role, as all the excitations can immediately enter the subsystem. Furthermore, one could argue that the finite chain effectively corresponds to an interval of size = L, which is the distance the quasiparticles have to cover before leaving the subsystem after reflection from the chain end. Clearly, the exact form of the plateau will not be the same in the two cases, but one expects the CFT results to be applicable in a regime vt L. Finally, there is a highly nontrivial symmetry s → −s displayed by all the terms (51)- (53) in the expression of ∆S 2 , corresponding to a change of the Luttinger parameter K → 1/K, which is expected to be observed also in the lattice calculations. Note that since K = 1 corresponds to the free-fermion point ∆ = 0, the symmetry relates interaction strengths of different sign. In Fig. 5 we show a comparison of ∆S 2 obtained from tDMRG calculations for a XXZ chain with L = 200 divided in the middle, to the CFT result (50) shown by the blue solid lines. For the latter we have set = L and d = 1 as discussed above, whereas the regularization was set by hand in order to achieve the best agreement with the numerical data. One indeed observes that the CFT result gives, up to oscillations, a good quantitative description of the XXZ numerics. Furthermore, for each ∆ = 0, we also performed the calculation for the conjugate ∆ corresponding to K = 1/K, leading to a remarkably good collapse of the curves.

V. ENTANGLEMENT DYNAMICS IN THE GAPPED PHASE
The CFT studies of the previous section give a rather good qualitative description of the entanglement spreading in the critical phase of the XXZ chain. To obtain a complete picture, in this section we shall study the dynamics in the gapped antiferromagnetic phase. For a physically motivated setting, we choose one of the symmetry-broken ground states |ψ ↑ from Eq. (8), with a nonvanishing staggered magnetization (9). We now consider local Majorana operators, defined in terms of the spin variables as Note that these operators differ from the ones in (22) discussed in the gapless phase by an interchange of the x and z spin components, but they also obey Majorana fermion statistics with anticommutation relations {m k ,m l } = 2δ kl . We focus on the case of a domain wall created bym 1 in the center of the chain, which is then time evolved by the XXZ Hamiltonian (1) Note that, in order to find the proper symmetry-broken ground state, in the DMRG simulation we add to the Hamiltonian a small staggered field in the z-direction, which is then decreased towards zero during the sweeps. First we have a look at the entropy growth ∆S for the half-chain r = 0 as a function of time, shown on the left of Fig. 6 for several values of the anisotropy ∆ > 1. One observes a clear saturation of the excess entropy for large times, which is reached very quickly for large values of ∆. The asymptotic value of ∆S decreases with ∆ and always exceeds ln (2). Remarkably, as shown on the right of Fig. 6, we find that the asymptotic excess entropy is well described by the formula ∆S = S(0) + ln(2), where S(0) is the ground-state entropy of the half-chain in the symmetry-broken state. Repeating the calculation for the excess Rényi entropy ∆S 2 , we find the exact same relation with S 2 (0). To gain a deeper understanding of the above relation, one should invoke the exact results for the reduced density matrix of the half-chain, which can can be found with the corner transfer matrix (CTM) method as [50] , where the single-particle eigenvalues are given by j = 2jφ with φ = acosh(∆), and n j = 0, 1 denotes fermionic occupation numbers. In other words, the entanglement Hamiltonian H CT M of the ground state is characterized by an equispaced single-particle entanglement spectrum. Strictly speaking, this result applies to a half-infinite chain, but in practice it holds also for finite chains of length much larger than the correlation length. Note also, that the result (58) applies for the symmetric ground state, whereas for the symmetry-broken state the term j = 0 is missing from the sum. In that case, the von Neumann and Rényi entropies can be simply expressed as [51] as well as It is easy to see that the inclusion of the term j = 0 with 0 = 0 simply yields an extra ln(2) contribution to the entropies. This change alone, however, would not explain our findings for the asymptotic excess entropy in Fig. 6, which seems to indicate that S(t) ≈ 2S(0) + ln(2) for t 1. Indeed, in order to obtain such a formula, one would have to add a double degeneracy for each j with j = 0. Let us now discuss how such a degeneracy is reflected in the eigenvalues λ l of the reduced density matrix. In fact, it is more convenient to introduce the scaled quantity where λ 0 denotes the maximal eigenvalue. For the initial symmetry-broken ground state, ν l are independent of ∆ and can only assume even integer values, with occasional multiplicities due to different integer partitions. The lowest lying λ l yield ν l = 0, 2, 4, 6, 6, . . . , i.e. the first degeneracy appears as 6 = 2 + 4. The inclusion of the 0 = 0 term simply gives an overall double degeneracy of the levels λ l . The doubling of the j for j = 0 further increases the degeneracies. Altogether, the combined effect would lead to the multiplicities (2, 4, 6) for ν l = 0, 2, 4.
To check these predictions, in Fig. 7 we have plotted the 12 lowest lying ν l calculated from the reduced density matrix eigenvalues, as obtained from tDMRG simulations after time evolving the state (57) to t = 100. One can see that the ν l lie indeed rather close to the expected even integer values, approximately reproducing the expected multiplicity structure. Interestingly, the largest deviation around ν l = 4 is found for ∆ = 5, where one actually finds the best agreement with the entropy formula, see Fig. 6. In fact, however, the contribution of these eigenvalues to the entropy is already negligible. Note that the situation for larger values of ν l is much less clear, as they correspond to very small eigenvalues λ l which are already seriously affected by tDMRG truncation errors. Although we find a nontrivial asymptotic behaviour of the half-chain entanglement, we expect that the full profile should still be described, up to a multiplicative factor, by the quasiparticle ansatz introduced in section III A, similarly to the Majorana excitation in the gapless phase in Fig. 2. Therefore, we put forward the ansatz and for the excess Rényi entropy we propose ∆S n = 1 + S n (0) ln 2 The quasiparticle fraction N must now be evaluated via (16) by using the spinon velocities (13) in the gapped phase. Note that the binary entropy functions are multiplied by a factor to reproduce our findings for the half-chain, where N = 1/2. The results of our numerical calculations for the profiles ∆S and ∆S 2 , plotted against the scaling variable ζ = r/t, are shown in Fig. 8. The solid lines show the respective ansatz (62) and (63), which give a very good description of the data for both ∆ values shown. In fact, we checked that the profiles are nicely reproduced even for ∆ = 1.5, which already corresponds to a relatively large correlation length.

A. Magnetization profiles
To conclude this section, we also investigate the spreading of the magnetization profiles for the antiferromagnetic domain wall excited bym 1 . This setting was studied previously with a focus on the edge behaviour of the profile [52]. In order to remove the dependence on the ground-state value (9) of the staggered magnetization, we consider the normalized profile which then varies between −1 ≤ M j (t) ≤ 1 along the chain. We are mainly interested in the quasiparticle description of the time-evolved profile. In fact, a very similar problem was studied for a ferromagnetic domain wall in the XY chain [28], by first expanding the excited state in the single-particle basis of the Hamiltonian, which can then be time evolved trivially.
Here we assume that the dominant weight for our simple domain wall is carried by single-spinon excitations |q . Strictly speaking, this is only possible if one considers antiperiodic or open boundary conditions on the spins, since for a periodic chain spinons are created in pairs (i.e. one actually has a pair of domain walls). The time evolved state can then be written as where ε s (q) is the spinon dispersion (10), while c(q) are the overlaps of the domain-wall excitation with the singlespinon states. Note that the momentum of a single spinon satisfies 0 ≤ q ≤ π, however, the total momentum of spinons above the quasidegenerate ground state is shifted by π. Since the domain wall is created by a strictly local fermionic operator, we assume that in the thermodynamic limit |c(q)| becomes a constant in momentum space, i.e. c(q) = e iα(q) / √ N is just a phase factor normalized by the number N of spinon states. Using this in (65), one obtains for the profile Clearly, the main difficulty of calculating (66) is due to the form factors p| σ z j |q . For the transverse Ising and XY chains, such form factors are known explicitly [53,54] and were used to obtain a double integral representation of the magnetization profile [26,28]. The hydrodynamic limit can then be obtained from the stationary-phase analysis of the integrals. Moreover, there exists a number of form factor results for the XXZ chain as well (see e.g. [55,56]), which were used in numerical studies of the magnetization profile after a spin-flip excitation [57]. Unfortunately, however, the expressions are typically rather involved or not in a representation that could be useful for our purposes. In fact, we are not aware of any results where the required single-spinon matrix elements are evaluated as a function of the spinon rapidity or momentum.
Nevertheless, based on the known results, we give a handwaving argument about how the main structure of the form factor should look like. Most importantly, we assume that it becomes singular for p → q and can be written as Here the only j-dependence is in the exponential factor that follows from the action of the translation operator, and the function F(q) denotes the slowly varying part of the form factor around its pole. The factor 1/N is required for a proper thermodynamic limit of (66). Converting the sums into integrals, one can proceed with the stationary phase analysis similarly to the XY case [28], by expanding the phases around Q = q − p = 0. Using a resolution of the pole and the definition of the step function one arrives at the following simple expression for the profile Note that the proper normalization of the profile for t = 0 requires to have π 0 dq π F(q) = 1 .
The result (69) is nothing else but the quasiparticle interpretation of the magnetization profile in the hydrodynamic limit. Indeed, the initial sharp domain wall is carried away by spinons of different momenta q and velocities v s (q), where F(q) gives the corresponding weight. Unfortunately, without an explicit analytical result on the form factor, one has to make a guess on the weight function. The simplest assumption is F(q) ≡ 1, which indeed holds true for the XY chain form factors [28]. With this simple choice one actually hasÑ = N , that is we recover the spinon fraction introduced in (16) for the description of the entropy profile. In Fig. 9 we show the comparison of this simple ansatz to the tDMRG data, with a rather good agreement for a large ∆ = 5. For ∆ = 2, however, one can already see the deviations from our simple ansatz, which fails completely for even smaller anisotropies. Thus, in sharp contrast to the case of the entanglement entropies, the spinon contributions to the magnetization cannot be taken to be equal, except for close to the Ising limit.

VI. SUMMARY AND DISCUSSION
We studied the entanglement spreading in the XXZ chain after excitations that are strictly local in terms of the fermion operators. In the gapless phase we found that the time evolution after a fermion creation operator yields an entropy profile that can be well described by a probabilistic quasiparticle ansatz for not too large ∆, assuming equal contributions from low-lying spinon excitations. On the other hand, for a local Majorana excitation we observe that the quasiparticle ansatz holds only up to a multiplicative factor, determined by the excess entropy at the operator insertion point. This is in agreement with our CFT calculations, which suggest that the excess entropy exceeds ln(2) for any ∆ = 0, with a very slow convergence towards the asymptotic value 2 ln(2). In the symmetry-broken gapped phase we considered a different Majorana excitation, creating an antiferromagnetic domain wall. For the entropy profile we find again a nontrivial prefactor, whereas our simple ansatz for the magnetization, assuming equal contributions from the spinons, holds only in the Ising limit ∆ → ∞.
The main limitation of our quasiparticle ansatz (17) is that it includes only the low-lying spinons. It is natural to ask how well such an assumption actually holds for our local excitations in the different regimes. A simple way to quantify the spectral weight of the spinons in the gapless regime is via the energy difference ∆E = ψ 0 | (m 1 Hm 1 − H) |ψ 0 of the Majorana excitation (equal to that of c † 1 by particle-hole symmetry) measured from the ground state, whereas in the gapped case we replace m 1 →m 1 . Our assumption in both regimes was that one can practically work with single-spinon states, whose energies above the ground state are given by the corresponding dispersions ε s (q) in (4) and (10), respectively. This yields the simple formula for the energy difference To test the validity of our assumption, in Fig. 10 we compare the energy difference obtained from DMRG to the formula (71) in both gapless and gapped phases. As expected, the result at the free-fermion point ∆ = 0 is exactly reproduced, while the error remains relatively small in the regime |∆| 0.5. However, not surprisingly, the overall behaviour of ∆E is not properly captured by the naive ansatz (71), especially for ∆ → −1, which is exactly what we observed for the entropy profiles in Fig. 1. On the other hand, in the gapped phase shown on the right of Fig. 10, one has a qualitatively good description in the entire regime, with the error decreasing for ∆ 1. This explains why we had a much better overall description of the entropy profiles for ∆ > 1 via the quasiparticle ansatz (62). Another feature that is not completely understood is the multiplicative factor of the spinon ansatz appearing for Majorana excitations. In the gapless phase this could be accounted for the mixing of the chiral boson modes and yields a factor 2 in the limit t → ∞ for any ∆ = 0. The exceptional behaviour of the XX chain can actually be also understood directly, using a duality transformation [58][59][60][61] that relates it to two independent and critical transverse Ising chains. Furthermore, as shown in [26], the Majorana excitation on the XX chain transforms under the dual map into a Majorana excitation acting only on a single Ising chain. Hence, the asymptotic excess entropy is given by ln (2) and there is no doubling in this case. On the other hand, in the gapped phase the prefactor in (62) seems to be nontrivially related to the ground-state entanglement entropy. Note that a similar observation was reported after a local quench in the non-critical transverse Ising chain [62], where the entanglement plateau was also found to be related to the ground-state value. A deeper understanding of these effects requires further studies.
Finally, let us comment about the case where the locality of the excitation is not imposed in the fermionic but rather in the spin picture. In other words, instead of the c † j excitation one could consider the spin operator σ + j by dropping the Jordan-Wigner string in (19). According to our tDMRG calculations carried out for this case, the entropy profiles change completely, becoming more flat in the center with a maximum that stays way below ln (2). In short, the fermionic nature of the local excitations turns out to be essential for the applicability of the quasiparticle description.
They are composed of chiral fermion fields which, after the Bogoliubov transformation (41), can be written as vertex operators (42) involving chiral boson fields. The holomorphic and anti-holomorhic components of the vertex operators are summarized in the table below, where c = cosh(ξ) and s = sinh(ξ).
ψ ψ †ψψ † α −c c s −s β −s s c −c We start by evaluating the two point function in the denominator of (34). Using the fact that vertex operators are primaries with conformal dimensions h = α 2 /2 andh = β 2 /2, one immediately obtains the nonvanishing two-point functions on the plane as From (32) we have w 1 − w 2 =w 1 −w 2 = 2 , thus we obtain for the two-point functions Let us now move to the four-point function on the Riemann surface Σ 2 . This is a sum of many terms, from which the nonvanishing contributions allowed by the neutrality conditions (44) are given by We first analyze the Jacobian of the transformation (38) from Σ 2 → Σ 1 . The derivatives of the mapping are given by Introducing the variable one obtains for the first four contributions in (A5) −2(c 2 +s 2 ) χ c 2 /2χs 2 /2 χ s 2 /2χc 2 /2 = −2(c 2 +s 2 ) |χ| c 2 +s 2 , whereas for the last two contributions we have, respectively −2(c 2 +s 2 ) χ c 2χ s 2 , −2(c 2 +s 2 ) χ s 2χ c 2 . (A9) In order to obtain an expression in terms of the cross-ratios, one can rewrite (A7) as Putting everything together, one arrives at the four-point function 2 (2 ) −2(c 2 +s 2 ) |η| (c+s) 2 + |1 − η| (c+s) 2 + 1 .
Evaluating the four-point function for the Majorana excitation (A2) is more cumbersome, since one has a large number of terms to consider. There are, however, some simple rules and symmetry arguments which make the task easier. First of all, one should clearly always have the same number of creation and annihilation operators, for the neutrality conditions (44) of the vertex correlation functions to be satisfied. This already drastically reduces the number of terms to consider. The remaining ones can be collected into families, some of them given by (A5).
Let us consider the family generated by the first term in (A5), by allowing permutations of the left-and right-moving operators separately (i.e. interchanging the first or last two operators). If only the first or last two are interchanged, the vertex correlator (A10) is modified by replacing whereas the correlator remains the same if both of them are interchanged. The next family is generated by the second term in (A5), which is actually related to the first one by Hermitian conjugation. Hence this just gives a factor of two. The same argument holds for the next two families, where interchanging only one pair modifies the correlator in (A11) as |η| 2cs → |η| −2cs .
Finally, the single interchange in the fifth family leads to whereas the last family follows by interchanging c and s above. There are, however, two additional families appearing where the left-and right-moving particles are intertwined. They are given by the representative correlators ψψ † ψ †ψ , ψ ψ †ψ † ψ .
Defining the variable the corresponding Jacobians contain the factors σ c 2σ s 2 and σ s 2σ c 2 , respectively. Furthermore, the vertex correlation functions yield and each term comes with a double multiplicity. Collecting all the terms, the four-point function takes the form 2 (2 ) −2(c 2 +s 2 ) (2A + B + C) , where the factors A, B and C are reported in (51)- (53).