Classical Simulations of Quantum Field Theory in Curved Spacetime I: Fermionic Hawking-Hartle Vacua from a Staggered Lattice Scheme

We numerically compute renormalized expectation values of quadratic operators in a quantum field theory (QFT) of free Dirac fermions in curved two-dimensional (Lorentzian) spacetime. First, we use a staggered-fermion discretization to generate a sequence of lattice theories yielding the desired QFT in the continuum limit. Numerically-computed lattice correlators are then used to approximate, through extrapolation, those in the continuum. Finally, we use so-called point-splitting regularization and Hadamard renormalization to remove divergences, and thus obtain finite, renormalized expectation values of quadratic operators in the continuum. As illustrative applications, we show how to recover the Unruh effect in flat spacetime and how to compute renormalized expectation values in the Hawking-Hartle vacuum of a Schwarzschild black hole and in the Bunch-Davies vacuum of an expanding universe described by de Sitter spacetime. Although here we address a non-interacting QFT using free fermion techniques, the framework described in this paper lays the groundwork for a series of subsequent studies involving simulation of interacting QFTs in curved spacetime by tensor network techniques.


I. INTRODUCTION
While by now it is generally expected that black holes will radiate and evaporate [1][2][3] in the presence of a quantum field, the details remain murky. The specific interplay between gravity and quantum fields during the early Universe's presumed inflationary [4][5][6][7] phase is similarly unclear. This is not (necessarily) because we lack a quantum theory of gravity: one expects that at least some of these effects can already be described acceptably in some sort of semiclassical regime. For example, one might consider the semiclassical field equations that relate the classical Einstein tensor G µν to the expectation value T ren µν of a quantum field's suitablyrenormalized stress-energy tensor T µν . Supposing such a regime indeed exists, the real problem is that (1) is hard to solve. On the left hand side, we have the fully non-linear Einstein equations, while on the right sits an out-of-equilibrium, perhaps strongly-interacting, quantum field theory in curved spacetime. These are intractable problems, even separately.
But either problem can be solved numerically. Since initial breakthroughs in 2005 [8,9] fully nonlinear simulations of general relativity are now routine. The quantum field theory is a bit trickier, due to the inability of classical hardware to represent highly entangled wavefunctions. In two spacetime dimensions, however, a tensor network ansatz known as the matrix product state [10][11][12][13][14][15][16][17][18][19][20][21] can efficiently simulate low-energy states of local Hamiltonians, as well as sufficiently short real-time evolutions of them. This is the first in a series of studies attacking (1) by combining numerical relativity with tensor network approaches. The first step is to write down the Hamiltonian evolving the quantum field forward in some time coordinate. This Hamiltonian is then mapped to the continuum limit of a sequence of lattice Hamiltonians. On each lattice, interesting states such as the ground state and other low energy states, as well as thermal states and their real-time evolutions, will be simulated with matrix product state techniques. From that sequence we extract T ren µν , the (renormalized) continuum limit of the analogous lattice expectation value. The expectation value T ren µν can then serve as the source term in a numerical relativity simulation.
A necessary first step towards the use of tensor networks such as matrix product states is a careful translation of the continuum QFT problem into one on the lattice. This paper documents that undertaking. Specifically, we describe the mapping between expectation values in a QFT on a fixed curved spacetime [22][23][24] and the aforementioned renormalized sequence of lattice expectation values.
We study Dirac fermions in a curved, two-dimensional Lorentzian spacetime. Dirac fermions in two spacetime dimensions are among the easiest fields to study on the lattice, because they map straightforwardly to an approximation by staggered lattice fermions [25][26][27] which in turn can mapped exactly to a quantum spin chain [28]. The quantum spin chain can then be simulated by matrix product state techniques without further ado, as we will demonstrate in [29]. In the current derivations, however, in order to not further complicate an already rather involved discussion, we consider a non-interacting Dirac field, which can be solved without resorting to matrix products states, using the free-fermion formalism.
In general spacetimes it is not always clear what states ought to be of physical interest. In the literature one often proceeds by making mode expansions of field operators, whose forms are somehow determined by the problem at hand. These mode expansions might not necessarily relate straightforwardly to the Hamiltonian formulation we use. To make comparisons with results based on them, we consider a class of states called "Hawking-Hartle vacua" [30][31][32][33][34][35][36]. They are defined in spacetimes containing static observers, bounded by a socalled "Killing horizon". In the Hawking-Hartle vacuum, that horizon radiates, but is at thermal equilibrium with its environment. It can be viewed as a thermal state of the Hamiltonian defined by static observers outside the horizon, at a certain "Unruh" temperature. We identify it (in the continuum limit) with thermal states of the analogous lattice Hamiltonian, at the relevant metric's Unruh temperature.
In a quantum field theory, the expectation value T µν of the stress-energy tensor T µν in the Hawking-Hartle vacuum is formally infinite even for free fields, and to extract a meaningful finite number it needs to be renormalized. This same problem occurs for any operator with quadratic terms in the fields. One possible solution, formulated directly in coordinate space and thus appropriate to curved spacetime, is Hadamard renormalization [37][38][39][40][41][42][43][44]. Here, the quadratic operator is replaced with the coincidence (zero-separation) limit of some correlation function, with any divergent terms subtracted before the limit is taken. In fact, the lattice two point functions have the same leading order form as their continuum limits in a short-distance expansion. Thus, substituting them for the continuum correlation functions within the coincidence limit yields the same result.
For illustrative purposes we present four applications in the continuum using lattice data. By order of complexity, they are: (i) computation of ground state expectation values as measured by inertial Minkowski observers; (ii) computation of thermal expectation values as measured by accelerated Minkowski observers (and thus reproduction of the Unruh effect [2]); (iii) computation of thermal expectation values as measured by static observers in the Hawking-Hartle vaccum of Schwarzschild black hole spacetime; and, finally, (iv) computation of thermal expectation values as measured by static observers in the Hawking-Hartle vacuum (often called the Bunch-Davies vacuum) of an expanding universe described by de Sitter spacetime [45][46][47][48][49][50].
This paper is divided into four main sections. The first is this introduction. The second develops the continuum problem to be numerically solved: it defines the free Dirac field in a curved two-dimensional spacetime, explaining its canonical quantization, and how point-splitting regularization and Hadamard renormalization allow quadratic expectation values to be defined. It also introduces examples of spacetimes with a bifurcate Killing horizon and defines the Hawking-Hartle vacuum. The third section outlines our lattice approximation: it introduces the lattice and then reviews the staggered fermion technique to map between its Hamiltonian and that of the continuous field theory, as well as the free fermion formalism to find correlation functions of non-interacting fermions on the lattice. Finally, the fourth section uses the above method to find Hadamardrenormalized expectation values in the four example cases (i)-(iv) previously listed.

