Notes on the Causal Structure in a Tensor Network

In this note we attempt to understand Lorentzian tensor networks, as a preparation for constructing tensor networks that can describe more exotic back- grounds such as black holes. We first compare tensor networks with the framework of algebraic quantum field theory (AQFT), and find that it is a natural arena for AQFT. We then construct simple examples of Lorentzian tensor networks, based on Gaussian models of fermions and integrable models. We discuss how salient features, such as the Unruh effect, well known in QFT in Minkowski space-times can be captured by these skeleton of spacetime, and find physical understanding of the Reeh-Schlieder theorem via these simple constructions. We find interesting subtleties when the same approach is directly applied to bosons-the operator algebra contains commutators that take the wrong sign-resembling bosons behind horizons.

1 Introduction relativistic QFT using operator algebra. This sets the stage for precisely identifying notions such as Cauchy surfaces, causal diamonds etc in a tensor network. Then in the second part of the paper, we study explicit models. We start with a quadratic fermionic model, and demonstrate, at least in some limits very explicitly, that notions such as Lorentz transformation and the Unruh effect can be described to very good approximation in the tensor network. Then we consider more general integrable models, and found that our notion of a boost operator is a close relative of the so called "corner transfer matrix" that has been considered in the literature of classical integrable lattice models.
This hopefully lays out some basic features of tensor networks describing Lorentzian space-time, and serves as preparation as we move on to a covariant construction of more general interesting spacetimes based on the tensor networks.

Axioms of AQFT
In this section we give a very brief summary of the basic axioms defining an AQFT, listing all the ingredients that are going to have a natural realization in a tensor network. Our discussion is heavily based on [5] which gives a relatively gentle review of the subject. We also find [6] a concise and physical exposition of the subject. The Kaag-Hastler axioms are motivated by incorporating locality and causality into an operator algebra that in turn defines a QFT.
The ingredients involved are therefore spacetime manifold M on the one hand, and some operator algebra A on the other. There are various conditions imposed on M. Of course traditionally, the discussion is based on smooth manifolds. We will list all the ingredients and discuss which can make direct contact with the tensor network. First, we need some constraints on the set of spacetimes concerned.
• M as a topological space is Hausdorff, connected and paracompact.
• M has a Riemannian metric g. This defines a causal structure. For a smooth curve γ(t) in M, it can either be space-like if g(γ,γ) < 0; causal (time-like or null like) if g(γ,γ) ≥ 0 • M is globally hyperbolic such that it does not contain closed causal curves and for any two points x and y, J + (x) ∩ J − (y) is compact. Here J + (x) denotes the collection of points that are in the "future" of x and that can be connected to x by a future directed causal curve. J − (x) corresponds to the causal past of x.
• Having a causal structure, it also means that we can define an order relation x ≺ y i.e. x precedes y if there exists a future directing causal curve connecting x to y. This relation can be generalized to subsets in M. Consider two sets O 1 and O 2 . They are space-like separated if they cannot be connected by a casual curve -if for all x ∈Ō 1 , J ± (x) has empty intersection with O 2 .
• There exist foliations by Cauchy surfaces, (locally) diffeomorphic to Σ × R • There exist "admissible embeddings" χ : M → N for globally hyperbolic spacetimes M, N , such that the map preserves the metric, orientations and causal structure.
• Collection B of bounded subsets of M forms a directed set. There exists a reflexive and transitive binary relation : for a pair O 1 , O 2 ∈ B, ∃O : These spacetimes form a "category" LOC -category of locally hyperbolic, oriented and time-oriented spacetimes.
With such a space-time as a starting point, the AQFT is a rule of assigning an operator algebra to hyperbolic spacetimes. Mathematically, an AQFT is a functor that maps between two categories. LOC on one hand, and OBS -the category of unital C* algebra corresponding to the algebra of operators of physical observableson the other.
In QFT, it is taken that observables form a C * -algebra. For completeness, we provide the definition of C * -algebra below. A C * -algebra A comes with a norm ||x|| and a map * that takes A → A. For x ∈ A, where A (λx) * =λx * (2.3) ||x * x|| = ||x||||x * ||. (2.4) The "C" stood for "(norm)-closed". The norm provides the algebra a topological structure. Let us make contact with quantum mechanics. Consider H to be a complex Hilbert space with inner product denoted ·, · . The collection of bounded linear operators on H, denoted by B(H), is a C * -algebra. The linear structure is clear. The product is by composition of operators. The * operation is the adjoint; for any operator a on H, its adjoint is defined by the equation a * ζ, η = ζ, aη , for all ζ and η in H. Finally, the norm is given by ||a|| = sup{||aζ|||ζ ∈ H, |ζ| ≤ 1}, for any a in B(H).
The axioms constraining the AQFT functor U are as follows: U assigns to each bounded subset O ∈ B a C * -algebra U(O). The algebra of spacetime M is defined as the inductive limit • time slice axiom. The solvability of the initial value problem is translated into a requirement of the operator algebra. The algebra U(N ) is isomorphic to U(M) for any causally convex neighbourhood of a Cauchy surface Σ. Note that a causally convex neighbourhood N is one in which no causal curve inside M meets N in a disconnected set.
• In a generally curved spacetime, it is more suitable to invoke local charts and they can be understood as the admissible maps χ discussed above. For each admissible embedding χ : N → M , there is an injective homomorphism α χ : U(N ) → U(M ). If χ 1 : M → N and χ 2 : N → L then we have i.e. altogether, the assignment of algebra via U has to be compatible with the structures of embeddings leading to maps between algebras. That makes an AQFT a covariant functor between LOC and OBS. In terms of these embedding maps, the Einstein Causality constraint can be phrased as follows: if χ 1 (M 1 ) ⊂ M and χ 2 (M 2 ) ⊂ M are causally disconnected, then [α χ 1 (U(M 1 )), α χ 2 (U(M 2 ))] = 0 (2.7) Similarly, the time-slice axiom can be restated if we introduce an embedding map χ : N → M , where N is the causally convex neighbourhood of a Cauchy surface, we have α χ an isomorphism.
• The above isomorphism describes general covariance. If we only have global isometries, then the statement is instead restricted to the existence of an isomorphism α L that maps α L : U(O) → U(LO), where LO is the region O transformed via some isometry L (such as Poincare symmetry in Minkowski spacetime).
• time evolution between any two Cauchy surfaces: given two Cauchy surfaces Σ 1,2 , since the neighbourhood N i of each can be associated an embedding map χ i : N i → M such that α χ i is an algebra isomorphism from U(N i ) → U(M ).

Realization of these objects in a tensor network
We would like to associate building blocks of a tensor network to ingredients discussed above.
A tensor network is a collection of tensors contracted with each other, with each tensor represented as a vertex and shared contracted indices between tensors represented by an edge connecting these vertices. Therefore in this simplest version, the tensor network associates a graph to a collection of tensors.

Topological space
A graph F can be viewed as a topological space naturally endowed with 0 simplices (vertices) and 1-simplices (edges). This topological space is the arena that plays the role of the spacetime manifold M . Clearly there are more structures to spacetimes than such a skeleton. We can follow which structure is directly made use of in the tensor network.
As a set of vertices and edges, the notion of O i ⊂ O j , where O i and O j are subgraphs of F is well defined. Therefore, this is also a directed set.

Hilbert space and local operator algebra
Each edge is associated to a contracted index of a tensor. Therefore each edge can be naturally associated with a Hilbert space. (Perhaps practically finite dimensional.) Operators acting on this Hilbert space would thus form a C* algebra, as discussed in the previous section. Associating an operator algebra to links do not constitute U(O). These operators are related to each other. Such relations will be laid out below.

Causal Structure-timelike vs spacelike separation
The graph does not generically possess any causal structure. The associated tensor network without any causal structure might well be associated to a Euclidean version of spacetime. To describe Lorentzian signature, the graph needs to acquire a causal structure.
One way of incorporating causal structure in the tensor network is based on local unitaries.
Consider the special case in which each tensor T has an even number of legs, each with the same bond dimension. We can split the indices into two equal groups g 1,2 . If T β j ∈g 2 α i ∈g 1 is a unitary from α i to β j , then we would consider the vertices connected via these edges to be causally connected. We can also assign a causal ordering.
Such an ordering can only be made consistent with each other however as follows. Suppose two indices α 1 ∈ g 1 and β 1 ∈ g 2 , where we have assigned the ordering g 1 ≺ g 2 . If there exists another split into two groups g 3 and g 4 such that α 1 and β 1 belong to g 3 , then at least one index in g 1 must now belong to g 4 and at least another one in g 2 now belong to g 3 . In such a case, requiring either g 3 ≺ g 4 or vice versa would be inconsistent with g 1 ≺ g 2 . Therefore a causal ordering assignment can be made only if there is a unique splitting g 1 , g 2 such that T is unitary. In which case, edges within the same group can be considered space-like separated. i.e. Each unitary evolution can be understood as defining a set of observers, or a frame.
For T being perfect tensors for example, then every pair of indices are time-like separated, and yet there is no ordering agreed by all frames.
Such an assignment is local. The above assignment makes it natural to include arrows in the edges to denote causal precedence. We can put in-going arrows among edges in g 1 and out-going arrows in edges in g 2 . As a unitary matrix, the number of arrows are conserved across each vertex.
For a global assignment of causal structure, one needs to pay special attention to how edges are contracted with each other. With the arrow assignment, we only need to ensure that when assembled together, an out-going edge proceeds to become an in-going edge in the other vertex it connects to.
This requirement is equivalent to the fact that the graph F is orientable. The requirement of the absence of closed time-like curves becomes the requirement that the graph is a directed acyclic graph. The causal structure is borne out by the building block being local unitaries. This is similar to the consideration in the causal-set approach to quantum gravity. (This is a huge subject in its own right. We refer interested readers to one of the original papers [7] and [8] and references therein for more recent discussions.) In the context of the tensor networks it is not clear whether graphs that are transitively closed/complete -graphs where there exists an edge connecting any 2 causally related points-play any special role.
To summarise, it appears that the conservation of in-going and out-going legs on a directed acyclic graph has replaced the notion of global hyperbolicity.

