Global propagator for the massless Dirac operator and spectral asymptotics

We construct the propagator of the massless Dirac operator $W$ on a closed Riemannian 3-manifold as the sum of two invariantly defined oscillatory integrals, global in space and in time, with distinguished complex-valued phase functions. The two oscillatory integrals -- the positive and the negative propagators -- correspond to positive and negative eigenvalues of $W$, respectively. This enables us to provide a global invariant definition of the full symbols of the propagators (scalar matrix-functions on the cotangent bundle), a closed formula for the principal symbols and an algorithm for the explicit calculation of all their homogeneous components. Furthermore, we obtain small time expansions for principal and subprincipal symbols of the propagators in terms of geometric invariants. Lastly, we use our results to compute the third local Weyl coefficients in the asymptotic expansion of the eigenvalue counting functions of $W$.

the Riemannian density. Let us clarify straight away why, when dealing with the massless Dirac operator, we restrict our analysis to the 3-dimensional case. The reason is twofold: on the one hand, dimension three is physically meaningful in that it represents the first step towards a potential future analysis of the relativistic 3+1-dimensional setting, and on the other hand, our method requires the eigenvalues of the principal symbol of our operator to be simple, cf. Section 3, which is not the case for the massless Dirac operator in higher dimensions.
Let e j , j = 1, 2, 3, be a positively oriented global framing, i.e. a set of three orthonormal smooth vector fields 1 whose orientation agrees with the orientation of the manifold.
In chosen local coordinates x α , α = 1, 2, 3, we will denote by e j α the α-th component of the j-th vector field. Throughout this paper we use Greek letters for holonomic (tensor) indices and Latin for anholonomic (frame) indices. We adopt Einstein's convention of summation over repeated indices. Let Here H 1 is the usual Sobolev space of functions which are square integrable together with their first partial derivatives.
In relativistic particle physics the massless Dirac equation is often referred to as the Weyl equation, which explains our notation. Our operator W appears as the result of separating out the time variable in the relativistic Weyl equation, see [19] for details. Henceforth, we refer to the massless Dirac operator simply as the Dirac operator, which conforms with the terminology adopted in differential geometry. The Dirac operator (1.4) is uniquely determined by the metric and spin structure modulo an SU(2) gauge transformation.
The Dirac operator is symmetric with respect to the L 2 inner product where dx = dx 1 dx 2 dx 3 . Furthermore, a simple calculation shows that it is elliptic 2 . It is well known that the Dirac operator is self-adjoint and its spectrum is discrete, accumulating to +∞ and to −∞. Let λ k be the eigenvalues of W and v k the corresponding orthonormal eigenfunctions, k ∈ Z. The choice of particular enumeration is irrelevant for our purposes, but what is important is that eigenvalues are enumerated with account of their multiplicity. Note that the Dirac operator has the special property that it commutes with the antilinear operator of charge conjugation see [24, Appendix A] for details, and this implies that eigenvalues have even multiplicity. (1.11) The Dirac propagator can be written as a sum of three operators (1.12c) We call the operators (1.12a), (1.12b) and (1.12c) positive, zero mode and negative propagators, respectively. These are time-dependent partial isometries. Note that the operator U 0 is nontrivial only if the Dirac operator has zero modes (i.e. if zero is an eigenvalue). We define the positive (+) and negative (−) local counting functions as N ± (y; λ) := 0 for λ ≤ 0, 0<±λ k <λ [v k (y)] * v k (y) for λ > 0. (1.14) The functions (1.14) are the 'global' counting functions, the only difference with the usual definition [45] being that we count the positive and negative eigenvalues separately. Letμ : R → C be a smooth function such thatμ = 1 in some neighbourhood of the origin and suppμ is sufficiently small. Here 'sufficiently small' means that suppμ ⊂ (−T 0 , T 0 ), where T 0 is the infimum of lengths of all the geodesic loops originating from all the points of the manifold.
Definition 1.4. We call local Weyl coefficients the smooth functions c ± k (y) appearing in the asymptotic expansions (1.17).
(v) An important topic in the spectral theory of first order elliptic systems is the issue of spectral asymmetry [1,2,3,4]. Let us mention that to observe spectral asymmetry for our Dirac operator one as to go as far as the sixth Weyl coefficients. This follows from the fact [12,31] that the eta function is holomorphic in the complex half-plane Re s > −2 and has its first pole at s = −2.
The value of the residue of the eta function at s = −2, which was computed explicitly by Branson and Gilkey [14], describes the difference between the sixth (global) Weyl coefficients.
Our paper has two main objectives.
Objective 1 Construct the propagators U ± (t) explicitly, modulo integral operators with infinitely smooth kernels, and do so as a single invariantly defined oscillatory integral global in space and in time.
Global propagator for the massless Dirac operator Remark 1.6. One cannot, in general, identify the third Weyl coefficient by looking at the asymptotic behaviour of the unmollified counting function. In order to illustrate this point, let us consider the 3-torus equipped with standard flat metric. Already in this simple case the mathematical statement is false. This fact can be established by writing down the eigenvalues explicitly as in [24, Appendix B] and using standard results [32] from analytic number theory.

Main results
The study of Dirac operators in curved space, arguably the most important operators from the point of view of physical applications alongside the Laplacian, has a long and noble history in the mathematical literature. Excellent introductions to the subject can be found in [39] and [30]. Due to the physical significance of the topic, numerous researchers have contributed to our current understanding of the spectrum of Dirac operators on Riemannian manifolds. One can ask, for example, how the eigenvalues behave under perturbations of the metric [13,28,25], how the spectrum depends on the spin structure [10], whether zero modes exist [8], et cetera.
Later in this paper we will be concerned with the study of the asymptotic behaviour of large (in modulus) eigenvalues of the Dirac operator on a closed 3-manifold. In the case of scalar elliptic operators, such as for example the Laplace-Beltrami operator, a wide range of classical techniques are available in the literature to compute spectral asymptotics. However, if one is interested in a first order system, whose spectrum is, in general, not semi-bounded, the heat kernel method can no longer be applied, at least in its original form, and even resolvent techniques require major modification [7]. A very natural approach in this case is the so-called wave method, going back to Levitan [40] and Avakumovic [5], which involves recovering information about the eigenvalue counting function from the behaviour of the wave propagator, see (1.11). How this can be done will be explained in greater detail later on. This partly motivates our interest in the Dirac propagator (1.9), which is also of interest on its own. Of course, the hyperbolic Cauchy problem (1.10) for W lies at the heart of relevant applications in theoretical physics (e.g., the mathematical description of neutrinos/antineutrinos in curved space).
In order to construct the propagator (1.9) precisely, one needs to know all eigenvalues and eigenfunctions of W , which is unrealistic for a generic Riemannian manifold. However, microlocal techniques allow one to construct the propagator (1.9) approximately, modulo an integral operator with smooth integral kernel. This fact is well-known and an extensive discussion can be found, for instance, in [33].
There are, however, several issues with this classical construction: (i) it is not invariant under changes of local coordinates, (ii) it is local in space and (iii) it is local in time. The last issue, locality in time, is especially serious: it is to do with obstructions associated with caustics. In practice, constructing a propagator locally in time means that for large times one has to use compositions