II. CONTINUUM PROBLEM
To obtain T ren µν in (1), we must understand how to compute renormalized expectation values of quadratic operators. The goal of this paper is to explain how to do so for so-called "Hawking-Hartle vacuum", describing a horizon in thermal equilibrium, of a theory of free Dirac fermions in two spacetime dimensions. In this section, we define our terms and outline this problem in detail.

A. Free Dirac Fields in Curved Spacetime
We work in a two-dimensional Lorentzian manifold (M, g µν (t, x)). Every two-dimensional spacetime is conformally flat [51], and thus can be covered by finite-size patches in which the metric g µν (t, x) differs from the Minkowski metric η µν ≡ diag −1, 1 by an overall function Ω 2 (t, x), called the conformal (Weyl, scaling) factor. We will always adopt coordinates in which this is manifestly true, The strategy to be outlined will not qualitatively depend upon this choice, though the specific equations we present will.
The manifold will shelter a "Dirac" field of two-component spinors ψ A (t, x), whose spinor components are indexed with capital-Latin letters. In order to simplify the notation, we will usually suppress these. When spinor indices are suppressed (ψ A (t, x) → ψ(t, x) and ψ †A (t, x) → ψ † (t, x)), ψ(t, x) and ψ † (t, x) are to be formally treated, apart from their transformation properties, respectively as two-component column and row vectors.
The spinors form a separate spinor representation of the Lorentz group at each point in M. In curved spacetime [24,52], along with the usual "flat spacetime" gamma matrices γ µ we need "curved spacetime" equivalents γ µ (t, x), respectively defined by their anticommutation relations, When spinor indices are suppressed (γ µ B A → γ µ ,γ µ B A →γ µ ), γ µ andγ µ interact with the Dirac spinors as 2 × 2 matrices. In the manifestly conformal chart (2), we will choosẽ so that (3a) imply (3b). It is convenient to define the Dirac adjoint, along with the inaccurately-named "fifth" gamma matrix, We will distinguish between two derivative operators. The first is the partial derivative, denoted interchangeably by either the usual ∂ µ or a comma, The partial derivative of a tensor generically does not transform as a tensor. To remedy this, one introduces the covariant derivative, denoted interchangeably by either ∇ µ or a semicolon, The covariant derivative satisfies the axioms of a derivative operator [53], notably including linearity and the product rule. It differs from the partial derivative by the subtraction of "connection" terms, such that the full expression transforms covariantly. The specific connection depends on the transformation properties of the operand. Details are available in e.g. [52,53], but in the case of the spinor fields we havē where the expressions in the right column (but not the left) are obviously specific to the manifestly conformal frame (2). We can work outγ µ ;ν (t, x) = 0 [52], which in fact holds generally. We will sometimes use the "arrow" notation Applying (9) we seeψ which is specific to 2 dimensions.
Dynamics will be set by the free Dirac action S and its corresponding Lagrangian density L(t, x), where g(t, x) is the determinant of the spacetime metric. In the manifestly conformal frame (2) we have In flat spacetime, one more commonly uses a Lagrangian density containing γ µ ∂ µ , as opposed to the 1 2 ↔ ∂ µ that appears in ours. The two Lagrangian densities differ only by a total derivative, but the present choice (13) is more convenient for our purposes due to (11), which lets us write The overall minus sign in (13) and (15) would flip if the opposite metric signature from ours were adopted.
The canonically conjugate momentum to ψ(t, x), π(t, x), is defined by Time evolution (Lie transport) of the fields between successive constant-t hypersurfaces of the spacetime is generated by the Hamiltonian, The stress-energy tensor T µν (t, x) of the Dirac field is [52] T µν (t, A field's stress-energy tensor measures its response to an infinitesimal diffeomorphism, in the sense that where dΣ µ (t, x) is the integration measure on a hypersurface Σ, generates Lie transport of the portion of that field from Σ along ξ ν (t, x). The Hamiltonian (17) can also be derived from (19) by fixing Σ to be constant-t hypersurfaces with spatial coordinate x and choosing ξ ν (t, x) to point along lines of constant x (see Appendix B). Unlike (17), however, (19) could also be used to propagate the fields along a different foliation than the constant-t one, an issue we will not explore further here.

B. The Quantum Field Theory
The theory can now be quantized e.g. canonically, by introducing the canonical anticommutation relations [24] {ψ which require ψ(t, x) to be understood as a field of operators. Inserting (16b) into (20), we find Note that the Dirac delta in (21a) makes the powers on the Ω(t, x) and Ω(t, x ) arbitrary, provided they sum to −1. We have chosen the present form to ease discretization later.
We take as given a Hilbert space of quantum states |Ψ(t) , such that the so-called Hadamard function (which is actually a 2×2 matrix containing 4 functions, see Appendix C) is finite, provided the primed and unprimed points not be lightlike separated: notably, provided they are not the same point. The primed index in (22) transforms as a spinor with respect to Lorentz transformations at (t , x ), rather than at (t, x) as does the unprimed index, but these indices will again usually be suppressed, G x, t , x ). The superscript (1) is a decorator often used in the literature to distinguish G (1) (t, x, t , x ) from other two-point functions such as the Feynman propagator (which will not appear here).
We will most often be interested in the equal-time Hadamard function As will be discussed below, the Hadamard function can be used to regulate and renormalize "quadratic" operators involving quadratic terms in the fields. Many interesting operators, including the stress-energy tensor (18), fall in this class. However, the stress-energy tensor must be further manicured even after renormalization in order for its expectation value to meet, for example, the divergence-free condition T µν (t, x) ;µ = 0. This poses a significant, though surmountable [44,54], complication that will not be addressed here. In this study we instead use as prototypes the condensate C I (t, x) , the pseudo-scalar C 5 (t, x) , and the current j µ (t, x) , Together, they span the four possible quadratic expressions of the form ψ †A ψ B that can be obtained with components of the Dirac spinors (see Appendix C).

C. Point-splitting regularization
The expectation values on the right column of (24) are formally infinite and thus ill-defined. We will replace them with well-defined quantities via so-called "point-splitting" [37][38][39][40][41][42][43][44]. Thus, we regulate our quadratic operators by explicitly treating them as (t , x ) → (t, x) "coincidence" limits of various traces of the Hadamard function, see Figure 1. For example, in flat spacetime, we would write where the trace is over the suppressed spinor indices, and the limit is along the unique geodesic connecting (t , x ) to (t, x).
In a curved spacetime, the point-split expressions will be more complicated, because as e.g.ψ(t , x ) is paralleldragged to (t, x), it will be deformed and rotated by the spacetime curvature. This is expressed symbolically by introducing the "spin parallel propagator" J A → 1 is the identity over bispinor indices. As always, the limit is understood to be along the coincident geodesic. These are, in fact, the parallel transport equations, so that multiplication ofψ(t , x ) by J (t, x, t , x ) correctly parallel-drags the suppressed spinor indices from (t , x ) to (t, x).
To point split, we now treat the formal expectation values in (24) as the replacements Notice that, had the spin parallel propagator J A B (t, x, t , x ) not been included in these expressions, the respectively unprimed and primed spinor indices of G x, t , x ) would transform at different points, and thus could not be traced over to yield a scalar.
To get expressions we can compare to lattice data, we need to specialize each of (29) to the manifestly conformal metric (2) with t = t . We will achieve this by expanding J (t, x, t, x ) in the spatial coordinate separation r ≡ x − x. As we will see in (33), the strongest divergences in (29) are proportional to 1/r. Thus recovery of the coincidence limit requires this expansion to linear order in r. In Appendix A we find

