A central limit theorem for integrals of random waves

TL;DR

This paper establishes a central limit theorem for the mean-square of high-frequency random waves on compact Riemannian manifolds, utilizing advanced integral estimates and the local Weyl law.

Contribution

It introduces a universal CLT for random waves on any compact Riemannian manifold, extending previous results to arbitrary dimensions and shrinking sets.

Findings

Proves a CLT for the mean-square of random waves in high-frequency limit.

Uses a novel estimate involving Bessel functions and Gegenbauer's addition formula.

Demonstrates universality across different manifold geometries.

Abstract

We derive a central limit theorem for the mean-square of random waves in the high-frequency limit over shrinking sets. Our proof applies to any compact Riemannian manifold of arbitrary dimension, thanks to the universality of the local Weyl law. The key technical step is an estimate capturing some cancellation in a triple integral of Bessel functions, which we achieve using Gegenbauer's addition formula.

Figures

Equations

$$\varphi (x)=\sum _j{c}_j{\varphi }_j(x)$$

$$Z=Z_B=\frac{{\int }_B{\varphi }^2-\mathbb{E}[{\int }_B{\varphi }^2]}{\sqrt{\mathrm{var}[{\int }_B{\varphi }^2]}},$$

$$\mathbb{E}[e^{itG}]=e^{-t^2/2}$$

$$\mathbb{E}\left[e^{itZ}\right]=e^{-t^{2}/2}\left(1+O\big{(}(rT)^{-1/2}\big{)}\right)$$

$$Z_B\to N(0,1).$$

$$f(X)=O\big{(}g(X)\big{)},\quad f(X)\lesssim g(X)$$

$$C_{1}g(X)\le f(X)\le C_{2}g(X).$$

$$\log X{\lesssim }_{\epsilon }{X}^{\epsilon }.$$

$$Z=\frac{{\mathfrak{z}}^{T}A\mathfrak{z}-\mathbb{E}\left[{\mathfrak{z}}^{T}A\mathfrak{z}\right]}{\sqrt{\mathrm{var}\left[{\mathfrak{z}}^{T}A\mathfrak{z}\right]}}.$$

$$\frac{\mathrm{tr}(A^3)}{\mathrm{tr}(A^2{)}^{3/2}}\to 0$$

$$\mathbb{E}\left[e^{itZ}\right]=e^{-t^2/2}\left(1+O\left(|t{|}^3\frac{\mathrm{tr}(A^3)}{\mathrm{tr}(A^2{)}^{3/2}}\right)\right).$$

$$K(x,y)=\sum _j{\varphi }_j(x){\varphi }_j(y).$$

$$\frac{{\int }_B{\int }_B{\int }_{B}K(x_1,x_2)K(x_2,x_3)K(x_3,x_1)d{x}_{3}d{x}_{2}d{x}_1}{{\left({\int }_B{\int }_{B}K(x,x^{\prime }{)}^{2}dxd{x}^{\prime }\right)}^{3/2}}\to 0$$

$$\mathbb{E}\left[e^{itZ}\right]=e^{-t^2/2}\left(1+O\left(|t{|}^3\frac{{\int }_B{\int }_B{\int }_{B}K(x_1,x_2)K(x_2,x_3)K(x_3,x_1)d{x}_{3}d{x}_{2}d{x}_1}{{\left({\int }_B{\int }_{B}K(x,x^{\prime }{)}^{2}dxd{x}^{\prime }\right)}^{3/2}}\right)\right).$$

$$\int_{B}\int_{B}K(x,x^{\prime})^{2}dxdx^{\prime}\asymp(rT)^{-n+1}\big{(}\operatorname{vol}(B)T^{n-1}\eta\big{)}^{2}.$$

$$\int_{B}\int_{B}\int_{B}K(x_{1},x_{2})K(x_{2},x_{3})K(x_{3},x_{1})dx_{3}dx_{2}dx_{1}\lesssim(rT)^{-3n/2+1}\big{(}\operatorname{vol}(B)T^{n-1}\eta\big{)}^{3}.$$

$${\int }_A{\varphi }_{\lambda }^{2}d\mathrm{vol}\to \mathrm{vol}(A)$$

$$\mathbb{P}(|X-\mathbb{E}[X]|\ge y\sqrt{\mathrm{var}[X]})\lesssim e^{-y^2/2}$$

$$\mathbb{P}(|X-\mathbb{E}[X]|\ge y\sqrt{\mathrm{var}[X]})\lesssim e^{-c{y}^2}$$

$${\int }_B{\varphi }^2=\sum _j\sum _k{c}_j{c}_k{\int }_B{\varphi }_j{\varphi }_k$$

$${\int }_B{\varphi }^2={\mathfrak{z}}^{T}A\mathfrak{z}$$

$$A_{jk}=\mathrm{var}[c]{\int }_B{\varphi }_j{\varphi }_k.$$

$${\mathfrak{z}}^{T}A\mathfrak{z}=(U\mathfrak{z}{)}^{T}D(U\mathfrak{z})=\sum _j{\lambda }_j{y}_j^2$$

$$\mathbb{E}\left[e^{s{\mathfrak{z}}^{T}A\mathfrak{z}}\right]$$

$$\mathbb{E}\left[e^{s\lambda y^2}\right]={\int }_{-\infty }^{\infty }{e}^{s\lambda y^2}{e}^{-y^2/2}\frac{dy}{\sqrt{2\pi }}=(1-2s\lambda {)}^{-1/2}.$$

$$\mathbb{E}\left[e^{s{\mathfrak{z}}^{T}A\mathfrak{z}}\right]=\prod _j{\left(1-2s{\lambda }_j\right)}^{-1/2}$$

$$\mathbb{E}[{\mathfrak{z}}^{T}A\mathfrak{z}]$$

$$\mathrm{var}[{\mathfrak{z}}^{T}A\mathfrak{z}]$$

$$g(s)=\log \mathbb{E}\left[e^{sX}\right]=\frac{1}{2}\sum _j-\log (1-2s{\lambda }_j)$$

$$\mathbb{E}[X{e}^{sX}]=\mathbb{E}[e^{sX}]\frac{1}{2}\sum _j\frac{2{\lambda }_j}{1-2s{\lambda }_j}$$

$$\mathbb{E}[X^2{e}^{sX}]=\mathbb{E}[X{e}^{sX}]\frac{1}{2}\sum _j\frac{2{\lambda }_j}{1-2s{\lambda }_j}+\mathbb{E}[e^{sX}]\frac{1}{2}\sum _j\frac{(2{\lambda }_j{)}^2}{(1-2s{\lambda }_j{)}^2}.$$

Full text

1. Introduction

The goal of this paper is to prove a central limit theorem for $\int_{B}\phi^{2}$, where $\phi:M\rightarrow\mathbb{R}$ is a random wave and the ball $B$ may shrink with the wavelength of $\phi$. On any compact manifold $M$ with a Riemannian metric and corresponding Laplace operator $\Delta$, the random functions we have in mind are given by

[TABLE]

where the eigenfunctions $\phi_{j}$ solve $(\Delta+t_{j}^{2})\phi_{j}=0$ with eigenvalues $t_{j}^{2}$ in a window $T-\eta(T)\leq t_{j}\leq T$. The coefficients $c_{j}$ are independent Gaussians of mean 0 and identical variance. The choice of variance will disappear when we pass to the standardized random variable

[TABLE]

which is invariant under scaling $\phi$ by a constant multiple or, what is the same, scaling the variance of the coefficients. Nevertheless, it is natural to choose the variance of each coefficient $c_{j}$ inversely proportional to the number of terms in the sum (1.1), so that $\mathbb{E}[\int_{M}\phi^{2}]=1$. A standard Gaussian $G$ has characteristic function

[TABLE]

and our main result is an estimate comparing $\mathbb{E}[e^{itZ}]$ to this.

Theorem 1.1**.**

For any compact Riemannian manifold $M$ of dimension $n\geq 2$, any window $\eta(T)\rightarrow\infty$ with $\eta(T)\max(T^{-1/2},r^{2}T)$ bounded above, and a geodesic ball $B=B_{r}(z)$ shrinking at such a rate that $rT\rightarrow\infty$ and $(rT)^{n-1}/\eta\rightarrow 0$, we have convergence of characteristic functions

[TABLE]