Global propagator for the massless Dirac operator
The propagator U(t) is a special case of a Fourier integral operator and it is known that handling compositions of such operators is a daunting task.
Our goal is to construct the Dirac propagator explicitly, in a global -i.e., as a single oscillatory integral -and invariant (under change of coordinates and gauge transformations) fashion. The key idea, originally proposed by Laptev, Safarov and Vassiliev in [38] and further developed in [45], is to use Fourier integral operators with complex-valued, as opposed to real valued, phase function. Crucially, this allows one to circumvent the topological obstructions due to caustics.
Our work partly builds upon [23] and [24]. In [23], using the wave method, Chervova, Downes and Vassiliev obtained an explicit formula for the second Weyl coefficient of an elliptic self-adjoint first order pseudodifferential matrix operator, fixing thirty years of incorrect or incomplete publications in the subject, see [23,Section 11]. In [24] the same authors applied the results from [23] to the Dirac operator. Unlike the current paper, the approach from [23] is not geometric in nature and the complexity of phase functions is not actually put to use. Note that some results from [24] were improved by Strohmaier and Li in [41], where the authors studied the second term of the mollified spectral counting function of Dirac type operators and characterised operators in this class with vanishing second Weyl coefficient.
A fully geometric global construction of the (scalar) wave propagator e −it √ −∆ on closed Riemannian manifolds, as a single oscillatory integral with complex-valued phase function, was recently proposed by the authors and Levitin in [17], and subsequently extended to the Lorentzian setting in [16]. The publication [17] is the starting point of the current paper. The extension of the results of the current paper to globally hyperbolic Lorentzian manifolds is carried out in [18].
Our main results are as follows.
1. We present a global construction of each of the two propagators, the positive propagator U + (t) and the negative propagator U − (t), as a single invariantly defined oscillatory integral, global in space and in time, with distinguished complex-valued phase function (Theorem 3.3, Definition 5.1, Definition 5.3). We provide a closed formula for the principal symbols of the propagators (Theorem 6.1) and an algorithm for the calculation of the subprincipal symbols and all asymptotic components of lower degree of homogeneity in momentum (subsection 3.3).
2. We give an explicit small time expansion of principal and subprincipal symbols of positive and negative propagators in terms of geometric invariants (Theorem 7.13).
3. We compute the third local Weyl coefficients in the asymptotic expansion of the two eigenvalue counting functions (1.13) (Theorem 8.1).
Along the way we prove a number of results about general first order elliptic systems and invariant representations of pseudodifferential operators on manifolds. Note that the third Weyl coefficients can, in principle and with some work, be also obtained by a different method using results available in the literature, see Remark 8.3.
Our paper is structured as follows. In Section 3 we explain how to construct explicitly positive and negative propagators for a general first order elliptic self-adjoint (pseudo)differential matrix operator, with a rigorous mathematical justification.

Global propagator for the massless Dirac operator
In Section 4 we deal with the delicate issue of invariant descriptions of pseudodifferential operators acting on scalar functions. In particular, we examine the relation between our g-subprincipal symbol and the standard notion of subprincipal symbol for operators acting on half-densities.
In Section 5 we apply the results from Section 3 to the Dirac operator. A formula for the principal symbol of positive and negative Dirac propagators is provided in Section 6, whereas small time expansions for principal and subprincipal symbols of positive and negative propagators are obtained in Section 7. Our final results are expressed in terms of geometric invariants: curvature of the Levi-Civita connection associated with the metric g and torsion of the Weitzenböck connection generated by the framing defining the Dirac operator.
In Section 8 we use the results from Section 7 to compute the third local Weyl coefficients for the Dirac operator.
Finally, in Section 9 we apply our techniques to two explicit examples: M = S 3 , where formulae are isotropic in momentum, and M = S 2 × S 1 , where they are not.
The paper is complemented by two appendices, containing background material and technical proofs.

Preliminary results for general first order systems
In this section we will consider a broader class of first order systems and we will prove fairly general results, which will be later applied to the special case of the Dirac operator. In doing so, we will need some of the technology developed in [23]. The setting of our analysis is somewhat different from that in [23], in that our operators are differential, as opposed to pseudodifferential (see also Remark 3.11), and act on scalar functions on a Riemannian manifold, as opposed to half-densities on a manifold with no metric structure. In particular, the change of the space in which the operator acts raises delicate issues concerning the invariance of the mathematical objects involved. For these reasons we provide here a modified version of some of the results from [23], adapted to the setting of our paper.
Throughout this section, M will be a smooth connected closed Riemannian manifold of dimension d ≥ 2.
Let A be an elliptic symmetric (with respect to (1.8)) first order m × m matrix differential operator acting on m-columns of smooth complex-valued scalar functions We denote by h (j) (x, ξ) the eigenvalues of A prin (x, ξ) enumerated in increasing order, with positive index j = 1, 2, . . . , m + for positive h (j) (x, ξ) and negative index j = −1, −2, . . . , −m − for negative h (j) (x, ξ). We assume that the eigenvalues of the principal symbol A prin are simple. Clearly, m = m + + m − , because the ellipticity condition det A prin (x, ξ) = 0 ensures that all eigenvalues are nonzero. In fact, as our operator is differential, one can show [29, Remark 2.1] that m can only be even and that we have Furthermore, the eigenvalues h (j) of the principal symbol and the corresponding normalised eigenvectors v (j) possess the symmetry Global propagator for the massless Dirac operator Under the above assumptions the spectrum of A is discrete and accumulates to +∞ and to −∞. We denote eigenvalues and orthonormalised (smooth) eigenfunctions of A by λ k and v k , respectively, enumerated with account of their multiplicity.
By replacing W with A, one can define the 'full' propagator U A (t) for A via (1.11), as well as the positive, zero mode and negative propagators via (1.12a)-(1.12c), which we denote by U + A (t), U 0 A and U − A (t), respectively. Each eigenvalue h (j) (x, ξ) of the principal symbol can be interpreted as a Hamiltonian on the cotangent bundle. The corresponding Hamiltonian flow (x (j) (t; y, η), ξ (j) (t; y, η)), i.e. the (global) solution to Hamilton's equationṡ with initial condition (x (j) (0; y, η), ξ (j) (0; y, η)) = (y, η), generates a Lagrangian manifold to which one can, in turn, associate a global Lagrangian distribution. See [17,Section 2] and references therein for details. In particular, the singularities of the solution to the initial value problem propagate along Hamiltonian trajectories generated by the eigenvalues of A prin .