D. Hadamard renormalization
The coincidence limits of these expressions (31) are still undefined, because via G (1) (t, x, x ) they contain terms that diverge as x → x. However, given certain smoothness assumptions, bundled together as the requirement that physical quantum states be "Hadamard", said divergences are completely determined by the local geometry and the mass. In particular, they are independent of the specific choice of "Hadamard" quantum state.
The locally determined contributions to (31), including all the divergences, made by any Hadamard state can then be computed as the first step of a procedure known as Hadamard renormalization. The next step is to subtract the locally-determined terms from each of (31) before the coincidence limits are taken, thus yielding finite answers.
We label the locally-determined terms with the superscript "loc", e.g. j loc x, x ). Subtracting each local contribution from its corresponding bare two-point function within the limit then yields a convergent remainder -which we label e.g. j ren µ (t, x, x ) -with a finite coincidence limit. Next, we identify each of (24) with the relevant coincidence limit: The locally-determined terms were computed in [44]. They are C loc j loc where µ is an undetermined dimensionful parameter that must be "measured". We will choose  where γ E = 0.5772 . . . is the Euler-Mascheroni constant. This arbitrary choice cancels a constant term that otherwise appears in a short-distance expansion of C I (t, x, x ) (see Appendix D), and is thus implicitly made by the standard practice of normal ordering. If a different choice µ were made, our results for the condensate would suffer the replacement Our results are otherwise fixed by the Hadamard renormalization procedure.