for any fixed parameter $t\in\mathbb{C}$.

The characteristic function \mathbb{E}\big{[}e^{itZ}\big{]} is defined for all real $t$ and, once $rT$ is sufficiently large, for any fixed complex $t$. In particular, Theorem 1.1 implies convergence in distribution.

Corollary 1.2**.**

(Central Limit Theorem) In the limit that $T\rightarrow\infty$ with $r\rightarrow 0$, $\eta\rightarrow\infty$ as above, the standardized mean squares $Z_{B}$ converge in distribution to a Gaussian of mean 0 and variance 1:

[TABLE]

For example, one could take $\eta=\log{T}$ and $r=\log(T)^{\varepsilon}/T$ for a sufficiently small $\varepsilon>0$. It is important for our proof that $r\rightarrow 0$, although one might expect the Central Limit Theorem to apply also for sufficiently small fixed radii in the limit $T\rightarrow\infty$. The scale $1/T$ is significant because it is the natural wavelength for Laplace eigenfunctions of frequency $T$, and hence for the random waves (1.1). We assume that $rT\rightarrow\infty$ so that a length scale $r$ is enough to contain many oscillations of $\phi$, and we have in mind that $rT\rightarrow\infty$ arbitrarily slowly so that one is almost at the wave scale. We do not expect the central limit theorem to hold if $rT$ remains bounded, as we explain in Section 10.3, so that the wave scale provides a natural barrier.

We now give an overview of our strategy for proving Theorem 1.1 and some of the notation we will employ. Throughout the paper, we write $\sum_{j}$ as an abbreviation for sums over a given window of frequencies $T-\eta(T)\leq t_{j}\leq T$. For example, we have already done so in (1.1). We use the notations

[TABLE]

with the same meaning: there is a constant $C>0$ such that $|f(X)|\leq Cg(X)$ for all sufficiently large values of the variable(s) $X$. We write $f(X)\asymp g(X)$ if both of the approximate inequalities $f(X)\lesssim g(X)$ and $g(X)\lesssim f(X)$ hold, that is, there are constants $C_{1},C_{2}>0$ such that

[TABLE]

If the constants $C,C_{1},C_{2},\ldots$ depend on parameters, we indicate this with subscripts, as in

[TABLE]

The function $f$ is not necessarily real when we write $f(X)\lesssim g(X)$, but we only write $f(X)\asymp g(X)$ when both functions are positive. In Section 2, we compute the characteristic function of a quadratic form in Gaussian random variables, which includes $\int_{B}\phi^{2}$ as a special case. We give a sufficient condition for such a quadratic form to obey the Central Limit Theorem.

Proposition 1.3**.**

Given a sequence of symmetric, positive definite $N\times N$ matrices $A=A_{N}$ and standard Gaussian vectors $\mathfrak{z}=(\mathfrak{z}_{1},\ldots,\mathfrak{z}_{N})$, let $Z$ be the standardized quadratic form

[TABLE]

If

[TABLE]

as $N\rightarrow\infty$, then the characteristic function of $Z$ obeys

[TABLE]

This can be seen as an instance of Lyapunov’s criterion for deriving a central limit theorem [5, p.362]: If $S=X_{1}+\ldots X_{N}$ is a sum of independent random variables, perhaps not with the same distribution, and $\mathbb{E}[|S|^{2+\delta}]$ is small compared to $\mathbb{E}[|S|^{2}]^{1+\delta/2}$ for some $\delta>0$, then $S$ converges to a Gaussian as $N\rightarrow\infty$. Proposition 1.3 takes $\delta=1$, and one diagonalizes the quadratic form to obtain a sum of independent random variables. For a quadratic form in Gaussians, one can replace the general argument of Lyapunov’s criterion with an explicit computation of the characteristic function.

In the case of $\int_{B}\phi^{2}$, the matrix $A$ has a special form that allows the traces $\operatorname{tr}(A^{p})$ to be computed in terms of the two-point function

[TABLE]

We have $\mathbb{E}[\phi(x)\phi(y)]=\operatorname{var}[c]K(x,y)$ so that, after choosing the variance of the coefficients, $K(x,y)$ gives the correlations between values of $\phi$ at different points. In Section 3, we reformulate the criterion of Proposition 1.3 as

Proposition 1.4**.**

If $\phi=\sum_{j}c_{j}\phi_{j}$ is a Gaussian random function with two-point function $K(x,x^{\prime})$ satisfying

[TABLE]

as the ball $B$ shrinks, then the local integral $\int_{B}\phi^{2}$ obeys the Central Limit Theorem in the form

[TABLE]

Theorem 1.1 is then deduced from the following two estimates:

Theorem 1.5**.**

In any dimension $n\geq 2$,

[TABLE]

Theorem 1.6**.**

In dimensions $n\geq 3$,

[TABLE]

In Section 4, we use the local Weyl law to approximate $K(x,y)$ and it is for this purpose that Theorem 1.1 assumes $\eta\lesssim T^{1/2}$. In Section 5, we prove Theorem 1.5, and in particular give a lower bound on the denominator $\int_{B}\int_{B}K^{2}$, which is essentially the variance of $\int_{B}\phi^{2}$. In Section 6, the hypothesis $\eta\lesssim r^{-2}T^{-1}$ comes into play as we approximate the triple integral in the numerator by a Euclidean version. For $n\geq 3$, we bound this integral from above in Sections 7 and 8, completing the proof of Theorem 1.6 in Section 9 for dimensions $n\geq 3$. We prove the two-dimensional case of Theorem 1.1 in Section 10.2 by a different argument: reducing to the case where $M$ is the round sphere.

For context, we recall Berry’s Random Wave Model [4], which uses monochromatic random waves of the form (1.1) as a stand-in to make predictions about non-random eigenfunctions. This is expected to be a good approximation for chaotic systems. The model applies to high-frequency Laplace eigenfunctions on a manifold of negative curvature and suggests that, at an appropriate scale, they should be uniformly distributed. Related to this is the Quantum Unique Ergodicity conjecture of Rudnick and Sarnak [24] that on a negatively curved manifold $M$,

[TABLE]

for any fixed measurable subset $A$ of $M$ and any sequence of Laplace eigenfunctions $\phi_{\lambda}$ with growing eigenvalue $\lambda\rightarrow\infty$. This is one of our motivations for considering local integrals $\int_{A}\phi^{2}$, and the limit $\operatorname{vol}(A)$ corresponds to the expected value of $\int_{A}\phi^{2}$ when $\phi$ is taken at random. The eigenvalue equation $(\Delta+\lambda)\phi=0$ imparts a quantum interpretation, $|\phi|^{2}$ being the probability density of a single quantum particle in $M$.

The quantum ergodicity theorem proved by Shnirelman [25, 26], Colin de Verdière [9], and Zelditch [29] shows that negative curvature implies convergence along a full subsequence of eigenfunctions, or equivalently on average over the eigenfunctions, but allows many other subsequential limits besides the uniform measure. Although it remains unknown whether the uniform measure is the only possibility, work of Anantharaman [1], Anantharaman-Nonnenmacher [2], Anantharaman-Silberman [7], and Dyatlov-Jin [12] places significant constraints on the measures that arise as quantum limits in general. In examples of arithmetic origin, QUE has been proved in work of Lindenstrauss [20, 21], and Bourgain-Lindenstrauss [6], Jakobson [19], Holowinsky [18], Holowinsky-Soundararajan [17].

The deterministic problem being so difficult, it is of interest to randomize the function $\phi$ and study the fluctuations in $\int_{A}\phi^{2}$. Our main result shows that for $A$ shrinking at a certain rate, one can expect central limit statistics.

The central limit theorem controls deviations of size $``\sqrt{\operatorname{var}}"$ in the sense that

[TABLE]

for any fixed $y$. We also obtain the following concentration estimate at other scales.

Theorem 1.7**.**

If $y(rT)^{-1/6}$ is sufficiently small in terms of the implicit constants in the previous theorems, then

[TABLE]

where $c>0$ is a numerical constant and $X=\int_{B}\phi^{2}$ as above. In particular, the conclusion holds if $y(rT)^{-1/6}\rightarrow 0$ as $rT\rightarrow\infty$.

Han and Tacy proved deviation bounds in a related setting in equations (4.4), (4.5) of [15], which was a motivation for the present work.

