Entanglement Growth after a Global Quench in Free Scalar Field Theory

We compute the entanglement and R\'enyi entropy growth after a global quench in various dimensions in free scalar field theory. We study two types of quenches: a boundary state quench and a global mass quench. Both of these quenches are investigated for a strip geometry in 1, 2, and 3 spatial dimensions, and for a spherical geometry in 2 and 3 spatial dimensions. We compare the numerical results for massless free scalars in these geometries with the predictions of the analytical quasiparticle model based on EPR pairs, and find excellent agreement in the limit of large region sizes. At subleading order in the region size, we observe an anomalous logarithmic growth of entanglement coming from the zero mode of the scalar.


Introduction
A global quench is a simple setting in which we can study thermalization in isolated quantum systems: at t = 0 we start with an atypical translationally invariant, short range entangled initial state |ψ 0 , and let the state evolve in time. 1 In a generic quantum system, during the process of thermalization all simple observables converge to the value they take in the Gibbs ensemble. A good characterization of thermalization is how close the reduced density matrix of a small subsystem A, ρ A [|ψ(t) ] is to the reduction of the thermal density matrix to the region A, ρ th A ∝ TrĀe −βH , where β is to be chosen such that the expectation value of the energy agrees between the two density matrices. One way to quantify the proximity to thermal behavior is to calculate the von Neumann entropy of ρ A (t), and follow as it evolves from an area law value to saturation at the thermal entropy.
In a free theory, because of the infinitely many conserved charges the above picture requires modification. The time evolution leads to simple observables converging to their values in the Generalized Gibbs Ensemble (GGE) [1] instead of the Gibbs ensemble. In this paper we will work with Gaussian states in free scalar field theories: for these states it is known that the set of charges that one has to include in the GGE are the particle numbers in each momentum mode [2]. 2 After the quench we focus on the case of massless fields.
We investigate entanglement entropy growth of these fields for geometric subregions in diverse dimensions. We discretize the theory on a lattice and use the correlator matrix approach to numerically compute entanglement and Rényi entropies. The continuum limit can be achieved by taking a scaling limit: where R is the characteristic size of the region A, t is the time measured from the quench, a is the lattice constant, and β k is the inverse of the effective temperature in the mode with wavenumber k. Let us introducê to get rid of the vacuum area law pieces in the entropy. 3 In the limit of large region sizes and times R, t β k , (1.3) it is expected that the entropy obeys a scaling form: where s is the entropy density in the GGE, 4 f (0) = 0 follows from the definition, and f (∞) = 1 assumes that the entropy reaches the equilibrium value predicted by GGE.
In the limits (1.1), (1.3) the finite area law pieces in the entropy are suppressed by the factor R/β. In summary, we want to work in the double scaling limit There is a useful toy model for entanglement growth introduced in [6,7] and generalized to higher dimensions in [8]. This model assumes that the quench creates quasiparticle EPR pairs 5 localized on length scales O(β). In the scaling limit (1.5) the pairs can be taken to be pointlike. In a massless theory, the pairs then propagate with the speed of light, 6 andŜ A (t) counts the number of pairs that have one member in A and the other inĀ. While the model reproduces the entropy of one interval in any 1+1 dimensional conformal field theory (CFT) [6], for more complicated geometries it only works in integrable CFTs [11]. In this work we find overwhelming evidence that the quasiparticle model reproduces the growth of entanglement in higher dimensional free massless scalar field theories in the scaling limit (1.5), by comparing the predictions of the quasiparticle picture to numerical computations in strip and sphere geometries, see Fig. 1. We study two types of quenches: the boundary state quench corresponds to starting the evolution from a regularized boundary state of the CFT, which leads to β k = β, while the mass quench corresponds to abruptly changing the Hamiltonian of the system by changing the mass parameter, and leads to a k-dependent effective temperature. The quasiparticle picture works for both quenches equally well.
We emphasize three key features of our findings. First, we find (in the two geometries we consider) that at early times the entropy exhibits linear growth of the form:Ŝ 4 In a generic system without any conserved quantity other than the energy, it would be the thermal entropy density. 5 In higher dimensions we can consider more complicated patterns of entanglement, as explored in [8]. Intuitively, however, for Gaussian states that we consider in this paper the bipartite entanglement structure encoded in EPR pairs seems to be the appropriate choice. 6 In massive integrable models, they follow a nontrivial dispersion relation [9,10]. Ā A R Figure 1: The two types of geometries we examine in this work. Regions A andĀ partition the system into two distinct regions. Starting with a pure state, we trace out regionĀ to obtain a reduced density matrix ρ A , from which we compute the entanglement and Rényi entropies. Left: The strip geometry with two sides separated by a distance 2R. Right: A spherical geometry of radius R.
where by area(A) we mean the area of the entangling surface, vol(∂A). The dependence on the shape only appears through area(A), and the entanglement velocity v E is shape independent. 7 Second, we comment on the saturation time t S . For spherical geometries entanglement saturates as fast as allowed by causality [15], (1.7) We find that t S is strongly shape dependent, 8 and for a strip geometry 9 We reiterate that the results (1.6), (1.7), and (1.8) are in agreement with the quasiparticle model. They are just simple properties of the functionŜ A (t) in the limit (1.5), which according to our findings is in complete agreement between the numerical computation in the free massless scalar field theory and the quasiparticle model. Third, we 7 In the regime β k t R the curvature of A should be irrelevant for the process, so (1.6) is intuitive. (1.6) is also known to hold in strongly coupled theories with a holographic dual [12,13,14]. 8 In chaotic (holographic) examples the shape dependence of t S is mild, but still non-trivial [13,14,16]. 9 The intuition behind (1.8) is that there are quasiparticle pairs propagating almost parallel to ∂A that take an arbitrary long time to start to contribute to the entropy. point out an unexpected aspect of our numerical results: we see a logarithmic growth of entropy even after the saturation time (1.7), which however is subleading in the limit (1.5), and therefore does not spoil the agreement with the quasiparticle model in the appropriate regime. 10 We identify the scalar zero mode as the source of this logarithmic growth, but the phenomenon deserves further investigation.
Besides the intrinsic interest in the study of equilibration in free field theory, we are also motivated by the scarcity of computations of entropy growth in field theories. Our results elevate the status of the higher dimensional quasiparticle model from a toy model to an actual description of entanglement growth in free massless scalar field theories. 11,12 The results then provide a useful benchmark for strongly coupled theories: the conclusion of [8] that in strongly interacting (holographic) theories entanglement spreads faster than allowed by free streaming, is reinforced. In general, collecting results on entanglement growth from various systems could lead to further insight into the workings of equilibration in quantum systems, both integrable and chaotic. For further discussion from this viewpoint see [16].
Using similar techniques, it is possible to study global quenches in free fermion theories. The analytical and numerical techniques for analyzing global quenches in free scalar fields could potentially be extended to interacting field theories either perturbatively [17] or non-perturbatively [18,19]. Such generalizations could shed new light on the dynamics of entanglement in interacting systems. In this paper we restrict our attention to instantenous quenches. It would be very interesting to extend our analysis to smooth quenches, where the duration of the quench δt introduces a new time scale. In the limit R, t, δt β, we expect the entropy to again obey a scaling form (1.4), but the scaling would become a function of two variables f (t/R, δt/R). Correlation functions obey universal scaling laws in this limit [20,21,22,23], and it would be interesting to explore, if those results carry over to the case of entanglement entropy. It would also be interesting to see, if a modification of the quasiparticle model could reproduce the scaling function f (t/R, δt/R). Perhaps, smearing the time of origin of the EPR pairs could be a useful starting point [8].
The plan of the paper is as follows. In Sec. 2 we provide an introduction to our setup: we review the correlation matrix approach of computing entropies, and we discuss the quenches considered, along with the quasiparticle model. Sec. 3 contains the numerical results for the entanglement and Rényi entropies, and a comparison with the quasiparticle model gives excellent agreement. A brief investigation into the logarithmically growing mode is also included. Some further details of the setup are relegated to the Appendices.

