Ballistic propagation of a local impact in the one-dimensional $XY$ model

We investigate propagation of a local impact in the one-dimensional $XY$ model with the anisotropy $\gamma$ in a magnetic field $h$. Applying a local and instantaneous unitary operation to the ground state, we numerically observe a light-cone-like propagation of the magnetization in the $z$ direction. In particular, we focus on the parameter region $0\leq\gamma\leq1$ and $0\leq h \leq2$ of the model. The dynamics in the light-cone region varies depending on the parameters. By combining numerical calculation with an asymptotic analysis, we find the following: (i) for $h\geq1-\gamma^{2}$ except for the case on the line $h=1$ with $0<\gamma<\sqrt{3}/2$, a wave front propagates with the maximum group velocity of quasiparticles, except for the case $\gamma=1$ and $0<h<1$, in which there is no clear wave front; (ii) for $h<1-\gamma ^{2}$ as well as on the line $h=1$ with $0<\gamma<\sqrt{3}/2$, a second wave front appears owing to multiple inflection points of the dispersion relation of quasiparticles. The propagation velocity of this second wave front is given by the amplitude of the group velocity at the second local extrema; (iii) for $h=1-\gamma^{2}$, the velocity of the second wave front vanishes, and as a result, the magnetization change forms a ridge at the impacted site. We find that the height of the wave front decays in a power law in time $t$ with various exponents depending on the model parameters, by an asymptotic analysis. We find a power-law decay $t^{-2/3}$ on the wave front except for the line $h=1$, on which the decay can be given by either $\sim t^{-3/5}$ or $\sim t^{-1}$. For $h=1-\gamma^{2}$, the magnetization at the impacted site relaxes in a power law $t^{-1/2}$ in time $t$ as opposed to the relaxation according to the power law $t^{-1}$ in other cases.


