Modular forms in the spectral action of Bianchi IX gravitational instantons

We prove a modularity property for the heat kernel and the Seeley-deWitt coefficients of the heat kernel expansion for the Dirac-Laplacian on the Bianchi IX gravitational instantons. We prove, via an isospectrality result for the Dirac operators, that each term in the expansion is a vector-valued modular form, with an associated ordinary (meromorphic) modular form of weight 2. We discuss explicit examples related to well known modular forms. Our results show the existence of arithmetic structures in Euclidean gravity models based on the spectral action functional.


Introduction
In this paper we prove that each Seeley-deWitt coefficient in the full heat kernel expansion for the Dirac-Laplacian on the Bianchi IX gravitational instantons is a vector-valued modular form, and we illustrate how they can be related to well known modular forms. The general form of the instantons involves a conformal factor that is a function of the cosmic time µ [4,25,39]. Thus, in section 2, we start by presenting an explicit formula for the Dirac operatorD of the Bianchi IX metrics with general cosmic expansion factors JHEP01(2019)234 w 1 (µ), w 2 (µ), w 3 (µ) and with a general conformal factor F (µ). We then prove a rationality result for the recursive structure of the Seeley-deWitt coefficients associated withD 2 , which illuminates the arithmetic nature of the coefficients, and stimulates the search for more explicit arithmetic structures in the expansion. The rationality result generalizes the analogous rationality results for Robertson-Walker metrics [8,19], and for the Bianchi IX case without a conformal factor [16]. We also recall the explicit Babich-Korotkin parameterization [4] of the Bianchi IX gravitational instantons in terms of theta functions with characteristics, which we shall use throughout the paper.
In section 3 we provide the necessary background material from the theory of modular forms that we will need for presenting our main results. Moreover, in order to prepare the ground for discovering the modular properties of the Seeley-deWitt coefficients of the Bianchi IX gravitational instantons, we identify the transformation properties of the anisotropy coefficients w 1 (µ), w 2 (µ), w 3 (µ) and the conformal factor F (µ) of the instantons. Since these functions are given in terms of theta functions, it makes sense to analytically continue the cosmic time variable µ to be a point in the right half-plane {z ∈ C : ℜ(z) > 0} so that the point iµ in the upper half-plane H = {z ∈ C : ℑ(z) > 0} is acted upon by the modular group PSL 2 (Z) = SL 2 (Z)/{±1}, by linear fractional transformations. It suffices to check the modular property with respect to the two generators of the modular group, commonly denoted by T and S, whose corresponding actions on iµ ∈ H are iµ → iµ + 1 and iµ → −1/(iµ), respectively. We comment briefly in section 1.1 on the role of this complexified time variable in physical models.
In section 4 we first derive the modular properties of the first two Seeley-deWitt coefficients, denoted byα 0 andα 2 , by performing explicit calculations. A general conceptual understanding of these modular properties, observed initially in the first two coefficients, is then achieved by showing that a particular modular action on the parameters of the Bianchi IX gravitational instantons gives rise to isospectral Dirac operators. We show that the isospectrality of the Dirac operators is indeed responsible for the modular properties of all the Seeley-deWitt coefficients in the full heat kernel expansion. In section 5 we prove that, for rational values of the parameters in the two-parameter family of Bianchi IX solutions, the Seeley-deWitt coefficients are vector-valued modular forms. We discuss the associated (meromorphic) modular form obtained by averaging over the orbit of the parameters, and we analyze its zeros and poles. Then, we illustrate in explicit examples how the resulting modular forms can be computed explicitly in terms of well known modular forms, namely the modular discriminant and Eisenstein series.
Finally, in section 6 we summarize our conclusions about the appearance of explicit arithmetic structures in spectral models of gravity and our approach to a conceptual understanding of the modularity of the Seeley-deWitt coefficients in the full expansion of the spectral action.
The modular properties presented in this article were initially identified through a series of lengthy explicit calculations of the first few Seeley-deWitt coefficients. This provided a clear picture of the overall pattern, on the basis of which we then derived the general proof of modularity presented in this paper. While we have not included here the explicit lengthy computations of the first few terms of the expansion, they remain available to the interested reader in the arXiv version of this article.