E. Hawking-Hartle Vacua
To specify a calculation we must identify a quantum state. In Minkowski spacetime, one "by default" considers the Minkowski vacuum, selected for example by the demand that G (1) (t, x, t , x ) be Poincaré invariant. But there is no general such prescription yielding an interesting state in an arbitrary spacetime. States must be identified in a manner suggested by the problem at hand. Our target in this study will be the Hawking-Hartle vacuum, which describes a (Killing) horizon in thermal equilibrium with its environment.
In spacetimes containing such Killing horizons, the Hawking-Hartle vacuum is, in a free theory, the unique state that is both stationary and possessed of a Hadamard function with no divergences in x independently of x . That such a state appears thermal to certain observers is the essential reason for the Unruh effect, the inflationary power spectrum, and black hole radiance, arguably the "poster children" of quantum field theory in curved spacetime. Given access to a direct simulation of the Hawking-Hartle vacuum, these effects can correspondingly all be studied, perturbed, and prodded in various ways. One could, for example, excite the state, and then study e.g. real time evolution of the von Neumann entropy [29].
In the specific spacetimes we study, the Hawking-Hartle vacuum is often referred to by a specific name. In Minkowski spacetime it is the Minkowski vacuum. In the prototypical case of Schwarzschild spacetime, it is the Hawking-Hartle vacuum, or sometimes the Kruskal vacuum. In de Sitter, it is the Bunch-Davies or Euclidean vacuum. These and similar comparisons are tabulated in Table I In all three spacetimes, we will first identify global observers at constant X in the global coordinates (T, X). Next we will identify static observers at constant x in the static coordinates (t, x). We will use these same labels for the global and static coordinates in all three spacetimes. The corresponding Weyl factors (2) will be written Ω i global (T, X) and Ω i static (x), where the superscript i labels the spacetime. That is, we will write Together, they divide the spacetime into four wedges: the right and left static wedges, and the top and bottom interior wedges. Global observers, whose time T runs upward in the diagram, can cross wedges. Static observers, whose time t runs hyperbolically in the diagram, are confined to one or the other static wedge. Analogous hyperbolic trajectories exist in all four wedges, but are spacelike in the top and bottom ones. In Schwarzschild spacetime, the interior regions would be capped by singular hyperbolae. In two-dimensional de Sitter, the T = const. lines here cover only half of a ring, and so covering the full spacetime requires a second copy of the diagram. The point where the horizons cross is called the bifurcation point.
causally disconnect the static observers from half of the manifold. The Killing horizon is "bifurcate" because of the reflection-like symmetry relating the top and right patches to the bottom and left ones. They meet at the origin of the diagram, called the bifurcation point, where a constant t hypersurface in either static wedge intersects the horizon. The important geometric features are thus the global observers, the static observers, their confinement by the horizon to the static patch, and the bifurcation point. Now, we will highlight these in concrete examples, starting with Minkowski spacetime.
Minkowski.-The global chart (T, X) on Minkowski spacetime is given by called "inertial" or "Minkowski" coordinates. Figure 2 shows an abstract bifurcate Killing horizon spacetime in global coordinates. Viewed as depicting Minkowski space, that diagram shows the inertial observers tracing straight vertical lines. The Minkowski metric is static because it does not depend on T . Vectors parallel to T translations form a symmetry of the spacetime, and are thus timelike Killing vectors.
Minkowski spacetime has another set of timelike Killing vectors, parallel to boosts. Observers parallel to these are called Rindler observers. They are at rest in, for example, the Lass chart which is also static because it does not depend on t.
The Lass coordinate x ranges from −∞ to ∞, but only covers the right static wedge in the diagram 2, where the Killing vector is future-pointing and timelike. In the left static wedge, it is past-pointing and timelike. In the top or bottom wedge, it is spacelike. At the boundary, the Killing horizon, it is lightlike.
The Killing horizon causally disconnects a particular group of timelike Killing observers from half of Minkowski spacetime. Which group depends on which point of Minkowski spacetime is identified with the origin, and upon the choice of acceleration parameter α. This parameter defines the rate at which accelerations diverge as observers closer to the horizons are considered, and can be viewed as characterizing the "strength" of the horizon. This can be quantified by defining the surface gravity κ, where k a is the Killing vector, and the expression is to evaluated on the horizon. In Minkowski spacetime we have κ Mink. = α.
Schwarzschild.-Now consider the Schwarzschild black hole, described by the familiar metric We use the term "Schwarzschild" black hole to describe the 2D manifold partially covered by (40), not the 4D solution to the Einstein equations. Since this line element is independent of t, the constant-ξ "Schwarzschild" observers are timelike Killing observers. The spatial coordinate ξ ranges from 2M (the event and Killing horizon) to ∞.
An analogous chart to the Lass metric (37) is found by defining the tortoise coordinate x to find which is again static. These coordinates cover the same region as the Schwarzschild coordinates, but have infinite range. The horizon is at x → −∞.
The region of Schwarzschild spacetime outside the ξ = 2M Killing horizon is thus analogous to the static patch of Minkowski spacetime. This can be made obvious by adopting the "Kruskal" chart, which covers all of Schwarzschild spacetime just as the inertial chart covered all of Minkowski. The Kruskal chart is defined by If Figure 2 depicted Schwarzschild spacetime, a singular hyperbola would need to be drawn at ξ = 0. Otherwise, Figure 2 might equally well depict Schwazrschild spacetime in Kruskal coordinates or Minkowski spacetime in inertial coordinates.
The const-ξ "Schwarzschild" observers in the patches are Killing, but they are not geodesic. Their local acceleration diverges towards the horizon at a rate set by the Schwarzschild mass M , which is in this way analogous to α in (37). However, while α parameterized a family of charts upon the same manifold, M parameterizes a family of physically inequivalent manifolds. The surface gravity is De Sitter.-Our final example will be the exponentially expanding de Sitter universe, often labelled dS 2 in the two-dimensional case. This spacetime is most straightforwardly discussed in terms of its embedding in higher dimensional flat spacetime, as in Figure 3. The global coordinates (T, X) with line element cover the full hyperboloid, and may be viewed as projections of the higher-dimensional Minkowski coordinates upon it. Constant-T slices are horizontal rings, which contract until T = 0 and then expand. The "de Sitter radius" α dS controls the rate of expansion.
Constant-X observers follow timelike geodesics, which appear hyperbolic and thus approach lightlike lines in the embedding diagram, their cosmological horizons. These horizons segment dS 2 into causally disconnected The corresponding timelike Killing observers, defining static coordinates (t, x), are not zero energy geodesics, but asymptote to them, as do all worldlines in dS 2 . Thus, each geodesic's cosmological horizon is also a Killing horizon, confining some set of static observers to some static patch, like the one shaded in the diagram. These static patches also wrap around behind the hyperboloid, and thus form spatially compact diamonds, enclosed on either side. The diagram in Figure 2 would need to be "doubled" and then periodically identified to reflect this.
patches. In fact, they are Killing horizons, and the geodesic observers are Killing observers, though this is not obvious from (45).
We can construct static coordinates bounded by a given Killing horizon as which cover half the circumference of the central ring. These observers are accelerated relative to the higherdimensional Minkowski coordinates. The surface gravity is

