Quench dynamics in lattices above one dimension: the free fermionic case

We begin a systematic investigation of quench dynamics in higher-dimensional lattice systems considering the case of non-interacting fermions with conserved particle number. We prepare the system in a translational-invariant non-equilibrium initial state -- the simplest example being a classical configuration with fermions at fixed positions on the lattice -- and let it to evolve in time. We characterise the system's dynamics by measuring the entanglement between a finite connected region and its complement. We observe the transmutation of entanglement entropy into thermodynamic entropy and investigate how this process depends on the shape and orientation of the region with respect to the underlying lattice. Interestingly, we find that irregular regions display a distinctive multi-slope entanglement growth, while the dependence on the orientation angle is generically fairly weak. This is particularly true for regions with a large (discrete) rotational symmetry group. The main tool of our analysis is the celebrated quasiparticle picture of Calabrese and Cardy, which we generalise to describe the case at hand. Specifically, we show that for generic initial configurations (even when restricting to classical ones) one has to allow for the production of multiplets involving ${n>2}$ quasiparticles and carrying non-diagonal correlations. We obtain quantitatively accurate predictions -- tested against exact numerics -- and propose an efficient Monte Carlo-based scheme to evaluate them for arbitrary connected regions of generic higher dimensional lattices.


I. INTRODUCTION
Finding an efficient description of quantum matter out of equilibrium is every bit as important and timely as it is difficult.Despite the first efforts to attack this problem dating back to the work of John von Neumann in the late 1920s [1], almost one century later we are still lacking general tools to characterise quantum many-body dynamics in an effective and systematic fashion.
These remarkable tools allowed us to understand several consequential concepts pertaining to the dynamics and eventual relaxation of quantum matter out of equilibrium.One key realisation has been that the equilibration process in quantum many-body systems works locally in space, i.e., local subsystems are eventually described by time-independent statistical ensembles even though the dynamics of the full system conserves probabilities [4].Another breakthrough has been to identify the quantum entanglement between a local subsystem and the rest as the "universal observable" able to characterise the full relaxation process in an elegant and basis-independent manner [18,[32][33][34].In essence, one can describe relaxation as the process of turning entanglement entropy into thermodynamic entropy [34][35][36].This process displays an astonishing universality across a huge spectrum of different locally interacting systems which has been explained as the result of a duality between space and time [37] (see also Refs.[38,39]).
Having now sharpened our theoretical tools it is natural to wonder whether we can move on from the onedimensional setting and start exploring higher dimensional cases.In this work we initiate this venture by studying the entanglement dynamics in a d > 1 dimensional lattice system of non-interacting fermions with conserved particle number (tight-binding model), which is driven out of equilibrium by means of a global quantum quench protocol.The main tool of our analysis is the quasiparticle picture of Ref. [17], which is based on the assumption that, after the quench, quantum correlations are transported throughout the system by pairs of correlated quasiparticles created by the quench.Supplemented with a few bits of microscopic data [40], this picture gives asymptotically exact predictions for the entanglement dynamics of free [17,41,42] and interacting integrable [40,43,44] theories, where quasiparticle excitations are infinitely stable.Our work parallels similar studies carried out in the context of continuum quantum field theory [45][46][47]; see also Refs.[48][49][50][51][52] for other field-theory studies of higher dimensional quenches.
We show that, to describe the dynamics from translational invariant initial configurations where the unit cell contains |ν| = ν 1 . . .ν d sites, the quasiparticle picture has to be generalised in the spirit of Refs.[53][54][55].Namely, one has to admit that, instead of pairs, the cor-related quasiparticles form an n-plet [53].Moreover, one has to account for the fact that the quasiparticles generically show complicated, off-diagonal correlations that can be determined by computing the "particle entanglement" [56] for a given bipartition of the multiplet [54,55].Proceeding in this way we obtain quantitatively accurate predictions, which show how the entanglement entropy of a region is transformed into thermodynamic entropy by the time evolution.Then we discuss how this process depends on shape and orientation of the region.We find that regions that are more irregular -characterised by different length scales -display a distinctive multi-slope entanglement growth.On the other hand, we see that the dependence on the orientation becomes increasingly weaker as we increase the discrete-rotation symmetry group of the region.We test our predictions against exact numerics, and propose an efficient Monte Carlo-based scheme to compute entanglement dynamics for arbitrary connected regions of generic d-dimensional lattices.
This paper is organised as follows.In Section II, we introduce the precise setting considered in this work.In Section III we introduce the quasiparticle description of entanglement growth and test its predictions against exact numerics for various d = 1 and d = 2 states.In Section IV we use the quasiparticle description to investigate the entanglement growth in d ≥ 1 and, in particular, how the latter depends on shape and orientation of the subsystem.Finally, in Section V we present our conclusions and outlook.Additional technical details are included in the appendix.

II. SETUP
In this paper we study a global quantum quench protocol in which a many-body quantum system is prepared in a non-equilibrium initial state |Ψ⟩ and let to evolve according to its own unitary dynamics.In this section we describe the specific system/initial state considered and define the observable of interest.

A. Hamiltonian
We consider a system of spinless fermions arranged, for convenience, on a square lattice in d spatial dimensions and linear size L. The dynamics are generated by the following Hamiltonian where J is the coupling strength, n ∈ aZ d L denotes a point on the d dimensional lattice with spacing a, ⟨n, m⟩ indicates that the sum is restricted to nearest neighbours, and finally c n denote canonical fermionic operators.In the following, the lattice spacing a will be set to 1 unless explicitly stated.
The Hamiltonian ( 1) is invariant under one-site translations and is diagonalised by Fourier transform where {c † k , ck } are the Fourier transformed fermions, k ∈ (2π/L)Z d L is a quasi-momentum in the d-dimensional Brillouin zone and k i denotes its i-th component.

B. Initial State
We focus on initial states that are Gaussian, lowentangled, invariant under ν-site translations, and with fixed particle number.Namely we consider states of the form where |ψ ν,j ⟩ is written in terms of fermionic operators within the j-th unit cell.In the following, L = (L, . . ., L) and is a d-dimensional vector, and the operations among d dimensional vectors are always intended element-wise L/ν = (L/ν 1 , . . ., L/ν d ), and, for a given vector n ∈ N d we set This paper will devote particular attention to the subset of these initial states, denoted by |ψ c ν ⟩, where the fermions are at fixed initial positions, i.e., they can be thought of as classical configurations.For these states we have Two concrete examples of these states, one in d = 1 and one in d = 2, are A diagrammatic representation of these classical configurations is provided in Fig. 1.  4) and (6).On the accompanying diagrams, both the lattice spacing a, and an arbitrary subsystem of dimension l, are indicated for later reference.a) the d = 1 state in Eq. ( 7); b) The d = 2 state in Eq. (8).