Time evolution of entanglement
We consider the time evolution of a Gaussian wave function in free scalar field theory through a quench. What enables us to do the computation is that the time evolution of an arbitrary Gaussian initial state remains Gaussian in a free theory. The computation simplifies in the global quench setup due to the preservation of translational and rotational symmetry: the kernel of the Gaussian remains diagonal for all times in momentum space. We can then apply the machinery developed for Gaussian states in free field theories [24,25,26,27,28] that we review below.

Gaussian wave function in free scalar field theory
Let us consider the Hamiltonian for a free massive scalar in d + 1 spacetime dimensions where π is the canonical momentum for φ. The Hamiltonian (2.1) can be discretized and written in a general form We will consider Gaussian wave functions 13 where N (t) includes the normalization of the wave function and an overall time dependent, but φ i -independent phase. If the wave function is of the form (2.3) at one instant in time, it remains of the same form for all times when evolved by (2.2). The ground state of the system is obtained by setting Ω = √ K.
Let us diagonalize K by making the orthogonal transformation O on the fields φ and canonical momenta π then K takes the form where K D is diagonal. K is the discretized version of the operator −∇ 2 + m 2 and O is the discrete Fourier transform. Then it is clear that to describe translationally and rotationally invariant quenches, it suffices to restrict to the case where Ω(t) commutes with K. In terms of the new q variables the Gaussian (2.3) is diagonal: A simple example to keep in mind is a mass quench: we prepare the initial state through an abrupt change of the Hamiltonian H m 2 → H m 2 =0 , which is implemented by the change K 0 → K = K 0 − m 2 1. All of our claims above are readily verified for this case.
As a warm-up problem for the time evolution, let us consider a quench in which we change the frequency of a harmonic oscillator abruptly at t = 0, from frequency ω before the quench toω after the quench. With initial state given by the pre-quench ground state the solution to the time-dependent Schrödinger equation reduces to a complex Riccati equation for the kernel Ω(t), which can easily be solved by standard methods to give the post-quench wave function 14 For the correlation matrix approach, we need to calculate all two-point functions. 15 A straightforward computation gives (2.9) Note that we do not lose any information by considering 1 2 {q, p} instead of qp, because An important physical quantity is how much energy we inject into the system when we quench this harmonic oscillator: (2.10) 14 There are two easy checks of this formula: at t = 0 it gives back the initial Gaussian, and for ω =ω we get the ground state wave function with trivial time dependence. 15 One point functions vanish by construction.
Based on the solution (2.8) for a single harmonic oscillator, we immediately see that for a collection of harmonic oscillators of the discretized scalar field theory, the initial and post-quench wave functions are where ψ(t, q) is given in (2.8). In (2.11) it is understood that one should make the replacement q → q i ,ω →ω i , and ω → ω i , where ω 2 i are the eigenvalues of Ω(0) that characterize the initial state, andω 2 i are the diagonal elements of K D . Recall that as discussed around (2.5), Ω(0) and K can be simultaneously diagonalized. To be completely explicit, we write out the φ-dependent part of the wave function as (2.12)

