Reduced density matrix and entanglement in interacting quantum field theory with Hamiltonian truncation

Entanglement is the fundamental difference between classical and quantum systems and has become one of the guiding principles in the exploration of high- and low-energy physics. The calculation of entanglement entropies in interacting quantum field theories, however, remains challenging. Here, we present the first method for the explicit computation of reduced density matrices of interacting quantum field theories using truncated Hamiltonian methods. The method is based on constructing an isomorphism between the Hilbert space of the full system and the tensor product of Hilbert spaces of sub-intervals. This naturally enables the computation of the von Neumann and arbitrary R\'enyi entanglement entropies as well as mutual information, logarithmic negativity and other measures of entanglement. Our method is applicable to equilibrium states and non-equilibrium evolution in real time. It is model independent and can be applied to any Hamiltonian truncation method that uses a free basis expansion. We benchmark the method on the free Klein-Gordon theory finding excellent agreement with the analytic results. We further demonstrate its potential on the interacting sine-Gordon model, studying the scaling of von Neumann entropy in ground states and real time dynamics following quenches of the model.

which is obtained by tracing out the degrees of freedom in B. We can define the von Neumann (vN) entropy of subsystem A which satisfies the properties of an entanglement monotone and for a pure state ρ acts as a measure of the entanglement between A and B. The logarithm of the density matrix is often difficult to compute. Thus, it is convenient to study also Rényi entropies from which the α → 1 limit recovers the vN entropy. More generally, the complete entanglement properties of a state are encoded in the entanglement Hamiltonian H A defined by ρ A = e −H A , and thus playing the role of a generating Hamiltonian for the reduced density matrix [8].
The computation of partial traces is very natural in lattice systems where the Hilbert space inherently has a local tensor product structure. The same calculation is much more challenging in the case of field theories. Here, the number of degrees of freedom is continuous and the reduced density matrix can only be defined through the path integral formalism while entropies only become meaningfully defined after the introduction of an ultraviolet (UV) cutoff. A widely used approach to compute the entanglement entropies is the replica trick [9]. This procedure can be carried out for free theories [10], conformal field theories (CFT) [11] and integrable theories [12]. It has been argued that some results apply also to the general non-integrable case [13], yet it remains an open question. Another field theoretical approach consists of the covariance matrix formalism where entanglement entropies are computed from the correlation functions of Gaussian states [10,14]. Finally, the third widely used approach is to approximate a field theory with a lattice system and use discrete methods like the density matrix renormalisation group (DMRG) to compute the entanglement measures. This has been extremely successful for computing equilibrium properties [3,4,15,16] but is suffering from severe limitations when computing nonequilibrium dynamics due to the exponentially increasing bond dimension [15]. In this work, we develop a method for computing reduced density matrices and entanglement measures in interacting field theories both in and out of equilibrium.
While HT is very successful at computing spectral properties and non-equilibrium time evolution, it has not been the most convenient choice for computing entanglement related quantities. In preceding works [43,44] HT was used to calculate matrix elements between higher excited states. By using analytic replica techniques, correlation functions of twist fields and other entanglement related objects were calculated for those states. While such approaches proved useful at computing low Rényi entropies, the calculation of vN entropies and entanglement negativity remained out of reach.
In this work, we develop a general way to construct reduced density matrices with HT. The output of our method is the density matrix of a state explicitly represented in a computational basis which is a tensor product of the basis of the left and right subsystems. This enables the direct computation of almost any entanglement related quantity: vN entropy, entanglement negativity, mutual information and the direct study of the entanglement Hamiltonian and reduced density matrix of an interacting field theory itself. Our method is general and can be widely used on top of any HT code that uses expansions in free bases (see Section II) which is a common choice in modern applications [19,20,22,24,26]. This enables us to take full advantage of the power of HT for real time evolution of a wide range of interacting QFT and study the whole spectrum of entanglement related quantities without needing to approximate the theory with a lattice system. It gives access to the entanglement properties of ground, excited and time dependent nonequilibrium states as well as thermodynamic entropies of thermal states. Additionally, the method has the potential to work also in D > 1 + 1.
The manuscript is organised as follows: in sec. II, we briefly introduce the basic concepts of HT. In sec. III we present our method for computing reduced density matrices. We begin in sec. III A with the theoretical construction and continue in sec. III B outlining an efficient algorithm for the numerical implementation. In sec. IV we introduce the QFT models that we test the method on. In sec. V we present the results and a comparison against analytical predictions. Sec. V A focuses on equilibrium states while sec. V B covers non-equilibrium dynamics. We conclude in sec VI with an overview, discussion and the scope for the future work. The appendices cover the more technical details of the method.

II. HAMILTONIAN TRUNCATION
HT is a numerical method for strongly interacting QFT first introduced in the 90's by Yurov and Zamolodchikov [45,46]. It is based on the Hamiltonian formalism and the idea is to represent the Hamiltonian of a field theory defined on a compact domain as where H solv is the solvable part of the Hamiltonian and Φ pert is a perturbing potential. Traditionally, the CFT of the UV fixed point of the theory was used as H solv but more modern approaches consist of using other solvable theories like free massless and massive theories. The perturbing potential Φ pert does not need to be small which gives HT the power to capture non-perturbative effects. The method proceeds with representing the potential Φ pert as a matrix in the Hilbert space of H solv , the space H solv . The crucial step of HT is to introduce a high energy cutoff, keeping only the low energy states of H solv which renders the matrices finite and enables numerical computation. The method converges if Φ pert does not mix significantly the low energy sector of H solv with the higher energy sectors. In case of an expansion around the CFT point, this is guaranteed by the renormalisation group theory for relevant perturbations Φ pert .
If computed for several high enough cutoffs the results of a HT simulation can often be extrapolated to obtain the infinite cutoff value. Alternatively, a numerical renormalisation group algorithm can be used [47]. Computing entanglement related quantities has been challenging for HT. The HT Hilbert space is usually spanned either by primary and descendant CFT states or free model eigenstates in momentum basis. It does not allow for an easy bipartition in position space. Earlier approaches [43,44] were based on mapping the problem to the CFT calculation of entanglement related CFT objects for descendant fields. While conceptually elegant, such calculations for higher descendant states are often tedious and are associated with several restrictions. They have so far been limited to the first few Rényi entropies. In this work we want to overcome this problem and construct a more general and robust approach which can be exploited for the calculation of almost any entanglement related quantity.

III. METHOD
Our main goal is the construction of reduced density matrices with Hamiltonian truncation. We consider a field theory defined on a finite interval F = [0, L] (full) with open boundary conditions and are interested in computing the entanglement between a subsystem L = [0, ] (left) and its complement R = [ , L] (right). Note that we use L for the full system size and as a label for the left subsystem. The interpretation is clear from the context. Using HT, we can compute the density matrix ρ of a ground, excited, thermal or non-equilibrium state of the theory expressed in the Hilbert space H F of the full interval F . However, tracing out a spatial part of the system is difficult in the momentum basis of H F . If we express ρ as ρ LR in a Hilbert space H L ⊗ H R built out of Hilbert spaces H L and H R on the intervals L and R, we can trace one part of the system directly. We thus need to construct H F , H L and H R and find the unitary transformation to compute The idea of the method is visualized in Figure 1.
A. Splitting the System