Hawking-Hartle Vacua
Now let us try to construct a physically reasonable "ground" state in a spacetime with a bifurcate Killing horizon. We will make two demands. First, we would like the equal-time Hadamard function to be independent of the static time t, so that we can interpret the state as "equilibrated". Second, we would like the state not to diverge in any special way at the horizon, since this would be very surprising to the global observers who can cross it.
These stipulations pick out, uniquely [36] in a free theory, the Hawking-Hartle vacuum. To reiterate, the Hawking-Hartle vacuum is the unique state on a bifurcate Killing horizon which is both smooth everywhere and static to the static observers.
In Figure 4  for example from the right half of Σ 0 to Σ s . Since H s is independent of t this translation can be used to define a stationary quantum state, but since the translation is only timelike within the static wedges, that state will generically be singular on the horizon. A nonsingular state that is still stationary with respect to H s can nevertheless be defined as a Euclidean path integral on the Euclidean manifold found by analytically continuing from e.g. one half of Σ 0 to the other. This Euclidean manifold is depicted as the half-cone in the inset. The edge of that half-cone overlaps with Σ 0 in the Lorentzian spacetime, and the path integral is from one edge to the other along the Euclidean time coordinate τ . To fix the singularity at the tip of the cone, τ must be made periodic with period β = 2π κ , resulting in a thermal state of H s at temperature (48).
forward in the static time t and the global time T . The first, H s , evolves the field forward along the worldlines of the static observers -e.g. from the right half of Σ 0 to Σ s -and is therefore independent of t.
The "ground state" of H s , called the static vacuum, is thus also independent of t, our first criterion. However, H s generates a lightlike rather than timelike translation at the horizon. Because of this, its Hadamard function diverges at the horizon independently of the distance between its arguments. It thus violates our second criterion.
Though unphysical, the static vacuum comes up often enough to have earned its own special names in each of Minkowski ("Rindler vacuum"), Schwarzschild ("Boulware vacuum") and de Sitter ("static vacuum").
The global Hamiltonian H g (T ) evolves the field forward in the global time T , e.g. from Σ g to Σ 0 , which is well behaved at the horizon. The "ground state" of H g (T ) at a particular time is thus well-behaved there. However, except in Minkowski spacetime, H g (T ) will depend on T , and through it, on t. This violates our first criterion. Though such a state is potentially physical, it is not "equilibrated" in any useful sense.
Let us now construct the Hawking-Hartle vacuum, which fulfills both criteria. Consider a hypersurface such as Σ 0 in Figure 4, formed by the union of constant-time hypersurfaces in the right and left static wedges. A stationary state of H s prepared everywhere on that hypersurface would be smooth. The difficulties in preparing one stem entirely from the bifurcation point O.
We will resolve them by analytically continuing past that point. That is, we will make the substitution t → τ = −it in the static chart. Because that chart is by definition time-independent, this only changes the signature of the metric from Lorentzian to Euclidean. For example, the hyperbolic trajectories of the static observers in the Lorentzian manifold become circles in the Euclidean one, as depicted in the inset of Figure 4. The "time" evolution operator e −τ Hs advances along those circles.
We can follow those circles to move from the part of Σ 0 in the right static patch to that in the left, without needing to cross the horizon. Of course the Euclidean circles are not physical trajectories, but path integrals along them can be used to define quantum states. These will be stationary with respect to the static Hamiltonian H s , the generator of motion along the circles.
One part of the horizon does lie upon the Euclidean cone: the bifurcation point O. This can be shown to be a literal conical singularity by going to polar coordinates in its vicinity. So long as that singularity persists, the quantum state defined by the path integral on the cone will be singular at the horizon.
Like any conical singularity, however, the singularity can be smoothed by choosing a certain periodicity β U of the relevant angular coordinate, τ , the inverse Unruh temperature [2]. States ρ defined as path integrals on closed Euclidean manifolds are thermal with respect to the Hamiltonian generating motion in the direction of the path integral, Thus the Hawking-Hartle vacuum describes the thermal equilibrium state of the horizon. Despite the apparent complexity of the above discussion, preparing it from a Hamiltonian formulation is simple: just find a thermal state of H s at the relevant Unruh temperature.