Operator pushing and local unitaries
The local unitaries T defines isomorphisms between operator algebras. Operators from in-going legs are related to operators in the out-going legs by conjugation (or equivalently operator pushing). An immediate consequence is that as we push an operator starting from a Cauchy slice across T 's, we sweep out a light-cone, either forward or backward in time, a feature already observed for example in the context of Multi-scale entanglement renormalization ansatz (MERA) tensor networks [9]. This ensures that information cannot flow faster than the speed of light, which is captured by the Einstein Causality condition below.

Operator algebra assigned to a connected subgraph O
Now consider O to be a connected subgraph of F . Then we can pick out a set of edges all mutually space-like separated. The operator algebra U(O) associated to this region O can be defined as the operator algebra on the maximal set S of mutually space-like separated edges. Operators acting on any other legs can be pulled back to operators acting on S via local unitaries. Or in other words, we will include in this operator algebra only those that can be pulled back to S. This defines U(O) that is generated by the operator algebra on S. To make subsequent discussion simple, when we discuss these subregions O, they should carry the structure of a causal diamond -the boundary of O should intersect ∂S. The boundary is thus separated into two pieces by ∂S. Each piece is related to S by a unitary map. Observation on the tensor network: if O 1 ⊂ O 2 , then S 1 ⊂ S 2 and thus U(O 1 ) ⊂ U(O 2 ). The set of algebra {U(O)} also forms a net of (C*) algebra.

Cauchy surface and the time slice axiom
For an acyclic directed graph, it admits a slicing which is a set of edges, such that not any two are related causally, and that the surface does not have a boundary. Such a surface can generically be chosen as a Cauchy slice Σ. A typical Cauchy surface on the network is illustrated in figure 2. A neighbourhood N of a Cauchy surface Σ can now be defined as picking up new edges connected to edges on the Cauchy surface by the tensors T . Since individual T are all local unitaries, this defines an algebra U(N ) which is isomorphic to the algebra generated by the C * -algebra on Σ via the unitary maps T . This is also isomorphic to the algebra U(F ), which is isomorphic to the C * -algebra at any Cauchy surface Σ.
Algebra on different Cauchy surfaces are related also by isomorphism. The isomorphic map α 1→2 : U(Σ 1 ) → U(Σ 2 ) is now provided by the sets of T contained between two Cauchy surfaces Σ 1,2 and this defines a unitary time evolution.
That any two choices of Cauchy surfaces can be related by a unitary map is probably equivalent to the Stone-von Neumann theorem in finite dimensional Hilbert spaces. (See for example [6] for an explanation of the Stone-von Neumann theorem. ) Figure 3 illustrates a unitary transformation between the horizontal surface and the slanted ladder-like surface. As we are going to see, in a homogenous network where every tensor is the same, this can be interpreted as an approximate Lorentz transformation. Figure 3. The set of network between the two Cauchy surfaces : the horizontal Cauchy surface marked black, and the ladder like Cauchy surface marked yellow, is a unitary matrix, and we define that as the discrete approximation of a Lorentz Transformation.
A set of observer defining a "frame" is characterized by an ordered set of Cauchy surfaces, and these Cauchy surfaces are related by unitary transformations. If there is time translation invariance, it corresponds to the fact that the unitary matrices that map one Cauchy surface to the next remain unchanged.

Einstein Locality
To make actual comparison between two operator algebra U(O 1,2 ) assigned to two different connected subgraphs O 1,2 , we embed them into U(Σ) for any choice of Cauchy surface Σ, which defines an embedding map α 1,2→F : With that, we can define Einstein Locality between two space-like separated regions O 1,2 using α 1,2→F : As noted above, these maps α i of operators are basically operator pushing across tensors T along unitary directions. The above commutation relations can be phrased equivalently as follows. Consider a Cauchy surface Σ 1 containing S 1 . Now pull the operators in S 2 to Σ 1 which we denote by the map α 2→Σ 1 . If α 2→Σ 1 (U(O 2 )) ⊂ U(P 1 ), where P 1 ⊂ Σ 1 and S 1 ∩ P 1 = 0, then (3.1) is satisfied.
Thus far, the tensor network falls short of being a "functor" mapping the category of graphs to OBS. The reason is that it is not obvious what is the physical data that goes into defining a functor that maps different graphs to different OBS that can be compared with a quantum field theory. 1 Note that to explicitly compute these commutators, we first construct Cauchy surfaces that contain each of these locations and then perform operator pushing of one operator in one of the Cauchy surfaces to the other surface. The precise choice of the Cauchy surface is immaterial, but existence of which is crucial. That would essentially rule out many potential confusion over whether two points are in fact connected causally.

Isotony
Isotony can be defined as follows in the tensor network. We will restrict our attention to subregions in spacetime describable by causal diamonds. The algebra attached to a causal diamond can be described as follows. We locate the space-like surface and the intersection of its causal future and causal past to define a causal diamond. The operator algebra associated to this region of space is defined as the operator algebra defined along the maximal space-like surface inside the causal diamond. Any other space-like surface are thus related to each other by isomorphisms. Now if a causal diamond is inside another causal diamond, it is always possible to map the algebra on the space-like surface in the smaller diamond to the algebra on the maximal space-like surface in the larger diamond. The operator algebra is sure to be a sub-algebra of the algebra defined on the larger causal diamond. Therefore, an algebra net is naturally recovered for causal diamonds.
Note that for a space-like region A with a causal diamond D(A), it follows that the entanglement entropy only depends on D(A) but not on the specific maximal 1 When defining a quantum field theory, a standard procedure is to define a Lagrangian for some given set of fields. The program based on category theory defines a quantum field theory without using a Lagrangian. For a CFT, one needs a set of primaries, their conformal dimensions and their OPE's to completely specify the CFT. It is not completely clear to the authors what is the full set of data needed to specify a generic QFT in this language. A tensor network probably has more data than is necessary. space-like surface A chosen inside D(A). This is because all these maximal space-like surfaces inside D(A) are all related by local unitaries which preserve entanglement entropy.

Illustration based on free fermions
In this section we would like to illustrate some of these ideas using a simple model. We would like to construct a unitary evolution that is explicitly expressed as a tensor network of local unitaries. Then we would like to construct notions familiar in a continuous space-time in the present context, and show that they can be approximated to some extent.

Tensor network evolution
First of all, we have to define a model that is inspired by the free fermions but whose time evolution takes the form of a tensor network built from local unitaries. In fact, breaking up a generic unitary evolution into a product of local unitaries is frequently employed in numerical simulations or actual experiments on quantum simulations [10,11]. The Lieb-Robinson bound is also based on such an approximation [12]. Here, however, we will take the view point that the tensor network defines the model.
Our tensor network is constructed as follows. Consider a set of fermion creation and annihilation operators a i , a † i , where i denotes the link where the fermion is located. These operators satisfy the usual anti-commutator.
To construct the simplest example of a unitary evolution, we consider the quadratic Hamiltonian, where 2L is the size of the spin chain. We have already taken lattice spacing a = 1 here.
To ensure that H is Hermitian, we take α =α 2i for some realα. The time evolution operator U (∆t) over the smallest unit of time ∆t is given by where We have, U a 2x U −1 = c 2 a 2x + c s a 2x−1 − c s a 2x+1 + s 2 a 2x+2 , (4.5) and where c = cos(iα∆t), s = sin(iα∆t), and x is the label of "site" on a given Cauchy surface -the vertex to which the links are connected. For the sake of convenience in later (numerical) computations of correlation functions, we have labeled these vertices by half-integers x, and the links by integers. We therefore have 2x + 1 labelling "even" links, and 2x labelling "odd" links. The labelling is illustrated also in figure 4. Time evolution is given by repeated application of U . There is thus a time translation invariance over t → t + ∆t. The simplest set of Cauchy surfaces are those that are "horizontal". The tensor network and the labelling are illustrated in figure 4. . The collection of horizontal Cauthy surfaces, defining a particular set of "inertial observers". We note that we adopt a rather odd labelling for convenience with numerics, and take "site" numbers as half-integers x, while links are labelled by integers 2x, 2x + 1.
The network is invariant under translation by two links i → i + 2. Therefore based on translation invariance, we expect that the eigen-modes should be given by plane waves.