Modularity in physics
The focus of this paper is primarily on deriving an explicit form of the Seeley-deWitt coefficients for the heat kernel expansion of the Dirac Laplacian on Bianchi IX gravitational instantons and identifying the modular properties of these expressions in terms of vectorvalued modular forms. We will discuss briefly in the Conclusions section the relevance of these Seeley-deWitt coefficients in the asymptotic expansion of the spectral action, which is an action functional for a model of modified (Euclidean) gravity, where the terms of the expansion represent higher derivative corrections to the usual Einstein-Hilbert gravity with cosmological term. A detailed treatment of models for gravity and gravity coupled to matter based on the spectral action functional can be found in [12,32,41].
We comment here briefly on the relevance of modularity properties in the context of various kinds of physical models. Some of the best known occurrences of modularity in physics include elliptic functions, modular forms and quasi-modular forms that arise as characters and partition functions associated to vertex operator algebras; generating functions in 4D gauge theory; modular forms occurring in parametric Feynman integrals in perturbative quantum field theory; various moonshine phenomena in String Theory; generating functions associated to mirror symmetry, Gromov-Witten invariants, and Eynard-Orantin recursion formulas; quantum degeneracies of black holes as Fourier coefficients of mock Jacobi forms.
The modular phenomena in physics most relevant for the topic discussed in this paper are those that arise in special solutions of the Einstein equations. A general introduction to modularity in the Einstein equations is given in Part II of [26]. A review including recent results on (quasi)modularity in the Einstein equations and in gravitational instantons is given in [36]: in this overview the relevance of modularity to String Theory is explained through the role of certain classes of gravitational instantons in compactifications of superstring models, or through dynamics of moduli of Supergravity and String Theory as sigma models with target a gravitational instanton. Physical applications of the modularity properties of gravitational instantons discussed in [36] include monopole scattering, Ricci flows, nonperturbative instanton corrections to string moduli spaces, and perturbative expansions.
As we discuss in detail in the following sections, the modularity that we analyze in the coefficients of the heat kernel expansion on Bianchi IX gravitational instantons corresponds to modular transformations of the time coordinate. The fact that this can be viewed as a complexified coordinate reflects the fact that the Euclidean Bianchi IX metrics can be seen as Wick rotations of Lorentzian signature metrics. The use of a complexified time coordinate in Bianchi IX metrics, which appears as a mathematical artifact, is in fact used for example in the context of tunneling geometries and caustics of classical solutions (see for instance [6]), or when considering complex-time contours in a minisuperspace approximation for convergence of Hartle-Hawking path integrals (see [35] for an overview), in topological effects on solutions of the Yang-Mills-Einstein equations, where the complexified time variable has a finite temperature interpretation [11]. A discussion of complexified time for gravitational instanton geometries in a different setting is also given in [30,31].
The physical context where we expect the modularity properties analyzed in this paper to play a role is mainly the setting of a Hartle-Hawking type quantum cosmology, based JHEP01(2019)234 on the spectral action functional for gravity rather than the Einstein-Hilbert action. A general setup for such a model is introduced in section 9 of [32]. In this setting the complexification of the time variable plays a role analogous to the original Hartle-Hawking setting (see [23,24,35]). The role of modularity in this setting can be investigated along similar lines as in [36], though adapted to a setting where the action functional for gravity is given by the spectral action. We will return to a more thorough discussion of these physical aspects in forthcoming work.

Bianchi IX gravitational instantons
This section covers necessary background material. We first present the explicit form of the Dirac operator on the Bianchi IX metrics with a time dependent conformal factor. We then extend to the case with non-trivial conformal factor the rationality result of [16] for the recursive form of the heat kernel Seeley-deWitt coefficients, originally obtained for Bianchi IX metrics with a trivial conformal factor. We also recall the Babich-Korotkin parameterization [4] of the Bianchi IX gravitational instantons in terms of theta functions, which we will be using throughout the paper.
The general form of the triaxial Euclidean Bianchi IX metrics that we consider in article is given by where the cosmic expansion factors (anisotropy coefficients) w 1 , w 2 , w 3 are functions of the cosmic time µ and the conformal factor F is also a function of µ. The σ i are left-invariant 1-forms on SU(2)-orbits satisfying The Bianchi IX gravitational instantons were discovered by first finding solutions for the expansion factors w 1 , w 2 , w 3 that make the Weyl tensor of the metric self-dual. Exploiting the fact that the Weyl tensor is conformally invariant, it was then discovered that a suitable choice of the conformal factor F can be used to obtain an Einstein metric [4,25,39].

The Dirac operator
Given a spin bundle S and a spin connection ∇ S on a Riemannian manifold M , the Dirac operator D acting on the smooth sections C ∞ (S) of S is, by definition, the operator Here, the second arrow # is essentially the musical isomorphism identifying the cotangent and tangent bundles using the Riemannian metric, and the third arrow c is obtained by considering the Clifford action of the tangent bundle T M on S. The spin connection ∇ S is a connection on the Spin(m)-principal bundle associated with S, which is obtained by

JHEP01(2019)234
lifting the Levi-Civita connection ∇ seen as a connection on the SO(m)-principal bundle of T M [29]. The calculation of the Levi-Civita connection and the Dirac operator D of the Bianchi IX metric with a trivial conformal factor, was performed in [16]. The explicit local formula for D, when the 3-dimensional sphere S 3 is parametrized by with the parameter ranges 0 ≤ η ≤ π, 0 ≤ φ < 2π, 0 ≤ ψ < 4π, is written as The gamma matrices γ 0 , γ 1 , γ 2 , γ 3 are explicitly given by Since the general form of the Bianchi IX gravitational instantons given by (2.1) involves a time dependent conformal factor F (µ), after performing calculations for the Levi-Civita connection and the Dirac operatorD of the metric (2.1), we find that: where D is given by (2.3). The Dirac-LaplacianD 2 is a Laplace-type operator and identification of the arithmetic structures hidden in the Seeley-deWitt Coefficients a 2n (D 2 ) associated withD 2 is the main focus of the present paper.

The heat kernel of the Dirac-Laplacian
We consider the Dirac-LaplacianD 2 for the 4-dimensional geometry (2.1). The Seeley-deWitt coefficients a 2n (D 2 ) appear in the small time heat kernel expansion of this operator.

JHEP01(2019)234
Namely, as t → 0 + , the trace of the heat operator exp(−tD 2 ) admits an asymptotic expansion of the from which means that for small t > 0, the trace goes to infinity in a controlled manner: for any non-negative integer N , we have As we shall explain in more detail, there are densities a 2n (x,D 2 ) d 4 x on the manifold M , which can be obtained in theory from the Riemann curvature tensor and its contractions and covariant derivatives [12,22], such that Intuitively, this is related to the fact that, long before reaching equilibrium, the diffusion of heat on a curved manifold depends heavily on the local curvature. Note that this type of result holds in any dimension. However, we write them for dimension 4 because of the focus of the present paper on the 4-dimensional Bianchi IX spacetimes. We now explain how the asymptotic expansion (2.5) can be derived by making use of parametric pseudodifferential calculus. This is explained in full detail in [22]. The starting point is the use of the Cauchy integral formula to write where the integration is over a contour γ in the complex plane that goes clockwise around the non-negative real numbers, where the eigenvalues ofD 2 are located. SinceD 2 is an elliptic differential operator, the operator (D 2 − λ) −1 in the integrand can be well approximated by pseudodifferential operators. The idea is to pass to Fourier modes and to represent the operatorD 2 by its pseudodifferential symbol σD 2 (x, ξ). For each coordinate x on the manifold and frequency ξ, the symbol is a 4 × 4 matrix whose entries are complex numbers. The explicit formula, for any spinor s, is whereŝ is the component-wise Fourier transform of s. SinceD 2 is a differential operator of order 2, we have σD 2 (x, ξ) = p 2 (x, ξ) + p 1 (x, ξ) + p 0 (x, ξ), where each p k is homogeneous of order k in ξ. Therefore the approximation of (D 2 − λ) −1 reduces to finding pseudodifferential symbols r j (x, ξ, λ) of order −2 − j, j = 0, 1, 2, . . . , such that the operator R λ with symbol σ R λ = ∞ j=0 r j is the inverse ofD 2 − λ in the algebra of pseudodifferential operators, modulo infinitely smoothing operators. In fact one finds recursively that