Positive and negative propagators: an abstract approach
Our aim is to show that U + A (t) and U − A (t) can be separately approximated by a finite sum of global oscillatory integrals. Before doing so, let us state and prove an abstract preparatory theorem.
if for every α > 0, every k ∈ N and every linear partial differential operator P with infinitely smooth coefficients of order k on M x × M y there exists a positive constant C α,P such that |P v| ≤ C α,P |λ| −α for ± λ > 1 , Theorem 3.2. Let (T − , T + ) ⊆ R be an open interval (possibly, the whole real line) and let u + (t, x, y), u − (t, x, y), u + (t, x, y) and u − (t, x, y) be elements of Furthermore, assume that for every ζ ∈ C ∞ 0 (T − , T + ) we have Global propagator for the massless Dirac operator Proof. Let ζ ∈ C ∞ 0 (T − , T + ). Multiplying (a) by ζ(t) we get (3.5) Applying the inverse Fourier transform F −1 t→λ to (3.5), letting λ → +∞ and using assumptions (b) and (c) we obtain Here, when dealing with the remainder from (3.5), we used the fact that the Fourier transform of a compactly supported smooth function is rapidly decreasing. The compactness of M ensures a uniform estimate in the spatial variables. Furthermore, (b) and (c) immediately imply Combining (3.6) and (3.7) we arrive at As ζ ∈ C ∞ 0 (T − , T + ) in the above formula is arbitrary, we conclude that A similar argument gives

Construction of positive and negative propagators
Theorem 3.3. Let A be an elliptic self-adjoint first order pseudodifferential operator acting on m-columns of scalar functions over M, whose principal symbol has simple eigenvalues. The positive and negative propagators U + A (t) and U − A (t) can be written, modulo an infinitely smoothing operator, as a finite sum of oscillatory integrals, global in space and in time. More precisely, we have where U (j) x;y,η) a (j) (t; y, η) χ (j) (t, x; y, η) w (j) (t, x; y, η) ( · ) ρ(y) dy dη (3.10) and Global propagator for the massless Dirac operator we mean that the operator on the LHS is equal to the operator on the RHS up to an integral operator with infinitely smooth integral kernel; • ( · ) is meant for insertion of f 0 (y) when computing (U (j) is an element in the class of polyhomogeneous symbols of order zero with values in m × m complex matrices, which means that a (j) admits an asymptotic expansion in components positively homogeneous in momentum, x; y, η) = 1 on the intersection of {(t, x; y, η) | |h (j) (y, η)| ≥ 1} with some conical neighbourhood of {(t, x (j) (t; y, η); y, η)}, where the smooth branch of the complex root is chosen in such a way that Remark 3.4. Note that the weight w (j) is the inverse of a smooth density in the variable y and a smooth scalar function in all other variables. The powers of the Riemannian density ρ in (3.12) are chosen in such a way that the symbol a (j) and the integral kernel of the operator (3.10) are scalar functions in all variables. The fact that the symbol is a genuine scalar function on R × T ′ M is a crucial feature of our construction. Taking the square and then extracting the fourth root in (3.12) serves the purpose of making the weight invariant under inversion of a single coordinate x α or a single coordinate y α . Note, however, that if one works on an orientable and oriented manifold, then one can simplify (3.12) to read Remark 3.5. The existence of a phase function satisfying conditions (i)-(iv) is a nontrivial matter. In fact, the space of phase function satisfying these conditions is nonempty and path-connected, see [38, Lemmata 1.4 and 1.7].
Remark 3.6. Let us emphasise that a phase function ϕ (j) satisfying conditions (i)-(iv) from Theorem 3.2 automatically satisfies see, e.g., [45,Subsection 2.4.1]. The equation is known in the literature as eikonal equation. Note that when x = x (j) (t; y, η) formula (3.15) turns into (3.14). In the classical approach to the construction of hyperbolic propagators, (3.15) is required to be satisfied in some open neighbourhood of This is a fundamental difference with the approach adopted in the current paper, where (3.15) is only required to be satisfied 'along the Hamiltonian flow', i.e., one only needs (3.14). Indeed, there is no open neighbourhood of (3.16) where the special phase functions that will be introduced and used in Section 5 -the Levi-Civita phase functionssatisfy (3.15). Relaxing the requirements on our phase functions is needed in order to accommodate an imaginary part and, consequently, circumvent obstructions arising from caustics.
Proof of Theorem 3.3. Suppose that we have constructed the symbols a (j) appearing in the oscillatory integrals (3.10) so that How to achieve this will be explained in subsection 3.3. Put so that the Schwartz kernel of the operator U A (t) reads

Global propagator for the massless Dirac operator
Let u(t, x, y), u + (t, x, y) and u − (t, x, y) be the Schwartz kernels of the operators U A (t), U + A (t), and U − A (t), respectively. Formulae (3.18a) and (3.18b) imply This fact can be established as follows. Let From the construction algorithm, we know that Here the superscript in A (x) indicates that the differential operator A acts in the variable x. Using functional calculus, we can write the functions u ∞ , f and ζ in terms of the eigenfunctions of A as Let ζ be the operator with integral kernel ζ(x, y), Then The operator A l ζA n is a pseudodifferential operator of order −∞, so formula (3.26) and the fact that λ k ∼ k 1/d when k → ∞ allow one to conclude that the c jk decay faster than any power of j and k as j, k → ∞. A similar argument shows that the b jk (t) and their time derivatives decay faster than any power of j and k as j, k → ∞ uniformly over any bounded open interval in R. Formula (3.25) now tells us that the same is true for the a jk (t). This, in turn, implies that the series on the RHS of (3.22) defines a function u ∞ (t, x, y) which is smooth in all variables. So we arrive at (3.19), which gives us assumption (a) in Theorem 3.2 with (T − , T + ) = R.
Resorting to standard stationary phase arguments -see, e.g., [45, Appendix C] -and using the properties (i)-(iv) of our phase functions, it is easy to see that u ± and u ± satisfy assumptions (b) and (c) of Theorem 3.2. Hence, Theorem 3.2 gives us (3.8) and (3.9).
The fact that the construction is global in time is guaranteed by [38, Lemma 1.2].
Remark 3.7. If one is prepared to give up globality in time, Theorem 3.3 and the corresponding proof can be adapted in a straightforward manner to the more customary case of real-valued -as opposed to complex-valued -phase functions. This is achieved by prescribing the phase functions to take values in R, dropping condition (iv) and replacing everywhere in the statement and in the proof the time domain R with the interval The values of T ± depend on the choice of particular real-valued phase functions, but we always have T − < 0 < T + . Observe that Theorem 3.2 was formulated in such a way that it covers both the case of real-valued and complex-valued phase functions.
The reader will have noticed that the zero mode propagator U 0 A does not appear in our construction. This is due to the fact that, clearly, We end this subsection with the observation that, thanks to the presence of the weight w (j) in formula (3.10), the scalar matrix-function a

