Entanglement Evolution in Lifshitz-type Scalar Theories

We study propagation of entanglement after a mass quench in free scalar Lifshitz theories. We show that entanglement entropy goes across three distinct growth regimes before relaxing to a generalized Gibbs ensemble, namely 'initial rapid growth', 'main linear growth' and 'tortoise saturation'. We show that although a wide spectrum of quasi-particles are responsible for entanglement propagation, as long as the occupation number of the zero mode is not divergent, the linear main growth regime is dominated by the fastest quasi-particle propagating on the edges of a widen light-cone. We present strong evidences in support of effective causality and therefore define an effective notion of saturation time in these theories. The larger the dynamical exponent is, the shorter the linear main growth regime becomes. Due to a pile of tortoise modes which become dominant after saturation of fast modes, exact saturation time is postponed to infinity.


Introduction
Propagation of entanglement in systems with a large number of degrees of freedom is of great importance to understand non-equilibrium physics in such systems [1,2]. This has been widely studied via considering a global quench and focusing on the evolution of a given initial state. The evolution of such a system due to the quench is generally a thermalization process. A typical measure for this thermalization process is how the reduced density matrix corresponding to a typical spatial subregion is close to a thermal density matrix.
The evolution of the system after the quench on the other hand is also an equilibration process. Since the post-quench evolution is a unitary evolution, the system will only reach local equilibrium though never a global one. In this also the relevant quantity is a density matrix corresponding to a subsystem which is expected to reach local equilibrium. At the end of the thermalization process when the system has reached local equilibrium, the local observables are given by the thermal (Gibbs) ensemble.
The story is different in case of integrable models (including free systems). In these systems the observables are much more constrained by an infinite number of conserved charges and the systems does not thermalize ending up in a Gibbs ensemble but relaxes to generalized Gibbs ensemble [3,4].
A typical measure to quantify the evolution of the system in a pure state after a global quench is to study entanglement entropy of a subsystem in the post-quench state. The entanglement entropy is defined as and ρ A := TrĀ|ψ ψ| where |ψ = e −iHt |ψ 0 . |ψ 0 is the pure pre-quench state and H is the post-quench Hamiltonian. Here by quench we mean a global quench which is defined by a sudden change of a parameter in the Hamiltonian of the system at t = 0. The dynamics of the system and specifically the spread of entanglement in such systems has been widely studied in a large number of papers (see [1] and [5] and references therein). The quasi-particle picture is a core of understanding how entanglement spreads over the system after a global quench. More recently this picture has been improved in [5] which makes it valid in a wider scope.
From a field theoretic point of view, the quasi-particle picture successfully describes propagation of entanglement in relativistic field theories. This is strongly based on the notion of causality in these theories. Of course the question of entanglement propagation as a probe to study how the system relaxes to a generalized Gibbs ensemble is a very important question in non-relativistic theories. The propagation of entanglement is specifically intertwined with the causality structure together with the equilibration process of the theory under consideration.
The goal of this paper is to investigate how entanglement propagates in a field theory with Lifshitz scaling symmetry [6,7] given by t → λ z t , x → λ x, (1.1) where z is the dynamical exponent and determines the anisotropy between space and time. This kind of scaling symmetry is typical at critical points where a continuous phase transition takes place. 1 We study evolution of entanglement in the vacuum state of a translational invariant system and try to understand the corresponding physics basically by means of a zero mode analysis. The rest of the paper is organized in the following order: As the rest of the introduction section we introduce a discrete version of Lifshitz scalar field theory which we will utilize in our numerical analysis. In section 2 we review propagation of entanglement in relativistic scalar theory after a global quench. In section 3 we give analytical arguments for quasi-particle picture in these theories followed by a numerical study of evolution of entanglement entropy. We analyse our numerical studies focusing on the spectrum of these theories, the key role of tortoise modes and the quasiparticle picture.