JHEP01(2019)234
and for any n ≥ 1, with the summations over all 4-tuples of non-negative integers α, and integers 0 ≤ j < n and 0 ≤ k ≤ 2 with |α| + j + 2 − k = n. More importantly, the densities that integrate to the Seeley-deWitt coefficients as in (2.6), can be expressed by the formula The integrals in the latter can be calculated in theory using complex analysis and ordinary integration methods. However, due to the very rapid growth in the length of the expressions for the r n (x, ξ, λ), even with computer assistance one can only calculate the first few terms. Moreover, it is nearly impossible to prove conceptual results about the densities in the full expansion directly from the recursion formula above.

A rationality result for the Seeley-deWitt coefficients
A rationality phenomenon for the coefficients of the heat kernel expansion was first observed computationally in [8] in the case of Robertson-Walker metrics, and proved rigorously in [19]. A rigorous rationality result was then proved in [16] for the Bianchi IX metrics of the form (2.2) (without a conformal factor). The proof for the Bianchi IX case required the devising of a new method in [16] for calculating the Seeley-deWitt coefficients, which is based on considering auxiliary flat tori and making use of the Wodzicki residue and its properties [43,44]. In this subsection we show that a similar rationality result holds for the general form of the Bianchi IX gravitational instantons given by (2.1). This rationality property is a consequence of the symmetries of the Robertson-Walker and the Bianchi IX metrics, which have the effect of cancelling out (additively) the terms with irrational coefficients in the calculation of the Seeley-deWitt coefficients. These rationality results are also a strong indication of existence of richer arithmetically representable structures in the coefficients. Indeed, the same rationality property for the Robertson-Walker metrics was shown in [21] to be part of a deeper arithmetic structure formulated as a realization of the coefficients of the heat kernel expansion as periods of mixed Tate motives. More importantly, it was shown in [20] that the coefficients of the heat kernel expansion in the Robertson-Walker case can be written explicitly in terms of Bell polynomials and Brownian bridge integrals, extending the results of [8]. A motivic interpretation of the coefficients in the Bianchi IX case is presented in [17]. In order to state the rationality result in the case of Bianchi IX metrics with conformal factor, it is convenient to introduce a notationα 2n (µ) for the time dependent terms. These are obtained by integrating each density coefficient a 2n (x,D 2 ) of the Seeley-deWitt coefficient given by (2.6), where x = (µ, y) with y ∈ S 3 , over the 3-dimensional sphere: (2.10)

JHEP01(2019)234
As we shall mention in the proof of Proposition 2.1, due to the spatial symmetries of the Bianchi IX metric, a 2n (µ, y),D 2 is independent of y ∈ S 3 and the above integration just gives a multiplicative volume factor in the calculation. We now state the rationality result.
Proof. The argument is similar to the proof of Theorem 5.1 in [16], so we only outline its main steps briefly. One can exploit the SU(2)-invariance of the 1-forms σ i appearing in the metric (2.1) to show that functions of the type (2.9), whose integrals give the Seeley-deWitt coefficientsã 2n , have no spatial dependence when computed with respect to this metric. Then, one uses the Wodzicki residue method as in Theorem 4.1 of [16] for the computation of the heat kernel coefficients, which gives Here ∆ −1 denotes the parametrix of the elliptic differential operator σ ∆ −1 ,−2n−2 is the homogeneous component of its symbol of order −2n−2, and ∆ T 2n−2 is the flat Laplacian on an auxiliary (2n − 2)-dimensional torus T 2n−2 = (R/Z) 2n−2 . Considering the explicit form of the pseudodifferential symbol ofD 2 , which can be computed from the symbol for D 2 of appendix A of [16] and the relation between D 2 andD 2 using (2.4), it can be checked inductively that where tr(P 2n (ζ)) depends polynomially on the variables ζ = (ζ 1 , . . . , ζ 2n+2 ) ∈ S 2n+1 with coefficients in the algebra over Q generated by trigonometric functions of the spatial coordinates and the functions w k , for k ∈ {1, 2, 3} and ℓ ∈ {0, . . . , 2n}. Using the fact that the termsα 2n have no dependence on the spatial S 3 -coordinates, and the explicit form in terms of Gamma functions of the integration on S 2n+1 of monomials ζ r 1 1 · · · ζ r 2n+2 2n+2 , one arrives at (2.11).
It is indeed surprising that only rational numbers occur in the coefficients of the terms α 2n , expressed as a function of w 1 , w 2 , w 3 , F , and their derivatives. Usually such rationality results point to the existence of an underlying arithmetic structure. In this paper we show that for the Bianchi IX gravitational instantons the arithmetic structure is related to the occurrence of modular forms.