Fields and Hilbert spaces
We start with the construction of the Hilbert space on the full interval H F . We expand the fields of the free theory (massless or massive) in terms of momentum modes. The Fock space generated by the mode creation operators serves as the computational basis. For concreteness, we choose an expansion around a massless free bosonic field theory with Dirichlet boundary conditions (φ(0) = φ(L) = 0) at the edges. The procedure is easily generalisable to expansions around massive free theories and other boundary conditions and could by construction be applied also in D > 1 + 1.
The field expansion can be written as with p k = k π L and A k the bosonic modes on the full interval fulfilling the commutation relations [A k , A l ] = A † k , A † l = 0 and A k , A † l = δ k,l . We refer to the modes A k as full modes in the rest of the text.
The full modes A k span the Hilbert space H F | n F ≡ |n 1 , n 2 , . . . with n k the bosonic occupation numbers, the normalisation N F = k>0 √ n k ! and |0 the vacuum of the massive free boson theory.
A cut at position divides the system into two subsystems L and R (left and right), as shown in Figure 1. In a similar fashion as H, we construct H L ⊗ H R by defining two fields φ L and φ R (split fields) living on the two subintervals and quantizing them.
The formulation of the fields on the intervals depends on the additional boundary conditions that we introduce at the cut. We choose to study Neumann (∂ x φ L ( ) = ∂ x φ R ( ) = 0) or Dirichlet (φ L ( ) = φ R ( ) = 0) boundary conditions at the cut. In the main text, we focus on Neumann boundary conditions at the cut. The treatment of Dirichlet boundary conditions at the cut is detailed in appendix B. The outer edges of the system are always fixed to be Dirichlet boundary conditions (φ(0) = φ(L) = 0).
For Neumann boundary conditions at the cut, the fields on the intervals are defined as where q and a σ m are the bosonic annihilation operators on the two partitions for σ ∈ {L, R}. Both, the fields φ L , φ R and the modes a L/R defined in (9) and (10) fulfill the respective bosonic commutation relations. Fields and modes on different subintervals commute. In analogy to the full fields, the modes on the sub-intervals span their respective Hilbert spaces H L and H R . The computational bases for the two sub-intervals are with the normalisation N σ = m>0 n m,σ !. The choice of mixed boundary conditions (Dirichlet on the edges and Neumann at the cut) has the advantage that no zero-modes appear in the system. The vacua of the sub-intervals are not equal to each other and in particular they are not equal to the full system vacuum |0 L = |0 R = |0 . The product space H L ⊗ H R is then generated by | n L , n R ≡ | n L ⊗| n R on top of the vacuum |0, 0 ≡ |0 L ⊗ |0 R = |0 . Before we continue, a couple of words on conventions and notation. Throughout the paper, we will use k as an index for the full modes A k and l, m as indices for partial modes a σ m . Greek indices always indicate a left or a right partition, σ ∈ {L, R}.

Bogoliubov transformation
At first glance it might not be obvious that the the descriptions of the system in terms of the full field and split fields are equivalent. From an intuitive point of view: for any given field configuration of φ, one can find a configuration of φ L and φ R that is arbitrarily close to φ and still respects the boundary condition at the cut. Indeed, we can choose φ L and φ R to be equal to φ everywhere except for a small neighborhood of the cut. There, they have to deviate in order to satisfy the boundary condition. But we can make this neighborhood arbitrarily small while still preserving the continuity of φ L and φ R and the boundary conditions. Later, we give a more detailed argument for the correspondence on the algebraic level.
We now formally construct the unitary mapping between H F and H L ⊗ H R proposed in (5) In order to compute the matrix elements (12), we take two distinct steps. Firstly, we express the full modes A k in terms of the partial modes {a σ m } m,σ and {a σ † n } n,σ . Secondly, we rewrite the full vacuum in terms of the partial vacua and partial modes. The latter is particularly important because neglecting that the Hilbert spaces are not defined on top of the same vacuum will lead to wrong results.
We rewrite the full modes as where the coefficients γ are to be determined. The coefficients γ are the result of identifying the fields on the full interval with the split fields This identification, the continuity condition, is the core of the unitary map between the full and the split Hilbert spaces. We can express the full modes A k in terms of the fields φ(x, t) and the canonical momentum operator π(x, t) = ∂ ∂t φ(x, t) Combining (15) with the continuity condition (14) links the full system modes A k to the partial modes a σ m . Evaluating the integral, we obtain the coefficients γ for Neumann boundary conditions at the cut as The special cases in equations (16) and (18) are divergences of the integrand. The detailed calculation as well as the expressions for the Dirichlet boundary conditions at the cut can be found in appendix B.

