text:
\begin{document}

\maketitle

\section{Introduction}

A  Weyl law for Toeplitz operators was proved
by Boutet de Monvel and Guillemin  (\cite{bg}, Theorem
13.1), using the theory of Fourier-Hermite distributions and their symbolic calculus based on
symplectic spinors. The result holds in the setting of so-called Toeplitz structures, a generalization
of strictly pseudoconvex domains, and the operators involved act on certain spaces
of half-forms. The circle bundle of a positive line bundle on a symplectic manifold, with an appropriate
connection, is an important
special case of a Toeplitz structure.

On the other hand, a circle of ideas has emerged recently
(\cite{z}, \cite{bsz}, \cite{sz}), in which asymptotic properties of positive line bundles
are studied by appealing directly to the microlocal description of Szeg\"{o} kernels
in \cite{bs}, and to its generalization to the almost complex symplectic setting in \cite{sz}.
These ideas and techniques lead to
several rather explicit local and global
asymptotic expansions related to the (generalized) Szeg\"{o} kernel of the line
bundle. It seems natural to ask whether the Weyl law can also be interpreted in light of this approach.

Accordingly, the aim of this paper is to provide an alternative,
perhaps relatively elementary proof of the
theorem of Boutet de Monvel and Guillemin, in the case of Toeplitz structures related to positive line bundles
on symplectic manifolds. As is well known, this setting is natural in algebraic and symplectic geometry,
and in the theory of geometric quantization.
For ease of exposition, we shall state our results in the context of
complex projective geometry, but in view of the theory of
\cite{sz} the present arguments apply to the more general
almost complex symplectic category.

In short, by pairing techniques from \cite{sz} with certain classical arguments used to estimate the spectral function
of a pseudodifferential operator \cite{h}, \cite{gs}, we shall obtain local Weyl
laws for the spectral function of a Toeplitz operator; the global Weyl law then follows by integration.
Unlike \cite{bg}, the operators will be understood to act on half-densities, rather that half-forms.

Before stating the result in point,
it is in order to fix some conventions, and to recall some basic notions from
\cite{bg}.

Suppose $M$ is a d-dimensional complex projective manifold, and let $A$ be an ample line
bundle on $M$. There exists an Hermitian metric on $A$, such that the
unique compatible connection has normalized curvature $\Theta=-2i\,\omega$, where
$\omega$ is a K\"{a}hler form. The induced volume form $dV_M=\frac{1}{\mathrm{d}!}\,\omega^{\wedge \mathrm{d}}$
on $M$ has
total volume $\mathrm{vol}(M)=\frac{\pi^\mathrm{d}}{\mathrm{d}!}\,\int _Mc_1(A)^\mathrm{d}$.

With this assumption, the unit disc bundle $D\subseteq A^*$ is a strictly
pseudoconvex domain, with boundary $X=\partial D$ the unit circle bundle in $A^*$.
The compatible connection on $A$ determines
the connection form $\alpha\in \Omega^1(X)$, and by positivity
$(X,\alpha)$ is a contact manifold. In particular, $X$ has the volume form
$\frac{1}{2\pi}\,dV_X=\frac{1}{2\pi}\,\alpha\wedge \pi^*(dV_M)$;
given this, we shall identify densities, half-densities, and functions on $X$.

Given the $L^2$-product on $X$ induced by
$d\mu_X=:\frac{1}{2\pi}\,dV_X$,
let $H(X)\subseteq L^2(X)$ denote  the
Hardy space of $X$, and $\Pi:L^2(X)\rightarrow H(X)$ the orthogonal projector. The Szeg\"{o} kernel of $X$ is the
distributional kernel of $\Pi$, which we shall denote by the same symbol, $\Pi\in \mathcal{D}'(X\times X)$.
If $\{s_j\}$ is an orthonormal basis of $H(X)$, then $\Pi(x,y)=\sum_{j=1}^{+\infty}s_j(x)\,\overline{s_j(y)}$ and
\begin{equation}
\label{eqn:pairing-density}
\Pi(f)(x)=\int _X\Pi(x,y)\,f(y)d\mu_X(y)\,\,\,\,\,\,\,\,\,\,\left(f\in L^2(X),\,x\in X\right).
\end{equation}

We shall deal with the following class of operators:

\begin{defn}
\label{defn:toeplitz}
If $m\in \mathbb{Z}$, a Toeplitz operator of order $m$ on $X$ is an operator
$T:\mathcal{D}'(X)\rightarrow \mathcal{D}'(X)$
of the form
$T=:\Pi\circ P\circ \Pi$, where $P$ is a pseudodifferential operator of classical type
and order $m$ on $X$.
\end{defn}

The composition in Definition \ref{defn:toeplitz}
is well-defined,
in view of the microlocal structure of $\Pi$ \cite{bs} and standard results on wave fronts \cite{d}.

Now $\alpha$ being a connection form on $X$ is equivalent to the condition that
the closed cone
\begin{equation}
\label{eqn:sigma-alpha}
\Sigma=:\left\{(x,r\alpha_x):x\in X,\,r>0\right\}\subseteq T^*X\setminus \{0\}
\end{equation}
be a symplectic submanifold of the cotangent bundle of $X$, equipped with its canonical symplectic structure.
By the theory in \S 1 and \S 2 of \cite{bg},
the following definition is well-posed:

\begin{defn}
\label{defn:toeplitz-symb}
Let $T:\mathcal{D}'(X)\rightarrow \mathcal{D}'(X)$ be a Toeplitz operator of order $m$,
and let $\Sigma$ be as in (\ref{eqn:sigma-alpha}).
If $P$ is as in Definition
\ref{defn:toeplitz}, the \textit{symbol} $\sigma_T$ of $T$ is the $\mathcal{C}^\infty$ function
on $\Sigma$ given by the restriction of the symbol of $P$.
Furthermore, $T$ is \textit{elliptic} if $\sigma_T$ is the restriction to $\Sigma$
of an elliptic symbol of degree $m$.
\end{defn}

Thus a Toeplitz operator $T$ determines a $\mathcal{C}^\infty$-function
$\varsigma_T$ on $X$, defined by
$\varsigma_T(x)=:\sigma_T(x,\alpha_x)$ ($x\in X$).

We shall say that the Toeplitz operator $T$ is self-adjoint
to mean that it is formally self-adjoint with respect to
the $L^2$-product associated to $d\mu_X$; in this case, $T$ is
an essentially self-adjoint operator with domain $\mathcal{C}^{\infty}(X)$,
and its symbol $\sigma_T$ is real valued.

Suppose that $T$ is a self-adjoint and elliptic Toeplitz operator
of degree $1$, with
$\sigma_T>0$. By the theory in \S 1 and \S 2 of \cite{bg},
$T\ge -C\,I$ for some constant $C\in \mathbb{R}$.
Furthermore, its spectrum $\mathrm{Spec}(T)$ is discrete, bounded from below and has
only $+\infty$ as an accumulation point (Proposition 2.14 of \cite{bg}); if we view $T$ as
an operator $T_H$ with domain $\mathcal{C}^{\infty}(X)\cap H(X)$,
every eigenvalue occurs with finite multiplicity.

Thus $\mathrm{Spec}\left(T_H \right)=\big\{\lambda_j\big\}\subseteq \mathbb{R}$, with $\lambda_1\le \lambda_2\le \cdots$,
where each $\lambda_j$ is repeated according to its multiplicity, and $\lambda_j\rightarrow +\infty$ for $j\rightarrow
+\infty$. For every $j=1,2,\ldots$, let $e_j\in H (X)\cap \mathcal{C}^\infty(X)$ be an eigenfunction of
$T$ corresponding to the eigenvalue $\lambda_j$ and of unit $L^2$-norm.
We can arrange that $\{e_j\}$ is an orthonormal basis of $H(X)$.

For every $\lambda\in \mathbb{R}$,
let us set
$$
\mathbb{T}_\lambda=:\mathrm{span}\Big\{e_j:\lambda_j\le \lambda\Big\},
$$
a finite dimensional subspace of $H(X)$, and denote by
$
\mathfrak{T}_\lambda :L^2(X)\rightarrow \mathbb{T}_\lambda
$ the orthogonal projector:
$\mathfrak{T}_\lambda$ is the smoothing operator with distributional kernel
the \textit{spectral function}
$$
\mathcal{T}(\lambda,x,y)=:\sum _{\lambda_j\le \lambda}e_j(x)\cdot \overline{e_j(y)}\,\,\,\,\,\,\,\,\,\,\,\,
(x,y\in X).
$$
We shall also let $\mathcal{T}(\lambda,x)=:\mathcal{T}(\lambda,x,x)$ ($x\in X,\,\lambda\in \mathbb{R}$).

\begin{thm}
\label{thm:local-weyl-law}
Let $T$ be a self-adjoint first order elliptic Toeplitz operator on $X$ with $\sigma_T>0$.
Then as $\lambda\rightarrow +\infty$ we have
$$
\mathcal{T}(\lambda,x)=\frac{\pi}{d+1}\cdot\left(\frac{\lambda}{\pi\, \varsigma_T(x)}\right)
^{\mathrm{d}+1}+O\left(\lambda^\mathrm{d}\right),
$$
uniformly in $x\in X$.
\end{thm}

The \textit{counting function} of $T$ is
\begin{equation}
\label{eqn:counting-fctn}
N_T(\lambda)=:\sharp \Big\{j\in \mathbb{N}\,:\,\lambda_j\le \lambda\Big\}\,\,\,\,\,\,\,\,\,\,\,\,
(\lambda\in \mathbb{R}).
\end{equation}
Clearly,
\begin{equation}
\label{eqn:counting-trace}
N_T(\lambda)=\mathrm{trace}(\mathfrak{T}_\lambda)=\int _X \mathcal{T}(\lambda,x)\,\mathrm{vol}_X(x).
\end{equation}

To check our normalizations, suppose
$T=\Pi\circ \left(\frac{1}{i}\,\frac{\partial}{\partial \theta}\right)\circ \Pi$, where $\frac{\partial}{\partial \theta}$
denotes the generator of the $S^1$-action.
Then by the Riemann-Roch Theorem
\begin{equation}
\label{eqn:riemann-roch}
N_T(\lambda)=\sum _{k\le \lambda}h^0\left (M,A^{\otimes k}\right)=\frac{1}{(\mathrm{d}+1)!}\,
\lambda^{\mathrm{d}+1}\,\int _M c_1(A)^\mathrm{d}+O\left(\lambda^\mathrm{d}\right)
\end{equation}
for $\lambda\rightarrow +\infty$.
Now Theorem \ref{thm:local-weyl-law} says that
$\mathcal{T}(\lambda,x)=\frac{\pi}{d+1}\cdot\left(\frac{\lambda}{\pi}\right)
^{\mathrm{d}+1}+O\left(\lambda^\mathrm{d}\right)$, and therefore predicts the estimate
$$N_T(\lambda)=\frac{\pi}{d+1}\cdot\left(\frac{\lambda}{\pi}\right)
^{\mathrm{d}+1}\cdot \mathrm{vol}(M)+O\left(\lambda^\mathrm{d}\right),
$$
in agreement with (\ref{eqn:riemann-roch}).

More generally, the Weyl law of \cite{bg} follows by using Theorem
\ref{thm:local-weyl-law} in (\ref{eqn:counting-trace}):
\begin{cor}
\label{cor:global-weyl-law}
Let
$\Sigma_1\subseteq \Sigma$ be the locus where $\sigma_T<1$.
In the situation of Theorem \ref{thm:local-weyl-law},
as $\lambda\rightarrow +\infty$ we have
$$
N_T(\lambda)=
\left(\frac{\lambda}{2\pi}\right)^{\mathrm{d}+1}\,\mathrm{vol}(\Sigma_1)+O\left(\lambda^\mathrm{d}\right).
$$
\end{cor}

\section{Proof of Theorem \ref{thm:local-weyl-law}.}

By Lemma 12.1 of \cite{bg}, there exists a
self-adjoint first-order pseudodifferential operator of classical type
$Q$ on $X$ such that
$[Q,\Pi]=0$, $q>0$ everywhere, and $T=\Pi\circ Q\circ \Pi$; here
$q:T^*X\setminus \{0\}\rightarrow \mathbb{R}$ is the symbol of
$Q$.

Let us consider the one-parameter group of $L^2$-unitary operators
$$
U(\tau)=:e^{i\tau Q}\,\,\,\,\,\,\,(\tau \in \mathbb{R})
$$
(\S 12 of \cite{gs}, \S 12 of \cite{bg}).

The $\mathcal{C}^\infty$-function $q$
generates a one-parameter group $\Phi_\tau$ ($\tau\in \mathbb{R}$)
of homogeneous symplectomorphisms
of $T^*X\setminus \{0\}$. For every $\tau\in \mathbb{R}$
the graph of $\Phi_\tau$, $\mathrm{graph}(\Phi_\tau)$, is therefore
a conic Lagrangian relation, and
$U(\tau)$ is a Fourier integral operator of degree $0$ associated to it.
In fact, there exists $\epsilon_0>0$ such that if $|\tau|<\epsilon_0$
then $U(\tau)=V(\tau)+R(\tau)$, where $V(\tau)$ and $R(\tau)$ are as follows.

First, $R(\tau)$ is a smoothing operator, and in fact
its kernel $R(\tau,x,y)$ is a $\mathcal{C}^\infty$ function on
$\big((-\epsilon_0,\epsilon_0)\times X\times X\big)$ (Proposition 12.3 of \cite{gs}).