Spectra of the model
From the translation symmetry of the network, we expect the eigen-operators to be expressible as for some coefficients f 2x (p) and g 2x+1 (p) given by where q 1,2 are independent of the site number x. They are to be determined by the following eigenvalue equations. The momentum is given by where n are integers satisfying −L ≤ n ≤ L − 1. In the limit L → +∞, we can take −π < p < π. The eigenvalue equation is obtained by demanding This implies For each given p, it gives two eigenmodes with eigenvalues We note that the eigenvalues are even in the momentum p. In the limit c → 0 the eigenvalue approaches lim c→0 λ ± (p) = exp(±i|p|), (4.12) which means that energy is linearly dependent on p. This recovers the dispersion of a relativistic free massless fermion. In the other limit where c → 1, we have λ(p) approaching a constant. The model becomes non-dispersive. Therefore while the graph on which the model is defined remains unchanged, the effective causal structure of the system changes dramatically as the unit of time ∆t in the model is varied. For λ − , the corresponding eigen-modes are given by, (4.13) Then for λ + the eigen-modes are, s sin(p) + (1 − cos(p)) (2c 2 + s 2 cos(p) + s 2 ) 2c sin p 2 (4.14) We can then define the operators corresponding to creation and annihilation of these eigen-modes.
We check that, which implies, Similarly for b p . This gives, On the other hand, the equation (4.10) implies which means that so b † p and a p behave as creation and annihilation operators respectively. We can now define the fermion field as where cos(χ) = c 2 + s 2 cos(p). Similarly, (4.23) We would like to define the vacuum state as, In a discrete spacetime the notion of "ground state" is clearly ill-defined. If the smallest unit of time is ∆t, then the energy is identified with a period i.e. E ∼ E + 2π/∆t. However, mimicking the continuous situation in search for a state that corresponds to a separating vector, and requiring that the p → 0 limit should recover our usual notion of low energy states, one can still identify a natural "vacuum" state, which we have picked above. We note that while a discrete evolution in Lorentzian signature does not identify a ground state without ambiguity, a Euclidean partition function would continue to project to a unique (unless otherwise degenerate) ground state. One wonders how such an analytic continuation can be defined if time were discrete. The folklore that a quantum model in d dimensions necessarily originate from a classical model in d+1 dimensions is by no means obvious when time becomes discrete. As we have seen above, variation of ∆t while keeping the Hamiltonian unchanged could lead to utterly different dispersion relations and subsequently causal structures. It appears that at least in the case of integrable models however, an analytic continuation between Lorentzian and Euclidean signatures remain well defined. We will postpone the discussion to the next section.

Lorentz transformation
Since the system has a well defined causal structure, it is very tempting to define the notion of boost -i.e. the set of (unitary) transformations between different observers. As already explained in the previous section, an observer is defined as a collection of Cauchy surfaces. The unitary transformation connecting observers are generically given by the set of tensors sandwiched between the respective Cauchy surfaces. In general, such a transformation does not preserve the "ground" state. For a general set of observers, there is not even any notion of time translation invariance, such that energy is completely ill-defined. For our model, there are different observers that enjoy some degree of time translation invariance. This is illustrated in figure 5 And if in the special case, that these different observers agree (perhaps only approximately) upon the notion of "positive energy" modes, it means that the vacuum state remains invariant under the corresponding transformation between these observers. These observers that (approximately) share the same ground state would be the closest analogue we have for "inertial observers" in a continuous Lorentz invariant spacetime.
In the current tensor network we are working with, there are natural families of Cauchy surfaces corresponding to inertial observers, inspired by results in continuous space-time. Each family of Cauchy surfaces are slanted surfaces with some given slope. We note however that in a discrete space-time, these slanted surfaces with given slopes are more accurately speaking ladders. For a given slope, the corresponding Cauchy "ladder" is not unique, and we understand them as slightly different approximations of these inertial observers which only become indistinguishable in the long wavelength limit.  Figure 5. This diagram illustrates two independent sets of "inertial observers", marked by sets of Cauchy surfaces with different colours. The tensors sandwiched between each pair of Cauchy surfaces of the same colour correspond to the corresponding "inertial time" evolution between the pair.
Boost transformation that relates these different observers are constructed below.

Approximating the boost operator
Consider the simplest scenario. There is a natural set of inertial observers with flat Cauchy surfaces. We consider another set of inertial observers over ladder surfaces. We would like to construct the unitary boost operator that transform the ladder into the flat surface. The most natural guess would correspond to the set of local unitaries sandwiched between the two Cauchy surfaces. However, after connecting the Cauchy surfaces, in general one can do further transformations that preserve the target Cauchy surface, such as translations on the surface. We therefore allow the boost operator to take the general form, corresponding to one that first translates along the time direction according to the sandwiched unitaries, which could be, depending on the situation at hand, followed by a set of translations along the target Cauchy surface. This is the same strategy taken by [13] in the Euclidean version when constructing transformations mapping between different (Cauchy) surfaces.
We can compute the effect these transformations have on our modes. To make the discussions most transparent and its similarity with the continuous case obvious, let us first consider the limit c → 0, in which case, we recall that the dispersion relations (4.11) reduce to a linear one. That is directly analogous to the discussion of the free fermions in the c → 0 limit in the main text. We note that in the limit c → 0, the fermion evolution (4.10) behaves like a (fermionic) SWAP, in which even links are translated to the left by two steps 2x + 1 → 2x − 1, and odd links are being translated to the right also by two steps 2x → 2x + 2. This suggests that the fermionic modes turn into a pair of chiral fermions propagating in opposite directions. In fact, one can directly check the expressions (4. 13-4.14) to confirm that modes with positive momenta has positive (negative) energies for the even (odd) sites. This is recapped below. We list the eigen-operators explicitly as c → 0. We note that the expression for the operator approaches a rather singular limit as c → 0 where special are has to be taken when taking limits of (4.13-4.14), the correlation functions are relatively straight forward to treat. Nonetheless, it is possible to show that the modes reduce to and then similarly The vacuum continues to be the one that is annihilated by a p and b p . The inverse transform would give In this case, we can inspect the effect the boost operator has on the operators a p and b p . The boost operator that we will illustrate in detail is shown in figure 6. The effect of the ladder operatorB has on a i is as follows. One can see that in For i = 0, the operator a 0 is invariant under the transformation of the ladder operator B. For i < 0, where In general, we could consider also the effect of a link-dependent translation that for some (odd) integer M which would preserve the current Cauchy surface. 2 But in the current illustration, let us do without further deformation of the Cauchy surface, which does not appear to alter the physics of the final result. Given the above expressions, one readily obtains the effect the boost operator has on the modes. For positive links, using (4.28), we obtain (4.33) Let us emphasize here that the apparent shift by w/4 above is an unfortunate result of our notation which takes x to be half-integers. The reader should be able to see from the figure 6 that the odd links particularly, are simply scaled. i.e. 1 → 3, 3 → 7 etc., with the steps scaling linearly with the link number. The even links behave similarly with a slight complication of some internal structure in the rescaling which is explained in equation (4.28) and easily recovered also from the same figure. Since our transformation only makes sense in the thermodynamic limit where the momentum also approaches a continuous limit, we replace the summation over p by an integral over the first Brillouin zone. Note that the lattice spacing is just a = 1 here. For negative links, using (4.30), we obtain .

