Entanglement spreading after local and extended excitations in a free-fermion chain

We study the time evolution of entanglement created by local or extended excitations upon the ground state of a free-fermion chain. A single particle or hole excitation produces a single bit of excess entropy for large times and subsystem lengths. In case of a double hole, some of the coherence between the excitations is preserved and the excess entropy becomes additive only for large hole separations. In contrast, the coherence is always lost for particle-hole excitations. Multiple hole excitations on a completely filled chain are also investigated. We find that for an extended contiguous hole the excess entropy scales logarithmically with the size, whereas the increase is linear for finite separations between the holes.


Introduction
Entanglement properties in one-dimensional quantum systems have been the topic of intensive research in the last decades [1,2,3]. The most emblematic feature uncovered is the existence of an area law, which bounds the scaling of entanglement entropy in ground states of local Hamiltonians [4]. Notable exceptions include quantum chains at criticality [5], where the logarithmic violation of the area law shows a universal character that can be understood via the underlying conformal field theory (CFT) [6]. Such a moderate growth of ground-state entanglement is the key ingredient behind the classical simulability of quantum chains via matrix product state methods [7].
The area law with eventual logarithmic corrections is, however, not restricted to the ground-state scenario. Indeed, under some reasonable assumptions, it has been shown to extend to low-lying excited states of generic local Hamiltonians [8]. The simplest case is the one of translational invariant systems possessing well defined quasiparticle excitations. The entanglement of such few-particle excitations have been studied in various integrable chains [9,10,11,12,13], in the context of quantum field theory [14,15,16,17], and recently also in some nonintegrable models [18,19]. In the context of free field theory, few-particle excitations carry a finite amount of excess entanglement, given by a simple binomial expression via the probabilities of finding the quasiparticles in the subsystem [14,15]. In particular, the quasiparticle contributions are equal and independent of the momenta. On the lattice, however, coherence effects between quasiparticles are retained, leading to corrections that are most dominant for nearby momenta [12,13].
Another class of low-lying excitations that was studied intensively is the one generated by primary operators in CFT [20,21]. In the simplest case, these correspond to low-lying excited eigenstates of the underlying critical quantum chain. More generally, one could consider excitations created by primaries that are local in space, and thus do not yield an eigenstate. The time evolution of entanglement after such local operator excitations was first studied in [22,23,24], finding a light-cone propagation and a finite excess entanglement given via the quantum dimension of the primary operator. The CFT approach was also extended to the case of multiple excitations [25,26,27], finding an additive behaviour of the excess entropy.
The lattice counterpart of local operator excitations has first been studied for a critical transverse Ising chain [28], confirming the CFT predictions. Entanglement evolution after local fermionic excitations, creating a domain-wall initial state in the ordered phase of the Ising and XY chains were studied in [29,30,31,32]. The resulting entanglement profiles along the chain could be understood by a simple quasiparticle picture. Namely, the excess entropy is determined via the density fraction of the quasiparticles that arrive to the entanglement cut by a binary expression [32], analogously to the case of quasiparticle eigenstates [14]. The quasiparticle description of the excess entropy has also been generalized to the interacting XXZ chain [33].
Here we address the question how the excess entropy changes when considering multiple local excitations in a free-fermion chain. More precisely, we consider excitations that can be written as a product of local fermion creation and annihilation operators and study the ensuing entanglement evolution via Gaussian techniques [34]. Starting with the case of a single particle or hole excitation, we derive an expression that puts the heuristic quasiparticle picture obtained in [32] on a firm ground. The situation is more complicated for a double hole excitation, as coherence effects enter the picture and the excess entropy becomes additive only for very large hole separations. For finite separations we obtain an asymptotic formula for large times. Interestingly, we find that particle-hole excitations always lose their coherence and contribute additively for large times. In case of multiple hole excitations on a completely filled background, we find that the coherence is essentially preserved only for an extended contiguous hole, which leads to a logarithmic scaling of the asymptotic excess entropy with the size of the excitation. In contrast, for a finite separation between the holes one obtains a linear scaling, as the individual hole excitations become effectively distinguishable.
The manuscript is structured as follows. The model and the general setting is introduced in Sec. 2, along with the main tools required for subsequent calculations. The following sections are devoted to the study of entanglement after single (Sec. 3) and double hole (Sec. 4), as well as particle-hole excitations (Sec. 5). Finally, the case of multiple hole excitations is discussed in Sec. 6. Our conclusions are presented in Sec. 7, followed by three Appendices on various details of the calculations. Entanglement spreading after local and extended excitations in a free-fermion chain 3