Next,
$V(\tau)$ can be locally
represented by a Fourier integral in the following manner.
On a coordinate patch, with abuse of notation let us identify $x,y,\cdots\in X$
with their local coordinates in $\mathbb{R}^{2\mathrm{d}+1}$.
Let
$\varphi (\tau,\cdot,\cdot)$ be a generating function
for $\Phi_\tau$. Since $\Phi_0$ is the identity and $\varphi$ satisfies the Hamilton-Jacobi equation,
\begin{equation}
\label{eqn:generating-for-small-tau}
\varphi(\tau,x,\eta)=x\cdot \eta +\tau \,q(x,\eta)+O\left(\tau^2\right)\cdot \|\eta\|.
\end{equation}
Then
\begin{equation}
\label{eqn:local-representation-V}
V(\tau)(x,y)=\frac{1}{(2\pi)^{2\mathrm{d}+1}}\,
\int _{\mathbb{R}^{2\mathrm{d}+1}}e^{i[\varphi(\tau,x,\eta)-y\cdot \eta]}\,a(\tau,x,y,\eta)\,d\eta.
\end{equation}
Here $a(\tau,\cdot,\cdot,\cdot)\in S^0_{\mathrm{cl}}$
for every $\tau$. To determine $a(0,\cdot,\cdot,\cdot)$, remark that
$U(0)=\mathrm{Id}$; therefore, if $V(\tau)$ was to act on functions by means of the Euclidean density $dx$
on $\mathbb{R}^{2\mathrm{d}+1}$, we would have $a(0,\cdot,\cdot,\cdot)=1$. Since however $V(\tau)$ acts
on functions by means of the density $d\mu_X$ (that is, as in (\ref{eqn:pairing-density}) with $V(\tau)$
in place of $\Pi$), we must have $a(0,\cdot,\cdot,\cdot)=1/\mathcal{V}(y)$, where
$\mathcal{V}\,dx$ is the local coordinate expression of $d\mu_X$.

By \cite{bs}, $\Pi$ is also a Fourier integral operator, with complex phase however. More precisely,
 $\Pi=\Pi _1+\Pi_2$, where $\Pi_2$ is smoothing, while locally
\begin{equation}
\label{eqn:local-representation-Pi}
\Pi_1(x,y)=\int _0^{+\infty}e^{it\psi(x,y)}\,s(t,x,y)\,dt,
\end{equation}
where $\Im \psi\ge 0$ and $s(t,x,y)\sim \sum _{j=0}^{+\infty}t^{\mathrm{d}-j}\,s_j(x,y)$.

If $\epsilon>0$ and
$\chi\in \mathcal{C}^\infty _0\big((-\epsilon,\epsilon)\big)$, for any suitable family of operators
$A(\tau)$ we set (see \S 12 of \cite{gs})
\begin{equation*}
\label{eqn:A-chi}
A_\chi=:\int _{-\epsilon}^\epsilon \chi(\tau)\,A(\tau)\,d\tau.
\end{equation*}

Let us pick a sufficiently small $\epsilon \in (0,\epsilon_0)$, and apply this to
$S(\tau)=:\Pi\circ U(\tau)\circ \Pi$ ($\tau\in \mathbb{R}$).
By the hypothesis, actually $S(\tau)=U(\tau)\circ \Pi$.
Thus,
\begin{equation}
\label{eqn:s-chi}
S_\chi=U_\chi\circ \Pi=V_\chi\circ \Pi+R_\chi\circ \Pi.
\end{equation}
Since $U_\chi$ is a smoothing operator \cite{gs}, so is $S_\chi$.

When $\chi$ is replaced by $\chi \cdot e^{-i\lambda \,(\cdot)}$, and (\ref{eqn:s-chi})
is rewritten with some abuse of language in terms of distributional kernels, we obtain
\begin{eqnarray}
\label{eqn:s-chi-kernel}
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,y)&=&\left( U_{\chi \,e^{-i\lambda \,(\cdot)}}\circ \Pi\right)(x,y)\\
&=&\left( V_{\chi \,e^{-i\lambda \,(\cdot)}}\circ \Pi\right)(x,y)
+\left( R_{\chi \,e^{-i\lambda \,(\cdot)}}\circ \Pi\right)(x,y)
\,\,\,\,\,\,\,(x,y\in X).\nonumber
\end{eqnarray}
We shall now study the asymptotics of (\ref{eqn:s-chi-kernel}) for $\lambda\rightarrow \infty$.

To begin with, the kernel $\mathcal{R}(t,x,y)$ of $R(t)\circ \Pi$ is a $\mathcal{C}^\infty$ function on
$(-\epsilon_0,\epsilon_0)\times X\times X$.
The second summand in (\ref{eqn:s-chi}) is then given by
\begin{equation}
\label{eqn:rapid-decay}
\int _{-\epsilon}^\epsilon e^{-i\lambda t}\,\chi(t)\,\mathcal{R}(t,x,y)\,dt=\widehat{\chi\,\mathcal{R}}(\lambda,x,y)=
O\left (\lambda ^{-\infty}\right)
\end{equation}
for $\lambda\rightarrow \infty$, uniformly in $x,y\in X$;
here $\,\widehat{}\,$ is the
Fourier transform with respect to $t$.

Let $\sim$ denote equality of asymptotic expansions.
By (\ref{eqn:s-chi}) and (\ref{eqn:rapid-decay}), for all $x\in X$
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag}
\lefteqn{
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim \left(V_{\chi \,e^{-i\lambda \,(\cdot)}}\circ \Pi\right)(x,x)}\\&=&\int _X\int_{-\epsilon}^\epsilon e^{-i\lambda \tau}\,\chi (\tau)\,
V(\tau)(x,z)\,\Pi(z,x)\,d\mu_X(z)\,d\tau.                                                                           \nonumber
\end{eqnarray}

Let us fix $x\in X$ and an arbitrarily small open neighborhood
$X_1\subseteq X$; let furthermore
$X_2\subseteq X$ be an open neighborhood of $X\setminus X_1$, such that $x\not\in \overline{X_2}$.
Let $\big\{\varrho_j\big\}$ be a a partition of unity
on $X$ subordinate to the open cover $\mathcal{U}=:\left\{X_1,X_2\right\}$ of $X$.
We obtain from (\ref{eqn:s-chi-lambda-diag}):
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part}
\lefteqn{
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim \left( V_{\chi \,e^{-i\lambda \,(\cdot)}}\circ \Pi\right)(x,x)}\\
&=&\sum _{j=1}^2\int _{X_j}\int_{-\epsilon}^\epsilon e^{-i\lambda \tau}\,\chi (\tau)\,\varrho _j(z)\,
V(\tau)(x,z)\Pi(z,x)\,d\mu_X(z)\,d\tau.                                                                           \nonumber
\end{eqnarray}

\begin{lem}
\label{lem:rapid-decay-second-term}
As $\lambda\rightarrow \infty$, we have:
$$
\int _{X_2}\int_{-\epsilon}^\epsilon e^{-i\lambda \tau}\,\chi (\tau)\,\varrho _2(z)
V(\tau)(x,z)\,\Pi(z,x)\,d\mu_X(z)\,d\tau=O\left(\lambda^{-\infty}\right).
$$
\end{lem}

\textit{Proof.} Since the singular support of
$\Pi$ is the diagonal in $X\times X$,
$\Pi(\cdot,x)$ is
$\mathcal{C}^\infty$ on $X_2$.
If $p_x(z)=:\varrho _2(z)\,\Pi(z,x)$, we have $p_x\in \mathcal{C}^\infty_0(X_2)$.

Now $V(\tau)$ is a $\mathcal{C}^\infty$ family of distributions on $X\times X$,
and by regularity
$V(\tau,x)=:V(\tau)(x,\cdot)$ is therefore a  well-defined $\mathcal{C}^\infty$-family of distributions on $X$,
parametrized by $\mathbb{R}\times X$.

Consequently,
$\gamma (\tau)=:\chi(\tau)\,\langle V(\tau,x),p_x\rangle$ is a compactly supported
$\mathcal{C}^\infty$ function of $\tau$, hence
its Fourier transform is rapidly decreasing.

\hfill Q.E.D.

\medskip

\begin{rem}
If $X_2$ is defined by an inequality on the distance
function from $x$ (using the natural metric on $X$), and the same
bound is used for every $x\in X$, the previous statement
holds uniformly in $x$.
\end{rem}

Recalling that $\Pi=\Pi_1+\Pi_2$, we conclude that

\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim }\\&=&\sum _{j=1}^2\int _{X_1}
\int_{-\epsilon}^\epsilon e^{-i\lambda \tau}\,\chi (\tau)\,\varrho _1(z)
V(\tau)(x,z)\,\Pi_j(z,x)\,d\mu_X(z)\,d\tau.\nonumber
\end{eqnarray}

Since $\Pi_2$ is smoothing, a similar argument to the proof of Lemma \ref{lem:rapid-decay-second-term}
shows that the second summand in (\ref{eqn:s-chi-lambda-diag-part-reduced}) is $O\left(\lambda^{-\infty}\right)$
as $\lambda\rightarrow \infty$.
To estimate asymptotically the first summand, we make use of (\ref{eqn:generating-for-small-tau}),
(\ref{eqn:local-representation-V}) and (\ref{eqn:local-representation-Pi}) and obtain

\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim }\\
&=&\frac{1}{(2\pi)^{2\mathrm{d}+1}}\,\int _{X_1}\int _0^{+\infty}\int_{-\epsilon}^\epsilon \int_{\mathbb{R}^{2\mathrm{d}+1}}
e^{i\,\Phi_1(x,z,t,\tau,\eta,\lambda)}\,A(x,z,t,\tau,\eta)\,\,d\mu_X(z)\,dt\,d\tau\,d\eta, \nonumber
\end{eqnarray}
where
\begin{eqnarray}
\label{eqn:phase-A}
\Phi_1(x,z,t,\tau,\eta,\lambda)&=:&\varphi (\tau,x,\eta)-z\cdot \eta+t\psi(z,x)-\lambda \tau\\
&=&
(x-z)\cdot \eta+\tau\,q(x,\eta)+t\psi(z,x) -\lambda \tau+O\left(\tau^2\right)\cdot \|\eta\|,\nonumber
\end{eqnarray}
and
\begin{eqnarray}
\label{eqn:amplitude-A}
A(x,z,t,\tau,\eta)=:\chi (\tau)\,\varrho_1(z)\,a(\tau,x,z,\eta)\,s(t,z,x).
\end{eqnarray}

Since $q$ is positively homogenous of degree one in $\eta$ and everywhere positive for
$\eta\neq \mathbf{0}$, (\ref{eqn:phase-A}) implies that
$\partial _\tau\Phi\sim \|\eta\|+|\lambda|$ for $\lambda\rightarrow -\infty$; integrating by parts in $\tau$,
$S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)=O\left(\lambda^{-\infty}\right)$ for $\lambda\rightarrow
-\infty$.

To study the asymptotics of (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier})
for $\lambda\rightarrow +\infty$, following \cite{gs} we choose a radial function
$F\in \mathcal{C}^\infty _0\left(\mathbb{R}^{2\mathrm{d}+1}\setminus \{\mathbf{0}\}\right)$
identically equal to $1$ for $\frac 1C\le \|\eta \|\le C$ for some $C\gg 0$,
and split (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier}) as:
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-spezzato}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim }\\
&=&\frac{1}{(2\pi)^{2\mathrm{d}+1}}\,\left\{\int _{X_1}\int _0^{+\infty}\int_{-\epsilon}^\epsilon \int_{\mathbb{R}^{2\mathrm{d}+1}}
e^{i\,\Phi_1(x,z,t,\tau,\eta,\lambda)}\,A(x,z,t,\tau,\eta)\,F\left (\frac \eta \lambda \right)\,d\mu_X(z)\,dt\,d\tau\,d\eta \right.\nonumber \\
&&+\left.\int _{X_1}\int _0^{+\infty}\int_{-\epsilon}^\epsilon \int_{\mathbb{R}^{2\mathrm{d}+1}}
e^{i\,\Phi_1(x,z,t,\tau,\eta,\lambda)}\,A(x,z,t,\tau,\eta)\,\left[1-F\left (\frac \eta \lambda \right)\right]\,d\mu_X(z)\,dt\,d\tau\,d\eta\right\}.\nonumber
\end{eqnarray}
Since $\|\partial _\tau \Phi\|\sim
\|\eta\|+\lambda$ where $1-F\left (\frac \eta \lambda \right)\neq 0$, the second term in
(\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-spezzato}) is
$O\left(\lambda ^{-\infty}\right)$ for $\lambda\rightarrow +\infty$.

To estimate the first term in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-spezzato}),
let us
perform the change of variables $t\rightsquigarrow \lambda t$,
$\eta \rightsquigarrow \lambda \eta$, and use the homogeneity of $\Phi$ to obtain:
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-riscalato}
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)
&\sim&\frac{\lambda^{2\mathrm{d}+2}}{(2\pi)^{2\mathrm{d}+1}}\,\int _{X_1}\int _0^{+\infty}\int_{-\epsilon}^\epsilon \int _{1/D}^D
\int_{S^{2\mathrm{d}}}
e^{i\,\lambda\,\Phi_1(x,z,t,\tau,r\omega,1)}\,A(x,z,\lambda t,\tau,\lambda r\omega)\nonumber\\
&&\cdot r^{2\mathrm{d}}\,F(r)\,
d\mu_X(z)\,dt\,d\tau\,dr\,d\omega,
\end{eqnarray}
where now integration in $dr$ takes place over the interval
$\frac 1D\le r\le D$ for some $D\gg 0$.

To proceed further, we shall need a more specific choice of
a system of local coordinates near $x$. It
will be useful to work in a system of \textit{Heisenberg local coordinates}
centered at $x$, as defined in \cite{sz}, for in these
systems of coordinates scaling limits of Szeg\"{o} kernel
exhibit their universal nature (the focus in \cite{sz} is
on the asymptotics for $k\rightarrow +\infty$, where $k$ indexes the
isotype with respect to the $S^1$-action, but we shall exploit some
formal analogies with the present computations). We refer to
\cite{sz} for a definition of Heisenberg local
coordinates.