The correlation matrix approach to quenches
For this wave function it is now easy to determine the two-point functions of the canonical variables. The generalization of the single harmonic oscillator results for the two-point functions (2.9) is Let us introduce a vector of canonical variables in region A (we trace overĀ) where i, j = 1, . . . , n are restricted to region A and I, J = 1, . . . , 2n. The canonical commutation relations read: 15) and the correlators can be collected into a 2n × 2n matrix [24,25,27]: Note that Γ IJ is a real symmetric positive definite matrix. Such matrices can be brought to Williamson normal form, i.e. there exists a symplectic matrix M that diagonalizes them. M implements a canonical transformation (that preserves (2.15)): (2.17) The easiest way of determining γ k is to obtain the eigenvalues of the matrix iJΓ, which are {±γ k }. We have now successfully mapped the problem to computing the entropy of n harmonic oscillators at finite temperatures: and the Rényi entropies (2.20) In the symmetric geometries A that we consider in this paper, the matrix Γ is block diagonal. In the case of the strip, the different blocks are labelled by the momenta parallel to the entangling surface; in the case of a sphere, the labels are the angular momentum quantum numbers. The matrix Γ is block diagonal as two point functions do not mix different linear or angular momenta. In these cases, the above steps can be performed block by block, and the entropy is just the sum of the contribution of each block. Details of the different coordinate systems can be found in the Appendices.