Model and setting
We consider an infinite hopping chain given by the Hamiltonian with fermionic creation/annihilation operators satisfying canonical anticommutation relations {c † m , c n } = δ mn . The Hamiltonian is diagonalized by a Fourier transform, leading to the dispersion ω(q) = − cos q − µ, where the chemical potential µ sets the filling of the chain via the Fermi wavenumber, satisfying ω(±q F ) = 0. The ground state |ψ 0 of the chain is characterized by its correlation matrix where the filling is given by C mm = q F /π, and all the higher-order correlation functions factorize according to Wick's theorem.
In the following sections we will consider excitations created upon the ground state of the chain by acting with an operator O, supported on a single or a few lattice sites, and then letting the system evolve, governed by the Hamiltonian in (1). The time-evolved state can be written as where the factor N O ensures normalization. We will focus on a class of excitations that preserve the Gaussianity of the state, i.e. Wick's theorem also applies to O |ψ 0 , see Appendix A for details. Since the time evolution is governed by a free-fermion Hamiltonian, |ψ(t) is just another Gaussian state characterized by the time-dependent correlation matrix In order to evaluate these matrix elements, one only needs to determine the effect of the local excitation. Denoting with C ′ the correlation matrix immediately after the excitation, the time evolution is obtained via the unitary transformation where the matrix elements of U are given via Bessel functions as Our goal is to calculate the entanglement entropy of the segment A = [1, L] which, owing to the Gaussianity of the problem, can be obtained via standard correlationmatrix techniques [34]. Indeed, the entropy can be calculated as Entanglement spreading after local and extended excitations in a free-fermion chain 4 from the eigenvalues ζ k (t) of the reduced correlation matrix C A (t), with matrix elements (4) restricted to the interval m, n ∈ A. In fact, we are rather interested in the excess entropy where S 0 denotes the ground-state entropy. One should stress that we are subtracting the entropy before the excitation takes place, thus ∆S(0) = 0 and may become negative. The quantity ∆S(t) provides information about the excess entanglement generated by the dynamics with respect to the ground state, and is studied for a range of single or few-particle excitations in the following sections.