C. Observable of Interest
We characterise the evolution of the system by studying the dynamics of quantum entanglement between a chosen subsystem A and the rest of the system Ā.Since the state of the entire system is pure, the entanglement is conveniently measured by computing the entanglement entropy, i.e., the von Neumann entropy of the reduced density matrix ρ A for the subsystem A [32].Namely, we consider In a free-fermionic system evolving from a Gaussian state all correlations are encoded in the fermionic two-point functions.In particular, in our case the entanglement entropy is expressed as [57] S(ρ (10) where C A is the correlation matrix of the subsystem A. The latter is obtained from the full correlation matrix by eliminating rows and columns indexing particles in Ā.As an example, in Appendix A we report the explicit form of C n,m for the states (6).In a translational invariant, non-interacting system evolving from a Gaussian state the correlation matrix at time t can be computed directly in the thermodynamic limit ( Ā → ∞) with an amount of resources scaling polynomially with the number of sites of the subsystem A. Therefore, Eq. ( 10) provides an efficient tool to characterise the entanglement dynamics.However, it does not provide direct insight into the relaxation process.To achieve the latter, in the next subsection we present a simple emergent description of the dynamics based on the propagation of stable quasiparticles [17].We will show that, upon supplementing it with a small set of microscopic data, this quasiparticle picture provides an exact asymptotic description.

III. QUASIPARTICLE PICTURE
The quasiparticle picture is based on the observation that non-interacting systems (but also interacting integrable ones [58,59]) feature stable quasiparticle excitations -in our case these are simply the momentum modes in Eq. ( 2).Following Ref. [17] one can then imagine that at t = 0 the quench produces a finite density of quasiparticle excitations, which, upon spreading through the system for t > 0, drive the relaxation process and the growth of entanglement.To elevate this idea to a quantitative description one needs to characterise the motion of the quasiparticles and how they "carry" quantum correlations through the system.This will be our task for the rest of this section.The final result is reported in Eqs. ( 21) and (23), while in Sec.III A we test it against exact numerical results and verify that its asymptotic value coincides with the thermodynamic entropy.
The motion of quasiparticles can be characterised straightforwardly.Since the system under examination is non-interacting, over large scales quasiparticles move like free classical particles and their trajectory is fully specified by their velocities [60].In our case the latter can be directly obtained by computing the group velocity of the momentum modes, i.e., Understanding how these modes contribute to the growth of entanglement, however, is far less straightforward and requires further physical insight.
The key assumption of the quasiparticle picture is that, while moving, the modes generate entanglement in position space but not in momentum space: that is, they merely propagate correlations already present in the initial state [17].Specifically, one assumes that the modes created at the same position are correlated as specified by the initial state and, while moving far apart, they spread this correlation through the system.This means that the entanglement between two regions can be obtained by finding all the multiplets of correlated modes shared between the two regions and summing up their contributions to the entanglement.Therefore, to find a quantitative prediction, one has to determine these quantities [53].
The task is particularly simple when the correlated modes come in pairs.Indeed, Ref. [40] showed that, in this case, the relevant contribution can be inferred from the stationary value of the entanglement entropy.Crucially, however, in our higher dimensional setting the initial states generically create correlations among more than 2 modes.This can be seen by expressing the states in terms of the fermionic operators in momentum space, For instance, considering the classical configurations ∆ (6) in Appendix B we find where |ν| = ν 1 . . .ν d denotes the volume of the unit cell.
As one can infer from this equation, these states generate correlations among |ν| modes.Requiring |ν| to be equal to 2 for d ≥ 2 forces all ν i but one to be equal to one.Namely, one is reduced to consider a state that is effectively one-dimensional.The same conclusion holds for all the states (4).Indeed, as shown in Appendix C, they can all be expressed as in (13) where | ψν,p ⟩ are written in terms of the fermionic operators with quasi momenta p + k and k ∈ (2π/ν)Z ν .Since, for these larger unit cells, one can no longer infer the entropy contributions from the stationary state as per Ref. [40], we instead follow Ref. [54] (see also [55]) and reconstruct the evolution of the initial correlations using a semiclassical approach.This produces a closed form expression for the entanglement entropy valid at large scales and gives a natural way to compute the entanglement contributions.
We begin by subdividing the system in hypercubic cells of linear size ∆, which is much larger than the lattice spacing (one for us) but much smaller than the linear size of A, i.e.
For convenience we take the hypercubic cells to be composed by an integer number of unit cells, i.e., ∆/ν ∈ N d .We then perform a partial Fourier transform to define the modes of the cell Intuitively, the idea is to define modes that are localised in both real and momentum space so that they can have well defined trajectories: see the schematic representation in Fig. 2. Writing the initial state in terms of the modes (16) we have where | ψν,(x,p) ⟩ coincides with | ψν,p ⟩ in Eq. ( 13) if one replaces cp+k with ĉx,p+k .For instance, for the classical configuration states (6) we have Assuming that the cell modes move classically in the limit (15), the reduced density matrix of the subsystem A can be computed by tracing out all the modes that are not in A at time t, namely where we introduced In words, Eq. ( 19) evaluates the reduced density matrix by tracing over the fermionic modes (j, p) that are out of the subsystem A.
Plugging the expression (19) into the definition (9) we have Eq. ( 21) represents the desired quasi-particle expression for the entanglement entropy at time t.
As we can see from Eq. ( 20), the contribution of the correlated multiplet represented by the momentum p ∈ [0, 2π/ν] is found by computing the entanglement between the modes in A and those out of it in the state | ψν,(x,p) ⟩.This quantity is a measure of entanglement between modes, or particles, rather than between regions of space and is referred to as particle entanglement [56].In fact, since the state | ψν,(x,p) ⟩ is Gaussian, the entanglement can be computed using the fermionic correlation matrix as described in Sec.II C. In particular, we define the |ν| × |ν| correlation matrix Ĉ(p, x) with matrix elements given by where k, k ′ ∈ (2π/ν)Z ν .We then define the submatrix ĈA (p, x, t) corresponding to the modes that are in A at time t by eliminating rows and columns corresponding to modes outside of A at time t, i.e., [ Ĉ(p, In terms of this submatrix we can finally write It is worth emphasising that, since the dynamics of modes depends only on the unit cell size ν, the exact structure of the initial state enters the entropy dynamics only through these entropy contributions. The integral in Eq. ( 21) is conveniently evaluated by tracing the motion of quasiparticles.Namely, for fixed (p, t), one can trace the backward light cone of each mode p + k being inside the subsystem at time t as a displacement of the subsystem by −v p+k t.The overlapping regions of these light cones indicate the origin of multiplets that have multiple modes inside the subsystem at time t, i.e. those corresponding to (x, p + k) / ∈ D A (p, x, t), see Fig. 4 for explicit examples.Proceeding in this way we can single out all possible splittings of correlated cell modes contributing to the entanglement -we call them cell mode bipartitions -and the spatial regions where they are produced.The entanglement contribution of each splitting is then evaluated via Eq.( 23).In summary, we can rewrite Eq. ( 21) as follows where B is the set of all cell mode bipartitions, s a (p) is the entanglement contribution of the bipartition a ∈ B, and A a (A, p, t) ≥ 0 is the area of the spatial region producing multiplets contributing to the bipartition a ∈ B at time t and for momentum p.If a specific bipartition a ∈ B does not appear for a given choice of (p, t) we have A a (A, p, t) = 0.
To gain a clearer understanding of how Eq. ( 24) is put into practice, we refer the reader to Section IV A, which considers the simplest irreducible form of these areas in d = 2 arising from a hypercubic lattice and subsystem and the choice ν x , ν y = 2.In this section, helpful illustrations are given for the areas A a (A, p, t) in both this and the analogous d = 1 case, as well as the explicit solution to one of these areas; the complete set of solutions is deferred to Appendix E.
For simple shapes, the areas A a (A, p, t) of the various regions can be determined analytically as a function of (p, t).Computing the corresponding entanglement contributions and integrating Eq. ( 24) numerically over p one can determine the time evolution of S A (t): Section IV A gives a practical example.Instead, for more general regions we evaluate both the areas A a (A, p, t) and the integral over p by means of a convenient Monte Carlo scheme that we detail in Appendix D.

A. Test of the Quasiparticle Formula
The quasiparticle formula (21, 23) may be tested against exact numerical results for finite subsystems A, whereby the real space correlation matrix of the subsystem is obtained in the thermodynamic limit and used in Eq. (10) -a representative example of the comparison is provided in Fig. 3.As shown in the figure, the finite size numerics approach the quasiparticle prediction in the scaling limit This is in agreement with the expectation that the quasiparticle description becomes asymptotically exact in this limit.
Another important test of the quasiparticle solution concerns the long-time behaviour.Indeed, the entanglement entropy is known to approach thermodynamic entropy as time increases [34][35][36].This means that the infinite time value of Eqs.(21,23) should coincide with the thermodynamic entropy of the stationary state reached by the subsystem A. In our case the latter is given by where the integral is over the full Brillouin zone and n(p) is the occupation number of the conserved momentum mode p, i.e., Eq. ( 26) is recovered by Eqs.(21,23) by noting that for t = ∞ the only cell mode bipartitions contributing are those where a single mode of the multiplet is in the system and all the others are outside.This is because, as the modes have different velocities, those starting at the same position are infinitely far from each other at t = ∞ and only one of them can be in A. In this limit the matrix ĈA (p, x, t) becomes 1 × 1 and coincides with the occupation number of the only mode in the system.Namely, where χ A (x) is the characteristic function of A, i.e., χ A (x ∈ A) = 1 and χ A (x / ∈ A) = 0. Plugging back into (21, 23) we then have where in the second step we combined the |ν| integrals over reduced Brillouin zones into a single integral over the whole zone.

IV. RESULTS
In this section we present the predictions for the growth of the entanglement entropy between a region A and the rest of the system after a quench from simple translation-invariant states in d = 1 and d = 2. First we consider the case of A being a simple hyper cubical region aligned with the lattice axes.Then, we investigate the effect of rotations with respect to the lattice.Finally, we study the entanglement growth for more general, irregular-shaped regions.
Note that the last two states of each row can be written with a smaller unit cell, ν = 2 and (ν x , ν y ) = (2, 1) re-spectively.This means that their entanglement dynamics can also be described by pairs of quasiparticles, which is not the case for the other four.In all these cases, the entanglement entropy can be efficiently computed by tracing the motion of the quasiparticles as described in Section III (cf.Eq. ( 24)).An explicit example of this is shown in Fig. 4. In particular, for all the states (30), the relevant cell bipartitions produce three distinct entanglement contributions, {s a (p)} a=1,2,3 .This effectively means that we have to specify only three areas {A a (A, p, t)} a=1,2,3 : see Figs. 4(a) and 4(b) for illustrations of the three areas in d = 1 and d = 2 respectively, where in the latter we consider the more general case of l x ̸ = l y to anticipate the section that follows.An example of the explicit form for one of these areas, A 2 (A, p, t) from Fig. 4

(b), is given by
where The explicit form for all areas A a (A, p, t) is reported in Appendix E, while the entanglement contributions are reported in Tables I and II.The resulting quasiparticle predictions are compared to the exact numerical solutions, in Fig. 3.The left and right panels correspond respectively to d = 1 and d = 2.   24) and Fig. 4b), for the three independent ν = 4 states of Eq. ( 30).
Apart from the agreement between quasiparticle solution and exact numerics for increasing system sizes, which we have already stressed in Sec.III A, these figures show two significant features.First, we see that the d = 2 case is clearly distinguished by the occurrence of a nonlinear initial regime (see the inset of Fig. 3).Indeed, using the explicit form of A a (A, p, t) we see that Eq. ( 24) contains a quadratic term in time proportional to (2s 1 − s 2 − s 3 ).From the entropy contributions presented in Table II, we see that the magnitude of this term is largest for , while it is zero for .This is consistent with the fact that the latter is effectively a one-dimensional setting.
Another key takeaway from Fig. 3 is the difference between the entropy plots of states with equal occupation numbers (27) (the last two states in both lines of Eq. ( 30) have n(p) = 1/2).Examples include the slower initial growth of versus in Fig. 3(a), and of versus in Fig. 3(b), which share the same saturation value.One can offer a heuristic explanation of this slower growth: states of equal charge density but less uniform site occupation impose more constraints on the initial site hoppings; the plots indicate that these constraints last until correlated quasiparticles first span the full width of the subsystem.Importantly, since occupation numbers fully specify the expectation value of all conserved charges, our result shows that the expectation values of all conserved charges are not enough to determine the full-time entanglement dynamics even in the scaling regime (25).
A stark feature of Tables I and II is that the classical configurations (30) turn out to give momentum independent entanglement contributions.To make sure that this property does not introduce any qualitative difference in the entanglement dynamics we also consider more general, non-classical states described by Eq. 4, where |ψ ν,j ⟩ is now a superposition of states in the |ν|-site unit cell that satisfies the Gaussianity condition given by Eq.F11.The entanglement contributions for these superposition states are reported in Appendix F while Fig. 5 reports some representative examples.Although in this case the entanglement contributions become momentum dependent (cf.App.F), we see that the plots are qualitatively similar to those in Fig. 3 and the agreement with the exact numerical solution is still excellent for large enough subsystems.

B. Rotations with respect to the lattice
Let us now consider the dependence of the entanglement growth on the orientation of A with respect to the underlying lattice.We begin considering the simple case of a rectangular region in d = 2 that is rotated by an angle θ with respect to the lattice.
In this case the explicit calculation of the areas in Eq. ( 24) becomes quite tedious and the quasiparticle prediction is more conveniently obtained integrating Eq. ( 21) via the Monte Carlo scheme discussed in App.D, which agrees with the explicit approach.The results for a representative initial state are reported in Fig. 6.From Fig. 6(a) we see that when the aspect ratio of the rect- angle r ≫ 1 the entanglement dynamics depends quite markedly on the orientation of the region, with the rectangle aligned with the lattice showing a slower relaxation.Interestingly, however, we see that the dependence on rotation angle decreases smoothly with the aspect ratio r such that whenever the edges l x and l y coincide the dependence on the rotation angle almost disappears.To exclude that this is not an artefact of the Monte Carlo integration routine we reproduced the result for the rotated square region using the quasiparticle tracing integration of Eq. ( 24) finding exact agreement.In this case the relevant areas depend on the rotation angle, see Fig. 10, and their explicit expression is reported in Appendix E while the entanglement contribution associated with each area is reported in Tabs.I and II of the previous section.
The behaviour in Fig. 6 can be expected as the square is "more rotationally symmetric" than the rectangle.More precisely, it is left invariant by greater number of discrete rotations: its cyclic group is of order 4 rather than 2. To highlight how the order of the cyclic group of A affects the orientation dependence, in Fig. 7 we consider polygons with cyclic group of order q = 1, 3, 5.We see that, as expected, the dependence on θ decreases with q.A surprising aspect, however, is how quickly it does so: the θ dependence is already negligible for q = 3.

C. General Shapes
Finally, we use our quasiparticle approach to investigate the the entanglement growth of irregular, connected, regions characterised by different cross sections.Interestingly, we observe that such regions display an entanglement dynamics that is qualitatively different from that of regular ones.
For instance, in Fig. 8 we show the entanglement dynamics of a region in the shape of five point star (see inset panel of Fig. 8 for an illustration) for different orientation angles θ with respect to the underlying lattice.We see that, as expected, there is essentially no dependence on θ, however, the entanglement evolution reported in the figure is quite peculiar: rather than the usual linear increase followed by saturation the entanglement shows a complicated multi-slope curve.This is due to the fact that each cross section of the figure corresponds to a non-analyticity of the quasiparticle prediction in time [61].These special moments correspond to the points in time where the backward light cones associated to the fastest quasiparticles separate through a cross section -in other words, when a sudden gap appears between the light cones.For instance, for a rectangular region and the states in the second line of Eq. ( 30), the quasiparticle solution obtained integrating Eq. ( 24) with the A a (A, p, t) in Appendix E reads 0 0.2 0.4 0.6 0. as where we introduced ζ = 2v max t/ l x l y , r = l x /l y , and This function has a non-analytic points (corresponding to a discontinuous second derivative) for z = 1: this means that (32) has non-analyticities at ζ = √ r and ζ = 1/ √ r.
An instance in which these two points are visible is the mid-blue solid curve of Fig. 6.More irregular regions have many of such non-analytic points and originate peculiar looking curves like the one in Fig. 8.

V. CONCLUSIONS
In this paper, we studied the spreading of entanglement in a free fermionic system defined on a lattice of dimension d ≥ 1 by generalising the quasiparticle picture of Calabrese and Cardy [17].In particular, we showed that if the initial state has a fixed number of particles, and is invariant under no less than ν j discrete lattice shifts in the direction j = 1, . . ., d, the quench produces a multiplet of ν 1 • • • ν d correlated quasiparticles.This means that only settings that are effectively one-dimensional can produce pairs: in all non-degenerate cases one has to consider larger multiplets.
Characterising the spreading of entanglement by generic multiplets of quasiparticles, we derived a general integral formula, Eq. ( 21), for the evolution of the entanglement entropy.Then, we studied its explicit predictions for d = 1, 2 in the case of a square lattice.In particular, we introduced an efficient Monte Carlo scheme for d ≥ 1 to study the entanglement of arbitrary connected regions in d = 2.
First, we showed that exact diagonalisation results recover the generalised quasiparticle description in the limit of large subsystems and times, i.e., when the quasiparticle picture is expected to apply.Then, we studied how the entanglement dynamics depends on shape and orientation of the subsystem with respect to the underlying lattice.We observed that, for subsystems with cyclic symmetry group of order larger than three, the dependence on the orientation is negligible.Moreover, we showed that irregular regions show a multi-slope entanglement growth.Interestingly, our results also provided simple examples showing that specifying the expectation values of all conserved charges is not enough to fully determine the entanglement dynamics, even at leading order.Namely, we found examples of two different initial states with the same expectation values for all local conserved charges that generate different entanglement dynamics.
Our general quasiparticle formulation and its Monte Carlo implementation provide a flexible and versatile method to study entanglement dynamics in noninteracting systems.The approach can directly be applied in charge conserving fermionic systems with arbitrary dispersion relation -including the case of anisotropic couplings -to study the entanglement of arbitrary regions in hypercubic lattices of generic dimensions d > 1.However, it can also be directly generalised to study systems without charge conservation, e.g.BCSlike ones, and to systems on arbitrary (regular) lattices.In all these cases, whenever the initial state is not onesite shift invariant, the entanglement will generically be transported by n-plets of correlated quasiparticles quasiparticles with n > 2.
A more immediate future direction for our work, however, is for us to use our generalised quasiparticle description to study the restoration of a discrete symmetry broken by the initial state and the possible occurrence of the quantum Mpemba effect [62,63].This can be efficiently done by using the recently introduced entanglement asymmetry [62,64], which can be treated using the quasiparticle picture.
Note Added.While this manuscript was being finalised, we became aware of the related work [65].The latter also studies entanglement dynamics in higher dimensional free fermionic systems (d = 2) but focuses on special regions that can be treated by the technique of dimensional reduction [66] (see also Refs.[67,68]).
This gives the time-dependent correlations e ip•(n−m) e ik•(m+ap) e it(ϵ(p−k)−ϵ(p)) .(A2) In the thermodynamic limit L → ∞ this becomes lim Appendix B: Fourier Transform of Classical Configurations The initial state defined by Eq. ( 6) may be written as The Fourier transform is defined for a square lattice of spatial dimension D and length L d along each spatial dimension.Provided each L d /ν d is an integer, the Fourier transform of c † ν•j−ap may be written as Applying Eq. (B2) to the operator b † ν,p of Eq. (B1) gives the Fourier transform of this operator as where (B † ν,p ) 2 = 0 is used to write as a sum over all possible permutations of n.Using that e ipα(ν d β−ap) is a Van der Monde matrix, we can then write its determinant as [69] det e ipα(ν (B4) Since this difference of elements is nonzero for all α ̸ = β, we see that the determinant of this matrix is also nonzero.The normalisation of B † ν,p then implies that this determinant must be equal to one.We therefore have Applying this result to the initial state of Eq. (B1) leads to the Fourier transform of this state given by Eq. ( 14), | ψν,p ⟩ .

Appendix C: Fourier Transform of General Initial States
Given that the initial state ( 4) is Gaussian, its density matrix has the exponential form In the end, we are left with a sum over all particle entanglement contributions, which we divide by the number of samples taken to get the approximation of the integral Eq. (D1).Note that, by also computing the average of the squares of the particle entanglement contributions, we can also keep track of the variance σ 2 of the Monte Carlo sampling, and so estimate the standard error of the mean σ/ √ number of samples.
Appendix E: Explicit form of Aj for rectangular regions in d = 1, 2 Here, we present the analytic solutions to the functions A j (A, p, t) defined by Eq. 24.Following sections IV A and IV B, we focus on d = 1 and d = 2 initial states, where A is a simple hypercubic subsystem.
We begin with the d = 1 states of IV A, namely with ν = 4 states in a subsystem of length l, whose three areas {A j (A, p, t)} j=1,2,3 are specified by Fig. 9(a).First, we write the solutions to A j (A, p, t) for arbitrary velocity ordering v a > v b > v c > v d .For convenience, we define the natural times τ ab = l/(v a − v b ) and lengths ∆ ab = (v a − v b )t, where we have applied this symmetry of mode velocities v c = −v b and v d = −v a , to write where H(x) is the Heaviside step function.These solutions can be rearranged into a full timewise solution of the quasiparticle dynamics as shown in Fig. E.

Time Interval Total Entropy
where the contributions {s j } j=1,2,3 for the d = 1 classical configurations of Section IV A are defined in Table I.Next, for these classical configurations, we introduce the integrated areas, We then combine Eqs.(E1) and (E3) to obtain the solution to the integrated areas {A j (A, t)} j=1,2,3 .It should be noted that, while these integrated areas yield the solution to all states of fixed s j (p) = s j , Eq. (E2) can easily be modified to a state dependent solution of general s j (p).
We now turn to the d = 2 states of IV A, namely to ν x = ν y = 2 states in a rectangular subsystem of lengths l x , l y aligned with the lattice, whose three areas {A j (A, p, t)} j=1,2,3 are specified by Fig. 9 (b).Our method closely follows the d = 1 case, where first we write the solutions to A j (A, p, t) for arbitrary velocity ordering v a > v b > v c > v d , ordered by projection onto the positive x-axis.For convenience, we define the natural times τ i ≡ l i /4 sin(k i ) and lengths X, Y ≡ 4 sin(k i )t for i = x, y, where we have applied this symmetry of mode velocities v c = −v a and v d = −v b , to write A 1 (A, p, t) = 4XY H(min(τ x , τ y ) − t) + 4l x Y H(t − τ x )H(τ y − t) + 4Xl y H(t − τ y )H(τ x − t) + 4l x l y (H max(τ x , τ y , t)) A 2 (A, p, t) = 2(l x − X)Y H(min(τ x , τ y ) − t) + 2(l x − X)l y H(t − τ y )H(τ x − t) A 3 (A, p, t) = 2(l y − Y )XH(min(τ x , τ y ) − t) + 2(l y − Y )l x H(t − τ x )H(τ y − t) .
In d = 1 we consider a superposition state of the form where we omitted an overall constant ensuring normalisation.The contribution to the entanglement when only one mode of the multiplet is in or out of the system depends on the mode.Specifically we have where and p = 0, • • • , 3 identifies the mode.Instead, when only two modes are in the system we have s p,q (k) = − λ p,q (k) log λ p,q (k) − (1 − λ p,q (k)) log(1 − λ p,q (k)), (F4) with The contribution of a single mode in (or out of) the system turns out to be mode-independent and is given by (F7) Concerning the contributions of two modes, the only ones produced by the quasiparticle dynamics are where g(k) = α cos(k)/(1+α 2 ).The saturation value for entanglement density is given by 4(α 2 +1) ), (F9) for both the d = 1 and d = 2 cases.In the limit α → 0, ∞ these contributions recover those presented in the second two columns of Tabs.I and II.We remark that the initial state should be Gaussian in order for our techniques to apply.This property is not immediately obvious for superposition states; however, in our translational invariant setting we can verify it by restricting to the unit cell.Specifically, we require a general superposition state such as to be annihilated by a new family of canonical fermions f i that is linearly related to c i and c † i .Namely, one can write This is always the case for the superpositions considered here.

FIG. 1 .
FIG.1.Examples of initial states defined by Eqs.(4) and(6).On the accompanying diagrams, both the lattice spacing a, and an arbitrary subsystem of dimension l, are indicated for later reference.a) the d = 1 state in Eq. (7); b) The d = 2 state in Eq. (8).