F. Summary
Everything needed for a well-defined calculation is now at hand. The steps to be followed are: 1. Choose a spacetime with a bifurcate Killing horizon. (17)

Compute renormalized quadratic expectation values via (32).
We will next show how to approximately follow this procedure using data obtained by numerically simulating a sequence of lattice theories.

III. STAGGERED FERMION DISCRETIZATION
One possible strategy to perform numerical simulations is to map the continuum quantum field theory onto a sequence of lattice models. In this section we outline a means to do this, based on a slight adaptation of the well-known "staggered fermion" discretization [25][26][27].

A. Staggered Fermions
We will numerically simulate the quantum field theory by extrapolating data from a sequence of lattice theories. Thus, consider a one-dimensional lattice whose sites, labelled by n, lie at points x n = na separated by constant lattice spacing a ≡ x n+1 − x n . The lattice may be finite or infinite; when finite, it will have an even number of sites N . We use a staggered-fermion discretization [25][26][27] to relate fermionic operators φ n on this lattice to components of the two-component continuum spinors ψ A (t 0 , x) at some reference time t 0 {φ † n , φ n } = δ n,n , (50a) As we move from site to site, these operators are identified with alternating components of the continuum spinors, The normalization recovers the canonical anticommutation relations (20) and (21) from the lattice anticommutation relations (50) in the continuum limit a → 0. In addition, it rescales φ n to be dimensionless: one of a or Ω(t 0 , x n ) has units of length, depending on whether in ds 2 = g µν (t, x)dx µ dx ν the coordinates dx µ or the metric components do.
Since this prescription involves individual components of the continuum spinors, we need to commit to a representation of the gamma matrices in order to use it. Our choice will be Quadratic operators are transcribed onto the lattice as in [26]. Those which couple differing components of ψ(t, x) are mapped to nearest-neighbour couplings such as φ † n+1 φ n . Those that couple like components are instead mapped to on-site couplings such as φ † n φ n , with oscillating signs introduced while moving from site to site as required by (52). Noting our slightly different conventions from [26], we have Per [26] and our differing normalization of φ n , the derivative term in the Hamiltonian (17) maps to Comparing (53) and (54) to the continuum Hamiltonian (17) yields where dx → a was adopted.
Thus regulated, such quadratic operators are well-defined on the lattice, but their continuum limits a → 0 will diverge. We extract defined quadratic expectation values by way of the point-split definitions made in Section II B. First, we map correlation functions onto the lattice, by first expanding their spinor dependence into components using the representation (52), and then applying (51), as illustrated in Appendix C. Starting from φ † n φ n , n, i same parity; n , j same parity 0 otherwise (56) we make the correspondences x, x , a), Tr x, x , a), x, x , a), The correlation functions in (57d) have been chosen so that the total number of lattice units between sites is the same in both terms. This choice is made so that the coincidence (x → x) limit is correctly renormalized by the locally determined terms (33). It is not necessary to do this in the case of (57b), because the coincidence limit is well-defined in that case. Notice also that the coordinate dependencies of the Ω(t, x) prefactors are located at the same point of the lattice correlator they are proportional to. This is how, for example, the derivative term in (33d) gets reproduced on the lattice.
As alluded to, the expressions (57) are finite in the continuum (a → 0) limit, but the continuum quantities they approximate diverge in the subsequent coincidence (x → x) limit. We can renormalize the latter divergence by subtracting the locally determined terms (33) after the continuum limit is taken. In fact, we find empirically that the continuum and coincidence limits can be renormalized simultaneously, by defining K ren 5 (t, x, ∆n, a) ≡ K 5 (t, x, x + a∆n, a), K ren 0 (t, x, ∆n, a) ≡ K 0 (t, x, x + a∆n, a) − The limit a → 0 with ∆n held fixed now also brings x and x together. The remaining terms may depend on ∆n, but do not diverge in the relevant limit ∆n → ∞. The limiting process is illustrated in Figure 5. In practice the ∆n dependence seems to be quite weak.
We can now make lattice correspondences to the Hadamard renormalized quadratic operators by combining (58) in the a → 0 and ∆n → ∞ limits with (32) and (31), to find FIG. 5: Extracting the coincidence limit from lattice correlators. The lattice correlator between sites at x and x + ∆na, with ∆n the number of lattice sites separating the two operators and a the lattice spacing, is identified with the continuum correlator ψ † i (t, x + ∆na)ψ 0 (t, x) , i = ∆n mod 2. Holding ∆n constant and comparing correlators on successively refined lattices, we get functions of r = ∆na with the same divergent terms in r as their continuum analogues. Subtracting those divergent terms (in practice by fitting them, up to a constant, to the numerical data) leaves a remainder whose difference from the Hadamard renormalized quadratic expectation value ψ † i (t, x)ψ 0 (t, x) ren falls off with ∆n.
These expressions are now both well-defined and straightforward to compute numerically. Indeed, we could now make a Jordan-Wigner transformation [28] to exactly map the lattice Hamiltonian (55) and lattice Hadamard functions (57) to expressions in terms of Pauli matrices while retaining their local character. This enables simulations of the ground state of (55) and other related states of interest using matrix product state techniques, as we explore in Ref. [29] in the context of simulating interacting QFTs.