Different types of quenches
As discussed in the introduction, in free theories there are infinitely many conserved charges, and equilibration only happens in the sense of the GGE. What this means for our purposes is that any mode can be quenched independently, and they have their own effective temperature. We shall focus on the case in which after the quench the mass of the scalar field is zero (soω(k) = k), so in the continuum limit (1.1) the time evolution is governed by a CFT.

Boundary state quench
From the point of view of effective temperatures mode by mode, a particularly nice state to consider is the conformal boundary state model of a quench [6] where |Dirichlet is the Dirichlet boundary state. This state will have finite energy in any dimensions, with mode-independent inverse temperature β k = β, and is specified by the relation (2.22) The requirement that the quench is described by the continuum field theory translates into β/a 1, as discussed around (1.5). 16 The Rényi entropy density arising from this quench is (in d spatial dimensions) leading to the entropy density .
(2.24) 16 Rewriting the inequality as 1 a 1 β , we intuitively want the energy scale of a thermal excitation, which is approximately 1/β, to be far less than the highest energy excitations which can be supported by our lattice theory, which go as 1/a.

Mass quench
In contrast to the boundary state quench, the mass quench -although it may seem more physical -has less favorable properties. In particular, the relation produces a mode-dependent temperature [29] i.e., high energy modes have diverging effective temperature. 17 Now the inequalities in (1.5) will not be satisfied for all k. Nevertheless, we can intuit that the weaker condition m a 1 (2.27) should guarantee that we stay close to the continuum limit, which corresponds to low energies. However, additional complications arise for d ≥ 3, where the mass quench does not produce a finite energy state, as the total injected energy (in excess of the vacuum energy) is where we used (2.10). In the continuum limit the change in energy density is where Λ U V m is some UV cutoff scale. For spatial dimensions d = 1, 2 the mass quench produces a state with finite energy density as Λ U V → ∞, while for d ≥ 3 we encounter ultraviolet divergences; in particular in d = 3 we find a logarithmic divergence. We may summarize these results as .
(2.30) 17 The average energy in a high energy mode is still low.
The Rényi entropy density for the mass quench is given by (2.23) with β → β k , leading to where γ E is the Euler-Mascheroni constant and ψ(z) is the digamma function. The above equations yield the entropy density and s = ∞ is divergent for spatial dimensions d ≥ 4. Note that in d = 3 the entropy density is finite, even though the energy density is infinite after a mass quench.
Furthermore, one also anticipates a divergent area law contribution to the entropy in any number of dimensions. The entropy difference from before to after the quench (1.2) is expected to result in an infinite area law correction for d ≥ 3. For example, for d = 3 there should be a log divergence in the change of the area law contribution [30] ∆S area = A m 2 24π log Λ U V m . (2.33) So we will only focus on mass quenches in 1 and 2 spatial dimensions.

Other quenches
The above formalism extends to any quench in which the kernel of Ψ 0 is diagonal in Fourier space, meaning it takes the form (2.11). We have discussed boundary state and mass quenches, since they are perhaps the most physical examples, but an arbitrary choice of ω(k) is allowed. If we want to preserve translation and rotation invariance, we only need to choose a function ω = ω(|k|). This can be parametrized by a type of mode dependent initial mass that we quench for the initial wave function. 18 Choosing a mass function that decays to zero (m 2 (k) → 0) fast enough for large k can make ∆E finite in any dimension.
Another generalization that we may consider is to replace the instantaneous quench with a smooth quench of duration δt, which would render the energy density injected into the system by a mass quench finite [23].