2. Characteristic function of a quadratic form in Gaussians

In this section, we prove Proposition 1.3. Our quantity of interest $Z$ is a quadratic form in the random Gaussian coefficients $c_{j}$, and therefore has an explicit characteristic function. Indeed, expanding the square in $\phi^{2}$ shows that

[TABLE]

Write $c_{j}=\sqrt{\operatorname{var}[c]}\mathfrak{z}_{j}$, where $\mathfrak{z}_{j}$ is a centered Gaussian of unit variance and $\operatorname{var}[c]$ denotes the common variance $\operatorname{var}[c_{j}]$ among all $j$. Writing $\mathfrak{z}$ for the vector with entries $\mathfrak{z}_{j}$, we have

[TABLE]

where the matrix $A$ has entries

[TABLE]

The moment generating function can be computed explicitly by diagonalizing this symmetric matrix $A$. Write $A=U^{T}DU$, so that

[TABLE]

where $\lambda_{j}$ are the eigenvalues of $A$. For an orthogonal matrix $U$ and a standard Gaussian vector $\mathfrak{z}$, the vector $y=U\mathfrak{z}$ also has a standard Gaussian distribution. For $s\geq 0$ small enough that $1-2s\lambda_{j}>0$ for all $j$, independence of the random variables $y_{j}$ implies that

[TABLE]

which is a product of Gaussian integrals. These can be evaluated as

[TABLE]

We summarize this as

Proposition 2.1**.**

If $A$ is a symmetric matrix with eigenvalues $\lambda_{j}\geq 0$, then the moment generating function of the quadratic form $\mathfrak{z}^{T}A\mathfrak{z}$ in Gaussian random variables $\mathfrak{z}_{j}$ is

[TABLE]

where both sides are defined for $s$ sufficiently small that $1-2s\lambda_{\max}>0$, $\lambda_{\max}$ being the largest eigenvalue of $A$.

Note that, by differentiating the moment generating function,

[TABLE]

Indeed, let $X=\mathfrak{z}^{T}A\mathfrak{z}$ and differentiate with respect to $s$ starting from

[TABLE]

to

[TABLE]

and then to

[TABLE]

Taking $s=0$ gives $\mathbb{E}[X]=\sum\lambda_{j}$ and then $\mathbb{E}[X^{2}]=(\mathbb{E}[X])^{2}+2\sum\lambda_{j}^{2}$. Thus, as claimed,

[TABLE]

Passing to the standardized random variable $Z$, write

[TABLE]

so that

[TABLE]

Taking logarithms,

[TABLE]

Expanding the logarithm $-\log(1-x)=\sum x^{p}/p$, the linear term cancels:

[TABLE]

Now we take $s=it$ and standardize. Shifting and scaling and starting from

[TABLE]

we find that the characteristic function of the standardized quantity

[TABLE]

is

[TABLE]

The factor $e^{-t^{2}/2}$ is the characteristic function of a standard Gaussian. To show that $Z$ converges to a Gaussian, we will estimate the traces $\operatorname{tr}(A^{p})$ and show that they are negligible compared to $\sigma^{p}=(2\operatorname{tr}(A^{2}))^{p/2}$.

We complete the proof of Proposition 1.3 by using the third moment to bound the higher ones. Since $p$-norms are monotone, we have an upper bound

[TABLE]

Therefore, summing a geometric series,

[TABLE]

where the series converges for any fixed $t$ as long as

[TABLE]

Here, as in the rest of the paper, we write $\lesssim$ to denote inequality up to a constant, and $\rightarrow$ refers to limits where $rT\rightarrow\infty$ as described in the Introduction. Thus an upper bound on the third moment and a lower bound on the second are enough to control all others. Combining this with (2.6), we obtain the following estimate for the characteristic function, as stated in the Introduction.

Proposition 1.3.

Given a sequence of symmetric, positive definite $N\times N$ matrices $A=A_{N}$, let $Z$ be the standardized quadratic form

[TABLE]

If

[TABLE]

then the characteristic function obeys

[TABLE]

The rest of our work is to bound the second and third moments.

3. Moments and the two-point function

We now show how Proposition 1.4 is a special case of Proposition 1.3. For any quadratic form $\mathfrak{z}^{T}A\mathfrak{z}$ and any $p\geq 2$, we have

[TABLE]

but for the particular one $\int_{B}\phi^{2}$, we can express the traces in terms of the kernel

[TABLE]

The matrix entries are $A_{jk}=\operatorname{var}[c]\int_{B}\phi_{j}\phi_{k}$. By writing the product of $p$ integrals over $B$ as a single integral over $B^{p}$, we obtain

[TABLE]

When we take the (finite) sum $\sum_{k}\sum_{j_{1}}\cdots\sum_{j_{p-1}}$ under the integral, it factors into $p$ copies of the two-point function:

[TABLE]

With the extra factor of $\operatorname{var}[c]$ from the entries $A_{jk}$, the trace becomes

[TABLE]

In particular,

[TABLE]

We summarize this discussion in the following restatement of Proposition 1.4, which may be helpful in other examples besides the monochromatic ensemble.

Proposition 3.1**.**

If $\phi=\sum_{j}c_{j}\phi_{j}$ is a Gaussian random function with two-point function satisfying

[TABLE]

as the ball $B$ shrinks, then the local integral $\int_{B}\phi^{2}$ obeys the Central Limit Theorem.

In our example, we will use semiclassics to estimate $K(x,x^{\prime})$ and verify the assumption that this ratio of integrals tends to 0 as $rT\rightarrow\infty$. The lower bound on the denominator is provided by

Theorem 1.5.

[TABLE]

while the upper bound on the numerator is

Theorem 1.6.

[TABLE]

4. Weyl’s law

The local Weyl law is an estimate for $K(x,x^{\prime})$ showing that it resembles a Bessel function at the scale $1/T$. Such results were obtained by Hörmander in the regime $d(x,y)\lesssim T^{-1}$, but we are interested in distances $d(x,y)$ shrinking arbitrarily slowly as $T\rightarrow\infty$. In this situation, we have the following by work of Canzani and Hanin.

Theorem 4.1**.**

(Canzani-Hanin) Let $M$ be a compact manifold of dimension $n\geq 2$, with smooth Riemannian metric $g$ and corresponding Laplace eigenfunctions obeying $(\Delta+t_{j}^{2})\phi_{j}=0$. Fix any point $x_{0}\in M$. For a large parameter $T$, let $B=B_{r}(x_{0})$ where $r=r(T)\rightarrow 0$ arbitrarily slowly as $T\rightarrow\infty$. Then the spectral function

[TABLE]

can be written

[TABLE]

where the remainder satisfies

[TABLE]

Hörmander had already proved this with a remainder at most $O(T^{n-1})$ for $d(x,y)$ of order $1/T$. In Theorem 2 from [8], Canzani and Hanin improved this in two ways: allowing $d(x,y)$ to shrink more slowly than $1/T$ and, assuming $x_{0}$ is not self-focal, achieving an error term $o(T^{n-1})$ instead of $O(T^{n-1})$. Without the assumption on $x_{0}$, their proof still gives $O(T^{n-1})$ for $d(x,y)$ shrinking arbitrarily slowly, as we proceed to sketch. Proposition 18 in [8] is the main technical estimate, and it allows $x_{0}$ to be self-focal. The role of the non-self-focal assumption is to guarantee the estimate (43) in [8], which says roughly that there is a constant $c>0$ such that when $x=y$, the remainder is at most $c\varepsilon T^{n-1}+O_{\varepsilon}(T^{n-2})$ for any $\varepsilon>0$. Without the (non)focal assumption, one still knows by Hörmander’s work that the remainder is $O(T^{n-1})$ on the diagonal and the same arguments from [8] (in particular Propositions 10 and 11) imply that $R(x,y)$ remains $O(T^{n-1})$ for distances $d(x,y)$ shrinking arbitrarily slowly.

The phase $\langle\exp_{y}^{-1}(x),\xi\rangle_{g_{y}}$ is an important aspect of [8]. It is defined as long as $d(x,y)$ is less than the injectivity radius of $M$, that is, even for distances of order 1. In contrast, Hörmander’s “adapted” phase functions $\psi(x,y,\xi)$ are obtained from the eikonal equation, a differential equation which is not guaranteed to have solutions when $Td(x,y)$ is unbounded.