Single hole or particle excitation
We start with the simplest case of a strictly local excitation, a single hole created by an annihilation operator inserted at site ℓ as where the normalization factor is given by the density N = q F /π. As pointed out in the previous section, the only ingredients we need are the correlation matrix elements C ′ mn = ψ 1h | c † m c n |ψ 1h after the excitation. Applying Wick's theorem one has The matrix C ′ has a simple form and is determined via the ground-state correlations. In particular, one has C ′ mℓ = C ′ ℓn = 0, i.e. the annihilation operator erases all the correlations along the ℓ-th row and column of the matrix. Moreover, the change in the density C ′ mm is non-local, decaying as a power-law away from the insertion site m = ℓ, even though the excitation is created by a local operator. However, the overall particle number decreases by one, Tr(C ′ − C) = −1, as it should for a single hole excitation, which follows from the fact that C 2 = C for a pure state.
We shall now apply the time evolution on C ′ as given by (5). As the first term in (10) is simply the ground-state correlation matrix, the time evolution affects only the difference term and one has where the elements of the difference matrix can be cast as Using the integral representations of the matrix elements of U and C, it is easy to show that the propagator is given by Entanglement spreading after local and extended excitations in a free-fermion chain 5 Note that the propagator depends only on the distance n−ℓ measured from the location of the excitation. It differs from the time-evolution operator in (6) only in the range of integration, which is now restricted to the Fermi sea. This has the clear physical meaning that a hole can only be excited from the initially filled modes. In order to evaluate the entropy, one needs to restrict the correlation matrix (11) to the interval A. From the structure of the matrix elements (12) one can immediately see, that ∆C A (t) is a rank-one matrix with a single nonzero eigenvalue In fact, ∆ζ(t) is nothing else but the decrease in the density within the interval at time t due to the hole excitation. However, the problem is that the difference matrix does not commute with the ground-state reduced correlation matrix, [C A , ∆C A (t)] = 0, thus one cannot immediately read off the spectrum of C A (t). The eigenvalues ζ k (t) in the middle of the spectrum are shown by the symbols in Fig. 1, for a hole excited at ℓ = 0 next to the left boundary of an interval with L = 50 in a half-filled chain. One observes that all but one of the eigenvalues move rapidly towards their ground-state values, indicated by the horizontal lines. The single anomalous eigenvalue is well captured by 1 − ∆ζ(t), shown by the red line, up to oscillations which tend to be stronger in the data. Thus, despite the non-commuting property, the rank-one update ∆C A (t) of C A modifies essentially only one of its eigenvalues lying exponentially close to one, accounting for the decrease of density in the subsystem. One should remark that the upwards motion of the anomalous eigenvalue is accompanied by avoided crossings in the spectrum, whenever one of the ground-state eigenvalues is reached. The first avoided crossing in Fig. 1 can be seen at t ≈ 60 (between light blue and dark blue symbols) while a second one just starts at t ≈ 200 (dark blue and yellow), and the same feature continues towards the upper edge of the spectrum.
The above approximate description of the spectra can now be used to understand the behaviour of the excess entropy, shown on the left of Fig. 2 for ℓ = 0 and two different sizes of the interval at half filling. Indeed, the change of entropy can mainly be accounted for the anomalous eigenvalue, and using (7) one has ∆S(t) ≈ s[∆ζ(t)], where we used the symmetry s(x) = s(1 − x). This approximation is shown by the lines which feature a plateau region for 0 < t < L, slightly overestimating the numerical data, followed by a very slow decay for t > L where the approximation becomes increasingly more accurate. Indeed, this is due to the behaviour of the anomalous eigenvalue in Fig. 1, which first tends towards the middle of the spectrum, giving a large contribution to the entropy, while for t > L moves slowly upwards, yielding the decay of the entropy. Similar features can also be observed for a hole excitation with ℓ < 0, created at some distance from the interval, as shown in the right of Fig. 2. The main difference is that the plateau region is shifted towards |ℓ| < t < |ℓ| + L, since the excitation must first reach the interval, and the plateau is also rounded off.  For the quantitative understanding of the excess entropy, one has to analyze the expression in (14), which is the density fraction of the propagating hole within the interval. From simple symmetry arguments one finds that for |ℓ| ≪ t ≪ |ℓ| + L one has ∆ζ(t) → 1/2, since half of the density travels to the right from the initial hole location. This immediately yields that the height of the entropy plateau is given by ln 2. Furthermore, one can even take into account the dispersion of the fermionic modes in a semi-classical picture. This is equivalent to calculating the stationary-phase approximation of (14) using the integral representation of the propagator (13). As shown in Appendix B, for t, L ≫ 1 and ℓ ≤ 0, this leads to the following expression where v q = dω dq = sin q are the single-particle velocities and Θ(x) is the Heaviside step function. This has the obvious semi-classical interpretation that the change of density is given by the fraction of the propagating modes located within the interval. Note that the 1/2 distance correction in the arguments of the step functions is due to the fact, that the initial hole density is not sharply localized at site ℓ. This correction only becomes important for small |ℓ| and at short times.  To check the validity of (15), in Fig. 3 we have calculated ∆ζ(t) for ℓ = 0 and various fillings, finding a good agreement with the semi-classical formula up to oscillations. One also observes a change of behaviour below and above half filling. Indeed, for q F < π/2 one has v max = v q F , and ∆ζ(t) starts to decay at v q F t ≈ L. In contrast, for q F > π/2 one has v max = 1, and the decay starts already at t ≈ L. However, there is a clear change of behaviour also at v q F t ≈ L, where the modes π/2 < q < q F having a pair π − q with the same velocity, are depleted. This changes the slope of decay and is observed as a cusp in the semi-classical approximation for q F = 0.8π in Fig. 3. The numerical data, however, interpolates smoothly between these two regimes. It is also interesting to have a look at the asymptotic decay. For v q F t ≫ L one has Entanglement spreading after local and extended excitations in a free-fermion chain 8 which yields for the excess entropy Thus ∆S(t) receives a logarithmic correction and its decay is slower than algebraic. Finally, it should be stressed that the above described relation between entropy and excess density improves for increasing filling. Indeed, in the case of complete filling q F = π, the correlation matrix becomes C ij = δ ij , and thus [C A , ∆C A (t)] = 0. Hence the formula ∆S(t) = s[∆ζ(t)] for the excess entropy becomes exact in this case.
To conclude this section, we also discuss the single particle excitation given by The correlation matrix immediately after the excitation reads whereC mn = δ mn − C mn is the ground-state correlation matrix after a particle-hole transformation. After time evolution one then has where Here the propagator is given by the integral whereF = [q F , π] ∪ [−π, −q F ] is the complement of the Fermi sea F , reflecting the fact that a single particle excitation can only create modes that were initially empty. It is easy to see, that the particle excitation with filling q F /π is trivially related to the hole excitation with filling 1 − q F /π via a particle-hole transformation.