The quasiparticle model for a quench
In this section, we will review the dynamics of entanglement using the quasiparticle picture in which entanglement is carried by a uniform density of noninteracting EPR pairs [6,7,8]. We assume that the two quasiparticles which comprise a pair travel in opposite directions at the speed of light, with an isotropic angular distribution.
In the quasiparticle model one can first fix a point x and time t, and determine the contribution to the entropy of region A coming from EPR pairs that originated from that point. These pairs will be positioned on a sphere of radius t with center x. Let us denote the part of this sphere incident in region A by L A (x, t), and by µ[L A (x, t)] the contribution to the entropy of this region. We have to sum over the points of origin to obtain the entropy of region A: Although not discussed in [8], the computation presented here applies both to the entanglement and Rényi entropies with the appropriate entropy density s (or s q ) used in (2.36).
In [8] the integral (2.35) was evaluated for symmetric entangling regions. Here we will need the results for S A (t) in d = 1, 2, 3 spatial dimensions for a strip of width L = 2R and a sphere of radius R. In d = 1, both cases degenerate to A being a finite interval of length 2R. For the the strip geometry in spatial dimension d, we have . (2.37) Likewise for spherical geometries in spatial dimension d, . (2.38) The saturation times (1.7) and (1.8) can be easily read off from these expressions. The expression for the entanglement velocity in d spatial dimensions is (2.39) In the next section, we will confirm the predictions of S A (t) (and hence for t S and v E ) with numerical simulations of global quenches of free scalar fields.
3 Numerical results for strips and spheres

