Partial Differential Equations and Quantum States in Curved Spacetimes

: In this review paper, we discuss the relation between recent advances in the theory of partial differential equations and their applications to quantum ﬁeld theory on curved spacetimes. In particular, we focus on hyperbolic propagators and the role they play in the construction of physically admissible quantum states—the so-called Hadamard states —on globally hyperbolic spacetimes. We will review the notion of a propagator and discuss how it can be constructed in an explicit and invariant fashion, ﬁrst on a Riemannian manifold and then on a Lorentzian spacetime. Finally, we will recall the notion of Hadamard state and relate the latter to hyperbolic propagators via the wavefront set, a subset of the cotangent bundle capturing the information about the singularities of a distribution.


Introduction
Partial differential equations or PDEs pervade both mathematics and physics so exhaustively that it is hard to tell which of these two disciplines they more naturally pertain to. It is thus only natural that PDEs serve as a channel of intensive transfer of knowledge and innovation between mathematics and physics. However, the growing specialization of research communities results in a divergence of terminology and perspectives, which begins to impede the free exchange of ideas and methods. The mathematical physics community therefore has the crucial role of catalyzing and managing the effective communication between various mathematical and physical subjects, including those relying on partial differential equations. The present expository paper is a modest contribution to these efforts, and addresses an interplay between microlocal analysis of PDEs on manifolds on one side and mathematically rigorous quantum field theory or QFT on the other side.
In general terms, the problem at hand is as follows. Solutions of PDEs on manifolds play an important role in QFT by representing the quantum states of a physical system. One of the remarkable difficulties of QFT is the fact that these solutions may not be differentiable or even continuous, but are only assumed to be distributions (generalized functions). The precise way in which these solutions fail to be infinitely differentiable-their singularity structure-becomes an important criterion of how physically relevant a given mathematical solution is. Constructing solutions with the desired singularity structure in a more or less explicit manner, e.g., as certain integrals over the phase space, is a very challenging mathematical task. Until recently, such constructions existed only locally, that is, the integral representations of these solutions were valid only in small portions of the entire space and time. There are, however, new developments in microlocal analysis that allow integral representations that are global, i.e., valid on the entire spacetime, and in the present exposition we will describe the procedure of constructing the desired solutions as global oscillatory integrals. This will be done step-by-step, beginning with propagators in pure mathematics and ending with the construction of physically relevant quantum states in QFT.
While PDEs come in all kinds and flavors, here we will be concerned only with very special linear elliptic and hyperbolic PDEs with variable coefficients on compact curved manifolds. More precisely, our subject of study will be Laplacian and wave-type second order partial differential operators associated with Riemannian metrics on closed manifolds. We will discuss the construction of approximate solutions to the Cauchy problem for these wave-type operators as a single Fourier integral operator global both in space and time, which is a rather new development in the subject, while the classical approaches to this problem were guaranteed to work only locally. Then, we will switch the perspective from integral operators to their kernels, viewed as distributions, and discuss the structure of their singularities. Finally, we will touch upon the problem of the choice of an appropriate vacuum state in quantum field theory in the presence of gravity and will show that it is related to the singularity structure of certain distributions on the spacetime. Our global construction of approximate solutions to the wave-type equation will conveniently yield physically reasonable quantum states. All terms and constructions mentioned above that go beyond the undergraduate mathematical curriculum will be duly defined or given references in the sequel.
The mathematical language used in this paper will mostly come from spectral theory and microlocal analysis, with elements of basic Riemannian and Lorentzian geometry necessary to introduce the geometrical setting. In the discussion of quantum field theory, reference to physics terminology will be inevitable, but we will do our best to remain reasonably self-contained. The paper is intended for the professional mathematician, whether a specialist in PDEs or mathematical QFT or neither, who is interested in a first acquaintance with some of the most recent interactions between the two subjects. The reader is thus assumed to be familiar with the standard upper-division undergraduate mathematical curriculum, such as the basics of functional analysis and differential geometry.
The structure of the paper is as follows. In Section 2, we introduce propagators for wave-type hyperbolic operators ∂ 2 t + D where D is an elliptic operator on a closed Riemannian manifold. We first provide motivation through very basic examples in Section 2.1, and then describe the construction of propagators with all details in Section 2.2. In Section 3, the notion of a wavefront set is explained together with some of the crucial properties and relevant facts. Finally, Section 4 introduces the basics of QFT in curved spacetimes. Section 4.1 gives the generalities of Hadamard states, while their construction is accomplished in Section 4.2.