For a window of frequencies $T-\eta<t_{j}\leq T$, the kernel is

[TABLE]

The remainders $R_{T-\eta}\lesssim(T-\eta)^{n-1}$ and $R_{T}\lesssim T^{n-1}$ are both $O(T^{n-1})$. Therefore Theorem 4.1 (Theorem 2 from [8]) implies an estimate for $K(x,y)$. Canzani-Hanin state the result as Theorem 3 in [8], but again they are interested in obtaining an improved remainder $o(T^{n-1})$ in the most crisp regime $\eta=1$ under the assumption that $x,y$ are close to a non-self-focal point $x_{0}$. In the present article, we avoid making this assumption at the price of taking a slowly growing window $\eta\rightarrow\infty$. Rescaling the integral for $E_{T-\eta}$ by $\xi\mapsto T\xi/(T-\eta)$, we have

[TABLE]

where we have suppressed the subscripts with the understanding that $|\xi|$ and $\langle\exp_{y}^{-1}(x),\xi\rangle$ are taken with respect to the metric at $y$. With the multiplicative factor $(\frac{T-\eta}{2\pi})^{n}$,

[TABLE]

It follows that

[TABLE]

We assume $\eta/T\rightarrow 0$, so the integration over the thin annulus $1-\eta/T<|\xi|<1$ can be appproximated by an integral over the sphere $|\xi|=1$ times the width $\eta/T$. In polar coordinates at $y$, $g_{y}$ is the identity matrix. Identifying $x$ and $y$ with the vectors in $\mathbb{R}^{n}$ representing them in coordinates,

[TABLE]

Finally, the Bessel function appears because of the integral formula

[TABLE]

We summarize this as follows, writing $d(x,y)$ for the distance between $x$ and $y$.

Proposition 4.2**.**

For any compact manifold $M$ and radius $r=r(T)\rightarrow 0$ as $T\rightarrow\infty$, the two-point function for $d(x,y)<r$ is given approximately by

[TABLE]

where $N$ is the number of frequencies $t_{j}$ in the interval $T-\eta<t_{j}\leq T$. We have

[TABLE]

for a positive constant $c_{n}>0$ depending only on the dimension.

We will often use Proposition 4.2 as an upper bound. The Bessel term is bounded, so we always have $K(x,y)\lesssim N\lesssim T^{n-1}\eta$. This can be improved when $d(x,y)$ is larger than $1/T$ since $J_{n/2-1}(u)$ decays as $u^{-1/2}$:

[TABLE]

To explain the normalization in Proposition 4.2, we review the local and global versions of Weyl’s law. Weyl’s law for counting eigenvalues can be deduced from the local version by integrating on the diagonal $x=y$. Indeed, for orthogonal eigenfunctions $\phi_{j}$, we have

[TABLE]

This leads to the usual Weyl’s law that the number of eigenvalues in our window is

[TABLE]

where $|B_{n}|$ is the volume of the Euclidean unit ball in dimension $n$. Recall that we write the Laplace eigenvalue as $t_{j}^{2}$ (largely to avoid confusion with the eigenvalues of $A$, already called $\lambda_{j}$). From the binomial expansion,

[TABLE]

If we choose $\eta(T)=T^{1/2}$ or smaller, then $T^{n-2}\eta(T)^{2}\leq T^{n-1}$ and we can absorb the higher order terms into the error $O(T^{n-1})$ already present in Weyl’s law. We need $T^{n-1}\eta(T)$ to be larger than $T^{n-1}$ in order for the main term to dominate. Thus for $\eta(T)$ diverging no faster than $T^{1/2}$, we have

[TABLE]

The variance of the coefficients should be inversely proportional to this in order for $\phi$ to be normalized in $L^{2}$ (on average). Indeed, for orthonormal eigenfunctions $\phi_{j}$, we have

[TABLE]

It is also of interest to study windows $\eta$ even larger than $T^{1/2}$. In this case, the error $O(T^{n-2}\eta^{2})$ from the binomial expansion overwhelms the error $O(T^{n-1})$ from Weyl’s law. If $\eta=(1-\alpha)T$ is proportional to $T$ instead of only $o(T)$, then there is no need to use the binomial formula and we instead have

[TABLE]

This is the case of band-limited random functions, synthesized from a band of frequencies between $(1-\alpha)T$ and $T$.

5. Proof of Theorem 1.5

We have Weyl’s law

[TABLE]

and hence

[TABLE]

We write $J$ as a shorthand for the Bessel function

[TABLE]

We integrate over $B\times B$ to get

[TABLE]

By the triangle inequality, the ball $B_{r}(z)$ contains a ball $B_{r-\xi}(x)$ with smaller radius and center shifted to $x$, where $\xi=d(x,z)$. The integrand $J^{2}$ is nonnegative, so we have a lower bound by integrating only over the smaller ball. Shifting to polar coordinates with respect to $x$, we have

[TABLE]

Indeed, the volume form in normal coordinates is given by

[TABLE]

since the volume form is obtained from the metric $g$ by $\sqrt{\det(g)}$ and we have the expansion

[TABLE]

We have

[TABLE]

Again, we achieve a lower bound by integrating only over a subset $B_{r}(z)$, namely $B_{r/2}(z)$. For $x$ in this smaller ball, we have $r-\xi\geq r/2$. Thus

[TABLE]

Change variables to $u=T\rho$ in the inner integral, so that $\rho d\rho=T^{-2}udu$. Then

[TABLE]

We have $\operatorname{vol}(B_{r/2})\asymp\operatorname{vol}(B_{r}(z))\asymp r^{n}$. Also, the function $uJ_{n/2-1}(u)^{2}$ is bounded because of the $1/\sqrt{u}$ decay of $J$-Bessel functions. It has a non-zero average value

[TABLE]

Hence

[TABLE]

In terms of $K$, this means

[TABLE]

To have a main term, we must assume that

[TABLE]

We have

[TABLE]

so the condition is simply

[TABLE]

which is one of the assumptions in Theorem 1.1.

One can obtain a comparable upper bound by a similar argument. Although we only need a one-sided bound to compare $\operatorname{tr}(A^{3})$ and $\operatorname{tr}(A^{2})$, it is worth knowing that the lower bound above is essentially of the correct order of magnitude. The ball $B_{r}(z)$ is contained in the re-centered ball $B_{2r}(x)$, so that

[TABLE]

As before, we shift to polar coordinates centered at $x$. Weyl’s law gives a pointwise upper bound on $K(x,x^{\prime})$:

[TABLE]

This diverges as $\rho\rightarrow 0$, since we would be better off using the trivial bound $K\lesssim N$ for $\rho<1/T$, but the singularity is integrable. We obtain

[TABLE]

Integrating over $x$ brings another factor of $\operatorname{vol}(B_{r})\asymp r^{n}$. The resulting upper bound for the integral $\int_{B}\int_{B}K(x,x^{\prime})^{2}dxdx^{\prime}$ is

[TABLE]

Under our assumption, $\eta$ is large compared to $(rT)^{n-1}$ so that the terms involving $\eta^{-1}$ or $\eta^{-2}$ are even smaller than $(rT)^{-n+1}$. Therefore the upper and lower bounds match up to a constant multiple:

[TABLE]

6. Euclidean integral for the third moment

We have stated our main theorem for any compact Riemannian manifold. This level of generality is possible because the radius of the ball $B$ is shrinking and the triple integral $\int_{B}\int_{B}\int_{B}K_{12}K_{23}K_{31}$ is therefore close to its Euclidean counterpart. Recall that we write $K_{12}=K(x_{1},x_{3}),K_{23}=K(x_{2},x_{3}),K_{31}(x_{3},x_{1})$ and similarly $J_{12},J_{23}J_{31}$ for their approximations by Bessel functions. Write Weyl’s law in the form

[TABLE]

where $J$ is the Bessel function from (5.1). It follows that

[TABLE]

Write $|x-y|$ for Euclidean distance and $d(x,y)$ for Riemannian distance between the points represented by $x$ and $y$ in coordinates. Because $J$ has bounded derivative,

[TABLE]

In a ball of radius $r\rightarrow 0$, the relative error in approximating Riemannian distances by Euclidean ones is a factor of $1+O(r)$. Hence

[TABLE]

