Modular forms in the spectral action of Bianchi IX gravitational instantons

Following our result on rationality of the spectral action for Bianchi type-IX cosmological models, which suggests the existence of a rich arithmetic structure in the action, we study the arithmetic and modular properties of the action for an especially interesting family of metrics, namely $SU(2)$-invariant Bianchi IX gravitational instantons. A variant of the previous rationality result is first presented for a time-dependent conformal perturbation of the triaxial Bianchi IX metric which serves as a general form of Bianchi IX gravitational instantons. By exploiting a parametrization of the gravitational instantons in terms of theta functions with characteristics, we show that each Seeley-de Witt coefficient $\tilde a_{2n}$ appearing in the expansion of their spectral action gives rise to vector-valued and ordinary modular functions of weight 2. We investigate the modular properties of the first three terms $\tilde a_0$, $\tilde a_2$ and $\tilde a_4$ by preforming direct calculations. This guides us to provide spectral arguments for general terms $\tilde a_{2n}$, which are based on isospectrality of the involved Dirac operators. Special attention is paid to the relation of the new modular functions with well-known modular forms by studying explicitly two different cases of metrics with pairs of rational parameters that belong to two different general families. In the first case, it is shown that by running a summation over a finite orbit of the rational parameters, up to multiplication by the cusp form $\Delta$ of weight 12, each term $\tilde a_{2n}$ lands in the 1-dimensional linear space of modular forms of weight 14. In the second case, it is proved that, after a summation over a finite orbit of the parameters and multiplication by $G_4^4$, where $G_4$ is the Eisenstein series of weight 4, each $\tilde a_{2n}$ lands in the 1-dimensional linear space of cusp forms of weight 18.