Suppose then that $(\theta,\mathrm{v}):X_1\rightarrow (-a,a)\times B_{2\mathrm{d}}\big(\mathbf{0},\delta\big)$
is a system of Heisenber local coordinates centered at $x$, for some $a,\delta>0$; here $B_{2\mathrm{d}}\big(\mathbf{0},\delta\big)\subseteq \mathbb{R}^{2\mathrm{d}}
\cong \mathbb{C}^\mathbf{d}$ is the open ball of radius $\delta$ centered at the origin.
Following \cite{sz},
let $z=x+(\theta,\mathbf{v})$ denote the point with coordinates $(\theta,\mathbf{v})$.
Furthermore, we shall write
$\omega =(\omega_0,\omega_1)$, where $\omega_0\in \mathbb{R}$, $\omega_1\in \mathbb{R}^{2\mathrm{d}}$ and
$\omega_0^2+\|\omega_1\|^2=1$.
On the upshot, the phase in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-riscalato})
becomes
\begin{eqnarray}
\label{eqn:phase-phi-heis}
\lefteqn{\Phi_2\big(x,\theta,\mathbf{v},t,\tau,r,\omega\big)=:\Phi_1(x,z,t,\tau,r\omega,1)}\\
&=:&-r\,\mathbf{v}\cdot\omega_1
+t\psi\Big(x+\left(\theta,\mathbf{v}\right),x\Big)-r\,\theta \,\omega_0
+r\tau\,q\big(x,\omega\big)-\tau+O\left(\tau^2\right)\,r.
\nonumber
\end{eqnarray}
The paring $\mathbf{v}\cdot\omega_1$ is the standard scalar product as vectors in $\mathbb{R}^{2\mathrm{d}}$.

Let $d\mu_X(z)=\mathcal{V}(\theta,\mathbf{v})\,d\theta\,d\mathbf{v}$
be the volume density expressed in local coordinates.
Integration in $d\mathbf{v}$ takes place over the open ball
$B_{2\mathrm{d}}(\mathbf{0},\delta)\subseteq \mathbb{C}^\mathrm{d}$.
By definition of Heisenberg local coordinates,
$\mathcal{V}(\theta,\mathbf{0})=1/2\pi$ for every $\theta$.

We next show that, up to a rapidly decreasing contribution, the integration in $d\omega$ can be
restricted to an arbitrarily small open neighborhood of $(1,\mathbf{0})$ in $S^{2\mathrm{d}}$.
Suppose $0<\delta_1<\delta_2<1$, and define
$$S_-=:\left\{\omega \in S^{2\mathrm{d}}:\omega _0<\delta_2\right\},\,\,\,\,\,
S_+=:\left\{\omega \in S^{2\mathrm{d}}:\omega _0>\delta_1\right\}.
$$
Let $\{\beta_-,\beta_+\}$ be a partition of unity subordinate to
the open cover $\mathcal{U}=\left\{S_-,S_+\right\}$
of $S^{2\mathrm{d}}$.
Inserting the identity $\beta_-+\beta_+=1$ in
(\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-riscalato}), we obtain
$$
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)_-+
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)_+,
$$
where in the former term the integrand in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-riscalato})
has been multiplied by $\beta_-(\omega)$, in the latter by $\beta_+(\omega)$.

\begin{lem}
\label{lem:reduction-on-sphere}
If $\epsilon$ is sufficiently small,
then $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)_-=O\left(\lambda ^{-\infty}\right)$ when
$\lambda\rightarrow +\infty$.
\end{lem}

\textit{Proof.}
If in the integral
representation (\ref{eqn:local-representation-V}) for $V(\tau)$ a cut-off is introduced in the neighborhood of
$\eta=0$, only a smoothing term is lost. The argument leading to (\ref{eqn:rapid-decay})
implies that the corresponding contribution to $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)$
is $O\left(\lambda^{-\infty}\right)$.

With this cut-off implicit, the identity
$\beta_-\big(\eta/\|\eta\|\big)+\beta_+\big(\eta/\|\eta\|\big)=1$ in (\ref{eqn:local-representation-V})
leads to
a decomposition
$V(\tau)=V(\tau)_-+V(\tau)_+$.
Clearly, $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)_-$ is
obtained by replacing $V(\tau)$ by $V(\tau)_-$ in the previous construction.

The wave front of $V(\tau)$ is locally parametrized by
$$
\left(x',\eta\right)\mapsto \left(x',\frac{\partial \varphi}{\partial x}\left(\tau,x',\eta\right),
\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),-\eta\right).
$$
Let us write $\eta=\big(\eta_0,\eta_1\big)\in \mathbb{R}^{2\mathrm{d}+1}
\cong\mathbb{R}\times \mathbb{R}^{2\mathrm{d}}$.
By construction the wave front of $V(\tau)_-$ satisfies
\begin{eqnarray}
\label{eqn:wave-front-of-V-a}
\lefteqn{\mathrm{WF}\left(V(\tau)_-\right)}\\
&\subseteq &\left\{\left(x',\frac{\partial \varphi}{\partial x}\left(\tau,x',\eta\right),
\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),-\eta\right):x'\in X_1,\,
\eta\in \mathbb{R}^{2\mathrm{d}+1}:
\eta_0/\|\eta\|<\delta_2\right\}\nonumber\\
&=&\left\{\left(\Phi_\tau\left(\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),\eta\right),
\left(\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),-\eta\right)\right):x'\in X_1,\,
\eta\in \mathbb{R}^{2\mathrm{d}+1}:
\eta_0/\|\eta\|<\delta_2\right\}.\nonumber
\end{eqnarray}
On the other hand, the wave front of $\Pi$ is
\begin{equation}
\label{eqn:wave-front-of-Pi}
  \mathrm{WF}(\Pi) =\big\{(x,r\alpha_x,x,-r\alpha_x):x\in X,\,r>0\big\}.
\end{equation}

By definition of Heisenberg local coordinates, the cotangent vector $(x,\alpha_x)$
corresponds to $\Big((0,\mathbf{0}),(1,\mathbf{0})\Big)$.
Therefore, if $X_1$ has been chosen sufficiently small, there exists $a>0$
such that (mixing intrinsic notation and expressions in local coordinates)
$\mathrm{dist}_X\left(\left (\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),\eta\right),(x,\alpha_x)\right)>a$,
if $x\in X_1$, $\|\eta\|=1$ and
$\eta_0<\delta_2$, where $\mathrm{dist}_X$ is the Riemannian distance on $X$.
Therefore, if $\epsilon>0$ is sufficiently small we also have
$\mathrm{dist}_X\left(\Phi_\tau\left (\frac{\partial \varphi}{\partial \eta}\left(\tau,x',\eta\right),\eta\right),(x,\alpha_x)\right)>\frac a2$ in the same range if $|\tau|<\epsilon$.

Given this, (\ref{eqn:wave-front-of-V-a}), and (\ref{eqn:wave-front-of-Pi}) we have
$\mathrm{WF}'\big(V(\tau)_-\circ \Pi\big)=\mathrm{WF}'\big(V(\tau)_-\big)\circ \mathrm{WF}'\big(\Pi\big)=\emptyset$;
here $\mathrm{WF}'$ is the image of $\mathrm{WF}$ under the map
$(x,\xi,y,\eta)\mapsto (x,\xi,y,-\eta)$. Hence, $V(\tau)_-\circ \Pi$
is a smoothing operator for every $\tau\in (-\epsilon,\epsilon)$. The claim
follows by arguing again as for (\ref{eqn:rapid-decay}).

\hfill Q.E.D.

\medskip

Given (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-riscalato})
and Lemma \ref{lem:reduction-on-sphere}, we have
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)}\\
&\sim&\frac{\lambda^{2\mathrm{d}+2}}{(2\pi)^{2\mathrm{d}+1}}\,\int _{-a}^a
\int_{B_{2\mathrm{d}}(\mathbf{0},\delta)}\int _0^{+\infty}\int_{-\epsilon}^\epsilon \int _{1/D}^D
\int_{S_+}
e^{i\,\lambda\,\Phi_2}\,r^{2\mathrm{d}}\,G\,d\theta\,d\mathbf{v}\,dt\,d\tau\,dr\,d\omega,\nonumber
\end{eqnarray}
where $\Phi_2$ is as in (\ref{eqn:phase-phi-heis}), and
\begin{equation}
G=:A\Big(x,x+(\theta,\mathbf{v}),\lambda t,\tau,\lambda r\omega\Big)\, F(r)\,\beta_+(\omega)\,
\mathcal{V}(\theta,\mathbf{v}).
\end{equation}

Our next task is to show that integration in $dt$ in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg})
may also be restricted to a suitable compact interval in $(0,+\infty)$.

To this end, suppose $D'\gg 0$ and let $\rho :\mathbb{R}\rightarrow
\mathbb{R}$ be a $\mathcal{C}^\infty$ function satisfying $\rho \ge
0$, $\rho (t)=1$ if $|t|<2\, D'$, $\rho (t)=0$ if $|t|>3\,D'$.
Using the identity $G=\rho\,G+(1-\rho)\,G$ in
(\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg}),
we get
\begin{equation}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-riscalato-tagliato-t}
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim S_{\chi \,e^{-i\lambda
\,(\cdot)}}^{(1)} (x,x)+S_{\chi \,e^{-i\lambda \,(\cdot)}}^{(2)}
(x,x),
\end{equation}
where in $S_{\chi \,e^{-i\lambda \,(\cdot)}}^{(1)} (x,x)$
(respectively, $S_{\chi \,e^{-i\lambda \,(\cdot)}}^{(2)} (x,x)$) $G$ is replaced by
$\rho \,G$ (respectively, by $(1-\rho)\,G$).
Thus, $\int _0^{+\infty}dt$ may be replaced by $\int _0^{3D'} dt$ in
$S_{\chi \,e^{-i\lambda \,(\cdot)}}^{(1)} (x,x)$, and by $\int
_{2D'}^{+\infty}dt$ in $S_{\chi \,e^{-i\lambda \,(\cdot)}}^{(2)} (x,x)$,
respectively.

Recall that $B_{2\mathrm{d}}(\mathbf{0},\delta)$ is the image
of the preferred local chart on $M$ centered at $\pi(x)$
that underlies the chosen Heisenberg local chart on $X$.
Thus $\|\mathbf{v}\|<\delta$ for any $z=x+(\theta,\mathbf{v})\in X_1$.

\begin{lem}
\label{lem:second-term-S-decays}
If $0<\delta\ll 1$ and $D'\gg 0$,
then
$
S_{\chi\,e^{-i\lambda \,(\cdot)}}^{(2)} (x,x)
=O\left(\lambda^{-\infty}\right)
$
as $\lambda\rightarrow +\infty$.
\end{lem}

\textit{Proof.} By construction,
$d_{(x,x)}\psi=\left(\alpha_x,-\alpha_x\right)$ for any $x\in X$
\cite{bs}, \cite{z}, and more generally, $
d_{(e^{i\theta}x,x)}\psi=\left(e^{i\theta}\,\alpha_{e^{i\theta}x},-e^{-i\theta}\,\alpha_{x}\right),
$
for any $x\in X$ and
$e^{i\theta}\in S^1$.

By construction, $(e^{i\theta}x,x)=\big(x+(\theta,\mathbf{0}),x\big)$ in local coordinates.
Thus, if
$\Upsilon:(\theta,\mathbf{v})\mapsto
\psi\big(x+\left(\theta,\mathbf{v}\right),x\big)$,
then $\partial
_\theta\Upsilon\big(\theta,\mathbf{0}\big)=
\,e^{i\theta}$ for every $\theta$. Hence, if $\delta$ is
sufficiently small then
\begin{equation}
\label{eqn:bound-on-dpsi}
2\ge \left|\partial
_\theta\Upsilon\big(\theta,\mathbf{v}\big)\right
| \ge \frac 12,
\end{equation}
for every $z=x+(\theta,\mathbf{v})\in X_1$.

Suppose $D'>4D$; thus, $t\ge 2D'>8D$ on $\mathrm{supp}\{1-\rho\}$.
In view of (\ref{eqn:phase-phi-heis}) and the lower bound in
(\ref{eqn:bound-on-dpsi}), since $|\omega_0|\le 1$ and $r\le D$ we have
\begin{eqnarray}
\Big|\partial _\theta
\Phi_2\big(x,\theta,\mathbf{v},t,\tau,r,\omega\big)\Big|&=&\left|t\partial
_\theta\Upsilon\Big(\theta,\mathbf{v}\Big)-r\,\omega_0\right|\\
&\ge &\frac{t}{2}-r\ge\frac{t}{4}+\left(\frac{t}{4}-D\right)\ge \frac{t}{4}+D\nonumber
\end{eqnarray}
for all $t\in \mathrm{supp}\{1-\rho\}$.
The statement follows by integrating by parts in $\theta$.

\hfill Q.E.D.

\medskip

Thus $\int _0^{+\infty}dt$ may be replaced by
$\int _0^{3D'}dt$ in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg}),
and the amplitude $G$ by $\rho\,G$.
We can similarly discard the contribution to $\int _0^{+\infty}dt$
coming from $(0,\varepsilon)$, if $0<\varepsilon\ll 1$.
Let $\upsilon :\mathbb{R}\rightarrow \mathbb{R}$
be $\mathcal{C}^\infty$, non-negative, and such that
$\upsilon (t)=0$ for $t\le \varepsilon$, $\upsilon (t)=1$ for $t>2\varepsilon$.
Thus
$S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)'+S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)''$,
where $\rho\,G$ has been replaced by $\upsilon \,\rho\,G$ in $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)'$, and
by $(1-\upsilon)\,\rho\,G$ in $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)''$, respectively.