Theta function parameterization of Bianchi IX gravitational instantons
The Bianchi IX gravitational instantons are metrics of the form (2.1) that are self-dual and Einstein, that is, they satisfy the self-duality of the Weyl tensor and the proportionality between the Ricci tensor and the metric.
An especially interesting property of these metrics is the fact that the self-duality and Einstein equations reduce to a classical system of ordinary differential equations with singularities, the Painlevé VI equation [4,25,39]. In turn, solving these equations in terms of elliptic theta functions [4,25,39], and using the formula for the τ -function of the Schlesinger equation [27] one obtains explicit parameterizations of the Bianchi IX gravitational instantons in terms of theta functions [4].
The resulting parameterization of [4] of the Bianchi IX gravitational instantons provides two classes of solutions: a two-parameter family of solutions that correspond to the case with non-vanishing cosmological constant, and a one-parameter family which gives the vanishing cosmological constant case.

Two-parameter family of gravitational instantons
We follow the same notation as in [4] for the theta functions with characteristics. Namely, for p, q, z ∈ C, iµ ∈ H, we let (2.12) We also write where Θ(z|τ ) is the Jacobi theta function defined by The two-parameter family of solutions of [4] (Bianchi IX gravitational instantons with non-vanishing cosmological constant Λ) with parameters p, q ∈ C is obtained by substitut-

One-parameter family of gravitational instantons
The one-parameter family of [4] (Bianchi IX gravitational instantons with vanishing cosmological constant) with parameter q 0 ∈ R is obtain by substituting in the metric (2.1) the functions where C is an arbitrary positive constant.

Asymptotics and singularities
It is shown in [30] that the Bianchi IX solutions of the gravitational instanton equations approximate for large µ gravitational instantons of Eguchi-Hanson type with w 1 = w 2 = w 3 , [13]. Moreover, it was already observed in [25] that these manifolds have singularities for certain special values of µ. In terms of the explicit parameterization of [4] the singularities can occur where at least one of the functions ϑ[p, q], ∂ q ϑ[p, q], ∂ q ϑ[p, q + 1 2 ], ∂ q ϑ[p + 1 2 , q] and ∂ q ϑ[p + 1 2 , q + 1 2 ] vanish. Since we are looking at the properties of the Seeley-deWitt coefficientsα 2n (µ) given by (2.10), our arguments will apply as long as µ is away from the singularities of the Bianchi IX gravitational instantons and we work with an interval of values around µ which also does not contain any singularities. We will discuss this more in detail in section 4. We will see that the modular forms we obtain by working on such a domain then extend meromorphically.

Modular forms and modular action on the Bianchi IX parameters
In this section we provide some background material from the theory of modular forms. We will also investigate the basic properties of the functions w 1 , w 2 , w 3 , F given by (2.14) and (2.15) under the action of the modular group PSL 2 (Z) = SL 2 (Z)/{±1} on any iµ in the upper half-plane H, by linear fractional transformations. Modular forms appear in many areas of mathematical physics, see [45]. One of their main advantages of modular forms is that, once the modularity property is detected, one can use algorithms to write them explicitly in terms of well known modular forms whose properties are well understood [38].

JHEP01(2019)234
3.1 Modular forms, Eisenstein series G 2k and the modular discriminant ∆ This subsection is dedicated to providing from [37] the basic definitions, properties, and explicit examples of modular forms that we will need in the sequel to present our results.
For an integer k ∈ Z, a function f : H → C is a modular function of weight 2k if it is meromorphic and satisfies the condition In order to illuminate the geometric nature of this definition, we point out that, since one has d(g · z)/dz = (cz + d) −2 , and one can equivalently write (3.1) as This means that the equation is representing an invariant differential form of weight k.
Since the matrices , in order to show that a meromorphic function f : H → C is modular of weight 2k, it suffices to check that A modular function is a modular form if it is holomorphic on H and holomorphic at infinity. The latter needs an explanation. The condition f (z) = f (z + 1) in (3.2) is clearly showing that a modular function is defined by its values on the strip Therefore, using the transformation Q = e 2πiz = e −2πy e 2πix , any modular function f can be identified with a meromorphic functionf : . It is customary to use the small letter q for e 2πiz . In this context, however, in order to avoid confusion with the parameter q used in the two parameter family of gravitational instantons given by (2.14), we use the capital letter Q for e 2πiz . Note that Q = e 2πiz = e −2πy e 2πix → 0 as y → ∞. Therefore a modular function f is said to be holomorphic at infinity if its associated functionf is holomorphic at Q = 0. This makes it clear how to define meromorhphicity of f at infinity as well. To summarize, one can say that a modular form of weight 2k is a JHEP01(2019)234 function f : H → C defined on the upper half-plane by a convergent series of the following form with coefficients a n ∈ C, We denote the linear space of modular forms of weight 2k by M 2k . Note that the cusp forms of weight 2k form a linear subspace, denoted by M 0 2k , of M 2k . In fact, M 0 2k is the kernel of the linear functional on M 2k that sends any f as in (3.3) to a 0 = f (∞). Therefore, in general, M 2k = M 0 2k or M 2k is the direct sum of M 0 2k and a 1-dimensional linear space. We will indicate shortly that the Eisenstein series (3.4) and the modular discriminant (3.5) are the explicit tools that one can use to provide a generator for the remaining 1-dimensional space in M 2k and to map onto M 0 2k , respectively. The Eisenstein series G 2k , which is defined for any integer k ≥ 2, is a modular form of weight 2k and it is given by The modular discriminant is a cusp form of weight 12 and it is given by or, equivalently by Now that we have explained rigorously how the holomorphic and meromorphic behavior of the modular functions at infinity are described, we can review a fundamental formula, called the valence formula, which we use crucially in this work. The formula is quite important for understanding the structure of the modular forms and thereby finding algorithms for identifying them. Before stating the formula, note that for any meromorphic function f defined on the upper half-plane H and P ∈ H, we denote the order of zero of f at P by ν P (f ). Therefore ν P (f ) is 0 if f is neither 0 nor has a pole at P , it is the positive integer equal to the order of 0 if f (P ) = 0, and it is the negative integer equal to minus the order of the pole if f has a pole at P . For a modular function f : H → C of weight 2k, the valence formula gives an important relation between the order of its zeros (and poles) and its weight. Due to the special orbifold nature of the points i and ρ = e 2πi/3 on the modular curve, the zeros occurring at these two points are singled out from the remaining zeros (and poles) in the fundamental domain H/PSL 2 (Z) and counted with the appropriate JHEP01(2019)234 fractional multiplicity. The resulting valence formula for a modular function of weight 2k is written as Using this formula, one can readily see that there are no modular form of weight 2k when k is a negative integer and when k = 1. Further arguments based on the above facts yield the following facts. First, for k = 0, 2, 3, 4, 5, the space M 2k of modular form of weight 2k is 1-dimensional with the basis 1 for k = 0 and G 2k , for k = 2, 3, 4, 5. Second, for k ≥ 2, we have M 2k = M 0 2k ⊕ CG 2k . Third, multiplication by the modular discriminant ∆ defines an isomorphism from M 2k to M 0 2k+12 . Finally, the space M 2k is generated linearly by the monomials G α 4 G β 6 where α and β are non-negative integers such that 4α + 6β = 2k. A practical advantage that makes modular forms very useful objects is that there are explicit algorithms for writing them in terms of monomials of the Eisenstein series G 4 and G 6 .
We end this subsection by providing the definition of a more general type of modular forms called vector-valued modular forms, which have been the subject of intensive research because of their appearance and applications in a number of fields [5,14,28,34]. One of the main results of the present paper is that each Seeley-deWitt coefficient in the expansion of the spectral action for the Bianchi IX gravitational instantons is a vector-valued modular function of weight 2. For the general definition one has to consider a representation π 0 : SL 2 (Z) → GL(V ) on a finite dimensional vector space V . A vector-valued modular function of weight 2k with respect to π 0 is a meromorphic function f : H → V such that The function is a vector-valued modular form of weight 2k if it is holomorphic everywhere including at infinity. In section 5, we will show that the Seeley-deWitt coefficients of the Bianchi IX gravitational instantons are vector-valued modular forms with respect to an explicit representation. Note that it is common to use the terminology meromorphic modular form for modular functions. In this paper, we clearly point to any existing poles in our functions to avoid any confusion.