Double hole excitation
Next we consider a double hole excitation given by where ℓ 1 = ℓ 2 . The correlation matrix immediately after the excitation can again be evaluated by Wick's theorem and reads Entanglement spreading after local and extended excitations in a free-fermion chain 9 After time evolution one has again the form (11) with the difference matrix given by where we have defined a = N 2 /N 2 and b = −N C ℓ 1 ℓ 2 /N 2 . The structure of (25) suggests that we have now to deal with a rank-two update of the correlation matrix. However, one cannot immediately write down the eigenvalues of ∆C(t) since the basis states I n−ℓα (t) with α = 1, 2 are not orthonormal. Nevertheless, introducing the notation as well as λ ± = (λ 1 ± λ 2 )/2, it is shown in Appendix C that the two nonzero eigenvalues of ∆C(t) can be obtained as Note that we have suppressed the argument t for notational simplicity. It should be stressed that the eigenvalues now depend, apart from the densities λ 1,2 , also on the complex overlap ∆ in (26). In the special case C ℓ 1 ℓ 2 = 0 one has a = 1 and b = 0 and thus (27) simplifies to Analogously to the case of a single hole excitation, we now expect that two anomalous eigenvalues should appear around 1 − ∆ζ 1,2 in the spectrum ζ k (t). This is illustrated in Fig. 4 for two different separations between the holes. On the left we have ℓ 1 = 0 and ℓ 2 = −1, such that the holes are located at neighbouring sites next to the boundary of the interval. For t < L one observes that there is indeed an eigenvalue well described, up to oscillations, by the smaller of the two 1 − ∆ζ 1,2 , whereas the larger one lies between two ζ k (t) that tend to converge toward it from both sides. The decay regime t > L is dominated by the smaller anomalous eigenvalue, which moves slowly upwards through the spectrum. On the right of Fig. 4 we also show the case ℓ 1 = 0 and ℓ 2 = −5, where both 1 − ∆ζ 1,2 seem to give a better description of the anomalous eigenvalues, featuring similar avoided crossings as in the single-hole case in Fig. 1.
The behaviour of the excess entropy is shown in Fig. 5 for ℓ 1 = 0 and various ℓ 2 . One expects that the entropy is now approximated by which is shown by the solid lines for L = 50 on the left of Fig. 5. While the approximation seems to improve for increasing separations between the holes, (29) always overestimates the excess entropy and fails especially for short times t ≪ L. This is most prominent in case of ℓ 2 = −1, which could already be anticipated from the behaviour of ∆ζ 1,2 in Fig. 5. Nevertheless, one also observes that the ansatz (29) improves for larger times, and it correctly gives the height of the plateau in the asymptotic regime |ℓ 1 |, |ℓ 2 | ≪ t ≪ L. This is illustrated on the right of Fig. 5, where the dashed lines indicate the ansatz (29) with the ∆ζ 1,2 evaluated via stationary-phase analysis. Indeed, for the densities one obtains λ 1,2 ≃ 1/2 as for the single hole, whereas the overlap in (26) is evaluated in Appendix B and yields Substituting into (27), one obtains for the asymptotic eigenvalues Entanglement spreading after local and extended excitations in a free-fermion chain 11 Interestingly, the plateau height is not monotonously increasing with the distance |ℓ 1 − ℓ 2 |. In particular, at half filling one has Im(∆) = 0 for mod (|ℓ 1 − ℓ 2 |, 4) = 0 and thus one obtains the maximal excess entropy ∆S = 2 ln 2. In general, ∆S also converges towards this maximum at large separations |ℓ 1 − ℓ 2 | ≫ 1, since from (30) one has ∆ → 0 as the coherence between the two holes is lost and they become effectively independent.