The algorithm
The integral kernel (3.13) of U A (t) can be constructed explicitly as follows.
Step 1. Choose a phase function ϕ (j) compatible with Theorem 3.3. We will see later on that for the special case of the Dirac operator we can identify a distinguished phase function, the Levi-Civita phase function. Furthermore, set χ (j) ≡ 1. In fact, the purpose of the cut-off is to localise integration in a neighbourhood of the h (j) -flow and away from the zero section: different choices of χ (j) result in oscillatory integrals differing by an infinitely smooth function.
Global propagator for the massless Dirac operator Step 2. Act with the operator −i∂ t + A (x) on the oscillatory integral (3.13). This produces a new oscillatory integral By making use of the fact that ϕ (j) and w (j) are positively homogeneous in momentum η of degree 1 and 0, respectively, one can write down an asymptotic expansion for the amplitude a (j) in components positively homogeneous in momentum: Step 3. As u (j) (t, x, y) is to be the (distributional) solution of the hyperbolic equation one would like to impose the condition a (j) (t, x, y, η) = 0. However, the amplitude a (j) , unlike the symbol a (j) , depends on x, and doing so would result in an unsolvable system of partial differential equations (PDEs). The current step consists in excluding the dependence of a (j) on x by means of a procedure known as reduction of the amplitude, to the end of reducing the system of PDEs to a system of ordinary differential equations instead. Put 3 Furthermore, the operators S are the pushforward of partial derivatives ∂/∂ x α along f −1 . Hence, the operators L α commute because the partial derivatives ∂/∂ x α commute. An adjustment of the above argument to our setting with f = ϕ (j) η (and with account of the fact that ϕ (j) is complexvalued) provides an alternative explanation for the commutation of the operators L The amplitude-to-symbol operator is defined as When acting on a function positively homogeneous in momentum, the operator S (j) −k excludes the dependence on x and decreases the degree of homogeneity by k.
The reduction of the amplitude is achieved by replacing the amplitude a (j) in (3.29) by The oscillatory integral (3.29) only by an infinitely smooth function. We refer the reader to [17, Appendix A] for further particulars and detailed proofs concerning the amplitude-to-symbol operator. Step Equations (3.32), combined with the initial conditions stemming from the constraint yield a hierarchy of (matrix) transport equations for the homogeneous components a Let us make a few remarks warranted by formula (3.33). The m oscillatory integrals appearing on the RHS of (3.8) and (3.9) are not independent of one another, but they 'mix' at t = 0 via the initial condition (3.33). Now, satisfying (3.33) involves representing the identity operator on C ∞ (M; C m ) in a somewhat nonstandard fashion, as ). In terms of the symbols a (j) , the initial condition (3.33) reads a (j) (0; y, η) = s (j) (y, η).

Global propagator for the massless Dirac operator
From the fact that the principal symbol of the identity operator is the identity matrix it follows that j a (j) However, obtaining formulae for subleading components s (j) −1 is already a challenging task, see [23, subsection 4.2]. In general, lower order components of s (j) depend in a nontrivial manner on the eigenvalues and eigenprojections of the matrix-function A prin (x, ξ) and on the choice of phase functions ϕ (j) .
The invariant representation of the identity operator -and, more generally, of pseudodifferential operators -on manifolds is not a well-studied subject. An initial analysis of the scalar case was carried out in [17,Section 6]. For the case of the Dirac operator a more detailed examination of (3.34) will be provided in subsection 5.2. A more extensive analysis of (3.34) for a general operator A is carried out in [20,21,15].
Remark 3.11. All statements and results presented in this section carry over verbatim to the case where A is an elliptic symmetric first order m × m matrix pseudodifferentialas opposed to differential -operator, with the following exceptions: • formulae (3.1) and (3.2) have to be dropped as they are no longer true; • 'Step 2.' in subsection 3.3 has to be modified to take into account the action of a pseudodifferential operator on an oscillatory integral in an invariant manner, along the lines of [11,Section 4.3].
Remark 3.12. Let us point out that in this section we did not use anywhere the fact that M carries a Riemannian structure. If one replaces the Riemannian density (1.1) with an arbitrary positive density, all statements and results stay the same.

Invariant description of pseudodifferential operators acting on scalar functions
In order to prepare ourselves to address the issue of initial conditions for our transport equations in the case of the Dirac operator, we need to discuss first the more general question of invariant representation of a pseudodifferential operator. We devote a separate section to this, as we believe this matter to be of independent interest. Note that we treat the case of a scalar operator merely for the sake of presentational convenience: all the formulae and arguments in this subsection remain unchanged for matrix pseudodifferential operators acting on m-columns of scalar functions.
Global propagator for the massless Dirac operator when x lies in a geodesic neighbourhood of y and continued smoothly elsewhere in such a way that Im φ ≥ 0. Here γ is the (unique) shortest geodesic connecting y to x, ζ is the parallel transport of η along γ, dist is the geodesic distance and ǫ is a positive parameter.
Let P be a pseudodifferential operator of order p acting on scalar functions over a Riemannian d-manifold. The operator P can be written, modulo an integral operator with smooth kernel, in the form where φ is the time-independent Levi-Civita phase function, p ∈ S m ph (T ′ M), χ 0 is a cut-off localising integration to a neighbourhood of the diagonal and away from the zero section (see also (I)-(III) in Theorem 3.2) and Here the smooth branch of the complex root is chosen in such a way that w 0 (y; y, η) = [ρ(y)] −1 .
Remark 4.2. Note that (4.3) is, effectively, a special case of (3.10) with t = 0. Formula (4.3) provides an invariant representation of the pseudodifferential operator P . Furthermore, we call the homogeneous functions p p and p p−1 the g-principal and gsubprincipal symbol, respectively 4 .
The notions of principal and subprincipal symbols of a pseudodifferential operator are nowadays standard concepts in microlocal analysis. The former makes sense for operators acting either on scalar functions or on half-densities, whereas the latter is only defined for operators acting on half-densities. We refer the reader to [33] for further details. Note that the concept of subprincipal symbol was introduced by Duistermaat  It is easy to see that the concept of principal symbol P prin and that of g-principal symbol p p coincide. As far as the subprincipal symbol is concerned, the situation is more complicated, in that before drawing a comparison we need to turn our operator into an operator acting on half-densities.
Put P 1/2 := ρ 1/2 P ρ −1/2 (4.5) and let P sub be the subprincipal symbol of the operator (4.5) defined in accordance with [ Theorem 4.4 implies that, in particular, the two notions of subprincipal symbol coincide when the principal symbol does not depend on η, i.e. when P is a pseudodifferential operator of the type "multiplication by a scalar function plus an operator of order −1". Note that the identity operator, whose invariant representation was investigated in [17,Section 6], falls into this class.
Remark 4.5. A tedious, yet straightforward, calculation shows that the RHS of (4.6) is a scalar function on the cotangent bundle. In fact, the second and third summands on the RHS of (4.6) admit an invariant representation in terms of the Laplace-Beltrami operator associated with the neutral metric n on the cotangent bundle T * M, which, in local coordinates (x 1 , . . . , x d , ξ 1 , . . . , ξ d ), reads The adjective 'neutral' refers to the fact that the metric n has signature (d, d). It turns out that the neutral metric is an effective tool in the development of an invariant theory of pseudodifferential operators on Riemannian manifolds. As the analysis of this matter requires a lengthy discussion and would take us away from the core subject of our paper, we plan to address it in detail elsewhere. See also [46].
Proof of Theorem 4.4. Consider the pseudodifferential operator P and turn it into an operator on half-densities P 1/2 via (4.5). In what follows we work in an arbitrary coordinate system, the same for x and y.
Dropping the cut-off, the integral kernel of P 1/2 now reads Our phase function (4.1) admits the expansion 9) which implies that Substituting (4.9) and (4.10) into (4.8), we get Global propagator for the massless Dirac operator Excluding the x-dependence from the amplitude in (4.11) by acting with the operator we arrive at Computing the subprincipal symbol of (4.13) and using the fact that p p = P prin = (P 1/2 ) prin , we obtain (4.6). Note that the sign in front of the correction term i 2 (P prin ) y α ηα is opposite to the usual one, see, for example, [19,Eqn. (A.3)]. This is due to the fact that in this paper we use the right -as opposed to left -quantization.