Modular action on the Bianchi IX parameters
In order to investigate the modularity properties of the heat kernel Seeley-deWitt coefficientsα 2n given by (2.10), we first need to compute explicitly how the anisotropic scaling factors w 1 , w 2 , w 3 , and the conformal factor F of the Bianchi IX gravitational instantons behave under the action of the group SL 2 (Z) acting on the variable iµ in the upper halfplane H by linear fractional transformations.

The two-parameter family
In the case of the two-parameter family (2.14), the anisotropic scaling factors w k [p, q](iµ), k = 1, 2, 3, and the conformal factor F [p, q](iµ) of the Bianchi IX gravitational instantons

JHEP01(2019)234
transform in the following way under the action of the generators T and S of the modular group. We recall that T (iµ) = iµ + 1 and S(iµ) = i/µ (we are assuming iµ ∈ H).
Theorem 3.1. The functions w j [p, q] and F of (2.14) transform under T as with their derivatives with respect to µ satisfying the same relations. The transformation under S of the w k , F and their derivatives is given by

JHEP01(2019)234
All these identities follow directly from the explicit expressions (2.12), (2.13), the identity and the Poisson summation formula, which gives We recall that the Poisson summation formula asserts that the sum on the integral lattice of the values of any rapidly decaying function is equal to the sum of the values of its Fourier transform on the lattice. For example, in the case of the last identity above one has The other cases are similar. We also see that the functions ϑ 2 , ϑ 3 , ϑ 4 , of (2.13) satisfy the transformation rules

JHEP01(2019)234
We then obtain The remaining cases are analogous and the relations for the derivatives are obtained directly by differentiating with respect to µ.

The one-parameter family
The argument for the one parameter case is analogous. The resulting transformation relations for the anisotropic scaling factors w k , k = 1, 2, 3, the conformal factor F , and their derivatives are as follows. (3.12) The derivatives with respect to µ satisfy the same transformation relation. The transformation under S is given by with derivatives transforming as (3.14)

Isospectral Dirac operators and modularity
Using the transformation properties of the functions w 1 , w 2 , w 3 , F under the generators S and T of SL 2 (Z), and explicit computation of the first two coefficientsα 0 andα 2 of the

JHEP01(2019)234
heat kernel expansion, we first identify the modularity property that theα 2n are expected to satisfy. We then prove it by a comparison of the Dirac operators and an isospectrality argument, which again relies directly on the transformation properties of the w 1 , w 2 , w 3 , F , which were discussed in section 3.2.
4.1 Modularity of the volume termα 0 and the scalar curvature termα 2 We compute here explicitly the transformation of the first two heat kernel Seeley-deWitt coefficientsα 0 andα 2 of the Bianchi IX gravitational instantons under the action of the generators T and S of the modular group.
Proof. These coefficients can be computed directly using either the parametrix method explained in section 2.1.1 or the method based on the Wodzicki residue which we elaborated on in section 2.1.2. The result is well known: theα 0 term is the volume form and theα 2 term is the Einstein-Hilbert action.
Theα 4 coefficient, which contains the Gauss-Bonnet and Weyl curvature terms, can also be computed explicitly, but the resulting expression is lengthier, so we do not report it here. Proposition 4.1. In the case of the two-parameter family of Bianchi IX gravitational instantons, the coefficientsα 0 andα 2 of (4.1) and (4.2) transform under the generators T (iµ) = iµ + 1 and S(iµ) = i/µ of the modular group as for µ a complex number in the right half-plane ℜ(µ) > 0. Similarly, in the one-parameter case the transformation under T and S is given bỹ