Introduction
The main result of this paper shows that the spectral action of Bianchi IX gravitational instantons is arithmetic, in the sense that its asymptotic expansion can be expressed in terms of rational combinations of (vector valued) modular forms, which in turn can be explicitly related to classical modular forms of weight 14 and 18. The rationality of the spectral action for a general triaxial Bianchi type-IX metric with an SU (2)-symmetry, obtained in [19], suggested the existence of a rich arithmetic structure in the Seeley-de Witt coefficients associated with the square of the Dirac operator of these cosmological models. Here the rationality assertion means that each coefficient is the time integral of an expression presented by a several variable polynomial with rational coefficients evaluated on the expansion factors and their derivatives, up to a certain order with respect to time. An earlier rationality result of a similar nature was obtained in [21] for the Robertson-Walker metrics, proving a conjecture formulated in [6].
The present article is intended to obtain a deeper understanding of the arithmetic properties of the spectral action, and to shed light on the role of the rational coefficients appearing in the expansion of the spectral action for Bianchi IX metrics. By imposing the condition of self-duality of the Weyl tensor and by employing a time-dependent conformal factor to obtain an Einstein metric from the Bianchi IX models, an especially interesting family of metrics called Bianchi IX gravitational instantons have been explored and well studied in the literature, see for example [48,45,27,35,2] and references therein. Interestingly, as explained in great detail in the latter, the differential equations for finding these metrics reduce to well understood equations such as the Halphen system and the Painlevé VI equation with particular rational parameters. In [2], following the work carried out in [48,27], these equations are solved by using the τ -function of the Schlesinger system formulated in terms of theta functions [30], and an explicit parametrization of the Bianchi IX gravitational instantons is given in terms of theta functions with characteristics. Considered along with our rationality result about the spectral action [19], this parametrization of the gravitational instantons will be the main ingredient in our construction of the modular expression for the terms appearing in the expansion of the spectral action in the energy scale. We will also describe an explicit connection between the modular functions that arise in the spectral action and well-known classical modular forms.
In the next section, we introduce and clarify our notation, and we briefly review all the necessary background material on the spectral action that we need to use throughout the paper. While most of the material we need to recall is standard, we prefer to maintain the paper as readable and self-contained as possible. The reader already familiar with these notions and notations can skip this section. We start by an explanation of the spectral action functional [5,9], which is based on the Dirac operator and provides a modified Euclidean gravity model. Implications of this action as a source of new early universe models and inflationary mechanisms in cosmology has been studied in recent years [3,31,38,39,40,41,42,43,44,18]. We recall briefly some basic facts about the Dirac operator, the heat kernel method and pseudo-differential calculus, and how they can be employed to compute the terms in the asymptotic expansion of the spectral action in the energy scale. The terms that appear in the expansion include the Einstein-Hilbert action and other modified gravity terms such as the Weyl curvature and Gauss-Bonnet terms. Indeed the latter expressions appear as the first few terms in the expansion. As a general problem it is highly desirable to achieve an understanding of the full expansion. In [19], we devised an efficient method for computing the terms appearing in the expansion by using the Wodzicki noncommutative residue [49,50], a powerful tool that, in addition to be very important for convenient calculations, yields an elegant proof of the rationality result presented in [19]. After briefly recalling, for comparison purposes, the traditional method of computing these coefficients, we end Section 2 by describing the final formulation of our noncommutative residue method, which prepares the ground for deriving a variant of the rationality result for the specific family of metrics that serve as a general form for the Bianchi IX gravitational instantons. Section 3 is devoted to the explicit computation of the Dirac operatorD of a general time-dependent conformal perturbation of the triaxial Bianchi IX metric. We also prove a rationality statement for its spectral action. The general form of the Bianchi IX gravitational instantons, which we mentioned earlier, is (1) ds 2 = F ds 2 = F w 1 w 2 w 3 dµ 2 + w 2 w 3 w 1 σ 2 1 + where the conformal factor F is a function of the cosmic time µ. Here, the metric ds 2 is the Bianchi IX model with general cosmic expansion factors w 1 , w 2 and w 3 , and SU (2)-invariant 1-forms σ 1 , σ 2 and σ 3 . We found it convenient to recall from [19] the expression for the Dirac operator D of the Bianchi IX model ds 2 and to use it for the presentation of the Dirac operatorD of the conformally equivalent metric ds 2 = F ds 2 . We then obtain a rationality statement for the general terms a 2n appearing in the small time asymptotic expansion of the heat kernel, The section ends by a presentation of explicit expressions forã 0 andã 2 , while the lengthy expression forã 4 is recorded in Appendix B. As general notation, a tilde is used on top of any symbol that represents an object associated with the conformally perturbed metric ds 2 = F ds 2 . We will follow this notational convention throughout the paper.
In Section 4 we recall briefly another result that we need to use essentially in the rest of the paper: the derivation of explicit formulas for the Bianchi IX gravitational instantons obtained in [48,27,2]. One starts by imposing the self-duality condition on the Weyl tensor of the Bianchi IX model and employing a conformal factor to obtain an Einstein metric. These conditions reduce to the Halphen system and the Painlevé VI equation with particular rational parameters, which can then be solved in terms of elliptic modular functions and theta functions. We have written explicitly, in separate subsections, the parametrization of the solutions in terms of theta functions with characteristics given in [2]: one subsection for a twoparametric family with non-vanishing cosmological constants and the other for a one-parametric family whose cosmological constants vanish.
One of our main goals is to understand the modular properties of the Seeleyde Witt coefficientsã 2n appearing in (2) when the solutions of the gravitational instantons in terms of the theta functions are substituted in the metric (1). That is, by extending the real time µ to the right half-plane ℜ(µ) > 0 in the complex plane, to the extent that the argument of the theta functions that belongs to the upper half-plane H can be replaced by iµ, we explore the changes in theã 2n under modular transformations on H. Therefore, in Section 5 we deal with the modular properties of the theta functions and their derivatives that appear in the termsã 0 , a 2 andã 4 . Using these properties in the explicit expressions (24), (25) and the one recorded in Appendix B, we show by direct calculations in Section 6 that, under the linear fractional transformations T 1 (iµ) = iµ + 1 and S(iµ) = i/µ, the termsã 0 ,ã 2 andã 4 satisfy interesting modular properties. What is especially interesting here is that this behavior reveals modular transformation properties that are encoded in the parameters of the metric and that are of the same nature as those of vectorvalued modular forms considered in the Eichler-Zagier theory of Jacobi forms [17].
It is then natural to expect that these surprising modular transformation properties obtained by direct calculations for the coefficientsã 0 ,ã 2 andã 4 , will continue to hold for allã 2n . Clearly it is impossible to resort to direct calculations to investigate similar properties for the general terms. However, the properties seen in the first few terms, which are reflected in the parameters of the metric, indicate that there is a deeper relation between the Dirac operators associated with the metrics whose parameters transform to each other. Indeed, in Section 7 we study the corresponding Dirac operators and we prove that they are all isospectral. By taking advantage of this fact and of uniqueness of the coefficients in the asymptotic expansion, we prove that indeed all of the Seeley-de Witt coefficientsã 2n of the Bianchi IX gravitational instantons satisfy modular transformation properties.
In Section 8, we first prove a periodicity property for the termsã 2n with respect to the parameters of the two parametric family of Bianchi IX gravitational instantons. Combining this with their modular transformation properties, we show that, for rational parameters, each termã 2n defines a vector-valued modular function of weight 2 with values in a finite-dimensional representation of the modular group P SL 2 (Z). Recall that the modular group is generated by the matrices corresponding to the linear fractional transformations T 1 (iµ) = iµ + 1 and S(iµ) = i/µ acting on the upper half-plane, iµ ∈ H. We then observe that, by running a summation over a finite orbit of the modular transformations on the rational parameters, each a 2n gives rise to a modular function of weight 2 with respect to P SL 2 (Z). This type of modular functions are sometimes called quasi-modular, weakly modular or meromorphic modular forms to indicate that they are allowed to possess poles, which is indeed the case for any modular function of weight 2.
In the second part of Section 8, we find an intimate connection between the modular functions that arise from the Seeley-de Witt coefficientsã 2n and wellknown modular forms. We consider two pairs of rational parameters that belong to two different general families. In the first case, we prove that such modular functions have only simple poles at infinity, hence multiplication by the cusp form of weight 12, ∆(q) = q ∞ n=1 (1 − q n ) 24 , q = exp(2πiz), z ∈ H, lands them in the 1-dimensional linear space of modular forms of weight 14, generated by a single Eisenstein series. In the second case, we show that the only poles of the resulting modular functions are of order 4 and located at the point ρ = e 2πi/3 . Moreover, by showing that they vanish at infinity, we prove that multiplication by G 4 4 , where G 4 (z) = noncommutative residues of Laplacians. This method is remarkably efficient from a computational point of view and provided an elegant proof of the rationality result in [19], a variant of which has more direct relevance to the subject of the present paper, as explained in Section 3.
A noncommutative geometric space is described by a spectral triple, which consists of an involutive algebra A represented by bounded operators on a Hilbert space H, and an unbounded self-adjoint operator D on H [10]. The metric information is encoded in the operator D, which is assumed to satisfy the main regularity properties of the Dirac operator. The spectral action [5,9] for a spectral triple (A, H, D) is an action functional that depends on the spectrum of the Dirac operator D. It employs a cutoff function f defined on the real line to consider where Λ is the energy scale. Its asymptotic expansion in the energy scale is usually, depending on the nature of the arising poles, of the form [13] (3) Trace The summation in the latter runs over the points β where the poles of the spectral zeta function of the Dirac operator D, ζ D (s), and the poles of other associated zeta functions are located. The set of such points is called the dimension spectrum of the spectral triple [14]. For classical manifolds, the poles of the spectral zeta functions are located at certain points on the real line. However, in general, the dimension spectrum of a noncommutative geometric space can contain points in the complex plane that do not belong to the real line. See for instance [32] for examples of geometric zeta functions whose poles are not necessarily located on the real line. The main commutative example of a spectral triple is given by a spin c manifold M [11], where the algebra of smooth functions A = C ∞ (M ) acts on the L 2 -spinors of M , and D is the Dirac operator. In this case, the coefficients in the expansion of the spectral action are determined by the Seeley-de Witt coefficients associated with D 2 , which are local invariants of the geometry [25]. That is, up to considerations that merely rely on momenta of the cutoff function f , the terms of the expansion (3) are determined by the coefficients appearing in a small time asymptotic expansion of the form Chapter 1 of the book [13] contains a detailed mathematical discussion of the asymptotic expansions related to the spectral action.
2.1. Dirac operator and its pseudo-differential symbol. Given a spin bundle S and a spin connection ∇ S on a Riemannian manifold M of dimension m, the Dirac operator D acting on smooth sections of S is defined by the following composition. It is given by composing the three maps where the second arrow is essentially the musical isomorphism # identifying the cotangent and tangent bundles, and the third arrow 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 lifting the Levi-Civita connection ∇ seen as a connection on the SO(m)-principal bundle of T M . Thus, finding an explicit local formula for the Levi-Civita connection is the first step in calculating the Dirac operator.
The Levi-Civita connection ∇ : is the unique connection on the tangent bundle that is compatible with the metric and torsionfree. In fact, these conditions characterize this connection uniquely and it can be computed by the method that will be explained shortly. Let us first explain that compatibility with the metric g means that Moreover, torsion-freeness refers to vanishing of the torsion tensor Since it is more convenient to work with coframes rather than with frames, for explicit calculations it is useful to take advantage of the identification of the tangent and cotangent bundles via the metric and write the Levi-Civita connection as a map from C ∞ (T * M ) to C ∞ (T * M ⊗T * M ). Then, if {θ a } is a local orthonormal coframe, i.e. an orthonormal basis for the local sections of T * M in a local chart U where it is trivialized, one can write where ω a b are local differential 1-forms. In this basis ∇ can be expressed as where d is the de Rham differential and ω = (ω a b ) is a matrix of 1-forms. In this picture, the compatibility with the metric and the torsion-freeness respectively read as The latter conditions yield a unique solution for the 1-forms ω a b and thereby one achieves an explicit calculation of ∇.
Having computed ∇, one can lift it to the spin connection as follows. The spin group Spin(m) is a double cover of the special orthogonal group SO(m) and there is an explicit isomorphism µ : so(m) → spin(m) identifying their Lie algebras, which is given by (see for example Lemma 4.8 in [46]), In the latter {e a } is the standard basis of R m viewed inside the corresponding Clifford algebra, in which the linear span of the elements {e a e b ; a < b} is the Lie algebra spin(m) of Spin(m). Combining this with uniqueness of the spin representation for the Clifford algebra, one can choose k × k matrices, k = rk(S), which satisfy the relations (γ a ) 2 = −I and γ a γ b + γ b γ a = 0 for a = b, and calculate the matrix of 1-forms ω S representing the spin connection That is, one can write Now we write the Dirac operator explicitly: it follows from its definition and the above explicit formulas for the ingredients of the definition that where {θ a } is a pre-dual of the coframe {θ a } and the ω b ac are defined by It is clear that D is a differential operator of order 1. Like any differential operator, or more generally, like any pseudo-differential operator, using the Fourier transform and the Fourier inversion formula, D can be expressed by its pseudodifferential symbol denoted by σ(D). The symbol σ(D) is defined locally from U × R m to M k (C) and allows one to write the action of D on a local section s as whereŝ is Fourier transform of s taken component-wisely. Note that the endomorphisms of the bundle are locally identified with M k (C). The expression given by (5) makes it clear that for x = (x 0 , x 1 , . . . , x m−1 ) ∈ U and ξ = (ξ 1 , ξ 2 , . . . , ξ m ) ∈ R m . For our purpose of studying the Seeley-de Witt coefficients appearing in the asymptotic expansion (4), we need to have the pseudo-differential symbol of D 2 . This can be achieved by finding an explicit formula for D 2 from (5), or more easily by using the composition rule for pseudo-differential symbols. In fact, pseudo-differential operators are closed under composition and there is an explicit and handy formula, which describes the symbol of the product of two such operators modulo infinitely smoothing operators. That is, if the operators P 1 and P 2 are associated with the symbols σ(P 1 ) and σ(P 2 ), then, the symbol of P 1 P 2 has the following asymptotic expansion: One can find precise technical discussions in Chapter 1 of the book [25] about the pseudo-differential calculus that we use. In the case of differential operators, since the symbols are polynomials in ξ with coefficients that are matrix-valued functions defined on a local chart U , the formula (8) gives a precise formula for the symbol of the composition. Returning to the Dirac operator D, one can use its symbol given by (7) and the composition formula (8) to calculate the symbol of D 2 : The latter is in general a polynomial of order 2 in ξ whose coefficients are matrixvalued functions defined on the local chart U . In Section 3, explicit formulas are presented for the Dirac operators of the cosmological models that we are interested in.