Particle-hole excitations
Our next example is the particle-hole excitation where, in order to ensure orthogonality ψ 0 |ψ ph = 0, we subtracted the ground-state expectation value from the operator. The particle-hole excitation defined this way still leads to a Gaussian state, as verified in Appendix A. The correlation matrix after the excitation then reads Thus the contributions from a single hole (10) as well as from a single particle (19) become additive, entering the expression with different signs. After time evolution this yields Following the arguments applied for a single excitation, one expects two anomalous eigenvalues to appear around 1 − λ 1 andλ 2 , respectively, where Consequently, the excess entropy should also be approximately additive The above ansatz is checked against numerical results in Fig. 6 for ℓ 1 = ℓ 2 = 0 at half filling. On the left one can clearly see that the anomalous eigenvalues are well captured in the ζ k (t) spectrum. The excess entropies are shown on the right of Fig. 6 for two interval sizes, showing convergence towards the predicted plateau with ∆S = 2 ln 2. The results are similar for ℓ 2 < ℓ 1 < 0, with the emergence of an intermediate plateau around ∆S = ln 2 in the regime |ℓ 1 | ≪ t ≪ |ℓ 2 |, in case of larger separations. However, in contrast to the double hole excitation, the asymptotic excess entropy for |ℓ 2 | ≪ t ≪ L is always given by ∆S = 2 ln 2.
In order to understand this, let us first remark that the two difference terms in (34) actually do not commute exactly. In fact, one could proceed similarly to the double hole case, treating the rank-2 update matrix as in (25). The eigenvalue problem could be dealt with in a very similar fashion, by introducing the overlap between the particle and hole propagator However, this overlap decays very quickly towards zero, which can be understood from a stationary-phase argument. Indeed, using the definitions (13) and (22), the phase of the resulting double integral could be made stationary by choosing equal momenta, see Appendix B. However, since the domains of the two integrals are now complementary in momentum space, the stationary phase condition cannot be satisfied. This explains why the coherence is not preserved for large times in case of a particle-hole excitation, which thus effectively behaves as two independent local excitations. Finally, one should remark that the plateau value of the excess entropy for ℓ 1 = ℓ 2 = ℓ could also be obtained via the CFT approach to local operator excitations [22,23,24]. Indeed, using the standard bosonization procedure [35], the operator applied to the ground state in (32) can be rewritten, up to an overall normalization factor, in terms of the chiral boson fields φ andφ as Here the first two terms are the chiral (right-and left-moving) currents, corresponding to primary fields with conformal dimensions (h,h) = (1, 0) and (0, 1), respectively. The last two terms correspond to Umklapp processes, removing a fermion at the right edge of the Fermi sea and placing it back at the left edge or vice-versa, and they both have conformal dimensions (1/2, 1/2). Thus the operator in (38) is a sum of chiral and non-chiral primaries with the same scaling dimension h +h = 1. This is exactly the situation considered in [36], finding the result that the excess entropy is simply given by the Shannon entropy calculated from the weights of the various operators. However, since (38) is an equally weighted sum, one recovers immediately the result ∆S = 2 ln 2.