JHEP01(2019)234
Proof. The relations (3.8) and (3.12) are satisfied by the derivatives w (n) k and F (n) with respect to µ of arbitrary order n ≥ 0. Explicit transformation relations under S can be computed for the derivatives w (n) k and F (n) directly by differentiation. For the second derivatives we have, for the two-parameter and the one-parameter cases, respectively and The other cases are obtained in a similar way.
A lengthier computation confirms that, in the two-parameter case, the termα 4 also satisfies the same transformation relations asα 0 andα 2 , while in the one-parameter case it satisfiesα 4

Isospectrality of the Dirac operators under the modular action
The result of Proposition 4.1 suggests the modularity property that we should expect to hold for each termα 2n in the full expansion of the heat kernel. In principle, one could analyze inductively the Seeley-deWitt coefficientsα 2n , using (2.9), (2.7) and (2.8), but such an approach would be computationally cumbersome and less transparent in meaning. We use a different method, based on comparing directly the Dirac-LaplaciansD 2 [p, q], D 2 [p, q +p+ 1 2 ] andD 2 [−q, p], for the two-parameter case, and the Dirac- with √ h the volume element on ∂M , satisfies γ 5 = χγ(ν), with χ 2 = 1 and χγ 5 + γ 5 χ = 0 and determines local elliptic boundary conditions B ± = 1 2 (1 ∓ χ), see [9,18,42]. The Dirac operatorD with a boundary condition B as above has discrete real spectrum {λ n }, unbounded above and below, and an orthonormal basis of eigenspinors u n = u n (µ, η, φ, ψ), with u n ∈ C ∞ (S), for the L 2 -completion of the space of u ∈ C ∞ (S) with B(u| ∂M ) = 0 and B((Du)| ∂M ) = 0, see [18] for an elementary proof based on the Weitzenböck formula, which does not require pseudodifferential arguments, or [7] for a proof based on pseudodifferential calculus.
The use of APS boundary conditions also leads to the aforementioned properties of the Dirac operator and Dirac-Laplacian. However, the local elliptic boundary conditions determined by boundary chirality operators are more suitable for applications to spectral action models of gravity, see [9], hence it is more natural in our setting to use those.
For simplicity, we suppress writing the dependence of u n on the spatial coordinates η, φ, and ψ on S 3 , since those will not be directly affected by the transformations we consider, and we assume it understood.
In a similar way, using the transformation (3.8) of w k [p, q](iµ) and F [p, q](iµ) under the modular generator T in the expression (2.4) for the Dirac operator, we see that, if D[p, q]u n = λ n u n , thenD which implies that the operatorsD 2 [p, q + p + 1 2 ] andD 2 [p, q] are also isospectral. Finally, it is easy to see that the transformations u n (µ) → u n (µ + i) and u n (µ) → γ 0 u n ( 1 µ ) are two-sided inverses of the two eigenspinor transformations considered above, so the eigenspinor transformations are indeed bijections between the sets of eigenspinors.
In the case of the one-parameter family of Bianchi IX gravitational instantons we have a similar result.
Theorem 4.2. For the one-parameter family of metrics characterized by q 0 , the operators Proof. The argument is similar to the two-parameter case, using the transformations (3.12), (3.13), and (3.14) in the expression (2.4) of the Dirac operator.

Modularity of the heat kernel and the Seeley-deWitt coefficientsα 2n
The heat kernel K t of the operator exp −tD 2 can be written as a sum that involves its eigenvalues and their corresponding eigenspinors. Thus, it follows immediately from Theorem 4.1 that the heat kernel for the Bianchi IX gravitational instantons with parameters (p, q) satisfies modular transformation properties. More precisely, we have the following result, where again we write explicitly only the dependence on the µ-coordinate and we

JHEP01(2019)234
suppress the dependence on the spatial coordinates in S 3 . We use the same notation M , M T , and M S for the manifolds considered in Theorem 4.1.