Global propagator for the Dirac operator
In this section we will start the analysis of the global propagator for the Dirac operator, specialising Theorem 3.3 to the case A = W . We denote by W prin (y, η) := σ α (y) η α (5.1) the principal symbol of W and by its zero order part, see Definition 1.1. The principal symbol W prin (y, η) has eigenvalues h ± = ±h, where h is given by (4.2), compare with (3.2). This fact, which can be easily established by writing down (5.1) in local coordinates, shows that the Dirac operator is indeed elliptic.
Definition 5.1. We call positive (+), resp. negative (−), Levi-Civita phase function the infinitely smooth function ϕ ± ∈ C ∞ (R × M × T ′ M; C) defined by Global propagator for the massless Dirac operator for x in a geodesic neighbourhood of x ± (t; y, η) and continued smoothly elsewhere in such a way that Im ϕ ± ≥ 0. Here dist is the Riemannian geodesic distance, the path of integration γ ± is the shortest geodesic connecting x ± to x, ζ ± is the result of parallel transport of ξ ± (t; y, η) along γ ± and ǫ is a positive parameter.
The positive and negative Levi-Civita phase functions are related as ϕ − (t, x; y, η) = −ϕ + (t, x; y, −η). It is easy to see that the positive and negative Levi-Civita phase functions satisfy conditions (i), (ii) and (iv) from Theorem 3.3. Furthermore, [45,Corollary 2.4.5] implies that condition (iii) is also satisfied. Hence, Theorem 3.3 ensures that the integral kernel of U ± can be written as a single oscillatory integral where ϕ ± is the positive/negative Levi-Civita phase function.
Definition 5.3. We define the full symbol of the positive (resp. negative) propagator to be the scalar matrix-function a + (resp. a − ), obtained through the algorithm described in Section 3.3 with Levi-Civita phase functions. We define the subprincipal symbol of the positive (resp. negative) propagator to be the scalar matrix-function a + −1 (resp. a − −1 ) obtained the same way. As to the principal symbol, this object was defined earlier, see Definition 3.8. We stress that the mathematical objects contained in the above definition are uniquely and invariantly defined. They only depend on the phase functions which, in turn, originate from the geometry of M in a coordinate-free covariant manner, cf. Definition 5.1.
To the best of our knowledge, there is no accepted definition of full symbol or subprincipal symbol for a Fourier integral operator available in the literature to date. The geometric nature of our construction allows us to provide invariant definitions of full and subprincipal symbol of the Dirac propagator, analyse them, and give explicit formulae. This paper, alongside [17], aims to build towards an invariant theory for pseudodifferential and Fourier integral operators on manifolds.
Before moving on to computing the principal and subprincipal symbols of the positive (resp. negative) Dirac propagator, an important remark is in order. In addition to what was discussed in Section 3 for the general case, the construction of the Dirac propagator has to be consistent with the gauge transformation (1.5), (1.6). In particular, the action of the gauge transformation needs to be carefully accounted for by the construction process.
in the oscillatory integral (5.7). Note that this introduces an x-dependence which has to be handled by means of amplitude-to-symbol reduction (3.31).

Transport equations
By acting with the Dirac operator W on (5.7) in the variable x and dropping the cut-off, we obtain 10) for k ≥ 0. Note that the a ± −k , k ≥ −1, are positively homogeneous in momentum of degree −k.
Recalling that v ± are the normalised eigenvectors of W prin corresponding to the eigenvalues ±h, denote by P ± (y, η) := v ± (y, η) [v ± (y, η)] * (5.14) the spectral projections along the eigenspaces spanned by v ± . Of course, Let us label the transport equations with nonnegative integer numbers in increasing order, so that (5.11) is the zeroth transport equation, (5.12) is the first transport equation and so on. Direct inspection of (5.9) and (5.10) reveals that Global propagator for the massless Dirac operator • multiplication of the n-th transport equation by P ∓ (x ± , ξ ± ) on the left allows one to determine P ∓ (x ± , ξ ± )a ± −n (t; y, η), n ≥ 0, (5.18) algebraically; • multiplication of the (n + 1)-th transport equation by P ± (x ± , ξ ± ) on the left and the use of (5.18) allows one to determine P ± (x ± , ξ ± )a ± −n (t; y, η), n ≥ 0, upon solving a matrix ordinary differential equation in the variable t.

Pseudodifferential operators U ± (0)
This subsection is devoted to the examination of operators U ± (0). We need to examine these operators because, as explained in subsection 3.3, their full symbols determine the initial conditions a ± −k (0; y, η) for our transport equations. We have where θ(λ) := 1 for λ > 0, 0 for λ ≤ 0.
We see that the operators U ± (0) are self-adjoint pseudodifferential operators of order zero, orthogonal projections onto the positive/negative eigenspaces of the operator W . The operator Id − U + (0) − U − (0) is the orthogonal projection onto the nullspace of the operator W , hence The principal symbols of the operators U ± (0) read [U ± (0)] prin = P ± (y, η), (5.21) where P ± are the orthogonal projections onto the positive/negative eigenspaces of the principal symbol of the operator W , see (5.14). The analysis of the full symbol of U ± (0) is a delicate task which was investigated, to a certain extent and in a somewhat different setting, in [23]. In order to develop the ideas from [23] we have to address a number of issues.
• We are now dealing with scalar fields as opposed to half-densities.
• We are now making full use of Riemannian structure.
• We are now working in the special setting of a system of two equations in dimension three with trace-free principal symbol.

Global propagator for the massless Dirac operator
In order to calculate the subprincipal symbols of the pseudodifferential operators U ± (0) we will need the following auxiliary result. Then where K (resp. K) is the contorsion tensor of the Weitzenböck connection (see Appendix A) associated with the framing {e j } 3 j=1 (resp. { e j } 3 j=1 ), the star stands for the Hodge dual applied in the first and third indices, see formula (A.7), and σ α (y) is defined by (1.3).
Proof. The proof is provided in Appendix B.1.
The following theorem is the main result of this subsection.
Theorem 5.6. We have