Hyperbolic Propagators
Let us begin with discussing wave-type hyperbolic equations and their solutions. Here, we will be concerned with equations of the kind (∂ 2 t + D) f = 0, where D is a time-independent second order positive elliptic operator (e.g., Laplacian, Schrödinger). The methods we discuss will work for all such equations under mild conditions, but for definiteness we will concentrate on the case D = −∆ g + m 2 , where ∆ g is the Laplace-Beltrami operator associated to a Riemannian metric g, and m 2 ≥ 0 is a constant. This equation is referred to as the Klein-Gordon equation on ultrastatic spacetimes in physics, and the parameter m is interpreted as the mass of rest of the quantum particle this equation describes. The reader interested in more general wave-type equations may consult [1][2][3][4].

Motivation
The simplest hyperbolic equation commonly used for a first introduction to the subject is the wave equation on the real line R, with initial conditions which admits the d'Alembert solution Unfortunately, the above explicit formula does not have analogues in higher dimensions, and the solution to the initial value (Cauchy) problem has to be found in less explicit ways using the Euclidean Fourier transform.
Consider the wave equation on R d with initial conditions Here x α , α = 1, . . . , d, are Euclidean coordinates in R d . The second order partial differential Equation (3) can be factorized as the product of two first order pseudo-differential equations Using Fourier transform methods, it is easy to see that one can represent the solution operators of the two equations corresponding to the factors on the RHS of (4) by means of oscillatory integrals as On account of the symmetry under time reflection, in order to solve the wave equation, it is sufficient to know the operator U(t) := e −it √ −∆ . The general solution to (3) can then be written as Of course, the above formula reduces to the d'Alembert solution (2) when d = 1. This connection helps to see that U(t) can be viewed as a propagator, which in d = 1 simply shifts the argument by t, but acts in a more complicated way for d > 1.
When working on a Riemannian manifold, the wave equation becomes a partial differential equation with variable coefficients, and the Fourier transform does no longer work. However, the spirit of the above argument still stands: solving the wave equation can be reduced to constructing the propagator U(t), and the latter can be expressed as an oscillatory integral that appropriately generalizes (5).