\begin{lem}
\label{lem:reduction-in-t-away-from-0}
If $\varepsilon$ is sufficiently small, then
$S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)''=O\left(\lambda^{-\infty}\right)$ for
$\lambda\rightarrow +\infty$.
\end{lem}

In the following, recall that $\mathrm{supp}(\chi)\subseteq (-\epsilon,\epsilon)$, and that $\epsilon$ has
been chosen suitably small.

\medskip

\textit{Proof.}
Given (\ref{eqn:phase-phi-heis}),
$\partial _\tau\Phi_2=r\cdot \Big [q\big(x,\omega\big)+O(\tau)\Big]-1$.
Let $M$ be the maximum reached by $q$ on the unit sphere bundle of
$X$, and suppose $r<1/(2M)$.
Then $\Big|\partial _\tau\Phi_2\Big|\ge 1-\frac{M+O(\epsilon)}{2M}\ge \frac 13$,
since $\epsilon$ is small.
Integration by parts in $\tau$ shows that only a rapidly decreasing
contribution in $\lambda$ is lost if we multiply the amplitude by a  cut-off in
$r$ so as to assume $r>1/(2M)$.

Now suppose that $\varepsilon<\delta_1/(16\,M)$, and argue as in Lemma \ref{lem:second-term-S-decays},
using the upper bound in
(\ref{eqn:bound-on-dpsi}), and the fact that $\omega_0>\delta_1$ on $S_+$: if $t<2\varepsilon$,
\begin{eqnarray*}
\Big|\partial _\theta
\Phi_2\big(x,\theta,\mathbf{v},t,\tau,r,\omega\big)\Big|&=&\left|t\partial
_\theta\Upsilon\Big(\theta,\mathbf{v}\Big)-r\,\omega_0\right|\\
&\ge &r\omega_0-2t\ge \frac{\delta_1}{2M}-4\varepsilon\ge \frac{\delta_1}{4M}.\nonumber
\end{eqnarray*}
The statement follows by integration by parts in $\theta$.

\hfill Q.E.D.

\medskip

On the upshot, for $\lambda\rightarrow +\infty$
we only lose a rapidly decreasing term in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg})
by replacing $\int_0^{+\infty}dt$ by $\int_{1/C}^{C}dt$ for some $C\gg 0$, and
$G$ by $\zeta (t)\,G $, where $\zeta \in \mathcal{C}^\infty (\mathbb{R})$
is supported in $(1/C,C)$, and identically one on
$(2/C,C/2)$.

As a further reduction, let $\gamma :\mathbb{C}^\mathrm{d}\rightarrow \mathbb{R}$ be
$\mathcal{C}^\infty$, nonnegative,
such that $\gamma (\mathbf{v})=1$ if $\|\mathbf{v}\|<2$ and
$\gamma (\mathbf{v})=0$ if $\|v\|>3$. Let us set
$\gamma_\lambda(\mathbf{v})=:\gamma \left(\lambda^{1/3}\mathbf{v}\right)$
($\lambda >0,\,\mathbf{v}\in \mathbb{C}^\mathrm{d}$).
We obtain $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)\sim S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)^{\prime}_1+
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)^{\prime}_2$, where the amplitude is $\gamma_\lambda \,G$
in $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)^{\prime}_1$,
and $(1-\gamma_\lambda )\,G$ in $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)^{\prime}_2$,
respectively.

\begin{lem}
\label{lem:riscalato-in-v-due-termini-decade}
$S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)^{\prime}_2=O\left(\lambda^{-\infty}\right)$ for
$\lambda\rightarrow +\infty$.
\end{lem}

\textit{Proof.} Let $\mathrm{dist}_X$ and $\mathrm{dist}_M$ be the Riemannian distance functions on
$X$ and $M$, respectively. If $m=:\pi(x)$, let
$m+\mathbf{v}$ be the point in $M$ with preferred coordinates $\mathbf{v}$.
By construction \cite{sz}, $m+\mathbf{v}=\pi\big(x+(\theta,\mathbf{v})\big)$
for every $\theta,\,\mathbf{v}$.

If $\mathbf{v}\in \mathrm{supp}\{1-\gamma_\lambda\}\cap B_{2\mathrm{d}}(\mathbf{0},\delta)$, then
$\delta>\|\mathbf{v}\|>2\,\lambda^{-1/3}$. By definition, in the preferred local coordinates
centered at $m$ the tangent space $T_mM$ is unitarily identified with $\mathbb{C}^\mathrm{d}$.
Therefore,
if $\delta$ is sufficiently small we may assume that
$\mathrm{dist}_M(m+\mathbf{v},m)\ge \frac 12\,\|\mathbf{v}\|$,
$\forall\,\mathbf{v}\in B_{2\mathrm{d}}(\mathbf{0},\delta)$. On the upshot,
\begin{eqnarray}
\label{eqn:bound-on-distance}
\mathrm{dist}_X\big(x+(\theta,\mathbf{v}),x\big)\ge\mathrm{dist}_M\big(m+\mathbf{v},m\big)\ge
\lambda^{-1/3},
\end{eqnarray}
$\forall \,\mathbf{v}\in \mathrm{supp}\{1-\gamma_\lambda\}\cap B_{2\mathrm{d}}(\mathbf{0},\delta),\,
\theta\in (-a,a)$.

By (1.4) of \cite{bs},
$\Im \psi (x,y)\ge C\,\mathrm{dist}_X(x,y)^2$ $\forall\,x,y\in X$ and a suitable
constant $C>0$.
In view of (\ref{eqn:bound-on-distance}),
we conclude that on the same support
\begin{eqnarray}
\label{eqn:bound-on-d-t-psi}
\Big|\partial _t\Phi_2\big(x,\theta,\mathbf{v},t,\tau,r,\omega\big)\Big|&=&\Big |\psi \big(x+(\theta,\mathbf{v}),x\big)\Big|
\nonumber \\
&\ge& \Im \psi \big(x+(\theta,\mathbf{v}),x\big)\ge C\lambda ^{-2/3}.
\end{eqnarray}
The statement then follows by integrating by parts in $t$, since each integration introduces a factor $\frac 1\lambda\cdot
\lambda ^{2/3}=\lambda ^{-1/3}$; notice that $t$ appears in the amplitude by a factor which has an asymptotic expansion
in $\lambda t$, with leading terms $\lambda ^\mathrm{d}t^\mathrm{d}$. Therefore, every
$t$-derivative of the amplitude
remains $O\left(\lambda ^\mathrm{d}\right)$.

\hfill Q.E.D.

\medskip

Lemma \ref{lem:riscalato-in-v-due-termini-decade} means that we only lose a rapidly decreasing term in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-solo-fourier-heisenberg}) if, in addition to the previous reductions,
$\int _{B_{2\mathrm{d}}(\mathbf{0},\delta)}d\mathbf{v}$ is replaced by
$\int _{B_{2\mathrm{d}}(\mathbf{0},3\lambda^{-1/3})}d\mathbf{v}$, and the amplitude is multiplied by
$\gamma_\lambda(\mathbf{v})$.

Let us introduce the change of variables $\mathbf{v}\rightsquigarrow \mathbf{v}/(r\sqrt{\lambda})$.
Since $D^{-1}\le r\le D$, integration in $d\mathbf{v}$ is now over a ball of radius $O\left(\lambda^{1/6}\right)$.
We obtain:
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-lots-of-reductions}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)}\\
&\sim&\frac{\lambda^{\mathrm{d}+2}}{(2\pi)^{2\mathrm{d}+1}}\,\int _{-a}^a
\int_{\mathbb{C}^\mathrm{d}}\int _{1/C}^{C}\int_{-\epsilon}^\epsilon \int _{1/D}^D
\int_{S_+}
e^{i\,\lambda\,\Phi_3}\,H\,d\theta\,d\mathbf{v}\,dt\,d\tau\,dr\,d\omega,\nonumber
\end{eqnarray}
where
\begin{eqnarray*}
\lefteqn{\Phi_3=:\Phi_2\left(x,\theta,\frac{\mathbf{v}}{r\sqrt{\lambda}},t,\tau,r,\omega\right)}\\
&=&-\frac{\mathbf{v}}{\sqrt{\lambda}}\cdot\omega_1
+t\psi\left(x+\left(\theta,\frac{\mathbf{v}}{r\sqrt{\lambda}}\right),x\right)-r\,\theta \,\omega_0
+r\tau\,q\big(x,\omega\big)-\tau+O\left(\tau^2\right)\,r
\end{eqnarray*}
and
$$
H=:\gamma\left(\frac{\mathbf{v}}{r\lambda^{1/6}}\right)\,
\zeta(t)\,G\left(x,\theta,\frac{\mathbf{v}}{r\sqrt{\lambda}},t,\tau,r,\omega\right).
$$
In view of (65) in \cite{sz}, we have
\begin{eqnarray}
\lefteqn{t\psi\left(x+\left(\theta,\frac{\mathbf{v}}{r\sqrt{\lambda}}\right),x\right)}\nonumber\\
&=&
it\,\left[1-e^{i\theta}\right]+\left[\frac{it}{2r^2\lambda}\,\|\mathbf{v}\|^2\,e^{i\theta}
+t\,R_3^\psi\left(\frac{\mathbf{v}}{r\sqrt{\lambda}}\right)\,e^{i\theta}\right],
\end{eqnarray}
where $R_3^\psi$ vanishes to third order at the origin.

The latter term gives rise to a bounded exponential for $\|\mathbf{v}\|\lesssim \lambda ^{1/6}$;
as in \cite{sz}, this exponential may then be incorporated as part of the amplitude in an appopriate oscillatory
integral.
More precisely, we have:
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-lots-and-lots-of-reductions}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)}\\
&\sim&\frac{\lambda^{\mathrm{d}+2}}{(2\pi)^{2\mathrm{d}+1}}\,\int_{\mathbb{C}^\mathrm{d}}\int_{S_+}
e^{-i\sqrt{\lambda}\,\mathbf{v}\cdot \omega_1}\cdot \left[
\int _{-a}^a
\int _{1/C}^{C}\int_{-\epsilon}^\epsilon \int _{1/D}^D
e^{i\,\lambda\,\Psi}\,\widetilde{H}\,d\theta\,dt\,d\tau\,dr\right]\,d\mathbf{v}\,d\omega,\nonumber
\end{eqnarray}
where
\begin{equation}
\Psi=:it\,\left[1-e^{i\theta}\right]-r\,\theta \,\omega_0
+r\tau\,q\big(x,\omega\big)-\tau+O\left(\tau^2\right)\,r,
\end{equation}
\begin{equation}
\widetilde{H}=:e^{-(t/2r^2)\,\|\mathbf{v}\|^2\,e^{i\theta}
+i\lambda\,t\,R_3^\psi\left(\frac{\mathbf{v}}{r\sqrt{\lambda}}\right)\,e^{i\theta}}\,H.
\end{equation}

Let us first consider the asymptotics of the inner integral in
(\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-lots-and-lots-of-reductions}),
an oscillatory integral $d\theta\,dt\,d\tau\,dr$ with phase $\Psi$, and amplitude
$\widetilde{H}$ that may be expanded in descending powers of $\lambda^{-1/2}$
for $\lambda\rightarrow +\infty$; here
$\mathbf{v}\in \mathbb{C}^\mathrm{d}$, $\omega\in S_+$ play the role of
parameters.
Using the the ellipticity of $q$ and the fact that $|\tau|<\epsilon$, $\epsilon$ small,
a straightforward computation yields the following.

\begin{lem}
\label{lem:stationary-point}
On the given domain of integration, $\Psi$ has a unique stationary point
given by
$(\theta_0,t_0,\tau_0,r_0)=\big(0,\omega_0/q(x,\omega),0,1/q(x,\omega)\big)$.
Furthermore, if $\Psi''$ is the Hessian of $\Psi$ at this critical point then
$$
\det\left(\frac{\lambda \Psi''}{2\pi\,i}\right)=\left(\frac{\lambda}{2\pi}\right)^4\,q(x,\omega)^2.
$$
\end{lem}

We are thus in a position to apply the stationary phase Lemma, following \S 5 of \cite{sz},
and conclude that for any integer $N\gg 0$ the inner integral in
(\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-lots-and-lots-of-reductions}) is given by
$S_N(x,\omega,\mathbf{v})+R_N(x,\omega,\mathbf{v})$, where
\begin{eqnarray}
\label{eqn:inner-integral-expanded}
\lefteqn{S_N(x,\omega,\mathbf{v})=\lambda^{\mathrm{d}-2}\,(2\pi)^{2}\,\chi(0)}\\
&&\cdot \frac{\omega_0^\mathrm{d}}{q(x,\omega)^{\mathrm{d}+1}}\,\,e^{-\frac 12\omega_0\,
q(x,\omega)\,\|\mathbf{v}\|^2}\beta_+(\omega)\,s(x,x)\,\left(1+\sum _{j=1}^NF_j(\omega,\mathbf{v})
\lambda^{-j/2}\right),\nonumber
\end{eqnarray}
where $F_j$ is polynomial in $\mathbf{v}$ and
$\big|R_N(x,\omega,\mathbf{v})\big|\le C_N\,\lambda^{\mathrm{d}-2-(N+1)/2}\,e^{-a\,\|\mathbf{v}\|^2}$, for some $a>0$.
Therefore, the asymptotic expansion may be integrated term by term in $d\mathbf{v}\,d\omega$.