Let $\xi=d(x_{1},z)$, $\rho=d(x_{1},x_{3})$, $s=d(x_{1},x_{2})$. We compare the integral to a Euclidean version, introducing a small error in order to pass to polar coordinates based at $x_{1}$. Take polar coordinates $(\rho,\omega)$ where the radial coordinate $\rho=d(x,x_{1})$ ranges from 0 to $2r$. As in Section 5, the volume form is given approximately by its Euclidean counterpart, namely

[TABLE]

where the error depends on the curvature of $M$. We integrate and conclude that

Proposition 6.1**.**

If $B=B_{r}(z)$ is a ball of radius $r\rightarrow 0$ in a compact manifold $M$ with a smooth Riemannian metric, then

[TABLE]

where $B_{0}$ is a Euclidean ball of radius $r$ centered at the coordinates for the basepoint $z\in M$ (or equally well, centered anywhere in $\mathbb{R}^{n}$ since the Euclidean integral is translation-invariant).

For the original kernel $K$, larger than $J$ by a factor of $T^{n-1}\eta$, this leads to an error $(\eta T^{n-1}\operatorname{vol}(B))^{3}Tr^{2}$. We assume $\eta\lesssim T/(rT)^{2}$ so that this can be absorbed into the error term $(\eta T^{n-1}\operatorname{vol}(B))^{3}\eta^{-1}$ already with us from Weyl’s law. Thus we may proceed with the Euclidean integral.

Let us first see what upper bound follows from crude estimates ignoring cancellation in the Bessel integrand. In higher dimensions, the angular variables $\alpha$ and $\beta$ range over $S^{n-1}$ instead of $S^{1}$. With $x_{1}$ as the origin in Euclidean space, write

[TABLE]

Then

[TABLE]

We may choose the coordinates to make $\zeta$ the north pole $(1,0,\ldots,0)$ and then account for the integration over $\zeta$ with an overall factor of $\operatorname{vol}(S^{n-1})$. The distance $\xi=d(z,x_{1})$ ranges from 0 to $r$. The distances $\rho=d(x_{1},x_{3})$ and $s=d(x_{1},x_{2})$ may be as large as $r+\xi$, with an angular range starting from a full sphere $S^{n-1}$ for distances less than $r-\xi$ and shrinking to nothing as the distance approaches $r+\xi$. More specifically, the condition $x_{2}\in B_{r}(z)$ is

[TABLE]

and similarly for $x_{3}$ with its angle $\beta$. Thus the angles $\alpha$ and $\beta$ obey

[TABLE]

and we may simply take the dot products $\alpha\cdot\zeta$ and $\beta\cdot\zeta$ as the first coordinates $\alpha_{1}$ and $\beta_{1}$, if we commit to the north pole. Let us write $\mathfrak{a}$ for the range of $\alpha$ defined by (6.2), and likewise $\mathfrak{b}$ for the range of $\beta$.

The integral is, up to a factor $\operatorname{vol}(S^{n-1})$,

[TABLE]

where we write $J(y)=c_{n}J_{n/2-1}(y)/y^{n/2-1}$ for the Bessel function. Change variables to $w=\xi/r$, $u=Ts$, $v=T\rho$, with change of measure

[TABLE]

The angles $\alpha$ and $\beta$ are unaffected by this scaling. The integral becomes

[TABLE]

Think of the prefactor $r^{n}T^{-2n}$ as $\operatorname{vol}(B)^{3}(rT)^{-2n}$, since $\operatorname{vol}(B)$ is comparable to $r^{n}$. We have

[TABLE]

which oscillates between values of order $u^{(n-1)/2}$. Simply observing that

[TABLE]

and ignoring cancellation gives an upper bound

[TABLE]

The exponent $n-1$ must be improved to something larger than $3(n-1)/2$ in order to beat the variance, since the denominator is $\operatorname{tr}(A^{2})^{3/2}\asymp(rT)^{-3(n-1)/2}$.

7. Gegenbauer’s addition formula

To handle the term $J(|u\alpha-v\beta|)$ in the inner integral, we have recourse to

Fact 7.1**.**

(Gegenbauer’s addition formula, [22, formula 10.23.8], [28, p.362-368]) If $\varpi=\sqrt{u^{2}+v^{2}-2uv\cos{\theta}}$, then for $\nu\neq 0,-1,-2,\ldots$,

[TABLE]

where $C_{k}^{\nu}$ are the Gegenbauer polynomials.

Recall that the polynomials $C_{k}^{\nu}(t)$ are orthogonal over the interval $[-1,1]$ with weight $(1-t^{2})^{\nu-1/2}$. Their maximum value is attained at

[TABLE]

The notation $(2\nu)_{k}$ is shorthand for $2\nu(2\nu+1)\cdots(2\nu+k-1)$. Note that when $\nu=1/2$ the Gegenbauer polynomials are bounded by 1, no matter how high the degree. For our application, $\nu=n/2-1$, so this bounded case arises when $n=3$. In higher dimensions $n>3$, $C_{k}^{\nu}$ at its largest grows as a power of $k$. Nevertheless, the series in Fact 7.1 converges because of the rapid decay in $J_{\nu+k}(u)J_{\nu+k}(v)$ as a function of $k$.

To bring Gegenbauer’s addition formula to bear on our inner integral, note that

[TABLE]

so we may take $\theta=\arccos(\alpha\cdot\beta)$. Write $\mathfrak{a}=\mathfrak{a}(u,v)$ and $\mathfrak{b}=\mathfrak{b}(u,v)$ for the region of integration over the angles $\alpha$ and $\beta$ from (6.2), suppressing the dependence on the outermost variables $w$ and $\zeta$. The integral is

[TABLE]

Next we integrate over $u$ and $v$. This can be interchanged with the sum over $k$, again thanks to the rapid decay of $J_{\nu+k}$ as $k$ grows. This leaves us with

[TABLE]

For many ranges of the parameters, $\mathfrak{a}(u,v)$ or $\mathfrak{b}(u,v)$ is a full sphere $S^{n-1}$, in which case $\int_{\mathfrak{a}}\int_{\mathfrak{b}}C_{k}^{\nu}(\alpha\cdot\beta)$ is exactly 0 by orthogonality of the Gegenbauer polynomials. In general, we make the following claim.

Proposition 7.2**.**

For $n\geq 3$, $\nu=n/2-1$, concentric arcs $\mathfrak{a}$ and $\mathfrak{b}$ as above, and degree $k\geq 1$

[TABLE]

Proof.

The angular ranges $\mathfrak{a}$ and $\mathfrak{b}$ are concentric spherical caps (centered at $\zeta$, the direction from $x_{1}$ to $z$). We may assume that $\mathfrak{a}\subseteq\mathfrak{b}$, reversing the order of $\alpha,\beta$ if need be. We integrate with respect to $\beta$ using polar coordinates $(\theta,\omega)$ based at $\alpha$. Thus $\cos{\theta}=\alpha\cdot\beta$ assumes some interval of values depending on the distance between $\alpha$ and $\zeta$, and also on the angle $\omega$. Because $\mathfrak{a}\subseteq\mathfrak{b}$, this interval contains $\theta=0$ since $\beta=\alpha$ automatically lies in the range $\mathfrak{b}$. The variable $\theta$ ranges over an interval $0\leq\theta\leq\theta_{+}$ where $\theta_{+}=\theta_{+}(\alpha,\omega)$, or this interval together with $\pi-\theta_{+}\leq\theta\leq\pi$ if $\mathfrak{a}$ and $\mathfrak{b}$ exceed a hemisphere. This only changes $\cos{\theta}$ by a sign and the argument below will apply equally well. Depending on the angle $\omega$, $\theta_{+}$ lies somewhere between the extreme cases ${\rm diam}(\mathfrak{b})/2-d(\alpha,\zeta)$ and ${\rm diam}(\mathfrak{b})/2+d(\alpha,\zeta)$. Note also that $\theta_{+}\leq\pi/2$ or else the hemispheres would merge. We write

[TABLE]

where $(\sin\theta)^{n-2}d\theta d\omega$ is the volume element on $S^{n-1}$ and the polynomials $C_{k}^{\nu}(\cos\theta)$ are orthogonal on the interval $0\leq\theta\leq\pi$ with respect to the weight $\sin(\theta)^{n-2}$. For the integral over $\theta$, note that with $x=\cos{\theta}$ and $\nu=n/2-1$,