The Wave Propagator on a Riemannian Manifold
The appropriate generalization of (5) to a Riemannian manifold (M, g) warrants a moment of thought. The natural generalization of the Laplacian operator ∆ is the Laplace-Beltrami operator ∆ g , but what can substitute Fourier transform? If in R d one interprets the Fourier transform as the spectral resolution of the self-adjoint operator ∆, then the spectral resolution of ∆ g should be the analogue of Fourier transform in the Riemannian setting. For a non-compact M, the spectrum of ∆ g will not be discrete, bearing many technical difficulties when dealing with eigenfunction expansions. Further sophistication will arise if we also allow M to have a boundary, on which appropriate boundary conditions should be set in order for ∆ g to be self-adjoint. For the sake of simplicity and completeness of arguments, in what follows, we will assume that M is compact (so that the spectrum of ∆ g is discrete) and without boundary (so as to avoid dealing with boundary conditions).
Let ∆ g be the Laplace-Beltrami operator on scalar functions over M corresponding to the metric g. The operator −∆ g is elliptic, self-adjoint and non-negative. Its spectrum is discrete and accumulates to +∞. We adopt the notation −∆ g v k = λ 2 k v k for eigenvalues and corresponding orthonormalized eigenfunctions of −∆ g , where eigenvalues are enumerated with account of their multiplicity as The operator is called the wave propagator (Note that in the literature the operator U(t) is sometimes referred to as a half -wave propagator). The propagator (6) is a one-parameter family of unitary operators that solves, in a distributional sense, the operator-valued hyperbolic Cauchy problem   where Id is the identity operator on scalar functions. Arguing as in the previous section, the operator (6) can be used to write down the solution for the full wave equation supplemented with two initial conditions. In order to construct the propagator (6) precisely, one needs to know all eigenvalues and eigenfunctions of −∆ g . Indeed, the functional calculus gives us where ρ g (x) = det g αβ (x) is the Riemannian density. For a general Riemannian manifold, eigenvalues and eigenfunctions of the Laplace-Beltrami are not known.
However, one can still construct the propagator approximately, modulo an integral operator with smooth integral kernel, using techniques from microlocal analysis.
The subject of microlocal analysis is another wonderful confluence of mathematics and physics, where the basic idea is that in order to fully understand a PDE, one has to look not only at coordinates, but also at momenta, that is, work in the phase space (cotangent bundle). The structure of singularities of solutions to PDEs in terms of microlocal analysis provides a connection with Hamiltonian mechanics and symplectic geometry, which can be seen as the materialization of the rather blurry correspondence between particle mechanics and field theory.
The microlocal construction of U(t), now classical, has been the subject of intense research since the middle of the twentieth century, see for example [5][6][7][8][9] (Volumes 3 & 4), [10][11][12]. The operator U(t) is written locally, in time and space, as the composition of oscillatory integrals whose phase functions and amplitudes are obtained by solving a hierarchy of ordinary differential equations. The classical construction is explicit, but presents a number of shortcomings: it is not invariant under change of local coordinates, it is local in space and it is local in time.
We will present here an alternative, geometric version of the classical construction, developed in [13] building upon the earlier results [14][15][16], which overcomes the abovementioned issues. The main idea is to use a distinguished complex-valued phase function with a non-negative imaginary part. The complexity allows one to achieve a global construction and circumvent topological obstructions, whereas prescribing the phase function makes the oscillatory integral invariantly defined. Put The function (7) is a nowhere vanishing smooth function on T M positively homogeneous in momentum ξ of degree 1. It can be viewed as a Hamiltonian on T M. The corresponding Hamilton's equations admit a unique solution defined for all times t ∈ R for each choice of initial condition (x * (0; y, η), ξ * (0; y, η)) = (y, η) ∈ T M.
The function (7) will play a crucial role in our construction, because it coincides with the principal symbol of −∆ g , This implies that singularities of solutions of the half-wave equation propagate along the flow defined by (8), see Section 3. It is not hard to see that the curves t → x * (t; y, η) are geodesics.
Definition 1 (Levi-Civita phase function). We define the Levi-Civita phase function to be the infinitely smooth function for x in a geodesic neighborhood of x * and continued smoothly elsewhere, in such a way that Im ϕ ≥ 0. Here, grad z stands for the gradient in the variable z and dist denotes the geodesic distance. Note that, equipped with a Riemannian metric, a manifold becomes a metric space.
The Levi-Civita phase function encodes information about the geometry of M and about the Hamiltonian flow of (7). This is formalized by the following lemma [13,15]. Lemma 1. The Levi-Civita phase functions (9) satisfies the following properties: This allows us to represent the propagator (6) explicitly as a single oscillatory integral, global in space and time, with phase function ϕ. Namely, χ is a cut-off function that serves the purpose of localizing the integration in a neighborhood of the orbit with initial condition (y, η) and away from the zero section, see [13] (Section 5); • the weight w is defined by where the branch of the complex root is chosen in such a way that The scalar function a(t; y, η) is called the amplitude of the oscillatory integral (10). We can represent a as an asymptotic expansion in components a −k positively homogeneous in momentum η of degree −k: Individual homogeneous components can be determined iteratively by solving a hierarchy of ordinary differential equations, as detailed in the algorithm below.
Step two. Construct a new oscillatory integral with x-independent amplitude b = b(t; y, η), coinciding with (10) up to an integral operator with infinitely smooth integral kernel. Such a procedure is called reduction of the amplitude. This can be done by means of special operators, as described below. Put All differentiations are applied to the whole expression to the right of them. The operator (12b) is well defined because the differential operators L α commute [13] (Lemma A.2). When applied to a function positively homogeneous in momentumn η, the operator S −k decreases the degree of homogeneity by k. Hence, denoting by a ∼ ∑ ∞ j=0 a 2−j the asymptotic expansion of the function a defined by (11), the homogeneous components of the symbol b are The operator S ∼ ∑ ∞ k=0 S −k is called the amplitude-to-symbol operator. It maps the x-dependent amplitude a to the x-independent symbol b.
Step three. Impose the condition that our oscillatory integral (10) satisfies the wave equation, modulo an integral operator with infinitely smooth kernel. This is achieved by solving transport equations obtained by equating to zero the homogeneous components of the reduced amplitude b, b l = 0, supplemented with initial conditions a −k (0; y, η; ) determined by imposing that at t = 0 our oscillatory integral (10) is, modulo C ∞ , the integral kernel of the identity operator. Formula (13) describes a hierarchy of ordinary differential equations in the variable t whose unknowns are the homogeneous components of the original amplitude a. (10) is to be interpreted in a distributional sense: if one multiplies the RHS of (10) by f (y) and integrates the result over M with respect to the variable y, one obtains [U(t) f ](x) modulo an infinitely smooth function of x and t.