In order to apply the stationary phase Lemma again, we set $\mu=\sqrt{\lambda}$. Since
$s_0(x,x)=\pi^{-\mathrm{d}}$, we have
\begin{eqnarray}
\label{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)}\\
&\sim&\frac{2\pi}{\pi^\mathrm{d}} \,\left(\frac{\mu^2}{2\pi}\right)^{2\mathrm{d}}\,\chi (0)\cdot \left(
\int_{\mathbb{C}^\mathrm{d}}\int_{S_+} e^{-i\mu\,\mathbf{v}\cdot \omega_1-\frac 12\omega_0\,
q(x,\omega)\,\|\mathbf{v}\|^2} \,\frac{\omega_0^\mathrm{d}}{q(x,\omega)^{\mathrm{d}+1}}
\,\beta_+(\omega)\,d\mathbf{v}\,d\omega \right.   \nonumber\\
&&\left.+\sum _{j=1}^{+\infty}\mu^{-j}\int_{\mathbb{C}^\mathrm{d}}\int_{S_+}
e^{-i\mu\,\mathbf{v}\cdot \omega_1-\frac 12\omega_0\,
q(x,\omega)\,\|\mathbf{v}\|^2} \,\frac{\omega_0^\mathrm{d}}{q(x,\omega)^{\mathrm{d}+1}}
\,\beta_+(\omega)\,F_j(\omega,\mathbf{v})
\,d\mathbf{v}\,d\omega \right).\nonumber
\end{eqnarray}

To obtain an asymptotic expansion for $S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)$ as $\lambda\rightarrow +\infty$,
we are thus reduced to
determining asymptotic expansions for each of
the summands in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}) as $\mu\rightarrow +\infty$.
To this end, we remark that $\omega_1$ is a system of local coordinates on $S_+\subseteq S^{2\mathrm{d}}$, and
$\omega_0=\left(1-\|\omega_1\|^2\right)^{1/2}$.
Each summand in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}) is an oscilatory integral whose phase is the quadratic form $(\mathbf{v},\omega_1)\mapsto -\mathbf{v}\cdot \omega_1$
on $\mathbb{R}^{\mathrm{4d}}\cong \mathbb{R}^{2\mathrm{d}}\times \mathbb{R}^{2\mathrm{d}}$.
The origin $\mathbf{v}=\mathbf{0},\,\omega_1=\mathbf{0}$ is the only stationary point, and the Hessian has determinant
one. Thus, the stationary phase Lemma applies and as $\mu\rightarrow +\infty$ we are left with an asymptotic expansion
\begin{eqnarray}
\label{eqn:asympt-exp-1st-summand}
\lefteqn{S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)}\\
&\sim&
\frac{2\pi}{\pi^\mathrm{d}} \,\left(\frac{\mu^2}{2\pi}\right)^{2\mathrm{d}}\,\chi (0)\cdot \left(\frac{2\pi}{\mu}\right)
^{2\mathrm{d}}\,q\big(x,(1,\mathbf{0})\big)^{-(\mathrm{d}+1)}\cdot \left(1+\sum_{j\ge 1}c_j\,\mu^{-j}\right)\nonumber\\
&=&2\pi\,\left(\frac{\lambda}{\pi}\right)^{\mathrm{d}}\,\chi (0)\,
\varsigma_T(x)^{-(\mathrm{d}+1)}\cdot \left(1+\sum_{j\ge 1}c_j\,\lambda^{-j/2}\right).\nonumber
\end{eqnarray}

\begin{lem}
\label{lem:c1vanishes}
$c_1=0$.
\end{lem}

\textit{Proof.} $c_1$ is the coefficient of $\mu^{2\mathrm{d}-1}$ in the second line of (\ref{eqn:asympt-exp-1st-summand}).
Looking at (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}), there are two possible contributions to
this coefficient.

One is the second term in the asymptotic expansion of the first summand in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}); by the stationary phase Lemma, this is multiple of
$$
\langle D_\mathbf{v},D_{\omega_1}\rangle \left(e^{-\frac 12\omega_0\,
q(x,\omega)\,\|\mathbf{v}\|^2} \,\frac{\omega_0^\mathrm{d}}{q(x,\omega)^{\mathrm{d}+1}}\,\beta_+(\omega)\right)(\mathbf{0},\mathbf{0})=0.
$$

The other possible contribution comes from the leading term in the asymptotic expansion of the second
summand in (\ref{eqn:s-chi-lambda-diag-part-reduced-P-2-stationary-phase-mu}), the one with $j=1$;
hence it is a multiple of the
evaluation at the critical point of the amplitude of that oscillatory integral. On the other hand, the appearance
of the fractional powers of $\lambda^{-1}$ in the amplitude originates from the Taylor expansion of $\widetilde{H}$
in $\mathbf{v}/(r\sqrt{\lambda})$; hence $F_1$ is linear in $\mathbf{v}$, and so it vanishes at the critical point.

\hfill Q.E.D.

\begin{cor}
\label{cor:summingup}
There exists $\epsilon>0$ such that for every $\chi \in \mathcal{C}^\infty_0\big((-\epsilon,\epsilon)\big)$

$$
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)=\left\{
\begin{array}{lcc}
  2\pi\,\left(\frac{\lambda}{\pi}\right)^{\mathrm{d}}\,\chi (0)\,\varsigma_T(x)^{-(\mathrm{d}+1)}
+O\left(\lambda^{\mathrm{d}-1}\right) & \mathrm{for} & \lambda\rightarrow +\infty \\
   &  &  \\
  O\left(\lambda^{-\infty}\right) & \mathrm{for} & \lambda\rightarrow -\infty
\end{array}\right.
.
$$
\end{cor}

To complete the proof of Theorem \ref{thm:local-weyl-law}, we shall now closely follow \S 12 of \cite{gs},
hence the exposition will be rather sketchy.

Le us fix once and forall $\epsilon>0$ such that the conclusion of Corollary \ref{cor:summingup}
holds true. One can see that $S_{\chi \,e^{-i\lambda \,(\cdot)}}$ is a smoothing operator, with kernel
$$
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,y)=\sum _j\widehat{\chi}(\lambda-\lambda_j)\,e_j(x)\,\overline{e_j(y)}
\,\,\,\,\,\,\,\,(x,y\in X),
$$
hence the Corollary implies that for $\lambda\rightarrow +\infty$
$$
\sum _j\widehat{\chi}(\lambda-\lambda_j)\,\big|e_j(x)\big|^2=
2\pi\,\left(\frac{\lambda}{\pi}\right)^{\mathrm{d}}\,\chi (0)\,\varsigma_T(x)^{-(\mathrm{d}+1)}
+O\left(\lambda^{\mathrm{d}-1}\right).
$$

Applying this to  $\chi \in \mathcal{C}^\infty_0\big((-\epsilon,\epsilon)\big)$ such that
$\chi (0)=1$ and $\widehat{\chi}\ge 0$ shows the following:

\begin{cor}
\label{cor:incremental-estimate}
There exists a constant
$E>0$ such that
$\mathcal{T}(\lambda+1,x)-\mathcal{T}(\lambda,x)\le E\,\lambda^{\mathrm{d}}$ for $\lambda\rightarrow +\infty$,
whence
$\mathcal{T}(\lambda,x)=O\left(\lambda^{\mathrm{d}+1}\right)$. These estimates are uniform in $x\in X$.
\end{cor}

Let us define measures $d\mathcal{T}_x=:\sum _j\big|e_j(x)\big|^2\,\delta_{\lambda_j}$ ($x\in X$)
on $\mathbb{R}$, so that
\begin{equation}
\label{eqn:measure-integrated}
\mathcal{T}(\lambda,x)=\int_{-\infty}^\lambda d\mathcal{T}_x(\eta),\,\,\,\,\,\,\,\,
S_{\chi \,e^{-i\lambda \,(\cdot)}} (x,x)=\int_{-\infty}^{+\infty}\widehat{\chi}(\lambda-\eta)\,d\mathcal{T}_x(\eta).
\end{equation}

Under the same assumption on $\chi$, set $G(\lambda)=:\int_{-\infty}^\lambda \widehat{\chi}(\tau)\,d\tau$; the integral
$\int _{-\infty}^{+\infty}G(\lambda-\eta)\,d\mathcal{T}_x(\eta)$ may be computed in two differen manners.

On the hand, we use the change of variable $\tau \rightsquigarrow \tau-\eta$ and the Tonelli-Fubini theorem
(recall that $\widehat{\chi}\ge 0$) to change the order of integration and obtain
$\int _{-\infty}^{+\infty}G(\lambda-\eta)\,d\mathcal{T}_x(\eta)=
\int_{-\infty}^\lambda S_{\chi \,e^{-i\tau \,(\cdot)}} (x,x)\,d\tau$. Using Corollary \ref{cor:summingup},
we get
\begin{equation}
\label{eqn:1stcomputation-of-int}
\int _{-\infty}^{+\infty}G(\lambda-\eta)\,d\mathcal{T}_x(\eta)
=\frac{2\,\pi^2}{(\mathrm{d}+1)}\cdot
\left(\frac{\lambda}{\pi\,\varsigma_T(x)}\right)^{\mathrm{d}+1}
+O\left(\lambda^{d}\right).
\end{equation}

On the other hand, if $H$ is the Heaviside function, we have the chain of equalities:
\begin{eqnarray}
\label{eqn:2ndcomputation-of-int}
\lefteqn{\int _{-\infty}^{+\infty}G(\lambda-\eta)\,d\mathcal{T}_x(\eta)=\sum_jG(\lambda-\lambda_j)\,\big|e_j(x)\big|^2}\\
&=&\sum _j\left(\int _{-\infty}^{\lambda-\lambda_j}\widehat{\chi}(\tau)\,d\tau\right)\cdot \big|e_j(x)\big|^2
=\sum_j\int _{-\infty}^{+\infty}H(\lambda-\lambda_j-\tau)\,\widehat{\chi}(\tau)
\cdot \big|e_j(x)\big|^2\,d\tau\nonumber\\
&=&\int _{-\infty}^{+\infty}\left(\sum_jH(\lambda-\lambda_j-\tau)\,
\, \big|e_j(x)\big|^2\right)\cdot \widehat{\chi}(\tau)\,d\tau =
\int _{-\infty}^{+\infty}\mathcal{T}(\lambda-\tau,x)\cdot \widehat{\chi}(\tau)\,d\tau \nonumber\\
&=&\mathcal{T}(\lambda,x)\cdot \int _{-\infty}^{+\infty}\widehat{\chi}(\tau)\,d\tau +
\int _{-\infty}^{+\infty}\Big(\mathcal{T}(\lambda-\tau,x)-\mathcal{T}(\lambda,x)\Big)\cdot \widehat{\chi}(\tau)\,d\tau\nonumber\\
&=&2\pi\,\mathcal{T}(\lambda,x)+O\left(\lambda^{\mathrm{d}}\right),\nonumber
\end{eqnarray}
where the estimate on the second summand is a consequence of Corollary \ref{cor:incremental-estimate} and the fact that
$\widehat{\chi}$ is of rapid decay.

Theorem \ref{thm:local-weyl-law} follows from (\ref{eqn:1stcomputation-of-int}) and
(\ref{eqn:2ndcomputation-of-int}).

\hfill Q.E.D.

\section{Proof of Corollary \ref{cor:global-weyl-law}}

In order to prove Corollary \ref{cor:global-weyl-law}, let us dwell on the symplectic structure
$\Omega_\Sigma$ of the cone $\Sigma$
induced by the canonical symplectic structure of $T^*X$.

The action $r$ of $S^1$ on $X$ lifts to an action on $\Sigma$, given by
$e^{i\theta}\cdot (x,\alpha_x)=:\left(e^{i\theta}\cdot x,\alpha_{e^{i\theta}\cdot x}\right)$;
denote by $\frac{\partial}{\partial \theta}$ the generator of this action.

$\mathbb{R}_+=(0,+\infty)$ also acts on the cone $\Sigma$, by $r\cdot (x,s\,\alpha_x)=:(x,rs\,\alpha_x)$;
let $\frac{\partial}{\partial r}$ be the generator of the latter action; it spans the vertical tangent bundle
of the projection  $p:\Sigma\rightarrow X$.

Thus $\frac{\partial}{\partial \theta}$ and $\frac{\partial}{\partial r}$
are smooth nowhere vanishing vector fields on $\Sigma$, and they are actually everywhere linearly independent.
Let $\mathfrak{V}\subseteq T\Sigma$ be the rank-2 vector subbundle of the tangent bundle of $\Sigma$
spanned by $\frac{\partial}{\partial \theta}$ and $\frac{\partial}{\partial r}$. On $\mathfrak{V}$ let us consider
the symplectic structure $\Omega _\mathfrak{V}$ uniquely determined by
$\Omega _\mathfrak{V}\left( \frac{\partial}{\partial \theta},\frac{\partial}{\partial r}\right)=1$.

For every $r>0$, the map
$\sigma_r:X\rightarrow \Sigma$, $x\mapsto (x,r\alpha_x)$, is a smooth section of $p$, smoothly varying with $r$.
Thus $d_x\sigma_r\big(T_xX\big)\subseteq T_{(x,r\alpha_x)}\Sigma$ is the fiber of a vector subbundle of rank
$2\mathrm{d}+1$ of $T\Sigma$;
therefore, $\mathfrak{H}_{(x,r\alpha_x)}=:d_x\sigma_r\big(\ker(\alpha_x)\big)$
is the fiber of a vector subbundle of rank $2\mathrm{d}$ of $T\Sigma$, isomorphic to the pull-back $q^*\big(TM\big)$,
where $q=:\pi\circ p:\Sigma\rightarrow M$. We have $T\Sigma=\mathfrak{H}\oplus \mathfrak{V}$.

On $\mathfrak{H}\cong q^*(TM)$ we have the symplectic structure $q^*(\omega)$; let $\Omega_\mathfrak{H}$ be the symplectic
structure on $\mathfrak{H}$ defined by
$\left(\Omega_\mathfrak{H}\right)_{(x,r\alpha_x)}=2\,r\,q^*(\omega)_{(x,r\alpha_x)}$.