Multimode squeezed coherent vacuum
When expressing the full modes as a superposition of partial modes, we also have to re-express the vacuum of the system. In order to find a formulation of the full system vacuum in terms of the partial modes, we identify Equation (13) as a multi-mode Bogoliubov transformation [48] with A = (A 1 , . . . , A N F ), a = (a L 1 , . . . , a L N L , . . . , a R N R ) and Note that M is not an operator here, but a matrix of numbers. Since all the coefficients γ in equations (16)- (19) are real, we focus on the case of real u and v. We use the same symbols as in (17)- (18) without the subscript indices to refer to matrices of coefficients. For ease of notation, we still express equations in terms of u = γ L,+ γ R,+ and v = γ L,− γ R,− . The transformation (21) expresses bosonic modes A in terms of different bosonic modes a. Thus, the transformation must preserve the commutation relations. These are encoded in the symplectic structure of M It can be verified that γ coefficients in eq. (16)- (19) obey the symplectic structure of the Bogoliubov transformation. The Bogoliubov transform (21) is equivalent to a unitary transformation [ with In contrast to M , U is an operator and not a matrix of numbers. Thus, the vacuum of the full modes |0 can be expressed in terms of the vacuum of the partial modes |0, 0 as because then A k |0 = U † aU U † |0, 0 = 0. U can be written in a more convenient form for actual computations, the so-called disentangling form with Since we have an expression of the full modes in terms of the partial modes and an expression of the full vacuum in terms of the partial vacua, we can compute the matrix elements of U T . The elements of U T are overlaps between states in the full basis | n F and the split basis | n L , n R . with The first bracket in (28) builds the occupation number state n L , n R | from the partial vacuum 0, 0|. The order of left and right creation operators does not matter here, since they commute as they act on different partitions. The second bracket represents the operators A † k which build | n F on top of the vacuum of the full modes |0 . We choose to express the full modes in terms of the partial modes [cf. Equation (13)]. The opposite way of expressing the partial modes in terms of full modes would work as well. The last bracket in (28) expresses the vacuum of full modes in terms of the split modes. Looking at formulation of the norm in (29), we recognize the first term from the transformation of the vacuum; it is constant for all matrix elements. The second and third term are the norms of the partial and the full occupation number states. The contributions of the exponentials with σ and τ in (26) vanish due to the order of operators in the exponentials and their action on the vacuum to the right.

Equivalence of Hilbert spaces
We have now fully exposed the unitary transformation between H and H L ⊗ H R . We are thus ready to justify why the two descriptions of a physical state, in terms of ρ and in terms of ρ LR , are equivalent despite the additional boundary condition at the cut.
The way to think of the transformation is in terms of Fourier analysis. The sine modes of the split intervals in eqs. (9) and (10) serve as a functional basis in which a field configuration of the full interval field is expanded. The fact that it is indeed a basis, is a simple consequence of Carleson's theorem for convergence of Fourier series [49,50].
However, we are dealing with quantum fields not just simple scalar valued functions. Upon quantization, the Fourier coefficients in the expansion become operator valued as indicated in eqs. (9) and (10). The symplectic structure of the Bogoliubov transform between the two algebras (22) guarantees that the two quantizations of the system, the one in terms of the full field in eq. (7) and the one in terms of the split fields, are equivalent.
Thus in the complete, infinite dimensional Hilbert spaces, the unitary map between H and H L ⊗ H R is exact. However, once we introduce a truncation, it becomes an approximation in the same fashion as HT is always an approximation of the quantum state using the low energy sector of the Hilbert space. Using the partial field expansion (10) then becomes in spirit very similar to using a truncated Fourier series to approximate a function. In section V we demonstrate that such an approximation indeed performs excellently at computing entanglement entropies.

Truncation
The mapping between the full system and the partitioned system is so far formulated without considering the truncation. In order to be able to implement the method with a finite amount of computer memory, we have to consider a finite dimensional approximation of the Hilbert spaces. In HT, it is often important to choose the most suitable truncation scheme for the problem. In our case, we have to introduce three truncations: of H, H L and H R . For convenience, we choose the cutoffs such that the Bogoliubov transformation (21) becomes a square matrix. This leads to the restriction s F = s L + s R , where s F , s L and s R are the number of momentum modes kept in the full system, the left and the right partition, respectively. We truncate all Hilbert spaces, H, H L and H R with an energy cutoff such that the cutoff energy is equal to the energy of the single excitation of the largest momentum mode kept, s F , s L and s R respectively. For a free massless boson Hamiltonian, the energy of the state | n F defined in eq. (8) above the ground state is ( n F ) − (0) = π L ∞ k=1 kn k . Similarly for the states in the split Hilbert spaces (11), where L is replaced by and L − , respectively.
In infinite dimensional Hilbert spaces, before the truncation, the symplectic structure (22) is fulfilled exactly. By introducing the cut-off this property can be compromised. We use the preservation of the symplectic structure as a guiding principle to select a cut-off scheme.
We investigate two different cut-off schemes. In the first one, we distribute the number of partial modes equally across the left and the right partition independent of the position of the cut (fixed cutoff): s L = s R = s F /2. This leads to an easy implementation but is questionable from a physical standpoint. A short interval with many modes leads to a greater resolution in position space than a large interval with the same number of modes. Furthermore, such a non-uniform UV cut-off leads to a position-dependent non-universal constant in the entanglement entropy that obscures the true functional dependence. The second cut-off scheme takes into account the length of the partition. The number of modes are distributed proportionally to the length of the interval s L = L s F and s R = s F − s L (constant mode density). This scheme keeps a constant density of momentum modes (s L / = s R /(L − )) and thus a homogeneous UV cut-off. Both cut-off schemes are displayed in Figure 2.
All further considerations use the constant mode density cut-off scheme. It indeed reproduces the bosonic commutation relations more faithfully as shown in appendix A where also further details of the cutoff scheme are discussed.

B. Algorithm
The evaluation of the overlap as given in (28) is difficult because the expression is not normal ordered and the sums in the expression are taken to the powers n k . It would require an overwhelming number of commutations of individual mode operators.

Generating functional formulation
The expression simplifies if we express it in the spirit of a generating functional. The repeated application of a mode a is equivalent to where J is a scalar variable.
Inserting the identity (30) in the formulation of the matrix elements of U T in (28), we obtain an expression resembling a generating functional Here, S are all terms related to the split modes, F is the term that is associated with the full system modes and V are the terms that build the vacuum. We introduced two kinds of additional, scalar variables. J k are the additional variables for the full modes and j m,σ are defined for the partial modes.
Using the Baker-Campbell-Hausdorff (BCH) relations [51], we can bring this expression into normal order and evaluate the expectation value. The series of commutators in the BCH relations terminate at most at the second order since the highest power of mode operators in the exponent is two. We can write the normal ordered expression as with The terms ComF, ComSF and ComAV are results of the BCH operations. They are defined as The commutator [A, V ] is linear in creation operators. The derivation of the commutators is detailed in appendix C 1.
The problem of computing the overlap reduces to computing multiple derivatives of a scalar expression if we express (28) as the derivative of a generating functional

Tackling the exponential complexity of differentiation
The next goal is the efficient calculation of all derivatives in (40). The pure symbolic evaluation of derivatives becomes prohibitively expensive with increasing cut-off. The number of terms grows as n!/ 2 n/2 (n/2)! with the number of derivatives n. For large n, this scales as exp [(n/2) (log(n/2) − 1/4)]. In the following, we show how the structure of the generating functional helps us to make the evaluation more efficient.
From an algorithmic point of view, there are two distinct exponentially scaling problems involved in the computation. On the one hand, the size of the unitary transformation grows with the size of the Hilbert space. We have to evaluate exponentially many terms in order to fill the matrix. It is impossible to circumvent this exponential since the method is based on exact diagonalization. On the other hand, each matrix element needs an increasing number of derivatives with increasing occupation numbers. The number of terms in the derivation scales also exponentially with the occupation number (and thus the cut-off). In the following section, we describe an algorithm to make the evaluation of the derivatives feasible for relevant cut-offs. We are not able to reduce the exponential growth of terms to a polynomial growth. However, we derive a procedure which largely reduces the exponential growth. Thus, it is possible to reach HT cut-offs that provide reasonable approximations of interesting physics.
A commonly used alternative to symbolic differentiation is automatic differentiation (AD) [52,53]. The algorithm of AD tracks the computation of the function and uses predefined derivatives of elementary functions to evaluate the derivative numerically. In our case, it is hard to use AD directly since we have to compute possibly very high derivatives of the function and that we do not need the actual function value. Furthermore, we can exploit the structure of the function to determine which terms must be 0 without computing them. Therefore, we take a more specialized approach and do not rely on AD.
By inspecting the structure of the expressions in Equation (40), we note that the derivative always acts on an expression of the form with T = ComF + ComSF + ComAV, a shorthand for all terms in the exponent of Z in (36). For the ensuing discussion, we introduce a shorthand notation for the derivatives of Z The expression T ≡ T [•, •] has two arguments because the commutators in T are always quadratic in J k and j m,σ . For the rest of the discussion of the algorithm, we will not distinguish j m,σ and J k . We can always write them in terms of a general J i by using i as a multi-index. This new notation helps us to demonstrate that many terms are 0 and we can drop them. Due to the commutativity of the derivatives and the step of setting J k = 0 in the end [cf. eq. (40)], we find All expressions that are not derived twice must be zero if we set all J k = 0 in the end because J k appears quadratically in each commutator. Thus, the number of T s for each derivative is given by N T = i ni 2 , where n i are the occupation numbers of the full-system state and the partitioned state.
The restrictions described above lead to a more efficient algorithm in comparison to symbolic derivation of the full expression. Considering the restrictions in equation (43), the result of the derivatives in (40) is heavily constrained. Every term must be derived twice (otherwise it is 0). Furthermore, we only sum over unique combinations since we can freely exchange the arguments of T and the order of the T [J l1 , J l2 ] in the product over l.
The input of the algorithm is a list of J i with corresponding powers n i and we only compute combinations of fully derived where c k are the multiplicities of the terms in T . The sum runs over all unique combinations of T . A combination is unique if it cannot be transformed into another combination of T s by swapping the arguments of T or commuting T s. This corresponds to iterating over all pairwise lexicographically ordered tuples of J . The exponents p kl are the powers of certain terms T if the same arguments (J l1 , J l2 ) appear multiple times in the same sequence. Thus, the number of terms in the product over l can vary depending on the number of individual combinations of (J l1 , J l2 ). Since we are not considering the full system and the split modes separately at the moment, the indices k, m, and l are used without further implications here. For a more detailed discussion of the prefactors c k , we refer to appendix C 3.
For concreteness, we consider a simple example of three modes J 1 , J 2 and J 3 to illustrate the procedure. We are interested in the second derivative with respect to each J i . In terms of occupation numbers, we can write the configuration as n = (2, 2, 2), where n is the vector of occupations numbers. The primed sum in (44) runs over all unique configurations of strings. In our example, there are five distinct configurations The tuples represent the arguments of T [•, •]. All other combinations other than those listed in eq. (45) can either be generated by swapping tuples or by exchanging the arguments inside of a tuple. A swap of two tuples is allowed due to the commutativity of multiplication in (44). The exchange of arguments is equivalent to exchanging the derivatives of a single T which corresponds to one of the identities in (43). Since there are six derivatives in total, we must have three distinct T terms in each string. We can express the five combinations in (45) more compactly with powers p kl Here, k is the index of the overall combination of all pairs J and l is the index of the tuple in the string. More concretely, p 2,1 = 2 because the second string contains (J 1 , J 2 ) 2 as first pair. Some of the configurations in (45) may appear multiple times during the application of the product rule in (44). Thus, we have to take care of the multiplicities in front of the terms. In our simple example, we can just list them as c = (1, 2, 2, 2, 8). Here, they are calculated by explicitly performing the derivatives on the left side of (44). The primed sum in (44) can be evaluated given all configurations in (45) and the vector c. All terms of the form T [J l1 , J l2 ] are numbers that can be evaluated by summing the derivatives of commutators in (36) explicitly.
As demonstrated in the example, the computation of (44) can be divided into two subproblems. Firstly, we have to determine all unique combinations of pairs (J l1 , J l2 ) for a given n. Secondly, we have to compute the coefficients c k given p kl and the tuples (J l1 , J l2 ).
The first task can be solved with a tree-based algorithm that is described in detail in appendix C 2. The idea is to build only the combinations of tuples (J l1 , J l2 ) that adhere to the uniqueness condition defined for the primed sum, i.e. lexicographical ordering of all index tuples. The condition can be checked locally at every node of the tree. Thus, only nodes that can still build valid configurations are expanded in subsequent operations. The trivial approach of listing all combinations of J i for a given n and filtering for the unique ones gets prohibitively costly already for low cut-offs.
The coefficients c k have a closed form expression and are given by where n i are the occupation numbers and N diag,k is the number of identical arguments for T in the string with index k. In our example of n = (2, 2, 2), N diag,1 = 3 and N diag,3 = 1. The proof of the equation is given in appendix C 3. Finally, we can put all the pieces together. An element of U T corresponds to the calculation of an overlap of the form n L , n R | n F . Each of the states is given as an occupation number vector. Equation (40) connects the occupation numbers to derivatives of a scalar function. These derivatives can be computed explicitly by first enumerating all unique configurations of the primed sum in (44). Each of the tuples in a configuration represents the arguments of T . The derivatives of T for some tuple (J l1 , J l2 ) can be evaluated explicitly. The product of all T values in a string is weighed by a factor (47) and summed to yield the final value of the matrix element. The explicit expressions for the derivatives of the commutators are given in appendix C 1.

A. Klein-Gordon model
For concreteness, we demonstrate the potential of our method on two well known QFT models. The first example is the Klein-Gordon (KG) model, the massive free boson theory, described by the Hamiltonian where φ(x) is a real scalar field and m is the boson mass. This free model serves as a perfect test bed for our method. Its entanglement properties are known analytically both from replica trick techniques [10,11] and from covariance matrix methods [10,14] including the equilibrium states and the non-equilibrium dynamics. For the massless case the following well known result has been derived [11]: where the central charge of the CFT c = 1, a is a UV cutoff, g is the Affleck-Ludwig boundary entropy [54] and U (a) is a non-universal constant dependent on the precise form of the cutoff.
Although being a free theory, the massive case of the KG model is the first nontrivial test example of our method. While the massless case is diagonal in our computational basis, the massive case is fully non-diagonal. Due to a finite correlation length ξ ∼ 1 m the entanglement for m > 1 L saturates to an area law plateau where S( ) = const. At distance closer than ξ to the boundaries, the curve S N interpolates smoothly to the zero value at the boundaries. For m < 1 L , there is a smooth crossover from a log law to an area law scaling of the entanglement entropy.
For thermal states, the vN entanglement entropy becomes the thermodynamic entropy and there is a smooth crossover with increasing temperature to a volume law S( ) ∝ . In non-equilibrium dynamics in the massless case, the vN entropy is expected to grow linearly in time [9]. In case of a finite system, the growth stops when excitations from the splitting point reach the boundaries of the system and one expects recurrent dynamics. In the massive case, the linear growth is superposed with an oscillatory component with a frequency given by the boson mass [55].
To generate analytical predictions for the KG model to compare our numerical method against, we employ the covariance matrix formalism [10,14]. It is a convenient framework because of its simplicity and because it also enables to model cutoff effects. Further details are outlined in appendix E.

B. Sine Gordon model
A paradigmatic model of strongly interacting QFT is the sine-Gordon (sG) model with the mass parameter m and the interaction parameter β. The sG model is one of the simplest models displaying confinement and is an integrable model solvable by S-matrix bootstrap techniques [56]. The model has solitonic topological excitations and a rich phase diagram. For β < √ 4π the interaction is attractive and the solitons form bound states -breathers. For √ 4π < β < √ 8π the interaction is repulsive, for the separating line β = √ 4π, the model can be mapped to a free Dirac fermion and at β ∼ √ 8π the model undergoes a Berezinskii-Kosterlitz-Thouless phase transition to a free model [56].
A convenient way to parameterize the sG interaction parameter in the attractive regime β < √ 4π is where the parameter λ is convenient because λ equals number of breathers present in the sG spectrum. The mass m n of the n-th breather is given by where M is the soliton mass. In particular, the mass of the lightest particle, the first breather, m 1 determines the gap of the system. In finite system size, these masses get modified and can be computed using the form factor and boundary bootstrap formalism [21,[57][58][59]. Each of the breathers has a tower of excited states as a result of acquiring a nonzero momentum. The set of allowed momentum values is discrete in finite volume. The expressions for finite volume breather energies are given in appendix F. The entanglement properties of the sG model in the repulsive regime have been studied by spectral form factor and corner transfer matrix techniques [60,61] and predict the height of the vN entropy area law plateau where M is the soliton mass which is a function of m and β. In the attractive regime, the entanglement properties are less understood. Based on general arguments for gapped systems, the vN entropy plateau is expected to follow the form [13]: where K 0 is the modified Bessel function, c is the central charge of the UV critical point, m α are the masses of the particles in the spectrum (breathers in the sG case), ξ 1 the correlation length corresponding to the lightest particle and U a constant.
Concerning the non-equilibrium dynamics, it has been recently shown using form factor techniques for small quenches of the sG model that the entanglement entropy exhibits non damped oscillations in time with frequencies corresponding to even breather masses [62].
Here we implement a HT for the sG as developed in [19,21]. We list all the HT matrix elements used in appendix D.

V. RESULTS
The following section is structured in two main parts. In the first part, we show results of the method for the Klein Gordon model in equilibrium. These results are compared to covariance matrix calculations and serve as a benchmark. Additionally, we show results for the interacting sine Gordon model to demonstrate that our method works beyond the free regime. The second part of the result section contains non-equilibrium evolution of the von Neumann entropy in real time for Klein Gordon and sine Gordon models.
All results shown in this section are computed for a finite cut-off s F = 18. This corresponds to 1597 states. The boundary conditions at the cut are chosen to be Neumann while the (physical) boundary conditions at the outer edges are Dirichlet. A constant mode density truncation scheme is used in all computations. Further details on the cut-off scheme are described in appendix A. Normalisations (like the prefactor in (29)) are enforced by normalising the reduced density matrix numerically.

A. Equilibrium
The Klein Gordon model (48) represents a non-trivial check for HT because its Hamiltonian is non-diagonal when expanded in the massless (CFT) basis for any mass m = 0.
As an initial check for the splitting procedure, we reproduce the correlations of the Klein Gordon theory in terms of the split modes. The correlations φ(x)φ(L − x) can be either calculated in terms of the full fields φ acting on the full interval density matrix ρ or in terms of the split fields of the left and right partition φ L and φ R acting on the partitioned density matrix ρ LR : Here, we use the notation φ L/R to refer to the field on the sub-interval that x belongs to. Figure 3 compares the correlations across the full range of the system for a cut at position /L = 1/3 for Neumann and Dirichlet boundary conditions at the cut. Both of the split field curves agree well with the correlations of the full system. In the case of Neumann boundary conditions at the cut, we only observe deviations at the cut. A plateau forms around the split at , since we impose ∂ x φ = 0. The Dirichlet boundary conditions enforce φ = 0 at and we notice that the correlations drop to zero as expected. The figure is symmetric around /L = 0.5 due to choice of arguments in the correlator. The overall wavy features in the curve for the full and the partial modes are a feature of the finite cut-off in HT. With an increase in the cut-off, we expect these features to reduce in amplitude.
Since correlations with Neumann boundary conditions at the cut agree better with the full correlations, we choose Neumann boundary conditions at the cut for all further entropy computations. We expect Dirichlet boundary conditions at the cut to be eventually equivalent to the choice of Neumann boundary conditions for higher cutoffs (cf. section III A 4). We continue checking the performance of our method by calculating the von Neumann entropy. We compare the von Neumann entropy with an analytic calculation using the covariance matrix approach (cf. Figure 4). The formalism is explained in detail in section E. All covariance matrix computations in this section are performed using 200 momentum modes. HT entropies are calculated at all points /L = n/s F , n = 1, . . . , s F − 1 since the bosonic commutation relations in the truncated split basis are fulfilled best at these points (cf. appendix A). The calculation of the entropy at other points is possible, but will result in more significant errors due to the truncation effect leading to worse preservation of the canonical commutation relation by the splitting procedure. The covariance matrix results (dashed lines) and the CFT results (solid lines) in the figure are shifted by a constant to coincide with the HT curves at /L = 0.5 for ease of comparison. This accounts for the non-universal cutoff dependent constant (see for example eq. (49)) which is slightly different in the analytic and the HT case due to the different truncation schemes.
In all cases, our method agrees excellently with the analytic predictions. The massless boson shows the expected logarithmic growth in entropy. This agrees perfectly with the CFT prediction [11], eq. (49).
With increasing mass, the curve develops a flat plateau in the central region, transitioning to the area law regime as expected for a massive boson. For distances less than a correlation length away from the boundary, the curve undergoes a non-linear behaviour before it reaches 0 at the boundaries due to finite size effects.
Similar data can be obtained for Dirichlet boundary conditions at the cut. We present the results for Neumann boundary conditions here since they show better agreement with the expectation. The Dirichlet data has slightly stronger deviations close to the boundaries.
HT provides access to the reduced density matrix at arbitrary temperatures so in addition to ground state properties, we can also access the von Neumann entropy of thermal states. For T > 0, the vN entropy coincides with the classical thermodynamic entropy. Figure 5 shows the entanglement entropy (T = 0) and the thermodynamic entropy (T > 0) of a massive (m = 5) free boson at different spatial positions. As before, the results obtained by our method agree well with the covariance matrix computation. The dashed curves are again shifted to coincide at /L = 0.5 to account for cut-off dependent constants. At T = 0, the we see the expected plateau of the area law of the entanglement entropy. At a finite temperature, the entropy becomes extensive and grows linear with the system size.
The basis transformation from a full to a split system is independent of the model. Thus, we apply the same methodology to the interacting sine-Gordon Hamiltonian (50) as shown in Figure 6. We compare the curves for two different values of the coupling parameter λ and different values of the soliton mass M . The cases of λ = 7, M = 25 and λ = 17, M = 60.29 are chosen such that the gap of the model (the mass of the first breather) matches. In comparison to the Klein-Gordon model, we do not see the onset of a plateau in the middle of the curve. For the matching breather mass case, the gap is m 1 = 11.13L meaning that the correlation length is less than one tenth of the system size. At such a short correlation length, an area law plateau would generally be expected. The log-like deviation from that could be indicative of longer range entanglement in the sG case which could be a consequence of the topological nature of solitons or a subtlety of the continuum missed by discrete calculations. It would be interesting to further understand this surprising scaling with analytical tools.
The perfect overlap of the curves λ = 7, M = 25 and λ = 17, M = 60.29 indicates that the vN entropy scaling in the attractive regime of the sG model is dominated solely by the first breather and not by the higher particles in the spectrum. This is consistent with the general expression (54). At large volumes, the K 0 corrections are highly suppressed, resulting in the value of vN entropy depending only on the correlation length.

B. Real time Dynamics
We continue by using our method to study the real time dynamics of the vN entropy following quenches.
We begin with the analytically tractable KG model. The mass quenches with several increasing post-quench masses are shown in Figure 7. We study the dynamics of the vN entropy between the left and the right half of the system ( = 0.5) and compare the HT results with the analytical results from the covariance matrix formalism. We displace the curves by a constant such that they start from the same point. This is to account for the nonuniversal cutoff dependent constant resulting form the difference in truncation schemes in the two methods.
In the quenches to the massless post-quench Hamiltonian, we observe the expected CFT linear growth of the vN entropy. The linear growth is interrupted at t = L/2 by a reflection when the quasiparticles from the cut reach the system boundaries. At t = L this results in a recurrence and the free dispersionless nature of the model leads to periodic dynamics. This is shown in panel (a).
At nonzero mass, the vN entropy develops an oscillatory component with a frequency proportional to the boson mass m. from the boundaries are still prominent. Thus, the linear growth becomes obscured by them which is what we see in panel (b). For intermediate masses for which the correlation length is an order of magnitude but not more smaller than the system size, the linear growth becomes visible as shown in panel (c). However, the plateau keeps undergoing significant oscillations due to reflections from the boundary.
The KG plots in Figure 7 expose the limitations of the truncated Hamiltonian approximation. At smaller masses [panels (a) and (b)], the HT results with our method match perfectly the analytic prediction up to times several times longer than the system size. This shows that the HT calculation of real time dynamics can be very reliable up to considerably long times. At higher masses [panel (c)], the real time dynamics start to deviate from the analytic curve for late times and the curve develops a phase shift. This is due to truncation effects -at higher masses the low energy part of the Hilbert space becomes too small to accommodate all the relevant modes for the dynamics. The quality of the time evolution depends also on the amplitude of the quench, the difference between the pre-and the postquench mass. For small quenches, the HT evolution is reliable even at large masses and for bigger quenches it gets less reliable also at smaller masses. This is because a large quench generates excitations high up in the spectrum, exceeding the HT truncation. The very high masses shown in panel (d) cannot be reliably simulated with our current implementation of the HT and we show only the analytic curve to support the discussions in the previous paragraphs. Such high masses could be implemented also with our methods, though, if a massive basis was chosen for the HT expansion instead of the massless basis. In this case, HT becomes exact also at nonzero masses.
We turn to studying the interacting sG dynamics. The vN entropy dynamics of the sG quench is shown in Figure 8 and compared with the KG quench. The comparison is done at such choices of the parameters that the gaps of the two systems agree. We observe an oscillatory motion as predicted by recent work by Castro- The spectrum is obtained using a discrete Fourier transform. The amplitude at frequencies ω is compared against energy levels of (moving) breathers. The breather energies are computed using reflection factors [21,[57][58][59], the expressions are listed in Appendix F. Due to the charge conjugation symmetry, only C-even breather state frequencies appear. The inset shows the original time evolution of the quench.
Alvaredo [62]. As known previously in the literature [55] and also demonstrated here in Figure 7, the oscillating dynamics is a generic consequence of the gap and in not a special feature of interaction. From our present results it is however not yet possible to determine whether the oscillations in sG quenches remain undamped at longer times as predicted by [62]. In Figure 9 we perform the Fourier analysis of the time series and compare the frequency spectra with breather energy levels. In order to have a more reliable time evolution at longer times, we study a small quench in mass -a quench generated by a moderate change of the soliton mass. The analytical breather energies for comparison are computed using the form factor and boundary bootstrap formalism in references [21,[57][58][59], for completeness the expressions are listed also in Appendix F. The sG ground states which are the prequench states are even under the charge conjugation C : φ → −φ which interchanges solitons with anti-solitons. Therefore, as predicted by [62] and confirmed also by our results, only C even states get populated during the quench. These include even breather states and even multiples of odd breathers.
Interesting questions remain whether the oscillations in the sG dynamics of vN entropy are damped and whether there is a linear growth superposed to oscillations as in the KG case or not. As shown in fig. 8, there seems to be a slight growth at early times but it is hard to determine reliably whether it is not just an oscillation of another slower frequency superposed on top of the higher frequency oscillations. In order to study both questions, the quenches would have to be computed at a much larger post quench soliton mass, allowing to ex-plore larger times before the reflection from the boundaries. Recently, an advanced HT implementation of the sG model has been developed [26] allowing for calculations with Hilbert space sizes of several hundred thousand states. It would be interesting to combine our method with such an approach to study sG quenches in large volume. This would, however, require even more efficient approaches to deal with the exponential complexity of derivatives discussed in Section III B 2.

VI. CONCLUSION
We presented a method to compute a reduced density matrix of a quantum field theory within the Hamiltonian truncation framework. Our method constructs a unitary transformation between the Hilbert space of the full system and a tensor product of Hilbert spaces corresponding to the subsystems. This maps the density matrix of a state into a form which is convenient for taking partial traces.
The method allows for the direct evaluation of a wide spectrum of entanglement related quantities, including von Neumann and Rényi entropies, mutual information, entanglement negativity and entanglement Hamiltonians. Our method makes it possible to study entanglement in ground, excited and thermal states as well as the real time evolution in non-equilibrium dynamics. Furthermore, our method is model-independent and can be applied for any HT that is based on an expansion around a free (massive or massless) theory, which is a common choice in modern implementations. By construction, this method could in principle be applied in dimensions D > 1 + 1.
We benchmarked the method using the massive free boson. Despite being a free theory, it represents a nontrivial test of the method because its Hamiltonian is a non-diagonal perturbation of the massless free theory. The exact solutions for this model can be obtained using covariance matrix methods [14] making it suitable as a benchmarking model. We found excellent agreement of the von Neumann entropy with theoretical predictions for ground and thermal states as well as for dynamics after quenches. We have demonstrated that the method is capable of a reliable time evolution up to times several times longer than the system size.
We proceeded by studying an interacting system, the sine-Gordon field theory in the attractive regime. For the scaling of the ground state von Neumann entropy, we found sG ground states to be much more long-range entangled than KG ground states, exhibiting a logarithmic scaling. In the large volume regime, the vN entropy depends only on the gap of the system but not the higher particle content. Studying the quench dynamics of the sG model, we found an oscillating behavior, as predicted by [62]. The resonances in the frequency spectrum of the time series matched the masses of the lowest breather states even under charge conjugation. The questions whether the oscillations are damped and superposed with a linear growth remains unanswered. In order to study that, a more sophisitcated implementation of the sG model would be required that would allow for higher cutoffs. The recently developed chirally factorised approach [26] could be a suitable candidate.
Our method opens the doors to many interesting explorations and can be extended in several directions. It would be interesting to explore further the oscillatory time dependence of the vN entropy dynamics following quenches. Here, a possible direction could be to explore the role of integrability and the effects of integrability breaking. To do that, non-integrable perturbations of the sG model could be considered, like the double sG model and the massive sG model. Another more fundamental possibility would be the φ 4 theory which is a canonical non-integrable QFT model and has already been successfully implemented in the HT framework [22,23]. Our method can be adapted for the φ 4 model with a straightforward step of computing the Bogoliubov coefficients for the massive field HT expansion.
Furthermore, it would be interesting to study the entanglement Hamiltonian and the Bisognano-Wichmann theorem [63,64]. Several interesting properties have been established for the CFT case [65][66][67][68] and it would be important to explore how they extend to the interacting gapped QFT [69,70]. The explicit representation of the reduced density matrix in a computational basis makes our method naturally suited to such a study.
Much attention has been recently devoted to symmetry resolved entanglement (see [71][72][73][74][75][76][77][78]). It would be an interesting extension of our method to make it sensitive to the symmetry charge of the subintervals and thus resolve the entanglement per sectors.
An implementation of our method in D = 2 + 1 would be a very interesting step because of the lack of methods in dimensions higher than D = 1 + 1. By construction, our method can be easily generalised to any dimension. The main obstacle would be the quickly growing size of the full and the split system Hilbert spaces. However, HT has already been successfully applied in higher dimensions [41] and since our dimension of the split Hilbert space in practice does not exceed the dimension of the full Hilbert space, such an undertaking seems possible.
Finally, in case of free theories, massless and massive, our construction yields an exact construction of the reduced density matrix of the theory. It could be a fruitful direction to use that to get further analytical insights into the entanglement structure of QFT.
Appendix A: Cut-off effects and symplectic structure The main approximation in our method is the representation of the full system modes in terms of a finite number of partial modes. We have to ensure that the approximation conserves basic properties of the system like the bosonic commutation relations. Equation (22) can be reformulated to test the transformation as Equation (A1) evaluates the commutation relations of the full modes A k expressed in terms of the partial modes. The structure of K on the diagonal of 1 on the first half of the diagonal and −1 on the second half reflects the anti-symmetric nature of the commutator upon exchanging its arguments. The commutation relations of the reverse transformation, partial modes expressed in full modes, are tested in (A2). The two equations provide us with an objective quality criterion of our truncated method. If the commutation relations of the bosonic modes are not fulfilled, the transformation is invalid. As described in the main text, we consider two cut-off schemes. The fixed cut-off scheme distributes the partial modes symmetrically across both intervals (s L = s R = s F 2 ). The second scheme, constant mode density, distributes the modes proportionally to the size of the intervals (s L = L s F , s R = s F − s L ). The quality of the two schemes can be assessed by checking the commutation relations of the transformed modes. Figure 10 shows the result of the calculation of (A1). The top row shows that the fixed cut-off scheme does not reproduce the bosonic commutation relations if the full modes are expressed in terms of partial modes for splits that are not at x = 0.5. A cut in the middle represents a special case. Here, the constant mode density cut-off scheme and the fixed cut-off scheme coincide. The bottom row illustrates that the constant mode density cut-off scheme reproduces the correct commutation relations for different cuts. All the computations in the main text are performed with this cut-off scheme. The constant mode density can be realised exactly only for a finite number of splittings . If we assume s F full modes, we can split the system at multiples of 1/s F , that is ∈ {1/s F , 2/s F , . . . , L − 1/s F } such that we have the exact same density of modes on the left and the right side of the cut (s L / = s R /(L− )). We call those values commensurate cuts. For other points, the mode densities cannot be chosen to be the same on the two partitions. A possible choice could be to pick a rounding scheme for the distribution of the partial cutoffs, for example s L = round( L s F ). Unfortunately, this leads to imperfect realisation of the symplectic structure (A2).
As expected the errors in the symplectic structure in non-commensurate points lead to incorrect values of the vN entropy. We show in Figure 11 the influence of different rounding schemes. The points in blue floor the number of modes in the left partition s L = L s F . This leads to increasingly bad results as we move to the right between commensurate cuts. If we round the number of left modes instead (s L = round( L s F ), depicted in orange), the problems get less severe and obtain a symmetric structure around the middle of the intervals. The method only gives correct results for commensurate cuts of the system which are drawn in green in Figure 11. All results presented in the main text are computed at commensurate splittings. For these points, the equilibrium state results for the Klein-Gordon model converge to the results predictions already at very modest cutoffs s F .

Appendix B: γ coefficients for Dirichlet boundary conditions at
In this Appendix, we first show into more detail how to recover (15) in the main text. For completeness we also outline the γ coefficients for the case of Dirichlet boundary conditions at the cut (φ L ( ) = φ R ( ) = 0). Via the Bogoliubov transform, eq. (20) in the main text, they relate full system modes A k in terms of the partial modes a R m and a L m The starting point is the mode expansion of the scalar field of the full system with p k = k π L and [A k , A l ] = A † k , A † l = 0 and A k , A † l = δ k,l . The expression can be inverted with the help of the canonical conjugate momentum field By taking a linear combination of the scalar field (B1) and its momentum (B2), we obtain for each k and t = 0: Projecting out by multiplying the expression on both sides with sin (p k x) and integrating completes the inversion of the mode expansion of the field Our aim is to express the field operator on the full interval in (B1) by the fields defined on the sub-intervals. For Dirichlet boundary conditions at the cut (φ L ( ) = φ R ( ) = 0), their mode expansion is given by: where we have defined p The relationship between the full fields (B1) and split fields (B5-B5) is given by the continuity condition, eq. (14) in the main text. Plugging it together with the split field expansions (B5-B6) into (B4) and performing the integrals gives the desired Bogoliubov transformation between the full and the split modes, eq. (20) in the main text.
The resulting γ coefficients for Dirichlet boundary conditions at the cut are The case of Neumann boundary conditions at the cut (∂ x φ L ( ) = ∂ x φ R ( ) = 0) which was used for most of the results presented in this work is outlined in the main text. The corresponding γ coefficients are given in eqs. (16)(17)(18)(19).
Appendix C: Algorithm

Derivation and Derivatives of Commutators
In order to use the scalar formulation for the matrix elements of U T in (40), we have to bring the terms into normal order. Our starting point is (31) a ξ, † m ρ ξ,λ m,n a λ, † n .
To avoid unnecessary jumping back and forth between the main text and the appendix, we will repeat some of the equations here. As mentioned in the main text, S creates the excitations of the split modes on the partial vacuum according to the occupation numbers in | n L , n R . F represents the creation operators of the full modes according to | n F expressed in the split modes. Finally, V transforms the full vacuum into a squeezed state on top of the split vacuum.
We normal order the expression in three steps. Firstly, we normal order the exponential e F which contains both creation and annihilation operators. Then, we commute e S , which consists of annihilation operators only, past the creation operators of e F . Finally, we commute all annihilation operators of e S and e F past the vacuum transformation e V .
Normal ordering or e F is achieved by the application of the Baker-Campbell-Hausdorff formula e X e Y = e X+Y + 1 where we have defined to be the parts containing creation/annihilation operators respectively. The commutator in (C1) evaluates to In a second step, we commute annihilation operators of partial modes in e S past e F + , using BCH in the form with the commutator Finally, we commute the annihilation operators in e S+F − and F − past the vacuum transformation e V . For convenience, we denote the annihilation operators in S and F − as A The computation of the matrix elements of U T in (40) does not depend on the form of T directly, but on the second derivatives T [J i , J k ]. All terms that are not derived twice will vanish once we set J k = 0. The second derivatives are: ComF, given in (C2): and other derivatives vanish. ComSF, given in (C4): and other derivatives vanish. ComAV, given in (C8):

Tree building algorithm
The primed sum in eq. (44) runs over all lexicographically unique configurations of the arguments (J 1 , J 2 ) of T . Lexicographically unique implies that all tuples are sorted internally J 1 < J 2 and the string of tuples is sorted as well. Two tuples are sorted by sorting them first by their first entry and then by second entry.
The generation of all unique pairs can be approached in at least two ways. We could take all J s independently and find all possible ways of distributing them as pairs. The available J s are determined by the occupation numbers in the states that are determining the matrix elements in the unitary transformation matrix. If the occupation number vector is n = (2, 2, 2) (as in the example), the available J s are (J 1 , J 1 , J 2 , J 2 , J 3 , J 3 ), i.e. we derive each twice with respect to each J . In the case of independent generation of all combinations, we would have to filter out all the repeated configurations due to ordering in the arguments of T and in the string.
Alternatively, we can incrementally create all orderings in a tree-like structure. By tracking the ordering as we progress, we can avoid the generation of forbidden configurations. We will only consider this second alternative since the generation of all permutations scales with n! where n = i n i , i.e. the number derivatives and the task of computing all combinations is unnecessary.
We start with vector n = (n 1 , n 2 , . . . , r N ). Since we consider only sorted tuples and a globally sorted string of tuples, we build all valid pairs (J i , J k ) with the first J i corresponding to the smallest non-zero n i . We proceed in a recursive manner and modify the vector n by subtracting one from n i and n k and select again all valid pairs in the next step of the tree. A pair is valid if the J i < J k for a pair (J i , J k ) and the pair is greater or equal to the previous selected pair.
In total, we build a tree of the form in figure 12. We start on the left with the full string of the example that is also used in the main text n = (2, 2, 2). Lists in round brackets describe the occupation numbers of the state. precisely constitute all the possible partitions of the set ν ∪ e i and the expression can be put in the exactly same form as (C25). This is written out in the last equality of the expression above.
Going back to evaluating the derivatives of the generating functional, eq. (C24), we can now use the above derived formula (C25) for f = exp(•) and g = T (J 1 , J 2 , . . . , J N ). Using the quadratic form of T , resulting in the derivative rules in eq. (C19), the nonzero contributions are only going to come from those partitions of ν where every subset contains exactly two elements: where we have copied the second line of eq. (C24) for comparison. It now remains to count how many partitions in the sum in the top line correspond to the same lexicographically ordered index set in the bottom line. Since on the partition side we distinguish singletons e However, not all the partitions that we obtain with this procedure are distinct. Because the partitions are insensitive to the order of elements in subsets ν j , we have over counted by a factor of 2 every diagonal ν j . These are subsets . This is indeed a manifestation of the commutativity of derivatives. We thus overall over counted the number of partitions by a factor of where N diag,k is the number of diagonal ν j (diagonal bi-indices (i, i)) in the index set I k . Furthermore, partitions are insensitive to the order of the subsets ν i . We have therefore over counted by the factor p kl ! for all those bi-indices (i l1 , i l2 ) that repeat p kl number of times in the index set I k . This gives an overall factor of l (p kl !) (C30) of overcounted partitions. Combining all these together, the number of distinct partitions that each lexicographically ordered index set I k corresponds to is This is the sought for combinatorial multiplicity.

Appendix D: Hamiltonian truncation matrix elements
We use the following matrix elements from the methods developed in [21,30] to perform the HT calculations of ground, thermal and non-equilibrium states.
For an operator O, we list matrix elements They are computed between states spanning the computational basis, the Hilbert space of the massless free boson: Here, A k , k = 1, 2, . . . are the bosonic modes fulfilling the canonical commutation relations [A k , A l ] = A † k , A † l = 0 and A k , A † l = δ k,l , the normalisation is N n = k>0 √ n k ! and A k |0 = 0 ∀k is the vacuum of the massless free boson theory.
The massless free boson Hamiltonian for Dirichlet boundary conditions is diagonal with matrix elements H n , n 0FB = π L ∞ k=1 kn k − 1 24 δ n , n .
The Hamiltonian of the massive free boson has the following matrix elements H n , n mFB = π L δ n , n ∞ k=1 1 + m 2 L 2 2π 2 k 2 kn k − j =k δ n j ,nj    1 k 2 n k k (n k − 1)k δ n k +2,n k + (n k + 2)k (n k + 1)k δ n k −2,n k The Hamiltonian of the sine-Gordon model can be expressed as Here V p (x) ≡ e iqφ(x) , p ∈ Z (D8) for q ≡ pβ is the vertex operator, M is the semi-classical soliton mass, the interaction related coefficient ∆ is defined as ∆ ≡ β 2 8π and the coupling-mass ratio κ(∆) [79] is The vertex operator can be written in normal ordered form as V p (z,z) = e iqφ(z,z) = |z −z| −q 2 /(4π) : e iqφ(z,z) : where for convenience we have introduced z ≡ e i π L x . The matrix elements are with the Heaviside step function Θ(•).
Appendix E: Covariance matrix formalism for entanglement entropies of free theories Eigen, thermal and time-evolved states of free theories are Gaussian and the computation of their entanglement properties can be simply achieved by means of the covariance matrix formalism [14].
A Gaussian state is completely determined by its covariance matrix where [φ m , φ n ] = [π m , π n ] = 0, [φ m , π n ] = iδ m,n are harmonic oscillator conjugate pairs. All higher order correlations are given by Wick's theorem.
In the case of a bosonic field theory, we have to introduce an IR (finite volume L) and UV (maximal momentum mode kept, s F ) cutoff in order to keep the covariance matrix finite. Then, the harmonic oscillators are finitely many and we can treat the covariance matrix either in momentum space or position space. For covariance matrix calculations it is convenient to go back and forth between those representations using a discrete sine transform. The field expansion can be written as with p k = kπ L and x n = na for k, n = 1, . . . , s F with the lattice spacing a = L s F +1 . The inverse discrete sine transform is achieved by π(x n , t) sin (p k x n ) .
Such definitions of fields correspond to approximating the field theory with a lattice system. We shall be keeping the relativistic dispersion, however. For the calculation of entanglement entropy it is most convenient to take the position space covariance matrix of the bosonic degrees multiplied by the lattice spacing a, that is φ n = aφ(x n ) and π n = aπ(x n ) to keep the correct dimensions. In position space, the covariance matrix of a reduced density matrix corresponding to a subsystem is achieved by taking only those matrix elements corresponding to the lattice points that lie within the subsystem. For example if we are interested in the entanglement between the interval [0, ] and the rest of the system, we take the covariance matrix of the lattice sites x n ∈ [0, ].
Then, the von Neumann entanglement entropy is computed by calculating the symplectic spectrum of the covariance matrix Γ. This is achieved by diagonalising iJΓ (E5) with the symplectic unit The eigenvalues appear in pairs ±γ k , k = 1, . . . , s F . This maps the problem to computing the entropy of s F harmonic oscillators at inverse temperatures The von Neumann entropy is then and the Rényi entropies are The covariance matrix approach to computing the entanglement entropies is convenient because it lets us model also the truncation effects by taking finite s F . The results are not exactly comparable to the HT cutoff at the same maximal momentum, because in the HT case, we have an energy cutoff which implies also a maximal occupation number for a mode. But it is the closest approximation of the cutoff effect that we can get. Taking s F large, we can recover the exact analytical results in the continuum limit. π k (t) = − k sin ( k t) φ k (0) + cos ( k t) π k (0) .
The procedure for a KG mass quench from the prequench mass m 0 to the postquench mass m is the following: take the momentum space representation of the thermal correlations (E11) for the prequench mass m 0 and propagate them using the equations of motion (E13) for the postquench mass m. The covariance matrix is transformed back to position space, the reduced covariance matrix corresponding to the subsystem taken as described above and the entanglement entropies computed.