25)
where T is the torsion tensor of the Weitzenböck connection (see Appendix A) associated with the framing {e j } 3 j=1 encoded within the Dirac operator W (see Definition 1.1) and the star stands for the Hodge dual applied in the second and third indices, see formula (A.6).
Proof. Let us fix a point y ∈ M and choose normal geodesic coordinates x centred at y such that e j α (y) = δ j α . Consider the (local) operator with constant coefficients where the s α are the standard Pauli matrices (1.2). Let us choose a smooth special unitary 2 × 2 matrix-function G such that compare with (5.24). It is easy to see that such a matrix-function G(x) exists and is defined uniquely modulo O( x 2 ).
Let us now compare the subprincipal symbols of the pseudodifferential operators θ(±W ) and θ(±G * W G), with G * W G understood as an operator acting in Euclidean space (constant metric tensor g αβ (x) = δ αβ ). It can be shown that at the origin we have Thus, the proof of the Theorem 5.6 has been reduced to the case when we are in Euclidean space and the operator W is given by formulae (5.24) and (5.26).
We have where Formulae (5.27) and (5.28) imply that Excluding the z-dependence from the amplitude Q ± by acting with the operator compare with (4.12), we arrive at Substituting (5.30) and (5.31) into (5.32) and setting x = 0, we get Global propagator for the massless Dirac operator Theorem 5.4 tells us that G x µ = i 2 * K µν s ν . Substituting this into (5.33), and using standard properties of Pauli matrices and (A.10), we get The above argument combined with (5.20) yields (5.25).
Observe that formula (5.25) implies where ζ ± (t; y, η) is the parallel transport of v ± (y, η) along x ± with respect to the spin connection, i.e.
Proof. It is known [44] [43,Subsection 3.4] that the principal symbols a ± 0 are independent of the choice of the phase function and read on matrix-functions on the cotangent bundle. In formula (6.5) the second term on the RHS is the result of switching to half-densities, see (4.5).
Introducing the shorthand q ± (t) := q ± (x ± (t; y, η), ξ ± (t; y, η)), the task at hand is to show that ζ ± (t; y, η) = e −i t 0 q ± (τ ) dτ v ± (x ± , ξ ± ). More explicitly, we need to show that where we premultiplied our expression by e i t 0 q ± (τ ) dτ for the sake of convenience. We shall prove (6.1) for a + 0 , which corresponds to the upper choice of signs in (6.6). The proof for a − 0 is analogous. Let us begin by computing To this end, let us choose geodesic normal coordinates centred at x + (t; y, η) = 0 and such that [ξ + (t; y, η)] α = δ 3α . Furthermore, up to a global rigid rotation of the framing, we can assume that e j α (0) = δ j α .

Explicit small time expansion of the symbol
Even though the presence of gauge degrees of freedom represents an additional challenge in the analysis of the propagator, one can put this freedom to use and exploit it to obtain a small time expansion for the propagator.
Our strategy goes as follows.
1. Compute the principal and subprincipal symbols of the positive (resp. negative) propagator for a conveniently chosen framing; 2. Using the gauge transformation (1.7), (1.6), switch to an arbitrary framing with the same orientation 5 ; 3. Express the final result in terms of geometric invariants.

Special framing
Let us fix an arbitrary point y ∈ M and let V j ∈ T y M, j = 1, 2, 3 be defined by V j := e j (y).
Definition 7.1 (Levi-Civita framing). Let U be a geodesic neighbourhood of y. For x ∈ U, let e loc j (x), j = 1, 2, 3, be the parallel transport of V j along the shortest geodesic connecting y to x. We define the Levi-Civita framing generated by {e j } 3 j=1 at y to be the equivalence class of framings coinciding with { e loc j } 3 j=1 in a neighbourhood of y.
With slight abuse of notation, in the following we will identify the Levi-Civita framing with one of its representatives, denoted by { e j } 3 j=1 . The choice of a particular representative does not affect our results.
Using the Levi-Civita framing is especially convenient due to the following property.
Lemma 7.2. In normal coordinates centred at y, the Levi-Civita framing admits the following expansion: where R is the Riemann curvature tensor 6 .
Proof. In normal geodesic coordinates centred at y, the unique geodesic connecting y to x can be written as where · E is the Euclidean norm, so that γ( x E ) = x. Assuming t and x E to be small and of the same order, let us perform an expansion in powers of t of e j . The parallel transport equation defining the framing { e j } 3 j=1 readṡ Since e j (0) = V j and Γ(0) = 0, at linear order in t we have˙ e j (γ(t)) = O(t), which implies and Formula (7.1) follows at once from (7.5) and the elementary identity The Riemann curvature tensor R has components R κ λµν defined in accordance with Corollary 7.3. In normal coordinates x centred at y, the Pauli matrices σ α (x) projected along the Levi-Civita framing (see (1.3)) satisfy Proof of Corollary 7.3. Formula (7.7) follows immediately from (7.1).
Corollary 7.4. Let W be the Dirac operator (1.4) corresponding to the choice of the Levi-Civita framing. Then, in normal coordinates centred at y, its zero order part W 0 (see formula (5.2)) admits the following expansion: Here Ric is the Ricci tensor, Ric αβ := R γ αγβ .
Proof. Formula (7.8) is obtained by expanding the RHS of (5.2) in powers of x in normal coordinates centred at y, substituting (7.6) and (7.7) in and performing a lengthy but straightforward calculation. It is a somewhat nontrivial fact that the coefficient of the linear term in (7.8) turns out to be trace-free.

Small time expansion of the principal symbols
The first step towards computing small time expansions for principal and subprincipal symbols of W is to obtain an expression for these objects in a neighbourhood of a given point y ∈ M for the choice of the Levi-Civita framing generated by our framing {e j } 3 j=1 at y. Observe that, as we are after a small time expansion of the symbols, it is enough to restrict our attention to a small open neighbourhood of y.
In the following, we will denote with a tilde objects associated with the Dirac operator W corresponding to the choice of the Levi-Civita framing.
Theorem 7.5. For the choice of the Levi-Civita framing, the positive and negative principal symbols are independent of t and read a ± 0 (t; y, η) = P ± (y, η). (7.9) Proof. In accordance with Theorem 6.1, the principal symbols are determined by the eigenvectors of W prin and their parallel transport with respect with the spin connection along the Hamiltonian trajectories. Hence, it suffices to show that ζ ± (t; y, η) = v ± (y, η). (7.10) Once this is achieved, (7.9) follows from the fact that W prin (y, η) = W prin (y, η).
In normal coordinates centred at y the parallel transport equation (6.2) reads We claim that Global propagator for the massless Dirac operator In fact, we have in view of Definition 7.1 and the properties of the Hamiltonian flows x ± , i.e. that x + ( · ; y, η) is geodesic and relation (5.3). By substituting (7.12) into (7.11) we arrive at (7.10).