[TABLE]

We have the following exact antiderivative

[TABLE]

This is a special case of formula (18.17.1) in [22], which gives the integral of Jacobi polynomials:

[TABLE]

The conversion between the ultraspherical notation $C_{k}^{\nu}$ and the Jacobi polynomial is

[TABLE]

where $(t)_{k}=t(t+1)\cdots(t+k-1)$ denotes a rising factorial. We have

[TABLE]

Converting between $C_{k}^{\nu}$ and $P_{k}^{(\nu-1/2,\nu-1/2)}$ under the integral, as well as $P_{k-1}^{(\nu+1/2,\nu+1/2)}$ and $C_{k-1}^{\nu+1}$ outside it, therefore brings a factor

[TABLE]

as stated in (7.3).

For the definite integral, note that when $\theta=0$, we have $1-x^{2}=0$, so that only the upper endpoint $\theta_{+}$ contributes. It follows that

[TABLE]

We recall Szegő’s asymptotic formula (8.21.17) in [27]: For large degree $k$, the Gegenbauer polynomial is given approximately by

[TABLE]

up to an additive error of $O(\theta^{1/2}k^{-\nu-1})$, and even smaller if $\theta\lesssim 1/k$. In particular, applying this with $\nu+1$ in place of $\nu$ gives

[TABLE]

We have already noted that $C_{k}^{\nu}(1)\lesssim k^{n-3}$, and changing $\nu=n/2-1$ to $\nu+1$ corresponds to increasing $n$ by 2. Thus

[TABLE]

The large-argument bound $J_{\nu+1/2}(y)\lesssim y^{-1/2}$ then gives

[TABLE]

We have used the fact that $n\geq 3$ to bound $\theta_{+}^{\nu}\lesssim 1$, the exponent being nonnegative. Taking account of the prefactor

[TABLE]

it follows that

[TABLE]

Integrating over $\alpha$ and $\omega$ brings only a bounded factor (volumes of spheres, depending only on $n$) and the claim follows. ∎

Using Proposition 7.2, the upper bound for the third moment becomes

[TABLE]

8. Beginning, middle, and end

In this section, we estimate the following one-variable integral arising from Section 7.

Proposition 8.1**.**

As $rT\rightarrow\infty$ with $m\leq 2rT$,

[TABLE]

whereas if $m>2rT$, then

[TABLE]

for any $b<1$, with $c(b)>0$ a correspondingly small constant. Both cases (8.1) and (8.2) are uniform with respect to $0\leq w\leq 1$ (or any bounded interval of $w$).

The proposition is easy to see for fixed $m$, because the integrand $uJ_{\nu}(u)J_{m}(u)$ is a bounded function of $u$ and the interval of integration has length at most $2rT$. The difficulty lies in allowing $m$ to grow with $rT$, but we will see that the same reasoning applies after subdividing the integral depending on the size of $m$. For our application, $m=k+\nu$ as before.

Fix a $\delta>0$ and let $m_{\pm}=m\pm m^{1/3+\delta}$. Then split the integral as

[TABLE]

where the latter two ranges are empty if $m$ is large enough that $m_{-}\geq(1+w)rT$. These ranges are chosen in light of the fact that $J_{m}(u)$ is very small when $u$ is much less than $m$, mildly decaying and oscillatory when $u$ is much larger than $m$, and with a transition between these regimes taking place when $|u-m|\lesssim m^{1/3}$. More precisely,

Proposition 8.2**.**

For $0<x\leq 1$ and $m\geq 0$ not necessarily an integer,

[TABLE]

For $u>m$,

[TABLE]

For any $u$, but of particular importance when $|u-m|<m^{1/3}$,

[TABLE]

The first of these is inequality (10.14.15) in [22]. The latter two bounds are crude consequences of much more detailed asymptotics available for $J_{m}(u)$. See, for instance, (10.20.4) of [22] or Chapter 7.4 of [13] for statements. For proofs, see chapter 8 of Watson’s book [28] or Olver’s article [23]. Note that the bounds are compatible in the sense that (8.4) recovers (8.5) when $u-m$ is of order $m^{1/3}$. Note also that (8.4) reduces to $J_{m}(u)\lesssim u^{-1/2}$ when $u$ is much larger than $m$, for instance if $m$ is fixed.

In the initial range $\int_{0}^{m_{-}}$, we apply (8.3) with $x\leq 1-m^{-2/3+\delta}$. First, we use

[TABLE]

to expand the upper bound

[TABLE]

as a power series in $\sqrt{1-x^{2}}$. Note the series expansion

[TABLE]

with $y=\sqrt{1-x^{2}}$ for our application. It follows that

[TABLE]

In the range $u\leq m_{-}$, we have $x\leq 1-m^{-2/3+\delta}$. Thus

[TABLE]

Now we can use $J_{\nu}(u)\lesssim u^{-1/2}$ to bound the initial contribution by an elementary integral:

[TABLE]

for a constant $c>0$.

For the middle contribution, we use the transitional bound $J_{m}(u)\lesssim m^{-1/3}$ together with $J_{\nu}(u)\lesssim u^{-1/2}$ for fixed $\nu$ to conclude that

[TABLE]

For the final contribution, we see from (8.4) that

[TABLE]

Changing variables to $x=u/m$ gives

[TABLE]

Because $x^{1/2}(x^{2}-1)^{-1/4}$ is integrable near $x=1$ and bounded for $x$ away from 1, we have

[TABLE]

If $m>2rT$, then for any value of $w$, only the initial segment $\int_{0}^{m_{-}}$ is present and so the integral is exponentially small, as claimed in (8.2). If $m\leq 2rT$, then in principle there may also be contributions as large as $m(1+rT/m)$ and $2m^{1/2+\delta}$, both of which are dominated by $rT$ for this range of $m$. Thus the integral is of the order $rT$ as claimed in (8.1). This proves Proposition 8.1.

9. Completing the proof

In Section 7, we used Gegenbauer’s addition formula to express the third moment as a sum obeying the upper bound

[TABLE]

up to a constant multiple. In Section 8, we estimated the integral appearing in the summands for each $k\geq 0$. In this section, we estimate the sum over $k$ and deduce that the third moment compares favourably to the variance. We write $m=\nu+k$.

We distinguish a bulk contribution, where $m\leq 2rT$, from a tail contribution where $m>2rT$. For $m$ in the bulk, combining the three ranges of integration gives the upper bound from (8.1):

[TABLE]

For $m$ in the tail, only the beginning of the integral is present:

[TABLE]

Since $m^{n/2-5/2}\exp(-cm^{3\delta/2})$ is summble, the tail contributes only a bounded amount in addition to the bulk sum:

[TABLE]

The upper bound on the third moment becomes

[TABLE]

This compares favourably to the variance:

[TABLE]

which converges to [math].

For any (positive definite) quadratic form in Gaussians, Proposition 1.3 gives

[TABLE]

We have shown that for our particular form $\int_{B}\phi^{2}$ in any dimension $n\geq 3$, the second and third moments obey

[TABLE]

This implies the central limit theorem in the form

[TABLE]

as $rT\rightarrow\infty$. Thus we have proved the central limit theorem for $n\geq 3$ using Gegenbauer’s addition formula.

10. Two-dimensional case

In dimension 2, we have $\nu=n/2-1=0$ so that Gegenbauer’s addition formula is not available and we must argue differently. Note that the variance bound is the same for $n=2$ as for higher $n$, namely $(rT)^{-n+1}$ or just $(rT)^{-1}$ in this case. A significant difference is that the exponent $\nu-1/2=(n-3)/2$ on the weight $(1-t^{2})^{\nu-1/2}$ becomes negative. The kernel $J_{0}(u)$ decays only as $u^{-1/2}$ instead of saving an additional power $u^{-(n/2-1)}$. The Gegenbauer polynomials reduce to Chebyshev polynomials $T_{n}(\cos{\theta})=\cos(n\theta)$.

As $r\rightarrow 0$, the local Weyl law leads to the same triple integral in Section 6 no matter the global geometry of $M$. As a result, if a direct calculation proves the limit theorem on any particular manifold, one can draw the same conclusion in general. Such a calculation is possible on the sphere, and we will deduce the two-dimensional case from the instance $M=S^{2}$. First we review the theory of spherical harmonics, including the case of higher dimensions. In any dimension, one could in principle use the same strategy of reducing to the sphere, but the two-dimensional case is simplest.