2.2.
Heat expansion using pseudo-differential calculus. Let D be the Dirac operator on a compact spin manifold of dimension m, as we described. For any t > 0, the operator exp(−tD 2 ) is an infinitely smoothing operator and thus can be represented by a smooth kernel. In particular it is a trace-class operator and as t → 0 + , the trace of exp(−tD 2 ) goes to infinity. However, it is quite remarkable that there is an asymptotic expansion with geometric coefficients that describe the rate of this divergence. That is, as t → 0 + , where the coefficients a 2n (D 2 ) are local invariants of the metric that encode geometric information obtained from the curvature tensor. In fact, typically they are integrals of certain expressions obtained from the Riemann curvature tensor and its covariant derivatives and contractions. It is evident that the left hand side of (9) depends only on the eigenvalues of D 2 . Since, except in rare cases, in general the eigenvalues of the Dirac operator are not known, it is significant that there are methods in the literature that allow one to express the coefficients a 2n (D 2 ) appearing on the right hand side of (9) as integrals of local expressions obtained from the metric. One of these methods, which is quite effective, starts with the Cauchy integral formula and employs parametric pseudodifferential calculus to approximate the kernel exp(−tD 2 ) and thereby accomplishes a recursive procedure for finding formulas for the coefficients a 2n (D 2 ).
Let us review this method briefly from Chapter 1 of the book [25]. The Dirac operator D is a self-adjoint unbounded operator, with respect to which the Hilbert space of L 2 -spinors admits a spectral decomposition. Invoking the Cauchy integral formula one can write (10) exp where the integration is over a contour γ in the complex plane that goes around the non-negative real numbers clockwise. Since D 2 − λ is an elliptic differential operator, it admits a parametrix which is the same as (D 2 − λ) −1 modulo an infinitely smoothing operator. Thus, the approximation of (D 2 − λ) −1 amounts to finding or approximating the parametrix R λ of D 2 − λ. This is achieved by exploiting the calculus of pseudo-differential symbols given by the composition rule (8). In order to compute the symbol of R λ , since D 2 − λ is of order 2, the leading symbol of R λ has to be of order −2 and one can write where each r j (x, ξ, λ) is a parametric pseudo-differential symbol of order −2 − j.
There is an important nuance that λ should be treated of order 2. Also, for the parametric pseudo-differential symbols depending on the complex parameter λ, being homogeneous of order −2 − j means that for any t > 0, As we discussed before, the square of the Dirac operator D 2 is a differential operator of order 2 and therefore for the symbol of D 2 − λ we have where each p k (x, ξ) is a polynomial in ξ whose coefficients are matrix-valued functions defined on the local chart. By passing to the symbols and using the composition rule (8), the solution of the equation (11) σ R λ (D 2 − λ) ∼ I yields the following recursive solution of the terms r j in the expansion of the symbol of the parametrix R λ . In fact, by a comparison of homogeneous terms on the two sides of (11), one finds that and for any n ≥ 1, where the summations are over all 4-tuples of non-negative integers α, and integers 0 ≤ j < n and 0 ≤ k ≤ 2 such that |α| + j + 2 − k = n.
Using this recursive procedure, one can choose a large enough N so that the operator corresponding to the symbol r 0 + · · · + r N gives a desired approximation of the parametrix R λ of D 2 − λ. By substituting the approximation in the Cauchy integral formula (10), one obtains an approximation of the kernel of exp(−tD 2 ). By integrating the approximation of the kernel over the diagonal of M × M against the volume form, one can derive the small time asymptotic expansion (9).
It is remarkable that this method shows instructively that each coefficient in the expansion is given by the integral, where a 2n (x, D 2 ) is an invariantly defined function on the manifold defined in the local chart by The analysis involved for deriving the above expansion and formula for the coefficient a 2n (D 2 ) is quite intricate and as we mentioned earlier, it is presented in detail in Chapter 1 of the book [25]. It should be stressed that the integrals involved in the expression for a 2n (D 2 ) are possible to work out since one can show by induction from the formulas (12) and (13) that tr r n (x, ξ, λ) = n=2j−|α|−2 |α|≤3n r n,j,α (x) ξ α tr r 0 (x, ξ, λ) j .
Therefore, using the method reviewed in this subsection, one can calculate the a 2n (D 2 ) explicitly in concrete examples. However, it should be noted that this method involves heavy calculations that are cumbersome even with computer assistance. In Subsection 2.3, we review an efficient method that we devised in [19] for computing the Seeley-de With coefficients a 2n (D 2 ) by making use of Wodzicki's noncommutative residue [49,50] which in general is the unique trace functional on the algebra of classical pseudo-differential on a vector bundle on M . This method is significantly more convenient from a computational point of view and yields elegant proofs for rationality statements of the type discussed in Section 3.

2.3.
Calculation of heat coefficients using the noncommutative residue.
The symbol of a classical pseudo-differential operator of order d acting on the smooth sections of a vector bundle on an m-dimensional manifold M admits by definition an expansion of the following form as ξ → ∞: The noncommutative residue [49,50] of the pseudo-differential operator P σ associated with a classical symbol σ of the above type is defined by (17) Res Some explanations are in order for the latter. First, note that m is the dimension of the manifold, which shows that the noncommutative residue vanishes on the classical operators of order less than −m = −dim(M ). In particular, it vanishes on infinitely smoothing operators, which allows one to view the noncommutative residue as a linear functional defined on the space of classical symbols with the following rule for composing to classical symbols σ 1 and σ 2 , inherited from (8): Second, S * M = {(x, ξ) ∈ T * M ; ||ξ|| g = 1} is the cosphere bundle of M and the formula (17) is the integral over M of a 1-density associated with the classical symbol σ, which is called the Wodzicki residue density. In each cotangent fibre where ξ belongs to, consider the volume form on the unit sphere |ξ| = 1, Then, the mentioned 1-density associated with the classical symbol σ with the expansion (16) is defined by Extensive discussions and alternative formulations of the noncommutative residue, which was first discovered for 1-dimensional symbols on the circle [1,34], are given in [49,50,28] and Section 7.3 of the book [26]. The spectral formulation of this residue plays an important role in noncommutative geometry since it is used in the local index formula for spectral triples, developed in [14]. Also, formulas similar to the one given by (17) are used in [23,22] to define noncommutative residues for the restriction of Connes' pseudo-differential calculus [8] to noncommutative tori, which are handy computational tools, see for example [20] for an application. Now let D be the Dirac operator on a 4-dimensional manifold. We are interested in computing and understanding the nature of the Seeley-de Witt coefficients a 2n (D 2 ) appearing in an asymptotic expansion of the form (9) with m = 4, when D is of this type, namely the Dirac operator of a geometry describing a cosmological model. Using an alternative spectral formulation of the noncommutative residue and the Künneth formula, we showed in [19] that for any integer n ≥ 1, a 2n (D 2 ) = 1 32 π n+3 Res(∆ −1 ), where ∆ −1 denotes the parametrix of the elliptic differential operator in which ∆ T 2n−2 is the flat Laplacian on the (2n − 2)-dimensional torus T 2n−2 = (R/Z) 2n−2 . Evidently ∆ acts on the smooth sections of the tensor product of the spin bundle of M and the 1-dimensional trivial bundle on the torus. Using the symbol of ∆, which is intimately related to the symbol of D 2 , we showed in Corollary 4.1 of [19] that is the homogeneous component of order −2n−2 in the asymptotic expansion of the symbol of the parametrix ∆ −1 of ∆ = D 2 ⊗ 1 + 1 ⊗ ∆ T 2n−2 . In the proof of the corollary, we explained in detail the calculation of a recursive formula for σ −2n−2 (∆ −1 ), and stressed that it has no dependence on coordinates of the torus, which is indicated in the formula (18).
As we mentioned in the beginning of this section, this method is used crucially for proving the rationality statements for Bianchi IX metrics, which are elaborated on in Section 3.