Multiple hole excitations
In the last part we shall consider the case of multiple hole excitations. For the sake of simplicity, we restrict our attention to the case of a completely filled chain, q F = π, and the excited state reads where X is an index set labeling the empty sites. Indeed, for a such a simple initial state the difference matrix can be found explicitly as where the propagator I n−ℓ (t) = i n−ℓ J n−ℓ (t) is given via the Bessel function. Furthermore, due to C mn = δ mn , the spectrum of the reduced correlation matrix is given exactly by ζ k (t) = 1 − ∆ζ k (t), where ∆ζ k (t) are the eigenvalues of (40) reduced to the subsystem, m, n ∈ A. Note also that S 0 = 0 and thus ∆S(t) = S(t).
In order to understand the entanglement evolution from such a multiple hole excitation, we shall apply a duality argument. First of all note, that thanks to the property the propagators corresponding to different lattice sites form a complete orthonormal basis on the whole chain. The difference matrix in (40) is built by choosing a subset ℓ ∈ X of this basis, and then reducing the matrix to the subsystem m, n ∈ A. Now the duality amounts to interchanging the roles of the basis and spatial indices, by defining the so-called overlap matrix where the indices are restricted to ℓ, ℓ ′ ∈ X. It is then easy to show, that the nonzero eigenvalues of the reduced matrices ∆C A (t) and M X (t) are identical. In fact, this is completely analogous to the duality between the reduced correlation matrix of a Fermi sea ground state and the overlap matrix indexed by the occupied momenta [37,38,39]. The duality relation is particularly useful if the number of holes is small compared to the size of the subsystem. Indeed, the result (14) for the single hole, as well as (28) for the double hole are immediately reproduced. Furthermore, one could consider directly the limit of a half-infinite chain, in which case the expression further simplifies to where is the so-called discrete Bessel kernel. This kernel is known to describe the time-evolved correlation matrix for a step-like initial state (i.e. the limit K → ∞), also known as the domain-wall quench [40,41,42,43,44,45].
We first consider the case where the excitation is given by a contiguous set of K holes located just next to the interval, X = [−K + 1, 0]. Then, by relation (43), we find that the half-chain entanglement for an extended excitation is exactly equal to that of a finite interval [1, K] in the domain-wall quench. In other words, the size of the subsystem and that of the excitation are dual to each other. In particular, the case K → ∞ corresponding to the domain wall is self-dual. In Fig. 7 we show the results for the entropy for various K. One can observe that up to t ≈ K, the entropy follows the domain-wall quench result which, up to oscillations, shows a logarithmic growth with c 0 ≈ 0.4785 [46,47,48]. This is indicated by the red dashed line. Around t ≈ K, where the last part of the extended hole propagates into the right half-chain, the entropy shows an abrupt jump. From here on, one can observe the convergence to a finite and K-dependent value. Using the asymptotics of the discrete Bessel kernel with m, n fixed one has [49] lim t→∞ B mn (t) = sin π 2 (m − n) π(m − n) , which is nothing else but the ground-state correlation matrix at half filling. Hence, we can conclude that the asymptotic entropy assumes the ground-state value for an interval of length K, shown by the horizontal dotted lines in Fig. 7.
The interpretation of the result is the following. If the K holes are located at neighbouring sites, they effectively behave as a composite excitation. Similarly to the case of a single hole, half of the density propagates to the right and the other half to the left. However, the coherence between parts of the composite excitation is preserved and, despite the density being spread out over large distances, the inhomogeneous halfinfinite chain effectively behaves as a finite homogeneous one of size K. The asymptotic entanglement thus scales as 1/3 ln K, i.e. the contribution of the holes is not additive. In the next example, we demonstrate what happens if the holes are located on nonadjacent sites, i.e. they become distinguishable. For simplicity, we shall check only the case where K holes are placed periodically at sites 0, −p, −2p, . . . , −(K − 1)p. In Fig. 8 we show the cases p = 2 (left) and p = 3 (right), respectively. One can immediately see that the situation is now completely different, with the increase being linear and the asymptotic entropy proportional to K. Indeed, due to the distinguishability of the excitations, they start to contribute independently to the entropy, and one can recognize a kind of step structure in the entropy increase. Regarding the asymptotic value, this is still given by the discrete sine kernel (46), however with a subsystem that includes only every p-th site. Such a non-contiguous subsystem is well known to have an extensive entropy [50,51]. In particular, the case p = 2 corresponds to the situation where the asymptotic overlap matrix becomes diagonal, due to the special checkerboard structure of the matrix elements in (46), and thus one has S(∞) = K ln 2. In the case p = 3, the subsystem involves odd distances, thus preserving some of the correlations between the holes which leads to S(∞) < K ln 2. The horizontal dotted lines indicate the actual values of the asymptotic entropy, and one can see a clear convergence towards them.
Finally one can check what happens when considering a finite subsystem, in which case we have to use the form (42) of the overlap matrix. The results for L = 50 are shown in Fig. 9, for both a contiguous set of holes (left) as well as for p = 2 (right).
In the latter case the entropy starts to decrease immediately after t ≈ L, when the first hole leaves the subsystem, and the decay continues in a monotonous fashion as the other holes follow. In contrast, for the extended excitation one first observes a jump at t ≈ L, followed by a second smaller jump at t ≈ L + K, when the extended excitation leaves the interval. The slow decay of the entropy sets in only after this transient, the quantitative description of which would require some further investigations.