Small time expansion of the subprincipal symbols
Let us now turn our attention to the subprincipal symbols a ± −1 . Unlike the principal symbols, the subprincipal symbols depend on the choice of phase functions. As here we are only interested in small time expansions and the injectivity radius Inj(M, g) is strictly positive, we can work, without loss of generality, in a neighbourhood of y with no conjugate points to y. The absence of conjugate points allows us to construct positive and negative propagators for small times by means of the algorithm described in subsection 3.3 for the choice of real-valued Levi-Civita phase functions In the remainder of this subsection we adopt the same coordinates for x and y and we choose normal geodesic coordinates centred at y. We remind the reader that, in such coordinates, and w ± (t, x; 0, η) = 1 + 1 12 Recall that the weight w is defined by (3.12). As explained in subsection 5.1, the subprincipal symbols are determined by the first and the second transport equations, (5.12) and (5.13). More precisely, if we are interested in expansions with remainder O(t 2 ), we need to determine (5.13) up to zeroth order in t and (5.12) up to first order in t.
To this end, we begin by observing that formulae (7.14) and (7.15), see also (3.30), imply that the differential evaluation operators S ± −2 and S ± −1 admit the following expansions in normal coordinates centred at y. Lemma 7.6. We have (a) Proof. (a) It is an immediate consequence of (7.13), (7.15) and (b) Substituting (7.15) into (3.30) with k = 1 and recalling that ϕ ± η x=x ± = 0, we get Formula (7.18) and the fact that In order to be able to compute the subprincipal symbols, we need to determine the initial condition a ± −1 | t=0 first. Proof. The subprincipal symbols are scalar functions, so it enough to establish (7.21) in one specific coordinate system. Let us choose normal coordinates centred at y = 0 such that e j α (0) = δ j α . We observe that the torsion of the Weitzenböck connection generated by the Levi-Civita framing at y vanishes at y, as a consequence of the fact that the first derivatives of the framing are zero, cf. A straightforward perturbation argument shows that Substituting (7.22) and (7.23) into (4.6) with P = U ± (0) and ǫ = 0 and using the fact that Christoffel symbols vanish at y, we arrive at (7.21).
We are now in a position to examine the first transport equation.
Lemma 7.8. The projection onto the negative (resp. positive) eigenspace of W prim of the subprincipal symbol of the positive (resp. negative) propagator is given by Ric αβ (y) η α P ± η β (y, η) Global propagator for the massless Dirac operator Proof. We will establish formula (7.24) by expanding the first transport equation (5.12) up to first order in t and then acting with P ∓ on the left. Recall that a ± −k is defined by (5.10).
Working in normal coordinates centred at y and using (7.14)-(7.15), we obtain (7.25) Furthermore, in view of Theorem 6.1 and Lemma 7.6(b), we have (7.26) Adding up (7.25) and (7.26) and projecting along P ∓ , we arrive at
Let us now move to the second transport equation.
Lemma 7.9. The projection onto the positive (resp. negative) eigenspace of W prin of the subprincipal symbol of the positive (resp. negative) propagator is given by Global propagator for the massless Dirac operator Proof. We will establish formula (7.35) by computing the second transport equation (5.13) up to zeroth order in t and then acting with P ± on the left. With account of Lemma 7.6, we have In carrying out the above calculations we used Theorem 7.5 and Lemma 7.7. Note that, when multiplying on the left by P ± , the terms containing a ± −2 disappear. Summing up (7.36), (7.37) and (7.38), and projecting along P ± , we obtain (7.40) Let us put (7.42) • B 2 : Differentiating (5.17) with respect to η β yields (7.43) Substituting (7.43) into B 2 in (7.41) we obtain Now, since P ± σ ρ (0)η ρ = P ± W prin (0, η) = ±h P ± and σ α σ ρ = − σ ρ σ α + 2 δ αρ Id, formula (7.45) implies Summing up (7.42), (7.44) and (7.46) we arrive at A straightforward calculation shows that Substituting the above expression into (7.47) and integrating in time with initial condition (7.21), we obtain (7.35).
The pieces of information from Lemma 7.8 and Lemma 7.9 can be combined to give the following result.
Theorem 7.10. For the choice of the Levi-Civita framing, the subprincipal symbols of the positive and negative propagators admit the following small time expansion: Proof. Summing up formulae (7.24) and (7.35), we obtain The substitution of (7.43) into the RHS of the above equation gives (7.48).

Invariant reformulation
In the previous subsections we derived the quite elegant and compact formulae (7.9) and (7.48), which were obtained under the assumption that the chosen framing is the Levi-Civita framing at y. Now the task at hand is to obtain similar formulae for the Dirac operator W corresponding to an arbitrary framing {e j } 3 j=1 . Given a framing {e j } 3 j=1 and a point y ∈ M, there exits a special unitary matrixfunction G, defined in a neighbourhood of y, such that {e j } 3 j=1 and the Levi-Civita framing { e j } 3 j=1 generated by {e j } 3 j=1 at y are related in accordance with e j α (x) = 1 2 tr(s j G * (x) s k G(x)) e k α (x), G(y) = Id, (7.49) cf. (5.22) and (1.7). The symbols a ± and a ± are related as cf. Section 5. Note that on the RHS of (7.50) the transformed symbol is acted upon by amplitude-to-symbol operators (3.31). The latter are needed because the gauge transformation G introduces an x-dependence in the amplitude, which has to be excluded. Working in normal coordinates centred at y, formula (7.50), combined with (7.13) and (7.9), implies Similarly, by means of (7.13) and Lemma 7.6, from (7.50) we get The last step towards expressing (7.51) and (7.52) invariantly is writing ∇G and ∇∇G in terms of geometric invariants. Theorem 5.4 tells us that The following theorem provides an expression for the second covariant derivatives of the gauge transformation.
Theorem 7.11. Let us fix a point y and let the special unitary matrix-function G be such that our framing {e j } 3 j=1 and the Levi-Civita framing { e j } 3 j=1 generated by {e j } 3 j=1 at y are related in accordance with (7.49) in a neighbourhood of y. Then we have where K is the contorsion tensor of the Weitzenböck connection (see Appendix A) associated with the framing {e j } 3 j=1 and the star stands for the Hodge dual applied in the first and third indices (see formula (A.7)).
Proof. The proof is given in Appendix B.2.
Remark 7.12. Note that, remarkably, the curvature of the Levi-Civita connection does not appear in the RHS of (7.54). Substituting (7.53) and (7.54) into (7.51) and (7.52) we arrive at the following result.

55)
Global propagator for the massless Dirac operator where * K denotes the Hodge dual in the first and third indices of the contorsion tensor of the Weitzenböck connection associated with the framing {e j } 3 j=1 .

An application: spectral asymptotics
In this section we will compute the third Weyl coefficient for the Dirac operator. In doing so we will use the same notation as in Section 1 -recall in particular formulae (1.17), (1.13) and the definition of the function µ.
where R is scalar curvature.
Proof. Let us fix a point y ∈ M and choose normal geodesic coordinates x centred at y. Let us also choose a Levi-Civita framing { e j } 3 j=1 , see Definition 7.1; here we make use of the fact that Weyl coefficients do not depend on the choice of framing.
We have where u ± is the Schwartz kernel of the propagator U ± and tr stands for the matrix trace. Note that at each point of the manifold the quantity tr u ± (t, y, y) is a distribution in the variable t and the construction presented in preceding sections allows us to write down this distribution explicitly, modulo a smooth function. Our task is to substitute (5.7) into the right-hand sides of (8.2) and (8.3) and expand the resulting quantities in powers of λ as λ → +∞. Thus, the problem reduces to the analysis of explicit integrals in four variables, η 1 , η 2 , η 3 and t, depending on the parameter λ . In what follows we drop the y in our intermediate calculations.
The construction presented in preceding sections tells us that the only singularity of the distribution tr u ± (t, y, y)μ(t) is at t = 0. Hence, in what follows, we can assume that the support ofμ is arbitrarily small. In particular, this allows us to use the real-valued (ǫ = 0) Levi-Civita phase functions ϕ ± . Theorems 7.5 and 7.10 imply that tr a ± 0 (t; η) = 1 ,