Remark 1. Formula
Remark 2. The construction presented above can be adapted to cover the case of more general scalar operator, as well as of systems of partial differential equations, under suitable assumptions, see [17][18][19][20][21][22]. We should mention that propagators are an important tool in abstract spectral theory, as they encode asymptotic information on the spectrum of the elliptic operator that generates them, cf. [14,23,24].

The Notion of Wavefont Set and Propagation of Singularities
In theoretical physics, field theory (classical and quantum) deals with solutions to hyperbolic wave-type equations, and the singularity structure of these solutions has profound ramifications on the consistency of the theory itself. Microlocal analysis provides the necessary tools for studying solutions of PDEs, and the object in the phase space that encodes the information about the singularities of a given solution is the wavefront set.
In this section, we will briefly review the notion of a wavefront set of a distribution. Recall that a distribution u ∈ D (M) is a continuous linear functional where C ∞ 0 (M) and C ∞ (M) are equipped with the standard Fréchet topology, see [25]. Let us now temporarily work in Euclidean space R d . Given a compactly supported distribution u ∈ E (R d ), we define its Fourier transform as We have the following standard result, relating the smoothness of a distribution and the decay of its Fourier transform [9] (Volume 1 Lemma 7.1.3).
This result allows one to identify directions in the Fourier space responsible for the non-smoothness of a distribution. Recall that a subset One can show that multiplying a compactly supported distribution by a smooth compactly supported function does not increase its set of singular directions, namely that By making the support of ϕ smaller and smaller, one can single out the set of singular directions accounting for the non-smoothness of u "at x". In this spirit, one defines the set of singular directions of u at x as We are now in a position to introduce the notion of a wavefront set.
Definition 3 (Wavefront set). Let u ∈ D (R d ). We define the wavefront set of u to be the closed The wavefront set captures the full structure of the singularities of a distribution, and it is the basic object upon which microlocal analysis is built. The adjective 'microlocal' refers to the fact that one localizes singularities of a distribution not only in the position variable x, but also in the dual variable (momentum) ξ. The latter effectively shows how the singularity looks in different directions.
The notion of a wavefront set can be suitably extended to manifolds as follows.
Theorem 2. Let ψ : R d → R d be a diffeomorphism. Then, the wavefront sets of u and ψ * u-the pull-back of u along ψ-are related as Let us choose an atlas {(U α , ψ α )} α for M and let {χ α } α , ∑ α χ α = 1, be a partition of unity subordinated to our atlas. Let u ∈ D (M) be a distribution on M. Then, χ α u is a compactly supported distribution on M, which can be viewed as a distribution u α := (ψ −1 α ) * (χ α u) on R d via the chart map ψ α , for which the wavefront set was defined above.
Then, one defines the wavefront set of u by The set (15) The distribution v is called the Schwartz (integral) kernel of V.
The operator (6) can be viewed as a linear map whose integral kernel u(t, x, y) ∈ D (R × M × M) is given by the RHS of (10). One can show that the wavefront set of u(t, x, y) is bounded above by the stationary points of ϕ, namely, the set C ϕ := {(t, x; y, η) | ϕ η (t, x; y, η) = 0} .
Direct inspection of (9) shows us that, for the Levi-Civita phase function, This agrees with the fact, guaranteed by the theorem of propagation of singularities [9] (Volume 1 Section 8.3), that the wavefront set of u(t, x, y) is the union of lifts of geodesics to the cotangent bundle. This piece of information was encoded within the Levi-Civita phase function, in that the latter was built so as to satisfy properties (i)-(iv) from Lemma 1, see also [15] (Section 1).