Note that the window $\eta$ plays a minor role. Under the conditions of Theorem 1.1, the central limit theorem reduces to bounding the Euclidean multiple integrals $\int\int\int K_{12}K_{23}K_{31}$ and $\int\int K^{2}$. These same integrals are obtained from the sphere with $\eta$ small enough that the window $T-\eta<t_{j}\leq T$ amounts to a single degree of spherical harmonics. Thus the general case of any manifold $M$ with a window as in Theorem 1.1 follows from the case of $M=S^{2}$ with, say, $\eta=1/2$.

10.1. Spherical harmonics

The standard basis of spherical harmonics of degree $l$ on $S^{p+1}$ is parametrized by integers

[TABLE]

whose sum is the total degree $\sum_{j}m_{j}=l$. For any such vector $m=(m_{0},\ldots,m_{p},\pm)$ together with a choice of sign $\pm$, one has a spherical harmonic

[TABLE]

where $C_{n}^{\lambda}$ denotes the Gegenbauer polynomial. As $m$ varies, these harmonics are orthonormal and span the space of all harmonics of degree $l$. For real-valued harmonics, it suffices to replace $e^{\pm im_{p}\varphi}$ by $\sin(m_{p}\varphi)$ and $\cos(m_{p}\varphi)$ (up to a $\sqrt{2}$ factor of normalization). These functions come from separation of variables in the coordinates

[TABLE]

on $\mathbb{R}^{p+2}$, where the angles $\theta_{j}$ range from [math] to $\pi$ while the final angle $\varphi$ ranges over a full circle $0\leq\varphi<2\pi$. We refer to [13], section 11.2 and especially formula 11.2(23).

Fix any point of $S^{p+1}$ as origin, say $(1,0,\ldots,0)$. The distance to the origin is then given by $\cos{\theta_{1}}$. We separate variables into a radial part and an angular part varying over the lower-dimensional sphere $S^{p}$. Thus

[TABLE]

where

[TABLE]

and the remaining factors are such that $Y_{m}$ is a spherical harmonic of degree $l-m_{0}$ on $S^{p}$. Note that

[TABLE]

and the corresponding coordinates on the lower-dimensional sphere $S^{p}$ are simply scaled by $\sin{\theta_{1}}$ in order to get $x_{2},\ldots x_{p+2}$. The spherical volume element on $S^{p+1}$ is

[TABLE]

where $d\alpha$ is the spherical volume element on $S^{p}$. In particular, two harmonics $f_{m}$ and $f_{m^{\prime}}$ with vectors $m\neq m^{\prime}$ are orthogonal on any ball $\{\theta_{1}<r\}$ because the angular parts are already orthogonal over $S^{p}$. Note that the condition $\sum_{j}m_{j}=l$ rules out the possibility that $f_{m}$ and $f_{m^{\prime}}$ agree on their angular coefficients but have $m_{0}\neq m^{\prime}_{0}$. Such harmonics would be of different degree.

For the two-dimensional sphere $S^{2}$, we have $p=1$. The coordinates are just $\theta_{1}$ and $\varphi$. Because of the constraint $m_{0}+m_{1}=l$, there is essentially just one component $m_{0}$ together with a choice of sign $\pm$. Thus we have the usual basis harmonics $Y_{l}^{m}$

[TABLE]

Now expand any given harmonic $\phi$ in the basis functions $f_{m}$ with respect to a chosen origin:

[TABLE]

where $m=(m_{0},\ldots,m_{p},\pm)$ varies over all vectors as above, that is, $\sum m_{j}=l$ and the entries of $m$ are decreasing from $m_{0}$ to $m_{p}$. We also specify a trigonometric factor $e^{\pm im_{p}\varphi}$ via the choice of sign. By orthogonality, integrating $f^{2}$ over a ball around the origin gives a diagonal quadratic form in the coefficients $c_{m}$:

[TABLE]

For our application, the functions $f_{m}$ must be normalized so that $\int_{S^{p+1}}f_{m}^{2}=1$. Thus

[TABLE]

since the angular integrals over $S^{p}$ are the same whether we integrate over the ball of radius $r$ or over the whole sphere $S^{p}$, that is, the ball of radius $\pi$. Note that the radial part $R_{m}$ depends only on $m_{0}$ and $m_{1}$. Given $m_{0}$ and $m_{1}$, the extra components $m_{2},\ldots,m_{p}$ and the sign $\pm$ weight the corresponding summand by a combinatorial factor. For example, if $m_{0}+m_{1}=l$, the weight is 2 because it must be that $m_{2}=\ldots m_{p}=0$, leaving only the two choices of sign. Likewise, if $m_{0}+m_{1}=l-1$, it must be that $m_{2}=1$ and the smaller components vanish, so the weight is again 2. But the larger the deficit between $m_{0}+m_{1}$ and the total degree $l$, the more combinations arise. A simplification available for $S^{2}$ compared to higher-dimensional spheres is that $m_{0},m_{1}$ are the only parameters and no further cases occur.

For large $n$, we have Szegő’s asymptotic for the Gegenbauer polynomial:

[TABLE]

where $N=n+a+1/2$. The error term is $O(\theta^{1/2}n^{-3/2})$, and even smaller if $\theta\lesssim 1/n$. See formula (8.21.17) in [27]. For us, $n=m_{0}-m_{1}$ and $a=m_{1}+(p-1)/2$. Thus

[TABLE]

10.2. Central limit theorem on the two-dimensional sphere

When $M=S^{2}$, the Laplace eigenvalues are $m(m+1)$ with multiplicity $2m+1$, where $m=0,1,2,\ldots$ is the degree (instead of $l$, as above in arbitrary dimension). If $T=\sqrt{m(m+1)}$ and $\eta$ is a small constant (for instance, $1/2$), then the window $T-\eta<t_{j}\leq T$ contains $2m+1$ copies of a single frequency equal to $\sqrt{m(m+1)}$. Note that $T=\sqrt{m(m+1)}\sim m+1/2$, so we may take $m$ rather than $T$ as the key asymptotic parameter. The eigenvalues of the quadratic form $\int_{B}\phi^{2}$ are given by

[TABLE]

This was used in [10] to estimate the $\lambda_{j}$ as follows. There are three cases, which we think of as “bulk”, “edge”, and “tail”.

Proposition 10.1**.**

Fix a point $z\in S^{2}$ and choose for our basis functions $\phi_{j}$ the standard ultraspherical polynomials rotated so that $z$ is at the North pole $(0,0,1)$. Then

[TABLE]

where the $\mathfrak{z}_{\nu}$ are independent standard Gaussians for $0\leq\nu\leq 2m$ and the coefficients $\lambda_{\nu}$ satisfy

[TABLE]

for any $a>0$ and a ratio $0\leq k/(rm)<1$ bounded away from $1$. For $k$ close to $rm$, we have an upper bound

[TABLE]

For $k$ so large that $k+k^{p}>rm$, where $p>1/3$, we have

[TABLE]

for a constant $c>0$.

Thus the eigenvalues of $A$ follow the semi-circle law. Perhaps one should expect some GOE behaviour because $A$ is given by integrals $\int_{B}\phi_{j}\phi_{k}$ and there is a large orthogonal group of symmetries acting by change of basis on the functions $\phi_{j}$.

As in Section 2, to show that $Z$ is approximately Gaussian, it is sufficient to show that for $p\geq 3$

[TABLE]

To verify (10.4), we can use the semicircle law from Proposition 10.1 to estimate $\sum\lambda_{j}^{p}$. First, we bound the contributions from $j$ near $rm$ or larger. Fix a value of $\delta>0$ and write $j\approx rm$ as an abbreviation for $|j-rm|\leq(rm)^{1/3+\delta}$. We have

[TABLE]

For larger $j$, we have

[TABLE]

The tail sum is very small. For example, comparing to an integral gives

[TABLE]

with $\varepsilon=3\delta/2$, $x\approx c^{2/(3\delta)}k$ so that $x^{\varepsilon}\approx ck^{3\delta/2}$, and $L=c^{1/\varepsilon}rm$. Meanwhile, the “bulk” contribution is

[TABLE]

For the main term, we have a Riemann sum approximation

[TABLE]

For the error term, we have

