\section*{Introduction}
If $(M,g)$ is a globally hyperbolic spacetime then the Cauchy problem for normally hyperbolic operators, i.e., operators of the form 
\begin{equation*}
	\Box = - g^{ij} \partial_i \partial_j + \textrm{lower order terms}, \quad i,j = 1, \ldots, \dim M 
\end{equation*}
(in local coordinates using Einstein's sum convention)
is well posed. 
This implies the existence of unique advanced and retarded fundamental solutions (propagators). 
The properties of these fundamental solutions are extremely important for the understanding of classical wave propagation, such as the electromagnetic waves. 
They appear naturally because of causality: the retarded fundamental solution propagates to the future, whereas the advanced fundamental solution propagates to the past. 
In quantum field theory the appearance of time-ordering and the enforcement of positivity of energy have led to the development of another type of fundamental solution, the Feynman propagator.  
It propagates positive energy solutions to the future and negative energy solutions to the past and thus combines causality with the notion of positivity of energy.
Explicit formulae for Lorentz invariant Feynman propagators for normally hyperbolic operators in Minkowski spacetime can be found in basic quantum field theory textbooks. 
The usual construction employs the Fourier transform. 
However, on a generic globally hyperbolic spacetime there is no Fourier transform or reasonable notion of energy in the absence of a global timelike Killing vector field. 
A priori, it is therefore unclear what a Feynman propagator should be.

In the theory of partial differential equations, the notion of a parametrix is often useful in the first stages on the  construction of a true fundamental solution. 
A parametrix is, by definition, an inverse of the operator $\Box$ modulo smoothing operators. 
Parametrices are considered equivalent if they differ by smoothing operators.
It was a deep insight of 
Duistermaat and H\"{o}rmander~\cite{Duistermaat_ActaMath_1972}   
that there is a well defined notion of Feynman parametrices and these parametrices are unique up to smoothing operators, in other words they are unique as parametrices. 
In fact, the notion of distinguished parametrices for general pseudodifferential operators of real principal type is given in their seminal article where they have identified a geometric notion of pseudoconvexity which allows to prove uniqueness of such parametrices. 
Roughly speaking, despite of not having any notion of energy due to the lack of a global timelike Killing vector field, there is still microlocal equivalence of positivity of energy and a corresponding flow on the cotangent bundle. 
This microlocal notion can be used to characterise Feynman parametrices: they are distinguished by their wavefront sets rather than the support properties of the Fourier transform.

Feynman parametrices play an extremely important role in quantum field theory on curved spacetimes, and the theory has actually been developed to a certain extent first in the physics literature.
Canonical quantisation of linear fields can be done in two stages. 
In the first step, one constructs a field algebra from the space of solutions of the 
respective equation of motion~\cite{Dimock_CMP_1980, Dimock_AMS_1982} 
and next, some state is required to construct a Hilbert space representation of this field algebra. 
Whereas, in Minkowski spacetime these two steps are usually combined into one owing to the existence of the vacuum state, it is more fruitful to separate them in curved spacetimes. 
The first step can be done without any problems, exactly the same way as in Minkowski spacetime but the second step necessitates a notion of a reasonable state. 
It has been realised that a state compels certain conditions in order to perform the usual operations in 
perturbative quantum field theory~\cite{Brunetti_CMP_2000, Hollands_CMP_2001, Hollands_CMP_2002}. 
One of the identified conditions is a restriction of the type of singularity that one obtains from the state, the so called Hadamard condition. 
Although Duistermaat and H\"ormander were certainly aware of the developments in physics, it was realised only much later by 
Radzikowski~\cite{Radzikowski_CMP_1996} 
that the expectation value of the time-ordered products with respect to states satisfying the Hadamard condition are Feynman propagators. 
In fact, the construction of a Hadamard state is equivalent to the construction of a Feynman propagator that satisfies a certain positivity property. The fact that this positivity property holds for parametrices was already shown by Duistermaat and H\"{o}rmander. The property for parametrices implies the existence of such Feynman propagators and hence of Hadamard states (see~\cite{Lewandowski}).

Apart from its applications in physics Feynman parametrices also play an increasingly important role in mathematics.
For instance, they appear in the context of Lorentzian index theory as inverses modulo compact operators of the Dirac operator with Atyiah-Patodi-Singer boundary conditions~\cite{Baer_AJM_2019} 
and Lorentzian local index theorem~\cite{Baer}.  
These concepts also arise in Vasy's treatment of asymptotically hyperbolic problems, see for example~\cite{GellRedman_CMP_2016}  
and references therein, where even a nonlinear problem is discussed in this context.
Moreover, the notion of Feynman propagator also comes up naturally in the Lorentzian generalisation of the Duistermaat-Guillemin-Gutzwiller trace formula~\cite{Strohmaier_AdvMath_2021}. 

So far, the Duistermaat and H\"ormander's direct construction of distinguished parametrices for pseudodifferential operators of real principal type has not been used much in the physics literature probably due to its great generality and the associated complex notation. The main idea is however compelling and simple. The operator can be conjugated microlocally to a vector field and it is therefore sufficient to construct a parametrix  for the operator of differentiation on the real line. There are two distinguished fundamental solutions for the operator of differentiation. A choice of parametrix for each connected component of the characteristic set of the operator results in a distinguished parametrix. Thus, there are $2^N$ distinguished parametrices if $N$ is the number of connected components of the characteristic set. For the wave operator this gives four distinguished parametrices in dimensions $n \geq 3$ and $16$ distinghuished parametrices in dimension two.

The aim of this article is three fold. 
First we would like to revise and simplify the construction of 
Duistermaat-H\"ormander~\cite{Duistermaat_ActaMath_1972} 
in the special case of normally hyperbolic operators. 
The second aim is to fill a gap in the literature:  the construction of distinguished parametrices in Duistermaat and H\"ormander is given only for scalar operators. 
Several constructions in index theory, in trace formualae and also in physics need the existence and uniqueness of Feynman parametrices for operators acting on vector bundles. 
It is known that most of the constructions carry over to the case of normally hyperbolic operators, since their principal symbols are scalars. 
There are, however, also important differences that appear on the level of subprincipal symbols. 
In this article we would like to provide a detailed construction of Feynman parametrices for vector bundles, with complete proofs for their uniqueness, and discuss the effect of curvature of the bundle connection. In case the vector bundle has a hermitian inner product with respect to which the operator is formally selfadjoint we show that the construction can be carried out in such a way that the above mentioned positivity property holds. 
Finally and thirdly, for Dirac operators we give a much more direct construction of Feynman parametrices satisfying a positivity property.
We also discuss some consequences, including the usual propagation of singularity theorems.
These are well known to hold for very general operators on vector bundles. 
For instance, the propagation of polarisation sets has been proven by Dencker using 
microlocalisation in the matrix setting~\cite{Dencker_JFA_1982}. 
 
As explained above the construction of Feynman parametrices satisfying the positivity property is equivalent to the construction of Hadamard states.
There are by now several constructions of Hadamard states, recently even in the 
analytic category~\cite{Gerard_CMP_2019}. 
Amongst the methods to construct them, there are direct ones using singularity expansions employing the 
Lorentzian distance function 
(the so-called Hadamard expansion)~\cite{Fulling_CMP_1978, Brown_JMP_1984, Sahlmann_RMP_2001, Marecki_master, Hollands_RMP_2008, Lewandowski},  
spectral methods that rely on frequency splitting 
(deformation method)~\cite{Fulling_AnnPhys_1981}, 
pseudodifferential methods~\cite{Junker_RMP_1996, Gerard_CMP_2014, Gerard_CMP_2015} 
(see also the recent monograph~\cite{Gerard_EMS_2019}), 
holography~\cite{Moretti_CMP_2008, Dappiaggi_JMP_2009, Dappiaggi_ATMP_2011, Gerard_AnalPDE_2016} 
(see also the exposition~\cite{Dappiaggi_Springer_2017}), 
and even 
global methods~\cite{GellRedman_CMP_2016, Vasy_AHP_2017, Derezinski_RMP_2018, Capoferri_JMAA_2020}. 




\section{The setup and the main results} \label{sec: setup}
Suppose that $(M,g)$ is an $n$-dimensional  globally hyperbolic spacetime, i.e., an oriented and timeoriented Lorentzian manifold of metric-signature $(-1,+1,\ldots,+1)$ that possesses a smooth Cauchy hypersurface.  Let $E \to M$ be a smooth complex vector bundle. 
A differential operator $\Box$ acting on sections of $E$ is called \textit{normally hyperbolic} if its principal symbol $\symb{\square}$ is scalar and equals to the metric on the cotangent bundle, i.e.,  
\begin{equation*}
	\symb{\square} (\xi) := g^{-1} (\xi, \xi) \, \mathbbm{1}_{\End{E}}, \quad \forall \xi \in C^{\infty} (M; T^{*} M).  
\end{equation*} 
Then there is a unique connection
$\nabla: C^\infty(M; E) \to C^\infty(M; E \otimes T^*M)$ and a unique potential $V \in C^\infty(M; \End{E})$ such that $\Box$ takes the form 
(see e.g.~\cite[Prop. 3.1]{MR1418003})
\begin{equation*}
	\Box = -\tr_{g} \nabla^2 + V.
\end{equation*}
This formula is sometimes called the Weitzenb\"{o}ck formula and we will refer to the connection $\nabla$ as the Weitzenb\"{o}ck connection.
If $E$ carries a bundle metric such that $\Box$ is formally selfadjoint then $\nabla$ is compatible with this metric.
We do not however assume this, unless stated otherwise.

It is a classical result that normally hyperbolic operators on globally hyperbolic manifolds have unique retarded and advanced Green's operators $G^{\mathrm{ret, adv}}: C_{\mathrm{c}}^{\infty} (M; E) \to C^{\infty} (M; E)$ with the characteristic property   
\begin{equation*}
	\mathrm{supp} \big( G^{\mathrm{ret, adv}} (u) \big) \subset J^{\pm} (\supp{u}), \quad \forall u \in C_{\mathrm{c}}^{\infty} (M; E),  
\end{equation*}
where $J^\pm(K)$ denotes the causal future/past of a set $K \subset M$. 
Constructions in great generality are due to  
Hadamard~\cite{Hadamard_ActaMath_1908, Hadamard_Dover_2003} 
and 
Riesz~\cite{Riesz_ActaMath_1949, Riesz_CPAM_1960}. 
Contemporary expositions include, for 
example~\cite{Friedlander_CUP_1975, Guenther_AP_1988, Baer_EMS_2007},  
and we refer here, for instance~\cite[Cor. 3.1.4]{Baer_EMS_2007} or~\cite[Prop. 4.1, Rem. 4.3 (b)]{Guenther_AP_1988}  
for normally hyperbolic operators. 

Throughout this paper we will not distinguish notationally between a continuous linear map $C_{\mathrm{c}}^{\infty} \to \mathcal{D}'$ and its distributional kernel. 
This means that the fundamental solutions $\advGreenOp, G^{\mathrm{ret}}$ are also to be understood as distributions in $\mathcal{D}'(M \times M; E \boxtimes E^*)$. 
On general manifolds operators map will be thought of as maps between half-densities. 
This allows pairings between functions without using a volume form. 
In this way the space of distributions $\mathcal{D}' (M; E \otimes \varOmega^{1/2})$ is defined as the topological dual of $C_{\mathrm{c}}^{\infty} (M; E^{*} \otimes \Omega^{1/2})$, where $\varOmega$ is the bundle of densities. 
In case a metric is given on $M$ we use the metric volume form to identify half-densities with functions, i.e., in this case the pairing between distributional sections and sections is achieved using the Lorentzian volume form. 
In the general case of a manifold without given metric we will be notationally suppressing $\Omega$, so all functions and sections will be considered as half-densities on the manifold.
We will denote the zero section removed cotangent bundle by $\dot{T}^{*} M$ and use Einstein's summation convention.  

Our convention of Fourier transform is 
\begin{equation*}
	\hat{f} (\theta) := \int_{\Rn} f (x) \, \mathrm{e}^{ - \mathrm{i} x \cdot \theta} \rd x, 
	\qquad 
	f (x) = \frac{1}{(2 \pi)^{n}} \int_{\Rn} \hat{f} (\theta) \, \mathrm{e}^{ \mathrm{i} x \cdot \theta} \rd \theta 
\end{equation*}
for any function $f \in L^{1} (\Rn, \rd x)$ and $\cdot$ is the Euclidean inner product. 
Suppose that $U \subset \Rn$ is a nonempty open set.
Recall that $\dot{T}^{*} U \ni (x_{0}, \xi^{0}) \notin \WF{u}$ is not in the wavefront set of a distribution $u \in \mathcal{D}' (U; E)$ if and only if there exists a compactly supported section $f$ on $U$ nonvanishing at $x_{0} \in U$ such that the Fourier transform of $\widehat{fu}(\xi)$ is rapidly decreasing in a conic neighbourhood of $ \xi^{0}$~\cite[Prop. 2.5.5]{Hoermander_ActaMath_1971}. 
This definition in fact makes sense on smooth manifolds as the wavefront set transforms covariantly under a change of coordinates. 
Details are available, for instance, in~\cite[Chap. VIII]{Hoermander_Springer_2003} 
or in the 
exposition~\cite{Strohmaier_Springer_2009} 
customised for quantum field theory on curved spacetimes. 
If $u \in \mathcal{D}' (M; E)$ then $\WF{u} \subset \dot T^* M$.
The \textit{twisted wavefront set} of a bidistribution $A \in \mathcal{D}'(M \times M ; E \boxtimes F^*)$ is given by 
\begin{equation*}
	\WFPrime{A} := \{ (x_{0}, \xi^{0}; y_{0}, - \eta^{0}) \in \dot{T}^{*} M \times \dot{T}^{*} M ~|~ (x_{0}, \xi^{0}; y_{0}, \eta^{0}) \in \WF{A} \}
\end{equation*} 
where $F \to M$ is a smooth complex vector bundle.  
The twisted wavefront  $\WFPrime{P}$ of the Schwartz kernel of a pseudodifferential operator $P$ is a subset of the diagonal in  $\dot T^* M \times  \dot T^* M$.  Then $\ES{P} := \{ (x_{0}, \xi^{0}) \in \dot{T}^{*} M \mid \; (x_{0}, \xi^{0}; x_{0}, \xi^{0}) \in  \WFPrime{P} \}$ is the 
\textit{essential support}\footnote{It 
	is also known as the \textit{microsupport} or \textit{wavefront set} of a pseudodifferential operator.
} 
of $P$. 
This is the smallest conic set such that $P$ is of order $-\infty$ in $\dot{T}^{*} M \setminus \ES{P}$ 
(see, e.g.~\cite[Prop. 18.1.26]{Hoermander_Springer_2007}). 

Throughout the article, only the polyhomogeneous symbol class $S^{m}$ 
(Definition~\ref{def: polyhomogeneous_symbol_mf}) 
will be used, so $\Psi \mathrm{DO}^{m}$ will be the set of pseudodifferential operators having polyhomogeneous total symbols.  
We will denote by $H_{\mathrm{loc}}^{s} (M; E)$, the space of sections that are locally in the Sobolev space $H^{s}$ with respect to any smooth chart and smooth bundle chart. The space of sections in  $H_{\mathrm{loc}}^{s} (M; E)$ of compact support is denoted by $H_{\mathrm{c}}^{s} (M; E)$ . As usual the space $H_{\mathrm{loc}}^{s} (M; E)$
is equipped with the locally convex topology of convergence locally in $H^s (M; E)$. The space $H_{\mathrm{c}}^{s} (M; E)$ is the  union $\bigcup_{K \Subset M} H_{\mathrm{c}}^{s} (K; E)$ and is equipped with the inductive limit topology, where the union runs over all compact subset $K$ of $M$. 
For details, we refer, for example~\cite[App. B1]{Hoermander_Springer_2007}. 
\subsection{Feynman parametrices and propagators}
A parametrix which is uniquely characterised by its wavefront set is called a 
distinguished parametrix~\cite[Sec. 6.6]{Duistermaat_ActaMath_1972}. 
\begin{definition}[Feynman parametrices] 
\label{def: Feynman_parametrix}
	Let $(M, g)$ be a globally hyperbolic spacetime whose punctured cotangent bundle is denoted by $\dot{T}^{*} M$. 
	The Feynman $G^{\mathrm{F}}$ and the anti-Feynman $\antiFeynParametrix$ parametrices for a normally hyperbolic operator acting on smooth sections of a vector bundle $E \to M$ are parametrices $G^{\mathrm{F}}, \antiFeynParametrix : C_{\mathrm{c}}^{\infty} (M; E) \to C^{\infty} (M; E)$ whose Schwartz kernel satisfies~\cite[p. 229]{Duistermaat_ActaMath_1972} 
	(see also~\cite[pp. 541-542]{Radzikowski_CMP_1996}) 
	\begin{equation*}
		\WFPrime{G^{\mathrm{F}}} \subset \varDelta^{*} \cup C^{+}, 
		\qquad 
		\WFPrime{\antiFeynParametrix} \subset \varDelta^{*} \cup C^{-}  
	\end{equation*}
	where $\varDelta^{*}  := \{ (x, \xi; x, \xi) \in \dot{T}^{*} M \times \dot{T}^{*} M \}$ is the diagonal in $\dot{T}^{*} M \times \dot{T}^{*} M$,  
	\begin{equation} \label{eq: def_Feynman_bicharacteristic_relation_NHOp}
		C^{\pm} := \big\{ (x, \xi; y, \eta) \in \dot{T}^{*} M \times \dot{T}^{*} M ~|~ g_{x}^{-1} (\xi, \xi) = 0, \exists s \in \mathbb{R}_{\gtrless 0} : (x, \xi) = \varPhi_{s} (y, \eta) \big\}  
	\end{equation}
	are the forward and backward geodesic relations on $\dot{T}^{*} M \times \dot{T}^{*} M$, and $\varPhi_{s}$ is the time $s$ geodesic flow on $\dot{T}^{*} M$. 
\end{definition}
In this exposition we give a vector bundle version of the classic result by 
Duistermaat-H\"{o}rmander~\cite[Thm. 6.5.3]{Duistermaat_ActaMath_1972} 
(see also~\cite[Thm. 26.1.14]{Hoermander_Springer_2009}) 
for normally hyperbolic operators. 
\begin{theorem}[Existence and uniqueness of Feynman parametrices] 
\label{thm: exist_unique_Feynman_parametrix_NHOp}
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$. 
	Then there exist unique (modulo smoothing kernels) Feynman $G^{\mathrm{F}}$ and anti-Feynman $\antiFeynParametrix$ parametrices for any normally hyperbolic operator $\square$ on $E$. 
	Furthermore, for every $s \in \mathbb{R}$, $G^{\mathrm{F}}$ and $\antiFeynParametrix$ extends  continuous maps from  $H_{\mathrm{c}}^{s} (M; E)$ to  $H_{\mathrm{loc}}^{s+1} (M; E)$ with   
	\begin{equation*}
		G^{\mathrm{F}} - \antiFeynParametrix \in I^{-3/2} \big( M \times M, C'; \End{E} \big), 
	\end{equation*} 
	where $I^{-3/2} \big( M \times M, C'; \End{E} \big)$ is the Lagrangian distribution (Definition~\ref{def: Lagrangian_distribution_bundle_local}) 
	and $G^{\mathrm{F}} - \antiFeynParametrix$ is noncharacteristic 
	(cf. Definition~\ref{def: elliptic_FIO}) 
	at every point of the geodesic relation 
	\begin{equation} \label{eq: def_bicharacteristic_relation_NHOp}
		C := \{ (x, \xi; y, \eta) \in \dot{T}^{*} M \times \dot{T}^{*} M ~|~ g_{x}^{-1} (\xi, \xi) = 0, \exists s \in \mathbb{R} : (x, \xi) = \varPhi_{s} (y, \eta) \}.
	\end{equation}
\end{theorem}
A special case of this result for massive scalar wave operator was given by 
Radzikowski~\cite[Prop. 4.2, Prop. 4.3, Prop. 4.4]{Radzikowski_CMP_1996} 
as a consequence of the aforementioned result by Duistermaat and H\"{o}rmander. 
Recently, 
Lewandowski~\cite[Prop. 3.11]{Lewandowski} 
has given a direct construction for real vector bundles employing the Hadamard series expansion along the presentation of~\cite{Baer_EMS_2007}. 
The proof described in our paper is in contrast purely microlocal as in the original 
treatment~\cite[Thm. 6.5.3]{Duistermaat_ActaMath_1972}.
This only requires a bundle version of microlocalisation, as described in Proposition~\ref{prop: microlocalisation_P} and Proposition~\ref{prop: microlocalisation_NHOp}. This is along the lines of the proof of propagation of polarization sets as given by  Dencker~\cite{Dencker_JFA_1982} 
but our presentation is inherently geometric. 
Global Fourier integrals for distinghuished parametrices for scalar wave operators with time independent smooth potential in spatially compact globally hyperbolic ultrastatic manifolds have also been constructed by Capoferri \textit{et al}.~\cite[Thm. 5.2]{Capoferri_JMAA_2020}.

In fact, the following positivity property can be shown.
\begin{restatable}{proposition}{positivityFeynmanMinusAdvNHOp} 
\label{thm: positivity_Feynman_minus_adv_NHOp}
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$ equipped with a nondegenerate sesquilinear form such that the normally hyperbolic operator $\square$ on $E$ is formally selfadjoint. 
	Then there exists a Feynman parametrix $G^{\mathrm{F}}$ for $\square$ such that $W := - \mathrm{i} (G^{\mathrm{F}} - \advGreenOp)$ is formally selfadjoint. 
	Additionally, if this form is positive-definite (hermitian) then $W$ can be chosen nonnegative.
\end{restatable}
This is essentially a bundle version of that in 
Duistermaat-H\"{o}rmander~\cite[Thm. 6.6.2]{Duistermaat_ActaMath_1972} 
which was proven by deploying a refined microlocalisation of scalar pseudodifferential operators~\cite[Lem. 6.6.4]{Duistermaat_ActaMath_1972}. 
We provide such a microlocalisation for normally hyperbolic operator~\eqref{eq: microlocalisation_NHOp_refined}. 

Finally, one can turn the Feynman parametrix $G^{\mathrm{F}}$ into a Feynman propagator $\FeynGreenOp$ utilising the well-posedness of Cauchy problem for normally hyperbolic operators on any globally hyperbolic manifold. 
\begin{restatable}[Existence of Hadamard bisolutions]{theorem}{existenceFeynmanGreenOpNHOp} 
\label{thm: existence_Feynman_propagator_NHOp}
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$ equipped with a nondegenerate sesequilinear form such that the normally hyperbolic operator $\square$ on $E$ is formally selfadjoint. 
	Then there exists a Feynman propagator $\FeynGreenOp$ for $\square$ such that $\omega := - \mathrm{i} (\FeynGreenOp - \advGreenOp)$ is formally selfadjoint.  
	In addition, if this form is positive-definite (hermitian) then the Feynman propagator can be chosen such that $\omega (u^{*} \otimes u) \geq 0$ for any $u \in C_{\mathrm{c}}^{\infty} (M; E)$. Thus $\omega$ defines a Hadamard state. 
	Here $u \mapsto u^{*}:  E \to E^{*}$ is the fibrewise antilinear mapping induced by the sesquilinear form.
\end{restatable}
Duistermaat-H\"{o}rmander~\cite[p. 229]{Duistermaat_ActaMath_1972} 
have actually considered a much wider class of operators, namely scalar pseudodifferential operators of real principal type on a smooth manifold. 
The pivotal step of determining the appropriate smoothing operators required to obtain Feynman propagators from respective parametrices was however left open.
This indeterminacy can be fixed in various ways on special spacetimes even in the absence of the timelike Killing vector field. Such constructions have appeared in the literature on microlocal analysis. 
For example, 
Gell-Redman \textit{et al}.~\cite[Thm. 3.6]{GellRedman_CMP_2016} 
treat the scalar wave operator in spaces with non-trapping Lorentzian scattering matrices.  Vasy constructs Feynman propagators 
by making assumptions on global dynamics for 
(i) scalar formally selfadjoint operators with real principal symbol in closed manifolds~\cite[Thm. 1]{Vasy_AHP_2017}, 
and in spaces of Lorentzian scattering matrices for 
(ii) wave operators in Melrose's b-pseudodifferential algebraic framework~\cite[Thm. 7]{Vasy_AHP_2017} 
and 
(iii) Klein-Gordon operators in Melrose's scattering pseudodifferential algebraic formalism~\cite[Thm. 10, Thm. 12]{Vasy_AHP_2017}.
The idea is to identify the appropriate spaces where these operators are invertible and then define the Feynman propagators as the inverse of those operators satisfying the required properties - a generalisation of Feynman's original "$\pm \mathrm{i} \epsilon$" prescription. 
As a consequence, he has also achieved respective positivity property for Feynman parametrices~\cite[Cor. 5, Cor. 9, Cor. 11, Cor. 13]{Vasy_AHP_2017},  
as orginally announced by 
Duistermaat-H\"{o}rmander~\cite[Thm. 6.6.2]{Duistermaat_ActaMath_1972}. 
A special case of Theorem~\ref{thm: existence_Feynman_propagator_NHOp} for Klein-Gordon operator minimally coupled to a static electromagnetic potential on a static spacetime has been proven by 
Derezi\'{n}ski-Siemssen~\cite[Thm. 7.7]{Derezinski_RMP_2018}. 
In the spirit of the limiting absorption principle, they have shown that the Feynman propagator can be considered as the boundary value of the resolvent of the Klein-Gordon operator.
Assuming Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} and ideas from~\cite[Sec. 3.3]{Gerard_CMP_2015}, 
Lewandowski~\cite[Thm. 4.9]{Lewandowski} 
has recently given a construction of Hadamard states for Riemannian vector bundles. 

Since squares of Dirac-type operators are normally hyperbolic, Theorem~\ref{thm: existence_Feynman_propagator_NHOp} can be used to show that Hadamard states exist for Dirac-type operators when $E$ is equipped with a positive-definite sesquilinear form. 
The existence of a hermitian inner product on which the Dirac operator is formally selfadjoint is however rather exceptional. 
The natural inner product on spinors on spacetimes is not positive-definite but indefinite and therefore the positivity property for the Feynman propagator cannot be concluded directly for the Dirac operator on a globally hyperbolic spacetime with spin structure.
We provide another construction of Feynman propagators with positivity for Dirac-type operators employing a direct microlocal decomposition~\eqref{eq: microlocal_decomposition_causal_propagator_Dirac_type_op} of Pauli-Jordan-Lichnerowicz operator in 
\begin{restatable}{theorem}{HadamardBisolutionDiracTypeOp}
	\label{thm: Hadamard_bisolution_Dirac_type_op}	
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$ endowed with a nondegenerate sesquilinear form $\langle \cdot | \cdot \rangle$ such that the Dirac-type operator $D$ on $E$ is formally selfadjoint. 
	Then there exists a Feynman propagator $\mathsf{S}^{\mathrm{F}}$ for $D$ such that $\omega := - \mathrm{i} (\mathsf{S}^{\mathrm{F}} - S^{\mathrm{adv}})$ is formally selfadjoint. 
	In addition, if there exists a positive-definite sesquilinear form on $E$ induced by Clifford multiplication with any future directed timelike covector (in the sense of~\eqref{eq: positive_definite_sesequilinear_form_Dirac_type_op}), then the Feynman propagator can be chosen such that $\omega$ is nonnegative with respect to $\langle \cdot|\cdot \rangle$ and hence defines Hadamard bisolution of $D$.
\end{restatable}

\section{Microlocalisation}
\subsection{Fourier integral operators}  
\label{sec: FIO}
Fourier integral operators originated from the study of the singularities of solutions of hyperbolic differential equations by 
Lax~\cite{Lax_DukeMathJ_1957} 
and in the context of geometrical optics by 
Maslov~\cite{Maslov_Moscow_1965}. 
Their local formulation has been later systematically developed and globalised in the seminal articles by    
H\"{o}rmander~\cite{Hoermander_ActaMath_1971} 
and by 
Duistermaat-H\"{o}rmander~\cite{Duistermaat_ActaMath_1972} 
for scalar operators and the vector bundle version is available in  
H\"{o}rmander's monograph~\cite{Hoermander_Springer_2009}.  
In this report, we will adopt a symbolic calculus viewpoint, as described below.  

Let $E \to M, F \to N$ be two complex smooth vector bundles over smooth manifolds $M, N$ and let $C$ be a canonical relation from the punctured cotangent bundle $\dot{T}^{*} N$ of $N$ to that $\dot{T}^{*} M$ of $M$. 
In a nutshell, a Fourier integral operator associated with $C$ is a continuous linear map from $C_{\mathrm{c}}^{\infty} (N; F)$ to $\mathcal{D}^{\prime} (M; E)$ whose Schwartz kernel is a Lagrangian distribution~\cite[Def. 25.2.1]{Hoermander_Springer_2009}. 
We will now explain briefly what this means. 
Suppose that $\varLambda$ is a smooth, closed and conic Lagrangian submanifold of the punctured cotangent bundle $\dot{T}^{*} M$.  
The space of Lagrangian distributions $I^m(M,\Lambda)$ can be defined with respect to local coordinate charts. 
To do this we first explain the concept of a clean phase function~\cite[p. 71]{Duistermaat_JInventMath_1975}. 
The reader is referred to the 
monograph~\cite[pp. 414 - 428]{Treves_Plenum_1980} 
for more details.
\begin{definition} \label{def: clean_phase_function}
	Let $M$ be a smooth manifold and $d$ a natural number not necessarily equal to the dimension of $M$. 
	A real-valued smooth function $\varphi$ on an open conic set $\mathcal{U} \subset M \times \dot{\mathbb{R}}^{d}$, homogeneous of degree $1$ in $\theta \in \dot{\mathbb{R}}^{d}$, is called a clean phase function if $\rd \varphi \neq 0$ and its fibre-critical set 
		\begin{equation*}
			\varSigma_{\varphi} := \big\{ (x_{0}, \theta^{0}) \in \mathcal{U} ~|~ \rd_{\theta} \varphi \,  (x_{0}, \theta^{0}) = 0 \big\}
		\end{equation*} 
		is a smooth manifold whose tangent space is given by the kernel of $\rd \partial_{\theta} \varphi$. 
		The number of linearly independent differentials $\rd (\partial \varphi / \partial \theta_{i}), i = 1, \dots, d$ is equal to $d - e$ on $\varSigma_{\varphi}$, where $e = \dim \varSigma_{\varphi} - \dim M$ is called the excess. 
		$\varphi$ is called a nondegenerate phase function if 
		$e = 0$    
		(see e.g.~\cite[Def. 21.2.15]{Hoermander_Springer_2007}).
\end{definition}

A clean phase function on an open set $U \subset \mathbb{R}^n$ defines an immersed conic Lagrangian submanifold $\varLambda_\varphi \subset \dot{T}^*U$ via the map
\begin{equation}
		\varSigma_{\varphi} \ni (x, \theta) \mapsto (x, \rd_{x} \varphi) \in \varLambda_{\varphi}.
\end{equation}
A Lagrangian distribution $u \in I^m(U, \varLambda_\varphi)$ is by definition a distribution given by an oscillatory integral of the form 
\begin{equation} \label{eq: def_IML_local}
		u(x) := (2\pi)^{(n + 2d - 2e) / 4} \int_{\mathbb{R}^{d}} \mathrm{e}^{\mathrm{i} \varphi (x, \theta)} \, a (x, \theta) \, \rd \theta.
\end{equation}
where $\rd \theta$ is the Lebesgue measure on $\mathbb{R}^{d}$, and $a \in S^{m + (n - 2d - 2e) /4} (U \times \dot{\mathbb{R}}^{d})$ is a symbol (Definition~\ref{def: symbol_mf}) with support in the interior of a sufficiently small  conic neighbourhood of
$\varSigma_{\varphi}$ contained in the domain of definition of the phase function. It then follows that $\WF{u} \subset \varLambda_\varphi$~\cite[Thm. 3.2.6]{Hoermander_ActaMath_1971}.

As usual, the oscillatory integral is understood as a formal expression that does not make sense for a fixed $x$ rather it defines a distribution in the sense that for any test function 
$f \in C^\infty_{\mathrm{c}} (U)$ we have
\begin{equation}
 u(f) = (2\pi)^{(n + 2d - 2e) / 4} \int_{\mathbb{R}^{d}} \int_{U} \mathrm{e}^{\mathrm{i} \varphi (x, \theta)} \, a (x, \theta) \, f(x) \, \rd x \, \rd \theta.  
\end{equation}
\begin{definition} \label{def: Lagrangian_distribution_bundle_local}
	Let $E \to M$ be a smooth complex vector bundle over a smooth manifold $M$. 
	Suppose that $\varLambda$ is a smooth, closed and conic Lagrangian submanifold of the punctured cotangent bundle $\dot{T}^{*} M$ of $M$.  
	A distribution $u \in \mathcal{D}' (M; E)$ is called an element of the Lagrangian distribution $I^{m} (M, \varLambda; E)$ of order $m \in \mathbb{R}$ if it can be written as 
	\begin{equation}
		u = \sum_{\alpha} u_{\alpha} 
	\end{equation}
	with locally finite $\supp{u_{\alpha}} \subset U_{\alpha} \subset M$ such that the components $(u^{1}, \ldots, u^{\rk E})$ of $u_{\alpha}$ with respect to a local chart and a bundle chart are in $I^{m} (U_{\alpha}, \varLambda_{\varphi_\alpha})$ with $\varLambda_{\varphi_\alpha} \subset \varLambda$.	
\end{definition}
The above of course implies that~\cite[Lem. 25.1.2]{Hoermander_Springer_2009} 
\begin{equation} \label{eq: WF_Lagrangian_dist}
	\WFPrime{A} \subset \varLambda.
\end{equation} 
The original 
definition~\cite[Sec. 1.4]{Hoermander_ActaMath_1971} 
was based on representations with nondegenerate phase functions. 
Subsequently, it has been 
globalised~\cite[Def. 3.2.2]{Hoermander_ActaMath_1971} 
and extended for sections of vector bundles~\cite[Def. 25.1.1]{Hoermander_Springer_2009}.  
We refer to~\cite[pp. 4 - 10]{Hoermander_Springer_2009} and~\cite[pp. 433 - 439]{Treves_Plenum_1980}  
for a systematic study.  
 
\begin{remark} \label{rem: parametrisation_FIO_canonical_relation}
	The parametrisation of a generic Lagrangian submanifold, $\varLambda = (x, \rd_{x} \varphi)$ by a phase function $\varphi$ cannot be performed globally due to topological restriction given by some cohomology class, known as the Maslov class~\cite[Thm. 1]{Lees_TransAMS_1979}. 
	Furthermore, a phase function does not uniquely characterise a Lagrangian submanifold, rather one has to consider an equivalence class of phase functions satisfying certain necessary and sufficient conditions, as originally analysed for 
	nondegenerate phase functions~\cite[Thm. 3.1.6]{Hoermander_ActaMath_1971} 
	albeit it holds true more generally for nonzero excess. 
	More precisely, if $\varphi$ and $\tilde{\varphi}$ are two clean phase functions in conic neighbourhoods of their fibre-critical points $(x_{0}, \theta^{0}) \in U \times \dot{\mathbb{R}}^{d}$ and $(x_{0}, \tilde{\theta}^{0}) \in U \times \dot{\mathbb{R}}^{\tilde{d}}$ respectively with the same excess $e$, then they are equivalent in some conic open neighbourhoods of these points under a diffeomorphism  from $(x_{0}, \theta^{0})$ to $(x_{0}, \tilde{\theta}^{0})$ if and only if~\cite[Prop. 1.5, p. 421]{Treves_Plenum_1980} 	
	\begin{subequations}
		\begin{eqnarray} 
			(x_{0}, \rd_{x} \varphi) 
			& = & 
			(x_{0}, \rd_{x} \tilde{\varphi}), 
			\\ 
			d 
			& = & 
			\tilde{d}, 
			\\ 
			\mathrm{sgn} \left( \frac{\partial^{2} \varphi}{\partial \theta_{i} \partial \theta_{j}} (x_{0}, \theta^{0}) \right) 
			& = & 
			\mathrm{sgn} \left( \frac{\partial^{2} \tilde{\varphi}}{\partial \tilde{\theta}_{i} \partial \tilde{\theta}_{j}} (x_{0}, \tilde{\theta}^{0}) \right),   
		\end{eqnarray}
	\end{subequations}
	where $\mathrm{sgn}$ denotes the signature of a matrix. 
	We note that a global parametrisation can be achieved, for instance, by an equivalence class of complex-valued nondegenerate phase functions in the particular case when Lagrangian submanifolds are generated by homogeneous symplectomorphisms~\cite[Lem. 1.2, Lem. 1.7]{Laptev_CPAM_1994}. 
\end{remark}
In order to define the class of Fourier integral operators, we substitute the manifold $M$ by a product manifold $M \times N$ and replace the Lagrangian submanifold $\varLambda$ by a homogeneous canonical relation $C$ from $\dot{T}^{*} N$ to $\dot{T}^{*} M$ which is closed in $\dot{T}^{*} (M \times N)$. 
In addition, the role of the vector bundle $E$ will be played by the vector bundle $E \boxtimes F^{*} \to M \times N$. 
We will usually abbreviate  $E \boxtimes F^{*} \to M \times N$ by $\Hom{F, E}$ when it is clear from the context what the base manifold is. 
The class of Fourier integral operators in $I^{m} \big( M \times N, C'; \Hom{F, E} \big)$ can be described in local coordinates as follows. Given a local chart and a local trivialisation we can identify an element
$A \in I^{m} \big( M \times N, C'; \Hom{F, E} \big)$ locally with an $\rk E \times \rk F$-matrix $(A_k^r)$
 with values in $I^{m} (U \times V, C'_{\varphi})$ where $U \subset \mathbb{R}^{\nM}$ and $V \subset  \mathbb{R}^{n_{\ms{N}}}$ are open subsets and $C_\varphi$ is of the form   
\begin{equation} \label{eq: local_twisted_canonical_relation_phase_function}
	C'_{\varphi} = \left\{ \left( x, \rd_{x} \varphi; y, - \rd_{y} \varphi \right) \in \dot{T}^{*} M \times \dot{T}^{*} N \right\}, 
	\quad 
	\rd_{\theta} \varphi \, (x_{0}, y_{0}; \theta^{0}) = 0
\end{equation}
on an open conic neighbourhood of each $(x_{0}, y_{0}; \theta^{0}) \in U \times V \times \dot{\mathbb{R}}^{d}$, generated by a clean phase function $\varphi$ with excess $e$ 
(see e.g.~\cite[(21.2.9)]{Hoermander_Springer_2007}). 
Each element $A$ of $I^{m} \big( M \times N, C'; \Hom{F, E} \big)$ 
has $\WFPrime{A}$ contained in $C'$ and $\WFPrime{A_k^r}$ has wavefront set in $C'_{\varphi}$.
Near $(x_{0}, \xi^{0}; y_{0}, \eta^{0}) \in \WFPrime{A_k^r}$ the singularity of the integral kernel of $A_k^r$ is described by a model of the form
\begin{equation} \label{eq: 25_1_3_clean_Hoemander}
	\tilde A_{k}^{r} (x, y) := (2\pi)^{-(\nM + n_{\ms{N}} + 2d - 2e) / 4} \int_{\mathbb{R}^{d}} \mathrm{e}^{\mathrm{i} \varphi (x, y; \theta)} a_{k}^{r} (x, y; \theta) \, \rd \theta 
\end{equation}
in the sense that in local coordinates $(x_{0}, \xi^{0}; y_{0}, \eta^{0}) \notin \WFPrime{A_k^r - \tilde A_k^r}$.
  
We list some rudimentary properties of Lagrangian distributions below. 
\subsubsection{Principal symbol} 
\label{sec: principal_symbol_FIO}
Loosely speaking, this is the highest order contribution in asymptotic sense of the Lagrangian distribution $I^{m} (M, \varLambda)$. The concept was introduced by 
H\"{o}rmander~\cite[Thm. 3.2.5]{Hoermander_ActaMath_1971} 
for vanishing excess. 
Subsequently its generalisation for nonzero excess was given by
Duistermaat and Guillemin~\cite[Lem. 7.2]{Duistermaat_JInventMath_1975}. 
Later, a geometric description of principal symbol of an arbitrary distribution on a manifold was provided by 
Weinstein~\cite{Weinstein_BullAMS_1976, Weinstein_TransAMS_1978} 
who has also extended H\"{o}rmander's analysis for vector bundles, $I^{m} (M, \varLambda; E)$. 
A compendium of these results is availbale 
in the 
monographs~\cite[pp. 13-16]{Hoermander_Springer_2009} 
and~\cite[pp. 439-450]{Treves_Plenum_1980}. 
We start with the notation 
\begin{equation*}
	I^{m - [1]} (\ldots) 
	:= 
	I^{m} (\ldots) / I^{m - 1} (\ldots).
\end{equation*}  
The principal symbol map will be defined as an
isomorphism ~\cite[Thm. 25.1.9]{Hoermander_Springer_2009} 
\begin{equation} \label{eq: def_symbol_map}
	\sigma : I^{m - [1]} \big( M \times N, C'; \Hom{F, E} \big) \to S^{m + \frac{\nM + n_{\ms{N}}}{4} - [1]} \big( C; \mathbb{M} \otimes \widetilde{\mathrm{Hom}} (F, E) \big), \quad A \mapsto  
	\symb{A}
\end{equation}
where $\widetilde{\mathrm{Hom}} (F, E)$ is the pullback of the bundle $\Hom{F, E} \to N \times M$ to $T^{*} (N \times M)$ followed by restriction to $C$, $\mathbb{M} \to C$ is the 
Keller-Maslov bundle~\cite[p. 155]{Hoermander_ActaMath_1971} 
(see also, e.g.~\cite[Def. 21.6.5]{Hoermander_Springer_2007}), 
and $S^{m + (\nM + n_{\ms{N}}) / 4 -[1]} \big( C; \mathbb{M} \otimes \widetilde{\mathrm{Hom}} (F, E) \big)$ is the $\mathbb{M} \otimes \widetilde{\mathrm{Hom}} (F, E)$-valued symbol on $C$ (see Appendix~\ref{sec: symbol} regarding this symbol class). 
The Keller-Maslov bundle is a smooth complex trivial line bundle whose structure group is $\mathbb{Z} / 4 \mathbb{Z}$.   
For its geometric details we refer the textbook treatments~\cite[pp. 408-412]{Treves_Plenum_1980} 
and~\cite[pp. 686-687]{Rudolph_Springer_2013}.   
The principal symbol is then given explicitly in local coordinates below and the corresponding bundles can be defined by the associated transition functions when changing phase functions and local coordinates. 

Since $A \equiv (A_{k}^{r})$, the principal symbol $\symb{A}$ is locally identified with the matrix $[\symb{A_{k}^{r}}]$. 
Therefore the result in~\cite[Prop. 25.1.5' and p. 15]{Hoermander_Springer_2009}  
for scalar Lagrangian distributions translates to vector bundles as   
\begin{equation}
	\sigma_{A_{k}^{r}} (x, y) = \sqrt{| \rd \xi |} \int_{\varSigma_{\xi^{0}}} a_{k}^{r} (x, y; \theta', \theta'') \, \dfrac{\mathrm{e}^{\mathrm{i} \pi \mathrm{sgn} \varPhi / 4}}{\sqrt{|\det \varPhi|}} \rd \theta'' 
	\in S^{m + \frac{\nM + n_{\ms{N}}}{4} - [1]} \big( C_{\varphi}; \mathbb{M}_{\varphi} \big), 
\end{equation}
where $C_{\varphi}$ is given by~\eqref{eq: local_twisted_canonical_relation_phase_function}, 
$\mathbb{M}_{\varphi}$ is obtained by a locally finite covering of $C$ by open cones with transition functions given by powers of $\mathrm{i}$, and 
\begin{equation}
	\varSigma_{\xi^{0}} 
	:= 
	\left\{ \big( C (y_{0}, \eta^{0}); y_{0}, \eta^{0} \big) \in U \times \mathbb{R}^{d_{\scriptscriptstyle M}} \times V \times \mathbb{R}^{d_{\scriptscriptstyle N}} | 
	\rd_{y} \varphi (y_{0}, \eta^{0}) = 0, \rd_{\eta} \varphi (y_{0}, \eta^{0}) =: \xi^{0} \right\}. 
\end{equation}
Here we split $\mathbb{R}^{d_{\scriptscriptstyle M}} \times \mathbb{R}^{d_{\scriptscriptstyle N}} = \mathbb{R}^{d} \in \theta \equiv (\theta', \theta'') \in \dot{\mathbb{R}}^{d - e} \times \mathbb{R}^{e}$ by a linear change of $\theta$ such that the projection $\varSigma_{\xi^{0}} \ni (x, y; \theta', \theta'') \mapsto \theta'' \in \mathbb{R}^{e}$ has bijective differential so that $T_{(x, y; \theta', \theta'')} \varSigma_{\xi^{0}} || T_{\theta''} \mathbb{R}^{e}$. 
Consequently, $\varphi$ is nondegenerate in a suitable conic open neighbourhood $\mathcal{U}'$ of $(x_{0}, y_{0}; \theta'_{0}, \theta''_{0}) \in U \times V \times \dot{\mathbb{R}}^{d - e}$ for a fixed $\theta''_{0} \in \mathbb{R}^{e}$ and the Hessian $\varPhi$ of $\varphi (x, y; \theta')$ at $(x_{0}, y_{0}; \theta^{0})$ is nonsingular.   
Thereby, we split $\mathcal{U}$ into $\mathcal{U}'$ and its complement (with respect to the Hessian considered as a bilinear form) $\mathcal{U}''$. 
Then $\rd \theta''$ is the induced measure on $\varSigma_{\xi^{0}} |_{\mathcal{U}''}$ in accord with the 
Thom splitting lemma
(see e.g.~\cite[App. C. 6]{Hoermander_Springer_2007} 
or~\cite[Thm. 12.5.2]{Rudolph_Springer_2013}) 
obtained from the pullback of Dirac measure on $\mathbb{R}^{d}$ by $\partial \varphi / \partial \theta_{i}$.  
\subsubsection{Adjoint} 
\label{sec: adjoint_FIO}
To describe the adjoint of a Fourier integral operator we recall some standard notions from linear algebra for our set up. 
Let $E^{*} \to M$ be the dual vector bundle of $E \to M$ and $x$ be any point in $M$. 
Then $\bar{E}_{x}$ is the complex conjugate of the vector space $E_{x}$, defined as the set of conjugate-linear maps from $E_{x}^{*}$ to $\mathbb{C}$ and the identity $E_{x} \to \bar{E}_{x}$ is conjugate-linear. 
Thus, we define a fibrewise sesquilinear form $\langle \cdot | \cdot \rangle_{x}$ on $\bar{E}^{*}_{x} \times E_{x}$ and consequently 
\begin{equation}
	\langle \cdot | \cdot \rangle_{E} := \int_{M} \langle \cdot | \cdot \rangle_{x} : C_{\mathrm{c}}^{\infty} (M; \bar{E}^{*}) \times C_{\mathrm{c}}^{\infty} (M; E) \to \mathbb{C}, 
\end{equation}  
where the first argument in conjugate-linear. 

Any Fourier integral operator $A : C_{\mathrm{c}}^{\infty} (N; F) \to  C^{\infty} (M; E)$ has a formal adjoint with respect to this sesquilinear form~\cite[Sec. 25.2]{Hoermander_Springer_2009}, 
i.e., there exists a Fourier integral operator $A^{*} : C_{\mathrm{c}}^{\infty} (M; \bar{E}^{*}) \to C^{\infty}  (N; \bar{F}^{*})$ with the property 
\begin{equation}
	\langle A^{*} \phi | v \rangle_{F} := \langle \phi | A v \rangle_{E}, 
	\quad \forall \phi \in C_{\mathrm{c}}^{\infty} (M; \bar{E}^{*}), \forall v \in C_{\mathrm{c}}^{\infty} (N; F). 
\end{equation}
In any orthonormal frames of $E \to M$ and $F \to N$ this yields 
\begin{equation}
	\big( A_{k}^{r} \big)^{*} (x, y) = (2\pi)^{-(\nM + n_{\ms{N}} + 2d - 2e) / 4} \int_{\mathbb{R}^{d}} \mathrm{e}^{-\mathrm{i} \varphi (x, y; \theta)} \bar{a}_{r}^{k} (y, x; \theta) \, \rd \theta 
\end{equation}
whenever $A$ is given by~\eqref{eq: 25_1_3_clean_Hoemander}. 
If $A \in I^{m} \big( M \times N, C'; \Hom{F, E} \big)$ then $A^{*} \in I^{m} \big( N \times M, (C^{-1})'$; $\Hom{\bar{E}^{*}, \bar{F}^{*}} \big)$. 
Moreover, if $\symb{A} \in S^{m + (\nM + n_{\ms{N}}) / 4 - [1]} \big( C; \mathbb{M}_{C} \otimes \widetilde{\mathrm{Hom}} (F, E) \big)$ then $\symb{A^{*}} = \mathsf{s} (\symb{A}^{*}) \in S^{m + (\nM + n_{\ms{N}}) / 4 - [1]} \big( C^{-1}; \mathbb{M}_{C^{-1}} \otimes \widetilde{\mathrm{Hom}} (\bar{E}^{*}, \bar{F}^{*}) \big)$ where $\mathsf{s} : N \times M \to M \times N$ is the interchanging map~\cite[Thm. 25.2.2]{Hoermander_Springer_2009}. 
\subsubsection{Algebra of Fourier integral operators} 
\label{sec: alg_FIO}
A first necessary assumption for the product of  two Fourier integral operators to be well-defined is, for the first operator to be properly supported. 
Then the defined composition may however still fail to be a Fourier integral operator.
For example, the composition of two canonical relations does not necessarily have the required properties to define another Fourier integral operator.  
In order to have a well-defined product that is again a Fourier integral operator one needs to assume that the Lagrangian distributions are properly supported and the composition of canonical relations is clean, proper and connected~\cite[Sec. 4.1]{Weinstein_Nice_1975},~\cite[Thm. 5.4]{Duistermaat_JInventMath_1975}.   
If 
$A \in I^{m} \big( M \times \tilde{N}, K'; \Hom{\tilde{F}, E} \big)$ 
and 
$B \in I^{r} \big( \tilde{N} \times N, L'; \Hom{F, \tilde{F}} \big)$ 
with the required restrictions, 
then~\cite[Thm. 25.2.3]{Hoermander_Springer_2009} 
(see also the monograph~\cite[Thm. 5.3, p. 461]{Treves_Plenum_1980})  
\begin{equation} \label{eq: product_Lagrangian_dist_bundle}
	AB := A \circ B \in I^{m + r + e/2} \big( M \times N, (K \circ L)'; \Hom{F, E} \big), 
\end{equation}
where $e$ is the excess of $K \circ L$ and its principal symbol is given by 
\begin{equation} \label{eq: product_symbol_Lagrangian_dist_bundle}
	\symb{AB} = \int_{C_{(x, \xi; z, \eta)}} \symb{A} \times \symb{B}, 
\end{equation} 
where $\sigma_{A} \times \sigma_{B}$ is the $\mathbb{M}$-valued density on $C_{(x, \xi; z, \eta)}$, the compact $e$ dimensional fibre over $(x, \xi; z, \eta) \in K \circ L$ of the intersection of $K \times L$ and $\dot{T}^{*} M \times \varDelta_{\tilde{N}}^{*} \times \dot{T}^{*} N$. 

The space of Fourier integral operators is not an algebra unless the canonical relation is symmetric and transitive~\cite[Ex. 1 and comment on it in the following page]{Guillemin_JFA_1993}. In case the composition of the canonical relation $C$ by itself is clean with excess zero, proper and connected the space of Fourier integral operators $\big( I^\bullet (M \times M, C'; \Hom{E, E}), \circ \big)$ is an associative graded algebra over the field $\mathbb{C}$.
The algebra $\big( I^{\bullet - [1]} \big( M \times N, C'; \Hom{E, E} \big), \eqref{eq: product_symbol_Lagrangian_dist_bundle} \big)$ of principal symbols constitute an algebra which is commutative only the scalar case but non-commutative in general. 
We remark that the product formula~\eqref{eq: product_symbol_Lagrangian_dist_bundle} becomes simpler (cf.~\eqref{eq: product_symbol_FIO_symplecto}) when both $K$ and $L$ are local canonical graphs and then the product is always defined; see Appendix~\ref{sec: FIO_symplecto_bundle} for details. 
\subsection{Connections and the subprincipal symbols of pseudodifferential operators}
The first important observation is that the subprincipal symbol of a pseudodifferential operator with scalar principal part transforms like a (partial-) connection $1$-form under change of bundle charts 
(see e.g.~\cite[Prop. 3.1]{Jakobson_CMP_2007}). 
This is perhaps not surprising given that the subprincipal symbol appears as a constant term in transport equations. 
Connections that are naturally defined from transport equations have appeared first in the work of Dencker on propagation of polarization sets~\cite[p. 365-366]{Dencker_JFA_1982}. 
They are often referred to as Dencker connections in the mathematical physics literature. 
We will now explain the precise relation between geometrically defined connections and the parallel transport induced by the subprincipal symbol in a systematic way. 

For a scalar $P \in \PsiDO{m}{M}$, the subprincipal symbol~\cite[p. 189]{Duistermaat_ActaMath_1972} 
(see also~\cite[Sec. 18.1]{Hoermander_Springer_2007})   
$\subSymb{P}$ is a well defined function on $\dot T^*M$, given by
\begin{equation} \label{eq: subprincipal_symbol}
	 \subSymb{P} = p_{m-1} - \frac{1}{2 \mathrm{i}} \frac{\partial^{2} \symb{P}}{\partial x^{i} \partial \xi_{i}} 
\end{equation}
where $\symb{P}$ is the leading order (principal symbol) and $p_{m-1}$ is the next to the leading order homogeneous term of the total symbol of $P$ (cf. Definition~\ref{def: polyhomogeneous_symbol_mf}) when expressed as an operator acting on half-densities.
If $P \in \PsiDO{m}{M; E}$, then this operator can locally be represented by a $\rk{E} \times \rk{E}$-matrix of pseudodifferential operators if we fix a local bundle frame
$(e_1,\ldots,e_{\rk E})$ of $E$. 
The subprincipal symbol of $P$ is then defined with respect to this local frame as the $\rk{E} \times \rk{E}$-matrix valued function on $\dot T^*M$ given by the subprincipal symbols of the matrix representing $P$ 
(see also~\cite[Def. 3.8]{Hintz_JST_2017}). 
We will see below that under a change of bundle frames this matrix will transform like a connection form along the Hamiltonian vector field in the bundle $\pi^* E \to \dot{T}^{*} M$, where $\pi: \dot T^*M \to M$ is the canonical projection.

First we summarise the transformation properties of the subprincipal symbol. 
In case $P \in \PsiDO{m}{M; E}$ has scalar principal symbol and $Q \in \PsiDO{k}{M; E}$ properly supported, we have the following multiplication formula~\cite[(1.4)]{Duistermaat_JInventMath_1975} 
\begin{equation} \label{eq: subprincipal_symbol_product}
	\subSymb{PQ} = \subSymb{P} \, \symb{Q} + \symb{P} \, \subSymb{Q} + \frac{1}{2 \mathrm{i}} X_{P} (\symb{Q}).  
\end{equation}
Here $X_P = \{ \symb{P},\cdot \}$ is the Hamiltonian vector field generated by the principal symbol $\symb{P}$ of $P$.
In case $P$ has a scalar principal symbol, we have~\cite[(1.3)]{Duistermaat_JInventMath_1975} 
\begin{equation} \label{eq: subprincipal_symbol_power}
	\subSymb{P^{k}} = k \, \symb{P}^{k-1} \, \subSymb{P}, \quad \forall k \in \mathbb{N}. 
\end{equation}
If, moreover $P$ is elliptic and $Q$ is a parametrix of $P$ then we have  
\begin{equation} \label{eq: subprincipal_symbol_inverse}
	\subSymb{Q} = - \, \symb{P}^{-2} \, \subSymb{P}  
\end{equation}
and in this sense the above formula also holds for negative integer $k$.
	
Now let us show that the subprincipal symbol indeed has the claimed transformation property under a change of bundle charts. In fact we show a microlocal version of this statement 
(see also~\cite[Prop. 3.11]{Hintz_JST_2017}). 
\begin{proposition} \label{prop: P_connection_transformation}
	Let $M$ be a manifold and let $m, k \in \mathbb{R}, N \in \mathbb{N}$. 
	Suppose that $P \in \PsiDO{m}{M,\mathrm{Mat} (N)}$ has a scalar principal symbol. 
	Assume that $Q \in \PsiDO{k}{M,\mathrm{Mat} (N)}$ is noncharacteristic at $(x_{0}, \xi^{0})$ and that $\tilde Q \in \PsiDO{-k}{M,\mathrm{Mat} (N)}$ is a microlocal parametrix for $Q$ in an open conic neighbourhood $\mathcal{U}$ of $(x_{0}, \xi^{0})$. 
	Then 
	\begin{equation}
		\subSymb{\tilde Q P Q}(x,\xi) = \symb{Q}^{-1} (x,\xi) \, \subSymb{P} (x,\xi) \, \symb{Q} (x,\xi) + \frac{1}{\mathrm{i}} \symb{Q}^{-1} (x,\xi) \, X_P(\symb{Q}) (x,\xi), \quad \forall (x, \xi) \in \mathcal{U}.
	\end{equation}
\end{proposition}
\begin{proof}
	First note that $(x_{0}, \xi^{0}) \notin \ES{\tilde Q P Q - P - \tilde Q [P,Q]}$. 
	Since $P$ has a scalar principal symbol, the order of $[P,Q]$ is $m+k-1$, and hence the order of $\tilde Q [P,Q]$ is $m-1$.
 	Thus, $\subSymb{\tilde Q P Q} = \subSymb{P}  + \symb{\tilde Q [P,Q]}$ on $\mathcal{U}$. 
 	It remains to compute the principal symbol of $[P,Q]$. Since the principal symbol of $P$ is scalar-valued, we obtain the commutator relation 
 	(see e.g.~\cite[(14)]{Jakobson_CMP_2007}) 
 	\begin{equation} \label{eq: symbol_commutator_PsiDO}
 		\symb{[P,Q]} = \frac{1}{\mathrm{i}} \{\symb{P},\symb{Q}\} + [\subSymb{P},\symb{Q}] 
 	\end{equation}
 	as a consequence of the multiplication formula~\eqref{eq: subprincipal_symbol_product}. 
	Therefore, on $\mathcal{U}$ we get $\symb{Q^{-1} [P,Q]} = - \mathrm{i} \symb{Q}^{-1} \{ \symb{P}, \symb{Q}\} + \symb{Q}^{-1}  \symb{P}  \symb{Q} -  \symb{P}$ and the result follows.
\end{proof}
As a simple consequence of the above product formula we also have the following proposition. 
\begin{proposition} \label{prop: P_orderonerescale}
	Let $M$ be a manifold and let $m, s \in \mathbb{R}, N \in \mathbb{N}$. 
	Suppose that $L \in \PsiDO{s}{M, \mathrm{Mat} (N)}$ has a scalar principal symbol $\symb{L}$ and vanishing subprincipal symbol.
	Assume that $P \in \PsiDO{m}{M,\mathrm{Mat} (N)}$ having a scalar principal symbol $\symb{P}$. 
	Then 
	\begin{equation}
		\subSymb{L P L} = \symb{L}^2 \cdot \subSymb{P}.
	\end{equation}
\end{proposition}
\begin{definition}\label{def: P_compatible_connection}
	Let $(E \to M, \nabla^{E})$ be a smooth complex vector bundle over a manifold $M$ endowed with a connection $\nabla^{E}$ and let $m \in \mathbb{R}$ be arbitrary.
	Suppose that $P \in \PsiDO{m}{M; E}$ whose principal symbol $\symb{P}$ is considered as an element of $C^{\infty} \big( \dot T^{*} M; \End{\pi^{*} E} \big)$ where $\pi : \dot T^{*} M \to M$ is the canonical projection. 
	$\nabla^{E}$ will be called $P$-compatible if and only if  $\subSymb{P} (x_{0}, \xi^{0}) = W (x_{0}, \xi^{0})$ for all $(x_{0}, \xi^{0}) \in \Char{P}$, where $W = \Gamma(\pi_*(X_P))$ is the connection $1$-form $\Gamma$ of $E$ evaluated on $\pi_*(X_P)$, i.e., 
	\begin{equation} \label{eq: def_Dencker_cov_deri}
		\nabla_{X_{P}}^{\End{\pi^{*} E}} = X_{P} + \mathrm{i} \, W
		: C^{\infty} \big( T^{*} M; \End{\pi^{*} E} \big)  \to C^{\infty} \big( T^{*} M; \End{\pi^{*} E} \big).
	\end{equation}
	Here $\nabla_{X_{P}}^{\End{\pi^{*} E}}$ is the covariant derivative along the Hamiltonian vector field $X_{P}$ of $\symb{P}$ on the bundle $\End{\pi^{*} E}$ induced by $\nabla^{E}$.
\end{definition}
This definition makes sense because both quantities, $W$ and $\subSymb{P}$ have the same transformation law under change of bundle frames. Hence $X_P + \mathrm{i} \, \subSymb{P}$ has an invariant meaning.
\begin{proposition} \label{prop: L_P_L_compatible_connection}
	Let $E \to M$ be a smooth complex vector bundle over a manifold $M$ and $m, s \in \mathbb{R}, N \in \mathbb{N}$.  
	Suppose that $P \in \PsiDO{m}{M; E}$ has a scalar principal symbol $\symb{P}$ and that $E$ admits a $P$-compatible connection $\nabla^{E}$. 
	Assume that $L \in \PsiDO{s}{M,\mathrm{Mat} (N)}$ having a scalar principal symbol $\symb{L}$ and vanishing subprincipal symbol. 
	Then $\nabla^{E}$ is $LPL$-compatible.
\end{proposition}
\begin{proof}
	We fix one bundle frame to check this.
	Define $\tilde P = L P L$. Then 
	\begin{equation}
		X_{\tilde P} = \{\symb{\tilde P} , \cdot \} = \{\symb{L}^2\symb{P} , \cdot \} =  \symb{L}^2 \{\symb{P} , \cdot \} + 2 \symb{P}\symb{L}  \{\symb{L} , \cdot \}.
	\end{equation}
	When restricted to covectors in $\Char{P}$ this equals $\symb{L}^2 \{\symb{P} , \cdot \} = \symb{L}^2 X_p$.
	Hence, on $\Char{P}$ we have 
	\begin{equation}
		\nabla_{X_{\tilde P}}^{\End{\pi^{*} E}} - X_{\tilde P} = \symb{L}^2  \Big( \nabla_{X_{p}}^{\End{\pi^{*} E}} - X_{p} \Big)=  \symb{L}^2 \subSymb{P} = \subSymb{L P L} = \subSymb{\tilde P}.
	\end{equation}
\end{proof}
The main observation is now that the Weitzenb\"{o}ck connection defined by a normally hyperbolic operator is compatible with this operator in the above sense. 
\begin{theorem} \label{thm: NHOp_compatible_connection}
	Let $E \to M$ be a smooth complex vector bundle over a spacetime $(M, g)$. 
	Suppose that $\square$ is a normally hyperbolic operator on $E$ and $\nabla^{E}$ is the Weitzenb\"ock connection.
	Then $\nabla^{E}$ is $\square$-compatible.
\end{theorem}
\begin{proof}
	Given a normally hyperbolic operator on a vector bundle over a spacetime, recall that, we can express it by the Weitzenb\"{o}ck formula: $\square = - \tr_{g} \big( \nabla^{T^{*} M \otimes E} \circ \nabla^{E} \big) +V$ where $V$ is a unique potential, $\nabla^{E}$ is the Weitzenb\"{o}ck connection and 
	\begin{equation}
		\nabla^{T^{*} M\otimes E}  =\nabla^{\mathrm{LC}} \otimes 1 + 1 \otimes \nabla^{E} 		
	\end{equation}
	is the connection induced by the Levi-Civita connection $\nabla^{\mathrm{LC}}$ and the Weitzenb\"{o}ck-connection. 
	Since the potential $V$ is of order zero, it does not contribute to the subprincipal symbol. 
	We will check compatibility in a local frame $(e_1,\ldots,e_{\rk E})$. 
	Then, we have
	$\nabla_{X}^{E} u = X^{\mu} (\partial_{\mu} u^{i} + \mathrm{i} \, \Gamma_{\mu j}^{i} u^{j}) e_{i}$ where $C^{\infty} (M; TM) \ni X = X^{\mu} \partial_{\mu}, \mu = 1, \ldots, n$ and $\Gamma$ is the connection $1$-form of Weitzenb\"{o}ck connection: $\nabla_{\partial_{\mu}}^{E} e_{j} = \mathrm{i} \, \Gamma_{\mu j}^{i} e_{i}$, and the Weitzenb\"{o}ck formula reads 
	\begin{equation}
		- \tr_{g} \big( \nabla^{T^{*} M \otimes E} \circ \nabla^{E} \big) + V  
		=
		- g^{\mu \nu} \nabla_{\mu}^{E} \nabla_{\nu}^{E}  + g^{\mu \nu} \Gamma_{\quad \mu \nu}^{\mathrm{LC} \, \rho} \partial_{\rho} +V_1
	\end{equation}
	for some potential $V_1$. 
	This expression is the formula for the operator acting on functions. 
	The formula~\eqref{eq: subprincipal_symbol} for the subprincipal symbol needs to be applied to the full symbol of the operator acting on half-densities. 
	Since we have the canonical Lorentzian volume form $\sqrt{|\det g|} \rd x^{1} \wedge \ldots \wedge \rd x^{n}$ on $(M, g)$, the formula for a normally hyperbolic operator on half-densities in local coordinates is given by
	\begin{equation}
		|\det g|^{\frac{1}{4}} \; \square \; |\det g|^{-\frac{1}{4}}
		= 
		- g^{\mu \nu} \nabla_{\mu}^{E} \nabla_{\nu}^{E} - \frac{\partial g^{\mu \nu}}{\partial x^\mu} \partial_\nu + V_3  
	\end{equation}
	for some potential $V_{2}$. 
	Using this representation we see that the subprincipal symbol of $\square$ is 
	\begin{equation} \label{eq: subprincipal_symbol_NHOp}
		\subSymb{\square} (x, \xi) \, e_j (x) = - 2 g^{\mu \nu} (x) \, \Gamma_{\nu j}^i (x) \, \xi_\mu \, e_i (x). 
	\end{equation}
	The Hamiltonian vector field of the principal symbol of $\square$ in local coordinates is  given by 
	\begin{equation} \label{eq: HVF_NHOp}
		X_{\square} = \frac{\partial g^{\mu \nu}}{\partial x^\alpha} \xi_\mu \xi_\nu \frac{\partial}{\partial \xi_\alpha} - 2 g^{\mu \nu} \xi_\mu \frac{\partial}{\partial x^\nu}.
	\end{equation}
	Hence, the pushforward $\pi_*(X_{\square})$ is $-2 g^{\mu \nu} \xi_\mu \partial_{x^{\nu}}$ and we can see that
	\begin{equation}
		\Gamma \big( \pi_*(X_{\square}) \big) e_j  = - 2 g^{\mu \nu} \Gamma_{\nu j}^i \xi_\mu e_i = \subSymb{\square} \, e_j. 
	\end{equation}
\end{proof}
As a simple application of this theorem and Proposition~\ref{prop: L_P_L_compatible_connection} we have
\begin{corollary} \label{cor: L_NHOp_L_compatible_connection}
	Let $E \to M$ be a smooth complex vector bundle over a spacetime $(M, g)$. 
	Suppose that $L \in \PsiDO{-1/2}{M,\mathrm{Mat} (N)}$ having a scalar principal symbol and vanishing subprincipal symbol. 
	Then $\nabla^{E}$ is $L \square L$-compatible, i.e. the Weitzenb\"ock connection is compatible with the first order operator $L \square L$.
\end{corollary}
\subsection{Microlocal conjugate of a normally hyperbolic operator} 
\label{sec: microlocalisation_NHOp}
Two pseudodifferential operators are called microlocally conjugate if they can be conjugated to one another by an elliptic Fourier integral operator once they have been appropriately localised in cotangent space. 
The key point is that under some natural assumptions a first order pseudodifferential operator can be microlocally conjugated to a vector field. 
This is due to 
Duistermat and H\"{o}rmander~\cite[Prop. 6.1.4]{Duistermaat_ActaMath_1972} 
for scalar operators which has been extended to vector bundles by 
Dencker~\cite{Dencker_JFA_1982} 
who formulated microlocal conjugate of a system of classical pseudodifferential operators locally of real principal type.  
More precisely, Dencker transformed the system of operators to a scalar pseudodifferential operator with vanishing subprincipal symbol $P$, modulo smoothing operators by conjugating with a system of elliptic pseudodifferential operators.  
Then he depicted the microlocal conjugation of $P$ with $D_{1} := - \mathrm{i} \mathbbm{1}_{\mathrm{Mat} (N)} \partial_{1}$ by an elliptic Fourier integral operator associated with the graph of the  symplectomorphism microlocally connecting $\symb{P}$ to $\symb{D_{1}}$.  

In this section, we will explain microlocalisation in an intrinsic geometric language and present Dencker's result in a slightly more general form. 
Let us now consider complex-matrix valued smooth functions $C^{\infty} \big( \Rn, \mathbb{C}^{N} \big)$ on Euclidean space $\Rn$ and let $D_{1} : C^{\infty} \big( \Rn, \mathbb{C}^{N} \big) \to C^{\infty} \big( \Rn, \mathbb{C}^{N} \big)$. 
We will show now that any first order pseudodifferential operator $P$ on $E$ having a scalar principal symbol is microlocally conjugate to $D_{1}$, in the sense of~\eqref{eq: microlocal_conjugate_P} in the proposition below. 
\begin{center}
	\begin{tikzpicture}
		\node (a) at (0,0) {$\mathcal{U}$};
		\node (b) at (4,0) {$\mathcal{U}'$}; 
		\node (c) at (0,2) {$\End{\pi^{*} E}_{\mathcal{U}}$};
		\node (d) at (4,2) {$\mathrm{Mat} (N)$}; 
		\node[below] at (2, 0) {$\varkappa$}; 
		\node[left] at (0,1) {$\symb{P}$}; 
		\node[above] at (2.1, 2) {$\alpha$}; 
		\node[right] at (4,1) {$\symb{D_{1}}$};  
	
		\draw[->] (b) -- (a); 
		\draw[->] (d) -- (c); 
		\draw[->] (a) -- (c); 
		\draw[->] (b) -- (d);
	\end{tikzpicture}
	\hfil
	\begin{tikzpicture}
		\node (a) at (0,0) {$C_{\textrm{c}}^{\infty} (M; E)$};
		\node (b) at (5,0) {$C_{\textrm{c}}^{\infty} (\Rn, \mathbb{C}^{N})$}; 
		\node (c) at (0,2) {$C_{\textrm{c}}^{\infty}(M; E)$};
		\node (d) at (5,2) {$C_{\textrm{c}}^{\infty} (\Rn, \mathbb{C}^{N})$}; 
		\node[below] at (2.2, 0) {$B$}; 
		\node[left] at (0,1) {$P$}; 
		\node[above] at (2.2, 2) {$\tilde{B}$}; 
		\node[right] at (5,1) {$D_{1}$}; 
	
		\draw[->] (b) -- (a); 
		\draw[->] (a) -- (c); 
		\draw[->] (c) -- (d); 
		\draw[->] (b) -- (d);
	\end{tikzpicture}
	\captionof{figure}{A schematic diagram of microlocalisation. The diagram on the right commutes in a microlocal sense near the point $(x_{0}, \xi^{0})$ and the map $\alpha$ is defined by  $\alpha(T) = \symb{B} T \symb{\tilde{B}}$.}
	\label{fig: microlocalisation}
\end{center}
\begin{proposition}[Microlocalisation] 
\label{prop: microlocalisation_P}
    Let $E \to M$ be a smooth complex vector bundle over an  $n$-dimensional manifold $M$.	
    Suppose that $P \in \PsiDO{1}{M; E}$ is properly supported with real scalar principal symbol $p$ such that: 
	\begin{enumerate}[label=(\alph*)]
		\item \label{con: p_x_xi_zero}
		$p (x_{0}, \xi^{0}) = 0$ for some $(x_{0}, \xi^{0}) \in \dot{T}^{*} M$ and  
		\item \label{con: HVF_radial}
		the Hamiltonian vector field $X_{P}$ of $p$ and the radial direction are linearly independent at $(x_{0}, \xi^{0})$.
        \end{enumerate}
	Then for any $m \in \mathbb{R}$ there exist\\
	{(i)}  a homogeneous symplectomorphism $\varkappa$ from an open conic neighbourhood $\mathcal{U}'$ of $(0, \eta_{1} \rd y^{1})$ in $\dot{T}^{*} \! \mathbb{R}^{n}$ to an open conic coordinate chart $\big( \mathcal{U}, (x^{1}, \ldots, x^{n}$; $\xi_{1}, \ldots, \xi_{n}) \big)$ of $(x_{0}, \xi^{0})$ in $\dot{T}^{*} M$ with graph  $\varGamma$ such that 
	\begin{equation} \label{eq: choice_symplecto}
		\varkappa^{*} p = \xi_{1} \mathbbm{1}_{\End{E}}, 
	\end{equation} and\\ \noindent
	{(ii)} properly supported Lagrangian distributions 
	$B \in I^{m} \big( M \times \Rn, \varGamma'; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$   
	and 
	$\tilde{B} \in I^{-m} \big( \Rn \times M, (\varGamma^{-1})'; \Hom{E, \Rn \times \mathbb{C}^{N}} \big)$ 
	such that $B \tilde{B}, \tilde{B} B$ both are zero order pseudodifferential operators and  
	\begin{subequations}
		\begin{eqnarray}
			&& 
			\WFPrime{B} \subset \mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_{1})}, 
			\qquad 
			\WFPrime{\tilde{B}} \subset \mathcal{U}'_{(0, \eta_{1} \rd y^{1}; x_{0}, \xi^{0})},
			\label{eq: microlocalisation_1}
			\\ 
			&& 
			(x_{0}, \xi^{0}) \notin \ES{B \tilde{B} - I_{E}}, 
			\qquad 
			(0, \eta_{1} \rd y^{1}) \notin \ES{\tilde{B} B - I}, 
			\label{eq: microlocalisation_2}
			\\ 
			&& 
			(x_{0}, \xi^{0}) \notin \ES{B D_{1} \tilde{B} - P}, 
			\quad 
			(0, \eta_{1} \rd y^{1}) \notin \WFPrime{\tilde{B} P B - D_{1}}.      
			\label{eq: microlocal_conjugate_P}
		\end{eqnarray}
	\end{subequations}
Here $D_{1} := - \mathrm{i} \mathbbm{1}_{\mathrm{Mat} (N)} \partial / \partial y^{1} : C^{\infty} (\Rn, \mathbb{C}^{N}) \to C^{\infty} (\Rn, \mathbb{C}^{N})$ and $\mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_{1}\rd y^{1})}, \mathcal{U}'_{(0, \eta_{1}\rd y^{1}; x_{0}, \xi^{0})}$ are small conic neighbourhoods of $(x_{0}, \xi^{0}; 0, \eta_{1} \rd y^{1}) \in \dot{T}^{*} M \times \dot{T}^{*} \! \mathbb{R}^{n}, (0, \eta_{1}\rd y^{1}; x_{0}, \xi^{0}) \in \dot{T}^{*} \! \mathbb{R}^{n} \times \dot{T}^{*} M$, respectively. 

	In addition, if $E$ is endowed with a nondegenerate sesquilinear form with respect to which $P$ is formally selfadjoint then $\tilde{B}$ can be chosen as the adjoint of $B$ if $\mathbb{C}^N$ is endowed with a standard sesquilinear form of the same signature as the sesquilinear form on $E$. 
	
	Moreover, if $E$ is equipped with a $P$-compatible connection $\nabla^{E}$ then the principal symbols of $B$ and $\tilde{B}$ can be chosen $\mathbbm{1}$ near $(x_{0}, \xi^{0}; 0, \eta_{1} \rd y^{1})$ and $(0, \eta_{1} \rd y^{1}; x_{0}, \xi^{0})$ with respect to a frame that is parallel along the Hamiltonian vector field $X_{P}$. 
	
	Furthermore, if $E$ is equipped with a nondegenerate sesquilinear form such that $P$ is formally selfadjoint and a $P$-compatible connection $\nabla^{E}$ then we can choose $B$ such that the principal symbol of $B$ equals $\mathbbm{1}$ near $(x_{0}, \xi^{0}; 0, \eta_{1} \rd y^{1})$ with respect to a frame that is unitary and parallel along the Hamiltonian vector field $X_{P}$ and  $\tilde{B}=B^*$.
\end{proposition}
A schematic of this notion has been portrayed in Figure~\ref{fig: microlocalisation}.  
\begin{proof}
	We will prove the proposition imitating the strategy used for the scalar version~\cite[Prop. 6.1.4 and Lem. 6.6.4]{Duistermaat_ActaMath_1972} 
	(see also the~\cite[Prop. 26.1.3]{Hoermander_Springer_2009}). 
	The existence of $\varkappa$ satisfying~\eqref{eq: choice_symplecto} is guaranteed by a 
	homogeneous Darboux theorem~\cite[Lem. 6.6.3]{Duistermaat_ActaMath_1972} 
	(see also~\cite[Thm. 21.3.1]{Hoermander_Springer_2007}) 
	which prerequisites our hypotheses~\ref{con: p_x_xi_zero} and~\ref{con: HVF_radial}.  
	Suppose that $b \in S^{m} \big( \varGamma; \mathbb{M} \otimes \widetilde{\mathrm{Hom}} (\Rn \times \mathbb{C}^{N}, E) \big)$ has an inverse in a conic neighbourhood of $(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1}) \in \varGamma$. 
	Then we can obtain a properly supported  $B_{1} \in I^{m} \big( M \times \Rn, \varGamma'; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$ such that $\WFPrime{B_{1}} \subset \mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1})}$ and $B_{1}$ is noncharacteristic (cf. Definition~\ref{def: elliptic_FIO}) at $(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1})$ from the construction given in 
	Appendix~\ref{sec: FIO_symplecto_bundle}. 



	  
	By Theorem~\ref{thm: existence_parametrix_FIO}, there exists a unique microlocal parametrix $\tilde{B}_{1} \in I^{- m} \big( \Rn \times M, (\varGamma^{-1})'$; $\Hom{E, \Rn \times \mathbb{C}^{N}} \big)$ so that~\eqref{eq: microlocalisation_2} is satisfied . 
	Since $B_{1}$ and $\tilde{B}_{1}$ have reciprocal principal symbols to each other on $\mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1})}$, so $\tilde{B}_{1} P B_{1}$ has the principal symbol $\eta_{1} \mathbbm{1}_{\mathrm{Mat} (N)}$ on $\mathcal{U}'$ due to~\eqref{eq: choice_symplecto} and the Egorov Theorem~\ref{thm: Egorov}. 
	Furthermore 
	\begin{equation} \label{eq: WF_B1_P_A1_D_Q}
		(0, \eta_1 \rd y^{1}) \notin \ES{\tilde{B}_{1} P B_{1} - D_{1} - Q}, 
		\quad 
		Q \in \PsiDON{0}.  
	\end{equation}



	
	To find the operators $B$ and $\tilde{B}$ with the claimed properties we construct properly supported elliptic $B_{2}, \tilde{B}_{2} \in \PsiDON{0}$ such that
	\begin{eqnarray}
		&& 
		\tilde{B}_{2} B_{2} - I \in \PsiDON{-\infty}, 
		\label{eq: 26_1_03_Hoermander}
		\\ 
		&& 
		\tilde{B}_{2} (D_{1} + Q) B_{2} - D_{1} \in \PsiDON{-\infty},    
		\label{eq: 26_1_3_Hoermander}
	\end{eqnarray}
	and we set $B := B_{1} B_{2}$ and $\tilde{B} := \tilde{B}_{2} \tilde{B}_{1}$. 
	Then~\eqref{eq: 26_1_03_Hoermander} and~\eqref{eq: 26_1_3_Hoermander} together with~\eqref{eq: WF_B1_P_A1_D_Q} imply 
	\begin{eqnarray}
		(0, \eta_1 \rd y^{1}) \notin \ES{\tilde{B}_{2} (\tilde{B}_{1} B_{1} - I) B_{2}} 
		= 
		\ES{\tilde{B} B - I},  
		\\ 
		(0, \eta_1 \rd y^{1}) \notin \ES{\tilde{B}_{2} (\tilde{B}_{1} P B_{1} - D_{1} - Q) B_{2}} 
		= 
		\ES{\tilde{B} P B - D_{1}}, 
	\end{eqnarray}
	which proves the second half of~\eqref{eq: microlocalisation_2} and~\eqref{eq: microlocal_conjugate_P} and the first half follows immediately after we multiply from left and right by $B$ and by $\tilde{B}$.



	
	It therefore remains to construct $B_2$ and $\tilde{B}_2$ such that~\eqref{eq: 26_1_03_Hoermander} and~\eqref{eq: 26_1_3_Hoermander} hold.
	By the existence of a parametrix 
	(see e.g.~\cite[Thm. 18.1.24]{Hoermander_Springer_2007}), 
	for every elliptic operator $B_{2} \in \PsiDON{0}$ there exists $\tilde{B}_{2} \in \PsiDON{0}$ such that $\tilde{B}_{2} B_{2} - I$ and $B_{2} \tilde{B}_{2} - I$ are smooth. 
	Multiplying~\eqref{eq: 26_1_3_Hoermander} by $B_{2}$ from the left we arrive at the equivalent condition for~\eqref{eq: 26_1_3_Hoermander} that 
	\begin{equation} \label{eq: 26_1_3_Hoermander_equiv} 
		[D_{1}, B_{2}] + Q B_{2} \in \PsiDON{-\infty}  
		\tag{\ref{eq: 26_1_3_Hoermander}$^{\prime}$} 
	\end{equation}
	for some elliptic $B_{2}$. 
	We will now construct a solution of~\eqref{eq: 26_1_3_Hoermander_equiv} order by order, starting with the principal symbol. 
	The principal symbol of~\eqref{eq: 26_1_3_Hoermander_equiv} vanishes provided  
	\begin{equation}
		- \mathrm{i} \left\{ \symb{D_{1}}, \symb{B_{2}} \right\} + \symb{Q} \symb{B_{2}} = 0,   
	\end{equation} 
	as the the subprincipal symbol of $D_{1}$ vanishes; cf.~\eqref{eq: symbol_commutator_PsiDO}.  
	If $q$ is the homogeneous principal symbols of $Q$, then the preceding equation yields for the homogeneous principal symbol $b_0$ of $B_2$ the equation
	\begin{equation} \label{eq: microlocalisation_first_transport}
		\parDeri{y^{1}}{b_0} = - \mathrm{i} q b_0.   
	\end{equation}
	This is a first order differential equation, hence a unique solution exists given the initial condition $(b_0)|_{y^{1}=0} = \mathbbm{1}_{\mathbb{C}^{N}}$ and this solution depends smoothly on $q$.
	%

	By construction  $b_0 (y, \eta)$ is homogeneous of degree zero and therefore $(b_0)^{-1}$ exists and is homogeneous of degree zero.  
	Moreover, $\det (b_0 (y, \eta))$ is nonvanishing.  
	Defining a properly supported $B_{2, 0} \in \PsiDON{0}$ with homogeneous principal symbol $b_0$, we now successively construct $B_{2, k} \in \PsiDON{-k}$ 
	so that, for every $k \in \mathbb{N}$:  
	\begin{equation} \label{eq: microlocalisation_k_transport} 
		[D_{1}, B_{2, 0} + \ldots + B_{2, k}] 
		+ 
		Q (B_{2, 0} + \ldots + B_{2, k}) = R_{k+1} \in \PsiDON{- (k+1)}.   
	\end{equation} 
	This is equivalent to the corresponding homogeneous principal symbols $b_{k}$ of $B_{2, k}$ and  $r_{k}$ of $R_{k}$, respectively, of degree $-k$ to satisfy	
	\begin{equation}
		- \mathrm{i} \parDeri{y^{1}}{b_{k}} + q b_{k} = - r_{k}. 
	\end{equation} 	
	This equation can be solved by the Duhamel principle and the solution reads 
	\begin{equation} \label{eq: microlocalisation_sol_k_transport}
		b_{k} (y, \eta) = - \mathrm{i}  \, b_0 (y, \eta) \int_{0}^{y^{1}} 
		b_0^{-1} (t, y^{2}, \ldots, y^{n}; \eta) \, r_{k} (t, y^{2}, \ldots, y^{n}; \eta) \, \rd t.   
	\end{equation}
	Then, using asymptotic summation of the symbols of $B_{2,k}$ (cf. Appendix~\ref{def: polyhomogeneous_symbol_mf}) we can now construct an operator $B_2$ satisfying \eqref{eq: 26_1_3_Hoermander_equiv}. 	\newline 



	
	\noindent 
	\textbf{The case when $P$ is formally selfadjoint:} 
	\\ 
	Next we prove that the operators $B$ and $\tilde{B}$ can be chosen microlocally unitary in case $P$ is formally selfadjoint with respect to a nondegenerate sesquilinear form on $E$.
	We endow the space $C_{\mathrm{c}}^{\infty} (\Rn, \mathbb{C}^{N})$ with a standard sesquilinear scalar product of the same signature as the sesquilinear form on $E$. The operator $D_{1}$ is then formally selfadjoint.  
	Futhermore, we make the choice $m = 0$. 
	Acting with $B \in I^{0} \big( M \times \Rn, \varGamma'; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$ from the left of~\eqref{eq: microlocal_conjugate_P} gives the equivalent microlocal conjugate relation between $P$ and $D_{1}$:  
	\begin{equation}
		(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1}) \notin \WFPrime{PB - B D_{1}}. 
	\end{equation}
	Taking its adjoint we have $(0, \eta_1 \rd y^{1}; x_{0}, \xi^{0}) \notin \WFPrime{B^{*} P - D_{1} B^{*}}$ where  $B^{*} \in I^{0} \big( \Rn \times M$, $(\varGamma^{-1})'; \Hom{E^{*}, \Rn \times \mathbb{C}^{N}} \big)$ and consequently   
	\begin{subequations}
		\begin{eqnarray}
			&& (0, \eta_1 \rd y^{1}) \notin \ES{B^{*} P B - D_{1} B^{*} B}, 
			\\ 
			&& (0, \eta_1 \rd y^{1}) \notin \ES{B^{*} P B - B^{*} B D_{1}}.
		\end{eqnarray}  
	\end{subequations}
	Thus $(0, \eta_1 \rd y^{1}) \notin \ES{[B^{*} B, D_{1}]}$ and $B$ is noncharacteristic at $(0, \eta_1 \rd y^{1})$. 
	Now the operator $B^{*} B$ is a pseudodifferential operator and since $\mathrm{i} D_{1}$ is differentiation with respect to $y^1$ the total symbol of the commutator $[B^{*} B, D_{1}]$ is the $-\mathrm{i}$ times the $y^1$-derivative of the total symbol of $B^{*} B$.
	Thus, on the the level of symbols, the above implies that the total symbol of $B^{*} B$ is the sum of a term independent of $y^1$ and a term that is rapidly decaying in a conic neighborhood $\mathcal{U}'$ of $\eta_1 \rd y^{1}$. 
	It is bounded below on $\mathcal{U}'$ because of the ellipticity of $B^{*} B$ there. 
	By Proposition~\ref{prop: FIO_2_2_2_matrix}, one can find a selfadjoint properly supported $\varPsi \in \PsiDONU{0}$ whose principal symbol is the same as that of $B$ and such that  
	\begin{subequations}
		\begin{eqnarray}
			&& 
			(0, \eta_1 \rd y^{1}) \notin \Char{\varPsi},
			\label{eq: microlocalisation_varPsi_1}
			\\ 
			&& 
			(0, \eta_1 \rd y^{1}) \notin \ES{\varPsi^{*} \varPsi - B^{*} B}, 
			\label{eq: microlocalisation_varPsi_2}
			\\ 
			&& 
			(0, \eta_1 \rd y^{1}) \notin \ES{[\varPsi, D_{1}]}, 
			\label{eq: microlocalisation_varPsi_3} 
		\end{eqnarray}
	\end{subequations}
	where $U'$ is the projection of $\mathcal{U}'$ on $\Rn$. 
	Note that Proposition~\ref{prop: FIO_2_2_2_matrix} (stated for operators acting on the same bundle) can be used only in a local trivialisation of $E$ by an orthonormal frame that identifies the fibre of $E$ with $\mathbb{C}^N$ in such a way that the sesquilinear forms are identified. It is at this stage that we must require the sesquilinear forms to have the same signature.
	
	Here the last property follows from the fact that the construction of the full symbol of $\varPsi$ in Proposition~\ref{prop: FIO_2_2_2_matrix} involves only multiplication and asymptotic summation of symbols. 
	The property of a symbol being a sum of two terms, one independent of $y_1$ and one rapidly decaying in a conic neighborhood of $\eta_1 \rd y^{1}$, is preserved under these operations. 
	The full symbol of $\varPsi$ is therefore also of this form.



	
	The ellipticity~\eqref{eq: microlocalisation_varPsi_1} entails a unique microlocal parametrix $\varPhi$ for $\varPsi$. 
	In other words, there exists a properly supported $\varPhi \in \PsiDONU{0}$ such that 
	\begin{subequations}
		\begin{eqnarray}
			&&  
			(0, \eta_1 \rd y^{1}) \notin \Char{\varPhi},
			\label{eq: microlocalisation_varPhi_1} 
			\\ 
			&& 
			(0, \eta_1 \rd y^{1}) \notin \ES{\varPhi \varPsi - I} 
			\Leftrightarrow (0, \eta_1 \rd y^{1}) \notin \ES{\varPsi \varPhi - I}.
			\label{eq: microlocalisation_varPhi_2} 
		\end{eqnarray}
	\end{subequations}
	Note,~\eqref{eq: microlocalisation_varPhi_1},~\eqref{eq: microlocalisation_varPsi_2} and~\eqref{eq: microlocalisation_varPsi_3} imply that 
	\begin{subequations}
		\begin{eqnarray}
			&& (0, \eta_1 \rd y^{1}) \notin \ES{\varPhi^{*} \varPsi^{*} \varPsi \varPhi - \varPhi^{*} B^{*} B \varPhi}, 
			\\ 
			&& (0, \eta_1 \rd y^{1}) \notin \ES{\varPhi [\varPsi, D_{1}] \varPhi}.   
		\end{eqnarray}
	\end{subequations}
	This, accounting~\eqref{eq: microlocalisation_varPhi_2} entails that 
	\begin{subequations}
		\begin{eqnarray}
			&& (0, \eta_1 \rd y^{1}) \notin \ES{I - (B \varPhi)^{*} B \varPhi}, 
			\\ 
			&& (0, \eta_1 \rd y^{1}) \notin \ES{[\varPhi, D_{1}]},  
		\end{eqnarray}
	\end{subequations} 
	which completes the proof since $(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1}) \notin \WFPrime{B [D_{1}, \varPhi]}$ and therefore with $B^{\backprime} := B \varPhi$ we have $(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1}) \notin \WFPrime{P B^{\backprime} - B^{\backprime} D_{1}}$. 
	\newline 


	%
	
	\noindent 
	\textbf{The case of connection $P$-compatibility:} 
	\\
	Now, we assume that $\nabla^{E}$ is $P$-compatible. 
	Since the construction of the symbols of $B$ and $\tilde{B}$ are local, we can fix a local frame and local coordinates. 
	We will reduce the general case to the case when the subprincipal symbol of $P$ vanishes near $(x_{0}, \xi^{0})$. 
	\\
	\textit{(i) The case of vanishing subprincipal symbol}:  \\
	We assume that $P$ is given in local coordinates with respect to a some bundle frame and that in this frame $\subSymb{P} = 0$ near $(x_{0}, \xi^{0})$. 
	Let $\tilde P$ be a scalar operator with the same principal symbol as $P$ and with vanishing subprincipal symbol.  
	Using Weyl-quantisation, the operators $B_1$ and $\tilde{B}_1$ can be constructed in such a way that Egorov's theorem holds up to the subprincipal symbol level~\cite[Thm. 1]{Silva_PublMat_2007}. We can choose these operators in such a way that their principal symbols are constant $\mathbbm{1}$ and such that the subprincipal symbol of $\tilde{B}_1 \tilde P B_1$ vanishes near $(0, \eta_1 \rd y^{1})$, and therefore so does the subprincipal symbol of  $\tilde{B}_1 P B_1$ near $(x_{0}, \xi^{0})$. We will now use the same construction as before starting with $B_1$ and $\tilde B_1$. Then the principal symbol $q$ of the remainder term $Q$ vanishes. The construction of $B_2$ and $\tilde{B}_2$ then yields operators with total symbols that are scalar and constant principal symbols equal to $\mathbbm{1}$. We conclude that in case the subprincipal symbol of $P$ vanishes near $(x_{0}, \xi^{0})$ the operators $B$ and $\tilde{B}$ can be chosen as scalar operators. The principal symbols that are constant along the flow lines of $X_P$ and $X_{D_1}$. 



	\\
	\textit{(ii) The case of nonvanishing subprincipal symbol}: \\
	We will microlocally transform $P$ to a scalar pseudodifferential operator $\tilde{P}$. 
	To be specific, for any properly supported $\tilde{P} \in \PsiDO{1}{M}$, we want to have a  $\hat{B} \in \PsiDO{0}{M; E}$ such that $\hat{B}$ is noncharacteristic at $(x_{0}, \xi^{0})$ and $(x_{0}, \xi^{0}) \notin \ES{P \hat{B} - \hat{B} \tilde{P} I}$. We construct this operator locally and therefore fix a local frame for the bundle $E$ near the point $x_0$. In the pullbacked bundle $\pi^* E \to \dot T^* M$ we can also construct a local frame that is parallel along $X_P$ with respect to $\nabla_{X_{P}}^{\pi^{*} E}$. This local parallel frame is constructed by choosing a local transverse to $X_P$ and use the original frame on this transverse. Parallel transport along the flow lines of $X_P$ then gives the desired frame. The change of frame matrix from the original frame to the parallel frame is then a local section $b$ of $\End{\pi^{*} E} \to \dot{T}^{*} M$. 
	By construction, this frame is homogeneous of degree zero. 
	We now choose an elliptic zero order pseudodifferential operator $\hat B$ whose principal symbol $\symb{\hat B}$ equals $b$ on $\mathcal{U}$. 
	Let $\check{B}$ be a parametrix for $\hat B$. 
	The subprincipal symbol of $\check{B} P \hat B$ equals
	$
	- \mathrm{i} X_{P} \, \symb{\hat{B}} + [\subSymb{P}, \symb{\hat{B}}] 
	$ 
	on $\mathcal{U}$, using~\eqref{eq: symbol_commutator_PsiDO}. By $P$-compatibility this is exactly the formula for the connection one form in the parallel local bundle frame and it therefore vanishes.
	By Proposition~\ref{prop: P_connection_transformation}, this is precisely the formula for the subprincipal symbol of $\check{B} P \hat{B}$ on $\mathcal{U}$ and hence, $\tilde P = \check{B} P \hat{B}$ (modulo smoothing operators) has vanishing subprincipal symbol on $\mathcal{U}$. 



	
	This reduces the problem to the case of vanishing subprincipal symbol discussed earlier and let $B_3$ and $\tilde{B}_3$ are the corresponding scalar Fourier integral operators. 
	This means $B$ and $\tilde{B}$ are of the form $B := \hat{B} B_3$ and $\tilde{B} := \tilde{B}_3 \check{B}$. 
	Since the principal symbol of $\hat{B}$ is the transition function to a parallel frame and $B_3$ is a scalar operator, these imply that the principal symbols of $B$ and $\tilde{B}$ can be chosen $\mathbbm{1}$ with respect to a parallel frame along $X_P$.  
	\newline  



	
	\noindent 
	\textbf{The case of selfadjoint $P$ with connection $P$-compatibility:} 
	\\
	Finally, suppose that $E$ is equipped with a sesquilinear form with respect to which $P$ is formally selfadjoint and $\nabla^{E}$ is $P$-compatible. 
	The above construction of $\hat B$ can be repeated using an orthonormal frame of $E$ and another orthonormal frame of $\pi^* E$ that is parallel along $X_P$. One obtains an operator that has principal symbol $\mathbbm{1}$ near $(x_0,\xi^0)$ with respect to the parallel orthonormal frame. Now one repeats the construction of $\varPsi$ as before and note that the principal symbol can be kept as $\mathbbm{1}$ in that way to turn $B$ into a microlocally unitary operator in the sense that $\tilde B = B^*$ with the desired properties.
\end{proof}
We will now show microlocalisation a normally hyperbolic operator $\square$ as a consequence of the preceding result by replacing the generic manifold $M$ by a Lorentzian (not necessarily globally hyperbolic) manifold $(M, g)$. 
In particular, the role of the characteristic point $(x_{0}, \xi^{0})$ of $P$ in Proposition~\ref{prop: microlocalisation_P}~\ref{con: p_x_xi_zero} will be played by a lightlike covector on $(M, g)$ and thus, the integral curves of the Hamiltonian vector field $X_{P}$ will be substituted by lifts of null-geodesics to the cotangent bundle (cf.~\eqref{eq: HVF_NHOp}). 
\begin{proposition} 
\label{prop: microlocalisation_NHOp}
	Let $E \to M$ be a smooth complex vector bundle over an $n$-dimensional Lorentzian manifold $(M, g)$ and let $\Box$ be a normally hyperbolic operator on $E$. Denote by $\nabla^{E}$,   the Weitzenb\"{o}ck connection associated with $\Box$ and by $\varGamma$, the graph of the homogeneous symplectomorphism from an open conic neighbourhood $\mathcal{U}'$ of $(0, \eta_1 \rd y^{1})$ in $\dot{T}^{*} \! \mathbb{R}^{n}$ to an open conic coordinate chart $\big( \mathcal{U}, (x, \xi) \big)$ centered at any lightlike covector $(x_{0}, \xi^{0})$ on $(M, g)$.
	Then for any $m \in \mathbb{R}$, one can find properly supported Lagrangian distributions 
	$A \in I^{-1/2 + m} \big( M \times \Rn, \varGamma'; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$  
	and 
	$\tilde{A} \in I^{-1/2 - m} \big( \Rn \times M, (\varGamma^{-1})'; \Hom{E, \Rn \times \mathbb{C}^{N}} \big)$ such that $A \tilde{A}, \tilde{A} A$ both are pseudodifferential operators of order $-1$ and 
	\begin{subequations}
		\begin{eqnarray}
			&& 
			\WFPrime{A} \subset \mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1})}, 
			\qquad 
			\WFPrime{\tilde{A}} \subset \mathcal{U}'_{(0, \eta_1 \rd y^{1}; x_{0}, \xi^{0})},
			\label{eq: microlocalisation_NHOp_1}
			\\ 
			&& 
			(x_{0}, \xi^{0}) \notin \ES{A \tilde{A} - I_{E}}, 
			\qquad 
			(0, \eta_1 \rd y^{1}) \notin \ES{\tilde{A} A - I}, 
			\label{eq: microlocalisation_NHOp_2}
			\\ 
			&& 
			(x_{0}, \xi^{0}) \notin \ES{A D_{1} \tilde{A} - \Box}, 
			\qquad 
			(0, \eta_1 \rd y^{1}) \notin \ES{\tilde{A} \Box A - D_{1}}.      
			\label{eq: microlocal_conjugate_NHOp}
		\end{eqnarray}
	\end{subequations}
	Here $D_{1} = - \mathrm{i} \mathbbm{1}_{\mathrm{Mat} (N)} \partial / \partial y^{1} : C^{\infty} \big( \Rn, \mathbb{C}^{N} \big) \to C^{\infty} \big( \Rn, \mathbb{C}^{N} \big)$ and $\mathcal{U}_{(x_{0}, \xi^{0}; 0, \eta_1 \rd y^{1})}, \mathcal{U}'_{(0, \eta_1 \rd y^{1}; x_{0}, \xi^{0})}$ are small conic neighbourhoods of $(x_{0}, \xi^{0}; 0, \eta_{0}) \in \dot{T}^{*} M \times \dot{T}^{*} \! \mathbb{R}^{n}, (0, \eta_1 \rd y^{1}; x_{0}, \xi^{0}) \in \dot{T}^{*} \! \mathbb{R}^{n} \times \dot{T}^{*} M$, respectively. 



	
	In addition, if $E$ is endowed with a nondegenerate sesquilinear form with respect to which $\square$ is formally selfadjoint then $A$ can be chosen as a scalar operator with respect to a unitary bundle frame that is parallel along the geodesic flow and microlocally $A = \tilde{A}^*$ if $\mathbb{C}^N$ is endowed with a standard sesquilinear form of the same signature as the sesquilinear form on $E$. . 
\end{proposition}
\begin{proof}
	The strategy is, as usual, to reduce $\square$ to a first order operator so that  Proposition~\ref{prop: microlocalisation_P} can be applied. 
	To do so choose a properly supported elliptic self-adjoint pseudodifferential operator on $E$ of order $-1/2$ having a homogeneous scalar principal symbol $l$ and vanishing sub-principal symbol. Such an operator always exists, as it can be constructed locally and then patched to a global operator using a suitable partition of unity.  Then $L \square L$ is a first order pseudodifferential operator whose principal symbol  $l^{2} (x, \xi) \, g_{x}^{-1} (\xi, \xi) \, \mathbbm{1}_{\End{E}}$ vanishes on lightlike covectors $(x_{0}, \xi^{0})$. 
	By Theorem~\ref{thm: NHOp_compatible_connection} and Corollary~\ref{cor: L_NHOp_L_compatible_connection}, $\nabla^{E}$ is compatible with both $\square$ and $L\square L$, respectively. 
	Therefore the hypotheses of Proposition~\ref{prop: microlocalisation_P} are satisfied and the conclusion entails by putting $A := LB, \tilde{A} := \tilde{B} L$ where $B, \tilde{B}$ are the Lagrangian distributions in Proposition~\ref{prop: microlocalisation_P}. In case $\square$ is formally self-adjoint with respect to a sesquilinear form on $E$ we choose $L$ self-adjoint.  
\end{proof}
\begin{remark}
\label{rem: global_microlocalisation_NHOp} 
	We have microlocally conjugated $\square$ to $D_{1}$ around a lightlike covector in the preceding proposition~\eqref{eq: microlocal_conjugate_NHOp}. 
	It is also possible to refine the construction as follows.
	Let $\mathbb{I}$ be a compact interval in $\mathbb{R}$ and $c : \mathbb{I} \to \dot{T}^{*} M$ an integral curve of $X_{P}$:    
	\begin{equation}
		p \circ c = 0. 
	\end{equation} 
	If the composition of $c$ and the projection $\dot{T}^{*} M \to \dot{T}^{*} M / \mathbb{R}_{+}$ is injective, then one can find a conic neighbourhood $\mathcal{V}'$ of $\mathbb{I} \times \{(0, \eta_{1} \rd y^{1}) \}$ and a smooth homogeneous symplectomorphism $\kappa$ from $\mathcal{V}'$ to an open conic neighbourhood $\kappa (\mathcal{V}') \subset \dot{T}^{*} M$ of $c (\mathbb{I})$ such that~\cite[Prop. 26.1.6]{Hoermander_Springer_2009}  
	\begin{equation}
		\kappa \big( \mathbb{I} \times \{ (0, \eta_{1} \rd y^{1}) \} \big) = c (\mathbb{I}), \qquad \kappa^{*} p = \xi_{1} \mathbbm{1}_{\End{E}}. 
	\end{equation}
	Imitating the proof of Proposition~\ref{prop: microlocalisation_P}, one can microlocalise $P$ to $D_{1}$ on $c (\mathbb{I})$ 
	(see~\cite[Prop. 26.1.3$'$]{Hoermander_Springer_2009} for the scalar version).  
	Consequently, if $\Gamma$ is the graph of a homogeneous symplectomorphism from a conic neighbourhood of $\mathbb{I} \times \{ (0, \eta_{1} \rd y^{1}) \}$ in $\dot{T}^{*} \! \mathbb{R}^{n}$ to a conic neighbourhood of lightlike geodesic (strip) $\gamma (\mathbb{I})$ on $\dot{T}^{*} M$ then the proof of Proposition~\ref{prop: microlocalisation_NHOp} shows that for any $m \in \mathbb{R}$, there exists Lagrangian distributions 
	$A \in I^{-1/2 + m} \big( M \times \Rn, \Gamma'; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$  
	and 
	$\tilde{A} \in I^{-1/2 - m} \big( \Rn \times M, (\Gamma^{-1})'; \Hom{E, \Rn \times \mathbb{C}^{N}} \big)$ such that $A \tilde{A}, \tilde{A} A$ both are pseudodifferential operators of order $-1$ and   
	\begin{subequations}
		\begin{eqnarray}
			&& 
			\WFPrime{A} \subset \mathcal{V}_{\Gamma}, 
			\qquad 
			\WFPrime{\tilde{A}} \subset \mathcal{V}'_{\Gamma^{-1}},
			\label{eq: microlocalisation_NHOp_1_global}
			\\ 
			&& 
			\gamma (\mathbb{I}) \cap \ES{A \tilde{A} - I_{E}} = \emptyset, 
			\qquad 
			\mathbb{I} \times \{ (0, \eta_{1} \rd y^{1}) \} \cap \ES{\tilde{A} A - I} = \emptyset, 
			\label{eq: microlocalisation_NHOp_2_global}
			\\ 
			&& 
			\gamma (\mathbb{I}) \cap \ES{A D_{1} \tilde{A} - \Box} = \emptyset, 
			\qquad 
			\mathbb{I} \times \{ (0, \eta_{1} \rd y^{1}) \} \cap \ES{\tilde{A} \Box A - D_{1}} = \emptyset,      
			\label{eq: microlocal_conjugate_NHOp_global}
		\end{eqnarray}
	\end{subequations}
	where $\mathcal{V}_{\Gamma}, \mathcal{V}'_{\Gamma^{-1}}$ are small conic neighbourhoods of $\Gamma$ restricted to $c (\mathbb{I})$ and its inverse, respectively. 
\end{remark}
We close this section by a simple application of Proposition~\ref{prop: microlocalisation_P} to derive a bundle version of 
H\"{o}rmander's propagation of singularity theorem~\cite{Hoermander_Nice_1970},~\cite[Thm. 6.1.1']{Duistermaat_ActaMath_1972}, 
as for example stated in  
Taylor~\cite[Thm. 4.1, p. 135]{Taylor_PUP_1981} (for Sobolev wavefront sets)
and
Dencker~\cite[Thm. 4.2]{Dencker_JFA_1982} (for the polarisation set).  
Since $C^{\infty} (M; E) = \bigcap_{s \in \mathbb{R}} H_{\mathrm{loc}}^{s} (M; E)$, such a refinement of the usual notion of (smooth) wavefront set is captured by the Sobolev wavefront set. 
For any $s \in \mathbb{R}$, the Sobolev wavefront set $\mathrm{WF}^{s} (u)$ of a distribution $u \in \mathcal{D}' (M; E)$ relative to Sobolev space $H_{\mathrm{loc}}^{s} (M; E)$ is defined by~\cite[p. 201]{Duistermaat_ActaMath_1972}  
\begin{equation} \label{eq: def_Sobolev_WF}
	\mathrm{WF}^{s} (u) 
	:= 
	\bigcap_{\substack{\Psi \in \PsiDO{0}{M; E} \\ \Psi u \in H_{\mathrm{loc}}^{s} (M; E)}} \Char{\Psi} 
	= 
	\bigcap_{\substack{\Psi \in \PsiDO{s}{M; E} \\ \Psi u \in L_{\mathrm{loc}}^{2} (M; E)}} \Char{\Psi} 
\end{equation}
where $\Char{\Psi}$ is the characteristic set of $\Psi$ and the intersection is over all properly supported $\Psi$. 
Locally this means, for any open subset $U$ of $\Rn$, $(x_{0}, \xi^{0}) \notin \mathrm{WF}^{s} (u)$ is not in the Sobolev wavefront set of a distribution $u \in \mathcal{D}' (U; E)$ if and only if there exists a compactly supported section $f$ on $U$ nonvanishing at $x_{0} \in U$ such that the weighted Fourier transform $(1 + \| \xi \|^{2})^{s} \, |\widehat{fu}(\xi)|^{2}$ is integrable in a conic neighbourhood of $ \xi^{0}$. 
An accumulation of important properties of this finer wavefront set is available, for example, in~\cite[App. B]{Junker_AHP_2002}. 
A propagation of $H^{s}$-regularity is now presented below. 
The result is similar to that of Taylor but formulated in a geometric fashion and proved in a different way employing microlocalisation on vector bundles developed in Proposition~\ref{prop: microlocalisation_P}.    
\begin{theorem}[Propagation of Sobolev regularity] 
\label{thm: propagation_singularity_Sobolev_WF}
	Let $E \to M$ be a smooth complex vector bundle over a manifold $M$ and $u \in \mathcal{D}' (M; E)$. 
	Suppose that, for some $m \in \mathbb{R}$, $\Phi \in \PsiDO{m}{M; E}$ satisfies the hypotheses~\ref{con: p_x_xi_zero} and~\ref{con: HVF_radial} of Proposition~\ref{prop: microlocalisation_P} and that $\mathcal{I}$ be an interval on an integral curve (on cotangent bundle) of the Hamiltonian vector field of the principal symbol of $\Phi$ such that $\mathcal{I} \cap \mathrm{WF}^{s} (\Phi u) = \emptyset$. 
	Then either $\mathcal{I} \cap \mathrm{WF}^{s + m -1} (u) = \emptyset$ or $\mathcal{I} \subset \mathrm{WF}^{s + m -1} (u)$. 
\end{theorem}
\begin{proof}
    Once Proposition~\ref{prop: microlocalisation_P} is at our disposal, the rest of the proof is the same as its scalar version~\cite[Thm. 6.1.1']{Duistermaat_ActaMath_1972}  
	(see also~\cite[Thm. 26.1.4]{Hoermander_Springer_2009}). 	For completeness we give the details here.



	 
	By conjugating $\Phi$ with appropriate Fourier integral operators we reduce the statement to the case $m = 1$. 
	In other words, one can just consider $P$ in Proposition~\ref{prop: microlocalisation_P} instead of $\Phi$ to conclude the assertion by utilising the Sobolev continuity properties of pseudodifferential 
	(see e.g.~\cite[p. 92]{Hoermander_Springer_2007}) 
	and Fourier integral operators 
	(see e.g.~\cite[Cor. 25.3.2]{Hoermander_Springer_2009}). 
	By Proposition~\ref{prop: microlocalisation_P}, the analysis further boils down to $P = D_{1}, s = 0$ and $(x_{0}, \xi^{0}) = (0, \eta_1 \rd y^{1})$ by choosing the order of the Fourier integral operator $B$ in Proposition~\ref{prop: microlocalisation_P} equals to $-s$. 
	Then the claim follows from the form of the advanced and retarded fundamental solutions of $D_{1}$ and as these map from $L_{\mathrm{c}}^{2} (\Rn; \Rn \times \mathbb{C}^{N})$ to $L_{\mathrm{loc}}^{2} (\Rn; \Rn \times \mathbb{C}^{N})$. 
\end{proof}
\subsection{Proof of Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp}}
We will make use of the strategy used by 
Duistermaat-H\"{o}rmander~\cite[Thm. 6.5.3]{Duistermaat_ActaMath_1972} 
(see also~\cite[Thm. 26.1.14]{Hoermander_Springer_2009}) 
for a very general class of scalar pseudodifferential operators and exploit the global hyperbolicity of manifold. 
To begin with, we show that the Feynman parametrix is unique, if it exists. 
\subsubsection{Uniqueness} 
\label{sec: uniqueness_Feynman_parametrix}
Suppose that $G_{\mathrm{L}}$ is a left and $G_{\mathrm{R}}$ is a right Feynman parametrix of $\Box$, i.e.,  
$\WFPrime{G_{\mathrm{L}}, G_{\mathrm{R}}} \subset \varDelta^{*} \cup C^{+}$ where $C^{+}$ is defined in Definition~\ref{def: Feynman_parametrix}. 
This implies that $G_{\mathrm{L}}, G_{\mathrm{R}} : C_{\mathrm{c}}^{\infty} (M; E) \to C^{\infty} (M; E)$. 
To prove that $G_{\mathrm{L}} - G_{\mathrm{R}}$ is a smoothing operator, we would like to argue that $G_{\mathrm{L}} \Box G_{\mathrm{R}}$ is congruent both to $G_{\mathrm{L}}$ and to $G_{\mathrm{R}}$ modulo smoothing operators. 
But $G_{\mathrm{L}}$ and $G_{\mathrm{R}}$ are not properly supported which makes this a nontrivial task. 
To circumvent this difficulty, we employ the fact that $G_{\mathrm{L}} H  G_{\mathrm{R}}$ is defined when $H$ is a pseudodifferential operator having Schwartz kernel of compact support in $M \times M$ and then $H : \mathcal{D}^{\prime} (M; E) \to \mathcal{E}^{\prime} (M; E)$. 

If $(x, \xi; y, \eta) \in \WFPrime{G_{\mathrm{L}} H G_{\mathrm{R}}}$ but $(x, \xi), (y, \eta)$ are in the complement of $\ES{H}$, then $(x, \xi; z, \zeta), (z, \zeta; y, \eta) \in C^{+}$ for some $(z, \zeta) \in \ES{H}$. This follows from from the behaviour of wavefront sets under composition of kernels,
Theorem 8.2.10 in~\cite{Hoermander_Springer_2003},
and the fact that $\WFPrime{H} \subset \varDelta^{*}$. 
This entails that $(x, \xi), (y, \eta)$ and $(z, \zeta)$ are on the same lightlike geodesic (strip) $\big( \gamma (s), \rd \gamma \, (s) \big)$ on $\dot{T}^{*} M$ with $(z, \zeta)$ in between the other two points. 

Since $(M, g)$ is assumed to be globally hyperbolic, $K' := J^{+} (K) \cap J^{-} (K)$ is compact for any compact $K \subset M$~\cite[Cor. Prop. 6.6.1]{Hawking_CUP_1973}. 
Therefore, if the endpoints of $\gamma$ are in $K$ then $\gamma (s) \in K'$, as shown in Figure~\ref{fig: uniqueness_Feynman_parametrix}.  
Consequently, if $\singsupp{H} \cap K' = \emptyset$ then $\WFPrime{G_{\mathrm{L}} H G_{\mathrm{R}}}$ cannot have any point over $K \times K$. 
Let $\chi \in C_{\mathrm{c}}^{\infty} (M; E)$ such that $\chi$ is identically $1$ on $K'$ and define $\mathrm{M}_{\chi}$ be the corresponding multiplication operator. 
We observe that $[\mathrm{M}_{\chi}, \square]$ vanishes identically on $K'$. 
Thus $\WFPrime{G_{\mathrm{L}} [\mathrm{M}_{\chi}, \square] G_{\mathrm{R}}}$ contains no  point over $K \times K$ and so it true for 
\begin{eqnarray}
	(\dot{T}^{*} K \times \dot{T}^{*} K) \cap \WFPrime{G_{\mathrm{L}} - G_{\mathrm{R}}}
	& = & 
	(\dot{T}^{*} K \times \dot{T}^{*} K) \cap \WFPrime{G_{\mathrm{L}} \mathrm{M}_{\chi} \square G_{\mathrm{R}} - G_{\mathrm{L}} \square \mathrm{M}_{\chi} G_{\mathrm{R}}} 
	\nonumber \\ 
	& = &    
	(\dot{T}^{*} K \times \dot{T}^{*} K) \cap \WFPrime{G_{\mathrm{L}} [\mathrm{M}_{\chi}, \square] G_{\mathrm{R}}} 
	\nonumber \\ 
	& = & 
	\emptyset. 
\end{eqnarray}
Since $K$ is arbitrary, we conclude that $G_{\mathrm{L}} - G_{\mathrm{R}}$ is a smoothing operator and similarly for the anti-Feynman parametrix.  
\begin{center}
	\begin{tikzpicture}
		\filldraw[color = cyan, fill = cyan!5!white] (0,0) -- (2,2) -- (4,0) -- (2,-2) -- (0,0); 
		\draw[blue, thick] (0,0) -- (4,0); 
		\draw[cyan, thick] (0,0) .. controls (2,1.5) .. (4,0); 
		\draw[cyan, thick] (0,0) .. controls (2,-1.5) .. (4,0);
	
		\node[left] at (0,0) {$\textcolor{blue}{K}$}; 
		\node[left] at (1.85,1.75) {$\textcolor{cyan}{K'}$}; 
		\node at (2,0.75) {\textcolor{cyan}{$\gamma (s)$}};
	\end{tikzpicture}
	\hfil
	\begin{tikzpicture}
		\draw[cyan, ultra thick] (1,0) -- (4,0); 
		\draw[red, ultra thick] (0,0) -- (1,0); 
		\draw[red, ultra thick] (4,0) -- (5,0); 
		\draw[thick] (1,2) -- (4,2);
		\draw[thick] (1,2) to [out=190, in=20] (0,0); 
		\draw[thick] (4,2) to [out=355, in=160] (5,0); 
		\draw[->] (0,0) -- (0, 2.5);
	
		\node[above] at (2.5, 2) {$\chi$};
		\node[below] at (2.5,0) {$\textcolor{cyan}{\gamma (s)}$}; 
		\node[left] at (0,2) {$1$};
	\end{tikzpicture}
	\captionof{figure}{A schematic visualisation of a consequence of global hyperbolicity of $(M, g)$ to entail the uniqueness of Feynman parametrices. 
	Here time flows from left to right.}
    \label{fig: uniqueness_Feynman_parametrix}	
\end{center}
\subsubsection{Construction} 
Notice that $\square G_{\mathrm{R}} = I + R$ is equivalent to $(G_{\mathrm{R}})^{*} \square^{*} = I + R^{*}$ and, since the adjoint $\square^{*}$ of $\square$ is also a normally hyperbolic operator, the existence of left parametrices with the listed properties in this theorem follows from the existence of right parametrices for $\square$. 
It is therefore sufficient to construct a right parametrix for $\square$ with the required regularities. 

To begin with, we note that the Schwartz kernels of the advanced, retarded and Pauli-Jordan fundamental solutions 
\begin{eqnarray} 
	&& 
	F^{\mathrm{ret}, \mathrm{adv}} = \pm \mathrm{i} \, \Theta (\pm y^{1} \mp z^{1}) \otimes \delta (y^{2} - z^{2}, \ldots, y^{n} - z^{n}), 
	\label{eq: adv_ret_fundamental_sol_partial_derivative} 
	\\ 
	&&  
	F := F^{\mathrm{ret}} - F^{\mathrm{adv}} = 
	\mathrm{i} (2 \pi)^{n-1} \int_{\mathbb{R}^{n-1}} \exp \Big( \mathrm{i} (y^{2} - z^{2}) \theta_{2} + \ldots + \mathrm{i} (y^{n} - z^{n}) \theta_{n} \Big) \, \rd \theta 
	\label{eq: causal_fundamental_sol_partial_derivative} 
\end{eqnarray}
of $D_{1}$ satisfy~\cite[Prop. 6.1.2]{Duistermaat_ActaMath_1972}  
(see also~\cite[Prop. 26.1.2]{Hoermander_Springer_2009})
\begin{eqnarray} 
	&& 
	\WFPrime{F^{\mathrm{ret}, \mathrm{adv}}} 
	= 
	\varDelta_{\mathbb{R}^{2n}}^{*} 
	\cup 
	\big\{ (x^{1}, x^{i}; 0, \xi_{i}; y^{1}, x^{i}; 0, \xi_{i}) | x^{1}, y^{1}, x^{i} \in \mathbb{R}, \xi_{i} \in \dot{\mathbb{R}}, x^{1} \gtrless y^{1} \big\}, \quad 
	\label{eq: Hoermander_Prop_26_1_2_i}
	\\
	&& 
	F, \vartheta F^{\mathrm{ret}, \mathrm{adv}} \in 
	I^{- \frac{1}{2}} \Big( \Rn \times \Rn, \big\{ (y^{1}, y^{i}; 0, \eta_{i}; z^{1}, y^{i}; 0, -\eta_{i}) ~|~ y^{1}, z^{1}, y^{i} \in \mathbb{R}, \eta_{i} \in \dot{\mathbb{R}} \big\} \Big) \quad 
	\label{eq: Hoermander_Prop_26_1_2_ii}
\end{eqnarray}
for $i = 2, \ldots, n$ and if $\vartheta$ is a smooth function on $\Rn \times \Rn$ vanishing near the diagonal. Again the integral \eqref{eq: causal_fundamental_sol_partial_derivative} is an oscillatory integral that needs to be understood in the sense of distributions.

We are now going to present a bundle version of the following result by 
Duistermaat and H\"{o}rmander~\cite[Lem. 6.5.4]{Duistermaat_ActaMath_1972}  
(see also~\cite[Lem. 26.1.15]{Hoermander_Springer_2009}). 
The proof simply carries over to the setting of bundles yet we include the details here for completeness. 
\begin{lemma} \label{lem: Lem_26_1_15_Hoermander}
	We use the set up and terminologies as in Proposition~\ref{prop: microlocalisation_NHOp}, that is, the normally hyperbolic operator $\square$ and $D_{1}$ are microlocally conjugate to each other, $A, \tilde{A}$ are the conjugating Lagrangian distributions with $m = 0$, and $\mathcal{U}$ is the conic neighbourhood of any lightlike covector $(x_{0}, \xi^{0})$ on a globally hyperbolic manifold $(M, g)$. 
	Suppose that $F^{\mathrm{ret}, \mathrm{adv}}$ are the Schwartz kernels of the retarded and the advanced fundamental solutions of $D_{1}$. 
	We set $F^{\pm} := \chi F^{\mathrm{ret}, \mathrm{adv}}$ where $\chi \in C^{\infty} (\Rn \times \Rn)$ is identically $1$ in a neighbourhood of the diagonal and vanishes outside another sufficiently small neighbourhood of the diagonal.  
	If $T \in \PsiDO{0}{M; E}$ with $\ES{T} \subset \mathcal{U}$ and $\tilde{F}^{\pm} := A F^{\pm} \tilde{A} T$, then  
	\begin{subequations}
		\begin{eqnarray}
			&& 
			\WFPrime{\tilde{F}^{\pm}} \subset \varDelta^{*} \cup C^{\pm}, 
			\label{eq: Hoermander_Lem_26_1_15_i}
			\\ 
			&& 
			\square \tilde{F}^{\pm} = T + R^{\pm}, 
			\quad 
			R^{\pm} \in I^{-3/2} \big( M \times M, (C^{\pm})', \End{E} \big), 
			\label{eq: Hoermander_Lem_26_1_15_ii}
			\\ 
			&& 
			\tilde{F}^{+} - \tilde{F}^{-} \in I^{-3/2} \big( M \times M, C', \End{E} \big).  
			\label{eq: Hoermander_Lem_26_1_15_iii} 
		\end{eqnarray}
	\end{subequations}
\end{lemma}
\begin{proof}
	The properties~\eqref{eq: Hoermander_Prop_26_1_2_i} and~\eqref{eq: Hoermander_Prop_26_1_2_ii} immediately entail~\eqref{eq: Hoermander_Lem_26_1_15_i} and~\eqref{eq: Hoermander_Lem_26_1_15_iii}. 
	To show~\eqref{eq: Hoermander_Lem_26_1_15_ii}, by definition:  
	\begin{equation} \label{eq: Hoermander_26_1_11}
		\square \tilde{F}^{\pm} 
		= \square A F^{\pm} \tilde{A} T 
		= (\square A - A D_{1}) F^{\pm} \tilde{A} T + A D_{1} F^{\pm} \tilde{A} T. 
	\end{equation} 
	Since $\square$ and $D_{1}$ are microlocal conjugates to one another,~\eqref{eq: microlocal_conjugate_NHOp} can be re-expressed as  
	\begin{equation} 
		(x_{0}, \xi^{0}; 0, \eta_{1} \rd y^{1}) \notin \WFPrime{\square A - A D_{1}}. 
		\tag{\ref{eq: microlocal_conjugate_NHOp}$^{\prime}$}   
	\end{equation}
	Then there is a conical neighbourhood $\mathcal{V}$ of $(0, \eta_{1} \rd y^{1})$ such that $(\square A - A D_{1}) v$ is smooth for any $v \in \mathcal{D}' (\Rn, \mathbb{C}^{N})$ when $\WF{v} \subset \mathcal{V}$. 
	Since $\WFPrime{F^{\pm}}$ can be made arbitrary close to the diagonal in $\dot{T}^{*} \! \mathbb{R}^{n} \times \dot{T}^{*} \! \mathbb{R}^{n}$ by choosing the support of $\chi$ close to the diagonal in $\Rn \times \Rn$, we can pick $\chi$ and a conic neighbourhood $\mathcal{W}$ of $(0, \eta_{1} \rd y^{1})$ such that $\WF{F^{\pm} v} \subset \mathcal{V}$ provided $\WF{v} \subset \mathcal{W}$. 
	If $\ES{T} \subset \varkappa (\mathcal{W})$ then the first term in the right-hand side of~\eqref{eq: Hoermander_26_1_11} is smooth and we are left with the last term of this equation. 
	By definition: $D_{1} F^{\pm} = I + \tilde{R}^{\pm}$ where $\tilde{R}^{\pm} := D_{1} (\chi) F^{\mathrm{adv}, \mathrm{ret}}$ and 
	\begin{equation}
		\tilde{R}^{\pm} 
		\in 
		I^{- \frac{1}{2}} \Big( \Rn \times \Rn, \big\{ (x^{1}, x^{i}; 0, \xi_{i}; y^{1}, x^{i}; 0, \xi_{i}) ~|~ x^{1}, y^{1}, x^{i} \in \mathbb{R}, \xi_{i} \in \dot{\mathbb{R}}, x^{1} \gtrless y^{1} \big\} \Big). 
	\end{equation} 
	Since $A \tilde{A} T - T = (A \tilde{A} - I) T$ is smooth as long as $\ES{T}$ is sufficiently close to $(x_{0}, \xi^{0})$, it follows that $\square \tilde{F}^{\pm} - T = R^{\pm}$ where $R^{\pm} - A \tilde{R}^{\pm} \tilde{A} T$ is smooth, which concludes the proof. 
\end{proof} 
To end the proof of Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp}, we choose a locally finite covering $\{ \mathscr{U}_{\alpha} \}$ of $\dot{T}^{*} M$ by open cones $\mathscr{U}_{\alpha}$ such that Lemma~\ref{lem: Lem_26_1_15_Hoermander} is applicable when $\ES{T} \subset \mathscr{U}_{\alpha}$. 
Denoting the projection of $\mathscr{U}_{\alpha}$ on $M$ by $U_{\alpha}$ and picking  $\mathscr{U}_{\alpha}$ so that $U_{\alpha}$ are also locally finite, we set 
\begin{equation}
	I = \sum_{\alpha} T_{\alpha}, 
	\quad 
	\ES{T_{\alpha}} \subset \mathscr{U}_{\alpha}, 
	\quad 
	\supp{T_{\alpha}} \subset U_{\alpha} \times U_{\alpha}. 
\end{equation}
Then for every $\alpha$, $A F^{\pm} \tilde{A} T_{\alpha}$ can be chosen according to the Lemma~\ref{lem: Lem_26_1_15_Hoermander} with $\mathrm{supp} (A F^{\pm} \tilde{A} T_{\alpha}) \subset U_{\alpha} \times U_{\alpha}$, and hence the sum 
\begin{equation} \label{eq: Feynman_parametrices_NHOp}
	G^{\mathrm{F, aF}} := \sum_{\alpha} A F^{\pm} \tilde{A} T_{\alpha} 
\end{equation}
is well-defined which satisfies the claimed properties in the statement of Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp}.  
Now note that $F^{\pm} := \chi F^{\mathrm{ret}, \, \mathrm{adv}}$ provided $\chi$ is taken as a function of $y-z$ and then, $F^{\pm}$ is a convolution by a measure of compact support and therefore continuously maps $H^s_{\mathrm{c}}(\mathbb{R}^n,\mathbb{C}^N)$ to $H^s_{\mathrm{c}}(\mathbb{R}^n,\mathbb{C}^N)$.
By the mapping properties of pseudodifferential and Fourier integral operators, ~\cite[Thm. 4.3.1]{Hoermander_ActaMath_1971}   
(see also~\cite[Cor. 25.3.2]{Hoermander_Springer_2009}), the factors $A, \tilde{A}, T_{\alpha}$ are 
continuous maps from $H^s_{\mathrm{c}}$ to $H^s_{\mathrm{c}}$. Since the sum of kernels is locally finite this shows that
 $G^{\mathrm{F, aF}}$ continuously maps $H_{\mathrm{c}}^{s} (M; E)$ into $H_{\mathrm{loc}}^{s} (M; E)$ for all $s \in \mathbb{R}$.  

Until now, we have just shown that $\square G^{\mathrm{F, aF}} = I + R^{\pm}$. 
By Lemma~\ref{lem: Hoermander_lem_26_1_16}, we can choose $G^{\pm} \in I^{-3/2} \big( M \times M, (C^{\pm})', \End{E} \big)$ so that $P G^{\pm} - R^{\pm} \in C^{\infty} (M \times M)$. 
The $G^{\pm}$ extend to continuous mappings from $H_{\mathrm{c}}^{s} (M; E)$ to $H_{\mathrm{loc}}^{s} (M; E)$ for every $s \in \mathbb{R}$. This follows from mapping properties of Fourier integral operators, for example, Theorem 25.3.8 in~\cite{Hoermander_Springer_2009} 
bearing in mind that the corank of the $2$-form $\omega_{\scriptscriptstyle \varGamma}$ on $\varGamma$ (Appendix~\ref{sec: FIO_symplecto_bundle}) 
is two. Since locally a Fourier integral operator on a bundle is a matrix of scalar Fourier integral operators this theorem can be applied here directly.
Hence, $\tilde{G}^{\pm} := G^{\mathrm{F, aF}} - G^{\pm}$ is a right parametrix which has this continuity property. 


Furthermore, the construction shows that $G^{\mathrm{F}} - \antiFeynParametrix$ and so $\tilde{G}^{+} - \tilde{G}^{-}$ is noncharacteristic at the diagonal of the set of lightlike covectors on $(M, g)$, as can be seen from the integral representation~\eqref{eq: causal_fundamental_sol_partial_derivative} of $F$ and wavefront set properties of $A, \tilde{A}$.  
Since $\square (\tilde{G}^{+} - \tilde{G}^{-})$ is smooth, it follows that the principal symbol satisfies a first order homogeneous differential equation along the lightlike geodesic (strip) on $T^{*} M$ by Theorem~\ref{thm: Hoermander_thm_25_2_4}.  
Therefore $\tilde{G}^{+} - \tilde{G}^{-}$ is noncharacteristic everywhere: $\mathrm{WF}' \big( \tilde{G}^{+} - \tilde{G}^{-} \big) = C$ and then it can be concluded that $\mathrm{WF}' \big( \tilde{G}^{+} - \tilde{G}^{-} \big) \supset C^{\pm}$ because of $\mathrm{WF}' \big( \tilde{G}^{\pm} \big) \subset \varDelta^{*} \cup C^{\pm}$. 
Finally, the fact $\varDelta^{*} = \WFPrime{I} = \WFPrime{\square \tilde{S}^{\pm}} \subset \WFPrime{\tilde{S}^{\pm}}$ with Lemma~\ref{lem: Hoermander_lem_26_1_16} complete the proof. 
\begin{remark} \label{rem: causal_propagator_NHOp_FIO}
	It is well-known that the wavefront sets of the retarded and the advanced Green's operators $G^{\mathrm{ret}, \mathrm{adv}}$ of a normally hyperbolic operator $\square : C^{\infty} (M; E) \to C^{\infty} (M; E)$ on a globally hyperbolic manifold $(M, g)$ are given by 
	\begin{equation} \label{eq: WF_advanced_retarded_Green_op_NHOp}
		\WFPrime{G^{\mathrm{ret}, \, \mathrm{adv}}} \subset \varDelta^{*} \cup C^{\mathrm{ret}, \, \mathrm{adv}}  
	\end{equation}
	where the retarded and the advanced geodesic relations are define by 
	\begin{equation} \label{eq: def_advanced_retarded_bicharacteristic_relation_NHOp}
		C^{\mathrm{ret}, \, \mathrm{adv}} := \{ (x, \xi; y, \eta) \in \dot{T}^{*} M \times \dot{T}^{*} M | g_{x}^{-1} (\xi, \xi) = 0, \exists s \in \mathbb{R} : (x, \xi) = \Phi_{s} (y, \eta), x \in J^{\pm} (y) \}.
	\end{equation}
	Another consequence of the proof of Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp}, is that the Pauli-Jordan-Lichnerowicz operator $G := G^{\mathrm{ret}} - \advGreenOp : C_{\mathrm{c}}^{\infty} (M; E) \to C_{\mathrm{sc}}^{\infty} (M; E)$ of $\square$ is a Fourier integral operator: 
	\begin{subequations}
		\begin{eqnarray}
			&& G \in I^{-3/2} \big( M \times M, C'; \End{E} \big), 
			\\ 
			&& \symb{G} = \frac{\mathrm{i}}{2}\sqrt{2 \pi |\rd_{C}|} u, 
			\\ 
			&& \Char{G} \cap C = \emptyset. 
		\end{eqnarray}
	\end{subequations}
	where $C$ is the geodesic relation~\eqref{eq: def_bicharacteristic_relation_NHOp}. Here $\rd_{C}$ the natural density on $C$ as for example explained in ~\cite[Thm. 7.1]{Strohmaier_AdvMath_2021}.
Moreover, $u$ is the unique section of the bundle $\End{E} \to C$ which is the identity map on the diagonal of the set of lightlike covectors on $(M, g)$ and parallel with respect to the Weitzenb\"ock connection.
This follows from Theorem~\ref{thm: Hoermander_thm_25_2_4} and the fact that the principal symbol of $\square$ vanishes on lightlike covectors 
(see~\cite[Thm. 6.6.1]{Duistermaat_ActaMath_1972} 
or e.g.,~\cite[Thm. 7.1]{Strohmaier_AdvMath_2021} 
for details given for scalar operator).  Note that the natural density we use on $C$ differs from ~\cite[p. 230]{Duistermaat_ActaMath_1972} by a factor of two because there the Hamiltonian flow of the principal symbol of $\square$ is used to parametrise lightlike geodesics, whereas we use the geodesic flow.
\end{remark}
\section{Hadamard states from Feynman parametrices}
In this section we will show that the unique Feynman parametrix $G^{\mathrm{F}}$~\eqref{eq: Feynman_parametrices_NHOp} constructed in Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp} can be turned into a fundamental solution $\FeynGreenOp$, the Feynman propagator for $\square$. 
Since Feynman propagator is state dependent, this is equivalent to construct a Hadamard $2$-point distribution.  
As remarked by 
Duistermaat-H\"{o}rmander~\cite[Thm. 6.6.2]{Duistermaat_ActaMath_1972},  
in the case of functions there is a remarkable positivity property. 
Namely, that one can modify $G^{\mathrm{F}}$ by a smoothing operator $R$ such that
\begin{equation}
	W + R := - \mathrm{i} (G^{\mathrm{F}} - \advGreenOp) + R
\end{equation}
is nonnegative in the sense that $(W + R) (\bar{u} \otimes u) \geq 0$ for any $u \in C_{\mathrm{c}}^{\infty} (M)$. 
This means, of course, that there exists a Feynman parametrix that has this property with $R = 0$. 
Thus, there exists a Feynman parametrix that satisfies $W \geq 0$ in the above sense. 
We state this here for the case of vector bundles.
\positivityFeynmanMinusAdvNHOp*
\begin{proof}

	Taking the adjoint $G^{\mathrm{F}} {}^{*} \square \equiv I ~\mathrm{mod}~ \PsiDO{-\infty}{M; E}$ of the Feynman parametrix $G^{\mathrm{F}}$ we find that $\WFPrime{G^{\mathrm{F}} {}^{*}} \subset \varDelta^{*} \cup C^{-}$ by Definition~\ref{def: Feynman_parametrix}.  
	Thus $G^{\mathrm{F}} {}^{*} - \antiFeynParametrix$ is smooth and so is $(G^{\mathrm{F}} - \advGreenOp) (G^{\mathrm{F}} - \advGreenOp)^{*}$. 
	Since $\mathrm{WF}' \big(- \mathrm{i} (G^{\mathrm{F}} - \advGreenOp) \big) \cap \mathrm{WF}' \big( \mathrm{i} (G^{\mathrm{F}} - \advGreenOp)^{*} \big) = \emptyset$, each of these is in $C^{\infty} (M; E)$. 
	Therefore $- \mathrm{i} (G^{\mathrm{F}} - \advGreenOp)$ can be normalised so that it is selfadjoint.  



	

	To address the positivity, first consider $D_{1}$ on $C^{\infty} (\Rn, \mathbb{C}^{N})$ (Section~\ref{sec: microlocalisation_NHOp}) and observe that $- \mathrm{i} F$ is positive by~\eqref{eq: causal_fundamental_sol_partial_derivative}:   
	\begin{equation} \label{eq: positivity_causal_fundamental_sol_partial_derivative}
		- \mathrm{i} F (u^{*} \otimes u) 
		= 
		\int u^{*} (y^{1}, y^{2}, \ldots, y^{n}) \, u (z^{1}, y^{2}, \ldots, y^{n}) \, \rd y \, \rd z^{1}  
		\geq 0. 
	\end{equation}



	
	To account $W$, we note that the Lagrangian distributions $A$ and $A^{*}$ in Proposition~\ref{prop: microlocalisation_NHOp} preserves wave front sets and thus, one obtains $(x_{0}, \xi^{0}; 0, \eta_{1} \rd y^{1}) \notin \WFPrime{W - A F A^{*}}$ from parametrix differences.  
	But it is not sufficient to construct a positive $W$ on a conic neighbourhood in $\dot{T}^{*} M \times \dot{T}^{*} M$ of the inverse image in $C$ of a neighbourhood of any point in $\Char{\square} / C$. 
	We need a refined version microlocalisation and hence respective homogeneous Darboux theorem. 



	
	On an $n$-dimensional globally hyperbolic manifold $(M, g)$ with Cauchy hypersurface $(\Sigma, h)$, the space of lightlike (unparametrised) geodesics is a $2n - 3$-dimensional smooth manifold which can be naturally identified with the cotangent unit sphere bundle $\mathbb{S}^{*} (\Sigma)$ of $\Sigma$ with respect to the induced metric $h$ on $\Sigma$~\cite{Low_NA_2001}. We choose $L$, an elliptic self-adjoint operator of order $-\frac{1}{2}$ with scalar principal symbol and vanishing subprincipal symbol as in Proposition~\ref{prop: microlocalisation_NHOp}. One can always choose the principal symbol of $L$ such that $X_{L \square L}$ is a complete vector field (see~\cite[Thm. 6.4.3 and p. 234]{Duistermaat_ActaMath_1972} 
	and also, e.g.,~\cite[Lem. 26.1.11 and Lem. 26.1.12]{Hoermander_Springer_2009} 
	and~\cite[Prop. 2.1]{Strohmaier_AdvMath_2021}).  
	Since the bundle of forward/backward lightlike covectors $T_{\pm}^{*} M \to M$ are connected components of the bundle of lightlike covectors $\Char{\square} \to M$, there exists a homogeneous diffeomorphism $\Char{\square} \to (\mathbb{S}^{*} (\Sigma) \oplus \mathbb{S}^{*} (\Sigma)) \times \mathbb{R} \times \mathbb{R}_{+}$ such that $X_{L \square L}$ is mapped to the vector field $\partial / \partial s$ when the variables in $\mathbb{S}^{*} (\Sigma) \times \mathbb{R} \times \mathbb{R}_{+}$ are denoted by $(\cdot, s, \cdot)$. 	In this setting, a refined Darboux theorem~\cite[Lem. 6.6.3]{Duistermaat_ActaMath_1972} 
	yields that, for every lightlike covectors $(x_{0}, \xi^{0})$ on $(M, g)$, there is an open conic coordinate chart $\big( \mathcal{U}, (x^{1}, \ldots, x^{n}; \xi_{1}, \ldots, \xi_{n}) \big)$ and a homogeneous symplectomorphism   
	\begin{equation} \label{eq: def_bijective_symplecto}
		\kappa : \mathcal{U} \to \mathcal{U}', \; (x, \xi) \mapsto 
		\kappa (x, \xi) := \big( y^{1} (x, \xi), \ldots, y^{n} (x, \xi); \eta_1  (x, \xi), \ldots, \eta_{n} (x, \xi) \big), 
	\end{equation}
	bijectively mapping $\mathcal{U}$ to an open conic neighbourhood $\mathcal{U}'$ of $\kappa (x_{0}, \xi^{0}) := (0, \eta^{0})$ in $\dot{T}^{*} \! \mathbb{R}^{n}$ such that 
	\begin{equation}
		\big( \kappa^{-1} \big)^{*} p = \xi_{1} \, \mathbbm{1}_{\End{E}}, 
	\end{equation}
	$\mathcal{U}' \cap \Char{D_{1}}$ is invariant under the translation along the $y^{1}$-axis, and $\Char{D_{1}} := \{ (y, \eta) \in \Rn \times \dot{\mathbb{R}}^{n} ~|~ y^{1} = 0 \}$ is symmetric with respect to the plane $y^{1} = 0$ and convex in the direction of $y^{1}$-axis. 
	Since $\big( \kappa^{-1} \big)_{*} X_{L^{*} \square L} = \mathbbm{1}_{\mathrm{Mat} (N)} \partial / \partial y^{1}$, it follows that $\mathcal{U} \cap \Char{\square}$ is invariant under $X_{\square}$. 



	Suppose that $\mathcal{V}_{1}, \mathcal{V}_{2}$ are closed conic neighbourhoods of $(x_{0}, \xi^{0})$ such that $\mathcal{V}_{1} \subset \mathcal{U}$ and $\mathcal{V}_{2} \subset \mathrm{int} (\mathcal{V}_{1})$ while $\mathcal{V}_{\iota} \cap \Char{\square}, \iota = 1, 2$ are invariant under the geodesic flow and that $\mathcal{V}'_{\iota} := \kappa (\mathcal{V}_{\iota})$ are convex in the $y^{1}$-direction and symmetric about the plane $y^{1} = 0$. 
	Remark~\ref{rem: global_microlocalisation_NHOp} then entails that there exists a properly supported Lagrangian distribution $A \in I^{-1/2} \big( M \times \Rn, \varGamma'_{1}; \Hom{\Rn \times \mathbb{C}^{N}, E} \big)$ such that $A^{*} A, A A^{*}$ both are pseudodifferential operators of order $-1$ and  
	\begin{subequations} \label{eq: microlocalisation_NHOp_refined}
		\begin{eqnarray}
			&& 
			\Char{A} \cap \varGamma'_{2} = \emptyset,
			\label{eq: microlocalisation_NHOp_1_refined} 
			\\ 
			&& 
			\ES{A^{*} A - I} \cap \mathcal{V}'_{2} = \emptyset  
			\Leftrightarrow
			\ES{A A^{*} - I_{E}} \cap \mathcal{V}_{2} = \emptyset, 
			\label{eq: microlocalisation_NHOp_2_refined}
			\\ 
			&& 
			\WFPrime{A^{*} \square A - D_{1}} \cap \varGamma_{2} = \emptyset, 
			\label{eq: microlocal_conjugate_NHOp_refined}
		\end{eqnarray}
	\end{subequations} 
	where $\varGamma_{1}$ and $\varGamma_{2}$ are the graphs of the restriction of $\kappa$ to $\mathcal{V}_{1}$ and to $\mathcal{V}_{2}$, respectively.



	
	It follows from Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp} and~\eqref{eq: microlocalisation_NHOp_refined} that we can microlocally conjugate $W$ to $F$ 
	(see~\cite[p. 237]{Duistermaat_ActaMath_1972} for the scalar version):  
	\begin{equation} \label{eq: FIO_6_6_14_bundle}
		\WFPrime{W - A F A^{*}} \cap (\mathcal{V}_{2} \times \mathcal{V}_{2}) = \emptyset, 
		\quad 
		\mathcal{V}_{2} \cap \Char{\square} \subset \big\{ \dot{T}_{\pm}^{*} M \big\}.  
	\end{equation} 




	\begin{lemma} 
	\label{lem: FIO_lem_6_6_5_bundle}
		For a normally hyperbolic operator $\square$ as in Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} and Lagrangian distributions $A$ given by~\eqref{eq: microlocalisation_NHOp_refined}, let $\varPsi \in \PsiDO{0}{M; E}$ such that $\ES{[\varPsi, \square]} \cap \Char{\square} = \emptyset$ and $\ES{\varPsi} \subset \mathcal{V}_{2}, 
			\Char{\square} \cap \mathcal{V}_{2} \subset  \dot{T}_{\pm}^{*} M $. Then 
		\begin{equation}
			W \varPsi \varPsi^{*} \equiv \varPsi A F A^{*} \varPsi^{*} ~\mathrm{mod}~\PsiDO{-\infty}{M; E}.
		\end{equation}
	\end{lemma}





	\begin{proof}
	 The idea of the proof is to use microlocalisation to reduce this to the special case $P=D_1$.
		To reach the conclusion our main task is then to compute the relevant commutators which is straightforward and exactly analogous to the scalar counterpart~\cite[Lem. 6.6.5]{Duistermaat_ActaMath_1972}. 
		So, we only sketch the main steps for completeness. 
		The ellipticity~\eqref{eq: microlocalisation_NHOp_1_refined} allows us to relate any properly supported $\varPsi \in \PsiDO{0}{M; E}$ and $\varPsi' \in \PsiDON{0}$ with $\ES{\varPsi} \subset \mathcal{V}_{2}$ and $\ES{\varPsi'} \subset \mathcal{V}'_{2}$ by 
		\begin{equation}
			\varPsi' \equiv A^{*} \varPsi A ~\mathrm{mod}~\PsiDON{-\infty}, 
			\quad 
			\varPsi \equiv A \varPsi' A^{*} ~\mathrm{mod}~\PsiDO{-\infty}{M; E}. 
		\end{equation} 
		A direct computation shows $\ES{[\varPsi', D_{1}]} \cap \Char{D_{1}} = \emptyset$ which implies that the derivative of the symbol of $\varPsi'$ with respect to $y^{1}$ is of order $-\infty$ in a neighbourhood of $\Char{D_{1}}$. 
		Denoting the convolution by the Dirac measure at $(h, 0, \ldots, 0)$ by $\tau_{h}$, we can rewrite~\eqref{eq: causal_fundamental_sol_partial_derivative} as $F = \int \tau_{h} \rd h$ and observe that $\tau_{h} \varPsi'  - \varPsi' \tau_{h} $ is of order $-\infty$ in a conic neighbourhood of $\Char{D_{1}}$. 
		Since $\Char{D_{1}}$ is invariant under translation in the $y^{1}$-direction, one has  
		\begin{equation}
			\WF{[\tau_{h}, \varPsi'] v} \cap \Char{D_{1}} = \emptyset, \quad \forall v \in \mathcal{E}'(\Rn, \mathbb{C}^{N}). 
		\end{equation}
		The result follows after the integration with respect to $h$ and a few algebraic manipulations.
	\end{proof}







	We observe that $\varPsi A F A^{*} \varPsi^{*}$ is nonnegative by~\eqref{eq: positivity_causal_fundamental_sol_partial_derivative} since $A A^{*}$ is nonnegative with respect to the hermitian form on $E$. 
	Thus the only thing left to conclude Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} is that the identity can be expressed as a sum of operators of the form $\varPsi \varPsi^{*}$ discussed in Lemma~\ref{lem: FIO_lem_6_6_5_bundle}, provided the prerequisites in this lemma are satisfied.  
	Aiming to show this, we consider a closed conic neighbourhood $\mathcal{V}_{3}$ of any lightlike covector $(x_{0}, \xi^{0})$ on $(M, g)$ where $\mathcal{V}_{3} \subset \mathrm{int} (\mathcal{V}_{2})$ and $\mathcal{V}_{3} \cap \Char{\square}$ is invariant under the geodesic flow.
One can prove that there exists $\varPsi \in \PsiDO{0}{M; E}$ with $\ES{\varPsi} \subset \mathcal{V}_{2}$ such that $\varPsi$ has a nonnegative (homogeneous) principal symbol equals to $1$ on $\mathcal{V}_{3}$ and $\ES{[\varPsi, \square]} \cap \Char{\square} = \emptyset$ (Lemma 6.6.6 in~\cite{Duistermaat_ActaMath_1972}  ).
	In other words, the hypotheses in Lemma~\ref{lem: FIO_lem_6_6_5_bundle} can be satisfied for $\varPsi$.
	Taking a suitable cover of the characteristic set one can use these operators to construct a family of operators $\varPsi_{\alpha} \in \PsiDO{0}{M; E}$ satisfying
	\begin{subequations}
		\begin{eqnarray}
			&& \ES{\varPsi_{\alpha}} \bigcap \Char{\square} \subset \dot{T}_{\pm}^{*} M , 
			\label{eq: WF_POU_1}
			\\ 
			&& \mathrm{ES} \bigg( I - \sum_{\alpha} \varPsi_{\alpha} \varPsi_{\alpha}^{*} \bigg) \bigcap  \dot{T}_{\pm}^{*} M  = \emptyset, 
			\label{eq: WF_POU_2}
			\\ 
			&& \ES{[\square, \varPsi_{\alpha}]} \bigcap \Char{\square} = \emptyset. 
			\label{eq: WF_commutator_NHOp_POU} 
		\end{eqnarray}
	\end{subequations}
	The construction of such a family that provides a microlocal partition of unity of the characteristic set of $\square$
	is carried out in detail in ~\cite{Duistermaat_ActaMath_1972}*{pp. 238}. The proof follows the usual stategy of summing the operators and multiplying from left and right by the root of a parametrix. One can verify that the hypotheses in Lemma~\ref{lem: FIO_lem_6_6_5_bundle} are preserved under this construction. 
	To finalise the proof of Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} we simply write 
	\begin{equation} \label{eq: Feynman_parametrix_P_mod_partial_deri}
		W 
		\equiv \sum_{\alpha} W \varPsi_{\alpha} \varPsi_{\alpha}^{*} 
		\equiv \sum_{\alpha} \varPsi_{\alpha} A_{\alpha} F A_{\alpha}^{*} \varPsi_{\alpha}^{*} 
		~\mathrm{mod}~ \PsiDO{-\infty}{M; E} 
	\end{equation}
	which terminates the proof as every term on the right hand side of~\eqref{eq: Feynman_parametrix_P_mod_partial_deri} is nonnegative with respect ot the hermitian form on $E$. 
\end{proof}
Of course, at the end one is interested in true solutions rather than parametrices. 
In fact there is a simple way, employing well-posedness of the Cauchy problem for normally hyperbolic operators, modifying the Feynman parametrix into a Feynman propagator preserving the aforementioned positivity. 
\existenceFeynmanGreenOpNHOp*
We postpone to proof this to make the following observations. 
Let $(\Sigma, h)$ be a Cauchy hypersurface of $(M, g)$ and assume $f_{0}, f_{1} \in C^\infty_{\mathrm{c}} (\Sigma, E |_{\Sigma})$. 
Since $M$ is globally hyperbolic, we will denote a point in $M$ as $(t, x) \in \mathbb{R} \times \Sigma$. 
Then we can define the distributions $f_{1} \otimes \delta_\Sigma, f_{0} \otimes \delta'_\Sigma \in \mathcal{E}'(M; E)$ by
\begin{equation} \label{eq: def_distributional_Cauchy_data_NHOp}
	f_{1} \otimes \delta_\Sigma (\phi)  
	:= 
	\int_\Sigma \phi (x) \, f_{1} (x) \, \dVolh (x), 
	\quad  
	f_{0} \otimes \delta'_\Sigma (\phi)  
	:= 
	\int_\Sigma (\partial_{N} \phi) (x) \, f_{0} (x) \, \dVolh (x), 
\end{equation}
where $\phi \in C_{\mathrm{c}}^{\infty} (M; E^{*})$ is any test section, $N$ is a future directed normal vector field on $\Sigma$, and $\dVolh$ is the volume form on $\Sigma$ induced by the Riemannian metric $h$. 
We denote by $\mathcal{W}_\Sigma \subset \mathcal{E}'(M; E)$ the span of the set of distributions of this form. By a duality argument both retarded and advanced fundamental solutions extend to continuous maps $\mathcal{E}'(M;E) \to \mathcal{D}'(M;E)$, and if $G = G^{\mathrm{ret}}- \advGreenOp$ and
$f = f_{0} \otimes \delta'_\Sigma + f_{1} \otimes \delta_\Sigma$, then $\phi = G(f)$ is a smooth solution of $\Box \phi=0$ with Cauchy data $(f_0,f_1)$ on $\Sigma$ 
(see e.g.~\cite[Lem 3.2.2]{Baer_EMS_2007}).

Assume that $w \Box$ and $\Box^{*} w$ are smooth for any bidistribution $w \in \mathcal{D}'(M \times M; E \boxtimes E^*)$.  
Let us denote by $\mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$, the set of compactly supported distributions with wavefront set contained in the conormal bundle $\mathrm{N}^*\Sigma$ of $\Sigma \subset M$. 
This space is endowed with a natural 
locally convex topology~\cite[p. 125]{Hoermander_ActaMath_1971} 
with respect to which $C_{\mathrm{c}}^{\infty} (M; E)$ is dense  
(see e.g.~\cite[Thm. 8.2.3]{Hoermander_Springer_2003} 
or the 
exposition~\cite[Sec. 4.3.1]{Strohmaier_Springer_2009} 
for details). 
The bidistribution $w$ can be defined as a sequentially continuous bilinear form on $\mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$.
Since $\mathcal{W}_\Sigma \subset \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$, the form $w$ is then also defined on $\mathcal{W}_\Sigma$.
\begin{lemma} \label{lem: positive_bisolution_NHOp}
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$ whose Cauchy hypersurface is $\Sigma$.
	Suppose that $E$ is endowed with a hermitian form and $w \in \mathcal{D}'(M \times M, E \boxtimes \bar{E}^*)$ be any bisolution of normally hyperbolic operators $\square$ on $E$. 
 	Then $w$ is nonnegative with respect to this form if and only if $w (\upphi^{*} \otimes \upphi) \geq 0$ for any $\upphi \in \mathcal{W}_\Sigma$. 
\end{lemma}
\begin{proof}
 First note that since $w$ is a bisolution it extends to a separately sequentially continuous bilinear form on $\mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E^*) \times \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$. We denote the extension by the same letter. Since this form is a bisolution it vanishes on $\Box^* C^\infty_{\mathrm{c}} (M; E^*) \times C^\infty_{\mathrm{c}} (M; E)$ and also on $C^\infty_{\mathrm{c}} (M; E^*) \times \Box C^\infty_{\mathrm{c}} (M; E)$. By sequential continuity and sequential density of $C^\infty_{\mathrm{c}}$ in $\mathcal{E}'_{\mathrm{N}^*\Sigma}$ it follows that $w$ vanishes on $\Box^* \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E^*) \times \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$ and on
 $\mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E^*) \times \Box \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$. We will show that any 
 $\phi \in C^\infty_0 (M; E)$ can be written as $\phi = \upphi + \Box h$, where $\upphi \in \mathcal{W}_\Sigma$ and $h \in \mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$, and the analogous statement holds for the dual bundle and the operator $\Box^*$. This obviously implies the statement. To see that this is indeed the case note that $\Box G \phi = 0$ and therefore, if $f_1,f_2$ are the Cauchy data of $G(\phi)$ on $\Sigma$ and $\upphi$ the corresponding element in 
 $\mathcal{W}_\Sigma$, then we have $G \phi = G(\upphi)$. This implies that $G(\phi - \upphi) = 0$ and therefore the retarded and advanced fundamental solutions applied to $\phi - \upphi$ are the same. Hence, $h= G^{\mathrm{ret}} (\phi - \upphi)$ must be a compactly supported distribution. By propagation of singularities this distribution must be again in $\mathcal{E}'_{\mathrm{N}^*\Sigma} (M; E)$ and we have $\phi - \upphi = \Box h$. The same proof works for the bundle $E^*$ and the operator $\Box^*$.
\end{proof}
Finally, we are prepared to show the existence of Feynman propagators of a normally hyperbolic operator on $E$. 
Moreover, if there exists a hermitian form on $E$ then Feynman propagators can be chosen nonnegative with respect to this form and hence implying existence of $2$-point bidistributions corresponding to Hadamard states. 
\begin{proof}[Proof of Theorem~\ref{thm: existence_Feynman_propagator_NHOp}]
 	Let $W$ be the parametrix for $\square$ constructed in Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp}. 
 	Since the wavefront set of $W$ does not intersect with the conormal bundle $\mathrm{N}^{*} \Sigma$ of any Cauchy hypersurface $\Sigma$ in $M$, we fix a $\Sigma$ and restrict $W$ and its normal derivative to create distributional Cauchy data on $\Sigma \times \Sigma$.
 	Suppose that $E$ is endowed with a hermitian form such that $W$ is nonnegative with respect to this form. 
 	Then the restriction of $W$ on $\Sigma$ is nonnegative as well. 
 	Next we construct $\omega \in \mathcal{D}' (M \times M; E \boxtimes \bar{E}^*)$ employing the solution operator $\mathcal{G}$ used in the proof of Lemma~\ref{lem: positive_bisolution_NHOp} with the same distributional Cauchy data as $W$ on $\Sigma \times \Sigma$. 
 	This will be a bisolution of $\square$ with the property 
 	\begin{equation}
 		\left( \omega (u^{*} \otimes u) \geq 0,  \quad \forall u \in \mathcal{E}'_{N^*\Sigma} (M; E) \right)
 \Leftrightarrow \omega \geq 0, 	\end{equation}
 	by the sequential density of $C_{\mathrm{c}}^{\infty} (\Sigma; E |_{\Sigma})$ in  $\mathcal{E}'_{\mathrm{N}^*(\Sigma)} (M; E)$ and Lemma~\ref{lem: positive_bisolution_NHOp}.
 	Note that $R = \omega - W$ is smooth since it solves the inhomogeneous problem with zero Cauchy data, which implies that $\FeynGreenOp := G^{\mathrm{F}} + \mathrm{i} R$ is the Feynman propagator with the required positivity property.
\end{proof}
\section{Feynman propagators for Dirac-type operators}
In this section we will construct Feynman propagators for Dirac-type operators. 
We recall that a first order linear partial differential operator $D$ on $E$ is called \textit{Dirac-type} if its principal symbol $\symb{D}$ satisfies the Clifford relation 
(see e.g.~\cite[Sec. 2.5 and Sec. 3.2]{Baer_Springer_2012} for details): 
\begin{equation}
	\big( \symb{D} (\xi) \big)^{2} = g^{-1} (\xi, \xi) \, \mathbbm{1}_{\End{E}}, \quad \forall \xi \in C^{\infty} (M; T^{*} M).  
\end{equation} 
Then the operator $ \square =  D^{2}$ is normally hyperbolic. 
Let $N$ be any future directed timelike tangent vector at any $x \in M$ and $\bullet$ denotes the Clifford multiplication. 
Suppose that $E$ is endowed with a nondegenerate sesequilinear form $\langle \cdot | \cdot \rangle$ such that 
\begin{equation} \label{eq: positive_definite_sesequilinear_form_Dirac_type_op}
	(\cdot | \cdot) := \big\langle - \mathrm{i} \symb{D} (N^{\flat}) \bullet \cdot \big| \cdot \big\rangle
\end{equation}
is fibrewise positive definite. Of course this positive definite inner product will depend on the choice of $N$.
If $D$ is formally selfadjoint with respect to $\langle \cdot | \cdot \rangle$ then the first condition in Theorem~\ref{thm: existence_Feynman_propagator_NHOp} is satisfied and we obtain Feynman propagator $\mathsf{S}^{\mathrm{F}}$ of $D$ immediately by setting $\mathsf{S}^{\mathrm{F}} := D \mathsf{G}^{\mathrm{F}}$ where $\mathsf{G}^{\mathrm{F}}$ is the Feynman propagator for $\square$. 
Nevertheless, the existence of Hadamard states cannot be deduced due to the positivity issue. 
One may wish to consider formally selfadjointness of $D$ with respect to $(\cdot | \cdot)$ so that the hypotheses of Theorem~\ref{thm: existence_Feynman_propagator_NHOp} are satisfied resulting existence of Hadamard states. 
But this turns out to be too restrictive to encompass all Dirac-type operators on globally hyperbolic manifolds (see Section~\ref{sec: twisted_Dirac_op} and Section~\ref{sec: Rarita_Schwinger_op} for concrete examples) and so we refrain to impose the condition. 
However, we will show below that there is a more direct way to construct a Feynman parametrix satisfying positivity, and hence a Hadamard state

To begin with, recall that the set of lightlike covectors $\Char{\square}$ on a globally hyperbolic manifold has only two connected components in three or higher  
dimensions\footnote{In 
	two dimensions, $\Char{\square}$ has four connected components.}:
the forward lightcone $\dot{T}_{+}^{*} M$ and the backward lightcone $\dot{T}_{-}^{*} M$.
Therefore, $C$ defined in Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp} has four different 
orientations\footnote{By 
	an orientation of a relation $C$, it is meant that any splitting of $C \setminus \Delta_{\Char{\square}}$ into a disjoint union of open $C^{+}, C^{-} \subset C$ which are inverse relations~\cite[p. 218]{Duistermaat_ActaMath_1972}.}
\cite[p. 218]{Duistermaat_ActaMath_1972} 
(see also~\cite[p. 540]{Radzikowski_CMP_1996}):  
$(C^{+}, C^{-}), (C^{-}, C^{+}), (C^{\mathrm{adv}}, C^{\mathrm{ret}}), (C^{\mathrm{ret}}, C^{\mathrm{adv}})$ corresponding to $\{ \Char{\square}, \dot{T}_{+}^{*} M, \dot{T}_{-}^{*} M, \emptyset \}$ where $C^{\pm}$ and $C^{\mathrm{ret}, \, \mathrm{adv}}$ are defined by~\eqref{eq: def_Feynman_bicharacteristic_relation_NHOp} and~\eqref{eq: def_advanced_retarded_bicharacteristic_relation_NHOp}, respectively. 
This enables to have a (nonunique) microlocal decomposition of Pauli-Jordan-Lichnerowicz operator $G$ as
\begin{subequations}
	\begin{eqnarray}
		&& 
		G \equiv G^{+} - G^{-} ~\mathrm{mod}~ \PsiDO{-\infty}{M; E}, 
		\\ 
		&& 
		G^{\pm} \square \equiv \square G^{\pm} \equiv 0 ~\mathrm{mod}~ \PsiDO{-\infty}{M; E}, 
		\\ 
		&& 
		G^{\pm} \in I^{-3/2} \big( M \times M, (C^{\pm})'; \End{E} \big), 
	\end{eqnarray}
\end{subequations}
Clearly, this is a slightly stronger version of Theorem~\ref{thm: exist_unique_Feynman_parametrix_NHOp} which can be proven exactly in the same fashion; cf Remark~\ref{rem: causal_propagator_NHOp_FIO}. 

We look for an analogue of the preceding results for Dirac-type operators $D$. 
In particular, $D$ admits unique advanced $S^{\mathrm{adv}} := D \advGreenOp$ and retarded $S^{\mathrm{ret}} := D G^{\mathrm{ret}}$ Green's operators on globally hyperbolic manifolds~\cite[Thm. 1]{Muehlhoff_JMP_2011}. 
We decompose the Pauli-Jordan-Lichnerowicz operator $S := S^{\mathrm{ret}} - S^{\mathrm{adv}}$ as $S^{+} + S^{-}$ up to smoothing operators where $S^{\pm} := \sum_{\alpha} Q_{\alpha}^{\pm} S Q_{\alpha}^{\pm}$ and $(Q_{\alpha}^{\pm})^{*} Q_{\alpha}^{\pm}$ are microlocal partition of unity such that $D$ commutes with $Q_{\alpha}^{\pm}$ up to smoothing operators and $\ES{Q_{\alpha}^{\pm}} \subset \{ (x, \xi; x, - \xi) ~|~ (x, \xi) \in \dot{T}_{\pm}^{*} M, g_{x}^{-1} (\xi, \xi) = 0 \}$. 
Consequently,  
\begin{subequations}
	\begin{eqnarray}
		&& 
		S \equiv S^{+} + S^{-} ~\mathrm{mod}~ \PsiDO{-\infty}{M; E}, 
		\label{eq: microlocal_decomposition_causal_propagator_Dirac_type_op}
		\\ 
		&& 
		S^{\pm} D \equiv D S^{\pm} \equiv 0 ~\mathrm{mod}~ \PsiDO{-\infty}{M; E},  
		\label{eq: distinghuished_parametrices_Dirac_type_op} 
		\\ 
		&& 
		S^{\pm} \in I^{-1/2} \big( M \times M, (C^{\pm})'; \End{E} \big), 
	\end{eqnarray}
\end{subequations}
by an application of~\eqref{eq: product_Lagrangian_dist_bundle} and the fact that $D$ is a local operator. 
The operator $Q_{\alpha}^{\pm}$ can be constructed by choosing any $q_{\alpha}^{\pm} \in S^{0 - [\infty]} \big( \dot{T}_{\pm}^{*} M; \End{E} \big)$ as its left total symbol such that $\mathrm{esssupp} (q_{\alpha}^{\pm}) = \{ (x, \xi) \in \dot{T}_{\pm}^{*} M ~|~ g_{x}^{-1} (\xi, \xi) = 0 \}$; in other words $q_{\alpha}^{\pm}$ is of order $-\infty$ near $y \in M$ and in a conic neighbourhood of $\eta$ for all $(y, \eta)$ in the complement of $\mathrm{esssupp} (q^{\pm})$. 
Furthermore, we pick $\supp{q_{\alpha}^{\pm}}$ slightly away from the projection of $\WF{Q_{\alpha}^{\pm}}$ on $M$ so that $[D, Q_{\alpha}^{\pm}]$ is smooth. 
We note that 

\begin{lemma}
	The Pauli-Jordan-Lichnerowicz operator $S$ of a Dirac-type operator on a globally hyperbolic manifold is nonnegative with respect to the nondegenerate sesquilinear form on $E$. 
\end{lemma}
\begin{proof} 
	Let $\psi, \phi \in C^{\infty} (M; E)$ having $\supp{\psi} \cap \supp{\phi}$ compact and let $\Omega$ be a compact subset of the globally hyperbolic manifold $(M, g)$ having smooth boundary $\partial \Omega$ with outward normal vector $N$ and volume element $\rd \mathsf{v}$. 
	Then the Green-Stokes formula yields 
	(see e.g.,~\cite[(1.7)]{Dimock_AMS_1982} or~\cite[(9)]{Baer_Springer_2012})
	\begin{equation}
		\int_{\Omega} \big( \langle D \psi | \phi \rangle - \langle \psi | D \phi \rangle \big) \mathrm{d} \mathsf{v}_{g} = \int_{\partial \Omega} (\psi | \phi) \, \rd \mathsf{v},  
	\end{equation} 
	where $\mathrm{d} \mathsf{v}_{g}$ is the Lorentzian volume form on $(M, g)$ and $(\cdot|\cdot)$ is given by~\eqref{eq: positive_definite_sesequilinear_form_Dirac_type_op}. 
	Setting $\psi := S^{\mathrm{ret}, \, \mathrm{adv}} u$ and $\phi := S u$ for any $u \in C_{\mathrm{c}}^{\infty} (M; E)$ the above formula entails $\langle u | S u \rangle \geq 0$ due to the compactness of $\supp{S^{\mathrm{ret}, \, \mathrm{adv}} u} \subset J^{\pm} (\supp{u})$. 
\end{proof}
Therefore, $\langle u | S^{\pm} u \rangle \geq 0$ by~\eqref{eq: microlocal_decomposition_causal_propagator_Dirac_type_op}. 
Since $D$ is a local operator, $\WFPrime{S^{\mathrm{adv}, \, \mathrm{ret}}} = \WFPrime{G^{\mathrm{adv}, \, \mathrm{ret}}} \subset \varDelta^{*} \cup C^{\mathrm{adv}, \, \mathrm{ret}}$ by~\eqref{eq: WF_advanced_retarded_Green_op_NHOp}. 
Note that, $C^{\pm} = C \cap (\dot{T}_{\pm}^{*} M \times \dot{T}_{\pm}^{*} M)$. 
Then inspecting the wavefront sets $C^{+}$~\eqref{eq: def_Feynman_bicharacteristic_relation_NHOp} and $C^{\mathrm{adv}}$~\eqref{eq: def_advanced_retarded_bicharacteristic_relation_NHOp} one can decompose 
$S^{-} = S^{\mathrm{adv}} - S^{\mathrm{F}}$ modulo smoothing operators~\cite[pp. 229-230]{Duistermaat_ActaMath_1972}, 
that is, $\langle u | W u \rangle$ is nonnegative where $W := - (S^{\mathrm{F}} - S^{\mathrm{adv}})$. 
Moreover, $S^{\mathrm{F}}$ is a unique parametrix because of $S^{\mathrm{adv}} D = D S^{\mathrm{adv}} = I$ and~\eqref{eq: distinghuished_parametrices_Dirac_type_op}. 
Hence, we arrive at the assertion which is the analogue of Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} for Dirac-type operators. 

As before, in order to turn $W$ into an exact distributional bisolution $\omega$, we are going employ the well-posedness of the Cauchy problem for Dirac-type operators on a globally hyperbolic manifold~\cite[Thm. 2]{Muehlhoff_JMP_2011}.  
For any $f \in C_{\mathrm{c}}^{\infty} (\Sigma; E |_{\Sigma})$, we define distributional Cauchy data $f \otimes \delta_{\Sigma}$ by 
(see~\cite[Prop. 2.4 (b)]{Dimock_AMS_1982} for the classical Dirac operator) 
\begin{equation} \label{eq: def_distributional_Cauchy_data_Dirac_type_op}
	f \otimes \delta_{\Sigma} (\phi) := \int_{\Sigma} \phi (t, x) \, ( N (x)  \bullet f (x)) \, \dVolh (x), 
\end{equation}
and by $\mathcal{X}_{\Sigma} \subset \mathcal{E}_{\mathrm{N}^{*} \Sigma}' (M; E)$, we denote the set of distributions of the above form where $\mathcal{E}_{\mathrm{N}^{*} \Sigma}' (M; E)$ is the set of compactly supported distributions whose wavefront set is contained in the conormal bundle $\mathrm{N}^{*} \Sigma$. 
As in the previous section, any bidistribution $w \in \mathcal{D}' (M \times M; E \boxtimes E^{*})$ such that $wD$ and $Dw$ are smooth in distributional sense, can be defined as a sequential continuous bilinear form on $\mathcal{X}_{\Sigma}$.  
\begin{lemma}
	Let $E \to M$ be a smooth complex vector bundle over a globally hyperbolic manifold $(M, g)$ whose Cauchy hypersurface is $\Sigma$. 
	Suppose that $w \in \mathcal{D}' (M \times M; E \boxtimes \bar{E}^{*})$ be a bisolution of Dirac-type operators $D$ on $E$ and that $E$ admits the positive-definite sesquilinear form~\eqref{eq: positive_definite_sesequilinear_form_Dirac_type_op}. 
	Then $\langle u | wu \rangle \geq 0$ if and only if $w (u^{*} \otimes u) \geq 0$ for all $u \in \mathcal{X}_{\Sigma}$. 
\end{lemma}
\begin{proof}
	The proof carries over its analogue Lemma~\ref{lem: positive_bisolution_NHOp} by obvious replace of $\square$, its Cauchy data~\eqref{eq: def_distributional_Cauchy_data_NHOp} and solution operator $\mathcal{G}$, by $D$, its Cauchy data~\eqref{eq: def_distributional_Cauchy_data_Dirac_type_op} and solution operator $\mathcal{S} : \mathcal{E}' (M; E |_{\Sigma}) \to \mathcal{D}' (M; E)$, respectively. 
\end{proof} 
Following exactly the route taken from Proposition~\ref{thm: positivity_Feynman_minus_adv_NHOp} to Theorem~\ref{thm: existence_Feynman_propagator_NHOp} the construction then yields a Feynman propagator that satisfies the positivity condition.
\HadamardBisolutionDiracTypeOp*
\section{Examples}
We list below several operators arising in the context of quantum field theories in globally hyperbolic manifolds $(M, g)$ to show that either they are of normally hyperbolic or Dirac-type. 
\subsection{Covariant Klein-Gordon operator}
Let $E \to M$ be a trivial line bundle, that is, its sections are just complex-valued smooth functions $C^{\infty} (M)$ on $M$. 
The covariant Klein-Gordon operator is defined by 
\begin{equation}
	\square := - \mathrm{div} \circ \mathrm{grad} + \mathrm{m}^{2} + \lambda \mathsf{R} : C^{\infty} (M) \to C^{\infty} (M), 
\end{equation}
where $\mathrm{m}^{2} \in \mathbb{R}_{+}$ is a parameter, physically interpreted as the mass-squared of a linear Klein-Gordon field, $\lambda$ is a coupling and $\mathsf{R}$ is the Ricci scalar of $M$. 
The special case of $\square$ when $\mathrm{m}, \lambda = 0$ is called the d'Alembert/relativistic wave operator. 
Since the natural inner product on $C^{\infty} (M)$ is positive-definite with respect to which $\square$ is formally selfadjoint 
(see e.g.~\cite[Example 2.13]{Baer_Springer_2012}), 
Theorem~\ref{thm: existence_Feynman_propagator_NHOp} applies and positive Feynman propagators can be constructed using this method.   
\subsection{Connection d'Alembert operator} 
Let $(E \to M, \nabla^{E})$ be a smooth vector bundle equipped with a connection $\nabla^{E}$ and $M$ is endowed with the Levi-Civita connection $\nabla^{T^{*} M}$ on its cotangent bundle $T^{*}M$. 
These connections induce a connection $\nabla$ on the vector bundle $T^{*} M \otimes E \to M$. 
The connected d'Alembert operator is then defined by the composition of following three maps  
\begin{equation}
	\square : C^{\infty} (M; E) 
	\stackrel{\nabla}{\to} C^{\infty} (M; T^{*} M \otimes E) 
	\xrightarrow{\nabla^{T^{*} M}} C^{\infty} (M; T^{*} M \otimes T^{*} M \otimes E) 
	\xrightarrow{\tr \otimes I} C^{\infty} (M; E), 
\end{equation} 
where $\tr_{g} : T^{*} M \otimes T^{*} M \to \mathbb{R}$ denotes the metric trace: $\tr_{g} \big( (x, \xi) \otimes (x, \eta) \big) := g_{x}^{-1} (\xi, \eta)$ and $I$ is the identity operator on $E$. 
A straightforward computation shows that this operator is normally hyperbolic 
(see e.g.~\cite[Example 2.14]{Baer_Springer_2012}). 
Then Theorem~\ref{thm: existence_Feynman_propagator_NHOp} entails that $\square$ admits a  Feynman propagator whereas the conclusion for positivity depends on the positivity of the inner product with respect to which $\square$ is formally selfadjoint. 
\subsection{Hodge-d'Alembert operator} 
\label{sec: Hodge_d_Alembert_op}
Let $E := \wedge^{k} T^{*} M$ be the bundle of $k$-forms. 
The Hodge-d'Alembert operator is defined by  
\begin{equation}
	\Box := - (\updelta \rd + \rd \updelta) : C^{\infty} (M; \wedge^{k} T^{*} M) \to C^{\infty} (M; \wedge^{k} T^{*} M),   
\end{equation}
where $\rd$ is the exterior derivative and $\updelta$ is coexterior differential. 
By the Bochner formula the Weitzenb\"ock connection is the Levi-Civita connection on $k$-forms and hence this connection is to be used for microlocalisation.
It can be shown that $\Box$ is a normally hyperbolic operator 
(see e.g.~\cite[Example 2.15]{Baer_Springer_2012}). 
There thus exists a Feynman propagator. Since the natural inner product with respect to which this operator is formally selfadjoint is not positive-definite unless $k=0$ or $k=n$, we cannot conclude nonnegativity directly. For example for $k=1$ nonnegativity is only expected on a subset as is usual for gauge theories. 
\subsection{Proca operator} 
\label{sec: Proca_op}
The Proca operator acting on covectors $A$ on $M$ is defined as   
\begin{equation}
	P := \updelta \rd + \mathrm{m}^{2} : C^{\infty} (M; T^{*}M) \to C^{\infty} (M; T^{*}M), 
	\quad \mathrm{m} \in \mathbb{R} \setminus \{0\},  
\end{equation}
which is not normally hyperbolic. 
Nevertheless, $P$ is equivalent to the normally hyperbolic operator $\updelta \rd + \rd \updelta + \mathrm{m}^{2}$ on $C^{\infty} (M; T^{*}M)$ together with the Lorenz constraint $\updelta A = 0$. 
Therefore, a Feynman propagator of $\updelta \rd + \rd \updelta + \mathrm{m}^{2}$ can be constructed by this method. The positivity issue is same as in Section~\ref{sec: Hodge_d_Alembert_op}. 
\subsection{Twisted Dirac operators} 
\label{sec: twisted_Dirac_op}
Let $S \to M$ be a spinor bundle over a globally hyperbolic spacetime $(M, g)$ 
admitting\footnote{A  
		spacetime $(M, g)$ admits a spin-structure if and only if the second Stiefel–Whitney class of $M$ vanishes. 
		We note that a spin structure always exists in a $4$-dimensional globally hyperbolic spacetime~\cite{Geroch_spin_JMP_1970}.
}  
a spin-structure. 
The Levi-Civita connection $\nabla^{\mathrm{LC}}$ on $(M, g)$ induces a spin connection $\nabla^{S}$ on $S$. 
Let $E \to M$ be a smooth vector bundle endowed with a connection $\nabla^{E}$.  
These two connections induce another connection $\nabla := \nabla^{S} \otimes \mathbbm{1}_{E} + \mathbbm{1}_{S} \otimes \nabla^{E}$ on the twisted spinor bundle $S \otimes E \to M$ and the twisted Dirac operator is defined by 
(see, e.g.~\cite{Baum_1981} for details) 
\begin{equation} \label{eq: def_twisted_Dirac_op}
	\slashed{D} := - \mathrm{i} \upgamma \circ g^{-1} \circ \nabla : C^{\infty} (M, S \otimes E) \to C^{\infty} (M, S \otimes E), 
\end{equation} 
where $\upgamma : C^{\infty} (M; TM) \to \End{S}$ is the Clifford mapping such that 
\begin{equation}
	\upgamma (X) \, \upgamma (Y) + \upgamma (Y) \, \upgamma (X) = 2 g (X, Y) \mathbbm{1}_{\End{S}}, \quad \forall X, Y \in C^{\infty} (M; TM)
\end{equation}
and the Clifford multiplication is given by $TM \otimes S \otimes E \ni X \otimes \psi \otimes u \mapsto (\upgamma (X) \psi) \otimes u \in S \otimes E$. 
The Schr\"{o}dinger~\cite{Schroedinger_GRG_2020}-Lichnwerowicz~\cite{Lichnerowicz_1963} 
formula entails  
(see e.g.~\cite[Example 1]{Muehlhoff_JMP_2011}) 
\begin{equation}
	\slashed{D}^{2} = \nabla^{*} \nabla + \frac{\mathsf{R}}{4} + \mathsf{F}
\end{equation}
where $\mathsf{R}$ is the scalar curvature of $M$ and $\mathsf{F}$ is the Clifford multiplied curvature of $\nabla^{E}$. 
Clearly, $\slashed{D}^{2}$ is a normally hyperbolic operator and thus $\slashed{D}$ is of  Dirac-type. This formula also shows that the Weitzenb\"ock connection is the twisted spin-connection on the bundle $S \otimes E$ and therefore this is the connection to be used for microlocalisation.
Since $\slashed{D}$ is formally selfadjoint with respect to the natural inner product $\langle \cdot | \cdot \rangle$ on $S \otimes E$, a Feynman propagator for $\slashed{D}$ exists by Theorem~\ref{thm: existence_Feynman_propagator_NHOp} yet the existence of Hadamard states cannot be concluded due to the fact that $\langle \cdot | \cdot \rangle$ is not positive-definite unless $(M, g)$ is a Riemannian manifold. 
However, $(M, g)$ is globally hyperbolic and thus, there exist smooth spacelike Cauchy hypersurfaces $\Sigma_{t} := \{t\} \times \Sigma$ with future-directed unit vector field $N$ normal to $\Sigma$ so that $\upgamma (N) \cdot$ defines the Clifford multiplication to achieve the positive-definite sesequilinear form~\eqref{eq: positive_definite_sesequilinear_form_Dirac_type_op}. 
Then Theorem~\ref{thm: Hadamard_bisolution_Dirac_type_op} is applicable and we have Feynman propagators with the desired positivity with respect to $\langle \cdot|\cdot \rangle$.  
\subsection{Rarita-Schwinger operator}
\label{sec: Rarita_Schwinger_op}
Let $S \to M$ be a complex spinor bundle over a spin manifold $M$ and $\upgamma$ be a the Clifford multiplication is given by 
$T^{*} M \otimes S \ni (\xi, u) \mapsto \upgamma (\xi^{\sharp}) u \in S$. 
Then, we have the representation theoretic splitting    
(see e.g.~\cite[Sec. 2.6]{Baer_Springer_2012}) 
\begin{equation}
	T^{*} M \otimes S = \iota (S) \oplus S^{3/2}, 
\end{equation}
where $S^{3/2}$ is defined as the kernel of the Clifford multiplication and the embedding $\iota$ of $S$ into $T^{*} M \otimes S$ is locally defined by 
$\iota (u) := - \frac{1}{n} \sum_{i=1}^{n} e_{i}^{*} \otimes \upgamma (e_{i}) u$ 
where $\{e_{i}\}_{1 \leq i \leq n}$ is an orthonormal basis of $TM$. 
Suppose that $\slashed{D} := \mathrm{i} (\mathbbm{1} \otimes \upgamma) \circ \nabla$ be the twisted Dirac operator on $T^{*} M \otimes S$ (cf.~\eqref{eq: def_twisted_Dirac_op}), then the Rarita-Schwinger operator is defined as~\cite[Def. 2.25]{Baer_Springer_2012} 
\begin{equation} \label{eq: def_Rarita_Schwinger_op}
	R := (I - \iota \circ \upgamma) \circ \slashed{D} : C^{\infty} (M, S^{3/2}) \to C^{\infty} (M, S^{3/2}). 
\end{equation}
$\Char{R}$ coincides with the set of lightlike covectors in dimensions $n \geq 3$  
(see e.g.~\cite[Lem. 2.26]{Baer_Springer_2012}), $R$ is formally selfadjoint whose Cauchy problem is well-posed when $(M, g)$ is globally hyperbolic,   
but $R^{2}$ is not a normally hyperbolic operator
(see e.g~\cite[Rem. 2.27]{Baer_Springer_2012}) 
and a positive-definite sesequilinear form~\eqref{eq: positive_definite_sesequilinear_form_Dirac_type_op} does not exist in $n \geq 3$~\cite[Example 3.16]{Baer_Springer_2012}. 

Originally, Rarita and Schwinger (in Minkowski spacetime) considered the twisted Dirac operator $D$ restricted to $S^{3/2}$ but not projected back to $S^{3/2}$~\cite[(1)]{Rarita_PR_1941} 
(see also~\cite[Sec. 2]{Homma_CMP_2019},~\cite[Prop. 2.7]{Wang_IndianaUMJ_1991} 
for Riemannian 
and~\cite[Rem. 2.28]{Baer_Springer_2012} 
for Lorentzian spin manifolds), 
that is,   
\begin{equation} \label{eq: constraint_Rarita_Schwinger}
	D |_{C^{\infty} (M; S^{3/2})} : C^{\infty} (M; S^{3/2}) \to C^{\infty} (M; T^{*} M \otimes S)
\end{equation}
to ensure the correct number of propagating degrees of freedom for spin-$3/2$ fields     
(see, for instance, the reviews~\cite{Sorokin_AIP_2005, Rahman_2015} for physical motivation and  different approaches used in Physics literature). 
This is an overdetermined system and so this constrained system limits possible curvatures of the spacetime  
(see e.g.~\cite[p. 856]{Homma_CMP_2019}): 
$(\mathsf{Ric} - \mathsf{R} g / n)^{*} u = 0$ 
where $\mathsf{Ric}$ and $\mathsf{R}$ are respectively the Ricci tensor and Ricci scalar curvature of the manifold $(M, g)$. 
In other words, Rarita-Schwinger operator~\eqref{eq: constraint_Rarita_Schwinger} exists only in Einstein spin manifolds in contrast to constraint Hodge-d'Alembert operator (cf. Section~\ref{sec: Hodge_d_Alembert_op}) and Proca operator (cf. section~\ref{sec: Proca_op}) discussed earlier. 
An application of the Weitzenb\"{o}ck formula~\cite{Branson_IJM_2002} shows that $R^{2}$ is normally hyperbolic on Einstein spin manifolds~\cite[Prop. 2.1]{Homma_CMP_2019}
\begin{equation}
	R^{2} = \Delta_{S^{3/2}} + \frac{n-8}{8n} \mathsf{R}, 
\end{equation} 
where $\Delta_{S^{3/2}}$ is the Bochner Laplacian. 
Therefore, Feynman propagators exist on a globally hyperbolic Einstein spin manifold by Theorem~\ref{thm: Hadamard_bisolution_Dirac_type_op} but the positivity issue is non-conclusive.
\subsection{Higher spin operators}
The straightforward attempts to generalise  
Dirac operator on Minkowski spacetime for arbitrary spin~\cite{Dirac_PRSA_1936} 
in 
curved spacetimes\footnote{Not 
	necessarily be globally hyperbolic.}
leads to difficulties
(see, e.g,~\cite[p. 324]{Illge_AnnPhys_1999} for a panoramic view 
and the reviews~\cite{Sorokin_AIP_2005, Rahman_2015}). 
A crucial advancement came through 
Buchdahl operator (in Riemannian manifold)~\cite{Buchdahl_JPA_1982} 
(see~\cite[Exam. 2.24]{Baer_Springer_2012} for a Lorentzian formulation) 
whose square turns out to be a normally hyperbolic operator~\cite[p. 8]{Baer_Springer_2012}, 
yet the minimum coupling principle seems violated and a "by hand" proposal is required in the original idea of Buchdahl. 
These minor imperfections were cured by 
W\"{u}nsch~\cite{Wuensch_GRG_1985},  
by 
Illege~\cite{Illge_ZAA_1992, Illge_CMP_1993} 
and by 
Illege and Schimming~\cite{Illge_AnnPhys_1999} 
for the massive case deploying the $2$-spinor formalism in $4$-dimensional curved spacetimes. 
In particular, the square of Buchdahl operator (as modified by W\"{u}nsch and Illge) is a normally hyperbolic operator on certain twisted bundles. 
In contrast, there are a few open questions for the massless case~\cite{Frauendiener_JGP_1999} 
(see also the reviews~\cite{Sorokin_AIP_2005, Rahman_2015} for the contemporary status and other formulations used in Physics literature). 
Hence, the existence of a Feynman propagator is evident either by Theorem~\ref{thm: existence_Feynman_propagator_NHOp} or Theorem~\ref{thm: Hadamard_bisolution_Dirac_type_op} but the issue of positivity is non-conclusive at this stage. 
\section*{Acknowledgement}  
OI is supported by a Leeds International Doctoral Studentship.