Rationality of spectral action for Bianchi IX metrics
The goal of this section is to present a variant of the rationality result proved in [19]. That is, we consider a time dependent conformal perturbation of the triaxial Bianchi IX metric treated in [19], and show that the terms appearing in the expansion of its spectral action in the energy scale are expressed by several variable polynomials with rational coefficients, evaluated on the expansion factors, the conformal factor and their time derivatives up to a certain order. The reason for considering this family of metrics is that they serve as a general form for the Bianchi IX gravitational instantons [48,27,2]. Indeed, combined with the parametrization of the latter in terms of theta functions with characteristics [2], the rationality result stimulates the construction of modular functions from the spectral action, which is carried out in the sequel.
For convenience, in passing, we recall the formalism and explicit calculation of the Dirac operator of the Bianchi IX metrics from [19], which can then be used for the presentation of the Dirac operator of the conformally perturbed metric, given by (1).

Triaxial Bianchi IX metrics. Euclidean Bianchi IX metrics are of the form
where the cosmic expansion factors w i are functions of the cosmic time µ. The σ i are left-invariant 1-forms on SU (2)-orbits satisfying In order to write this metric explicitly in a local chart, in [19], we parametrized the 3-dimensional sphere S 3 by (η, φ, ψ) → cos(η/2)e i(φ+ψ)/2 , sin(η/2)e i(φ−ψ)/2 , with the parameter ranges 0 ≤ η ≤ π, 0 ≤ φ < 2π, 0 ≤ ψ < 4π. We then wrote the metric (19) in the local coordinates x = (µ, η, φ, ψ), the expression of which was found to be Going through the definitions and the construction reviewed in Subsection 2.1, the Dirac operator of this metric was then explicitly calculated: The pseudo-differential symbol associated with the latter has the following expression: We also note that in our calculations we use the following gamma matrices γ 0 , γ 1 , γ 2 and γ 3 , which are respectively written as Using the methods reviewed in Subsection 2.2 and Subsection 2.3, the terms a 0 , a 2 and a 4 in the expansion of the spectral action for the metric (19) were computed. Moreover, the main result of [19] is that a general term a 2n in the expansion, modulo an integration with respect to µ, is of the form where Q 2n is a polynomial of several variables with rational coefficients. The rationality result was proved by exploiting the SU (2) invariance of the 1-forms σ i appearing in (19) and by making use of the method reviewed in Subsection 2.3.

3.2.
Time dependent conformal perturbations of Bianchi IX metrics. By making a correct choice of a conformal factor, an especially interesting family of Bianchi IX metrics called Bianchi IX gravitational instantons have been explicitly expressed in [2], which are Einstein metrics while having self-dual Weyl tensors.
It is remarkable that starting from writing a gravitational instanton with SU (2) symmetry in the general form, where, like the w i , F is also a function of the cosmic time µ, the solutions of the equations for the self-duality of the Weyl tenor and proportionality of the Ricci tensor to the metric are classified completely in terms of solutions to Painlevé VI integrable systems [27,45,48]. In turn, the latter can be solved [2] by using the τ -function of the Schlesinger system formulated in terms of theta functions [30]. We will review the explicit parametrization of the Bianchi IX gravitational instantons in Section 4. In this subsection we present a rationality statement for the spectral action of the metric (21). This result indicates the existence of an arithmetic structure in the spectral action of these metrics, which, combined with the parametrization in terms of theta functions with characteristics [2], leads to a construction of modular functions, which is one of our main objectives in this paper.
By an explicit calculation following the notions and the construction described in Subsection 2.1, we find that the Dirac operatorD of the metric (21) is given by where D is the Dirac operator given by (20), the Dirac operator of the metric (19).
We also find that Therefore, we have the pseudo-differential symbol ofD 2 , since that of D 2 was calculated in [19]. Now, by following a quite similar approach to the one taken in [19] for the rationality result, we present a variant of that result for the metric (21). That is, considering the asymptotic expansion eachã 2n is of the general form written in the following statement.
Theorem 3.1. The termã 2n in the above asymptotic expansion, modulo an integration with respect to µ, is of the form whereQ 2n is a polynomial of several variables with rational coefficients.
Proof. We provide an outline. One can exploit the SU (2)-invariance of the 1forms σ i appearing in the metric (21) to show that functions of the type (15), whose integrals give the coefficientsã 2n , have no spatial dependence when the metric is given by (21). Then, one employs the formula (18) along with the pseudodifferential symbol ofD 2 , and continues with similar arguments to those of Theorem 5.1 in [19].
Let us end this section by recording explicit expressions for the first few coefficients appearing in (23), which were computed in two different ways to confirm their validity. In fact we first computed them by the method reviewed in Subsection 2.2 leading to the formulas (14) and (15), and then confirmed that the expressions match precisely with the outcome of our calculations based on the formula (18) which used the noncommutative residue.
The first coefficient, which is the volume term, is given by The next term, which is the Einstein-Hilbert action term, has the following rather short expression, which indicates occurrence of remarkable simplifications in the final formula: The termã 4 , which is the Gauss-Bonnet term, also enjoys remarkable simplifications in its final formula, however, since it has a lengthier expression, we present it in Appendix B.

Bianchi IX gravitational instantons
There is an especially interesting family of metrics called Bianchi IX gravitational instantons, which have been explored in the literature by imposing the self-duality condition on triaxial Bianchi IX metrics and by employing a time-dependent conformal factor F (µ) to obtain an Einstein metric, see [48,27,2] and references therein. Let us provide an outline of some of the main ideas and steps for deriving these metrics, from the literature. In particular, we will then present the explicit parametrization of the solutions from [2] in terms of theta functions with characteristics.
The sought after gravitational instantons can be written in the general form Then, the differential equations derived from imposing the self-duality of the Weyl tensor and the condition of being an Einstein metric are solved [48,27,2] by turning them into well-studied systems of differential equations as follows. One can start by considering a basis of anti-self-dual 2-forms where all cyclic permutations (j, k, l) of (1, 2, 3) are understood to be considered. Consider the connection 1-forms α j k appearing in dϕ j = α j k ∧ ϕ k . This leads to writing α j k = A k w k σ k , where the functions A k satisfy the system of equations It can then be seen that the self-duality condition on the Weyl tensor yields the classical Halphen system Remarkably, the latter has well-known solutions in terms of the theta functions that will be defined shortly by (31).
One can define a new variable x = A1−A2 A3−A2 , which has an evident dependence on µ. Then the Halphen system (28) reduces to the following differential equation which is satisfied by the reciprocal of the elliptic modular function, see [48] and references therein: Therefore, one can solve the equation (28) and substitute the solution in (27). The latter can then be solved more conveniently by setting and by viewing x as the independent variable, which yields: It is well-known that these equations reduce to the Painlevé VI equation with particular parameters, which along with the condition of making the metric Einstein in the conformal class using a time-dependent conformal factor, one reduces to the following rational parameters [27,45,48]: .
It should be noted that in general a Painlevé VI equation with parameters (α, β, γ, δ) is of the form Going through the process outlined above and solving the involved equations in terms of elliptic theta functions [48,27,2], and using the formula for the τfunction of the Schlesinger equation [30] with an additional elegant calculation the square root of some expressions in [2], cf. [27], a parametrization of the Bianchi IX gravitational instantons can be given as follows.
The final solutions in [2] are written in terms of theta functions with characteristics, which for p, q, z, ∈ C, iµ ∈ H, are given by (29) ϑ Considering Jacobi's theta function defined by and by using the notation which sets z = 0 in (29), the following functions are also necessary to be introduced: We are now ready to write the explicit formulas for the Bianchi IX gravitational instantons presented in [2], in the following two subsections. The first family is two-parametric, which consists of the case of non-vanishing cosmological constants, and the second is a one-parametric family whose cosmological constants vanish. By studying the asymptotic behavior of these solutions as µ → ∞, it is shown in [36] that they approximate Eguchi-Hanson type gravitational instantons with w 1 = w 2 = w 3 [16] for large µ.
4.1. The two-parametric family with non-vanishing cosmological constants. The two-parametric family of solutions with parameters p, q ∈ C is given by substituting the following functions in the metric (26): 4.2. The one-parametric family with vanishing cosmological constants. The one-parametric family of the metrics with the parameter q 0 ∈ R is given by the following solutions that need to be substituted in the metric (26): where C is an arbitrary positive constant.