(4.34)
To obtain the effect it has on individual a p≥0 , we take another inverse Fourier transform to get (4.35) The above modes live only on the "odd" links. Similar expressions can be obtained for b p<0 and a p<0 which depend on the "even" links. Although one can expect from (4.28) that the result would be less clean.
The interesting observation is that in the thermodynamic limit where the total number of sites approach infinity, the sum over sites gives us a delta function δ(q − where we have denoted Λ = (1 + w/2). We note that in the current example at hand Λ = 2. We thus recover the expected transformation of the modes under a Lorentz boost. For Λp > π, one would have mapped a positive energy mode to a negative energy one. In other words, under the operation of our guess of a "Boost", it (almost, up to some lattice effect which disappears as L → ∞ and |p| sufficiently small i.e. continuous limit) preserves the "positive" and "negative" energy modes, which is a necessary ingredient of a Lorentz transformation in a Lorentz invariant theory that preserves the ground state. Due to our choice of the form of the ladders, it does not treat the p > 0 and p < 0 modes in a completely symmetric manner, because the ladder breaks the symmetry between even and odd sites. However, we associate such deviations with lattice effects.
As we are going to see when we work with integrable models, where some exact results of the modular Hamiltonian are available, we reckon our guess turns out to be a very sensible one.

Spectra of the boost operator
In the previous subsection, we made a crude guess of the form of the boost operator. To proceed with understanding the Unruh effect, we would like to solve for the spectrum of the boost operator, and obtain eigen-modes in the x < 0 and x > 0 halves of the lattice.
Since the boost operator does not act on the link at i = 0, it essentially breaks up the space into two halves, the eigen-modes on the two halves of the lattice decouple. Therefore, the form of the eigen-modes can be considered separately on the right and left half of the lattice.
We will continue to work with c = 0 to illustrate the point, where even and odd sites remain decoupled.
We begin with solving for eigen-operators on the right side of the lattice. We first write down an ansatz for the eigenmodes. Since Then we require thatB where we are denoting σ p≥0 = a p , σ p<0 = −b † p to avoid clutter. Therefore the eigen wavefunction must take the form Strictly speaking, the sum over x leads to a delta function in the momenta which is only defined up to 2π. As we already noted, the boost operator only preserves the modes with long wavelength compared to the lattice spacing a. In the limit a → 0, the limits of the p integral would be replaced by ±π/a → ±∞. Solutions to the equation are then given by the following: Both leads to the same eigenvalue Note that the e ip/2 is from our fractional labels of sites x = 1 2 , 3 2 , · · · . (In the a → 0 limit this phase is really e ipa/2 → 1. We will however keep this factor at present which will not make any material difference. ) This is a solution in the rhs of the lattice where x > 0. We note that p can be both positive and negative here. Therefore, to recover the wavefunction in configuration space, we can in fact take an inverse transform, which gives Now, one can check that Here, it is necessary that −1 < Re κ < 0 while assuming x is real. This computation is not regulated based on giving a small imaginary part to x which would not have been a good regulator if p < 0. We would also like to check the result for ψ − κ (x > 0). It gives We can work out the eigen-modes in x < 0 in a similar manner. As we mentioned above, the x < 0 sector is the mirror image of the x > 0 sector. (See figure 6.) This immediately suggests that we have a set of eigen-modes given by where χ ξ (x) satisfies an entirely analogous set of eigenvalue equations (4.42). We would then end up with the solutions Now this recovers An important question then arises. What values should κ and ξ takes? Recall that we are taking these infinite ladders as a unitary evolution. Therefore, κ and ξ should be so chosen such that we have a pure phase.
We therefore would like to have (4.51) In other words, for positive would define positive energy modes.
Here for the χ ξ (x) modes, by comparison with the discussion in continuous field theory, would require that we associate for positive to correspond to positive energy modes. We interpret here that the unit time evolution is evolved backward, and so ∆t has to take a negative value, thus inverting the definition of positive energies. Here, we labor further upon some ambiguity that creeps in. There is not an exact translation symmetry with respect to the evolution between ladders of different slopes -the analogue of evolution of Rindler observers here. As a result, there is not a precise concept of conserved energy (within the "energy Brillouin zone") here. But nonetheless, there is an approximate symmetry, so that we can still define 0 ≤ 1 to be approximately conserved positive energy modes in this analogue "Rindler" frame. Correspondingly −1 ≤ 0 defines negative energy modes. That these "Rindler Hamiltonian" actually does approximate the entanglement Hamiltonian and thus bear actual resemblance to the continuous scenario with exact killing vectors would be discussed in our numerical computations and subsequently in the discussions of integrable models in the next section.
Summarizing, we have four sets of wavefunctions, two independent ones for the x > 0 region and another two for the x < 0 regions.
where they correspond to pairs of positive and negative wavefunctions in each region for each . We use these to define four sets of eigen-operators (4.55) The complex conjugates of each operator can be obtained from the above. Entirely analogous expressions, albeit looking less clean and simple, can be obtained for modes on the even sites, which we will not dwell on in detail here. Before we move on to a detailed discussion of the Unruh effect in the current context, we note that the boost operator constructed from an ever rising sequence of stairs are crucial towards the decoupling between left and right moving modes which can be thought of as arising from a "Rindler horizon" near the origin. Initially we attempted to approximate this effect with a stair with only 1 step. This is equivalent to solving a semi-infinite system with a fixed boundary condition at one end. One can show that there is always reflection at x = 0, leading to very different physics. This is discussed in detail in C.2 in the general context of integrable models.

Approximating the Unruh effect
Having worked out the eigenfunctions, we can follow a very similar path as the continuous case and look for the Bogoliubov transformation that connects eigenmodes of the "inertial observers" and the "Rindler observers".
Equation (4.42) is actually supplying the Bogoliubov transformations between modes in "inertial" frames and "Rindler" frames.
We would like to express a p in terms of A L,R ±i . We note that i.e. In the last equality, the mode correspond to ψ + i , but analytically extended in the region x < 0. This gives Here, one has to make a choice of the branch-cut. If we make a choice that is one would recover the standard result that the temperature for the half-space observer is given by 3T This can be compared with the standard result for relativistic field theories (in 1+1 dimensions) at vanishing mass. A detailed computation can be found for example in [14] where we have, where α R l are (positive energy bosonic) Rindler modes, and a k are inertial modes. The speed of light and acceleration has been set to 1. The solution for f l (k) is, up to normalization, precisely what we have found in (4.42). The fact that the modes are bosonic only changes the analysis concerning normalizations. Otherwise the relationships between the modes are a direct result of Fourier transforms, and the plane wave solutions are shared between the 1+1 dimensional fermions and bosons.
The even sites can be treated in a similar manner, although the result would not have looked as clean due to our breaking of symmetries between even and odd sites.

Comparison with half-space entanglement
What we have demonstrated in the previous subsection is that we can construct an analogue of a boost operator. Its effect on the modes in the inertial frame is very much analogous to the expected form in the continuous case -a p → a Λ −1 p , and this action (approximately) preserves the notion of "positive" and "negative modes" that we have defined in the discrete space-time.
Then we showed that eigen-modes of this boost operator have very similar forms to the Rindler modes in continuous field theories. We demonstrated this explicitly at special values of the parameters, and showed that they closely approximate modes solved in Rindler space, and subsequently define an analogous Bogoliubov transformation between "inertial" and "Rindler" modes.
What we have not demonstrated however, is that the boost operator that we have defined is indeed the same as the entanglement Hamiltonian controlling halfspace entanglement. In the case of continuous field theory with Lorentz invariance, this is guaranteed by the Bisognano-Wichmann theorem [15] ( see a review for example in [16]). The construction relies heavily on the analytic continuation of the boost operator that rotates t → −t and x → −x. As we are going to discuss in later sections inspecting more generic integrable models, such an analytic continued operator indeed exists and can be readily constructed -it is the Euclidean version of the "corner-transfer-matrix".
In this section however, we will compute explicitly the half space entanglement Hamiltonian in the fermionic model for different values of c, and demonstrate that the modular Hamiltonian does resemble the guess of a boost operator we made in the previous subsection, therefore adding extra evidence that the boost operator that we have constructed does indeed control the half-space entanglement. In that case, the evolution generated by our boost does correspond to evolution of observers that approximately have no access to half of the space-time, closely resembling Rindler observers.

The computation of the entanglement Hamiltonian
The entanglement Hamiltonian, which is also called modular Hamiltonian, is defined by taking the log of the reduced density matrix. In the case of free fermions, given that the vacuum state we have chosen has zero fermion number, the reduced density matrix takes the following form: where ρ V is the reduced density matrix of region V , H is the entanglement Hamiltonian of region V . The normalization constant is specified as K = det(1 + e −H ) −1 . In [17], the entanglement Hamiltonian is calculated by correlation functions of local creation and annihilation operators a † i , a j , which satisfy {a i , a † j } = δ ij . The correlation functions are defined as The entanglement Hamiltonian in matrix form is given as [17] where C is the matrix of correlation functions C ij in region V . In our analysis, we take region V as the links connected with sites labeled from 1/2 to L − 1/2, i.e. links from 1 to 2L. Here we take link a 1 as an even link, because a 1 is a link of a 2x for x = 1/2 and similarly, we take link a 2 as an odd link, because a 2 is a link of a 2x+1 for x = 1/2. We do the numerical calculation of entanglement Hamiltonian. We take L = 11 and take the parameter c = 0.000001, 1/ √ 2, 1 respectively. We plot them in the following figures 7, 8 and 9 for three different values of the parameter c. We explicitly checked that our entanglement Hamiltonian H is Hermitian, which is expected for a positive definite reduced density matrix. We can see that the entanglement Hamiltonian is quite local for different parameters, since H ij takes the form of localized peaks close to the diagonal of the plot, meaning that the two sites involved, i, j in a i and a † j are very closely spaced. The real part of H ij is much smaller than its imaginary part by about 10 times. Focussing on the imaginary parts, one can read off all the positions and heights of these dominant localized peaks.

Peaks of the Imaginary Part of Modular Hamiltonian
We consider the peaks of the imaginary part of modular hamiltonian of half space in the fermion model. In the following figures of this section, we take L = 11.
We note that in the case of very small c, H ij is non-vanishing only when j = i±2. This is because when c ∼ 0, the time evolution is basically decoupled between the even and odd sites. As a result, the two sublattices do not entangle. As c increases, the two sublattices begin to couple. At c = 1, H ij becomes non-vanishing for j = i±1. The entanglement Hamiltonian is dominated by these local contributions. There are some very much suppressed non-local terms involving i, j further apart. We would like to understand these dominant local terms.
The shape of the plots above suggests that H ij might carry a quadratic dependence on i. We fit it with a parabola ai 2 + bi + f . One can see that the fit works very well when c → 0 and c → 1, although in the former we have to separate the even and   the two sub-lattices are beginning to mix. For the imaginary part of modular hamiltonian with parameters c = s = 1/ √ 2, there are two set of peaks. One are located at |i − j| = 1. The other set are located at |i − j| = 2.  These plots suggests that the entanglement Hamiltonian takes the form for c → 0, and that for c → 1 for some overall normalizations N . For other values of c between 0 and 1, it is interpolating between (4.65) and (4.66). This operator in (4.65) or (4.66) have been discussed before in the literature. In the computation of entanglement entropy of lattice fermionic Gaussian states, it is observed that the entanglement Hamiltonian always commute with an operator of the above form [18,19]. Here, we are observing that it is itself a fair approximation of the entanglement Hamiltonian, at least in the c → 0 or c → 1 limit. The intermediate values of c has extra complications as the translation symmetry of the "ground state" changes. We also note that as c increases from 0 to 1, the lattice size has essentially doubled. But ν 1 and ν 0 moved a little. We take that as evidence that it does not scale with L.
When we are sufficiently close to the entanglement boundary, one could see that it is approximately taking the form for some constant σ and for n sufficiently small and close to the boundary of the entanglement cut. We reckon this closely resembles the construction of the finite boost operator that we constructed, where the amount of time evolution grows linearly with the distance from the entanglement cut located at x = 0. In a continuous Lorentz invariant theory, the vacuum entanglement Hamiltonian of half-space is famously known to be given by the boost operator which is given by K = dx xT 00 at the t = 0 slice. Therefore, the tensor network is naturally recovering this results albeit only approximately. Our naive guess of the boost appears to match the actual entanglement Hamiltonian reasonably well, for different values of c, not restricted only to c → 0, even though we solved the spectra of the boost operator only in that limit. This form of the boost generator (4.67) is confirmed yet again when we encounter integrable models, where exact results there are again available. We note however, that the dispersion relation for generic c is non linear in p and it is not immediately clear how to decipher the invariance of the ground state under the boost as a symmetry transformation in the space-time graph. This would be further illustrated in the correlation function, where the c → 0 and c → 1 results admit immediate interpretations, but it is more obscure for generic values of c.
In the discussion below, we will inspect the commutation relations between the entanglement Hamiltonian and the creation/annihilation operators a p , and show that the algebra closely resembles our naive guess. It also closely resembles the results of the corner transfer matrix to be discussed in section 5.2.

Computing the algebra of the entanglement Hamiltonian
The half-space entanglement Hamiltonian is only defined within the positive half line. We can complete the operator into one that acts on all of space by considering H −H whereH is the entanglement Hamiltonian of the complement of the region corresponding to H, in this case therefore, the "negative half region". We consider the commutator [H −H, a p ], i.e. links connected on sites x = −L + 1 2 , · · · , − 1 2 . To make comparison with the boost operator that we constructed explicitly, we will particularly try to obtain the algebra of the entanglement Hamiltonian obtained above for c → 0. Also, in the limit that L is very large, we assume that we can ignore the corrections to the entanglement Hamiltonian that arises from the other entangling boundary in a finite system. 4 In this limit, the entanglement Hamiltonian therefore takes the form for some appropriate normalization N , and ν is some constant.
We thus have The x sin p term is already very suggestive that the action of Γ involves ∂ p a p , which is the infinitesimal version of the boost transformation that we introduced in (4.36), although it does contain other terms. We will show that there is a similar issue in the context of integrable models, and yet it is known how appropriate massage using properties of ∞ lead to the (appearance) of the correct algebra. This is a subtle issue that perhaps could only be clarified on more precise control of boundary conditions at infinity. This issue will be discussed in detail in section 5.2.2.

Correlation functions -visualizing the light cone
In this section, we would like to inspect the two point correlation functions evaluated wrt the "ground state" that we have defined. We would like to inspect the causal structures intrinsic to these tensor network and their manifestation in the correlation function. It is also instructive to see how the boost symmetry manifests themselves in the correlation functions, at least in some limits such as c → 0 and c → 1.
The correlation functions can be readily obtained. We relegate some of the details to the appendix, while keeping only the main results. The complete expressions of the correlation functions are presented in the appendix. We reproduce here the limit where c → 0.  There are various special limits in which the results are particularly illuminating. First, one could take the c → 0 limit, where each tensor behaves like a swap between the neighbouring fermionic sites. It is noteworthy that in (4.70), only half of the modes contribute. It is a direct consequence that in this limit, half of the sites are decoupled from each other. Half of the sites contribute to modes moving right i.e. have positive energy for p > 0; the other half contribute to left moving modes, with p < 0 having positive energies. At the other end, where c → 1, the time dependence completely drops out. The ground state however, is chosen to have non-trivial entanglement, since it is still annihilated by half of the momentum eigenmodes taken as annihilation operators, as their eigen-energies tend to 0 − in the non-dispersive limit.
We note also that when we further take the limit L → ∞ in the c → 0 regime, the correlation functions reduce to the result of a free massless fermion in 1 + 1 d.
In figure 14, 15 and 16 we present the correlation functions corresponding to c = 0.1, c = s = 1/ √ 2 and c = 1 respectively. We use the notations C 00 ≡ 0|a 2x (x, 0)a † 2y (0, t)|0 and C 01 ≡ 0|a 2x (x, 0)a † 2y+1 (0, t)|0 . The lattice is finite, where we take L = 200. We can see the light cone clearly in the figures, demonstrating the causal structure that follows immediately from the construction of the tensor network by local unitary. Moreover, as already alluded to in the previous subsections, a given topology of the tensor network graph gives only the upper bound of the size of the light cone. The precise values that the tensor takes control the actual size of the light cone. At c = 1 the light cone has shrunk to vanishing size as expected of a non-dispersive theory. At c = 0, the dispersion relation is linear, and the invariance under boost translates into an invariance of the correlation functions under a simple transformation in (x, t) coordinates: namely the Lorentz transformation. The hyperbolas marking "equipotential" lines are clearly visible in figure 14.

Comments on Bosons
This section might have conveyed the message that the network was built in from the beginning a discretization of the continuous path-integral and it is always possible, at least with some notion of a long-wavelength limit, to recover the physics of the continuous field theory. There is however a reason why the paper focused on a discussion of Gaussian models of fermions rather than bosons. The authors have begun the journey with the hope of obtaining a tensor network of free bosons. Following almost exactly the same strategy as we have described for the fermions in this section, it is found that we invariably end up with a pair of eigen-modes, one satisfying the usual commutator [a p , a † q ] = δ p,q , while the other pair satisfying the commutator with the wrong sign. This is somehow reminiscent of the issue encountered when quantizing bosons behind the black hole horizon. (See for example [20].) We have not understood the physical reasons for such a behaviour other than the technical reasons why it ended up that way. It may be suggesting that there is potential obstruction to supporting particular kinds of operator algebra in a given network. This is not unheard of, as in the case of chiral fermions which are known to be impossible to be simulated on a discrete (spatial) lattice. Discretising both space and time in the case of a tensor network might suffer additional obstruction which is an interesting subject in its own right. We relegate the details of our (failed) attempt to model free bosons in a tensor network in the appendix.

Generalization to Integrable models
In the previous section, we have worked with a tensor network that corresponds to a system of effectively free fermions i.e. the time evolution is based on a quadratic Hamiltonian. In this section, we would like to generalize our consideration beyond completely free theories. To retain some degree of analytic control, we will focus on a family of integrable models in 1+ 1 dimensions, namely the XYZ model. Much of our discussion however continues to hold for more general (integrable) models. Integrable model is a vast subject. It is impossible to give a complete account of this subject. Our perspective is one based heavily on the set of classical statistical models which will be taken as Euclidean continuation of our quantum tensor network models. One classic reference on these classical lattice model is [21]. Our notation is mainly inherited from the review of [22]. It takes a somewhat more modern perspective compared to [21] and sets up the model directly using Lax operators L n,f (v) as building blocks. The Lax operators satisfy where subscripts f i denote auxiliary spaces on which the respective matrices act, and n denotes the physical Hilbert space at site n. The parameters a, b are the spectral parameters. R denotes the "R matrix" satisfying the Yang-Baxter equation.  Generically, one can take R = L as a specific solution. The solutions are obtained through the (homogenous) matrices T N,f (v), which is defined as 2) where the lattice has N physical sites. Taking the trace wrt the auxiliary space f defines the transfer matrix T (v) which is equivalent to the "row-to-row" transfer matrices in the classical models (in periodic boundary condition) reviewed for example in [21]. Since T (v) evaluated at different v commutes, they define a set of commuting operators, which underlines integrability. The quantum Hamiltonian H is often taken as The energy eigenstates are constructed using T N,f (u). The above relations are generic. In the particular families of XYZ models, each physical site accommodates a spin 1/2. The auxilliary space f can be chosen to be two dimensional, matching the dimensionality of the physical Hilbert space at each site. In this case it was shown that and B(v) can be used to constructing eigenstates. The generic XYZ algebraic Bethe ansatz looks somewhat cumbersome. These can be found for example in [21] . The equations simplify significantly if we focus on the XXZ models. To get a feel of the schematics, we have then where |0 is the reference state corresponding to all spin up in the physical sites [21][22][23]. The important point is that the parameter v i here controls the momenta and energy of the eigenstates. In fact where p(u) is the momentum carried by the "excitation" B(u i ). The parameter λ is a constant depending on the model itself. The meaning of these parameters defining the XYZ families of integrable models are explained in the appendix C.