Quantum Field Theory on Curved Spacetimes: Hadamard States
While a consistent theory of quantum gravity is still far from our reach, the interaction of classical gravity with the rest of quantum matter can be studied relatively reliably. Quantum field theory in curved spacetimes is a mathematically rigorous framework that describes the propagation of quantum fields on a curved spacetime governed by what are called the semiclassical Einstein's equations. Following the principles of General Relativity, the spacetime is described by a smooth four-dimensional (orientable and time-orientable, globally hyperbolic) Lorentzian manifold. Classical fields are described by smooth sections in vector bundles over the spacetime, while quantum fields are given by operator-valued distributional sections in the same bundles. The dynamics or equations of motion are described by hyperbolic (at least in some sense) PDEs that the fields must satisfy. The most widely studied hyperbolic equations are either of wave-type or Dirac-type, corresponding to integer-spin and half-integer-spin particles, respectively. Below, we will concentrate on the Klein-Gordon equation ( + m 2 ) f = 0, where f is a section in the line bundle (function or distribution) over the spacetime, and is the d'Alembert operator associated with the Lorentzian metric (the analogue of the Laplace-Beltrami operator in the Lorentzian setting). This equation describes the field associated with scalar (spin zero) particles. Much of what follows works for general wave-type and Dirac-type operators [1,26].

Hadamard States
Hadamard states play a distinguished role in the algebraic formulation of quantum field theory, see, e.g., [26][27][28] for a recent review.
Much of the theory behind Hadamard states, including our constructions below, depend on an important property of the spacetime (M, g) called global hyperbolicty. Going back to the archetypical wave Equation (1), a well-posed Cauchy problem can be set up if the initial data are given on a curve that intersects every inextendible causal curve (x(τ), t(τ)) exactly once (the curve is causal if |ẋ(τ)| ≤ |ṫ(τ)|). Such an initial curve will be called a Cauchy curve. If the initial curve on which initial data are given fails to be Cauchy, then either of the following two problems arise: (i) Some causal curves will intersect the initial curve more than once, in which case, for generic initial data, a single-valued solution will not exist. (ii) Some causal curves will not intersect the initial curve at all, in which case a portion of the plane will not be in the domain of dependence of the initial curve, and the values of the solution in that region will not depend on the initial conditions, resulting in non-uniqueness.
The situation in more general Lorentzian manifolds with wave-type equations is qualitatively the same. In order for the Cauchy problem to be well-posed, the spacetime needs to possess a codimension one hypersurface that intersects every inextendible (i.e., with maximal interval of existence) causal curve exactly once. Such a hypersurface is called a Cauchy hypersurface, and a spacetime admitting a Cauchy hypersurface is called globally hyperbolic. There are several other equivalent definitions of global hyperbolicity, such as the compactness of causal diamonds, which we will not discuss here [1,29]. It was recognized early on that a globally hyperbolic manifold M is topologically homeomorphic to R × M, where M is the Cauchy surface. However, that M can always be chosen to be an embedded submanifold, and that (M, g) is isometrically diffeomorphic to (R × M, −β 2 (t, x)dt 2 + h t ), where β is a positive smooth function and h t is a smoothly time-dependent family of Riemannian metrics, is a relatively new development [30]. The converse problem of establishing which β and h t yield a globally hyperbolic spacetime in this manner has been settled in [31].
The algebraic (i.e., by means of operator algebras, which is a branch of analysis) approach to quantum field theory on curved spacetime can be succinctly formulated for the Klein-Gordon field as follows. Consider a globally hyperbolic spacetime (M ∼ = R × M, g) and let C ∞ 0 (M) be the space of test functions (smooth, compactly supported) on M. The hyperbolic operator + m 2 is symmetric with respect to the natural inner product and possesses unique advanced and retarded Green's functions E ± : C ∞ 0 (M) → D which satisfy is the causal future/past of the spacetime region N ⊂ M [1] (Section 3.4). Let A be a complex unital * -algebra generated by the images of a continuous linear map . It is clear that ω 2 is a weak bi-solution of the Klein-Gordon equation, It is customary in physics literature to use function notation φ(x) for the operator valued distribution φ by formally requiring φ( f ) = φ, f for all f ∈ C ∞ 0 (M). In these notations, one writes ω 2 (x, y) = ω(φ(x)φ(y)).
We will work here with a particular class of quantum states for which the 2-point function fully determines the state.