[TABLE]

Thus the bulk contribution is

[TABLE]

Compare this to the edge contribution, bounded by $(rm)^{(-4/3+a)p+1/3}$, and the tail contribution, which obeys the even stronger bound $\exp(-c(rm)^{3\delta/2})$. For any chosen $a<1/3$, these are negligible compared to the main term of order $(rm)^{-p+1}$ and even to the error term of order $(rm)^{-p+1-1/3+a}$. Thus we have, for the sum over all three ranges,

[TABLE]

with an implicit constant proportional to $p\int_{0}^{1}(1-u^{2})^{p/2}u^{2/3+a}du$. In particular,

[TABLE]

From this, we obtain

[TABLE]

which is negligible as $rm\rightarrow\infty$ provided $p\geq 3$. Even when we sum over all $p\geq 3$, we obtain a geometric series:

[TABLE]

A final detail remains: We took $s=it/\sigma$, which we now justify by using the semicircle law to show that the power series for $\log(1-2s\lambda_{j})$ does converge at that point. From the semicircle law Proposition 10.1, the largest $\lambda_{j}$ are of order $1/(rm)$, whereas $\sigma$ is of order $1/\sqrt{rm}$. It follows that, for any given $t\in\mathbb{R}$, once $rm$ is sufficiently large we do have

[TABLE]

for all $j$. The series will converge once $rm$ is so large that

[TABLE]

Combining these estimates, we have

[TABLE]

and hence

[TABLE]

For any fixed $t$, this converges to $e^{-t^{2}/2}$ as $rm\rightarrow\infty$, which completes our proof of Theorem 1.1.

10.3. Failure of the CLT at the wave scale

Suppose that $rT=c$ for a fixed constant $c$, or more generally that $rT$ remains bounded as $T\rightarrow\infty$. At this scale, we do not expect the central limit theorem to hold, and a direct calculation on $S^{2}$ supports this. From the variance formula (3.1) with $p=2$,

[TABLE]

We have $\lambda_{j}>0$ because they are the eigenvalues of the positive quadratic form $\int_{B}\phi^{2}$. It follows that

[TABLE]

where $\lambda_{0}$ is the eigenvalue corresponding to the zonal spherical harmonic. As in the proof of the semicircle law Proposition 10.1, we have

[TABLE]

But unlike Proposition 10.1, we now have $rm=c$ of order 1, so that

[TABLE]

It follows that

[TABLE]

Thus the criterion of Proposition 1.3 fails, and we believe that instead of converging to a Gaussian, the standardized local integral $Z$ converges to a weighted sum of finitely many squares of Gaussians (with a number of degrees of freedom proportional to $c$, the higher values of $k$ being exponentially suppressed).

11. Tail probability

One consequence of Theorem 1.1 is that the random variable $X=\int_{B}\phi^{2}$ has tails satisfying

[TABLE]

for any fixed $y>0$. Note that $\operatorname{var}[X]\asymp(rT)^{-n+1}\rightarrow 0$, and another approach is needed to control deviations where $|X-\mathbb{E}[X]|$ exceeds a given $\varepsilon>0$ of larger size. One obtains such a bound using the foregoing estimates together with a Chernoff bound. Namely, for any $s\geq 0$ for which $\mathbb{E}[e^{sX}]$ is defined, the upper tail probability obeys

[TABLE]

The lower tail event $X-\mathbb{E}[X]<-\varepsilon$ can be treated with the same procedure applied to $-X$ instead of $X$. The parameter $s$ is at our disposal, and the optimal choice is to minimize $\log\mathbb{E}\left[e^{sX}\right]-s\mathbb{E}[X]-s\varepsilon$. As we saw in Section 2,

[TABLE]

with both sides defined for $1-2s\lambda_{\max}>0$, where $\lambda_{\max}$ is the largest $\lambda_{j}$. Differentiating with respect to $s$, we find that the equation for the optimal $s$ is

[TABLE]

Upon expanding $\lambda_{j}/(1-2s\lambda_{j})$ in a geometric series, the first term $p=1$ cancels $\sum_{j}\lambda_{j}$ and we obtain

[TABLE]

Truncating at $p=2$ suggests the choice

[TABLE]

but it is not clear whether $1-2s\lambda_{\max}>0$, that is, whether this choice of $s$ is valid. Using our bound on the third moment $\operatorname{tr}(A^{3})\lesssim(rT)^{-3n/2+1}$ we conclude that

[TABLE]

With the desired choice of $s$,

[TABLE]

by the third-moment bound for $\lambda_{\max}$ and the fact that $\sum\lambda_{j}^{2}$ is bounded above and below by constant multiples of $(rT)^{-n+1}$. It follows that the choice (11.2) is valid provided

[TABLE]

or more generally when $\varepsilon$ is a sufficiently small constant multiple of $(rT)^{-n/2+2/3}$.

When $\varepsilon$ is small enough that (11.2) is a valid choice of $s$, the Chernoff bound gives

[TABLE]

assuming we can truncate the Taylor expansion at second order. This truncation is justified under the same condition allowing the choice of $s$ to begin with, since the higher terms can be bounded by $\varepsilon^{3}\left(\sum_{j}\lambda_{j}^{2}\right)^{-3}\sum_{j}\lambda_{j}^{3}.$ This completes the proof of Theorem 1.7. ∎

Even when (11.2) is not valid, one can obtain an upper bound by choosing a different $s$ (and (11.2) is itself only an approximation to the optimal choice). The workaround from [11] is to instead take $s$ to be a small multiple of \big{(}\sum\lambda_{j}^{2}\big{)}^{-1/2}, the validity of which follows from the bound $\lambda_{\max}\leq\sqrt{\sum\lambda_{j}^{2}}$. This allows one to take, for example, $\varepsilon$ of order 1 at the cost of replacing $(rT)^{n-1}$ by $(rT)^{(n-1)/2}$ in Theorem 1.7.

12. Conclusion and outlook

As explained in Section 5, the exponent $(rT)^{-n+1}$ in our lower bound on the variance is optimal. Less clear is the true size of the third moment. One approach to this would be to give an explicit treatment of higher-dimensional spheres along the lines of the semicircle law for $S^{2}$.

One of the aspects of the Central Limit Theorem is its universality: The normal distribution appears as a limit for a wide class of underlying distributions for the coefficients. In contrast, we have assumed Gaussian coefficients from the beginning, so it is perhaps less surprising that one recovers a Gaussian in the limit. An interesting problem would be to prove the CLT for random waves with other distributions for the coefficients $c_{j}$, not necessarily Gaussian. We expect that the central limit theorem will hold for a quite broad class of random coefficients. Indeed, it might not be necessary to randomize at all if the underlying manifold $M$ has favourable dynamical properties.

Some further examples of coefficients seem to be within reach, namely if (a) the vector of coefficients is rotation-invariant, or (b) several moments of the coefficient distribution agree with those of a Gaussian. Assuming (a), the vector of coefficients is rotation-invariant and as above one can diagonalize the quadratic form $X=\sum_{j}\sum_{k}c_{j}c_{k}\int_{B}\phi_{j}\phi_{k}$ to obtain a sum of independent random variables. Lyapunov’s criterion then applies and reduces the central limit theorem to estimates for the second and third moments as above. Assuming (b) up to moments of order six, these estimates follow from the direct calculation with Gaussians carried out here. Note that the $p$-th moment of the quadratic form $X=\int_{B}\phi^{2}$ only involves moments up to order $2p$ of the coefficients of $\phi$. For an example of rotation-invariant but non-Gaussian coefficients, one can take the vector of coefficients uniformly from the sphere $S^{N}$. However, the coefficients in this example are not independent.

We have given an upper bound on the probability that $\int_{B}\phi^{2}$ deviates from its average, at least in some ranges. An interesting problem would be to use the theory of large deviations to give a comparable lower bound on this probability. It could also be possible to prove the Central Limit Theorem directly from this machinery.

Acknowledgments

We thank Yaiza Canzani for helpful discussions about Weyl’s law. MDCI thanks Peter Sarnak for his advice, encouragement, and support over the course of this work and thanks the Natural Sciences and Engineering Research Council of Canada for a PGS D grant. ML thanks Rupert L. Frank for multiple stimulating discussions on the topic, and thanks the Institute for Advanced Study for its hospitality during the 2017-2018 academic year.