Analytic continuation of the spectral parameter
Caution has to be taken here about the spectral parameter. Here, we adopt the choice of parametrization in Faddeev, so that T (v) defines a unitary evolution. The parametrization in [21,23,24] is based on the statistical model however. One such example is the 6-vertex model, which is related to the quantum XXZ model is reviewed in the appendix. Generically for square lattices one can define the so called row-to-row transfer matrix, which is denoted as T E (v). It is related to the quantum transfer matrix defined above by When we work with such "Euclidean" versions of any operators we will include the superscript explicitly. This connection between the Lorentzian and Euclidean signatures in a discretized world is noteworthy. The usual procedure asserts that a quantum Lorentzian theory in continuous space-time can be connected to a classical theory in one higher dimension by analytic continuation -the Lorentzian time t is continued to a Euclidean time τ E via However, in a discrete space-time, how such a continuation should be defined to connect the Lorentzian quantum model with the classical model becomes less clear. At least in the context of integrable models, the spectral parameter takes up the role of continuous time, allowing one to continue between signatures. It is worth understanding in greater depth whether a continuation can be defined more generically.

Inhomogenous Algebraic Bethe ansatz vs the tensor network
To make connections with tensor networks, we need to construct a tensor network with a well defined causal structure. The row-to-row transfer matrix makes the connection with local unitaries somewhat obscure. We would like to construct an integrable tensor network made up of local unitaries that takes a similar form as the fermions considered in the previous section. It turns out that this problem has in fact been considered in the literature [22,25], although that was well before the dawn of the notion of tensor networks. Operators known as the inhomogeneous transfer matrix were considered. It takes the following form [22]: where there are 2N sites here, and f denotes an auxiliary space. L i,f (λ) are the Lax operators that satisfy the relation R a 1 ,a 2 (w − u)L n,a 1 (w)L n,a 2 (u) = L n,a 2 (u)L n,a 1 (w)R a 1 ,a 2 (w − u), (5.10) where R a 1 a 2 is a set of R matrix that acts in the auxiliary spaces a 1 and a 2 . This is the star-triangle relation that underlies integrability. Now define Suppose one identifies where L n,f (u) = P n,f l n,f (u), (5.14) for P n,f corresponding to a SWAP operator between the spaces n and f , and V is a shift operator that takes n → n + 1. One can see that (5.13) indeed recovers a network that is exactly analogous to the free fermion model that we considered in the previous section. Compare for example with the evolution of "horizontal" Cauchy surfaces illustrated in figure 4. The eigen-modes of this evolution can be solved again by the Algebraic Bethe Ansatz explained in detail in [22], similar to the homogenous case . In the following however, this is not the path that we are going to follow. The reason being that as already demonstrated in the previous section, to discuss Unruh effect we find it useful to obtain the form of the eigenwavefunctions as well, and the Algebraic Bethe Ansatz is not particularly convenient for the purpose. We would instead like to work with coordinate Bethe Ansatz. These models that we work with, particularly the XYZ family of models, are well known classical lattice statistical models whose partition function has been studied in depth in the literature. A comprehensive review can be found in [21]. The XYZ model can be recovered from the classical 8-vertex lattice model. The definition of the 8 vertex model is reviewed in the appendix. In particular, the definition of the model is summarized in figure 23. We also note that the inhomogenous transfer matrix defined above is essentially generating the so called "diagonal-to-diagonal" transfer matrix, as opposed to the usual "row-to-row" transfer matrix. This corresponds to rotating the classical integrable lattice model by 45 degrees (which can be readily illustrated in the simplest limit of the XXX model). We will explore the change of basis between the row-to-row and diagonal-to-diagonal picture in section 5.2.3. We solved the eigen-wavefunctions in a way directly analogous to the case of the fermions near the ferromagnetic ground state |0 . Since the actual computation is very similar to the previous section, we have relegated the details to appendix C, where eigen-functions and the form of correlation functions for the XXZ model are explicitly computed. It is noteworthy that the excitations above the reference state |0 cannot be made "positive definite" in the way of the fermions discussed in the previous section. It necessrily contains both "positive" and "negative" excitations. (See the discussion near (C.13). ) On the other hand, boost operators can be constructed in an exactly analogous manner as in the case of fermions. All the computations at least for single spinon states have a direct analogue in the integrable model. We note that there are limits of the parameters where the dispersion relations simplify, and the correlation functions again resemble the case in the previous sections recovering the usual invariance along hyperbolas. (See figures 24 and 25 and the discussion nearby.) As we are going to discuss in the next section, the way we have constructed the boost operator is a (Lorentzian) approximation to the boost operator that can be constructed in an integrable model based on the corner-transfer-matrix. The ferromagnetic ground state is indeed invariant under such a boost. That it is a direct product state however, can be traced to the observation we made above -that excitations around the state cannot be made positive definite (in the "first Brillouin zone" for energies). We will discuss the relation of this fact with the Reeh-Schlieder theorem in section 5.4.