Arithmetics of Bianchi IX gravitational instantons
This section is devoted to the investigation of modular properties of the functions appearing in the formulas (32) and (33) and that of their derivatives. When the functions w 1 , w 2 , w 3 and F are substituted from the latter identities in the metric (26), Theorem 3.1 implies that the Seeley-de Witt coefficientsã 2n are rational functions of ϑ 2 , ϑ 3 , ϑ 4 , ϑ[p, q], ∂ q ϑ[p, q], e iπp and their derivatives with rational coefficients. Therefore, finding out modular properties of these theta functions and consequently that of the functions w 1 , w 2 , w 3 , F and their derivatives, will help us to investigate modular transformation laws of theã 2n , under modular transformations on iµ belonging to the upper half-plane H.
We begin studying the two-parametric case (32). First, let us note that for the derivatives of the function ϑ[p, q](iµ) given by (30), we have We also need to prove the following lemma, which will be of crucial use for proving the transformation properties investigated in this section.
Proof. It follows from the following identity and the Poisson summation formula: For convenience, we need to use the constants which arise naturally in exploring the following transformation properties.
The following lemma shows that the function ϑ[p, q](iµ) and its derivatives with respect to µ possess periodic and quasi-period properties with respect to the variables p and q.
Lemma 5.2. The functions ϑ[p, q] are holomorphic in the half-plane ℜ(µ) > 0 and satisfy the following properties: . Proof. One can start from the definition given by (30) and (29) to write the following for proving the first identity.
The second identity can be seen to hold by writing This proves the quasi-periodicity property of the theta function and its derivative with the quasi-period 1 in the q-variable. The periodicity properties with period 1 in the p-variable can be similarly investigated by writing We now focus on modular transformations. For convenience, we use the following notation for the linear fractional transformations corresponding to the generators of the modular group acting on the upper half-plane H in the complex plane: In the following lemma, we present the transformation properties of the ϑ[p, q](iµ) and its derivatives under the modular transformations T 1 , T 2 , and S on iµ ∈ H. Lemma 5.3. Let µ be a complex number in the right half-plane ℜ(µ) > 0. We have Proof. Considering the formulas (30) and (29) we can write This establishes the first two identities for the transformation properties of ϑ[p, q] and its derivatives with respect to the modular action T 1 given by (34); the third and the fourth identity follows immediately from the first and the second by applying them twice.
In order to investigate the transformation properties with respect to the action of S given in (34), we need to use Lemma 5.1, which allows us to write Also we have: Now we investigate the transformation properties of the functions ϑ 2 , ϑ 3 , ϑ 4 , given by (31) and their derivatives, under the same modular actions as the ones considered above, namely T 1 and S given by (34) transforming iµ in the upper half-plane. 22 Lemma 5.4. Let ℜ(µ) > 0. The functions ϑ 2 , ϑ 3 , ϑ 4 and their derivative satisfy the following modular transformation properties: Proof. The first identity follows from considering the formula (31) and writing For the next identity, which involves the action S(iµ) = i/µ, one needs to use Lemma 5.1: The analogous identities for the functions ϑ 3 , ϑ 4 , and their derivatives can be proved in a similar manner. In the case of ϑ 3 we have ∂ n µ ϑ 3 (iµ + 1) = ∂ n µ ϑ[0, 0](iµ + 1) = ∂ n µ ϑ[0, and for ϑ 4 we have 23 5.1. The case of the two-parametric family of metrics. Using the above lemmas, we proceed to work out the modular transformation rules for the functions w j [p, q] given by (32) and their derivatives. First, let us deal with the transformation T 1 given by (34).
Lemma 5.5. The functions w j [p, q] and their derivatives of an arbitrary order n ≥ 1 with respect to µ satisfy the following identities: Proof. Using Lemma 5.3 and Lemma 5.4, we have Similarly for w 2 and w 3 we can write: The identities for the arbitrary derivatives follow easily from differentiating the above equalities with respect to µ.
The following lemmas show that the functions w j and their derivatives satisfy some transformation laws with respect to the modular action S given by (34) as well. Let us start by the properties of w 1 .
Lemma 5.6. Assume that the variable µ belongs to the right half-plane ℜ(µ) > 0. The function w 1 [p, q] and its derivatives up to order 4 with respect to µ satisfy the following identities: Proof. Using lemmas 5.3 and 5.4 we have By taking a derivative with respect to µ, it follows from the latter that One can then continue by taking another derivative to obtain For the next higher derivatives we find that: We record in the next two lemmas the transformation rules for the functions w 2 [p, q], w 3 [p, q] and their derivatives up to order 4 with respect to the modular action S(iµ) = i µ . Their proofs are presented in Appendix A, which are similar to the proof of Lemma 5.6. First, we treat w 2 [p, q] and its derivatives with respect to µ noting that the action of S on iµ yields expressions in terms of w 2 [−q, p] and its derivatives.
Lemma 5.7. The function w 2 [p, q] and its derivatives up to order 4 with respect to µ, ℜ(µ) > 0, satisfy the following identities: Proof. It is given in Appendix A.
Now we present the transformation laws for w 3 [p, q] and its derivatives under the modular transformation S(iµ) = i µ . We note that w 3 [p, q] behaves similarly to w 1 [p, q] under the action of S. In fact, the statement of the following lemma can be obtained from that of Lemma 5.6 by swapping the indices 1 and 3. 26 Lemma 5.8. Assuming ℜ(µ) > 0, the function w 3 [p, q] and its derivatives up to order 4 with respect to µ satisfy the following identities: Proof. See Appendix A.
We also need to know how the function F [p, q] given in (32) and its derivatives, transform under the modular transformations on iµ in the upper half-plane. First, their properties with respect to the action of T 1 given by (34) are presented. Lemma 5.9. For µ in the right half-plane ℜ(µ) > 0, the function F [p, q] and its derivative of an arbitrary order n ≥ 1 satisfy the following properties: Proof. Using Lemma 5.3, we have The identity for the derivatives of F [p, q] follows easily from differentiating the latter with respect to µ.
The following lemma gives the transformation properties of F [p, q] and its derivatives with respect to the action S given by (34) on the upper half-plane. 27 Lemma 5.10. If ℜ(µ) > 0, then Proof. We can use Lemma 5.3 to write Taking the derivate of the latter, we have: We then continue by taking another derivative to obtain Continuing this process, we find the properties of the higher derivative of F [p, q]:

5.2.
The case of the one-parametric family of metrics. We now turn to the one-parameteric case (33), for which, similarly to the case of the two-parametric case (32) analysed in Subsection 5.1, the necessary modular properties can be investigated by using the lemmas proved in the beginning of this section. Since the proofs are based on direct calculations combined with making use of the lemmas 5.3 and 5.4, we just present the statements.
First, let us deal with the functions w j [q 0 ] given in (33). With respect to the modular transformation T 1 given by (34) acting on iµ in the upper half-plane, these functions and their derivatives satisfy the following properties.
Lemma 5.11. The functions w j [q 0 ] and their derivative of an arbitrary order n ≥ 0 with respect to µ satisfy these properties when ℜ(µ) > 0: The next three lemma presents the transformation properties of the functions w j [q 0 ] and their derivatives up to order 4 with respect to the modular action of S given by (34) on iµ.
Lemma 5.12. The function w 1 [q 0 ] and its derivatives up to order 4 with respect to µ satisfy the following identities provided that ℜ(µ) > 0: While in the latter lemma for w 1 [q 0 ], the expressions on the right hand sides of the identities involve w 3 [ 1 q0 ], in the following lemma we observe that for w 2 [q 0 ] similar properties hold, however, the index 2 does not change in the sense that the properties are expressed in terms of w 2 [ 1 q0 ].
The function w 3 [q 0 ] and its derivatives behave similarly to the case of w 1 [q 0 ]: the statement given in the following lemma can be obtained from swapping the indices 1 and 3 in Lemma 5.12.
Lemma 5.14. The function w 3 [q 0 ] and its derivatives up to order 4 with respect to µ satisfy the following identities provided that ℜ(µ) > 0: Finally we present the modular transformation properties of the function F [q 0 ] given in (33) and its derivatives. Lemma 5.15. Provided that the complex number µ belongs to the right half-plane ℜ(µ) > 0, the function F [q 0 ] and its derivatives with respect to µ satisfy the following identities: 6. Modular properties ofã 0 ,ã 2 andã 4 by direct calculations Our aim in this section is to investigate modular properties of the Seeley-de Witt coefficientã 0 ,ã 2 , andã 4 , which are respectively given by (24), (25), and in Appendix B, in the case of the gravitational instantons parametrized by (32) and (33). The termsã 2n appear in general in the asymptotic expansion (23) of the heat kernel ofD 2 , whereD with the expression (22) is the Dirac operator of a time dependent conformal perturbation of a triaxial Bianchi IX metric, given by (21). The latter serves as a general form of Bianchi IX gravitational instantons, as we explained and provided references about it in Section 4. Since the expressions in (32) are parametrized by a pair of parameters [p, q], we denote the corresponding Seeley-de Witt coefficients byã 2n [p, q]. Similarly, the termsã 2n associated with the one-parametric case (33) with the parameter q 0 will be denoted byã 2n [q 0 ].
Proof. The identities can be seen to hold by applying directly lemmas 5.6, 5.7, 5.8 and 5.10 to the explicit expressions forã 0 ,ã 2 andã 4 given respectively by (24), (25) and in Appendix B. In fact, for the first term we can writẽ For the next term, first we just consider its expression, replace iµ by S(iµ) = i/µ and use the mentioned lemmas to write: Now we turn our focus to the modular properties of the termsã 0 [q 0 ],ã 2 [q 0 ] and a 0 [q 0 ] associated with the one parametric case (33).
Proof. The identities follow directly from applying lemmas 5.11, 5.12, 5.13, 5.14 and 5.15 to the explicit expressions forã The direct calculations carried out in this section show that the first three Seeleyde Witt coefficients behave similarly to vector-valued modular forms [17] since the modular transformations act accordingly on the parameters of the metrics. This strong evidence indicates that there should be an interesting connection between the Dirac operators of the metrics that are related to each other by the modular transformations. Indeed, by taking a close look at the spectral properties of these Dirac operators, we prove in Section 7 that the modular transformation properties proved for the first three terms hold for all of the termsã 2n .

Modular properties ofã 2n using identities for Dirac operators
Motivated by the modular transformation properties presented in theorems 6.1, 6.2 and 6.3 for the termsã 0 ,ã 2 andã 4 , the main goal in this section is to prove that all of the Seeley-de Wittã 2n satisfy the same properties. Clearly, this cannot be achieved by resorting to direct calculations. However, the direct calculations performed in Section 6 for the first three terms reveal that the effect of the modular transformations on the Seeley-de Witt coefficients is reflected in a certain way on the parameters of the metric, and by studying the corresponding Dirac operators carefully in this section we can prove that the transformation properties hold for the general terms.
Note that the the gravitational instantons parametrized by (32) and (33) are given in terms of theta functions and their derivatives. Thus in the two-parametric case (32) with the parameters (p, q) as well as in the one-parametric case (33) with the parameter q 0 , we are allowed to consider an extension of the definitions to p, q, q 0 ∈ C with iµ ∈ H, excluding possible zeros of the theta functions and their derivatives that appear in the denominators. Thus with (p, q) and q 0 fixed, we may consider the operatorsD 2 [p, q] andD 2 [ 1 q0 ] acting on a spin bundle over the base manifold M 0 = I (a,b) × S 3 , where I (a,b) is an arbitrary horizontal path in the upper-half complex plane H. We may also consider the operatorsD 2 [p, q + p In the following theorem, we prove the isospectrality of the Dirac operators corresponding to the metrics whose parameters are related to each other via the effect of the modular transformations on the Seeley-de Witt coefficientsã 0 ,ã 2 and a 4 , presented in theorems 6.1 and 6.2. Proof. For the two-parametric family of metrics, suppose u n (µ, η, φ, ψ) is a section of the spin bundle such thatD[p, q]u n (µ, η, φ, ψ) = λ n u n (µ, η, φ, ψ). For simplicity, let us suppress writing the dependence of u n on η, φ, and ψ and assume that it is understood. Using the explicit expression (22) for the Dirac operator we have: ) .
Now we can use in the latter lemmas 5.6, 5.7, 5.8 and 5.10, which show how the modular transformation S(iµ) = i/µ affects the functions w j [p, q] and F [p, q]. This yields Thus we showed that, for any eigenspinor u n (µ) ofD[p, q] with eigenvalue λ n , the spinor −γ 0 u n ( 1 µ ) is an eigenspinor of −D[−q, p] with the same eigenvalue λ n . From the above equations, one can also verify immediately that for any eigenspinor u n (µ) of −D[−q, p] with eigenvalue λ n , the spinor γ 0 u n ( 1 µ ) is an eigenspinor ofD[p, q] with the same eigenvalue λ n . Consequently, we see that there is a bijection λ n → −λ n between the eigenvalues ofD[p, q] and those of −D[−q, p]. Since the eigenspinors of D 2 are exactly the eigenspinors ofD with the corresponding eigenvalues squared, it follows thatD 2 [p, q] andD 2 [−q, p] are isospectral.
In a similar manner, we also have that Therefore using lemmas 5.5 and 5.9 in the latter, which explain how the parameters of the functions w j and F are affected by the modular transformation T 1 (iµ) = iµ + 1, we can write: Thus the following identity is proved which implies that the operatorsD 2 [p, q + p + 1 2 ] andD 2 [p, q] are isospectral. Finally, it is easy to see that the transformations u n (µ) → u n (µ+i) and u n (µ) → γ 1 u n ( 1 µ ) are two-sided inverses of the two eigenspinor transformations defined in the theorem, so the eigenspinor transformations are indeed bijections between the sets of eigenspinors.
Since the kernel K t [p, q] of the operator exp −tD 2 [p, q] can be written as a sum that involves its eigenvalues and their corresponding eigenspinors, it follows immediately from Theorem 7.1 that indeed the heat kernels corresponding to the parameters (p, q) satisfy modular transformation properties.
Corollary 7.1. With the spatial dependence suppressed, we have the following modular transformation properties for the heat kernel K t [p, q](iµ 1 , iµ 2 ) of the operator D 2 [p, q]: Proof. Recall that the defining equation for the heat kernel is for any eigenspinor u n (µ) ofD 2 [p, q] acting 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 1 according to Theorem 7.1.
Thus we see that K t [p, q](iµ 1 +1, iµ 2 +1) and − 1 ) satisfy the defining equations of K t [p, q + p + 1 2 ](iµ 1 , iµ 2 ) and K t [−q, p](iµ 1 , iµ 2 ) respectively. Now, the uniqueness of the heat kernel implies that Having established the modular transformation properties for the heat kernel K t [p, q], we can now show that all of the Seeley-de Witt coefficientsã 2n [p, q] inherit the same properties from the kernel.
Corollary 7.2. For any non-negative integer n we havẽ In addition, at any point iµ = P ∈ H, let v P (f ) be the order of zero in Q of the function f (Q), where Q = e −πµ . Then we have So in particular, all of the above Seeley-de Witt coefficients have the same zeros with the same orders in H.
Proof. The Seeley-de Witt coefficientsã 2n [p, q](iµ) can be uniquely defined by the asymptotic expansion or equivalently, So, using Corollary 7.1, 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 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 the book [25] for the latter inequality and the definition of the norm | · | ∞,k . So it follows 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,t [p, q] to be finite and nonzero, we need to have v n = v, which proves that the Seeley-de Witt coefficients have the same zeros of the same orders.
As an important remark, it should be mentioned that the proof of the latter corollary covers the case of poles of meromorphic functions considered as zeros of negative orders.
We now turn our focus to the one-parametric case. We proved in Theorem 6.3 that for the one-parametric Bianchi IX gravitational instantons (33) as well the termsã 0 [q 0 ],ã 2 [q 0 ] andã 4 [q 0 ] satisfy modular transformation properties. In order to show that all of the termsã 2n [q 0 ] satisfy the properties mentioned in this theorem, similarly to the two-parametric case treated so far in this section, we can use the isospectrality of the involved Dirac operators. An immediate corollary of the latter theorem is the following statement for all of the Seeley-de Witt coefficientsã 2n [q 0 ]. Theorem 7.3. For any non-negative integer n, we have: ] are isospectral, they have the same small time heat kernel expansions, namely that, for all non-negative integers n we haveã Therefore, for arbitrary real numbers a and b we have This show that which is equivalent to the statement of this theorem.