\begin{lem}
\label{lem:symp-str-cone}
$\big(T\Sigma,\Omega_\Sigma)\cong \big(\mathfrak{H},\Omega_\mathfrak{H}\big)\oplus \big(\mathfrak{V},\Omega_\mathfrak{V}\big)$
as symplectic vector bundles.
\end{lem}

\textit{Proof.} Suppose $x\in X$, and let $U\subseteq X$ be an open neighborhood of $x$,
on which a local coordinate chart $\mathbf{t}=(t_i):U\rightarrow B$ is defined;
here $B$ is an open ball in $\mathbb{R}^{2\mathrm{d}+1}$, and $\mathbf{t}$ is a diffeomorfism.

Using the basis $\{dt_i\}_i$,
a differential 1-form $\widetilde{a}=\sum_i a_i\,dt_i$ on $U$ may be represented
as a smooth function $a:B\rightarrow \mathbb{R}^{2\mathrm{d}+1}$, $\mathbf{t}\mapsto
\big(a_i(\mathbf{t})\big)$.

Using the basis $\{dt_i\wedge dt_j\}_{i<j}$,
a differential 2-form $\widehat{\nu}=\sum _{i<j}\nu _{ij}\,dt_i\wedge dt_j$ may be represented
as a smooth function $\nu:B\rightarrow A_{2\mathrm{d}+1}$,
$y\mapsto \nu(\mathbf{t})=\big[\nu_{ij}(\mathbf{t})\big]$; here
$A_{2\mathrm{d}+1}$ is the vector space of skew-symmetric $(2\mathrm{d}+1)\times (2\mathrm{d}+1)$ real matrices.
For any $y\in U$, $\nu(y)$ is the matrix representing $\widehat{\nu}_y$
as a skew-symmetric bilinear pairing on $T_yX$.

If $g:B\rightarrow \mathbb{R}^{2\mathrm{d}+1}$ is smooth, let $\mathrm{Jac}(g)$ be the its Jacobian matrix;
then $d\widetilde{g}=\widehat{S(g)}$, where $S(g)=:\mathrm{Jac}(g)^t-\mathrm{Jac}(g)$.

If $\alpha =\widetilde{a}$ on $U$, and $y\in X$ has local coordinates $\mathbf{t}\in B$,
the horizontal subspace $\ker (\alpha_y)\subseteq T_yX$
corresponds to the hyperplane $\big\{\mathbf{v}:a(\mathbf{t})\cdot \mathbf{v}
=0\big\}\subseteq \mathbb{R}^{2\mathrm{d}+1}$. Let $\kappa (\mathbf{t})\in \mathbb{R}^{2\mathrm{d}+1}$
correspond to $\left.\frac{\partial}{\partial \theta}\right|_y$.

On $\mathbb{R}^{2\mathrm{d}+1}\times \mathbb{R}^{2\mathrm{d}+1}$, with linear coordinates $(\mathbf{q},\mathbf{p})$,
consider the symplectic structure $d\mathbf{p}\wedge d\mathbf{q}$; the local chart $\mathbf{t}$ determines
a the symplectic coordinate chart $T^*(U)\cong  B\times \mathbb{R}^{2\mathrm{d}+1}$.
Locally in this chart, $\Sigma$ is parametrized as $\big\{\big(\mathbf{t},r\,a(\mathbf{t})\big):\mathbf{t}\in B,\,r>0\big\}$.
Therefore, if $r>0$ and
$z=:(y,r\,\alpha_y)$ then the tangent space $T_{z}\Sigma$
corresponds to the vector subspace
\begin{eqnarray}
\label{eqn:tangent-space-sigma}
\lefteqn{\left\{\left(
                  \begin{array}{c}
                    \mathbf{v} \\
                    r\,\mathrm{Jac}_\mathbf{t}(a)\,\mathbf{v} \\
                  \end{array}
                \right):\mathbf{v}\in \mathbb{R}^{2\mathrm{d}+1}\right\}\oplus \mathrm{span}\left\{\left(
                  \begin{array}{c}
                    0 \\
                    a(\mathbf{t}) \\
                  \end{array}
                \right)\right\}}\\
                &=&\left\{\left(
                  \begin{array}{c}
                    \mathbf{v} \\
                    r\,\mathrm{Jac}_\mathbf{t}(a)\,\mathbf{v} \\
                  \end{array}
                \right):\mathbf{v}\in \mathbb{R}^{2\mathrm{d}+1},\,a(\mathbf{v})=0\right\}\nonumber \\
                &&\oplus
                \mathrm{span}\left\{\left(
                  \begin{array}{c}
                    \kappa (\mathbf{t})\\
                    r\,\mathrm{Jac}_\mathbf{t}(a)\,\kappa (\mathbf{t})\\
                  \end{array}
                \right),\,\left(
                  \begin{array}{c}
                    0 \\
                    a(\mathbf{t}) \\
                  \end{array}
                \right)\right\}.\nonumber
                \end{eqnarray}

Since the chart is symplectic, the statement follows from (\ref{eqn:tangent-space-sigma}) recalling that
$d\alpha=2\,\omega$.

\hfill Q.E.D.

\medskip

We can now prove Corollary \ref{cor:global-weyl-law}.
Given Lemma \ref{lem:symp-str-cone}, the volume form on $\Sigma$ is
\begin{eqnarray}
\label{eqn:vol-form-cone}
\lefteqn{dV_\Sigma=\frac{1}{(d+1)!}\,\Omega_\Sigma^{\wedge (\mathrm{d}+1)}=
\frac{1}{\mathrm{d!}}\,\Omega_\mathfrak{H}^{\wedge \mathrm{d}}
\wedge \Omega_\mathfrak{V}}\\
&=&2^{\mathrm{d}}\,r^{\mathrm{d}}\,q^*\big(dV_M\big)\wedge \big(d\theta\wedge dr\big)
=2^{\mathrm{d}+1}\,\pi\,r^{\mathrm{d}}\,p^*(d\mu_X\big)\wedge  dr.\nonumber
\end{eqnarray}

By (\ref{eqn:vol-form-cone}), the symplectic volume of $\Sigma_1=\big\{(x,r\,\alpha_x)\in \Sigma: r\le 1/\varsigma_T(x)\big\}$
is
\begin{eqnarray}
\mathrm{Vol}(\Sigma_1)&=&2^{\mathrm{d}+1}\,\pi\,\int _X\left(\int _0^{1/\varsigma_T(x)}r^{\mathrm{d}}\,dr\right)\,d\mu_X\\
&=&2^{\mathrm{d}+1}\,\frac{\pi}{\mathrm{d}+1}\,\int _X\left(\frac{1}{\varsigma_T(x)}\right)^{\mathrm{d}+1}\,d\mu_X.\nonumber
\end{eqnarray}

Given this, Theorem \ref{thm:local-weyl-law} and (\ref{eqn:counting-trace}) imply
\begin{eqnarray}
N_T(\lambda)&=&\left(\frac{\lambda}{\pi}\right)^{\mathrm{d}+1}\,\frac{\pi}{d+1}\,
\int_X\left(\frac{1}{\varsigma_T(x)}\right)^{\mathrm{d}+1}\,d\mu_X
+O\left(\lambda^\mathrm{d}\right)\\
&=&\left(\frac{\lambda}{2\pi}\right)^{\mathrm{d}+1}\,\mathrm{Vol}(\Sigma_1)
+O\left(\lambda^\mathrm{d}\right).\nonumber
\end{eqnarray}

\hfill Q.E.D.