Corner transfer matrix
The corner transfer matrix (CTM) was introduced by Baxter [21]. The corner transfer matrix is illustrated in the picture figure 18 It was observed that the corner transfer matrix has a discrete spectrum even in the thermodynamic limit, and that in that limit its spectrum can be exactly solved in various models, such as the XYZ model, making it a very powerful tool. It has been noted that the CTM can be used to compute entanglement entropies of integrable model. The partition function Z of the statistical model can be expressed in terms of the CTM's. It is shown that the reduced density matrix of half-space can be expressed as [21]   The partition function Z is given by tr(ρ half ) = tr(A.B.C.D). Spectra of the reduced density matrix and their corresponding entanglement entropy, of the XYZ model for example, have been discussed in detail [21]. So now we have two loose ends to be tied together at least in the context of integrable model in which exact solutions are more readily available.. On the one hand, we have half-space reduced density matrix ρ half which can be solved exactly. On the other hand, we have constructed naive versions of boost operators in the previous section, which appears, at least in some analytically controllable limits, to illustrate the physics of the Unruh effect. This in turn should be a manifestation of half-space entanglement. Here, we would like to discuss the connections between the boost operator and half-space entanglement in the broader context of integrable models, and demonstrate that such a connection is generic. The boost operator can indeed be defined, and whose algebraic structure has been explored. In fact, our naive construction is a generalization of the boost operator that has been constructed in the past that preserves lattice symmetries, allowing for some sensible approximations when the boost is not an exact symmetry of the space-time lattice. The CTM defines rotation in Euclidean signature by 90 degrees in the row to row basis. Up to a basis change, it is indeed the same constructions of the boost operator as we have, apart from the fact that we allow for more general angles of rotation/boost by choosing the step size of ladders. Our construction does not correspond to exact symmetries of the space-time lattice, but only approximate ones. Moreover, our construction is based on Lorentzian signature in the "diagonalto-diagonal" basis, which also differs from the usual "row-to-row" basis. We will discuss how the analytic continuation to a Lorentzian signature is performed, and how the basis transformation is performed, at least perturbatively in the spectral parameter λ.

CTM and Lorentz transformation
It was observed in the early days that the CTM is related to the notion of a boost operator [24,26]. They have constructed the generator of a boost operator explicitly and obtained the commutation relations that they satisfy. For concreteness, we will explore the families of 8-vertex models, which is related to the XYZ models, and reduces to the XXZ and XXX models in some limits explained in the appendix near equation C.3. By setting the boundary conditions of the Lax operator L E n,f to satisfy the boundary conditions L E n,f (u = 0) = P n , (5. 16) it is readily shown graphically (see figure 19), that nH XY Z (n, n + 1), (5.17) where H XY Z (n, n + 1) = − 1 2 (J x σ x n σ x n+1 + J y σ y n σ y n+1 + J z σ z n σ z n+1 ). (5.18) Using properties of the CTM, where Figure 19. Expansion of the corner transfer matrix around u=0. Picture courtesy [24].
and that where λ is a parameter defining the 8-vertex model (see figure 21) and ogether with the small spectral parameter expansion (5.17), it implies that [24,26] A E (u) = exp(−uK).
HereÂ E is customarily defined as A E (u)/a 0 , where a 0 is the largest eigenvalue of A E . The matrix is known to generate a 90 degree rotation in the Euclidean plane. This tells us that A E are precisely the Euclidean version of the boost operator. To obtain inversion of (x, t) → (−x, −t) which is a 180 degree rotation in the Euclidean plane and which is the ingredient necessary in the Bisognano-Wightman theorem, it is given by (A E ).(B E ) t , using figure 18. It is shown that the Lax operator and H XY Z (n, n + 1) satisfy the following commutation relations (5.22) where L E n,f (v) denotes derivative of the (Euclidean) Lax operator wrt v. This relation follows from the Yang-Baxter equation, and can thus be generalized to generic integrable model, by identifying the Hamiltonian with H(n, n + 1) = L n,n+1 (0)P n,n+1 , where the derivative is taken wrt to the spectral parameter v. Now similar to the case of the free fermions, one can define a "completed" modular Hamiltonian (4.68), and get nH XY Z (n, n + 1). (5.24) i.e. We simply complete the boost operator by extending the sum from n = −∞ → +∞. Using (5.22) gives (5.25) where T N,f (v) are the transfer matrices defined in (5.2). Therefore, the operator K c shifts the spectral parameter in B(u) defined in (5.4). Using (5.6), one sees that the shift of the spectral parameter v under the effect of the generator K c is to vary p.
There is a generalization of the Lorentz algebra [24]. By defining and that C 0 = P i.e. the translation operator shifting globally by 1 lattice site, and C 1 = H. It would reduce to the usual Lorentz algebra when C 2n = P and C 2n+1 = H. This is argued to occur in appropriate continuum limit [24]. We also note that, as already alluded to in section 5.0.1, the rotation operator A E that generates rotation in the Euclidean plane can be continued to a Lorentzian boostÂ by taking u → iu.
There is an interesting lesson here.Â(u) as a tensor network has been fixed, as shown in figures 18 and 19 . Therefore, for a partition function defined at some fixed value of u, A E (u) generates a rotation on the network by 90 degrees. On the other hand, the value of u also plays the role of the angle that is actually rotated. If we consider A E (v) at some other value of v, such as v = u/2, then A E (u/2) is approximating the operator that generates a rotation by 45 degrees on the network defined at u. This can be readily checked in the small u limit discussed in section 5.2.3. In other words, any given effective angle that is actually rotated can either be approximated by the ladder structure that we have introduced, or varying the variable u, allowing it to deviate from the value that defines the tensor network.