Conclusions
We have studied the spreading of entanglement for a family of excitations that can be written as a product of a few local fermion operators in a hopping chain. In case of a single particle or hole, the excess entanglement is essentially determined by the density fraction of the excitation located within the subsystem. This can be well approximated by a semi-classical formula involving the single-particle velocities of the fermionic modes. The result is thus completely analogous to the ones found for local fermionic excitations in the XY [31,32] as well as the XXZ chain [33]. The situation is more complicated for a double hole excitation, where some of the coherence between the holes is preserved by the time evolution. The two holes thus do not behave independently and the additivity of the excess entropy is recovered only in the limit of large separations. In fact, this behaviour is completely analogous to the case of two-particle eigenstates, where the additivity is recovered only for large differences between the momenta [12,13].
In sharp contrast, the coherence of particle-hole excitations decays rapidly in time and thus the entropy becomes additive for arbitrary separations. It should be stressed, however, that it is crucial to consider a particle-hole excitation that is orthogonal to the ground state. If the vacuum expectation value in (32) is not subtracted, then offdiagonal contributions appear in the difference matrix, similarly to the case of the double hole in (25). In the particle-hole case, however, we observed that the eigenvalues ∆ζ 1,2 of the difference matrix cannot be directly connected to the anomalous eigenvalues in the ζ k (t) spectrum. This is probably due to the nonvanishing overlap between the ground and excited states. In any case, we found that the asymptotic excess entropy is given by a nontrivial value that is strictly less then 2 ln 2, converging towards it only for large separations. Understanding this mechanism requires some further investigations.
The case of multiple hole excitations was studied only at complete filling, and could be related to the domain-wall problem by a duality argument. A contiguous set of holes was found to behave as a composite excitation, with the asymptotic entropy scaling logarithmically with the number of holes. On the other hand, a finite spacing between the holes makes them distinguishable, and leads to a linear scaling of the entropy. In fact, this latter scenario is more reminiscent to the one observed in a global or inhomogeneous quench [52,53,54,55]. It would be interesting to check how these results for multiple excitations generalize to arbitrary fillings.
The characteristic feature of a local excitation is the appearance of an anomalous eigenvalue in the reduced correlation matrix spectrum, which is responsible for the excess entropy. One should note that this behaviour is very reminiscent of the one found earlier in the local quench problem, for a subsystem chosen symmetrically around the initial cut [56]. Indeed, this mechanism was found to lead to the slow asymptotic decay (17) of the entropy in both cases. However, while for a local excitation the anomalous eigenvalue is related to the change of density in the subsystem, in a local quench the density is unchanged and thus the origin of the mechanism must be more subtle. On the other hand, the half-chain entropy in the local quench shows a logarithmic increase [57,58], in contrast to a finite excess entropy for the local operator excitation. The proper understanding of the similarities and differences between these two nonequilibrium protocols is an interesting open question for future research.
We show that the state O |ψ 0 is Gaussian when In other words, acting with an arbitrary linear combination of Majorana operators on the ground state leads to a Gaussian excited state. We shall first consider the case where the coefficients w m are real, thus the operator is unitary, O † O = 1. In this case the excited state is clearly Gaussian. Indeed, working in the Heisenberg picture, one obtains a linear combination of the original operators An arbitrary p-point correlation can then be rewritten as