Modular forms arising fromã 2n [p, q]
In this section, we use the modular transformation properties of the Seeley-de Witt coefficients studied in the previous sections to show that when the parameters of the metric are rational in the two-parametric case, they give rise to vectorvalued modular forms. Then we show that by a summation over a finite orbit of the parameters, they give rise to ordinary modular functions. At the end we investigate their connection with well-known modular forms. Indeed, we show that, in examples of two general cases, one with poles at infinity and the other with no poles at infinity, the modular functions corresponding to the Seeley-de Witt coefficients land in a direct way in the modular forms of weight 14 or in the cusp forms of weight 18. 8.1. Vector-valued and ordinary modular functions fromã 2n [p, q]. The following lemma shows the periodicity of all Seeley-de Witt 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 modular group P 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. Lemma 8.1. For any non-negative integer n and any parameters (p, q) of the Bianchi IX gravitational instantons we have: Proof. Recalling the ingredients of the explicit formulas for the metric from Subsection 4.1 and using Lemma 5.2 we know 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 is invariant under q → q + 1, so the periodicity in q also follows. Now, relying on Lemma 8.1, we can correspond the following maps to the generators of P SL 2 (Z) acting on the ordered pair (p, q) ∈ S = [0, 1) 2 = (R/Z) 2 : S(p, q) = (−q, p), in both of which the parameters are considered modulo 1. For rational p, q with N being a common multiple of 2 and the denominators of p and q, it follows immediately from the definitions above that the orbit O (p,q) of (p, q) under the action of P SL 2 (Z) consists of (p, q) pairs with p, q ∈ N = {0, 2/N, . . . , and thus O (p,q) is finite, with any element of P SL 2 (Z) acting as a permutation on O (p,q) .
Theorem 8.1. Starting from a pair of rational numbers (p, q), for any non-negative integer n, the termã 2n , is a vector-valued modular function of weight 2 for the modular group P SL 2 (Z).
Proof. Since the orbit O (p,q) is finite in this case, we can arrange the functions a 2n [p ′ , q ′ ](iµ) with (p ′ , q ′ ) ∈ O (p,q) into a finite column vectorÃ 2n iµ; O (p,q) of some dimension d. Since any M ∈ P SL 2 (Z) acts as a permutation on O (p,q) , we may denote by ρ : S d → GL(p, C) the natural permutation representation of S d which acts onÃ 2n iµ; O (p,q) by permuting its components in the corresponding way.

From Corollary 7.2 we know that
for any M ∈ P SL 2 (Z). So 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.
Corollary 8.1. For any pair of rational number (p, q), and any non-negative integer n, the sumã defines a modular function of weight 2 for the modular group P SL 2 (Z).
Proof. Summing up all the components of the column vectorÃ 2n iµ; O (p,q) in Theorem 8.1, we find that the sumã 2n iµ; O (p,q) satisfies for any M ∈ P SL 2 (Z) which acts on the variable iµ as a Möbius transformation, and on (p, q) as defined previously. Namely,ã 2n iµ; O (p,q) is a modular function of weight 2 with respect to P SL 2 (Z).

8.2.
Connection betweenã 2n [p, q] and well-known modular forms. Here we investigate the connection between modular functions of type constructed in Corollary 8.1 and well-known modular forms. It was seen in this corollary that when the parameters (p, q) of the metric are rational, eachã 2n iµ; O (p,q) has the standard modular transformation properties of a modular function of weight 2. However, 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µ. For a detailed discussion of holomorphic modular forms, Eisenstein series and some related fundamental results used in this subsection, one can for example refer to [47].
By using Corollary 7.2, 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/P SL 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µ) when the parameters p and q are both real. In fact, according to Lemma 8.1, p and q are defined modulo 1, so in what is to be presented we will restrict to (p, q) ∈ S = [0, 1) 2 .