A comment on the commutator [K c , T N,f (u)]
In the previous section, we reviewed the construction of the boost generator K c , and that the Yang-Baxter equation (5.1) directly provides a way to construct a boost generator that shifts the spectral parameter appearing in the transfer matrix T E . This can be compared with the effect of our boost operator in equation (4.36), which looks very similar. There, it appears that while the commutator is producing the desired term ∂ p a p , it contains an extra term. We would like to inspect the mechanism of equation (5.25) that generates a simple shift in the spectral parameter. As we are going to see -this is the subtlety that comes with an infinite space. It involves and infinite v dependent normalization, and also a push of the discrepancy towards infinity. We do not know as yet whether this makes physical sense, but we present the detailed mechanism.
In order to see this, we will again restrict our attention to the XXX model where the expressions simplify.
In the XXX limit which is reviewed also in the appendix C, we take λ = π + , keeping /u fixed as u, → 0. We define the new spectral parameter to beũ = u/ .
We note that σ + annihilates the reference state |0 . Also σ z (n)|0 = |0 . Therefore for any given lattice of N sites, the single spinon state is given by where σ − appears at exactly one site n, and n is summed over all sites. In this limit we have and that (ũ + i/2) N /(ũ − i/2) behaves like normalization for any given N . Now naively we find where Formally the corner transfer matrix is defined for N → ∞, where N is the number of sites. We keep N explicit for now. One can see that very similar to the case of free fermions considered in the previous section, a naive computation of the commutation relation does not appear to be consistent with (5.25) -the term that is not linearly dependent on n does not have the right form so that together this is an over derivative in p. On the other hand, the result (5.22) should guarantee that this is true.
On the other hand, the left over terms, namely can be further combined with addition and subtraction of This can be continued indefinitely. Then as we collect all the terms involving σ − (n), we would obtain a better and better approximation of the term in the limit N → ∞. Of course here there are two subtleties. First the normalization N (p) defined in (5.34) plays a crucial role in the algebra, and yet its value is not well-defined in the large N limit. Secondly, we are adding and subtracting terms to push the discrepancy of the result from a total derivative to infinity. At present, we only present the mechanism that led to (5.25) at the spinon level. This mechanism should also be equally applicable to the free fermion case. This is evidence that the our boost operator, which generates finite rescaling in p is indeed consistent with the infinitesimal transformation generated by the entanglement Hamiltonian, despite appearance. Whether this is consistent with sensible and physical boundary conditions are to be investigated and clarified.