B. The Free Fermion Method
At present, however, since here we work with a quadratic fermion Hamiltonian H φ (t), we can exploit the free fermion formalism to completely and efficiently characterize its (spontaneous) ground state, as well as thermal states and more general Gaussian states, in terms of the equal-time two-point correlation matrix where the possible time dependence comes from the Gaussian state ρ(t) under consideration.
For instance, let us briefly recall how to characterize the instantaneous ground state and thermal states of a more generic quadratic Hamiltonian H(t) ≡ n,k A nn (t)φ † n φ n (with A nn (t) = A n n (t) * ). We proceed by first identifying a unitary transformation is a diagonal matrix with single-particle energies λ p (t) in its diagonal, that is D pp (t) = δ pp λ p (t). The unitary O(t) in turn defines a unitary U (t) acting on the Hilbert space of the lattice such that it diagonalizes H(t), where c p (t) are a new set of fermionic operators. Then the two point correlator of the Gaussian thermal state with inverse temperature β, is given by which implies In our applications below, we will consider the case of a time independent Hamiltonian H φ . Correspondingly, the above expressions simplify to We can now compare (66) with e.g. (57) to estimate traces of G (1) (x, x ). In the simplest case of Ω(t, x) = 1, zero mass m = 0, zero temperature (β = ∞) and infinite lattice length N → ∞, we can easily compute the two point lattice correlator exactly to obtain where r ≡ a∆n. Comparing with (57) and then (31), we see that the Hadamard divergences (33) are indeed obtained also on the lattice by sending a to 0 with ∆n held constant. More generally, we find empirically that this is also true at finite mass m > 0, finite temperature (0 < β < ∞), and finite lattice length, as well as for space-dependent Ω(x) (given that the normalizations (51) are taken into account).

IV. APPLICATIONS
In this section we apply the formalism described above to compute, in the Hawking-Hartle vacua of the two-dimensional Minkowski, Schwarzschild, and de Sitter spacetimes, expectation values of the Hadamardrenormalized quadratic operators defined in (32), namely those of the condensate C ren I (x) , the current j ren i (x), and the pseudoscalar C ren 5 (x) . Note that since all of the metrics we work in will be manifestly static, our notation here differs from that preceding by the suppression of temporal dependencies, e.g. C ren I (t, x) → C ren I (x) . We will compute the Hawking-Hartle vacuum in four different coordinate charts. First, we will find it in Minkowski spacetime, both as the β → ∞ state set by the inertial observers (Ω Mink. global = 1) and as the β U = 2π α state of the Rindler observers, Ω Mink. static (x) = e αx . We will find equivalent results in both cases despite the quite different lattice Hamiltonians, and in particular will find smooth behaviour when approaching the horizon in the Rindler chart, but only when β = β U .
Next we will find it in Schwarzschild spacetime, (42). We will find that at the "Hawking" temperature β H = 8πM results approach those of Minkowski far from the black hole, and yield smooth behaviour approaching it. Finally, we will find the "Bunch-Davies" vacuum in de Sitter spacetime, Ω dS static (x) = sech(x/α dS ), and again will show smooth behaviour approaching the horizon at the Unruh temperature β U = 2πα dS .  (31) in the ground state within the inertial Minkowski chart, computed on a finite and discrete lattice via the discretizations (58), plotted against the inertial coordinate X. We use m = 0.6. Darkening colours indicate simultaneously decreasing point-separation r and lattice spacing a. Red lines marked with crosses show the results of the infinite size (L → ∞) and coincidence r → 0 extrapolations, themselves detailed in the insets to the top left and bottom right plot. J ren 0 (X, X + r) (bottom left) and C ren 5 (X, X + r) (top right) vanish essentially to machine precision even at finite separation and lattice size. J ren 1 (X, X + r) (bottom right) is nonzero on the lattice at finite separation, but vanishes in the continuum/coincidence limit. C ren I (X, X + r) (top left) approaches a finite value at nonzero mass. The analytically calculated valued (see Appendix D) is plotted as a dotted black lines.