Lemma 8.2.
Let v ∞ (F ) be the order of zero of any function F (iµ) at infinity, and where p is defined as the number p ≡ p mod 1 such that p ∈ [− 1 2 , 1 2 ). Proof. These results can be obtained directly by keeping only the leading order terms in the defining formula (29) and its q derivative.
The latter can be used to obtain information about the order of zero ofã 0 [p, q] at infinity as follows.
Proof. Both of the above statements can be proved by the substitution of the identities given in Lemma 8.2 into the formula Now we can prove the following statements which will play crucial roles in studying the poles of the modular functions that are studied in detail in the sequel and have been related to well-known modular forms. Then n is necessarily even, and with v τ [p, q] denoting the order of zero of the function ∂ q ϑ[p, q](iµ) at iµ = τ , we have the equation where ρ = e 2πi 3 , n 0 is the number of points (p ′ , q ′ ) ∈ O (p,q) such that p ′ = 0, and the starred summation excludes the points i and ρ.
Proof. First notice that if (p, q) and (−p, −q) correspond to the same point in S = [0, 1) 2 then we have p ≡ −p mod 1 and q ≡ −q mod 1, i.e., p ≡ 0 mod 1 2 and q ≡ 0 mod 1 2 . So in S, the only possibility is that .
, it has a partner (−p ′ , −q ′ ) = (p ′ , q ′ ) in S that is also in the orbit. It is evident that only identical points in S can have the same partner, so it follows that in this case n is necessarily even.
Recall from Lemma 5.3 that ∂ n q ϑ[p, q](iµ + 1) = e −πip(p+1) ∂ n q ϑ[p, q + p + Consequently, with n an even number, if we define the following function then we have 45 which shows that f (iµ) is a modular function of weight n. Therefore using the valence formula 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 in H. Also, 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 v P (f ) = (p ′ ,q ′ )∈O (p,q) v P [p ′ , q ′ ], with v P [p ′ , q ′ ] ≥ 0 at any P in the fundamental domain. At infinity, we know from 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 Proof. This follows directly from the non-negativity of v P [p, q].
When n is small, for instance when n = n 0 , one can solve for all the v P [p ′ , q ′ ] by simple arithmetic arguments. For example, when n = 8, one has where (p ′ , q ′ ) and (−p ′ , −q ′ ) correspond to different points in S. Therefore, the following can be the only non-negative solution: v ρ [p 0 , q 0 ]= v ρ [−p 0 , −q 0 ] = 1 for some (p 0 , q 0 ) ∈ O (p,q) , and v P [p ′ , q ′ ] = 0 for any other P in the fundamental domain and (p ′ , q ′ ) ∈ O (p,q) .
We now work out different examples explicitly. For instance, we look at the following orbit generated by (p, q) = (0, 1 3 ): , , .
For this case we have the following statement.
Theorem 8.2. For any non-negative integer n,ã 2n iµ; where ∆ is the modular discriminant (a cusp form of weight 12), and G 14 is the Eisenstein series of weight 14.
Proof. From the orbit above one observes that n = 24 = 6n 0 in this case, so Corollary 8.3 shows that ∂ q [p, q](iµ) has no zeros in the upper-half complex plane for any (p, q) in the orbit. Therefore,ã 0 iµ; O (0, 1 3 ) and thusã 2n iµ; O (0, 1 3 ) are holomorphic on H as we have argued previously.
In addition, by either Corollary 8.2 or direct expansion in Q = e −2πµ we know thatã 0 iµ; O (0, 1 3 ) has a simple pole at infinity, iµ = ∞, Q = 0. Since the modular discriminant ∆ has a zero of order 1 at iµ = ∞, it follows that ∆(iµ)·ã 0 iµ; is a modular function holomorphic on H as well as at iµ = ∞. Namely, ∆(iµ) · a 0 iµ; O (0, 1 3 ) is a modular form of weight 12 + 2 = 14. Since the space of modular forms of weight 14 is a one-dimensional space generated by the Eisenstein series G 14 , the desired result follows.
Indeed, by writing their Q-expansions explicitly, we confirm that forã 0 iµ; Using this approach and taking advantage of the lemmas proved here, one can prove similar results for many other orbits O (p,q) . As another example, we can also look at the following orbit generated by (p, q) = ( 1 6 , 5 6 ) containing 8 points: ) .
In this case, the statement is as follows, in which G 6 denotes the Eisenstein series of weight 6. Indeed, explicit Q-expansions in the latter case also confirm that:

Conclusions
The results obtained in this paper present a novel occurrence in quantum cosmology of modular functions and the vector-valued modular forms considered in the Eichler-Zagier theory of Jacobi forms [17]. This was indeed suggested to us by a combination of two different sources: the existence of an explicit parametrization of Bianchi IX gravitational instantons in terms of theta functions with characteristics [2] (see also [48,27]) and our rationality result about the Seeley-de Witt coefficients in the asymptotic expansion of the spectral action for triaxial Bianchi IX metrics [19]. These two results combined reveal that each Seeley-de Witt coefficient in the expansion is a rational function, with rational coefficients, in the theta functions ϑ 2 , ϑ 3 , ϑ 4 , ϑ[p, q], ∂ q ϑ[p, q], e iπp and their derivatives (the latter theta functions are written explicitly in Section 4).
Bianchi IX gravitational instantons are especially interesting since they admit an explicit parametrization in terms of elliptic modular functions and theta functions with characteristics [2], see also [48,45,27] and references therein. These are obtained by imposing the self-duality condition on the Weyl tensor of Bianchi IX metrics, and reducing the corresponding partial differential equations to well known ordinary differential equations, namely the Halphen system and the Painlevé VI equation. This result is followed by still another crucial step aimed at making the result an Einstein metric. That is, a correct choice of a time-dependent conformal factor is essential for making the Ricci tensor proportional to the metric. This fact is relevant to our present work in the following interesting ways. On the one hand, we have explained in this paper that a similar rationality result holds for a general time-dependent conformal perturbation of the triaxial Bianchi IX metric treated in [19]. The result is proved by employing our method based on Wodzicki's noncommutative residue [49,50] and the Künneth formula. On the other hand, it is necessary to involve the correct conformal factor in our calculations, in order to obtain the modular transformation properties that we discussed. These properties add to the many interesting and special features of the Bianchi IX gravitational instantons.
Modular forms appear in a variety of areas in mathematics and physics. Since modular forms of a certain weight form a finite dimensional linear space and can be computed with algorithmic methods, they have a wide range of applications. Thus, it is of great importance in general to find an explicit way of relating any modular function or modular form that arise from a mathematical structure or from physical problems to well-known modular forms, whose Fourier expansion, for example, is known. We have accomplished this task for the modular functions arising from the spectral action for Bianchi IX metrics, in this paper, by exploring their intimate connection with modular forms of weight 14 and cusp forms of weight 18, both of which form 1-dimensional linear spaces. That is, we have shown that, when the two parameters of a gravitational instanton are rational, belonging to two different general families, there is a finite orbit of the parameters, for each case, such that summation over the orbits leads to the following. In the first case, after multiplication by the cusp form ∆ of weight 12, each modular function arising from the Seeley-de Witt coefficientã 2n lands in the space of modular forms of weight 14. This indicates that each modular function arising in this case has only one simple pole, which is located at infinity. In the second case, after multiplication by G 4 4 , where G 4 is the Eisenstein series of weight 4, the modular functions arising from the Seeley-de Witt coefficients land in the 1-dimensional space of cusp forms of weight 18.
In order to illustrate how the present work fits in the general panorama of other rich arithmetic and number theoretic structures in theoretical physics, let us mention the following examples. A first example is the setting in which Feynman integrals are interpreted as periods, see [37] for an overview. In this case, the relevant amplitude forms and domains of integration are algebraic over the rationals or integers and this fact has direct implications on the class of numbers that arise as periods. In particular, an interesting connection to modular forms also arises in this setting [4]. A second example in which rational coefficients play an important role is in the zero temperature KMS states of quantum statistical mechanical systems. For example, in the construction in [12], which is explained also in Chapter 3 of [13], an arithmetic algebra of observables over the rationals is constructed, whose link to modular functions allows to have KMS states that take their values in the modular field. There are many occurrences of modular forms in physics, especially in the context of String Theory. The literature on the subject is extensive and we cannot mention all the relevant results here, so we only point the reader to a couple of significant recent examples, such as [7,15]. The setting we considered here is very different, as modular forms arise in the gravity action functional (the spectral action) of a specific class of gravitational instantons, rather than in settings such as superstring amplitudes, or counting functions for BPS states, or mirror symmetry. There are many other examples in the literature of arithmetic structures arising in physics, see for example the contributions collected in the volume [29]. As it is noticeable from the present work as well, it is in general a challenging and promising task to further explore the hidden arithmetic structures in different areas of physics, including gravity and quantum cosmology.
Appendix A. Proofs of Lemma 5.7 and Lemma 5.8 For the sake of completeness, we provide here the proofs of the modular transformation properties of the functions w 2 , w 3 and their derivatives that were stated in lemmas 5.7 and 5.8. In fact the proofs are very similar to that of Lemma 5.6. For the function w 2 , using lemmas 5.3 and 5.4, we can write By taking derivatives of the latter with respect to µ consecutively we obtain:  Similarly, for the function w 3 we have: The expressions for the termsã 0 andã 2 appearing in the asymptotic expansion (23) in whichD is the Dirac operator of the metric (21), were given at the end of Section 3. Since the termã 4 has a lengthy expression, it is recorded in this appendix. It is clear that this term also, in accordance with Theorem 3.1, possesses only rational coefficients: ) .
Now we can use lemmas 5.12, 5.13, 5.14 and 5.15, which explain the modular transformation properties of the functions w j and F in the one-parametric case (33) with respect to the modular transformation S(iµ) = i/µ. Using these lemmas in the above we expression, we can write Therefore, for any eigenspinor u n (µ) ofD[q 0 ] with the eigenvalue λ n , the spinor −γ 0 u n ( 1 µ ) is an eigenspinor of iq −1 0D [ 1 q0 ] with the same eigenvalue λ n . HenceD 2 [q 0 ] and −q −2