TITLE: When is $\sum_{n=0}^{\infty}\frac{n^k}{3^n}$ an integer?
QUESTION [16 upvotes]: The motivation for this was to find some nice expression for the sum
$$\sum_{n=0}^\infty\frac{P(n)}{b^n}$$
where $P(n)$ is a polynomial and $b$ is a positive integer greater than 1.
Clearly, it suffices to find sums of the form $\sum_{n=0}^\infty\frac{n^k}{b^n}$ for integer k.
If we set $b=2$, a simple translation and subtraction of the terms in the series shows that $f(k)=\sum_{n=0}^\infty\frac{n^k}{2^n}$ satisfies the recurrence relation $f(k)=\sum_{i=0}^{k-1}{k \choose i}f(i)$, producing the series $2,2,6,26,...$ (A076726 in OEIS).
If $b>2$, a similar argument shows that the recurrence is $$f(k)=\sum_{i=0}^{k-1}\frac{{k \choose i}f(i)}{b-1}$$
When b=3, the series is as follows:
$\frac{3}{2},\frac{3}{4},\frac{3}{2},\frac{33}{8},15,\frac{273}{4},\frac{1491}{4},\frac{38001}{16},17295,...$
The series takes on integer values at strange places; the first several such values are 4, 8, 12, 13, 16, 20, 24, 32, ... (not in OEIS)
Some unproved notes about the integer positions: They're usually multiples of 4, and most multiples of 4 seem to be on the list except powers of 2 minus 4. All powers of 2 larger than 2 that I've computed so far (up to 1024) are integers. The odd numbers that take on integer values are frequently separated by powers of 2; if I recall correctly, the next odd number after 13 is 77=13+64.
Curiously, there appear to be $no$ integer values for $b>3$.
Does anyone know why these patterns should occur? If nothing else, more data on these values would be great (the few programming languages I'm competent in weren't able to handle large numbers very well, and I've lost the data up to 1024 that I used to have).

REPLY [1 votes]: I tested whether $f(k) = \sum_{k=0}^\infty n^k \,3^{-n}$ is an integer for $k = 0$ to $50000$.  In this range, $f(k) \in \mathbb{Z}$ if and only if $k$ has one of the following forms (for $j \ge 0$):

$4+4j$ (* with exceptions—see below)
$13+64j$
$122+256j$
$379+1024j$
$1529+2048j$
$8181+32768j$

There are two sequences of exceptions to the first case (i.e., when $k$ is a positive multiple of $4$).  When $k$ has either of the forms

$28+32j$
$1016+1024j$

then $f(k) \notin \mathbb{Z}$.
How does this extend to all $k$?  A natural conjecture is that there is an infinite list of sequences like the ones above, as well as an infinite list of exception sequences to the multiple-of-$4$ case.  We would expect that they all have of a power-of-$2$ stride and that they begin at an integer smaller than their stride (except for the first sequence).


TITLE: If a series is convergent is it possible to perform mathematical operations with its sum?
QUESTION [0 upvotes]: Suppose the series given by $$\left(\displaystyle{\sum_{k=1}^{\infty}|x_{k}|^{2}}\right)^{\frac{1}{2}}$$ is convergent to a number S. Where $x_{k}\in\mathbb{C}$ and $|.|$ is the complex modulus. I want to know if it's possible to know if $\displaystyle{\sum_{k=1}^{\infty}|x_{k}|^{2}}$ is a convergent series. I suppose it is a convergent series since S is a finite number and square a finite number is also finite, am I right?

REPLY [1 votes]: $S$ is not a series; it is the square root of a series. It turns out you are (in a roundabout way) correct though. $S^2$ is the square of the square root of a series, so it is a series itself.
To answer your general question, yes you can perform arithmetic with a convergent series. Like you said, you can treat it as a finite number.


TITLE: $P\cong P^\ast$ iff $P$ is a f.g projective module?
QUESTION [0 upvotes]: Is it true that for a noncommutative $R$, a module $P$ is f.g projective iff $\mathsf{hom}(P,R)=P^\ast \cong P$?
Here's what I thought of as a proof: Since $(-)^\ast$ is additive, it preserves biproduct diagrams. The projective modules are the direct summands of free modules, and f.g projective modules can be taken as direct summands of f.g free modules. For finitely generated free modules the result is known. Then write $F=P\oplus Q$ for $F$ f.g free and apply the dualizing functor. Since its additive we have $P\oplus Q=F\cong F^\ast =P^\ast \oplus Q^\ast$ and this means (I think) that $P^\ast \cong P$ as desired. So it seems the only difference from the $P\cong P^{\ast\ast}$ is just that this isomorphism is not natural.
Is the above proof correct? What's the justification for inferring $P\cong P^\ast$ from the biproduct diagrams?

REPLY [5 votes]: If $R$ is noncommutative and $P$ is a right $R$-module then $\text{Hom}_R(P, R)$ naturally has the structure of a left $R$-module (and vice versa if $P$ is a left $R$-module), so you can't even ask for this isomorphism because the two objects belong to different categories.
But this is still false if $R$ is commutative. Take $R = \mathcal{O}_K$ to be the ring of integers of a number field $K$ whose ideal class group has a nontrivial element of order anything other than two. The corresponding ideal $I$ is f.g. projective and $I^{\ast}$ corresponds to the inverse of $I$ in the ideal class group; moreover, isomorphism of modules corresponds to equality of elements in the ideal class group.
Your proof fails because you don't ask for any compatibility between the isomorphism $F \cong P \oplus Q$ and the isomorphism $F \cong F^{\ast}$.


TITLE: A function $y$ satisfies the differential equation $y'=x^2y^2$. Find a formula for $y$ in the solution curve passing through $(1,1)$.
QUESTION [0 upvotes]: A function $y$ satisfies the differential equation $y'=x^2y^2$.  Find a formula for $y$ in the solution curve passing through $(1,1)$.

How would I set up an equation?  How should I start?  Thank you.

REPLY [4 votes]: This DE is separable. Another way of writing "passing through point (1,1)" is $y(1)=1$. As the DE is separable, we can integrate
$$\int\frac{dy}{y^2}=\int x^2dx$$
to obtain
$$-\frac{1}{y}=\frac{1}{3}x^3+c$$
where $c$ is the constant of integration. Now apply the condition $y(1)=1$ to obtain the constant.
$$-1=\frac{1}{3}+c\implies c=-\frac{4}{3}$$
So the solution is
$$y=\frac{3}{4-x^3}$$


TITLE: Why is $\exists$ used in the ZFC Axiom of Power Set instead of $\exists !$
QUESTION [0 upvotes]: I recently obtained an old textbook on the theme Mathematical Structures (in my mother language, Slovenian). The content is truly fascinating, so naturally, after I dug through some of the main concepts, I sought universal applicability, or in simpler words, how the notation ultimately corresponds, regardless of the language in which the different textbooks are written.
To my horror, I quickly found a damming difference between my textbook and the ZFC axioms as presented in this Wikipedia website. The Definition of Power Set in my textbook is such:

To put it in LaTeX context:
$$\mathcal{P}A = \{ X; \; X \subseteq A \}$$
Which we could read as: A Power Set of a chosen set A is a set which consists of all the subsets of set A. As you can see, an example is given in the photo I provided as well. Up to this point, everything is normal. But for the existence of a Power Set, we have to propose a special axiom, which ensures a non-contradictory framework (so Cantor sets won't be a problem):

Once again, a LaTeX equivalence:
$$(\forall A)(\exists ! C)(C = \mathcal{P}A)$$
We read: For a given set A, there exists precisely one set C such that set C is the Power Set of set A. I find the little detail in this axiom that is $\exists !$ absolutely necessary, as there mustn't exist multiple modules of one set that derive from one definition. However, the Wikipedia page about ZFC axioms proposes this:
$$\forall X \exists Y \forall Z[Z \subseteq X \Rightarrow Z \in Y]$$
Which, put in words, gives: For every set X, there exists a set Y, such that for every set Z follows that if Z is a subset of X, then Z is an element of Y. At first glance, this axiom is even more elegant than the one found in my textbook, as it uses three sets and doesn't immediately assume their origin. But I am very bothered by the fact that the operand is not a $\exists !$, rather a $\exists$, thus not excluding a possibility of multiple Power Sets, at least it seems so.
Please help me. Though I have a hypothesis as to why this is so, it is NOT straight-forward: My initial thought was that the operand $\exists !$ was not needed in the Wikipedia page, because the definition of a power set was wider, therefore making the specification obsolete. Answers, criticisms and remarks are welcome. Cheers!
EDIT: Please ignore the fact that the author used old mathematical notation (e.g.: the old equivalence of $\exists !$ is $\mathsf{E} !$). Rest assured that all of the symbols mean the same as new notation.

REPLY [4 votes]: There are two reasons:

$\exists!$ is not a standard quantifier. It is a shorthand, $\exists!x\varphi(x)$ is really $\exists x(\varphi(x)\land\forall y(\varphi(y)\rightarrow x=y))$. This means that any statement where $\exists!$ is used, is inherently more complicated than a statement where only $\exists$ is used.
The axiom of extensionality ensures that definitions give us unique sets. In particular, if there are two sets whose elements are exactly the subsets of $x$, then the two sets are equal. Therefore, we can prove that $\exists y$ is really $\exists!y$, when we define $y$ to be the power set of $x$.


text:
\begin{document}

\newcommand{\eqd}{\stackrel{d}{=}}
\renewcommand{\r}{\mathbb R}
\renewcommand{\Pr}{\mathbb P}

\renewcommand{\refname}{References}

\title{\vspace{-1.8cm}
A functional approach to estimation of the parameters of generalized negative binomial and gamma distributions}

\author{
A.\,K.~Gorshenin\textsuperscript{1},
V.\,Yu.~Korolev\textsuperscript{2}}

\date{}

\maketitle

\footnotetext[1]{Institute of Informatics Problems, Federal Research
Center ''Computer Science and Control'' of Russian Academy of
Sciences, Russia; \url{agorshenin@frccsc.ru}}

\footnotetext[2]{Faculty of Computational Mathematics and
Cybernetics, Lomonosov Moscow State University, Russia; Institute of
Informatics Problems, Federal Research Center ``Computer Science and
Control'' of Russian Academy of Sciences, Russia; Hangzhou Dianzi
University, China; \url{vkorolev@cs.msu.su}}

\maketitle

\begin{abstract}
The generalized negative binomial distribution (GNB) is a new
flexible family of discrete distributions that are mixed Poisson
laws with the mixing generalized gamma (GG) distributions. This
family of discrete distributions is very wide and embraces Poisson
distributions, negative binomial distributions, Sichel
distributions, Weibull--Poisson distributions and many other types
of distributions supplying descriptive statistics with many flexible
models. These distributions seem to be very promising for the
statistical description of many real phenomena. GG distributions are widely applied in signal and image processing and other practical problems. The statistical estimation of the parameters of GNB and GG distributions is quite complicated. To find estimates,
the methods of moments or maximum likelihood can be used as well as
two-stage grid EM-algorithms. The paper presents a methodology based
on the search for the best distribution using the minimization of
$\ell^p$-distances and $L^p$-metrics for GNB and GG distributions,
respectively. This approach, first, allows to obtain parameter
estimates without using grid methods and solving systems of
nonlinear equations and, second, yields not point estimates as the
methods of moments or maximum likelihood do, but the estimate for
the density function. In other words, within this approach the set
of decisions is not a Euclidean space, but a functional space.
\end{abstract}

\section{Introduction}
The generalized negative binomial distribution (GNB) is a new
flexible family of discrete distributions that are mixed Poisson
laws with the mixing generalized gamma (GG) distributions. The GNB
distributions were introduced and studied
in~\cite{Korolev2017} under the name of GG mixed Poisson
distributions. This family of discrete distributions is very wide
and embraces Poisson distributions (as limit points corresponding to
a degenerate mixing distribution), negative binomial (Polya)
distributions including geometric distributions (corresponding to
the gamma mixing distribution, see~\cite{Greenwood1920}), Sichel
distributions (corresponding to the inverse gamma mixing
distributions, see~\cite{Sichel1971}), Weibull--Poisson
distributions (corresponding to the Weibull mixing distributions,
see~\cite{Korolev2016}) and many other types supplying
descriptive statistics with many flexible models. These
distributions seem to be very promising for the statistical
description of many real phenomena being very convenient and almost
universal models. It is quite natural to expect that, having
introduced one more free parameter into the pure negative binomial
model, namely, the power parameter in the exponent of the original
gamma mixing distribution, instead of the negative binomial model
one might obtain a more flexible GNB model that provides even better
fit with the statistical data. For example, GNB distributions can be
successfully applied to modeling statistical regularities in
duration of specific periods in data.

The GG distributions are proposed in order to have a flexible
Bayesian model with a mixing (prior) distribution which is
``responsible'' for the description of statistical regularities of
the manifestation of external stochastic factors. The class of
GG distributions was first described as a unitary family in $1962$
by E. Stacy~\cite{Stacy1962}. The family of GG distributions
contains practically all the most popular absolutely continuous
distributions concentrated on the non-negative half-line including
Weibull and gamma distributions.

GG distributions are widely applied in many practical problems.
There are dozens of papers dealing with the application of
GG distributions as models of regularities observed in practice.
As an example, the following research areas involving
models based on GG distributions can be mentioned:
\begin{itemize}
\item climatic and hydrological problems: drop size distributions~\cite{Maur2001},
drought data~\cite{Nadarajah2007},  phenomena in warm clouds~\cite{Xie2009};
\item synthetic-aperture radar (SAR) image processing and various applications:
distribution for the real and imaginary parts of the complex SAR backscattered signal~\cite{Li2010},
flexible model for the SAR images with different land-cover typologies~\cite{Li2011},
statistical modeling of SAR images~\cite{Qin2015,Sportouche2017};
\item astrophysical problems, for example, new galaxy luminosity functions~\cite{Zaninetti2010};
\item speech signal processing: parametric characterization of speech spectra~\cite{Shin2005},
modelling speech samples~\cite{Li2008}, real-time implementations of algorithms~\cite{Song2008}.
\end{itemize}

Apparently, the popularity of GG distributions is due to that most
of them can serve as adequate asymptotic approximations, since all
the representatives of the class of GG distributions listed above
appear as limit laws in various limit theorems of probability theory
in rather simple limit schemes.

The problem of statistical estimation of the parameters
of GNB and GG distributions (for example, the search for maximum
likelihood (ML) estimates) is quite complicated. To find the
estimates of the parameters, the method of moments or ML
method~\cite{Huang2006} for the GG distribution as well as the
two-stage grid EM-algorithm for the GNB, can be used. It should be
noted that the implementations of the methods of moments and ML
method  for GG distribution are difficult computational tasks,
moreover, the efficiency depends on the sample size (ML method is
better for large volumes).

The paper presents a methodology based on finding the best
distribution using minimization of $\ell^1$-, $\ell^2$- and
$\ell^{\infty}$-distances and $L^1$-, $L^2$- and
$L^{\infty}$-metrics for GNB and GG distributions, respectively.
This approach, first, allows to obtain parameter estimates without
using grid methods and solving systems of nonlinear equations and,
second, yields not point estimates as the methods of moments or
maximum likelihood do, but the estimate for the density function. In
other words, within this approach the set of decisions is not a
Euclidean space, but a functional space.

\section{The GNB and GG distributions}

It will be assumed that all the random variables are defined on the
same probability space $(\Omega,\,\mathfrak{F},\,{\mathbb P})$.

A random variable having the gamma distribution with shape parameter
$r>0$ and scale parameter $\mu>0$ will be denoted $G_{r,\mu}$,
\begin{equation}
\label{Gamma}
\Pr (G_{r,\mu}<x)=\int_{0}^{x}g(z;r,\mu)dz,\ \ \text{with}\ \
g(x;r,\mu)=\frac{\mu^r}{\Gamma(r)}x^{r-1}e^{-\mu x},\ x\ge0,
\end{equation}
where $\Gamma(r)$ is Euler's gamma-function,
$\Gamma(r)=\int_{0}^{\infty}x^{r-1}e^{-x}dx$, $r>0$.

A GG distribution is the absolutely continuous distribution defined
by the density
\begin{equation}
\label{GG}
g^*(x;r,\gamma,\mu)=\frac{|\gamma|\mu^r}{\Gamma(r)}x^{\gamma
r-1}e^{-\mu x^{\gamma}},\ \ \ \ x\ge0,
\end{equation}
with $\gamma\in\mathbb{R}$, $\mu>0$, $r>0$. The distribution
function corresponding to the density $g^*(x;r,\gamma,\mu)$ can be
denoted $F^*(x;r,\gamma,\mu)$.

The properties of GG distributions were described
in~\cite{Stacy1962,Korolev2013}. A random variable with the
density $g^*(x;r,\gamma,\mu)$ will be denoted
$\overline{G}_{r,\gamma,\mu}$. It can be easily made sure that
\begin{equation}\label{GG1}
\overline{G}_{r,\gamma,\mu}\eqd G_{r,\mu}^{1/\gamma},
\end{equation}
and hence,
\begin{equation}\label{GG2}
(\overline{G}_{r,\gamma,\mu})^{\gamma}\eqd G_{r,\mu}.
\end{equation}

The symbol $\eqd$ in~\eqref{GG1} and~\eqref{GG2} denotes the coincidence of distributions.

A random variable $N_{r,p}$ is said to have the negative binomial
(NB) distribution with parameters $r>0$ (``shape'') and $p\in(0,1)$
(``success probability''), if
\begin{equation}
\label{NB}
{\mathbb P}(N_{r,p}=k)=\frac{\Gamma(r+k)}{k!\Gamma(r)}\cdot p^r(1-p)^k,\
\ \ \ k=0,1,2,...
\end{equation}

Let $r>0$, $\gamma\in\r$ and $\mu>0$. We say that the random
variable $N_{r,\gamma,\mu}$ has the GNB distribution, if
\begin{equation}
\label{GNB}
{\mathbb P}(N_{r,\gamma,\mu}=k)=\frac{1}{k!}\int_{0}^{\infty}e^{-z}z^kg^*(z;r,\gamma,\mu)dz,\
\ \ \ k=0,1,2...,
\end{equation}
and $g^*(z;r,\gamma,\mu)$ is determined by formula~\eqref{GG}.

\section{A functional approach to estimation of the parameters of GNB distributions}

The problem of statistical estimation of the parameters of GNB
distribution (for example, the search for maximum likelihood
estimates) is extremely complicated. To find estimators, the
two-stage grid EM-algorithm for the GNB distribution
$F(x;r,\gamma,\mu)$ can be used. At the first stage, the main part
of the support of the mixing distribution is determined. That is, a
bounded interval is determined such that the probability of a
GG distributed mixing random variable to fall into this interval is
insignificantly less than one. This interval is covered by a finite
grid containing (possibly, a very large number) $K\in\mathbb{N}$ of
known nodes $\lambda_1,...,\lambda_K$. The GNB distribution under
study is approximated by the finite mixture of Poisson
distributions:
\begin{equation}\label{Poismixt}
F(x;r,\gamma,\mu)(x+0)\approx\sum_{j=0}^{[x]}\frac{1}{j!}\sum_{i=1}^Kp_ie^{-\lambda_i}\lambda_i^j,\ \ \ x\in\mathbb{R}.
\end{equation}

In the mixture on the right-hand side of \eqref{Poismixt}, only the
parameters $p_1,...,p_{K}$ are unknown. At the second stage, it
remains to use some standard method for fitting the GG distribution
to the histogram-type data $(\mu_1, p_1),..., (\mu_K, p_K)$,
obtained at the first stage. For example, the parameters $r$,
$\gamma$ and $\mu$ can be determined as the point minimizing the
corresponding chi-square statistic or some special least squares
problem.

However, with a fixed grid, the two-stage method yields only
approximate estimates of the parameters of GG distributions.
Moreover, the accuracy of the approximation depends on the choice of
the grid. The estimates can be consistent in the traditional sense
only if the grid mesh becomes infinitely small as the sample size
infinitely increases in an appropriate way. Moreover, the conditions
unifying the rate of decrease of the grid mesh with the rate of
increase of the sample size that provide the statistical consistency
of the estimators are very cumbersome and practically unverifiable.

In this section we present an alternative methodology based on
finding the best GNB distribution using minimization of $\ell^1$-,
$\ell^2$- and $\ell^{\infty}$-distances (they correspond to the
spaces of sequences whose series are absolutely convergent, the
space of square-summable sequences and the space of bounded
sequences, respectively). Namely, the histogram of the initial data
should be obtained. The integer rule is used as bining algorithm
(due to that the observations in the sample are integer), so bins
are created for each value. Let $N_{b}$ be the number of histogram
bins (with a uniform width that equals $1$), $\bf h$ be the vector
of bar heights ($h_i\in [0,1]$ for all $i=1,\ldots,N_{b}$). The
value of each component $h_i$ is equal to the ratio of a number of
observations in the bin to a total number of observations, the sum
of the bar areas is $1$. So, the bars of empirical distribution can
be approximated by ones of GNB. For finding estimations of unknown
parameters of generalized negative binomial distributions the
following optimization problems should be solved (the density
$g^*(x;r,\gamma,\mu)$ is determined by~\eqref{GG} and the
probability ${\mathbb P}(N_{r,\gamma,\mu}=k)$ is determined
by~\eqref{GNB}).
\begin{itemize}
\item If the target function is based on $\ell^1$-distance:
\begin{equation}
\label{l1}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \sum\limits_{k=1}^{N_{b}} \left| \frac{1}{k!}\int_{0}^{\infty}e^{-z}z^kg^*(z;r,\gamma,\mu)dz - h_k\right|.
\end{equation}
\item If the target function is based on $\ell^2$-distance:
\begin{equation}
\label{l2}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \sqrt{\sum\limits_{k=1}^{N_{b}} \left( \frac{1}{k!}\int_{0}^{\infty}e^{-z}z^kg^*(z;r,\gamma,\mu)dz - h_k\right)^2}.
\end{equation}
\item If the target function is based on $\ell^{\infty}$-distance:
\begin{equation}
\label{lInf}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \max\limits_{k=\overline{1,N_{b}}} \left| \frac{1}{k!}\int_{0}^{\infty}e^{-z}z^kg^*(z;r,\gamma,\mu)dz - h_k\right|.
\end{equation}
\end{itemize}

Formulas~\eqref{l1}--\eqref{lInf} allow to obtain parameter
estimates without using grid methods. It should be noted that this
methodology can also be used for the classical negative binomial
distribution~\eqref{NB} (the ratio~\eqref{NB} should be used in
formulas~\eqref{l1}--\eqref{lInf} instead of~\eqref{GNB}).

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{NB_GNB_l1}
\end{center}
\caption{\label{FigNB_GNB_l1} Approximation of the initial data distribution by optimization of $\ell^1$-distance.}
\end{figure}

A special~\verb"MATLAB" program is implemented for
finding GNB approximations and plotting figures. The numerical
optimization is based on the simplex search
method~\cite{Lagarias1998}. The functions for estimating the values
of all three unknown parameters of the GNB distribution or two
parameters provided the shape parameter $r$ estimate based on NB
distribution is given are created.

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{NB_GNB_l2}
\end{center}
\caption{\label{FigNB_GNB_l2} Approximation of the initial data distribution by optimization of $\ell^2$-distance.}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{NB_GNB_lInf}
\end{center}
\caption{\label{FigNB_GNB_l1Inf} Approximation of the initial data distribution by optimization of $\ell^{\infty}$-distance.}
\end{figure}

The histogram and approximating graphs of the NB and GNB
distributions as well as errors in the corresponding metrics are
plotted. The examples of results are shown on
Figs.~\ref{FigNB_GNB_l1}--\ref{FigNB_GNB_l1Inf}. They demonstrate a
high quality of the approximation of the histogram of the initial
data by each type of distributions. Table~\ref{Tab_NB_GNB}
represents approximation errors, the parameters are estimated by
each of the metrics. Obviously, the results for GNB distributions
are better (see the bold marked items).

\begin{table}
\caption{Approximation errors for negative binomial  and generalized
negative binomial  distributions, test sample.}\label{Tab_NB_GNB}
\centering

\begin{tabular}{|c|c|c|c|}
\hline
{\bf Distribution}&{\bf Error ($\ell^1$)}&{\bf Error ($\ell^2$)}&{\bf Error ($\ell^{\infty}$)} \\
\hline
NB ($\ell^1$-optimization)&$0{,}047$&$0{,}022$&$0{,}018$\\
\hline
GNB ($\ell^1$-optimization)&$\bf 0{,}043$&$0{,}018$&$0{,}014$\\
\hline
NB ($\ell^2$-optimization)&$0{,}051$&$0{,}0195$&$0{,}015$\\
\hline
GNB ($\ell^2$-optimization)&$0{,}05$&$\bf 0{,}015$&$0{,}0097$\\
\hline
NB ($\ell^{\infty}$-optimization)&$0{,}061$&$0{,}022$&$0{,}012$\\
\hline
GNB ($\ell^{\infty}$-optimization)&$0{,}061$&$0{,}017$&$\bf 0{,}007$\\
\hline
\end{tabular}
\end{table}

\section{Recurrence formulas for GNB distributions}

Using formulas~\eqref{GG} and~\eqref{GNB} we can obtain the following results:
\begin{gather*}
\Pr(N_{r,\gamma,\mu}=k)=\frac{|\gamma|\mu^r}{\Gamma(r)k!}\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma}\, z^{\gamma r+k-1}\,dz= \frac{|\gamma|\mu^r}{\Gamma(r)k!}\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma} \,d\frac{z^{\gamma r+k}}{\gamma r+k}=\\
=\frac{|\gamma|\mu^r}{\Gamma(r)k!}\times\left[\frac{z^{\gamma r+k}}{\gamma r+k}\,e^{-z-\mu z^\gamma} \bigg|_{0}^{\infty}-
\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma} \frac{z^{\gamma r+k}}{\gamma r+k} (-1-\mu\gamma z^{\gamma-1})\,dz\right]=\\
=\frac{|\gamma|\mu^r}{\Gamma(r)k!}\times\left[\frac1{\gamma r+k}\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma} z^{\gamma r+k}\,dz +\frac{\mu\gamma}{\gamma r+k}\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma} z^{\gamma r+\gamma+k-1}\,dz\right]=\\
=\frac{k+1}{\gamma r+k}\,\Pr(N_{r,\gamma,\mu}=k+1)+\frac{\gamma^2\mu^{r+1}}{(\gamma r+k)\Gamma(r)k!}\int\limits_{0}^{\infty}
e^{-z-\mu z^\gamma} z^{\gamma r+k-1+\gamma}\,dz=\\
=\frac{k+1}{\gamma r+k}\,\Pr(N_{r,\gamma,\mu}=k+1)+\frac{|\gamma|\mu}{\gamma r+k}\Pr(N_{r+1,\,\gamma,\mu}=k).
\end{gather*}