FIG. 3 .
FIG. 3. Plots of the entropy over linear dimension SA(t)/|A| against rescaled time for (a) d = 1 states with ν = 4, and (b) d = 2 states with νx = νy = 4 and a square subsystem.The states that are displayed reflect the number of 'independent" states that do not map onto one other by single-site shifts of the unit cell, i.e. the number of states with unique entropy dynamics for a hypercubic subsystem.To compute the integral we employ an inverse FFT with (a) 10000 subdivisions, and (b) 250 × 250 subdivisions.For each state, the quasiparticle solution (QP) is plotted against the finite size numeric solution with two different values of linear subsystem length |A| = l in order to illustrate the rate of convergence.The dashed lines indicate the saturation values obtained from the stationary state solutions.The inset of (b) focuses on the initial regime, comparing each quasiparticle solution against a straight dashed line tangent to the solution at t = 0.

A. d = 1
and d = 2 Square Subsystems We begin our discussion considering the case of A being either a connected segment in d = 1 or a rectangular region aligned with the lattice (with edges in the x and y directions) in d = 2.To fix the ideas we consider cases producing four correlated modes: (i) d = 1 and ν = 4; (ii) d = 2 and ν x = ν y = 2; and look at the following classical configurations

A 1 FIG. 4 .
FIG. 4. Diagram of the quasiparticle dynamics for a) a d = 1 state with ν = 4, and b) a d = 2 state with νx = νy = 2 and rectangular subsystem, where lx > ly.The red dashed line marks the boundary of the subsystem and the three areas {Aj(A, p, t)}j=1,2,3 are depicted for a) 0 < p < π/4 and l/2(sin(p) + cos(p)) < t < l/4 sin(p), and b) px = py and t < lx/4 sin(px).These areas are obtained by tracing the motion of quasiparticles as outlined in Section III, and their explicit solutions are given in Appendix E.
the momentum integral of Eq. (E2) into intervals for which the ordering of mode velocities is fixed, A j (A, t) = A, p, t) .(E3)