I. INTRODUCTION
Non-equilibrium dynamics of quantum many-body systems have been of great interest theoretically and experimentally. Recent experimental and numerical advance in simulating and examining quantum dynamics has motivated a wide range of studies on dynamics of isolated quantum systems [1,2]. The problems of thermalization and information propagation in isolated quantum systems are fundamental issues in this field. Important questions include how and under what conditions a pure initial state approaches to thermal equilibrium under unitary time evolution. Intensive studies in the last two decades have made remarkable progress in understanding the condition and mechanism of the thermalization in isolated quantum systems [3][4][5][6][7]. A large number of theoretical and experimental studies, including the early investigation by von Neumann [8], have shown that local observables generally relax to their steady values, which in most cases are described by a thermal ensemble [9][10][11][12][13][14][15]. On the other hand, understanding of non-equilibrium dynamics towards the seemingly thermal state has not been well established yet since the way of equilibration varies considerably among the systems, and even the generic equilibration timescale has not been well understood [16][17][18].
Among the studies of non-equilibrium dynamics in isolated quantum systems, ballistic spreading of a signal, namely the light-cone dynamics, is a widely observed phenomenon. Such dynamics has been studied in various ways. The celebrated Lieb-Robinson bound [19] provides an upper limit on the velocity of propagation of a local disturbance in systems with short-range interactions, and several important problems have been solved by its application [20][21][22]. While this rigorous result and seminal works offer an intuitive explanation for the * yoshi9d@iis.u-tokyo.ac.jp light-cone behavior [23][24][25], the propagation dynamics exhibits phenomena of a wide variety depending on the situation. For instance, the actual velocity of information propagation depends not only on the local Hamiltonian as the Lieb-Robinson velocity does, but also on the band structure of the total Hamiltonian and the initial state [26,27]. Indeed, we show below that a local impact propagates much slower than in the velocity given by the Lieb-Robinson bound.
The most common setup of the Hamiltonian and the initial state to study information propagation is a protocol that we refer to as the global quench, in which one prepares the ground state of a given Hamiltonian and suddenly and forever changes (namely "quenches") global system parameters, such as the interaction strength and a magnetic field [5,28,29]. Light-cone-like propagation of information has been observed under this protocol in a wide range of systems regardless of their integrability, mainly by calculating two-point correlation functions, entanglement entropy, and out-of-time ordered correlations, and the importance of information propagation in the relaxation process has been discussed [15,[30][31][32][33][34]. The global-quench protocol is also used to explain the dynamics and the speed of propagation in integrable systems, in terms of a quasiparticle picture in which a pair of correlated particles are emitted from each point on the chain after a quench and propagate with the maximum group velocity of the quasiparticles [32,35,36].
Considering inhomogeneous initial states is another way of investigating propagation dynamics. A protocol that is often considered in these decades is to connect the edge of two chains in different phases and producing an initial state with a domain wall. In this case, the energy and magnetization as well as the correlation functions exhibit propagation dynamics unlike in systems with homogeneous initial states. In the XX spin chain and the transverse-field Ising chain, a universal behavior characterized by a scaling t −1/3 of time dependence has been revealed [37][38][39] as well as a staircase structure of magnetization profiles [40]. In fact, similar scaling has been also found in the asymptotic behavior around the light-cone edge of the correlation functions in the global-quench setting [27,41,42], since the Airy function, which is characterized by the time dependence t −1/3 , is used in both cases. Another inhomogeneous initial conditions considered in this context is a local excitation, or a so-called droplet-like initial state, in which a few sites of a homogeneous system is perturbed [43,44]. In the XXZ model, a light-cone dynamics with multiple wave fronts has been found under local excitation [45,46]. While it can provide a further insight in understanding propagation properties of quantum systems, such a locally excited situation has been addressed in few studies so far [47][48][49], compared to the usual global quenches.
In the present paper, we investigate propagation dynamics of perturbation in an integrable quantum spin chain after locally disturbing the system. We first prepare the ground state of a given system, and then apply a unitary operation U loc which acts on the state only locally and instantaneously. The state after our protocol is therefore written as |ψ(t) = e −iHt U loc |ψ , with the initial state |ψ being an eigenstate of H. We refer to the localized unitary operation U loc as a local impact since it can be described by an instantaneous change of local parameters of the Hamiltonian, i.e., H(t) = H + V δ n,0 δ(t) with U loc = e −iV δn,0 , where we define the position of the disturbed site as the origin n = 0.
Important characteristics of this protocol as opposed to other quench protocols include the following two points: (i) we can observe a light-cone-like propagation of quasiparticles in terms of local observables unlike in global quenches which generally yield no transport in the dynamics; (ii) and that the translational invariance and the integrability of the Hamiltonian are conserved after the local impact in contrast to a local quench where the Hamiltonian is locally changed forever [50].
We specifically consider the spin-1/2 anisotropic XY chain in a magnetic field, and calculate the dynamics of the magnetization [28,51]. We find that the model exhibits rich propagation dynamics of the wave front, such as the existence of a second wave front and power-law decay with several exponents depending on the model parameters. For the asymptotic behavior of the wave fronts, the Airy function has been widely used in the previous studies for integrable systems. We perform an asymptotic analysis by generalizing the Airy scaling techniques, and demonstrate that it successfully captures the long-time behavior of the wave fronts in most cases. We also show that this technique fails when the model reduces to the Ising model, or when the system is on the Ising transition line.
This paper is organized as follows. In Sec. II, after introducing the model, we derive an integral form of the magnetization change under the local impact protocol and perform an asymptotic analysis to find the velocity of the propagation. In Sec. III, we present the phase diagram according to the inflection points of the dispersion relation, or the local extrema of the group velocity of quasiparticles and investigate the propagation dynamics by numerical integration as well as analytical calculation. We conclude the paper in Sec. IV, summarizing our findings and proposing future research. We also provide appendices to show details of calculation.

II. ANALYTIC CALCULATION OF THE TIME EVOLUTION OF THE MAGNETIZATION CHANGE
We consider the one-dimensional spin-1/2 antiferromagnetic XY model described by the Hamiltonian where {S x n , S y n , S z n } are the spin-1/2 operators, L denotes the system size, γ denotes the XY anisotropy, and h denotes a magnetic field. We require the periodic boundary conditions S L+1 = S 1 and take the system size L to be an even number in the diagonalization below. In this study, we particularly investigate the dynamics in the parameter region of 0 ≤ h ≤ 2 and 0 ≤ γ ≤ 1.
In our local-impact protocol, we here use the ground state |GS for the initial state and specifically give the local impact U loc = e iθS z 0 , namely a rotation over the angle θ of S 0 around the z-axis. Then we analyze the spacial propagation of the effect of the local impact on the state by calculating the dynamics of the magnetization in the z direction at each site n, according to the original Hamiltonian (1). We focus on the change of the local magnetization where S z n (t) denotes e iHt S z n e −iHt . (We set = 1 here and hereafter.) The Jordan-Wigner transformation, which we introduce later, makes Eq. (2) equal to the change of the fermion density at site n. In the calculation of the propagation dynamics, we take the thermodynamic limit L → ∞.
In this section we first give a brief summary of the diagonalization of the XY model in one dimension under the periodic boundary condition. After that we derive an integral expression of the magnetization change ∆(S z ; n, t) and perform an asymptotic analysis in order to discuss the velocity of the propagating wave fronts. All the results on the propagation dynamics in this study also hold for the ferromagnetic XY model.

A. Diagonalisation of the XY model in one dimension
For the diagonalization of the Hamiltonian (1), we rewrite it in terms of spinless fermions by using the Jordan-Wigner transformation [52], which is defined by where the operators c n , c † n obey the fermionic anticommutation relations {c n , c m } = {c † n , c † m } = 0, {c † n , c m } = δ nm . We thereby obtain where the boundary condition is given by The operator e iπN has the eigenvalues ±1 and commutes with the Hamiltonian (4).
We can therefore block-diagonalize the Hamiltonian as with where P ± ≡ 1 2 (1 ± e iπN ) with P + + P − = 1 (8) are the projection operators onto the respective blocks, which commute with H and H ± as well as S z n for all n. The blocks given by P ± are sometimes referred to as the Neveu-Schwarz sector and the Ramond sector, respectively [53].
The dispersion relation ε(p) of the quasiparticles in Eq. (9) is given by for the anisotropic case γ = 0. For the isotropic case γ = 0, it reduces to and hence we have s p = 1 and t p = 0. The dispersion relation (13) can have a multimodal shape as we show in Fig. 1(a). The group velocity of the quasiparticles is given by for the anisotropic case γ = 0, whereas v g (p) = − sin p for the isotropic case γ = 0. We here use the ground state of the Hamiltonian (1) for the initial state of our local-impact protocol. The ground state of Eq. (1) is given by either or both of the ground states |GS ± of the Hamiltonians (9), where the sign of the subscript of the ground states corresponds to the one of the superscript of the Hamiltonians. In fact, the choice of the ground state of Eq. (1) depends on L, γ, and h as discussed in Ref. [54]. Nevertheless, whether we choose |GS + , |GS − , or a superposition of them as the ground state of H, is not relevant in the calculation of ∆(S z ; n, t) for L 1, which we show in Appendix A.
For brevity, we here describe the derivation of only |GS + . For the anisotropic case, the ground state of H + is given by the vacuum of η p since ε(p) > 0 for a finite even L: For the isotropic case, the ground state is the state in which only the levels with negative energies are filled with fermions η p : where we assumed for simplicity that no momentum p satisfies ε(p) = cos p+h = 0. If there is a value of p with ε(p) = 0 and the ground state has degeneracy owing to this zero-energy excitation, it would only make difference of O (1/L) in the magnetization change (2). We use Eqs. (16) and (17) in deriving Eqs. (25)-(32) from Eqs. (21) and (24) in Sec. II B.
Finally we discuss the role of the local impact regarding the quasiparticle excitation on the ground state. The local impact that we use here is expressed as follows in terms of η p : The third term on the right-hand side of this expression represents the creation and annihilation of the quasiparticles. When this term is applied to the ground state (16) for the anisotropic case γ = 0, it excites all possible pairs of quasiparticles with momentum −p and q since the ground state is the vacuum. For the isotropic case, this term is reduced to In this case, the local impact excites all possible pairs of quasiparticle excitation and hole on the ground state (17) since it is occupied by the quasiparticles below the Fermi level.
In both cases, the local impact excites quasiparticles with broad range of energies. As a consequence, information of the local impact is ballistically transferred by quasiparticles of all possible momenta, and the fastest quasiparticles form propagating wave fronts of a light cone, regardless of the detail of the local impact. Even if the impact is weak in the sense that U loc ∼ I, i.e. θ ∼ 0, the above picture holds, and the excitation is not limited to low-energy quasiparticles.

B. Time evolution of the magnetization change
Now we present an analytical expression of the magnetization change (2): where the angular brackets · · · denote the expectation value with respect to the ground state of our choice and with We here note that the anti-commutation relations on the righthand side of Eqs. (21) and (22) are actually c-numbers; see Appendix A.
For the anisotropic case γ = 0, we obtain the analytic expressions of the functions F, G, Q, and W as follows by using the quasiparticle expression (10) in the thermodynamic limit: where For the isotropic case γ = 0, the functions G and W vanish, and the function F is represented by the Bessel function of the first kind J n (z): Q(n, t) = {p : cos p+h≤0} dp 2π e i(cos p+h)t+ipn .
We provide an outline of the derivation of these expressions in Appendix A. The derivation of Eq. (19) can be generalized to other spinchain Hamiltonians which are mapped into quadratic fermion systems by the Jordan-Wigner transformation as well as for initial states other than the ground state, including a finitetemperature thermal equilibrium state. For a thermal initial state with the temperature β, we replace the integrals π −π dp in Eqs. (27) and (28), and {p : cos p+h≤0} dp in Eq. (32) with π −π dp (1 + exp(βε(p))) −1 .

C. Asymptotic analysis and the velocity of propagation
The propagation velocity of the magnetization change is well characterized by the group velocity of the quasiparticles that are emitted from the impacted site. From Eqs. (19) and (25)-(32), we can expect that the dominant component of the wave front propagates with the group velocity at the local extrema. Here we roughly explain it by approximating the integrals F, G, Q, and W in the space-time scaling limit, that is, for a large time The functions in Eqs. (25)-(28) have the following integral form in common: where g(p) is a continuous function for p ∈ [−π, π]. Since the magnetization change ∆(S z ; n, t) is expressed by a quadratic sum of the integrals F, G, Q and W in Eqs. (25)-(32), we can estimate the behavior of ∆(S z ; n, t) by investigating the asymptotic behavior of Eq. (33). For a large t with v = x/t fixed, the leading contribution is obtained from the integral around a stationary point p * at which or ε (p * ) = v holds [55]. Then we expand ε(p) around p * as where and we assumed g(p * ) = 0 and v ≤ max p v g (p) so that the stationary point p * may exist. We perform the Fresnel integral to obtain (We present more precise approximation in Sec. III B 1 and Appendix C.) This shows that the integral (37) generally decays as t −1/2 except that it decays slower than t −1/2 when we choose v to be the group velocity v g (p) at one of its local extrema, where the corresponding stationary point p * satisfies ε (p * ) = 0, and thereby κ ≥ 3. Therefore the integral (33) yields wave fronts which propagate with the group velocity at its local extrema, forming the profile of a light cone and standing out from the bulk inside the light cone.
For the anisotropic case γ = 0, the dispersion relation (13) can have two inflection points in 0 < p ≤ π for some parameter regions, which means that v g (p) can have two local extrema in 0 < p ≤ π. (We only describe the inflection points in 0 < p ≤ π hereafter since the dispersions (13) and (14) are even functions of p.) In this case, there generally appears two wave fronts propagating with the velocities V 1 ≡ |v g (p * )| and V 2 ≡ |v g (p * * )|, where p * and p * * denote the inflection points as in ε (p * ) = ε (p * * ) = 0, and we assumed V 1 ≥ V 2 without loss of generality. The second velocity V 2 is defined only when the dispersion relation ε(p) has two inflection points in 0 < p ≤ π.

III. LIGHT-CONE DYNAMICS IN VARIOUS PHASES
In this section, we calculate the magnetization change (19) by numerical integration of Eqs. (25)-(27) and (32), and investigate the propagation dynamics under the local-impact protocol analytically. For the model parameters, we mainly investigate the region 0 ≤ h ≤ 2 and 0 ≤ γ ≤ 1. We particularly present the results for the local impact U loc = e iθS z 0 with θ = 2π/3, while the choice of θ makes only subtle change in the propagation dynamics because quasiparticles with any p are excited anyway as we stressed at the end of Sec. II A.
We observe that the local impact creates a ballistically propagating wave fronts, forming a light cone, except for the case of γ = 0, h ≥ 1, in which no dynamics is obtained since the ground state becomes an eigenstate of the local impact U loc = e iθS z 0 as well as S z 0 (see the expectation value of S z 0 for γ = 0 in Fig. 1(b)), and for the case of γ = 1 and h = 0, namely when the model reduces to the trivial Ising model H = L n=1 S x n S x n+1 , in which the local impact only causes an oscillation in S z 0 and does not spatially propagate, i.e. ∆(S z ; n = 0, t) = 0. We do not consider these exceptional cases hereafter.  (13) in the parameter region 0 ≤ h ≤ 2, 0 < γ ≤ 1. The region A represents the parameter region h > |1 − γ 2 |, excluding the orange thick line for h = 1, on which ε(p) has only one inflection point in 0 < p < π. The regions B, C, and D represent the parameter regions h = 1 − γ 2 (the black thick carve), h < 1 − γ 2 (the region with orange vertical lines), and h = 1 with 0 < γ < √ 3/2 = γc (the orange thick line), respectively. In B, C, and D, the dispersion ε(p) has two inflection points in 0 < p ≤ π. On the chain line E (the region γ = 1 with 0 < h < 1) and on the Ising transition line h = 1 (the broken line), a relation cos p * + h = 0 holds for one of the inflection point p * in 0 < p ≤ π (see Sec. III B). On the isotropic line γ = 0 (the dotted line), the dispersion relation is given by Eq. (14), and its inflection points are p = ±π/2.
First, we provide a phase diagram according to the number of inflection points, which is relevant in investigating the propagation of quasiparticles, and then present some results obtained by numerically integrating the functions (25)- (28) with particular interest in the speed of the propagating wave fronts. After that we investigate the relation between the long-time behavior of the amplitude of the wave front and the derivatives of the dispersion, and show that the wave front decays in a power law in time with exponents depending on which point of the phase diagram the model is located at.

A. Phase diagram and the propagation dynamics
As we have explained in Sec. II C, the number of inflection points of the dispersion relation generally corresponds to the number of propagating wave fronts. First we show in Fig. 2 the phase diagram according to the number of inflection points in 0 < p ≤ π. We also provide in Fig. 3 plots of the group velocities v g (p) = ε (p) for some parameter sets in the regions in the phase diagram. Note that the inflection points of the dispersion relation correspond to the local extrema of the group velocity. For |h| ≤ |1 − γ 2 | with γ = 0 (the regions B and C in Fig. 2), and for h = 1 with 0 < γ < √ 3/2 = γ c (the region D in Fig. 2), the dispersion ε(p) has two inflection points in 0 < p ≤ π, whereas in the other cases it has only one inflection point. Now we provide results of the dynamics of the magnetiza-  tion change ∆(S z ; n, t) under the local-impact protocol. We obtained the dynamics of ∆(S z ; n, t) in the thermodynamic limit by numerical integration of Eqs. (25)- (32). Figure 4 shows the dynamics of |∆(S z ; n, t)| for four sets of the model parameters in the regions of A, B, C, and E in the phase diagram. In all panels, the effect of the local impact spreads linearly in time, forming a light cone, and the magnetization change is exponentially suppressed at the sites n outside of the light cone. We also provide the profiles of the magnetization change ∆(S z ; n, t) at time t = 400 in Fig. 5 for eight parameter sets. The vertical line which indicate V 1 t and V 2 t agree well in most cases with the peak positions of the magnetization change that we obtained from numerical integration. This confirms the validity of our analysis in Sec. II C.
In most cases, a pair of wave fronts propagates ballistically with a clear peak, forming a light cone, as is exemplified in Figs. 4(a)-4(c). All panels of Fig. 5 except for Fig. 5(f) demonstrate that the first peak position agrees well with the red vertical line, which indicates V 1 t from the analysis in Sec. II C.
In addition, there emerges another pair of wave fronts inside the light cone when the parameter set is located in the region C in the phase diagram, as is exemplified in Fig. 4(c).  On the line B, the second set of wave fronts on the right and left sides merge in the middle to create a ridge at n = 0 as in Fig. 4(b), which we refer to as a "frozen" wave front. This is consistent with the fact that V 2 = 0 on the line B. We show below that the frozen wave front decays slower than the first wave front, as we can observe in Fig. 4(b).
When the parameter set is located on the line E in Fig. 2, where the XY model reduces to the transverse Ising model with the magnetic filed 0 < h < 1, we observe no clear peak around the wave front, as is shown in Fig. 5(f), whereas a peak appears around the wave front for γ = 1 with h ≥ 1 as shown in Figs. 5(g) and 5(h). We will reconsider the behavior in Fig. 5(f) below in Sec. III B 3.
On the line D in Fig. 2, a second local extremum of v g (p) emerges, and hence we would expect the appearance of a second wave front as is the case of the region C, but it is in fact hard to identify in Figs. 5(c) and 5(d). In Sec. III B 3 we discuss the origin of this behavior analytically. The second extremum disappears at the upper end of the line D, and hence we obtain a single wave front in Fig. 5(e).
The maximum group velocity V 1 and the group velocity at its second local maximum V 2 are shown in Fig. 6. As we observed above, they mostly give good estimates of the location of the wave fronts. We obtained the velocities by numerically searching the inflection points of ε(p) in the parameter region 0.05 ≤ γ ≤ 1, 0 ≤ h ≤ 2. In the isotropic case γ = 0, namely on the isotropic line in Fig. 2, the maximum group velocity is always unity and there are no second local maxima as in Fig. 5(a) because the dispersion relation (14) is ε(p) = cos p + h. For h = 0 with 0 < γ < 1, the second velocity V 2 coincides with the first one V 1 and hence there appears only one wave front even though the number of inflection points in 0 < p ≤ π is two. Incidently, we show in Appendix B that the Lieb-Robinson velocity is much faster than V 1 and V 2 . The parameter dependence is also essentially different from the one in Fig. 6.

B. The profile of the wave fronts
We now focus on the long-time behavior of the wave fronts. Figure 7 shows the time dependence of the amplitude of a wave front of the magnetization change for six modelparameter sets. They show power-law decay with various exponents. Below we discuss the origin of these exponents by using a stationary phase analysis. 1. The decay t −2/3 on the wave front As we have explained in Sec. II C, the integral (33) decays as a power law in time t as t −1/κ in the space-time scaling limit with the integer κ determined by Eq. (36). We can estimate the decay of the wave front of the magnetization change |∆(S z ; n, t)| as t −2/κ since it has a quadratic form of the integrals of the form (33). Focusing on the parameter region 0 < γ ≤ 1, 0 ≤ h ≤ 2, we find that κ is three except for the case of h = 1 − γ 2 , in which κ becomes four with p * = π, and for the case of γ = γ c , h = 1, in which κ becomes five with p * = π.
We can approximate the profile of the wave front for large t by extending the asymptotic analysis in Eq. (37). Around the wave front i.e. n ∼ vt with v = V 1 and V 2 , the integral (33) can be approximated by (38) as long as (see appendix C for the derivation), where we define when ε (n) (p * ) > 0. When ε (n) (p * ) < 0, on the other hand, we change B n (X n ) in Eq. (40)  ε (p * 2 ) = ..., we add up all the contributions from these points, i.e. A n (p * 1 , x, t) + A n (p * 1 , x, t) + ... .) Using the approximation (38) for Eqs. (25)-(28), we obtain as the leading behavior of Eq. (20) for large t with v = V 1 and V 2 , while the next-order term in this approximation is estimated at O t −1/κ × O t −2/κ = O t −3/κ as a crossing term from the first and second terms in Eq. (38). The approximations (43) and (44) are useful as long as cos p * + h = 0. Not only that the expressions (38) and (40) show that the integral I(x, t) decays as t −1/κ for large t with n ∼ vt as we derived in Sec. II C but also that they well reproduce the profile of the magnetization change off the wave front. The integral B n (X) can be seen as a generalization of the Airy function of the first kind since B 3 (X) = Ai(X). Figure 8 demonstrates a good agreement between the numerical calculation of the magnetization change profile and the approximation obtained from Eq. (38) with κ = 3 for Eqs. (25)- (28).
Although the agreement becomes worse for the second wave front compared to the fastest wave front due to the interference from the tail of the fastest one, the approximation succeeds in capturing the location and the height of the two fronts.
This kind of analysis has been performed in several studies for the XX model (γ = 0), the Ising model (γ = 1) and the Bose-Hubbard model, for instance. In Refs. [27,41,42], the wave front of correlation functions after global quenches are argued to be well described in terms of the Airy function Ai(X) = B 3 (X). In Refs. [37,38], the Airy function is also used to characterize the wave fronts after quenches from steplike inhomogeneous initial states. On the other hand, to the best of our knowledge, the asymptotic behavior of wave fronts has not been carefully investigated for the XY model with 0 < γ < 1 so far.
2. The decay t −1/2 at the origin When the model-parameter set is located on the phase boundary B with h = 1 − γ 2 , there emerges a "frozen" wave front (see Fig. 4(b)) which decays as t −1/2 , in addition to the propagating wave front, which is well described by the use of the Airy function in Eqs. (43) and (44) with κ = 3. In this case, the dispersion relation ε(p) has an inflection point at p * = π (see local extrema of ε (p) in Fig. 3(a) with h = 0.75), at which the group velocity v g (p) = ε (p) as well as the third derivative of the dispersion vanish, while the forth derivative of the dispersion is given by d 4 ε(p = π)/dp 4 = 3/γ 2 − 3. Therefore, ε (4) (p = π) is finite and κ for this inflection point is four except for the case of γ = 1 and h = 1 − γ 2 = 0, at which the dispersion becomes constant, i.e., ε(p) = 1. We thereby observe that the magnetization change on the frozen wave front at the impacted site n = 0 remains for a long time when the model-parameter set is on the line h = 1 − γ 2 .
The result in Fig. 7(e) demonstrates that the amplitude of the frozen wave front decays as t −1/2 with an oscillation, owing to the decay t −1/4 of κ = 4 at n = 0 of the functions F, G, Q, and W . Since the decay t −1/2 is slower than t −2/3 of the first wave front, the frozen wavefront stands out as in Fig. 4(b). In the other parameter regions in the phase diagram, the magnetization change decays as t −1 at n = 0 as is exemplified in Fig. 7(d) (at which the model is located in the region A), owing to the decay t −1/2 of κ = 2 with an oscillation of the same functions.
3. The decay t −1 and t −3/5 in special cases So far we have discussed the cases in which the long-time dynamics of the wave fronts can be well described by the approximation (43) and (44). The coefficient cos p * + h in Eqs. (43) and (44) vanishes when the model-parameter set is located on the line E and on the Ising transition line h = 1 in the phase diagram. In these parameter regions, the wave fronts show anomalous behavior.
For h = 1, the dispersion relation ε(p) has an inflection point at p * = π. Since the right-hand side of Eqs. (43) and  Fig. 2. The blue points represent numerical result, while the orange line represents the approximation using Eq. (38) with κ = 3. The magnetization change ∆(S z ; n, t) has two wave fronts for n > 0 owing to two inflection points p * and p * * of ε(p) with vg(p * ) > vg(p * * ) > 0. In the approximation, we estimate each integral of (25)- (28) as the sum of the integrals around p * and p * * , and then used the the approximation (38) with κ = 3 for each inflection point to calculate the approximated value of ∆(S z ; n, t).
(44) vanish as in (cos π + 1)/ε(π) = (cos π + 1)/2γ 2 = 0, the long-time behavior of the wave front owing to this inflection point is given by a higher-order term in the approximation. We numerically find the decay t −1 for h = γ = 1 as we show in Fig. 7(c), which is consistent with our estimation on the time dependence of the second-order term in the approxi-mation (43) and (44), namely O t −3/κ with κ = 3.
For 0 < γ < γ c = √ 3/2, namely on the line D, the dispersion ε(p) has another inflection point in 0 < p < π; see the local extrema of ε (p) in Fig. 3(b). The wave front corresponding to this new inflection point propagates faster than that of p * = π and decays as t −2/3 , forming a light cone (see Figs. 5(c) and 5(d)). In this case, the second wave front inside this light cone region is harder to identify compared to the case in the region C (see Fig. 5(b)) because the amplitude of this wave front is expected to decay as t −1 which is faster than the decay t −2/3 of the fastest wave front. Note that ∆(S z ; n = vt, t) ∼ t −1 typically holds since the integral I(vt, t) in Eq. (33) behaves as I(vt, t) ∼ t −1/2 when |v| is neither V 1 nor V 2 .
At the point h = 1 and γ = γ c , namely at the upper edge of D in Fig. 2, another inflection point collapses with the inflection point at p = π. At this point, the third and fourth derivatives of the dispersion vanish at p * = π, while the fifth derivative ε (5) (p = π) survives. Therefore, despite the leading behavior of the integral Eq. (38) is expected to decay as ∼ t −1/5 , the wave front of the magnetization change shows the decay ∼ t −3/5 as we have observed in Fig. 7(b), since the coefficients in Eq. (43) and (44) vanish at p * = π. Again this decay is consistent with our estimation of the next order of (43) and (44), O t −3/κ with κ = 5.
In the case with γ = 1, 0 < h < 1, namely when the XY model reduces to the transverse Ising model, the coefficients in Eqs. (44) and (43) again vanish. In this case, we observe that there is no clear peak around the wave front as is exemplified in Fig. 5(f), whereas there appears a clear peak for h ≥ 1 as is exemplified in Figs. 5(g)-5(h), and in Fig 7(c). The behavior for γ = 1, 0 < h < 1 as well as h = 1 is considered to be described by higher-order terms in the approximation (44) and (43), whose exact form we have not succeed in obtaining analytically.