So, the recurrence formulas for GNB distributions can be represented as follows:
\begin{gather}
(\gamma r+k)\,\Pr(N_{r,\gamma,\mu}=k)=(k+1)\,\Pr(N_{r,\gamma,\mu}=k+1)+|\gamma|\mu\,\Pr(N_{r+1,\,\gamma,\mu}=k),\notag\\
\intertext{or}
\Pr(N_{r,\gamma,\mu}=k+1)=\frac{\gamma r+k}{k+1}\,\Pr(N_{r,\gamma,\mu}=k)-\frac{|\gamma|\mu}{k+1}\,\Pr(N_{r+1,\,\gamma,\mu}=k).
\label{Reccurence}
\end{gather}

Unfortunately, the representation~\eqref{Reccurence}
does not significantly simplify the computational process, since in
addition to the value $\Pr(N_{r,\gamma,\mu}=k)$ a value
$\Pr(N_{r+1,\gamma,\mu}=k)$ should be known.

\section{A functional approach to the estimation of the parameters of GG distributions}

In this section we present a methodology based on the search for the
best GG distribution using minimization of $L^1$-, $L^2$- and
$L^{\infty}$-metrics (they correspond to the spaces of functions for
which the $p^{th}$ power of the absolute value is Lebesgue
integrable, where functions that agree almost everywhere are
identified).

The histogram of the initial data should be obtained. The
Freedman--Diaconis rule~\cite{Freedman1981} is used as bining
algorithm due to its suitableness for data with heavy-tailed
distributions. It uses a bin width of
\begin{equation}
\label{Freedman}
2 \,\frac{x_{0{.}75}-x_{0{.}25}}{\sqrt[3]{n}},
\end{equation}
where $x_{0{.}25}$, $x_{0{.}75}$ are $0{.}25$- and
$0{.}75$-quantiles, numerator of fraction~\eqref{Freedman}
represents an interquartile range and $n$ is a sample size.

Let $N_{b}$ be the number of histogram bins, $\bf h$ be the vector
of bar heights ($h_i\in [0,1]$ for all $i=1,\ldots,N_{b}$). The
value of each component $h_i$ is equal to the ratio of the number of
observations in the bin and the total number of observations, the
sum of the bar areas is $1$. Let $\bf b$ be the vector of bin edges.
The bars of empirical distribution should be approximated by
GG distribution.

To find the estimates of unknown parameters of GG distributions the
following optimization problems should be solved (the density
$g^*(x;r,\gamma,\mu)$ is determined by~\eqref{GG}).
\begin{itemize}
\item If the target function is based on $L^1$-metric:
\begin{equation}
\label{L1}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \sum\limits_{k=1}^{N_{b}-1} \int\limits_{b_k}^{b_{k+1}}\left|g^*(z;r,\gamma,\mu) - h_k\right|\,dz.
\end{equation}
\item If the target function is based on $L^2$-metric:
\begin{equation}
\label{L2}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \sqrt{\sum\limits_{k=1}^{N_{b}-1} \int\limits_{b_k}^{b_{k+1}}\left(g^*(z;r,\gamma,\mu) - h_k\right)^2}\,dz.
\end{equation}
\item If the target function is based on $L^{\infty}$-metric:
\begin{equation}
\label{LInf}
(r^*,\gamma^*,\mu^*)=\arg\min_{r,\gamma,\mu} \max\limits_{k\in[1,N_b-1]}\int\limits_{b_k}^{b_{k+1}}\left|g^*(z;r,\gamma,\mu) - h_k\right|\,dz.
\end{equation}
\end{itemize}

Formulas~\eqref{L1}--\eqref{LInf} allow to obtain parameter
estimates without using grid methods. It should be noted that this
methodology can also be used for the classical gamma
distribution~\eqref{Gamma}.

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{GG_L1}
\end{center}
\caption{\label{FigGG_L1} Approximation of the initial data distribution by optimization of $L^1$-metric.}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{GG_L2}
\end{center}
\caption{\label{FigGG_L2} Approximation of the initial data distribution by optimization of $L^2$-metric.}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{GG_LInf}
\end{center}
\caption{\label{FigGG_LInf} Approximation of the initial data distribution by optimization of $L^{\infty}$-metric.}
\end{figure}

A special \verb"MATLAB" program is implemented for finding
GG approximations and plotting figures. The numerical optimization
is based on the simplex search method~\cite{Lagarias1998}. The
functions for estimating values of all three unknown parameters of
GG distribution or two parameters provided the shape parameter $r$
is given based on the gamma distribution model are created. The
histogram and approximating probability density functions of gamma
and GG distributions as well as errors in the corresponding metrics
are plotted. The examples of results are shown on
Figures~\ref{FigGG_L1}--\ref{FigGG_LInf}.

They demonstrate a high quality of the approximation of the histogram of the initial data by each type of distributions. Table~\ref{Tab_Gamma_GG} represents approximation
errors, the parameters are estimated by each of the metrics. Obviously, the results for GG distributions are better (see the bold marked items).

\begin{table}
\caption{Approximation errors for gamma and generalized gamma distributions, test sample.}\label{Tab_Gamma_GG}
\centering

\begin{tabular}{|c|c|c|c|}
\hline
{\bf Distribution}&{\bf Error ($L^1$)}&{\bf Error ($L^2$)}&{\bf Error ($L^{\infty}$)} \\
\hline
Gamma &$0{,}212$&$0{,}095$&$0{,}044$\\
\hline
GG &$\bf 0{,}155$&$\bf 0{,}064$&$\bf 0{,}035$\\
\hline
GG, fixed $r$&$0{,}211$&$0{,}095$&$0{,}044$\\
\hline
\end{tabular}
\end{table}

\section{Conclusion}

The classical negative binomial distribution was successfully used
as a model for the number of subsequent wet days in precipitation
problems for the data registered in climatically different points
(see, for
example,~\cite{Gorshenin2017a,Gorshenin2017b,Gorshenin2017c}).
It was demonstrated that the fluctuations of the data with very
high confidence fit the negative binomial distribution. Obviously, a
more flexible GNB model could provide even better fit with the
statistical data. Herewith the GG distribution can be
effectively used to model aggregated data (for example, volumes
accumulated over a period) and can be useful for statistical testing
of hypotheses about their extremality.

Moreover, such types of mixed probability models are
quite adequate for information systems (for example, in
insurance~\cite{Grandell1997,Bening2002}, financial
mathematics~\cite{Gorshenin2013a,Korolev2015},
physics~\cite{Korolev2006,Gorshenin2012,Gorshenin2013b}, data
flows~\cite{Gorshenin2013c} and many other fields).
The developed functional methods for the estimation of
the unknown distribution parameters can be implemented as numerical
procedures in the research support system for stochastic data
processing~\cite{Gorshenin2016,Gorshenin2017d} to analyze events in
various information flows.

\section*{Acknowledgments.}

The research is partially supported by the Russian Foundation for
Basic Research (project~17-07-00851) and the RF Presidential
scholarship program (No.~538.2018.5).

config: arxiv