CTM in the diagonal-to-diagonal-basis
The CTM is expressed in terms of the "row-to-row" basis in the above. However, our construction of the tensor network so that the time evolution is explicitly made up of local unitaries are more conveniently expressed in terms of the "diagonal-to-diagonal" basis, which is the same lattice rotated by 45 degrees. The story developed above for row-to-row transfer matrices can be translated to the "diagonal-to-diagonal" basis. To that end, we need to obtain a rotation matrix J E (u) that rotates by 45 degrees. This can be constructed. This is basically given by the tensor network sandwiched between the Cauchy surface at 45 degrees to the horizontal Cauchy surface in the diagonal-to-diagonal basis. This is shown in figure  20.
Small u limit also note that J E would be responsible for transformation between the "row-to-row" basis and the "diagonal-to-diagonal" basis. Given any operator O rr in the "row-to-row" basis we have The differential equation ( (5.39) The effective boost operator in the diagonal-to-diagonal basis is thus given by J E −1 (u)KJ(u).
In the limit that u is small, one can readily show that J is given by the identity matrix to leading order. Therefore, We can also inspect the large u limit. This can be computed very readily in the case of the XXX model. Things simplify tremendously because the Lax operator already discussed in (5.30) and also the l n,f operators defined in (5.14) simplifies to L n,f (ũ) =ũI n,f + iP n,f , l n,f = iI n,f + uP n,f . (5.41) As shown in (5.13), each cross in the tensor network is simply given by l n,f up to an overall normalization. The boost operator K in (5.17) reduces to Note that we have dropped a constant term since that only changes the overall normalization of the CTM and do not contribute to any of the commutation relations. The boost operator in the diagonal-to-diagonal basis can be computed by conjugation by J defined in (20) as in (5.38). In the largeũ limit, this can be computed readily. The leading term in the largeũ expansion is contributed by the replacement l n,f → P n,f . This immediately gives This appears like a singular transformation. Despite building the tensor network by local unitaries, it appears that half of the degrees of freedom (the odd sites) are missing. This is because odd sites are thrown towards infinity n → ∞. In fact u → ∞ limit is the infinite boost limit which appears to decouple the left and right moving modes.
One could readily obtain a few sub-leading terms, expanding in 1/ũ and find that they involve growingly non-local terms for higher orders of 1/ũ. In both limits however, the boost generator takes on a simple local form in the diagonal-to-diagonal basis.

The eigen-modes of the corner transfer matrix
In the discussion of free fermions, we present the eigen-wavefunctions of the candidate boost operator that we proposed. Here, we would like to discuss the eigen-modes of the boost operators and compare them with the results following from our guesses. There has been considerable amount of work exploring the entanglement entropy following from the reduced density matrices constructed from the CTM as in (5.15). [27][28][29][30]. It is found for example, that near the critical points, the resultant entanglement entropies recovers a log divergence and whose coefficient matches with the expected central charge of the corresponding CFT [29].
Here, we will particularly focus on the XXZ model, in which the eigen-wavefunctions have been solved explicitly [23,31], in addition to the eigenvalues.
The set of solutions take a very simple form by noting the following. Eigenstates are given by where B(u) is the operator defined in (5.4).
Here, the rapidity parameter u is related to the actual lattice momentum by where λ, u are those that appeared in equation (5.6). It is observed in [24,26], and recollected in (5.25) that the boost operator K generates a shift in u. Therefore the eigenstates transform as This suggests that the correct eigenstates of K is given by Fourier transforming in the rapidity z. For the "single spinon state" for example, We would like to comment on the relation of this solution with the fermionic solution. Now, consider the critical limit which is explained in the appendix B, in which |∆| → 1. In that case, λ → 0, and we can expand (5.45). If we in addition also assume that z is small, we get This means that p = 1/i ln(−1(1 − λ − 2λz)). This says that p has a background value determined by the coupling of the model λ, and a dependence on the rapidity which is a variable that changes the momentum δp(z) = 2iλz. (5.50) One notes that δp depends linearly on z. Therefore the eigenstates (5.47) would take the form |l = dz z z l B(z)|0 .
which is a power in δp(z). This can be compared with (4.42), where we find l ↔ κ. 5 Here we note the following. The reference state of the ferromagnetic phase |0 as introduced in (5.5) is a direct product state in configuration space. However, it is an exact eigenstate of the boost operator K. Since K can be defined up to some constants, it means that this state can be taken as invariant under boost, even though it contains no entanglement.

A comment on the Reeh-Schlieder theorem
This can be contrasted with the case of the fermions, where we demonstrated that the ground state is highly entangled, and at the same time (approximately) invariant since the boost operator (approximately) preserves the notion of positive/negative energy.
In AQFT, there is an important theorem, namely the Reeh-Schlieder theorem, that guarantees that in a translation invariant theory with a positive definite energy operator, the ground state is cyclic and separating with respect to any algebra A(O) associated to any subregion O. In particular, that means that no local operator can annihilate the state [16]. This implies that the ground state is highly entangled, which is the crucial ingredient of the Unruh effect. This is clearly not the case for the reference state |0 in the ferromagnetic phase for the 6-vertex model. i.e. This state |0 is not a cyclic and separating vector in the Hilbert space. In fact, it is a direct product state with no entanglement at all. At the same time when we solve for "excitations" around this reference state, we find that the spectra always involves both one eigenvalue and its complex conjugate. (The details of the states can be found in the appendix.) This means that there is no natural notion of "positive energy" states in this case, which can be contrasted with the fermion ground state that we constructed.
Of course, in the tensor network, the notion of "positive energy" is itself ambiguous. Since the time evolution is discrete in units of ∆t, energy is only defined up to 2π/∆t. Moreover, in any finite (spin) system such as the 6-vertex model considered here, the Hilbert space is finite dimensional and so energy is lower bounded by definition. In such cases, is there anything we can learn from the Reeh-Schlieder theorem as guidance to the entanglement structure of the ground state?
Here, we note that the proof of the Reeh-Schlieder theorem made use of the following. The assumption of a positive energy operator allows one to decompose any (local) operators into the following [16] Q(t, x) = Q + + Q − + Q 0 , Q − |Ω = Q + * |Ω = 0, (5.52) where Q ±,0 corresponds to positive, negative and zero modes under a Fourier transfrom wrt time. i.e. Translation invariance requires that eigen-operators can be further decomposed as sum over operators with definite momenta. Therefore These Q i (p) cannot be local in space as operators with definite momentum. Therefore, (5.53) implies that the ground state is annihilated by roughly half of all the operators built with definite momenta. This guarantees that the vacuum is highly entangled, and conceivably (although we haven't produced a rigorous proof, examples are easily constructed) that the reduced density matrix can be inverted. This boundedness of the spectrum turns into properties of correlation functions, in which it is observed that correlation functions involving any local operator Q's F (x 1 , x 2 · · · ) = ΨQ(x 1 )Q(x 2 ) · · · |Ω , (5.54) can be extended to an analytic function over a tubular region in complex coordinates Here we denote x a = (t a , x i a ), and η i are d-vectors lying in the forward cone. The computation at real x i becomes the boundary values of this analytic function. This analytic extension to an analytic function over an extended region was crucial towards showing that |Ω is cyclic and separating. As we see above, this analyticity of correlation functions is directly related to entanglement of the ground state, although in a discrete system the procedure of such analytic extension would become obscure even though the cause of entanglement is still very much applicable.

A comment on the anti-ferromagnetic case
Before we end, let us comment on the antiferromagnetic case. As it is well known the antiferromagnetic case (corresponding to ∆ < −1) has a continuous limit near the critical point where the model admits a description as a non-linear sigma model that is Lorentz invariant [21,22]. The ground state of the anti-ferromagnetic state is built up from the |0 state by populating the reference state by spinons B(z i ) until exactly half of the spins are flipped. Naively, such a state has no hope of being invariant under the boost operator K, since K shifts all the rapidities z i uniformly when commuting through the B(z i ). The rescue, as proposed in [31], is that such a shift is not consistent with the Bethe ansatz, assuming that we first take periodic boundary conditions before taking the thermodynamic limit. The Bethe ansatz takes the form of a non-trivial integral equation constraining the density of states which would be violated if we shift the rapidities z i uniformly, thus violating the boundary condition we started out with. It is believed that some non-trivial interplay of restoring the Bethe ansatz and the shift of z i should eventually leave the ground state invariant. This has not been shown directly, although [31] took an alternative route in demonstrating that the true ground state is a K eigenstate perturbatively order by order in 1/|∆|. We believe this amazing interplay is closely related to restoring the cyclic and separating property of the ground state in the antiferromagnetic state. A thorough exploration however is beyond the scope of the current paper.

Conclusions and Outlook
In this paper, our goal is to explore the description of a Lorentzian space-time via tensor networks. This is the first step towards building a controllable tensor network description of more exotic backgrounds, such as more generic curved backgrounds and perhaps ones mimicking AdS spaces and AdS black holes, which would shed light on a covariant understanding of the error correcting code/ tensor network description of the AdS/CFT correspondence which has achieved a lot of successes in static spacetimes [2,3].
We approached the problem by first comparing the framework of algebraic quantum field theory with the tensor network, and find that the tensor network can very naturally fit into that picture, allowing one to formulate interesting questions about the tensor networks in very much the same language as is used in the AQFT. Some results of AQFT, such as the Reeh Schlieder theorem, or modular theory could find explicit physical interpretations in the setting of the tensor network.
In the second half of the paper, we consider explicit toy models. In particular, as a first exercise, we would like to illustrate that physics of the Minkowski space can already be (partially) captured in these simple settings, allowing one to explore questions such as the Unruh effect, which is a close relative of Hawking temperature and Hawking radiation in black holes. This would serve as preparation for future explorations of more exotic features of quantum field theories in curved spacetimes.
Finally, we generalize these constructions to integrable models. We find that our naive guesses of the boost operators in the free fermion system is basically an approximation of an operator known to the literature as the corner transfer matrix, which is found to be a close analogue of the Lorentz boost operator in lattice models. This gives extra support to the methods pursued and we hope to generalize our constructions to curved backgrounds, and to higher dimensions in a future publication.

A Free bosons
In this appendix, we will construct the tensor network that consist of free bosons.
As the case of free fermions, we consider a set of bosonic creation and annihilation operators a † n and a n . which satisfy the usual commutation relation [a n , a † n ] = δ n,n . (A.1) The Hamiltonian that generates the time evolution is 2i a † n−1 a n − a † n a n−1 .
Again, over sufficiently small unit of time ∆t, the time evolution operator U (∆t) can be well approximated as follows, Following the same procedure as in the fermion case, the eigen-operators are given by which follow from the eigen-equations U † a p U = Ea p and U † b p U = E * a p . Furthermore, the explicit form for these eigen-equations are where c and s denote cos ∆t/2 and sin ∆t/2 respectively . It is straightforward to see that In the limit c → 0, s → 1, we have When p is summed, the results become (B.10) When we take the limit L → ∞, the above correlation functions become (B.13) In the limit c → 1, s → 0, we have There is no time dependence in the correlation functions in the limit c → 1, s → 0, as expected of a theory with trivial dispersion relation.
Summing p, the results become 0|a 2x (x, 0)a † 2y (y, t)|0 = When we take the limit L → ∞, the above correlation functions become

C More details on integrable models
We will illustrate in detail how to obtain eigen-wavefunctions in the 6-vertex model, and compute some of the correlation functions in detail illustrating the lightcone effect.
We consider an integrable model. The classical statistical model is defined as in figure 21, which gives the assigned weights to each local configuration. The partition sum is the weighted sum of all configurations. The transfer matrices defined for example in (5.4, 5.7) can be read-off from the statistical model as follows. In those cases, they correspond to the "row-to-row" transfer matrices. It is illustrated in Figure 22. "row-to-row" transfer matrices.
The incoming index correspond to configurations of links at the bottom of the layer, and the out-going index correspond to the configurations of links at the top layer. Each element of the matrix is obtained by doing the weighted sum over links connecting vertices in the middle for given fixed boundary conditions at the top and bottom of the layer.
The condition for integrability is that transfer matrices T and T commute, which leads to a parametrisation of the weights given as the following: k and λ are fixed constants associated with the model and u is variable. These are elliptic functions defined as where H and Θ are theta functions. The details of these functions can be found in [21], chapter 15. When k → 0, snh(u) → sin u. In this limit, the 8-vertex (or XYZ) model then reduces to the 6 vertex (XXZ) model. i.e. (C.1) becomes If we further replace λ = π + while taking , u → 0 and /u fixed, we recover the XXX model. In our tensor network construction based on local unitaries, we are interested not in the row-to-row transfer matrix, but the "diagonal-to-diagonal" transfer matrix, obtained by rotating the square lattice by 45 degrees. Each individual vertex can now be viewed as a matrix with incoming indices from the bottom pair of links, and outgoing indices from the top pair of links. Now, requiring that each such local transformation to be a unitary matrix, we have Here we put a = 1, b a real positive number and c a pure imaginary number i|c|.
We consider the lattice as in figure 23. The lattice repeats itself every two layers of evolution. We would like to solve for eigen wavefunctions based on the coordinate Bethe Ansatz. Eigen modes of the diagonal-to-diagonal transfer matrix obtained via the coordinate Bethe Ansatz can be found in [32]. Our parametrization of the lattice is somewhat different from his, and we will solve it from scratch, borrowing heavily his strategy. We will then show that indeed our solution could have been obtained from his via appropriate reparametrization. It is well known that the 6-vertex partition function satisfies a conversation of arrows. Namely for each of the 6 vertices shown in figure  21, it preserves the number of "down" arrow across the vertex. Therefore, one can take the convention that a down arrow represents a particle, and consider "particle excitations" over a reference state with all arrows up. We have the eigen equations Solving these two equations, we have eigenvalues Λ ± = b 2 cos p + c 2 ± i 2b 2 cos 2 p 2 (b 2 − b 2 cos p − 2c 2 ), (C.13) which we note is again an even function in p. For Λ + , we get ψ L (I) = α + e ipI , (C.14) ψ R (I) = β + e ipI , (C. 15) where Here, we note some important difference from the free fermion case which underlies the fact that the reference state here is a direct product state with no entanglement. One can see that in the limit c → 0, the dispersion relation (C.13) also approaches the linear one as in the case of the free fermion. This is not surprising. It is well known that the critical point occurs at |∆| = 1. Comparing with (C.3), we find that it matches precisely with the c → 0 limit where we recover linear dispersion. Relativistc feature would naively be recovered as well, as in the case of the fermions. But this is not so. We note that the reference state is a direct product state and so it could not possibly resemble the highly entangled structure of a relativistic ground state. This is captured by the fact that the excitations around the reference "ground state" contain both positive and negative energy ones, unlike the case of fermions.

C.1 Correlation functions
We 2i sin(π(y − x − t)) e i π L (y−x−t) − 1 . (C.45) With (C.34) to (C.37), we plot the correlation functions of integrable model with different parameters. We use the notations C 00 ≡ 0|σ(0, 0)σ(2x, t)|0 and C 01 ≡ 0|σ(0, 0)σ(2x + 1, t)|0 . In the following figures, we take L = 200. We can see the light cone clearly in the figures. The light cone is very clear in the figures. It becomes small as the parameter b becomes small. One major difference that led to a final answer that does not preserve Lorentz invariance arises from the sum over all p, whereas half of the p modes were canceled out in the free fermion theory.  Here we would like to discuss also solutions of integrable models with boundaries. In the main text, we have taken the "boost" operator as a Hamiltonian, and discussed the corresponding eigenstates. There is something quite interesting, reflecting the fact that there is some sort of a horizon at the boundary. We noted, particularly We can again take the XXZ model as an example. To construct a tensor network here, we can adopt exactly the same strategy as described in (5.9) where we introduce an inhomogeneous transfer matrix. In the presence of a boundary, it would then take the form assuming that there are 2N physical lattice sites. Following similar routes as in the main text and substituting (5.14) into t, putting w = v, t(v, v) becomes a tensor network with boundaries. Explicitly, suppose N = 2 we have t(v, v) = tr f [K f + l 4,f ]l 23 K 1 − l 12 l 34 . (C.50) In the case of a semi-infinite lattice, the right boundary charecterized by K + would be taken off to infinity.
Choose a half-space Hamiltonian which is illustrated in figure C.2. This is a special case, in which we are simply taking K − to be proportional to the identity, which is a well known case satisfying the algebraic constraint described above.
In the case of the XXZ or 6-vertex model eigenmodes can be solved in the same manner as discussed in the previous section. The recursion relations (C.11) continue to apply, except that we have to include new relations that apply only at the boundary.
Suppose the boundary link is located at I = 1. Let the boundary be K = hI. Then the extra boundary recursion relation is given by Λψ L (1) = h(b ψ L (1) + c ψ L (2)), (C.51) where Λ has to be the same eigenvalue as in the bulk of the network determined in (C.11). This can be solved by taking the modified ansatz ψ(I) = m e ipI α + (p) β + (p) + n e −ipI α − (p) β − (p) , (C.52) where α ± , β ± have been determined in (C. 16) and (C. 19), and m, n are constants that can be readily solved using (C.51) and overall normalization of the wavefunction analogous to (C.23). These are simply typical solutions where the plane-waves are reflected at the boundary. Multiple spinon solutions can be solved similarly. For semi-infinite lattice, there is no extra constraints that follow from the boundary.