np
A m 1 n 1 A m 2 n 2 . . . A mpnp a n 1 a n 2 . . . a np (A.5) and Gaussianity follows from the multilinearity of Wick's theorem.
The general case will be handled as follows. For an arbitrary operator O, we shall construct a unitary U such that O |ψ 0 = U |ψ 0 is satisfied, and thus the Heisenberg operators can be defined via a ′ m = U † a m U. First note, that for an arbitrary quadratic Hamiltonian there exists an orthogonal transformation that brings the skew-symmetric Hamiltonian matrix in the canonical block-diagonal form Rewriting the Majorana fermions as one arrives at the diagonal form Importantly, since the ǫ k are chosen to be all positive, the ground state of H is the vacuum of the b k operators.
Let us now consider an operator as in (A.3) with general complex weights w m . Its action on the ground state can be written as where we used b k |ψ 0 = 0 and defined the real weights (A.11) Now we can replace again the operators b † k by either d 2k−1 or id 2k , respectively, which allows us to rewrite which is now an excitation with purely real coefficients. Hence, transforming back to the original Majorana variables, we obtain a unitary operator defined as Let us now consider an operator of product form O = O 2 O 1 , where the procedure described above can be iterated. Indeed, one can rewrite where the Heisenberg operators are defined as and U 1 is the unitary associated to O 1 . In the next step the Heisenberg operator O ′ 2 is written as a linear combination of the a ′ m operators, and one follows the exact same procedure to find the unitary operator U ′ 2 associated to it. Clearly, the procedure can be iterated for an arbitrary number of such operators in a product form.

Appendix B. Stationary phase calculation
Here we derive the stationary phase approximation of (14). We need to evaluate the following oscillatory double integral ∆ζ(t) = π q F L n=1 q F −q F dp 2π It is easy to see, that the stationary phase condition yields q = p, one can thus introduce the variables Q = q − p and P = (q + p)/2 and expand the phases around Q = 0. Using and setting v P = sin(P ), one arrives at The only difference is the phase factor e −iqℓ 1 +ipℓ 2 = e iP (ℓ 2 −ℓ 1 ) e −iQ(ℓ 1 +ℓ 2 )/2 , (B.7) which then yields the semi-classical result