Proposition 4.2.
In the case of the two-parameter family of Bianchi IX gravitational instantons, the heat kernel K t [p, q](iµ 1 , iµ 2 ) of the operator D 2 [p, q] satisfies the modular transformations Proof. Recall that the defining equation for the heat kernel is for any eigenspinor u n (µ) ofD 2 [p, q] on M . As a result, for (µ 1 , where u n (µ − i) is the general form of an eigenspinor ofD 2 [p, q + p + 1 2 ] acting on M T , according to Theorem 4.1. Similarly, for (µ 1 , µ 2 ) ∈ M S ×M S , we have ( 1 where −γ 0 u n ( 1 µ ) is the general form of an eigenspinor ofD 2 [−q, p] acting on M S according to Theorem 4.1. Thus, we see that K t [p, q](iµ 1 +1, iµ 2 +1) and − 1 The uniqueness of the heat kernel then implies that the modularity relations in the statement are satisfied.
Having established the modular transformation properties for the heat kernel K t [p, q], we can now show that all of the Seeley-deWitt coefficientsα 2n [p, q] inherit the same transformation properties.
Before stating the next theorem we recall from section 3.1 the notation ν P (f ) for the order of the zero of a function f at a point P (which is an integer).

JHEP01(2019)234
Theorem 4.3. For the two-parameter family of Bianchi IX gravitational instantons, and for any non-negative integer n, the Seeley-deWitt coefficients of the Dirac-Laplacian heat kernel satisfỹ In particular, all the Seeley-deWitt coefficients have the same zeros and poles with the same orders in H.
Proof. The Seeley-deWitt coefficientsα 2n [p, q](iµ) are uniquely determined by the asymptotic expansion or equivalently, Thus, using Proposition 4.2, we have Since the left and the right hand side of the latter have the following small time asymptotic expansions respectively, it follows from the uniqueness of the asymptotic expansion that Also, in a similar manner, we can write

JHEP01(2019)234
which implies thatα If v is the order of zero in Q of the function K t [p, q](Q) at Q 0 , then for t ∈ (0, 1) and some finite C t [p, q] = 0. On the other hand, the asymptotic expansion means that for any integer k there is some N (k) such that See sections 1.1 and 1.7 of [22] for the latter inequality and the definition of the norm | · | ∞,k . So it follows that where the convergence is uniform. Consequently, suppose that for some finite C n,t [p, q] = 0. Then we can switch the order of the two limits below and obtain As a result, for C n [p, q] to be finite and nonzero, we need to have v n = v, which proves that the Seeley-deWitt coefficients have the same zeros and poles of the same orders.
The case of the one-parameter family of Bianchi IX gravitational instantons can be approached similarly, using Theorem 4.2. We obtain the following modularity property for the Seeley-deWitt coefficients.

JHEP01(2019)234
Proof. Since the operatorsD 2 [q 0 ],D 2 [q 0 − i] and −q −2 0D 2 [ 1 q 0 ] are isospectral by Theorem 4.2, they have the same small time heat kernel expansions. Namely, for all non-negative integers n we haveα Therefore, for arbitrary real numbers a and b, such that the interval (a, b) does not contain any singularity of the metric, we have This show that which is equivalent to the statement of this theorem.

Modular forms in the full expansion of the spectral action
In order to properly interpret the modularity conditions for the Seeley-deWitt coefficients α 2n [p, q](µ) andα[q 0 ](µ) of Theorem 4.3 and 4.4, we show here that, for rational values of the parameters (p, q) the Seeley-deWitt coefficientsα 2n [p, q] give rise to vector-valued modular functions, in the sense of [14]. We then show that, by summation over a finite orbit of the parameters, these vectorvalued modular functions give rise to ordinary modular functions, and we investigate their relation to well known modular forms. In particular, we discuss two specific examples, one with poles at infinity and the other with no poles at infinity, where the modular functions corresponding to the Seeley-deWitt coefficients belong to the space of modular forms of weight 14 and to the cusp forms of weight 18, respectively.

The Seeley-deWitt coefficients as vector-valued modular forms
The following lemma shows the periodicity of all Seeley-deWitt coefficientsα 2n [p, q] in both of the parameters of the metric with period 1. This is a crucial step for showing that eachα 2n is defining a vector-valued modular form with respect to a finite dimensional representation of the group SL 2 (Z). More importantly, it also allows one to construct ordinary modular functions from eachα 2n [p, q], which can then be related to well known modular forms.

JHEP01(2019)234
Lemma 5.1. For any non-negative integer n and any parameters (p, q) of the Bianchi IX gravitational instantons we havẽ Proof. We know from the explicit form of the metric (2.14) and the transformations (3.11) that all the involved functions are invariant under p → p + 1, so the periodicity in p follows trivially. In addition, we have Consequently, we see from the equations above that the metric (2.1), for the two-parameter family of Bianchi IX gravitational instantons is invariant under q → q +1, so the periodicity in q also follows.
Thus, the orbit O (p,q) is finite, with any element of SL 2 (Z) acting as a permutation on O (p,q) . We can now show that each Seeley-deWitt coefficient is a vector-valued modular function as in the definition given by (3.7).
where on the left-hand-side M is acting on iµ by linear fractional transformations and on the right-hand-side M is acting on the parameters (p, q) according to the representation of (5.1). Therefore, by definition,α 2n [p, q](iµ) is a vector-valued modular function of weight 2.
Since the orbit is finite for a rational choice of parameters, by a summation over the orbit we obtain ordinary modular functions as follows. Recall that modular forms of weight 2 are meromorphic differentials on H/PSL 2 (Z). As we will discuss more in detail below, the modular forms we obtain in this way are indeed meromorphic and we provide more information later about their zeros and poles.
Corollary 5.1. For any pair of rational number (p, q), and any non-negative integer n, the sumα is a modular function of weight 2 for the modular group PSL 2 (Z).
Proof. Summing the components of the column vectorÃ 2n iµ; O (p,q) in Theorem 5.1, we find that the sumα 2n iµ; O (p,q) satisfies for any M ∈ PSL 2 (Z) which acts on the variable iµ as a linear fractional transformation, and on (p, q) as defined previously.
In the following we refer to the functionsα 2n (iµ; O (p,q) ) as averaged Seeley-deWitt coefficients.

Zeros and poles of the modular Seeley-deWitt coefficients
We have seen in Corollary 5.1 that, when the parameters (p, q) of the metric are rational, eachα 2n iµ; O (p,q) satisfies the transformation properties of a modular function of weight JHEP01(2019)234 2. Since there are no non-trivial holomorphic modular forms of weight 2, it is necessary forα 2n iµ; O (p,q) to have poles in the variable iµ.
By Theorem 4.3, in order to find the locations and multiplicities of the zeros of α 2n [p, q](iµ), it is enough for us to investigate the poles ofα 0 [p, q](iµ), and the result would apply to allα 2n [p, q](iµ). Recall that 4 .
Since all the theta functions and theta derivatives are holomorphic for iµ ∈ H, the singularities ofα may appear only at the zeros of the function ∂ q ϑ[p, q](iµ). In addition, because of the modular properties, it would be enough for us to look for poles in the fundamental domain H/PSL 2 (Z) and at infinity. Consequently, in principle, we only need to know the locations and multiplicities of the zeros of ∂ q ϑ[p, q](iµ) in order to figure out the space of modular forms that can be constructed fromα 2n iµ; O (p,q) by removing the poles. So we proceed by proving the following lemmas concerning the zeros of ϑ[p, q](iµ) and ∂ q ϑ[p, q](iµ) with parameters (p, q) ∈ [0, 1) 2 .
Proof. This follows directly from the leading order terms in (2.12) and its q derivative.
Lemma 5.2 can be used to obtain information about the order of zero ofα 0 [p, q] at infinity as follows.

JHEP01(2019)234
Then the function where n is an even number, hence, f (iµ) is a modular function of weight n. Therefore using the valence formula (3.6) we can write Recall that the zeros of ϑ[p, q](iµ) satisfy the equation for some m, k ∈ Z. So, for any real (p, q) = ( 1 2 , 1 2 ), ϑ[p, q](iµ) does not have any zeros (or poles) in H. Moreover, notice that the theta functions and the theta derivatives are all holomorphic in iµ ∈ H. So it follows that the function f (iµ) has no pole away from infinity, and for P ∈ H we have exactly where n 0 is the number of points in the orbit whose first coordinates are zero, p ′ = 0. Finally it follows from the valence formula for f that
In this case also we can compute the explicit Q-expansion as we find that indeed for the first three averaged Seeley-deWitt coefficients we havẽ

Conclusions
Unlike the round metric of the Robertson-Walker spacetimes, the Bianchi IX metrics are only left-invariant under the action of SU(2) but not right-invariant. This SU(2) symmetry is encoded in the invariant 1-forms of (2.1). We have shown that this SU(2) symmetry of the Bianchi IX gravity models provides the ground for the existence of arithmetic structures in the resulting spectral model of gravity, which can be identified explicitly with mathematical methods. More precisely, we have taken the rationality phenomenon presented in Proposition 2.1, which is due to the symmetries, as a first strong indication of the presence of deeper arithmetic structures hidden in the Seeley-deWitt coefficients of the small time heat kernel expansions for these geometries. We have then identified explicitly this arithmetic structure by showing that all the Seeley-deWitt coefficients associated with the Dirac-Laplacian of Bianchi IX gravitational instantons are vector valued modular forms. By studying their zeros and poles we have illustrated how they are related to well known modular forms. We should stress that in general it is formidable to have explicit formulas for the Seeley-deWitt coefficients beyond the first few terms [1,3,22,40,42], as it is quite difficult to keep pace with the rapid growth in complexity of the terms. However, by exploiting the fact that the solutions of the Bianchi IX instantons are given by theta functions, as written in section 2.2 from [4], we have used the action of the modular group and thereby JHEP01(2019)234 the theory of modular forms to obtain conceptual results about the structure of the full heat kernel expansion for the Dirac-Laplacian of the instantons. A crucial step in our work was to perform lengthy calculations with the explicit formulas for the first few terms in the expansion, which are computable, to identify the existing modular properties. We then sought to identify a conceptual reason for the observed modularity, and we proved that it holds because of the isospectrality of the Dirac operators of the metrics related by the action of the modular group on the parameters defining the family of instantons. With this method we proved modularity of all the Seeley-deWitt coefficients in the expansion. Different areas of mathematics and physics are familiar with this type of extraordinary performance of modular forms in identification of the structure of complicated expressions [45], as modular forms typically live in low dimensional vector spaces, and there are explicit algorithms for relating them to well understood modular forms [38].
Part of the motivation for the results described in this paper about the heat kernel of the Dirac-Laplacian on the Bianchi IX gravitational instantons comes from the construction of models of Euclidean gravity based on the spectral action functional [10]. This approach extends to noncommutative spaces and was used to construct models of gravity coupled to matter, as well as for applications to cosmological models, see [12,32] for an overview. The spectral action functional is defined as a regularized trace of the Dirac operator, where the regularization is given by a smooth approximation f to a cutoff function on the real line, Trace f (D/Λ) , and Λ is an energy scale. It is shown in [10] that the spectral action has an asymptotic expansion for large Λ of the form Trace f (D/Λ) ∼ β∈Π f β Λ β − |D| −β + f (0)ζ D (0) + · · · , (6.1) where the terms in the expansion are associated to poles β of the zeta function ζ D (s) of the Dirac operator, and a more general family of associated zeta functions, and the coefficients in each term depend on the momenta of the test function f and on the residues of the zeta function. Mellin transform, in turn, relates these residues to the Seeley-deWitt coefficients of the heat kernel of the Dirac-Laplacian. The spectral action can be regarded as the action functional of a modified gravity model, since the leading terms of its asymptotic expansion recover the usual Einstein-Hilbert action for gravity with cosmological constant, as well as additional higher derivative terms such as conformal gravity (Weyl curvature) and Gauss-Bonnet gravity. More precisely, when written non-perturbatively as Trace(f (D/Λ)) the spectral action is viewed as a functional on a configuration space of Dirac operators D, and a corresponding functional integral should be taken over this space of possible Dirac geometries. When seen perturbatively through the asymptotic expansion obtained from the heat kernel, one sees the individual coefficients as higher derivative correction terms to the classical action functional of general relativity. The results obtained in this paper on the Seeley-deWitt coefficients of the Bianchi IX gravitational instantons being vector-valued modular forms translate into JHEP01(2019)234 an analogous statement for the terms in the asymptotic expansion of the spectral action functional on these geometries. The presence of arithmetic structures in spectral models of gravity is observed in [17,21] from a different point of view, where the coefficients of the heat kernel expansion (hence those of the spectral action functional) are expressed in terms of periods of mixed Tate motives. This type of result is in a sense analogous to the well known occurrence of motives and periods in the asymptotic expansions of perturbative quantum field theory, see [33] for an overview. The results of the present paper show that arithmetic structures are also present in the form of modular forms in the spectral action of the Bianchi IX gravitational instantons.
Finally, we mention that considering the approach to quantum gravity based on the Wick rotation and thereby the heat kernel expansion in the realm of elliptic operators (cf. [2]), our results have important implications in gravity and cosmology models, which will be discussed elsewhere (see section 1.1 for a brief outline). Also, the geometries considered in the present article are especially interesting for a theory of Euclidean quantum gravity and quantum cosmology because they provide the minisuperspace model approximation in Hartle-Hawking gravity, see [15].