IV. DISCUSSION
In this paper we have investigated the propagation dynamics in the one-dimensional XY model under a magnetic field. We introduced the local-impact protocol, which is described by local and instantaneous unitary operation U loc applied to a steady state, and focused on the velocity of the propagation and the asymptotic behavior of the amplitude of the propagating wave front. We found distinctive features of the profile of the magnetization in the XY model, which mediates two prototypical integrable models, the XX chain (γ = 0) and the transverse field Ising chain (γ = 1), particularly in the antiferromagnetic phase |h| < 1 as well as in the critical phase |h| = 1.
Using numerical calculation and analytical computation, we demonstrated that the model exhibits a second wave front inside the light-cone region for |h| ≤ 1 − γ 2 , namely in the regions B and C in the phase diagram; see Figs. 4(b) and 4(c). This second wave front only emerges for γ = 0 since it originates from multiple local extrema of the group velocity of the quasiparticles, which can appear only for γ = 0 (more specifically, |h| ≤ 1 − γ 2 ). In the global-quench protocol, which has been a major protocol for investigating nonequilibrium dynamics and information propagation in isolated quantum systems, the second wave front would be blurred by all waves from other points. Nevertheless, Ref. [56] numerically observed that quasiparticles with the mode at the second local extremum of v q (p) can carry a dominant part of information which contributes to the entanglement growth in a global-quench setting.
We also found that the profile of the magnetization change exhibits drastic difference, that is the absence of a peak around the wave front (see Fig. 5(f)) for γ = 1, namely on the line E in Fig. 2. While we provide an analytical description for the origin of this behavior in Sec. III B 3, it will be interesting to find a physically relevant explanation using a quasiparticle picture, as well as investigating the universality of this difference in terms of other observables.
The transition line |h| = |1 − γ 2 | from the viewpoint of the dynamical behavior of the XY model has been identified in some other studies including work on a non-equilibrium steady state [57] and on the relaxation of the magnetization after a global quench [28]. In our protocol, on the other hand, we can capture this transition line by simply observing the dynamics of the "frozen" wave front around the impacted site after applying a local unitary operation to the system. Our results suggest that observing propagation dynamics of the local disturbance in terms of a local spin magnetization can solely show a rich and nontrivial behavior of dynamical properties of quantum systems.
Several studies have investigated long-time behavior of correlation functions around the light-cone edge under quench protocols. The Airy function associated with the scaling t −1/3 has been used to describe the dynamics around the wave front in order to discuss the asymptote of its height, width and velocity [27,37,38,41,42], and the scaling t −1/3 and its square t −2/3 appeared universally in light-cone dynamics. In the present paper, in contrast, we have revealed using the localimpact protocol that the scaling for the height of the wave front around the light-cone edge can be given not only by t −2/3 by t −1 and t −3/5 (see Fig. 7), by carefully investigating the dispersion relation ε(p) and the coefficient for the approximation. In particular, we found that the leading terms (44) and (43) in the approximation of ∆(S z ; n, t) vanish when the model is on the line h = 1 or on the line 0 < h < 1 and γ = 1, for which the relation cos p * + h = 0 holds for a local extremum p = p * of v g (p).
The local-impact protocol which we introduced in this paper may provide a new insight into the study of dynamics in isolated quantum systems. It will be important to investigate the propagation dynamics in this protocol in terms of other observables, such as the magnetization in the x and y directions and the entanglement entropy as are studied in [47,49] for the transverse Ising model. Studying a relaxation process after applying the local impacts for all sites is an interesting direction for future research. Therefore, we can parallelly calculate ∆(S z ; n, t) for the ground state in the two sectors. We note that the ground states of the XY model can be |GS + , |GS − , or a superposition of them depending on the size L, the anisotropy γ, and the field h; see Ref. [54].
We first show that anti-commutators {c † n (t), c m } and {c n (t), c m } are c-numbers, and then derive the expressions of F , G, Q, and W in Eqs. (25)- (32). From the equation of motion of quasiparticles η p , i.e., (d/dt)η p (t) = i[H a , η p (t)] = −iε(p)η p (t) with respect to each sector a = ±, we obtain η p (t) = η p exp [(−iε(p))t], and thereby obtain the expression of the Jordan-Wigner fermions c n (t) in terms of η p as from Eq. (10). Using Eq. (A3) and the relations s p = s −p = s * p , ε(p) = ε(−p) and t p = −t −p = −t * p , we obtain Equations (A4)-(A8) clearly show that the anti-commutators {c † n (t), c 0 } and {c n (t), c 0 } are c-numbers, which we denote F = F (n, t) and G = G(n, t), respectively, as in Eqs. (21) and (22). In Eqs. (A7) and (A8), we used an expression c n (t) = ( 1/L) p a η p exp [−iε(p)t − ipn] for γ = 0 since s p ≡ 1 and t p ≡ 0 from the definitions (11) and (12). We ob-tain the expression in (19) with (20) by utilizing the fact that Eqs. (A4) and (A5) are c-numbers. We calculate Q and W in a similar way, additionaly by taking into account Eqs. (16) for γ = 0, and Eq. (17) for γ = 0 in calculating the expectation values with respect to the ground state: a GS| c † n (t)c 0 |GS a = 1 L p,q a s p s q e iε(p)t a GS| η † p η q |GS a − t p t q e −iε(p)t a GS| η −p η † −q |GS a e ipn = 1 L p a − t 2 p e −iε(p)t e ipn = 1 L p a |t p | 2 Φ * p (n, t), a GS| c n (t)c m |GS a = 1 L p,q a s p t q e −iε(p)t a GS| η p η † −q |GS a + t p s q e iε(p)t a GS| η † −p η q |GS a e −ipn for the isotropic case γ = 0.
We arrive at the expressions (25)-(32) by taking the thermodynamic limit L → 1 to replace the sum p a over p = (2πj/L with j = −(L − 1)/2, ..., −1/2, 1/2, ..., (L − 1)/2 for the Neveu-Schwarz sector a = + and j = −L/2 − 1, ..., −1, 0, 1, ..., L/2 for the Ramond sector a = −) with the integral π −π dp over p. Now we show that the choice of the ground state is not relevant to the calculation of the magnetization change ∆(S z ; n, t), i.e. the difference between and can be ignored in the thermodynamic limit. The difference only comes from the way in which we take sums over p before taking the thermodynamic limit in order to obtain the integral representations in Eqs. (25)- (28), (32), and (30). Since the correction for replacing a discrete sum over p for the integral over [−π, π] is estimated at O L −1 , the difference is irrelevant in the thermodynamic limit. Therefore, we do not have to specify which sector the ground state of (1) belongs to in the calculation of the magnetization change ∆(S z ; n, t). As we mentioned in Introduction, the Lieb-Robinson bound [19] provides a bound for the velocity of the information propagation in lattice spin systems with local interactions. However, the Lieb-Robinson velocity depends only on the operator norm of the local Hamiltonian, particularly in onedimensional systems with nearest-neighbor interaction [58]. The characteristic velocity for the propagation dynamics in a given system can generally depend nontrivially on the property of the system. The Lieb-Robinson velocity specifically for the XY model is given by v LR = 2e H n = max e, e h 2 + γ 2 , where the local Hamiltonian H n (i.e., H = n H n ) of the XY model is H n = (1 + γ)S x n S x n+1 + (1 − γ)S y n S y n+1 + (h/2)(S z n + S z n+1 ), and e is the constant e = exp(1) ∼ 2.718... In Fig. 9, we compare the Lieb-Robinson velocity v LR and the maximum group velocity V 1 (see Fig. 6(a)) in terms of the dependence on the model parameters. We can see that V 1 is much less than v LR and the dependence on the model parameters is also different.