Definition 4 (Quasifree state). A state ω :
satisfy ω n ( f 1 , . . . , f n ) = 0 for n odd, where Π denotes all possible partitions of the set {1, 2, . . . , n} into pairs While all quantum states defined as above make perfect sense from a mathematical viewpoint, an arbitrary state is, in general, too singular to be a 'good physical states'. When manipulating states to describe the physics of the system, the type of operations that one can perform is constrained by the singular structure of the states themselves. Now, the singular structure of a generic state can be pretty wild, and this quickly becomes a serious technical limitation. It is worth stressing that the singular structure of a state cannot be just anything, in that the wavefront set of its 2-point function is constrained (from above) by condition (16). The latter, however, still leaves quite some freedom in the choice of ω 2 .
Let M denote the Minkowski space (the space R 4 equipped with the Lorentzian metric −dt 2 + dr 2 1 + dr 2 2 + dr 2 3 ) and let us adopt the notation x = (t, r) ∈ M, t ∈ R, r ∈ R 3 . It is well known [54] (Equation (4)) that the 2-point function of the Minkowski vacuum is where · denotes the Minkowski product, J 1 is the Bessel function of first kind, and with ψ(z) := d dz ln(Γ(z)) being the digamma function. The key idea underpinning the definition of Hadamard states is to prescribe that their 2-point functions possess the same singularity structure as (17). In order to turn this into a mathematically precise statement, we need to introduce further definitions and notation.
T is any local time coordinate increasing towards the future, is a reference length scale making the argument of the logarithm dimensionless, u, v n ∈ C ∞ (O × O) are uniquely determined by the geometry and the equation of motion (independently of ω) and w ∈ C ∞ (O × O) is determined by the state ω.
Definition 5 warrants a number of remarks.

Remark 3.
(a) The limit in the RHS of (18) has to be understood in the sense of distributions: first, one integrates against a test function, then one takes the limit for → 0 + . (b) The smooth functions u, v n ∈ C ∞ (O × O) are known as Hadamard coefficients. They are obtained as unique solutions of a hierarchy of differential equations that arise by imposing that the RHS of (18) solves the Klein-Gordon equation in the variable x, interpreting y as a parameter and setting w = 0. See, e.g., [38] (Appendix A) for further details. (c) Definition 5 immediately raises the question: does the series on the RHS of (18) converge?
The answer, unfortunately, is negative. The convergence of the series is only guaranteed when (M, g) is analytic. In the general smooth case, the series appearing in (18) has to be understood as an asymptotic expansion 'in smoothness', namely, the identity (18) means for every N = 1, 2, . . . . However, if one wants to work with a uniformly convergent series, the issue of non-convergence can be circumvented as follows. Choose a smooth cut-off χ : R → [0, 1], Then, there exists a real sequence  Although Definition 5 is rooted in good physical intuition, it is not very easy to handle or convenient to check. An alternative equivalent definition, based on the notion of wavefront set, was given in the now classical paper [56] by Radzikowski.
Here, (x 1 , ξ 1 ) ∼ (x 2 , ξ 2 ) means that there exists a lightlike geodesic γ connecting x 1 and x 2 , such that ξ 1 is cotangent to γ at x 1 and ξ 2 is the parallel transport of ξ 1 from x 1 to x 2 along γ with respect to the Levi-Civita connection. The notation ξ 1 0 means that ξ 1 is future directed.