Lifshitz-type QFTs on Square Lattice
A quantum field theory that is invariant under a Lifshitz scaling transformation, introduced in Eq. (1.1), is what we call a Lifshitz-type QFT. As the simplest example we consider a free scalar 1 Various properties of entanglement entropy is such theories in static cases has been previously studied in [8][9][10][11][12][13][14]. field in d + 1-dimensions with the following action [15] (1. 2) The above action has a Lifshitz scaling symmetry in the massless limit (m = 0). Although in this manuscript we will study both m = 0 and m > 0 theories. This is clearly a generalization of Klein-Gordon theory (z = 1) to generic z with non-relativistic scaling symmetry. The realization of Klein-Gordon theory on a lattice is the well-known "harmonic lattice". There also exists a family of models which is generalization of harmonic model on a (hyper)square lattice in generic spatial dimensions. These models are known as Lifshitz harmonic models [11,12] and realize the action given in Eq.(1.2) on a square lattice.
To have a intuitive picture of the model lets focus on 1 + 1-dimensions. In this case the model is a chain of coupled harmonic oscillators, which each oscillator is coupled to z other oscillators from each side where z is the dynamical exponent. Obviously at this level z should be considered as a positive integer.
The Hamiltonian of these models is given as follows One can easily check that this Hamiltonian reduces to the well-known harmonic model for the case of z = 1. It is straightforward to generalize this Hamiltonian to higher dimensions which we are not interested in this paper (see [11] for 3d generalization). One can also show that after a canonical transformation (q n , p n ) → ( this is a discretized version of a free Lifshitz theory on a square lattice with mass m and lattice spacing = M/K, 2 introduced in Eq.(1.2).
The Hamiltonian of this model in any number of dimensions on a (hyper)square lattice (more precisely on a d dimensional torus) can be diagonalized in a standard way which we are not going to review here leading to the following dispersion relation [11,12] where k = {k 1 , k 2 , · · · , k d } is the momentum vector, k i 's refers to the momentum components in all spatial directions and N x i refers to the number of sites in the i-th spatial direction. Different aspects of quantum entanglement in the vacuum state of these models has been studied previously [11][12][13] 3 . We will utilize the covariant correlator method to study entanglement propagation in these models. The method is briefly reviewed in the appendix section.

Entanglement Propagation in Relativistic Theories
Before getting into the question of how the relaxation of these systems is affected by non-Lorentzian dynamical exponent, in this section we briefly review the corresponding physics in relativistic field theories. The issue has been studied for 2d CFTs in a series of papers starting by [1]. Lets first focus on the simplest case which is a global quantum quenches in the context of 2d CFTs. The post-quench Hamiltonian is scale invariant while the pre-quench one is not. Thanks to the simplification due to the conformal symmetry of the post-quench Hamiltonian, EE can be worked out analytically as a function of time. It is well-known that EE exhibits a linear growth behavior up to a saturation time when a sudden transition to a saturation regime happens. The saturation time is proportional to the length of the entangling region (with a proportionality factor of 1/2) and the saturation value of EE depends on the details of the initial state. These universal features have been verified in various scale invariant models including transverse Ising spin chain in the same reference above.
In figure 1 we have summarized mainly the results of [1], where the numerical results are obtained using a harmonic lattice model. Note that although in the case of a CFT the saturation happens instantaneously, lattice simulation shows a mild transition. In order to clarify different scaling regimes of the process, in the right panel we have also numerically plotted the time deriva-tive of EE. At the very beginning there is quadratic growth regime 4,5 after which there is the well-known linear growth. The linear growth is followed by an extremely slower growth which is argued in [18] to logarithmically depend on time. 6 We will denote this regime by tortoise saturation in what follows and will also derive an analytic time dependence for this regime in a certain limit of mass quenches. Due to the mild transition, i.e. the existence of the tortoise saturation regime, saturation time is postponed [1,18]. We will elaborate on this point both in Lorentzian and in Lifshitz theories.
There are technical problems comparing between harmonic lattice and CFT results. Utilizing the harmonic lattice model, the straightforward choices is to impose periodic boundary condition which leads to a translational invariant model, although a discrete one. But this model in general suffers from IR divergence due to the existence of a zero mode. This zero mode and a pile of other slow modes which come into the game as soon as the model needs an IR cut-off are actually responsible for this mild transition. On the other hand one can impose a Dirichlet boundary condition on both ends of the lattice and sticking the entangling region to one of them in order to get around the problem. Although this makes the transition sharper since the zero mode is killed, but does not solve the problem and the price is loosing translational invariance.
The qualitative behavior of S A (t) during this relaxation process is well-known to be understood by the quasi-particle picture [1]. In this picture the entanglement between a subregion A and its complementĀ is carried by a uniform density of free streaming quasi-particles. A pair of entangled quasi-particles is created at each spatial point and the entropy between A andĀ is measured by the number of quasi-particles which a pair is in A and the other pair is inĀ. As a consequence of causality the propagation of quasi-particles constrained to be inside the light-cone such that the massless quasi-particles move along the null rays. Using this simple scenario it has been shown that the saturation time is given by t s = /(2v g ), where v g denotes the group velocity of entanglement propagation quasi-particles which in a critical theory is v g = 1.
Recently an analytic description for the time evolution of EE after a quantum quench based on the quasi-particle picture was proposed in [2,5]. 7 It was shown that considering a quasi-particle picture for spread of entanglement and also knowing the late time stationary state provided by the Bethe ansatz, one can find an analytic description for time evolution of EE. This construction works for several integrable models including non-interacting bosonic and fermionic systems. Since we will utilize this picture to understand our results, we will shortly review its main features in the following.
According to this description the general prediction for time dependence of EE in a 2d theory Figure 2: Alba-Calabrese quasi-particle picture for a relativistic free boson and three distinct modes together with the zero mode. The saturation time for each mode is shown with {t * 1 , t * 2 , t * 3 }.
for a subregion with length is given by [5] S(t) = 2t where s(k) is the entropy density and v g ≡ dω k /dk is the group velocity of the corresponding quasi-particles with momentum k.
To have a better intuitive understanding of what Eq.(2.1) means, we have illustrated the physics in figure 2. In general there is an infinite number of modes in the model and in this figure we have shown four of them to describe the physics in a relativistic model. The fastest mode is shown in green and two slower ones in orange and blue together with the zero mode in red. When a global quench happens at any spatial point, corresponding to each mode a pair of quasi-particle starts to move around back to back. Different quasi-particles have different group velocities bounded by |v g | ≤ 1. Each mode has a saturation time with regard to its group velocity given by /(2v g ). The saturation time of each mode (for a single interval) is by definition the point where the corresponding rays from the two ends of the region intersect. These are shown by {t * The zero mode saturation time is infinite. Due to this intuitive picture one can understand the physics of Eq.(2.1). At any time t, there are a number of modes which are fast enough to be saturated some time before t, these modes no more contribute to the time evolution of EE and form the second term in Eq.(2.1). At the same time the slower modes which satisfy /(2v g ) > t, still contribute to the time evolution and form the first term in Eq.(2.1). The smaller t is, the larger the number of modes contributing to the first term.
Indeed in [2] it was shown that for certain integrable models, s(k) can be fixed employing the equivalence between the entanglement and the thermodynamic entropy in the stationary state. Note that the thermodynamic entropy can be calculated from the generalized Gibbs ensemble describing the stationary state in terms of the expectation value of mode occupation number as follows where n k = n k GGE . Also note that becausen k is an integral of motion, n k can be found using the initial state, i.e., In the case of a free scalar theory one finds [28] n k = 1 4 where ω 0,k (ω k ) is the frequency before (after) the quantum quench. Now we are equipped with all we need to calculate the time dependence of EE using Eq.(2.1). In figure 3 we demonstrate the evolution of EE predicted by Eq.(2.1) and compare it with the numerical results for different values of once with m 0 = 1, m = 2 and once with m 0 = 1, m = 0. Note that as we increase the length of the entangling region, we find better agreement. For the first family where m = 0 the results for 200 number of sites is almost the same as the thermodynamics limit but for massless post-quench Hamiltonian it does not coincide with the quasi-particle picture which we believe is due to the zero mode effect. We will try to analyse this point in details in what follows.
Remarkably, there exists an interesting relation between the spectrum of the quasi-particles and the saturation time of the EE. As the concluding note for this section we explain this relation in the simplest case, i.e., the harmonic lattice. In a 2d QFT with relativistic scaling the lattice dispersion relation for the corresponding quasi-particles is given by Eq.(1.4) (setting z = 1) as follows 8 The group velocity of the quasiparticles from the above equation is Now one can find the mode with the maximum group velocity as follows where κ = 2 sin 2 kmax 2 . Regarding Eq.(2.6) and Eq.(2.7) the following comments are in order • To trace the effect of tortoise modes (those with small momenta) more precisely, we look at the k → 0 limit of the group velocity. The behaviour is as follows One should note that since we are considering a translational invariant lattice, which the dispersion relation is given in Eq.(2.5), there is an IR regulator in the model. Thus to track the effect of the zero mode in a scale invariant post-quench state one should always be careful about the order of the m → 0 and k → 0 limits. Because of the existence of the IR regulator, one should take the k → 0 limit first. Note that k = 0 is a permissible mode in this model thus whatever the IR regulator (the mass of the post-quench state) is, the zero mode has vanishing group velocity and modes with small momenta (k m) also have very small group velocities. These modes are responsible for tortoise saturation after all other fast modes are saturated.
• For a massive quasiparticle we have κ = − m 2 2 + m 1 + m 2 2 and the maximum group velocity is given by (2.8) 8 Note that in the following without loss of generality we consider the continuum limit of the dispersion relation.
which shows that the maximum velocity keeps decreasing as a function of the mass parameter.
• From the above picture in principle it is straightforward to workout the analytical behaviour of EE after saturation time. Here there are tortoise modes which do not saturate at finite time, but one can still define an effective saturation time regarding to the saturation of the fastest mode t max s . In general there is a maximum momentum which we denote by k max corresponding to each mode. In principle the time dependent part of entropy gives the time dependence of entropy after t max s . Even for the free bosonic case this is not analytically computable but in a certain limit which one considers a quench from a very heavy state (m 1) to a gapless state (m 1) the above integral can be performed analytically which leads to (for t > /2) where S 0 is the value of EE at t = /2, c 1 and c 2 are constants depending on m 0 , m and .
The key point to work out this behaviour is that the above intuitive picture lets us to think about the time dependence of EE as the time dependence of k max .
It is worth to mention that the above behaviour, although is found for a quench from highly gapped to a gapless Hamiltonian, but numerical data show that it works well even in the regime of our study where m 0 ∼ O(1) and m = 10 −6 .
• Thinking for a while about the continuum scalar field theory, the qualitative features of the previous results do not change. In this case the group velocity can be obtained from the Once again the existence of massive particles with vanishingly small momentum (tortoise modes), leads to a tortoise saturation regime for EE at late times.
• To get around these tortoise modes, one approach is to consider mode-dependent mass quench, i.e., m(k) [18]. In this case the group velocity becomes where we should impose v g (k ∼ 0) → v 0 (> 0) to prevent the excitation of tortoise modes.
Any mass function with the specific behavior of m(k ∼ 0) = m 0 + v 0 k + · · · satisfies this condition and removes these modes from the spectrum of quasi-particles 9 . Considering these family of quenches, one no longer finds a finite saturation time with no tortoise saturation regime. It is important to note that according to Eq.(2.11) the massless quasi-particles with linear dispersion relation (ω ∼ k) move along the null rays with the maximum momentum independent velocity, v g (k) = 1.

Entanglement Propagation in Lifshitz Theories
In this section we study how EE propagates in theories with Lifshitz scaling. To be more precise, we study post-quench states both with m = 0 and m > 0. We first study analytically the quasiparticle picture for Lifshitz theories, modelled by Lifshitz harmonic lattice, after which we present numerical studies of EE and discuss about the physics of propagation of entanglement in different scaling regimes.

Analytic Description
In principle since our Lifshitz field theory of interest is a generalization of Klein-Gordon theory but the spatial correlations are stretched out via the dynamical exponent, one would naturally expect a similar analysis to what we have reviewed in the previous section from [2,5] is valid in this theory. The intuitive description we discussed in the previous section, simply shows that similar to the relativistic case, the general time dependence of EE predicted in a 2d free scalar theory with Lifshitz sacling is given by Eq.(2.1). So what we need is to workout the exact dependence of n(k), s(k) and v g on the dynamical exponent. Using the dispersion relation Eq.(1.4) a strightforward computation gives Also similar to the relativistic case the entropy density in terms of the expectation value of mode occupation numbers in the initial state is given by Eq.(2.4) where we should consider the dependence of ω k and ω 0,k on z as follows where m 0 and m are the mass parameters before and after the quantum quench respectively. In 9 Note that in order to have a finite injected energy during the quench scenario we should impose m(k ∼ ∞) → 0. general the value of n(k) and s(k) increase at a given momentum as we increase the dynamical exponent. In other words the contribution of individual quasi-particles to the thermodynamic entropy of the generalized Gibbs ensemble describing the steady state increases due to a z > 1.
In the following we will present some results which the role of tortoise modes becomes extremely important. As an example we have plotted the occupation number and the entropy density in figure 4 for quenches from a fixed pre-quench state (m 0 = 1) to different values of smaller postquench mass parameters. We have considered these parameters to figure out what happens toward quenching to gapless models (m → 0). One can easily check from the expression of n(k = 0) which is n(k = 0) = 1 4 that the occupation number (and also the entropy density) diverges in the following three cases: m 0 /m 1, m/m 0 1, and z 1. In the case of scale invariant system in principle the postquench mass vanishes. But here due to the IR cut-off, this case is a special case of m/m 0 1. We will show in what follows that the numerical results deviate from the quasi-particle picture in these three regimes, although the picture works perfectly in other regimes.
In the following sections after describing the quasi-particle picture we will present numerical results where our main focus is comparing the vanishing and non-vanishing mass parameters in the post-quench state to inquire our understanding of the quasi-particle picture.

Quasi-particle Picture
As we discussed in section 2, the spectrum of quasi-particles together with the notion of causality interestingly describes the linear growth and tortoise saturation regimes. Here we explore the role of z in this scenario focusing on 2d theories in order to avoid any complication arising in case of more than one velocity component. In general we would expect the qualitative picture should be straightforwardly generalized to higher dimensions. Although the notion of causality is a totally subtle notion in these theories, we show that (except in certain cases which the tortoise mode contributions become dominant) different scaling regimes during the relaxation of the system can be perfectly described by Alba-Calabrese quasi-particle picture in presence of z > 1 dynamical exponents. Comparing to the z = 1 case there are interesting features in z > 1 which we will describe in the following.
The group velocity for quasi-particles in a Lifshitz harmonic lattice is given by Eq. the maximum group velocity can be derived by first solving dv g /dk = 0 for the maximum momenta. This can be done analytically for small post-quench mass (for z > 1) parameter as According to the above result a few comments are in order: • For z = 1 in the massless limit we have v max g = 1 which is consistent with the group velocity of quasi-particles which we reviewed in section 2. Also note that the mass correction is negative independent of z which shows that as expected intuitively, the velocity is a decreasing function of the mass parameter.
• Similar to what we discussed in Eq.(2.9), here also it is possible to extract an analytical temporal behaviour of EE for z > 1. The most straightforward case is z = 2 which the same analysis again in a limit which a Hamiltonian with a large mass is quenched to a tiny mass Hamiltonian one finds where S 0 is the value of EE at t = /(2v max g ) and v max g is given by Eq.(3.6) for z = 2, c 1 and c 2 are constants depending on m 0 , m and . The analysis is a bit hard to be extended for higher values of z because of a technical problem which one cannot easily solve for the generic time dependence of k max (t).
• For z → ∞ where the corresponding scalar theory given by Eq.(1.2) becomes strongly nonlocal the maximum group velocity diverges, i.e., • For z < 1, v max g becomes pure imaginary. Based on this we conclude that there may not be a self-consistent analytical continuation for this model for z < 1, although our results show that there should be such a continuation for non-integer z > 1. Based on this we have not considered this range of parameter space of these theories. 12 It is also worth to mention that due to null energy condition, the same constraint on the dynamical exponent arises in the holographic theories with Lifshitz scaling symmetry [37].
• In figure 5 we have plotted v max g both as a function of the dynamical exponent (z > 1) and 10 The existence of Lieb-Robinson bounds in discrete systems but with infinite dimensional Hilbert space at each site, such as Harmonic lattice model, has been proved in [30,31]. The same proof works for Lifshitz harmonic model replacing the corresponding dispersion relation. We would like to emphasize that these bounds are not strong enough to lead to a physically reasonable Lieb-Robinson velocity. 11 For a related study in non-relativistic theories see [33]. 12 Entanglement dynamics in some long-range models which resembles our mode for 0 < z < 1 has been studied in [35,36].  In the left panel there is perfect agreement between Alba-Calabrese quasi-particle picture and numerics. In the middle panel one can see that due to the zero mode effect even for = 400 numerics deviate from Alba-Calabrese picture around the offset of tortoise saturation regime. The same is correct for higher values of z although one must wait much longer to see this.
mass parameter. Obviously there is no divergence in the m → 0 limit but as expected it diverges as z → ∞.

Numerical Results: A Zero Mode Analysis
In this section we present numerical results for propagation of entanglement after a mass quench in Lifshitz theories. We consider a mass quench while the dynamical exponent as another parameter in the dispersion relation is held fix. We consider two family of quenches which basically differ in the post-quench state. In both families the state prior to the quench is the vacuum state of an infinite Lifshitz harmonic lattice with parameters (m 0 , z) where m 0 = 0. The post-quench state is again the vacuum state of an infinite Lifshitz harmonic lattice with parameters changed to (m, z). By two families we mean the post-quench state is either chosen to be massless which the system has Lifshitz scaling symmetry or massive which does not have such a symmetry. In figure 6 we present numerical results regarding to both families and compare them with the quasi-particle picture. In the left panel which the post quench mass is finite one can see perfect agreement between the quasi-particle picture and our numerical results. In the middle and right panel we have shown results for the case where the post-quench mass vanishes. On can see that in these cases since the zero mode effect become more important specifically after the short linear growth regime is finished, there is a deviation from the quasi-particle picture although the deviation is suppressed as one gets closer to the thermodynamic limit. The latter gets much harder to reach as the dynamical exponent is increased because of what we have explained previously about the occupation number of tortoise modes (see figure 4).
In these curves one can see that the larger the value of the dynamical exponent is, the faster EE grows in time. Thus the saturation value of EE is an increasing function of z as expected.  Figure 7: Numerical versus analytic results for time evolution of EE for z = 1, 2, 3. Here we set m 0 = 1 and m = 2. There is a very good matching between quasi-particle analytic results and numerics. The small difference at the early times is due to numerical instability and the mild transition in numerical results is due to the zero mode which is not captured by the quasi-particle picture. In these plots we have set the subregion size to be 100 in units of lattice spacing.
This behaviour is expected due to the enhancement of spatial correlations for larger values of the dynamical exponent (for a detailed discussion on the relation between the strength of spatial correlation and z see [11]). We have checked this for higher values of z up to z = 5 but we do not present the results here since the curves have the same qualitative behaviour in different scales. We believe that the physics is perfectly captured by these values of z. Also it is notable that we have not studied z > 5, since in order to interpret the results in the thermodynamic limit and avoid lattice effects for these higher values of z, one needs to consider much larger subregions which requires a much higher numerical cost.
In figure 7 we focus on numerical results regarding to the first family which is a mass quench from (m 0 = 1, z) to (m = 2, z) for z = 1, 2, 3. In figure 8 we do the same analysis for the second family that is a mass quench from (m 0 = 1, z) to (m = 0, z) for the same values of dynamical exponent. In both figures 7 and 8 we have presented the time derivative of EE. In figure 6 it is not easy to distinguish between different growth regimes specifically in the second family while this can be clearly seen in figure 7. We would say that figure 7, 8 carry the most important results of this paper since all essential features of propagation of entanglement in these theories can be extracted from them.
There is an initial rapid growth regime which rapidly disappears and is very hard to follow from these cures. During the rapid growth regime the system has not even reached a local-equilibrium and the scaling of entropy is expected to be understood from the only well-defined physical quantity during this regime which is the energy density of the system. We will not present a careful study of this regime here but will make some comments on the z-dependent scaling of EE in this regime in the discussion section.
The system rapidly enters its main growth regime which in general longs up to a certain point (see figure 7 and figure 8). After this point a tortoise saturation regime starts which in principle carries on up to infinite time. Both families share this property.
In general one can obviously mention that in our first family which is the generic case of  Figure 8: Numerical versus analytic results for time evolution of EE for z = 1, 2, 3. Here we set m 0 = 1 and m = 10 −6 . The matching between quasi-particle analytic results and numerics is not as good as the m 0 = 1 and m = 2 case. In this case which tortoise modes are highly occupied, the thermodynamic it is much harder to numerically reach the thermodynamic limit. In these plots we have set the subregion size to be 200 in units of lattice spacing. In this case the difference between numerics and analytical results at the early times is both due to highly occupation of tortoise modes and also numerical instability. these theories, the prediction of quasi-particle picture applies here perfectly (figure 7) while in the second family which the contribution of the tortoise modes is increased due to the increase in their occupation number there is a serious non-concurrence with the quasi-particle picture ( figure  8). In the following subsections we will interpret these results and argue that these results can be understood if there is an effective notion of causality in these theories while this is totally non-trivial specially from the field theoretic point of view.

Physical Interpretation
The quasi-particle picture which works perfectly in many cases including 2d CFT and other massive 2d theories also offers a nice understanding of entanglement propagation in Lifshitz theories. Due to our analysis in section 3 there is a large spectrum of quasi-particles with different group velocities which are responsible for propagation of entanglement.
The simplest case understood by this picture was the case of 2d CFT which there was a monochromatic quasi-particle with v g = 1. In this case the sharp transition from linear growth to saturation regime is totally clear. At the instance that the quench happens, a pair of monochromatic quasi-particles start to propagate back-to-back to each other at all spatial points and they can pass through the entangling region. At early times those pairs initiated from the spatial points near the boundary of the entangling region have one mate inside the region and the other outside. These pairs are responsible for the linear growth. After a while proportional to the length of the region, the phenomena of ingoing and outgoing of quasi-particles through the entangling region equilibrates and thus EE saturates.
In the case of massive theories where there exists a spectrum of quasi-particles with different group velocities, the situation is a bit more complex. One should consider the behaviour of different types of quasi-particles to understand the propagation of entanglement. The group velocity in this case is constrained to 0 ≤ v g ≤ 1. All quasi-particles propagate inside the light-cone. Those with maximum group velocity propagate on the null directions and the zero modes depending on the value of the mass propagate extremely slow. In presence of these extremely slow modes, the above scenario still works for all types of quasi-particles. After the effect of the fast quasi-particles, specifically the one with maximum group velocity is equilibrated, the slow modes take over the role. Due to existence of these modes there always exists pairs which one mate is inside the entangling region and the other mate is outside. Although the number of these pairs is decreasing in time, but the decrease rate is infinitely slow thus the tortoise saturation regime extend over infinite time.
The very interesting thing about Lifshitz theories is that even in the massless case, there exists a spectrum of quasi-particles with 0 ≤ v g (z) ≤ v max g (z), where v g (z) is given in Eq.(3.1) and v max g (z) is given in Eq. (3.4) and v max g (1) = 1. The quasi-particles propagate inside a widen light-cone which its maximum velocity rays are given by v max g . In figure 9 we have shown this behaviour for an entangling region for values of z = {1, 2, 3}. The slope of the blue, orange and green rays are respectively given by {±v max g (1), ±v max g (2), ±v max g (3)}. The three coloured points are the saturation times found from the fastest monochromatic quasi-particle. For a single interval entangling region, the saturation time for a monochromatic quasi-particle is found to be defined as the intersection of the world-line of two of these quasi-particles starting from the end points of the region. This is simply found to be t max s (z) = /(2v max g (z)). The interesting point here is that in Lifshitz theories there still exists an effective notion of saturation time which we denote by t * s (z). The saturation time is defined to be the instance that evolution of EE experiences a sudden transition from the main growth regime to the tortoise saturation regime. This notion is also well-defined for z = 1 case, which is specifically important when the system is massive and a spectrum of quasi-particles are propagating around. Our numerical data hints to interesting physics in Lifshitz theories. The picture is that before EE enters the tortoise saturation regime, the monochromatic quasi-particle with v max g (z) dominates the entanglement propagation. This the physical reason why figure 9 is a very good approximate for the physics of entanglement propagation in the main growth regime. In other words it is clear that in principle t max s (z) is not the same as the effective saturation time we have defined above by t * s (z). But here what we find is that these two are actually the same. This picture clearly shows that in presence of tortoise modes the precise saturation time goes to infinity, meanwhile the effective saturation time t * s (z) goes to zero (following the green, orange and blue points in figure 9 toward the t = 0 axis). In other words as the dynamical exponent increases, the dominant regime in time is the tortoise saturation regime.
The other interesting behaviour which is also understandable with the quasi-particle picture is comparing the behaviour of EE for different values of z in the tortoise saturation. As we have argued the tortoise modes are mainly responsible for this behaviour. We have shown in 3.1 that the occupation number increases while z is increased. This is the reason why the tortoise saturation regime becomes slower for higher values of z. The EE saturates at infinite time for z > 1, and these infinite times comparing to each other is larger for higher values of z.

Lattice versus Continuum
An interesting question which was one of the main motivations for this study is how propagation of entanglement is related to the causality structure of the theory. It is well-known that propagation of entanglement is governed by the dispersion relation (spectrum of the quasi-particles) and the causality structure of the theory of interest. This was clearly understood at least in the case of 2d CFTs and massive scalar theory. On the other hand the picture is also consistent with the lattice version of massive scalar theory, i.e. the harmonic lattice.
Lets more precisely compare what is going on in our lattice model and how these results are supported by a simple analysis in the continuum limit. The dispersion relation and group velocity of propagating modes in this theory is given by Clearly this group velocity does not have an extremum at least for z > 1. This makes it highly counter-intuitive to have a causal structure in these theories. On the other hand our results which are strongly supported by the quasi-particle picture (modulo the states which the zero modes are highly occupied) show that there is at least an effective notion of causal structure due to propagation of entanglement in these theories. It is worth to note that we are comparing the discrete model in the thermodynamic limit with the continuum theory. This is actually what is well-known in the context of many-body systems by the Lieb-Robinson bound [29]. There are several well known models which in the continuum correspond to local theories which do not have a Lorentzian causal structure but due to locality the correlations between points decay exponentially with their distance and there is a Lieb-Robinson velocity defined by the bound which information cannot spread faster it.
To be more concrete in the field theory language, there should be a light-cone like structure in the theory which measurements inside the cone should not be affected from outside the cone and vice versa. In the language of our theory in the simple case of d = 1, this means that the simplest thing is to look at the commutator of fields at different space-time points. In the Lorentzian case this is going to vanish for any spacelike separated points. In the case of Lifshitz scalar theory we present a numerical study of this quantity for different values of dynamical exponent and mass parameter in figure 10. It can be shown from the plots that there exists an effective widen cone that the commutator vanishes outside it. This is in agreement from what we expect from Lieb-Robinson bound in the lattice version and also with what we have found from the behaviour of entanglement entropy. We postpone a more concrete study of causality in these theories to [32]. 13

Conclusions and Discussions
We have studied the relaxation process of quenched states in free Lifshitz scalar theories to generalized Gibbs ensemble. Our study was mainly in 2d, which is the simplest case to utilize the quasi-particle picture to understand the process physically. A wide range of dynamical exponents is studied although we have presented results for a few number of them which was enough to focus on the physical picture. An important feature of these theories is a momentum-dependent group velocity of its propagating modes. We have utilized a discrete version of these theories, i.e., Lifshitz harmonic lattice models introduced in [11]. Different regimes of growth of entanglement entropy in these theories after a sudden quench is studied. Our results are mainly captured by an improved version of the quasiparticle picture introduced in [2]. In the following we would like to summarize our main results: • Comparing two specific values of dynamical exponent, say z 1 and z 2 which z 2 > z 1 , as expected from the enhancement of equal-time correlation functions we have shown that the value of EE is larger for z 2 . The growth rate of EE is also grater for the z 2 case, except in a short period of time after z 2 has entered the tortoise saturation regime while z 1 is still linearly increasing. This short period is expected due to two different phenomena. On one hand as the dynamical exponent increases the linear regime is shortened, and on the other hand the occupation number of slow modes is increased. Thus one would expect that the growth rate of z 1 is larger than z 2 during the end of the late linear regime of z 1 .
• We have shown that the larger the dynamical exponent is, the slower EE saturates. Actually exact saturation for these theories (and for any theory admitting a pile of slow modes) is postponed to infinity. For larger values of dynamical exponent, this process becomes slower and slower. This effect is due to the increase of occupation number of slow modes with the dynamical exponent. Based on this, even in the scale invariant case, i.e., the mass parameter of the post-quench Hamiltonian vanishes, no sudden saturation happens in these theories in contrast with scale invariant 2d relativistic theories.
• We have shown that except the cases which the contribution of tortoise modes becomes dominant (in other words the occupation number of these modes is much larger than that for fast modes), the Alba-Calabrese quasi-particle picture works perfectly in these theories. We have shown this for EE and its time derivative for mass quenches. Thanks to Alba-Calabrese, we have an analytic expression for propagation of entanglement in these theories.
As a byproduct we have worked out the analytical time dependence of EE during the tortoise saturation regime for z = 1 and z = 2.
• Our result which was explained in sections 3.4 and 3.5 are showing that the propagation of entanglement is well understood with the existence of an effective notion of widen light-cone in these theories. The larger the dynamical exponent, the wider the light-cone is. Our study is a strong, although still indirect, representative of existence of at least an effective causal structure in Lifshitz theories. Of course from the dispersion relation of these theories this is not obvious at all. We postpone a rigorous study of causal structure in these theories to a later work [32].
• We would like to compare our results with holographic studies of evolution of EE in Lifshitz space times [40,41]. This is mainly studied by considering a Vaidya geometry with an asymptotically Lifshitz spacetime leading to a finite saturation time. But we have shown that based on the existence of tortoise modes the tortoise saturation regime is prolonged and due to that the saturation time goes to infinity. This is another sign for non-robustness of considering asymptotically Lifshitz geometries in a relativistic theory as a dual to a state in a Lifshitz-type theory. 14 There are some other results/comments regarding to our study which we would like to mention in the following: We think that a interesting question which becomes more significant after this work is how is it possible to generalize the quasi-particle picture to be able to capture the strange effects of the tortoise modes in the dynamics?
In this work we have considered a sudden quench since it is the simplest framework to study propagation of entanglement turning off other interesting features of the theory and only focusing on propagation of entanglement. The more general question is what are the critical exponents of these theories under a Kibble-Zurek phase transition which our study is a fast quench limit of that more general setup (for similar studies in Lorentzian theories see [43][44][45][46]). We postpone reporting results regarding to this question to a feature work [47].
The numerical analysis of the early time rapid growth in these theories show that the growth of EE in this regime seems to be have a very good fit with S ∼ t 1+ 1 z . This behaviour was previously found in the holographic context [40,41]. This kind of scaling with t seems a bit confusing. The reason is that one would expect EE in the very early region after a sudden quench, which even no local equilibrium state is reached, to scale with E · A · t α , where E the energy density of the system, A is the area of the entangling region and α is fixed by dimensional analysis which in a Lifshitz theory terns out to be 1 + z. This behaviour is clearly expected to be the case for any number of spatial dimensions. It would be interesting to figure out what is the correct scaling of entanglement in this regime.
Another interesting future direction of this work would be considering a new family of quenches regarding to the other parameter in the Hamiltonian. The dynamical exponent is a parameter in the dispersion relation very similar to the mass parameter. It would be interesting to study the pattern of EE after a dynamical exponent quench specifically from a scale invariant theory (m 0 = 0 and z = 1) to another scale invariant theory (m 0 = 0 and z > 1).

A Time Dependent Correlator Method
We are now interested in a case that the frequency of the Hamiltonian is suddenly changed from its initial value ω 0,k to ω k . After this sudden change the vacuum state of the initial Hamiltonian evolves unitarily with the new Hamiltonian. If we denote the vacuum state of the initial state as |0 we need to compute the following correlators in order to study the time evolution of entanglement and Renyi entropies.
The explicit form of these correlators is given by where X k = 1 ω k ω k ω 0,k cos 2 ω k t + ω 0,k ω k sin 2 ω k t P k = ω k ω k ω 0,k sin 2 ω k t + ω 0,k ω k cos 2 ω k t It is worth to note that in the numerical calculations of this paper we have used a continuous version of these correlators, that is we have used the integral version of Eq.(A.2), that is N x 1 → ∞. Having these we are almost equipped to compute the entropy via the eigenvalues of iJ · Γ which we denote by {ν k (t)} where The entropies are given by where n A is the number sites in region A.