Inertial Frame
We first discuss the somewhat trivial case of Minkowski spacetime in Minkowski coordinates, Ω Mink. global = 1. The results in this case ( Figure 6) provide a baseline to compare our other results against. Here, we have computed each of (31), by identifying each with its lattice counterpart in (58). The superscript ren indicates that the locally-determined terms (33)  left panel. We have chosen m = 0.6 here. The dependence on X is artificial due to the periodic boundary conditions, but we maintain it nevertheless for easier comparison with subsequent results.
Since the locally determined Hadamard terms have been subtracted, the r → 0 extrapolation depicted in the bottom inset of the top left plot in Figure 6 yields an estimate of the associated quadratic expectation value in the field theory, (32). We have not depicted the ∆n → ∞ limit since it has no visible effect. The results here and throughout are at ∆n = 2.
C 5 and J 0 already vanish even on the lattice. J 1 , which measures the current density of fermions in the continuum limit, also appears to approach zero, though slowly; the actual extrapolation in this case yields about 10 −4 , but can be made smaller by extrapolating from smaller values of r.
The only nonzero expectation value in the continuum is the condensate C I , which takes a finite, massdependent continuum value. This quantity measures the degree to which so-called "chiral" symmetry [55] is broken, as it is by the mass term in (12).

Unruh Effect
Thermal states defined by the accelerated Lass chart Ω Mink.
static (x) = e αx have the same two-point functions as zero-temperature states defined by the inertial Minkowski observers: this is the famous "Unruh effect". We reproduce this on the lattice by making the same calculations and extrapolations using that metric as we did I (x) again shows results both at and slightly away from the Hawking increasing temperature β H , with brightening colours indicating increasing temperature (and thus decreasing β). Approaching the black hole, the results away from β H diverge, but those at β H converge to a finite horizon value.
for the inertial observers. The results are depicted as Figure 7, where we indeed qualitatively recover the same straight lines as we did in Figure 6. Note that in this case the lattice is not translationally invariant, as reflected in the very different curves away from the continuum limit.
The most interesting example is the condensate C ren I (x) , depicted in the top left plot, since it does not limit to zero. In this case we have plotted results both at the Unruh temperature β U and at temperatures differing from the latter by factors of 0.9 and 1.0, which thus do not yield the Hawking-Hartle vacuum. These, indeed, visibly diverge as the horizon is approached with x → −∞, while the β U curve remains smooth, and closely agrees with the value calculated in the Minkowski metric.

B. Schwarzschild: Hawking-Hartle Vacuum
Now we repeat the calculation in Schwarzschild spacetime in tortoise coordinates (42). The Hawking-Hartle vacuum is now no longer the Minkowski vacuum, although it approaches the latter far from the black hole x → ∞. Again, all operators but the condensate vanish in the continuum limit, albeit nontrivially (the small bump in J 0 (x) around x = 0 shrinks with the smallest value of r used in the extrapolation). Note the full Here we take α dS = 1 2π . Once again, all operators except C ren I (x) fall off to zero, and once again, the top-left plot of the latter shows results both at and slightly away from the Unruh temperature β U . Here we used δn = 2. Notice the slight asymmetry about x = 0 of the unconverged data, due to r always being added to x (rather than subtracted). current is conserved, J µ (x) ;µ = 0. For a time-independent metric the stipulation that time derivatives ought to vanish as they must in the Hawking-Hartle vacuum, J 0 (x) ,0 , implies J 1 (x) = 0 as well.
The top-left plot again shows the condensate at temperatures at and near the Hawking temperature β H = 8πM . Far from the black hole, results approach the inertial Minkowski results at this temperature. Close to it, they approach the Rindler results for a horizon with the same surface gravity. A x around 0, the Hawking-Hartle vacuum smoothly interpolates between these asymptotes.

C. de Sitter: Bunch-Davies Vacuum
We finally turn to the exponentially expanding de Sitter universe, whose static patch is covered by Ω dS static (x) = sech x α . The Hawking-Hartle vacuum in this case is called the "Bunch-Davies vacuum". In this case results are symmetric (antisymmetric, for J 0 ), so we show only those from x → −∞ to x = 0. Again, everything except the condensate vanishes in the field theory limit, whereas the condensate smoothly approaches the horizon only at the Unruh temperature.

V. CONCLUSION
Let us review before parting ways. We have illustrated, by simple computations using free Dirac fermions, a method to estimate two-point functions and quadratic expectation values of quantum expectation values upon curved 1+1 dimensional backgrounds. To get finite expectation values we made use of Hadamard renormalization, which perhaps unexpectedly can be applied directly to lattice data, provided the latter be computed holding the lattice length and separation constant. We used these results to demonstrate the Unruh effect, and to compute expectation values in de Sitter spacetime. This work opens the way for numerous future studies: one could quite straightforwardly study time-dependent physics, for example, as well as interacting theories once matrix product states are employed.

VI. ACKNOWLEDGMENTS
A. G. M. Lewis is supported by the Tensor Network Initiative at Perimeter Institute. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. G. Vidal is a CIFAR fellow in the Quantum Information Science Program. X is formerly known as Google[x] and is part of the Alphabet family of companies, which includes Google, Verily, Waymo, and others (www.x.company).
In light of (A7), taking the coincidence limit of both sides of (A6) yields which, inserted back into (A1), yields the desired result (30).

Appendix C: Spinor Components
For the sake of simplicity, in the main text we have mostly suppressed the spinor indices to lighten an otherwise already heavy notation. In this Appendix we will rewrite a few key expressions explicitly in terms of spinor components, given the important role the two spinor components play in the staggered lattice discretization.