Remark 4.
Loosely speaking, the microlocal spectrum condition means that-pointwise in the cotangent space-the 'Fourier transform' of the 2-point function is rapidly decaying everywhere, but in the future light cone, and this holds in a consistent manner as we move the base point around. This effectively implements on curved spacetimes the usual positive frequency condition of QFT. The generalization of this equivalence result to more general vector-valued fields (together with a minor technical fix in the original proof) was given in [57].
The two equivalent definitions of Hadamard states-Definition 5 and Theorem 4are inherently local. While there is no clear way of defining a global object starting from Definition 5, one can use the microlocal spectrum condition to do so, by simply dropping the restriction to a convex neighborhood in (19). We call the resulting object a global Hadamard state.
The relation between local and global Hadamard states is clarified by the following two theorems.
Theorem 5 (Radzikowski's local-to-global theorem [58]). If a quasifree state ω on a globally hyperbolic spacetime (M, g) is locally of Hadamard form, then it is a global Hadamard state.
Theorem 6 (Propagation of the Hadamard property [59]). Let (M, g) be a globally hyperbolic spacetime and let ω be a quasifree state. Suppose that the microlocal spectrum condition (19) holds for ω in a neighborhood of a Cauchy surface M. Then it holds globally.
Finally, we conclude this section with a result due to Fulling, Narcowich and Wald [60], which shows that one always has Hadamard states to work with. It is important to point out that Hadamard states are far from being unique on any globally hyperbolic spacetime. As it is apparent from the definition, the Hadamard condition only constrains the singularity structure of the 2-point function ω 2 of a state ω, and any other state ω whose 2-point function ω 2 differs from ω 2 by a smooth function has the same singularity structure and is therefore another Hadamard state. In fact, on every globally hyperbolic spacetime, there are infinitely many Hadamard states, and in a sense, all reasonable quantum states are expected to be of the Hadamard form. Moreover, Hadamard states are closed under quantum field theoretical operations (see, e.g., [61,62]); once your chosen vacuum state is Hadamard, then the entire corresponding folium of states is such. Thus, the actual physical content of the state is, roughly speaking, encoded in the smooth rather than in the singular part. The Hadamard condition, essentially, rules out pathological states that would prevent the reasonable application of physically motivated mathematical operations. In particular, vacuum states are widely assumed to be Hadamard. While there have been alternative proposals for the definition of the class of physically relevant quantum states, Hadamard states seem to be the most popular and best justified of all [63]. Now, the selection of a vacuum state among infinitely many Hadamard states is another big problem in quantum field theory on curved spacetimes. By physical intuition, a vacuum state should minimize the energy in some sense, but in quantum field theory, the energy density is not bounded from below pointwise. In the Minkowski spacetime, the Minkowski vacuum defined earlier is singled out as the canonical choice of a vacuum state due to its unique feature of being invariant under all spacetime symmetries. A generic curved spacetime, on the other hand, may not have any symmetries, and thus infinitely many candidate vacuum states will stand on equal footing. There have been proposals in the literature for singling out a unique vacuum state, such as the so-called states of low energy, which minimize the energy density as measured by a hypothetical observer [45,64]. However, this definition is very observation-dependent and the general question of the choice of a vacuum state in a curved spacetime remains open.

Construction of Hadamard States
The construction of the wave propagator can be suitably extended to Lorentzian manifolds with compact Cauchy surface. When the manifold carries a globally hyperbolic Lorentzian, as opposed to Riemannian, metric, greater care is needed, as time can no longer be treated as an external parameter and the components of the metric tensor explicitly depend on it. As a result, constructing an analogue of the propagator (6) in a global, invariant fashion presents some additional challenges arising from the spacetime geometry.
One can show [1,27] that singularities of solutions of the wave equation on globally hyperbolic spacetimes propagate along light-like geodesics. Now, the four-dimensional version of our Hamiltonian (7) vanishes identically on light-like covectors. Nevertheless, one can still construct a distinguished geometric phase function in the spirit of (1) as follows.
Let (M, g), M R × M, be a globally hyperbolic spacetime of dimension d. Let us denote by ι s : M s → M the embedding of M s {s} × M into M. Let Y = (s, y) ∈ M. For any η ∈ T * y \ {0}, let η + be the unique future pointing null covector in T * Y M such that ι * s η + = η. Here, and further, the upper * denotes the pull-back along the map it is applied to. Furthermore, let us define where η h s = h s αβ (y)η α η β .
The Levi-Civita flow is designed in such a way that it encodes information about propagation of singularities and X * and Ξ * are positively homogeneous in η of degree 0 and 1, respectively.
It is easy to see [19] (Lemma 4.3) that the Levi-Civita flow can be reparameterized using the global time coordinate t. With slight abuse of notation, we denote by t → ( X * (t; s, y, η), Ξ * (t; s, y, η)).
Starting from the Lorentzian Levi-Civita phase function, one can set in motion machinery similar to that of Section 2.2. We refrain from discussing the technical details here, and we refer the interested reader to [19].
In the remainder of this section, we will explain how one can construct Hadamard states using these techniques in the special case of ultrastatic spacetimes. Observe that a conformal transformation turns a static spacetime into an ultrastatic one. Therefore, the following immediately applies to the class of static spacetimes.
An explicit formula for the scalar symbol f can be obtained by constructing e −i(t−s) √ −∆ g first, and then composing the outcome with the parametrix −∆ g . Alternatively, one can construct (20) in one go by replacing the initial condition U(0) = Id from Section 2.2 with the condition Ω(s, s) = 1 2 −∆ g . Remark 5. From the above discussion, the reader may be led to think that there are no fundamental differences between studying partial differential equations on Riemannian manifolds and doing so on Lorentzian ones, provided one is prepared to take care of the amount of additional technicalities that the latter brings about. While this is true in some instances, it is very much not the case in general. Indeed, fully relativistic equations of mathematical physics are not always associated with a natural inner product, not even an indefinite non-degenerate one. As a result, in the Lorentzian setting one is often forced to work with equations, as opposed to operators, see [68,69].
Whilst a systematic review of the subject goes beyond the scope of this paper, we should conclude by mentioning that a number of different approaches to construct Hadamard states, alternative to that presented here, have been proposed over the years, see, for example, [60] or [70,71]. The closest to the above in terms of techniques employed are perhaps those of Gérard and Wrochna [72,73], heavily reliant on pseudo-differential calculus.

Conclusions
The interplay between partial differential equations and mathematical physics has, over the years, proved extremely fruitful in both directions. In this review, we have discussed a very specific aspect of this relationship, namely, how the construction of propagators for hyperbolic PDEs as global oscillatory integrals is useful in describing physically meaningful quantum states in curved spacetimes, the so-called Hadamard states. The results presented here are relatively recent, and we believe that this expository paper, which brings together theory and applications in a self-contained manner, will help to make some progress towards establishing a common language and, ultimately, bridging different research communities working in partial differential equations and mathematical physics.
The topics presented here open the way to future exciting avenues of research. Interesting questions are, for example: Can one perform a similar global construction for systems-say, for the Dirac field? Can one perform a similar global construction on globally hyperbolic manifolds with non-compact Cauchy surface and with boundary? Or, finally, can one apply the above techniques to construct Feynman propagators? These generalizations are expected to be related to remarkable technical complications and may necessitate the introduction of additional mathematical tools.
Author Contributions: All authors contributed equally. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.