4)
Global propagator for the massless Dirac operator Formula [17, (B.11)] reads ϕ + (t, η) = − η t + O(t 4 ) , which, in view of (5.5), implies Using formulae (8.2)-(8.6) and arguing as in [17, Appendix B], we conclude that where S 2 = 4π is the surface area of the 2-sphere. But 1 2π R 2 r m e i(λ−r)tμ (t) dr dt = λ m , m = 0, 1, 2, . . . , so (8.7) can be rewritten as We see that the large (in modulus) eigenvalues of the Dirac operator are distributed approximately the same way as the eigenvalues of the operator √ −∆ , differing only in the third Weyl coefficient. which are related to the counting functions via the Mellin transform, as in [26,14]. See also [22].

Examples
In this section we present two explicit examples, which show how our constructions work in practice and which give us an opportunity to double-check our formulae. The specific choice of examples is motivated by the fact that the first, M = S 3 , is isotropic in momentum whereas the second, M = S 2 × S 1 , is anisotropic in momentum.

The case M = S 3
Let R 4 be Euclidean space equipped with Cartesian coordinates x α , α = 1, 2, 3, 4, and put Consider the 3-sphere 7 with orientation prescribed in accordance with [28, Appendix A], equipped with the standard round metric g and with the global framings {V ±,k } 3 k=1 defined as the restriction to S 3 of the vector fields in R 4 It is easy to check that the vector fields (9.1) are tangent to S 3 , so that they restrict to smooth vector fields on the 3-sphere. Note that (9.1) is an adaptation of [28, Eqn. (C.1)] to the case at hand. Let us introduce coordinates on S 3 with the north pole excised by stereographically projecting it onto the hyperplane tangent to the 3-sphere at the south pole. The stereographic map is given by It is easy to see that the coordinate system (u, v, w) has positive orientation. In stereographic coordinates the metric reads and our framings are given by 2V ±,2 = (uv ± 2w) ∂ ∂u (9.3) 7 We shifted the sphere so as to place the south pole at the origin.
A straightforward calculation shows that {V ±,k } 3 k=1 are positively oriented framings formed by (orthonormal) smooth Killing vector fields with respect to the metric g.
The framings {V ±,k } 3 k=1 define, via (1.4), two Dirac operators W ± related in accordance with is the SU(2) gauge transformation relating the two framings via (7.49) with e k = V +,k and e k = V −,k . Let us deal with W + first. On account of the symmetries of the 3-sphere, we will write formulae for principal and subprincipal symbols of the propagator of W + at the south pole (y = (0, 0, 0)) for the choice of momentum η = (0, 0, 1).
Global propagator for the massless Dirac operator Taking the Fourier transform of the RHS of (9.18) we get Let M = S 2 × S 1 be endowed with the metric g = g S 2 + dϕ 2 , where g S 2 is the round metric on the 2-sphere. Let y ∈ M be given. In this subsection we shall compute the small time expansion for the subprincipal symbols of the Dirac propagator W associated with a Levi-Civita framing at y. In this case, the result will not be isotropic in momentum η, because, unlike the previous example, (S 2 × S 1 , g) is not an Einstein manifold. Without loss of generality, we assume that y coincides with the north pole when projected onto S 2 . The exponential map exp y : T y M → M is realised explicitly by (u, v, w) → (θ = √ u 2 + v 2 , φ = arctan(v/u), ϕ = w), (9.21) where (θ, φ) are standard spherical coordinates on S 2 . Formula (9.21) defines geodesic normal coordinates (u, v, w) in a neighbourhood of y. In such coordinates, the metric g reads g(u, v, w) = 1 We will assume that normal coordinates are chosen so that the Levi-Civita framing satisfies e j α (y) = δ j α . In this case, the Hamiltonian flows generated by the eigenvalues of W prin read, simply, z ± (t; 0, η) = ±t η η , ξ ± (t; 0, η) = η.
The Ricci curvature of g in normal coordinates (u, v, w) is given by Hence, Theorem 7.10 tells us that a ± −1 (t; y, η) = ∓it 1 12 η P ± (y, η) − 1 8 η 2 (η 1 σ 1 (y) + η 2 σ 2 (y)) + O(t 2 ) = ∓ it 24 η 2 η Id +(−3 ± 1) η α σ α (y) + 3 s 3 η β e 3 β (y) + O(t 2 ), (9.24) where the s j and the σ α are defined by formulae (1.2) and (1.3) respectively, and e 3 is the vector field ∂/∂ϕ (unit vector field along the positive direction of the circle S 1 ). Let us stress once again that, even though the intermediate steps depend on the choice of coordinates, the final result (9.24) is a scalar matrix-function, thus independent of the choice of coordinates. The only assumption involved in the derivation of formula (9.24) is that we used a particular Levi-Civita framing at the point y, one which respects the product structure of the manifold. The presence of the vector field e 3 in formula (9.24) is a manifestation of anisotropy.

Acknowledgements
We are grateful to Yiannis Petridis for providing useful references. Furthermore, we would like to thank an anonymous referee for insightful comments, in particular, for suggesting the argument in Remark 3.9. and the curvature tensor vanishes identically. The Weitzenböck connection coefficients and the Christoffel symbols are related via the identity is called contorsion of ∇ W . Note that the torsion tensor is antisymmetric in the second and third indices, T α βγ = −T α γβ , whereas the contorsion tensor is antisymmetric in the first and third ones, K αβγ = −K γβα (the first index was lowered using the metric). Torsion and contorsion can be expressed one in terms of the other and capture the geometric information encoded within the framing.
In dimension three antisymmetric tensors of order two are equivalent to vectors. Therefore, we define where F kα := [A k ] x α (0). Now, differentiating (5.22) with respect to x and evaluating the result at 0, we obtain Substitution of (B.7) into (B.2) gives (5.23).

B.2 Proof of Theorem 7.11
Recall that according to formula (7.49) we have G(y) = Id. In the following we work in a sufficiently small neighbourhood U of y and we choose normal coordinates centred at y = 0 such that e j α (0) = e j α (0) = δ j α . Since G ∈ C ∞ (M; SU(2)) and G(0) = Id, there exist smooth real-valued functions A k , k = 1, 2, 3, such that v k (0) = 0 and in a neighbourhood of y = 0. Differentiating (B.8) twice with respect to x and evaluating the result at zero we obtain =is k H kαβ − δ jk Id F jα F kβ . The RHS of (B.19) is symmetric in α and γ, whereas ε α γ ρ is antisymmetric in the same indices, so (B.18) follows.