Intervals in 1 spatial dimension
In d = 1 spatial dimension the results for the entropy for intervals of different sizes in a boundary state and a mass quench can be found in Fig. 2. For convenience we impose periodic boundary conditions at the ends of the 1 dimensional region. In this figure we have used the subtracted entropy (1.2). Mass quenches for the similar systems were analyzed in [31,27]. Related analytical and numerical results for local and global quenches in [26] and references therein. In Fig. 2 we see that the two types of quenches, the boundary quench and the mass quench, give results that are nearly indistinguishable to the eye. The quasiparticle prediction [6,7] matches closely for t < R with the correct entropy density s (2.24). However, instead of sharp saturation at t = R, we see that the entropy keeps increasing as a logarithm of time. To understand this deviation, we reproduce the same quenches in Fig. 2 (bottom), but plotted in logarithmic time. From fitting we find that the coefficient of the logarithm is independent of R, and appears to equal 1/2. 19 Since the coefficient of the logarithm is independent of R, in the limit (1.5) this 1 2 log t behavior is subleading. Hence the prediction of the quasiparticle model is obeyed in the limit of large region size and time. 20 Nevertheless, it is interesting to understand the origin of this 1 2 log t behavior, since naïvely we would have expected saturation in finite time (possibly with 1/t # power law behavior), and corrections to the quasiparticle model to be suppressed by β/R. The massless free scalar theory is known to exhibit peculiar logarithmic corrections in the entanglement entropy in static situations due to the presence of a zero mode. 21 This motivated us to modify our setup in an attempt at getting rid of the contribution of the soft modes. We observe that the 1 2 log t growth disappears if the quench leaves the soft modes in their ground state (i.e. we choose an appropriate m 2 (k) in (2.34)), 22 or if we take the harmonic chain to be finite and consider the interval to be at the end of the chain. Both cases are analyzed in Fig. 3 and Fig. 4, and the logarithmic growth is clearly gone.
Based on these results, we can give a heuristic explanation of the 1 2 log t behavior based on the dynamics of the zero mode. The following argument was suggested to us by A. Wall. After all the other modes have saturated, we can concentrate on the noncompact zero mode of the scalar field. Its wave function is initially localized, and it spreads under time evolution. The width of the wave function should go as √ t. Regarding the entropy as the number of available states we immediately obtain the contribution log # √ t = 1 2 log t + . . . to the entropy. That the zero mode contributes the logarithm of its target space volume to the entropy was discussed before in [35,36]. In the smooth mass quench analyzed in Fig. 3, the zero mode is not excited, while 19 In more detail, we have fittedŜ for the data points with t > R. (3.1) diverges for t = R, so we started the fitting procedure a couple of lattice units later in time. One can also consider introducing a time shift as an additional fitting parameter, but this hardly changes the value of S and c. We found that c = 1/2 within 2% accuracy. On Fig. 2 (bottom) we have plotted (3.1) with c = 1/2 and S fitted. The match is excellent. 20 Unless we extrapolate this growth to exponentially large times. 21 It was suggested to us by P. Calabrese that the behavior we observe here may be related to the log log R a correction to the one interval entropy in the vacuum discussed in [32,33,34]. 22 This is somewhat subtle, as the zero mode does not have a normalizable ground state.
for the finite chain analyzed in Fig. 4, the zero mode is absent due to the Dirichlet boundary condition, see also [34]. While the detailed exploration of complicated entangling regions and finite volume systems is outside the scope of this work, in Fig. 4 we follow the time evolution for long times on a finite harmonic chain, where the interval is at one end of the chain. We have included this geometry to demonstrate that the quasiparticle picture continues to hold in more complicated setups. The entropy exhibits exact revivals with profile exactly in agreement with the quasiparticle model, which we obtain by mirroring the chain at each end infinite amount of times. conditions on both ends. The region at the end of the chain is of size 2R/a = 300, and the full chain is chosen to be 2L/a = 1050 long in order to avoid any special ratio L/R. We follow the time evolution for a long time and find exact revivals. Note that because the region is at the end of the chain, the slope of the curves is half of that in (2.37).

Strips in 2 and 3 spatial dimensions
For the strip geometry in spatial dimensions d ≥ 2 we can decompose the fields in momentum modes transverse to the entangling surface. Let us denote the entanglement entropy of a massive scalar field of an interval by S I (R, t, β, m), where β is the effective temperature in the quench, 2R is the width of the interval, and t is the time. Then in d spatial dimensions (d ⊥ = d − 1 transverse dimensions) the entropy of the strip is (for details see Appendix A.1): where A ⊥ is the cross-sectional area of the sides of the strip and k ⊥ is the transverse momentum running parallel to the sides of the strip. We use this formula to compute the entropy numerically. To get a quantity with a well-defined continuum limit we make the subtaction (1.2).
The results are collected in Fig. 5. In the figures we have also plotted the time evolution of Rényi entropies, in addition to the von Neumann entropies. The results are compared to the predictions of the quasiparticle model (2.37), and precise agreement is found. By precise agreement, we mean up to area law terms subleading in the large region limit, which are not accounted for in the quasiparticle model. In the graphs below we have allowed ourselves to shift the numerical data points by an arbitrary constant to match the quasiparticle prediction.
We have checked that as we increase the region size this shift scales as the area, and thus is negligible in the limit of large region sizes, see Fig. 9 for a demonstration of this in the particular case of a boundary state quench for a spherical geometry, which are discussed in a following section. We expect similar results for the strip geometries. 23 The attentive reader may notice some deviation from the quasiparticle curve at early times, t ∼ β in Fig. 5. Such times do not obey the double scaling limit (1.5), hence we do not expect a precise match between the numerical results and quasiparticle curve. In particular until t ∼ β the entropy grows quadratically [6,13,14,31], while the quasiparticle curve exhibits linear growth (1.6). By smearing the time of origin of the EPR pairs, one can incorporate this quadratic growth into the quasiparticle model [8], but we chose to work with the simplest version of the model, which does not involve any adjustable parameters. 24 Next we focus on an important aspect of the entropy growth, the entanglement velocity defined in (1.6). Because the quasiparticle model predicts exact linear growth until t = R, and because at early times we observe more deviation from linearity, we extract v E from the slope of the curve at t = R: Numerical results are given in Fig. 6 based on this equation, and they show very good agreement with the quasiparticle value (2.39) even for fractional dimensions.

Spheres in 2 and 3 spatial dimensions
We complete the presentation of the numerical results with the sphere geometry. As discussed in the Introduction, it is important to consider different geometries for the entangling region to confirm the quasiparticle model. The sphere is a particularly nice case to analyze because of its symmetries, which make the numerical computations for large regions possible. See Appendix A.2 for details of the setup.
In Fig. 7 we plot the results of a boundary state quench in spatial dimensions d = 2 and d = 3, and in Fig. 8 we plot the results of a mass quench in spatial dimensions d = 2. All of these closely match the quasiparticle expectations for all times. Two highlights are that the linear regime is governed by the entanglement velocity (2.39), and the saturation time is t S = R (1.7). We note that at late times we again see a logarithmic rise of the entropy after t S , as in the one interval case in d = 1. This growth is most pronounced on Fig. 8, but the volume law in d = 2 provides more suppression than in d = 1. 25 Finally, for the example of a boundary state quench in d = 2, in Fig. 9 we demonstrate that the additional shift we apply to the numerical data points to get a closer fit with the quasiparticle model curves obeys the area law, hence it is subleading for large regions.   : Area law subtraction for a boundary state quench in a spherical (disk) geometry in 2+1 dimensions for different radii R of the disk. Because for spherical geometries t S = R, we require a match between the numerical data points and the quasiparticle curve at t = R. Then the shift that we apply to the numerical data isŜ(t = R) − S qp (t = R). It is linear in the radius R of the disk, as expected. relation when k ⊥ = k ⊥ . When we take k ⊥ = k ⊥ the delta function gives an infrared divergence that is regulated by assuming a finite size transverse region area A ⊥ , i.e., This motivates re-scaling the fields as follows The corresponding commutation relation appears canonically normalized for a field dependent only on the x variable, with k ⊥ just an external parameter The Hamiltonian may then be re-written in terms of these new fields as where H 1 is the Hamiltonian of a massive scalar field (of mass where in this last equation we have suppressed writing the dependence on the transverse momenta k ⊥ to emphasize that this is really just a Hamiltonian defined in 1 spatial dimension with axis running from one side of strip to the other.
The discretization of modes in 1 spatial dimension is rather straightforward, as we now explain. Slightly more complicated discretization that are relevant to the disk or sphere will be discussed in the next subsection. With lattice spacing a, the Hamiltonian in 1 spatial dimension is where the physical lattice points are x = j a, an IR cutoff is Λ IR = N a, and we have defined the dimensionless mass M a ≡ M a. We can also impose the conditions q(N + 1) = q(1) and p(N + 1) = p(1) to impose periodic boundary conditions where necessary. The non-zero elements of the K matrix of eq. (2.2) are given by K jj = 2 + M 2 a K j,j+1 =K j+1,j = −1 , (A.11) and if periodic boundary conditions are imposed K 1N = K N 1 = −1.
Since the different k ⊥ -modes decouple, and we are only tracing over the strip in the x-direction, rather than the x ⊥ direction, then just as the Hamiltonian decomposes as in (A.8) then so too does the entanglement entropy Similar results go through for the boundary state quench, where the entropy depends on the inverse temperature β, as given in (3.2).
We note that since the entropy S I is only a function of the magnitude and not the direction of the transverse momenta k 2 ⊥ in (A.12), then the angular integration of k ⊥ is trivial. This gives a factor of the area of the d − 2 dimensional unit sphere, i.e.,

A.2 Spherical coordinates and spherical harmonics
In 2+1 dimensions, for a spherical (disk) geometry, we use the Fourier expansions of the field φ and conjugate momentum π φ(r, θ) = Since φ is real, φ * = φ − , and similarly for π. Note that the harmonics φ , π satisfy the canonical commutation relations [φ (r), π (r )] = i 2πr δ + δ(r − r ) . We may now discretize this Hamiltonian with a uniform lattice in the radial direction: where a is the lattice spacing and r = ja, we introduced an IR cutoff Λ IR = N a, and we have defined the dimensionless mass m a ≡ m a. The radius of the disk is taken to be: So in 2+1 dimensions the non-zero elements of the K matrix for the discrete Hamiltonian, which was defined earlier in (2.2), are K 11 = 3 2 + 2 + m 2 a K jj = 2 + 2 j 2 + m 2 a K j,j+1 = K j+1,j = − j + 1/2 j(j + 1) .

(A.23)
This matrix and its eigenvalues form the basis for the numerical computations in the correlator method for the entanglement entropy. In the counting for the different modes in the entropy calculations, we must sum over all ≥ 0, with the = 0 mode getting a factor of 1 and the other modes a factor of 2: In 3+1 dimensions, for a spherical geometry, the development is similar to 2+1 dimensions, with some important differences. We use the expansions of the field φ and conjugate momentum π in terms of spherical harmonics As in 2+1 dimensions, this matrix and its eigenvalues form the basis for the numerical computations in the correlator method for the entanglement entropy. The entropy is finally given by the sum over the entropies coming from each angular momentum sector: (2 + 1) S . (A.32)