Q}[K^{4+\varepsilon}(\tilde{Z},\tilde{Z})]<+\infty$ for any probability
distribution $Q$. Now from [21, Theorem 3.1] it follows that
$\tilde{\eta}\to\eta$ almost surely.

Recall that ${\mathcal{G}}_{n}$ is the graph constructed on the pooled samples

$(Z_{1},\ldots,Z_{n})=\left(X^{(1)}_{1},\ldots,X^{(1)}_{{n}_{1}},X^{(2)}_{1},\ldots,X^{(2)}_{{n}_{2}},\ldots,X_{1}^{(M)},\ldots,X^{(M)}_{{n}_{M}}\right),$

with $n_{1},\ldots,n_{M}$ being nonrandom, and

$\hat{\eta}=\frac{\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in{\mathcal{G}}_{n}}K({\Delta}_{i},{\Delta}_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K({\Delta}_{i},{\Delta}_{j})}{\frac{1}{n}\sum_{i=1}^{n}K({\Delta}_{i},{\Delta}_{i})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K({\Delta}_{i},{\Delta}_{j})},$

with $(\Delta_{1},\ldots,\Delta_{n})=(1,\ldots,1,\ldots,M,\ldots,M)$ being the
class labels such that $n_{m}$ of them is $m$, i.e.,
$\sum_{i=1}^{n}I(\Delta_{i}=m)=n_{m}$, for $m=1,\ldots,M$. The goal is to show
that $\hat{\eta}-\tilde{\eta}\overset{a.s.}{\longrightarrow}0$.

By the strong law of large numbers,

$\frac{1}{n}\sum_{i=1}^{n}K(\tilde{\Delta}_{i},\tilde{\Delta}_{i})\overset{a.s.}{\longrightarrow}\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})]=\sum_{i=1}^{M}\pi_{i}K(i,i)=\lim_{\frac{n_{i}}{n}\to\pi_{i}}\frac{1}{n}\sum_{i=1}^{n}K({\Delta}_{i},{\Delta}_{i}).$

By the strong law of large numbers for U-statistics, we have

$\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})\overset{a.s.}{\longrightarrow}\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]=\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)=\lim_{\frac{n_{i}}{n}\to\pi_{i}}\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K({\Delta}_{i},{\Delta}_{j}).$

Since $K$ is characteristic, whenever
$\tilde{\Delta}_{1}\neq\tilde{\Delta}_{2}$,
$K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})+K(\tilde{\Delta}_{2},\tilde{\Delta}_{2})-2K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})=\|K(\tilde{\Delta}_{1},\cdot)-K(\tilde{\Delta}_{2},\cdot)\|_{\mathcal{H}}^{2}>0$.
This implies
$\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})-\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})>0$,
so the denominator of $\hat{\eta}$ has a nonzero deterministic limit. Hence we
only need to show

$\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in{\mathcal{G}}_{n}}K({\Delta}_{i},{\Delta}_{j})-\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{\tilde{d}_{i}}\sum\limits_{j:(\tilde{Z}_{i},\tilde{Z}_{j})\in\tilde{\mathcal{G}}_{n}}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})\overset{a.s.}{\longrightarrow}0.$

Recall that ${\mathcal{G}}_{n}$ is the graph constructed on
$X^{(1)}_{1},\ldots,X^{(1)}_{{n}_{1}},X^{(2)}_{1},\ldots,X^{(2)}_{{n}_{2}},\ldots,X_{1}^{(M)},\ldots,X^{(M)}_{{n}_{M}}$,
and $\tilde{\mathcal{G}}_{n}$ is the graph constructed on
$X^{(1)}_{1},\ldots,X^{(1)}_{\tilde{n}_{1}},X^{(2)}_{1},\ldots,X^{(2)}_{\tilde{n}_{2}},\ldots,X_{1}^{(M)},\ldots,X^{(M)}_{\tilde{n}_{M}}$.
Hence $\tilde{\mathcal{G}}_{n}$ can be obtained from ${\mathcal{G}}_{n}$ be
replacing $\sum_{i=1}^{M}|n_{i}-\tilde{n}_{i}|/2$ points. Note that replacing
one point changes at most $q_{n}+t_{n}$ edges. More specifically, if
$\tilde{\mathcal{G}}_{n}$ can be obtained from $\mathcal{G}_{n}$ by replacing
one point with another point (possibly elsewhere), then

$\\#\\{(j,k):(j,k)\in\left(\mathcal{E}(\mathcal{G}_{n})\setminus\mathcal{E}(\tilde{\mathcal{G}}_{n})\right)\cup\left(\mathcal{E}(\tilde{\mathcal{G}}_{n})\setminus\mathcal{E}(\mathcal{G}_{n})\right){\rm\
or\ }i\in\\{j,k\\}\\}\leq q_{n}+t_{n}.$

Hence:

$\left|\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{\tilde{d}_{i}}\sum\limits_{j:(\tilde{Z}_{i},\tilde{Z}_{j})\in\tilde{\mathcal{G}}_{n}}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})\right|\leq\frac{2(q_{n}+t_{n})\|K\|_{\infty}}{nr_{n}}\leq\frac{C}{n},$

for some $C>0$. Therefore, for general $\tilde{\mathcal{G}}_{n}$, by replacing
$\sum_{i=1}^{M}|n_{i}-\tilde{n}_{i}|/2$ points,

$\left|\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{\tilde{d}_{i}}\sum\limits_{j:(\tilde{Z}_{i},\tilde{Z}_{j})\in\tilde{\mathcal{G}}_{n}}K(\tilde{\Delta}_{i},\tilde{\Delta}_{j})\right|\leq\frac{C\sum_{i=1}^{M}|n_{i}-\tilde{n}_{i}|}{2n}\overset{a.s.}{\longrightarrow}0.$

This completes the proof of the theorem. ∎

###  C.8 Proof of Theorem 3

We will use the following result (proved in Appendix C.15) which is an
extension of Pham et al. [58, Theorem 3.2].

######  Theorem C.1.

Suppose we have real numbers $\\{a_{ij}\\}_{i,j=1}^{n}$,
$\\{b_{ij}\\}_{i,j=1}^{n}$ ($a_{ij}$, $b_{ij}$ depend on $n$) satisfying
$a_{ii}=b_{ii}=0$, $a_{ij}=a_{ji}$, $b_{ij}=b_{ji}$ for all $i,j$, and $\sum
a_{ij}=\sum b_{ij}=0$. Let $(R_{1},\ldots,R_{n})$ be a random permutation of
$(1,\ldots,n)$, with all permutations being equally likely. Set
$B_{ij}:=b_{R_{i}R_{j}}$. Suppose $t_{n}$ is a deterministic sequence such
that $\frac{t_{n}^{r}}{n}\to 0$, as $n\to\infty$, for all $r\in\mathbb{N}$,
and

  * (A1)

$\max_{i}\sum_{j}|a_{ij}|=O(\max|a_{ij}|\cdot t_{n})$,

  * (A2)

$\liminf_{n}\sum a_{ij}^{2}/(n\max a_{ij}^{2})>0$,

  * (B1)

$\limsup_{n}\left(\frac{2\sum b_{i+}^{2}}{n}\right)/\sum b_{ij}^{2}<1$,  where
$b_{i+}:=\sum_{j=1}^{n}b_{ij}$,

  * (B2)

$\sum|b_{ij}|^{r}/n^{2}=O\left((\sum b_{ij}^{2}/n^{2})^{\frac{r}{2}}\right)$,
for $r=3,4,\ldots$.

Then

$\frac{\sum a_{ij}B_{ij}}{w_{n}}\overset{d}{\to}N(0,1)$

where
$w_{n}^{2}:=4\sum^{\prime}a_{ij}a_{ik}\sum^{\prime}b_{ij}b_{ik}/n^{3}+2\sum
a_{ij}^{2}\sum b_{ij}^{2}/n^{2}$, and $\sum^{\prime}$ means summing over all
possible distinct indices (so here it means summing over all distinct triplets
$(i,j,k)$).

######  Remark C.1 (On the normalizing constant $w_{n}$).

Note that $w_{n}^{2}$ is not the variance of $\sum a_{ij}B_{ij}$ but they are
asymptotically equivalent in the sense that $w_{n}^{2}/{\rm Var}(\sum
a_{ij}B_{ij})\to 1$. We write $w_{n}$ here because it arises naturally in the
proof, and it can be replaced by ${\rm Var}(\sum a_{ij}B_{ij})$ in the
statement of Theorem C.1.

We now show how Theorem C.1 can be used to prove Theorem 3. Note that given
$\mathcal{F}_{n}$, the denominator of $\hat{\eta}$ is deterministic, so we
only need to focus on the numerator of $\hat{\eta}$. For $i\neq j$, let

$\displaystyle a_{ij}$
$\displaystyle:=\frac{1}{d_{i}}I\\{(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})\\}+\frac{1}{d_{j}}I\\{(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})\\}-\frac{2}{n-1},$
(C.3) $\displaystyle b_{ij}$
$\displaystyle=K(\Delta_{i},\Delta_{j})-\frac{1}{n(n-1)}\sum_{p\neq
q}K(\Delta_{p},\Delta_{q}).$

Then, $a_{ii}=b_{ii}=0$ for all $i$, and $\sum a_{ij}=\sum b_{ij}=0$. Note
that $\sum a_{ij}B_{ij}$ has the same distribution as the permutation
distribution of (the numerator of) $\hat{\eta}$, up to scaling by a constant.
Recall Assumption 3 on the boundedness of $\frac{t_{n}}{r_{n}}$, where $t_{n}$
and $r_{n}$ are the upper bound and lower bound on vertex degrees in
$\mathcal{G}_{n}$. To verify (A1), note that

$\displaystyle\sum_{j}|a_{ij}|$ $\displaystyle=\sum_{j:(Z_{i},Z_{j})\ {\rm or\
}(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}|a_{ij}|+\sum_{j:(Z_{i},Z_{j}),(Z_{j},Z_{i})\notin\mathcal{E}(\mathcal{G}_{n})}|a_{ij}|$
$\displaystyle\leq 2+\sum_{j:(Z_{i},Z_{j})\ {\rm or\
}(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{2}{r_{n}}$
$\displaystyle\leq 2+\frac{2t_{n}}{r_{n}}\;=\;O(1)$
$\displaystyle=O(\max|a_{ij}|\cdot t_{n}),$

where in the last line we have used the fact that $\max|a_{ij}|\cdot
t_{n}\geq\left(\frac{1}{r_{n}}-\frac{2}{n-1}\right)\cdot t_{n}\geq
1-\frac{2t_{n}}{n-1}\to 1$ as $n\to\infty$. To check (A2), observe that

$\displaystyle\sum a_{ij}^{2}$
$\displaystyle\geq\sum_{i}\sum_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}a_{ij}^{2}$
(C.4)
$\displaystyle\geq\sum_{i}\sum_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\left(\frac{1}{d_{i}}-\frac{2}{n-1}\right)^{2}\;=\;\sum_{i}\left(\frac{1}{d_{i}}-\frac{2}{n-1}\right)^{2}d_{i}$
$\displaystyle\gtrsim\sum_{i}\frac{1}{d_{i}^{2}}\cdot
d_{i}\;\geq\;\frac{n}{t_{n}}\;\geq\;\frac{n}{t_{n}}\left(\frac{\max|a_{ij}|}{\frac{2}{r_{n}}}\right)^{2}$
$\displaystyle\gtrsim\frac{r_{n}}{4}n\,\max a_{ij}^{2}\;\geq\;\frac{1}{4}n\max
a_{ij}^{2},$

where $x\gtrsim y$ means $x\geq Cy$ for some fixed $C>0$. The third line
follows as $\max|a_{ij}|\leq\frac{2}{r_{n}}$, and the second $\gtrsim$ follows
from $\frac{t_{n}}{r_{n}}$ being bounded (by Assumption 3).

To check (B1), we first show that $\frac{1}{n^{2}}\sum
b_{ij}^{2}\overset{a.s.}{\longrightarrow}{\rm
Var}\left[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right]$, where
$\tilde{\Delta}_{i}\overset{i.i.d.}{\sim}{\rm Multinoulli}(\pi)$.

$\displaystyle\frac{1}{n^{2}}\sum b_{ij}^{2}=$ $\displaystyle\sum_{i,j=1,i\neq
j}^{M}\frac{n_{i}n_{j}}{n^{2}}\left(K(i,j)-\frac{1}{n(n-1)}\sum_{p\neq
q}K(\Delta_{p},\Delta_{q})\right)^{2}$
$\displaystyle+\sum_{i=1}^{M}\frac{n_{i}(n_{i}-1)}{n^{2}}\left(K(i,i)-\frac{1}{n(n-1)}\sum_{p\neq
q}K(\Delta_{p},\Delta_{q})\right)^{2}.$

Since $\frac{n_{i}}{n}\to\pi_{i}$ and $\frac{1}{n(n-1)}\sum_{p\neq
q}K(\Delta_{p},\Delta_{q})\overset{a.s.}{\longrightarrow}\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]$,
we have:

$\frac{1}{n^{2}}\sum
b_{ij}^{2}\overset{a.s.}{\longrightarrow}\sum_{i,j=1}^{M}\pi_{i}\pi_{j}\left(K(i,j)-\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]\right)^{2}={\rm
Var}\left[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right].$ (C.5)

Next, we show $\frac{\sum
b_{i+}^{2}}{n^{3}}\overset{a.s.}{\longrightarrow}\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})K(\tilde{\Delta}_{1},\tilde{\Delta}_{3})-\left(\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right)^{2}$.
Since $b_{ii}=0$,

$\displaystyle\frac{\sum b_{i+}^{2}}{n^{3}}$
$\displaystyle=\frac{\sum^{\prime}b_{ij}b_{ik}+\sum b_{ij}^{2}}{n^{3}}.$

Similar to the previous argument, we have

$\displaystyle\frac{\sum^{\prime}b_{ij}b_{ik}}{n^{3}}$
$\displaystyle\overset{a.s.}{\longrightarrow}\sum_{i,j,k=1}^{M}\pi_{i}\pi_{j}\pi_{k}\left(K(i,j)-\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]\right)\left(K(i,k)-\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]\right)$
$\displaystyle=\mathbb{E}\left[\left(K(\tilde{\Delta}_{3},\tilde{\Delta}_{4})-\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]\right)\left(K(\tilde{\Delta}_{3},\tilde{\Delta}_{5})-\mathbb{E}[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})]\right)\right]$
$\displaystyle=\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})K(\tilde{\Delta}_{1},\tilde{\Delta}_{3})-\left(\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right)^{2}.$

Together with $\frac{1}{n^{2}}\sum
b_{ij}^{2}\overset{a.s.}{\longrightarrow}{\rm
Var}\left[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right]$ which implies
$\frac{1}{n^{3}}\sum b_{ij}^{2}\overset{a.s.}{\longrightarrow}0$, we have:

$\frac{\sum
b_{i+}^{2}}{n^{3}}\overset{a.s.}{\longrightarrow}\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})K(\tilde{\Delta}_{1},\tilde{\Delta}_{3})-\left(\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right)^{2}.$
(C.6)

To verify (B1), it remains to check that

$2\cdot\frac{\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})K(\tilde{\Delta}_{1},\tilde{\Delta}_{3})-\left(\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right)^{2}}{{\rm
Var}\left[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right]}<1.$ (C.7)

This follows from the fact

$\mathbb{E}\left(K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})-\mathbb{E}_{3}K(\tilde{\Delta}_{1},\tilde{\Delta}_{3})-\mathbb{E}_{3}K(\tilde{\Delta}_{3},\tilde{\Delta}_{2})+\mathbb{E}_{3,4}K(\tilde{\Delta}_{3},\tilde{\Delta}_{4})\right)^{2}>0.$

The inequality is strict because otherwise
$K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})=f(\tilde{\Delta}_{1})+g(\tilde{\Delta}_{2})$
almost surely, which implies $K(p,q)=f(p)+g(q)$ for all
$p,q\in\\{1,\ldots,M\\}$, and consequently
$\sum_{p,q=1}^{M}K(p,q)x_{p}x_{q}=0$ as long as $x_{1}+\ldots+x_{M}=0$,
contradicting the fact that $K(\cdot,\cdot)$ is characteristic.

To verify (B2), for any $r=3,4,\ldots$, since the kernel is bounded,
$|b_{ij}|^{r}$ is also bounded. Together with the fact that
$\frac{1}{n^{2}}\sum b_{ij}^{2}\overset{a.s.}{\longrightarrow}{\rm
Var}\left[K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})\right]>0$, we have

$\sum|b_{ij}|^{r}/n^{2}=O(1)=O\left((\sum
b_{ij}^{2}/n^{2})^{\frac{r}{2}}\right).$

Hence all the conditions needed for Theorem C.1 are satisfied and we have the
asymptotic normality of the numerator of $\hat{\eta}$.

Next, we simplify the permutation variance ${\rm
Var}(\hat{\eta}|\mathcal{F}_{n})$. Since if the pooled sample is given, the
graph $\mathcal{G}_{n}$ is known up to a permutation. We can arbitrarily fix a
labeling of the vertices as $1,\ldots,n$ with out-degree $d_{i}$ (which is a
slight abuse of notation, since $d_{i}$ was originally defined as the degree
of $Z_{i}$ whose label is $\Delta_{i}$). Let
$(\check{\Delta}_{1},\ldots,\check{\Delta}_{n})$ be a uniformly random
permutation of $({\Delta}_{1},\ldots,{\Delta}_{n})$. Since the denominator of
$\hat{\eta}$ is constant, it suffices to compute the variance of

$N_{n}=\sqrt{n}\left(\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\check{\Delta}_{i},\check{\Delta}_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\right).$

Since $\mathcal{G}_{n}$ does not contain any self-loop, for any $i\neq
j\in\\{1,\ldots,n\\}$,
$\mathbb{E}[K(\check{\Delta}_{i},\check{\Delta}_{j})]=\frac{1}{n(n-1)}\sum\limits_{l\neq
m}K(\Delta_{l},\Delta_{m})$. This implies $\mathbb{E}[N_{n}]=0$ and

${\rm
Var}(N_{n})=\frac{1}{n}\mathbb{E}\left[\left(\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\check{\Delta}_{i},\check{\Delta}_{j})\right)^{2}\right]-\frac{1}{n(n-1)^{2}}\left(\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\right)^{2}.$

Let

$\displaystyle\tilde{a}$
$\displaystyle:=\mathbb{E}[K^{2}(\check{\Delta}_{1},\check{\Delta}_{2})]=\frac{1}{n(n-1)}{\sum}^{\prime}K^{2}(\Delta_{i},\Delta_{j}),$
$\displaystyle\tilde{b}$
$\displaystyle:=\mathbb{E}[K(\check{\Delta}_{1},\check{\Delta}_{2})K(\check{\Delta}_{1},\check{\Delta}_{3})]=\frac{1}{n(n-1)(n-2)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{i},\Delta_{l}),$
$\displaystyle\tilde{c}$
$\displaystyle:=\mathbb{E}[K(\check{\Delta}_{1},\check{\Delta}_{2})K(\check{\Delta}_{3},\check{\Delta}_{4})]=\frac{1}{n(n-1)(n-2)(n-3)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{l},\Delta_{m}),$

where $\sum^{\prime}$ means the summation indices are required to be distinct.
Then

$\displaystyle\quad\mathbb{E}\left[\left(\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\check{\Delta}_{i},\check{\Delta}_{j})\right)^{2}\right]$
$\displaystyle=\tilde{a}\left(\sum_{i=1}^{n}\frac{d_{i}}{d_{i}^{2}}+\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{j}}\right)$
$\displaystyle+\tilde{b}\left(\sum_{i=1}^{n}\frac{d_{i}(d_{i}-1)}{d_{i}^{2}}+\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{s},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}+\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{s}),(Z_{s},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}+\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{s},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}\right)$
$\displaystyle+\tilde{c}\cdot{\sum}^{\prime}\frac{\left(d_{i}-1_{(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\right)\left(d_{j}-1_{(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\right)-T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}},$

where
$T^{\mathcal{G}_{n}}(Z_{i},Z_{j})=\sum_{k=1}^{n}I\\{(Z_{i},Z_{k}),(Z_{j},Z_{k})\in\mathcal{E}(\mathcal{G}_{n})\\}$
is the number of common out-neighbors of $Z_{i}$ and $Z_{j}$. Note that

$\displaystyle\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{s},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}$
$\displaystyle=\operatorname*{{\sum}^{\prime}}_{(Z_{s},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{d_{i}-1_{(Z_{i},Z_{s})\in\mathcal{E}(\mathcal{G}_{n})}}{d_{i}d_{s}},$
$\displaystyle\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{s}),(Z_{s},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}$
$\displaystyle=\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{s})\in\mathcal{E}(\mathcal{G}_{n})}\frac{d_{s}-1_{(Z_{s},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}}{d_{i}d_{s}},$
$\displaystyle\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{s},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{s}}$
$\displaystyle={\sum}^{\prime}\frac{T^{\mathcal{G}_{n}}(Z_{i},Z_{s})}{d_{i}d_{s}}.$

Recall that $\tilde{g}_{1}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$,
$\tilde{g}_{2}=\frac{1}{n}\sum_{i,j=1}^{n}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}=\frac{1}{n}\sum^{\prime}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}+\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$,
and
$\tilde{g}_{3}=\frac{1}{n}\operatorname*{{\sum}^{\prime}}_{(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{j}}$.
Hence we have:

$\frac{1}{n}\mathbb{E}\Bigg{[}\Bigg{(}\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\check{\Delta}_{i},\check{\Delta}_{j})\Bigg{)}^{2}\Bigg{]}=\tilde{a}(\tilde{g}_{1}+\tilde{g}_{3})+\tilde{b}(3-2\tilde{g}_{1}+\tilde{g}_{2}-2\tilde{g}_{3})+\tilde{c}(n-3+\tilde{g}_{1}-\tilde{g}_{2}+\tilde{g}_{3}).$

The second term of ${\rm Var}(N_{n})$ is easier to handle:

$\left(\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\right)^{2}=n(n-1)(n-2)(n-3)\cdot\tilde{c}+4n(n-1)(n-2)\cdot\tilde{b}+2n(n-1)\cdot\tilde{a}.$

Combining the two, we get

${\rm
Var}(N_{n})=\tilde{a}\Bigg{(}\tilde{g}_{1}+\tilde{g}_{3}-\frac{2}{n-1}\Bigg{)}+\tilde{b}\Bigg{(}\tilde{g}_{2}-2\tilde{g}_{1}-2\tilde{g}_{3}-1+\frac{4}{n-1}\Bigg{)}+\tilde{c}\Bigg{(}\tilde{g}_{1}-\tilde{g}_{2}+\tilde{g}_{3}+\frac{n-3}{n-1}\Bigg{)}.$
(C.8)

Hence

$\frac{\hat{\eta}}{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}=\frac{\sqrt{n}\left(\frac{1}{n}\sum\limits_{i=1}^{n}\frac{1}{{d}_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K({\Delta}_{i},{\Delta}_{j})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\right)}{\tilde{a}\left(\tilde{g}_{1}+\tilde{g}_{3}-\frac{2}{n-1}\right)+\tilde{b}\left(\tilde{g}_{2}-2\tilde{g}_{1}-2\tilde{g}_{3}-1+\frac{4}{n-1}\right)+\tilde{c}\left(\tilde{g}_{1}-\tilde{g}_{2}+\tilde{g}_{3}+\frac{n-3}{n-1}\right)}.$

From the conditional CLT, the unconditional CLT also follows:

$\mathbb{P}\left(\frac{\hat{\eta}}{{\rm Var}(\hat{\eta}|\mathcal{F}_{n})}\leq
x\right)=\mathbb{E}\left[\mathbb{P}\left(\frac{\hat{\eta}}{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}\leq
x|\mathcal{F}_{n}\right)\right]\to\Phi(x),$

by the dominated convergence, where $\Phi(x)$ is the cumulative distribution
function of $N(0,1)$. ∎

###  C.9 Proof of Corollary 2

Note that $\tilde{a},\tilde{b},\tilde{c}$ are bounded as the kernel is
bounded. Further, $\tilde{g}_{1}\leq\frac{1}{r_{n}}\leq 1$,

$\tilde{g}_{2}=\frac{1}{n}\operatorname*{{\sum}^{\prime}}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}+\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}\leq\frac{t_{n}^{2}}{r_{n}^{2}}+\frac{1}{r_{n}}=O(1)$

since $t_{n}$ is an upper bound for the in-degree, and

$\tilde{g}_{3}\leq\frac{1}{n}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{j}}\leq\frac{2t_{n}}{r_{n}^{2}}=O(1).$

Hence ${\rm Var}(N_{n})=O(1)$ and therefore ${\rm
Var}(\hat{\eta}|\mathcal{F}_{n})=O\left(\frac{1}{n}\right)$. If the $M$
distributions are not equal, then $\eta>0$, and
$\hat{\eta}\overset{a.s.}{\longrightarrow}\eta$ since Assumption 1 holds.
Hence:

$\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}\overset{a.s.}{\longrightarrow}+\infty,$

which implies the universal consistency of the test based on $\hat{\eta}$.

###  C.10 Proof of Theorem 4

The convergence of the denominator of $\hat{\eta}$ (in (4)) is easy to see, as
$n\to\infty$,

$\frac{1}{n}\sum_{i=1}^{n}K(\Delta_{i},\Delta_{i})-\frac{1}{n(n-1)}\sum\limits_{i\neq
j}K(\Delta_{i},\Delta_{j})\overset{a.s.}{\longrightarrow}\sum_{i=1}^{M}\pi_{i}K(i,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)>0.$

To see why the limit is positive, we write
$\tilde{\Delta}_{i}\overset{i.i.d.}{\sim}{\rm Multinoulli}(\pi)$, for
$i=1,\ldots,n$, and it is equivalent to showing that
$\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})-\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})>0$.
Since $K(\cdot,\cdot)$ is characteristic, the distance between Dirac measures
${\rm MMD}(\delta_{i},\delta_{j})>0$ whenever $i\neq j\in\\{1,\ldots,M\\}$, so
equivalently
$K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})+K(\tilde{\Delta}_{2},\tilde{\Delta}_{2})-2K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})>0$
whenever $\tilde{\Delta}_{1}\neq\tilde{\Delta}_{2}$. After taking the
expectation, we have
$\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{1})-\mathbb{E}K(\tilde{\Delta}_{1},\tilde{\Delta}_{2})>0$.

The (conditional) variance of $N_{n}$ (in Equation C.8), i.e., the numerator
of $\hat{\eta}$, involves $\tilde{a},\tilde{b},\tilde{c}$ and
$\tilde{g}_{1},\tilde{g}_{2},\tilde{g}_{3}$. It is clear that:

$\displaystyle\tilde{a}$
$\displaystyle=\frac{1}{n(n-1)}{\sum}^{\prime}K^{2}(\Delta_{i},\Delta_{j})\overset{a.s.}{\longrightarrow}\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K^{2}(i,j)=a,$
$\displaystyle\tilde{b}$
$\displaystyle=\frac{1}{n(n-1)(n-2)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{i},\Delta_{l})\overset{a.s.}{\longrightarrow}\sum_{i,j,l=1}^{M}\pi_{i}\pi_{j}\pi_{l}K(i,j)K(i,l)=b,$
$\displaystyle\tilde{c}$
$\displaystyle=\frac{1}{n(n-1)(n-2)(n-3)}{\sum}^{\prime}K(\Delta_{i},\Delta_{j})K(\Delta_{l},\Delta_{m})$
$\displaystyle\overset{a.s.}{\longrightarrow}\sum_{i,j,l,m=1}^{M}\pi_{i}\pi_{j}\pi_{l}\pi_{m}K(i,j)K(l,m)=\left(\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)\right)^{2}=c,\
{\rm as\ }n\to\infty.$

To show the convergence of $\tilde{g}_{1},\tilde{g}_{2},\tilde{g}_{3}$, note
that they are functions of $\mathcal{G}(Z_{1},\ldots,Z_{n})$, with
$\\{Z_{1},\ldots,Z_{n}\\}$ being the pooled sample, and $\mathcal{G}$ being a
stabilizing graph as defined in Section 4.2. Further, they can all be
expressed in terms of $\sum_{i=1}^{n}\xi(Z_{i},\\{Z_{1},\ldots,Z_{n}\\})$, and
thus can be analyzed using the tool of stabilizing functions to be introduced
next.

######  Definition C.1 (Stabilizing functions [56]).

Let $\xi(x,\mathcal{X})$ be a measurable $\mathbb{R}^{+}$-valued function
defined for all finite set $\mathcal{X}\subset\mathbb{R}^{d}$ and
$x\in\mathcal{X}$. If $x\notin\mathcal{X}$, define
$\xi(x,\mathcal{X}):=\xi(x,\mathcal{X}\cup\\{x\\})$. $\xi$ is said to be
translation invariant, if $\xi(x+y,\mathcal{X}+y)=\xi(x,\mathcal{X})$ for all
finite set $\mathcal{X}\subset\mathbb{R}^{d}$ and $x,y\in\mathbb{R}^{d}$. For
a locally finite131313$\mathcal{X}\subset\mathbb{R}^{d}$ is said to be locally
finite if its intersection with any compact set is finite. set $\mathcal{X}$,
if

$\displaystyle\liminf_{m\to\infty}\underbrace{\inf_{n\in\mathbb{N}}\left(\mathop{\rm
ess\ inf}_{\begin{subarray}{c}\mathcal{A}\subset\mathbb{R}^{d}\backslash
B(\mathbf{0},m)\\\
|\mathcal{A}|=n\end{subarray}}\xi\left(\mathbf{0},\big{(}\mathcal{X}\cap
B(\mathbf{0},m)\big{)}\cup\mathcal{A}\right)\right)}_{{\rm denoted\ by\
}\underline{\xi}(\mathcal{X},m)}$ $\displaystyle=$
$\displaystyle\limsup_{m\to\infty}\underbrace{\sup_{n\in\mathbb{N}}\left(\mathop{\rm
ess\ sup}_{\begin{subarray}{c}\mathcal{A}\subset\mathbb{R}^{d}\backslash
B(\mathbf{0},m)\\\
|\mathcal{A}|=n\end{subarray}}\xi\left(\mathbf{0},\big{(}\mathcal{X}\cap
B(\mathbf{0},m)\big{)}\cup\mathcal{A}\right)\right)}_{{\rm denoted\ by\
}\bar{\xi}(\mathcal{X},m)},$

where the essential supremum/infimum is taken with respect to the Lebesgue
measure on $\mathbb{R}^{dn}$, then $\xi$ is said to stabilize on
$\mathcal{X}$. In such a case, we define $\xi(\mathbf{0},\mathcal{X})$ as the
limit of the above quantity.

We will be interested in functions that stabilize on the homogeneous Poisson
process $\mathcal{P}_{\lambda}$, as the pooled sample
$\\{Z_{1},\ldots,Z_{n}\\}$ is locally close to a homogeneous Poisson process.
These functions arise naturally from stabilizing graphs. For example, if
$\mathcal{G}$ is stabilizing on $\mathcal{P}_{\lambda}$, and $\xi$ satisfies,
for any finite set $\mathcal{X}$ containing $\mathbf{0}$,
$\xi(\mathbf{0},\mathcal{X})$ only depends on the edge set
$\mathcal{E}\left(\mathbf{0},\mathcal{G}(\mathcal{X})\right)$, then $\xi$
stabilizes on $\mathcal{P}_{\lambda}$.

######  Theorem C.2 ([56, Theorem 2.1]).

Suppose $q=1$ or $q=2$. Let $X_{1},X_{2},\ldots$ be i.i.d. $d$-dimensional
random variables with common density $f$ and
$\mathcal{X}_{n}:=\\{X_{1},\ldots,X_{n}\\}$. Suppose $\xi$ is translation
invariant and almost surely stabilizing on the homogeneous Poisson process
$\mathcal{P}_{\lambda}$ for all $\lambda>0$. If $\xi$ satisfies the moment
condition

$\sup_{n\in\mathbb{N}}\mathbb{E}\left[\xi(n^{1/d}X_{1},n^{1/d}\mathcal{X}_{n})^{p}\right]<\infty$

for some $p>q$, then as $n\to\infty$,

$\frac{1}{n}\sum_{x\in\mathcal{X}_{n}}\xi(n^{1/d}x,n^{1/d}\mathcal{X}_{n})\to\int_{\mathbb{R}^{d}}\mathbb{E}[\xi(\mathbf{0},\mathcal{P}_{f(x)})]f(x){\rm
d}x,\quad in\ L^{q}$

and the right-hand side above is finite. Here $\mathcal{P}_{f(x)}$ is a
homogeneous Poisson process with intensity $f(x)$.

In particular, if $\xi(x,\mathcal{X})$ is also scale invariant (i.e.,
$\xi(ax,a\mathcal{X})=\xi(x,\mathcal{X})$ for all $a>0$), then

$\frac{1}{n}\sum_{x\in\mathcal{X}_{n}}\xi(x,\mathcal{X}_{n})\to\mathbb{E}[\xi(\mathbf{0},\mathcal{P}_{1})],\quad
in\ L^{q}.$

Let us now get back to the proof of our Theorem 4. Since the graph
$\mathcal{G}$ is stabilizing,
$\xi(x,\mathcal{X}):=\frac{1}{d(x,\mathcal{G}(\mathcal{X}))}$ is a bounded
stabilizing function. Using the above theorem, by the same coupling used in
Appendix C.7 (the proof of the consistency theorem for $\hat{\eta}$), the
limit of $\tilde{g}_{1}$ is unchanged when we replace $Z_{1},\ldots,Z_{n}$ by
i.i.d. observations from $\sum_{i=1}^{M}\pi_{i}P_{i}$:

$\tilde{g}_{1}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\frac{1}{d\big{(}\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\big{)}}\right]=g_{1}.$

By Penrose et al. [56, Lemma 3.3], $\xi(x,\mathcal{X}):=\sum_{y\neq
z:(y,x),(z,x)\in\mathcal{E}(\mathcal{G}(\mathcal{X}))}\frac{1}{d\left(y,\mathcal{G}(\mathcal{X})\right)d\left(z,\mathcal{G}(\mathcal{X})\right)}$
is also bounded and stabilizing on $\mathcal{P}_{\lambda}$ for any
$\lambda>0$. Hence

$\displaystyle\tilde{g}_{2}-\tilde{g}_{1}$
$\displaystyle=\frac{1}{n}\sum_{i\neq
j}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}=\frac{1}{n}\sum_{i=1}^{n}\sum_{j\neq
k:(Z_{j},Z_{i}),(Z_{k},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}d_{k}}$
$\displaystyle\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{y\neq
z:(y,\mathbf{0}),(z,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(z,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]=g_{2}-g_{1}.$

Similarly, by considering
$\xi(x,\mathcal{X}):=\sum_{y:(y,x),(x,y)\in\mathcal{E}(\mathcal{G}(\mathcal{X}))}\frac{1}{d\left(x,\mathcal{G}(\mathcal{X})\right)d\left(y,\mathcal{G}(\mathcal{X})\right)}$,
we have

$\tilde{g}_{3}=\frac{1}{n}{\sum_{(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}}\frac{1}{d_{i}d_{j}}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{y:(y,\mathbf{0}),(\mathbf{0},y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]=g_{3}.$

Hence the (conditional) variance of the numerator of $\hat{\eta}$, scaled by
$\sqrt{n}$, converges in $L^{2}$ to

$\displaystyle\quad{a}\left({g}_{1}+{g}_{3}\right)+{b}\left({g}_{2}-2{g}_{1}-2{g}_{3}-1\right)+{c}\left({g}_{1}-{g}_{2}+{g}_{3}+1\right)$
$\displaystyle=(a-2b+c)(g_{1}+g_{3})+(g_{2}-1)(b-c).$

To show the limiting variance is positive, note that from (C.7), $a-2b+c>0$,
and from (C.6), $b\geq c$. Hence it remains to show $g_{2}\geq 1$. Note that

$\displaystyle\tilde{g}_{2}$ $\displaystyle=\frac{1}{n}\sum_{i\neq
j}\frac{T^{\mathcal{G}_{n}}(i,j)}{d_{i}d_{j}}+\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\sum_{j\neq
l:(Z_{j},Z_{i}),(Z_{l},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}d_{l}}+\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{j:(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}}\right)^{2}-\frac{1}{n}\sum_{i=1}^{n}\sum_{j:(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}^{2}}+\frac{1}{n}\sum_{i=1}^{n}\frac{1}{d_{i}}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\sum_{j:(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}}\right)^{2}$
$\displaystyle\geq\frac{1}{n^{2}}\left(\sum_{i=1}^{n}\sum_{j:(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{j}}\right)^{2}=1,$

where in the last step we have used the Cauchy-Schwarz inequality. Hence
$g_{2}$, as the $L^{2}$ limit of $\tilde{g}_{2}$, is greater than or equal to
1. ∎

###  C.11 Proof of Proposition B.1

Recall $\phi=\sum_{i=1}^{M}\pi_{i}f_{i}$. Let $\mathcal{O}$ be the interior of
${\rm supp}(\phi)$ — a set with probability 1, under $\phi$. Let
$\mathcal{C}_{i}$, $1\leq i\leq I$, be a finite collection of infinite open
cones covering $\mathbb{R}^{d}$ with angular radius $\pi/12$ and apex at
$\mathbf{0}$. For $x\in\mathcal{O}$ and $1\leq i\leq I$, let
$\mathcal{C}_{i}(x)$ be the translate of $\mathcal{C}_{i}$ with apex at $x$.
Let $\mathcal{C}_{i}^{+}(x)$ be the open cone concentric to $\mathcal{C}_{i}$
with apex $x$ and angular radius $\pi/6$. For $1\leq i\leq I$, let
$N^{-1/d}R_{i}(x)$ be the distance from $x$ to its $k$-th nearest neighbor in
$\mathcal{P}_{N\phi_{N}}\cap\mathcal{C}_{i}^{+}(x)$, if this $k$-th nearest
neighbor exists at a distance less than ${\rm
diam}(\mathcal{C}_{i}(x)\cap\mathcal{O})$, and otherwise set
$N^{-1/d}R_{i}(x)={\rm diam}(\mathcal{C}_{i}(x)\cap\mathcal{O})$. Then
$\max_{1\leq i\leq I}R_{i}(x)$ is a radius of stabilization. Note that
$\mathbb{P}\left[R_{i}(x)>t\right]$ is 0 unless ${\rm
diam}(\mathcal{C}_{i}(x)\cap\mathcal{O})>N^{-1/d}t$. Now, consider ${\rm
diam}(\mathcal{C}_{i}(x)\cap\mathcal{O})\geq N^{-1/d}t$. Then there exists a
$y\in\mathcal{C}_{i}(x)\cap\mathcal{O}$ such that $|y-x|=N^{-1/d}t$. By
convexity, $\frac{x+y}{2}\in\mathcal{C}_{i}(x)\cap\mathcal{O}$. With
$\gamma:=\frac{1}{2}\sin(\pi/12)$, we have $B(\frac{x+y}{2},\gamma
N^{-1/d}t)\subset\mathcal{C}_{i}^{+}(x)\cap\mathcal{O}$. $R_{i}(x)>t$ implies
that $B(\frac{x+y}{2},\gamma N^{-1/d}t)$ contains less than $k$ points. Since
$\phi_{N}$ is bounded below, the number of points in $B(\frac{x+y}{2},\gamma
N^{-1/d}t)$ is a Poisson random variable with mean at least
$N\int_{B(\frac{x+y}{2},\gamma N^{-1/d}t)}\inf\\{\phi_{N}\\}{\rm d}z=\delta
t^{d}$ for some $\delta>0$ independent of $i,x,t,N$. Hence

$\displaystyle\mathbb{P}[R_{i}>t]$ $\displaystyle\leq\mathbb{P}\left[{\rm
Poisson}\left(\delta t^{d}\right)<k\right]$
$\displaystyle=\sum_{i=0}^{k-1}\frac{(\delta t^{d})^{i}}{i!}e^{-\delta
t^{d}}\leq C_{1}e^{-C_{2}t^{d}},$

where $C_{1},C_{2}>0$ are independent of $i,x,t,N$. Hence,

$\displaystyle\tau(t)$ $\displaystyle=\sup_{N\geq
1,x\in\mathcal{O}}\mathbb{P}\left[\max_{1\leq i\leq I}R_{i}>t\right]\leq
IC_{1}e^{-C_{2}t^{d}}.$

which yields the desired result. ∎

###  C.12 Proof of Theorem B.1

We first provide a roadmap for the proof. Recall the definitions of $H_{n}$,
$\mathcal{F}_{n}$, $\\{\kappa_{i}\\}_{i=1}^{2}$ in Appendix B.1. In the
following, we first provide some lemmas that describe the local behavior of
$\mathcal{P}_{N\phi_{N}}$. Next, in Step 1 of the main proof, we show
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\mid\mathcal{F}_{n}\overset{d}{\to}N(0,\kappa_{1}^{2})$,
where the provided lemmas are crucial in establishing the convergence of the
variance, and the CLT is established using Stein’s method. Next, in Step 2, we
show that
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{2}^{2})$.
This is the step where we use $\sqrt{N}\left(\frac{N_{p}}{N}-\pi_{p}\right)\to
0$ and the power-law stabilization. To establish the convergence of variance,
we need power-law stabilization of order $q>d$ (recall that
$\mathcal{Z}\subset\mathbb{R}^{d}$), and to establish the CLT using Stein’s
method, power-law stabilization of order $q>16d$ is needed. Finally, combining
Step 1 and Step 2, we can show that
$H_{n}-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{1}^{2}+\kappa_{2}^{2})$.

####  C.12.1 Preliminaries

The following four lemmas are modifications of Lemma 3.1, Lemma 3.2, and
Proposition 3.1 from Penrose et al. [56], which describe that the local
behavior of $\mathcal{P}_{N\phi_{N}}$ at $z$ is similar to a homogeneous
Poisson process $\mathcal{P}_{\phi(z)}$, where
$\phi_{N}(z):=\sum_{i=1}^{M}\frac{N_{i}}{N}f_{i}(z)$, for
$z\in\mathcal{Z}\subset\mathbb{R}^{d}$, with $\frac{N_{i}}{N}\to\pi_{i}$ as
$N\to\infty$, and $\phi(z):=\sum_{i=1}^{M}\pi_{i}f_{i}(z)$.

We say that $x\in\mathbb{R}^{d}$ is a Lebesgue point of a function
$f:\mathbb{R}^{d}\to\mathbb{R}$, if

$\lim_{r\to 0^{+}}\frac{1}{|B(x,r)|}\int_{B(x,r)}|f(y)-f(x)|{\rm d}y=0.$

If $x$ is a continuity point of $f$, then $x$ is a Lebesgue point of $f$. More
generally, the Lebesgue differentiation theorem states that given $f\in
L^{1}(\mathbb{R}^{d})$, a.e. $x$ is a Lebesgue point of $f$.

######  Lemma C.1.

Recall that $\phi(z)=\sum_{i=1}^{M}\pi_{i}f_{i}(z)$, for
$z\in\mathcal{Z}\subset\mathbb{R}^{d}$. Suppose $x_{0}$ is a Lebesgue point of
$f_{1},\ldots,f_{M}$. Then there exists a homogeneous Poisson process
$\mathcal{P}_{\phi(x_{0})}$, coupled with $\mathcal{P}_{N\phi_{N}}$, such that
for all $K>0$,

$\mathbb{P}\left[N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\cap
B(\mathbf{0},K)=\mathcal{P}_{\phi(x_{0})}\cap B(\mathbf{0},K)\right]\to
1,\quad{\rm as\ }N\to\infty.$

###### Proof.

Let $\mathcal{P}$ be a homogeneous Poisson process of rate 1 on
$\mathbb{R}^{d}\times[0,\infty)$. Let $\mathcal{P}_{N\phi_{N}}$ be the
projection of the set

$\left\\{(x,t)\in\mathcal{P}:t\leq N\phi_{N}(x)\right\\}$

onto the $x$ space, i.e., $(x,t)\mapsto x$. Then $\mathcal{P}_{N\phi_{N}}$ is
a non-homogeneous Poisson process with intensity function $N\phi_{N}$. Let
$\mathcal{P}_{\phi(x_{0})}$ be the image of the point set

$\\{(x,t)\in\mathcal{P}:t\leq N\phi(x_{0})\\}$

under the mapping

$(x,t)\mapsto N^{1/d}(x-x_{0}).$

Note that $\mathcal{P}_{\phi(x_{0})}$ is a homogeneous Poisson process on
$\mathbb{R}^{d}$ with intensity $\phi(x_{0})$ for all $N$. The number of
points in
$\left(N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\triangle\mathcal{P}_{\phi(x_{0})}\right)\cap
B(\mathbf{0},K)$ equals the number of points $(x,t)\in\mathcal{P}$ such that
$x\in B(x_{0},N^{-1/d}K)$ and $t$ is between $\phi_{N}(x)$ and $\phi(x_{0})$,
and follows a Poisson distribution with mean

$N\int_{B(x_{0},N^{-1/d}K)}|\phi_{N}(x)-\phi(x_{0})|{\rm d}x\leq
N\int_{B(x_{0},N^{-1/d}K)}(|\phi_{N}(x)-\phi(x)|+|\phi(x)-\phi(x_{0})|){\rm
d}x,$ (C.9)

which converges to 0 as $x_{0}$ is also a Lebesgue point of $\phi$. ∎

######  Lemma C.2.

Suppose $x_{0}$ is a Lebesgue point of $f_{1},\ldots,f_{M}$ and
$\phi(x_{0})>0$. Let $\xi(x,\mathcal{X})$ be a translation invariant function
that almost surely stabilizes on $\mathcal{P}_{\phi(x_{0})}$ (see Definition
C.1). Suppose the moment condition

$\sup_{N\geq
1}\mathbb{E}\left[\left|\xi\left(\mathbf{0},N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\right)\right|^{p}\right]<\infty$

is satisfied for some $p>1$. Then

$\xi\left(\mathbf{0},N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\right)\to\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)$

in distribution and in expectation.

###### Proof.

Consider the coupling of $\mathcal{P}_{N\phi_{N}}$ and
$\mathcal{P}_{\phi(x_{0})}$ in the previous lemma. We have

$\displaystyle\quad\mathbb{P}\left[\left|\xi\left(\mathbf{0},N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\right)-\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)\right|>\varepsilon\right]$
$\displaystyle\leq\mathbb{P}\big{[}N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\cap
B(\mathbf{0},K)\neq\mathcal{P}_{\phi(x_{0})}\cap
B(\mathbf{0},K)\big{]}+\mathbb{P}\big{[}\bar{\xi}(\mathcal{P}_{\phi(x_{0})},K)-\underline{\xi}(\mathcal{P}_{\phi(x_{0})},K)\geq\varepsilon\big{]},$

where $\bar{\xi}$ and $\underline{\xi}$ are defined in Definition C.1. By the
stabilization assumption, we can choose $K>0$ so that second term is less than
$\delta$, for any $\delta>0$, by taking $K$ large. By the previous lemma, the
first term converges to 0 as $N\to\infty$. Since $\delta>0$ is arbitrary, it
follows that

$\xi\left(\mathbf{0},N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\right)\overset{d}{\longrightarrow}\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right).$

The assumption on the boundedness of the $p$-th moment implies that
$\\{\xi\left(\mathbf{0},N^{1/d}(\mathcal{P}_{N\phi_{N}}-x_{0})\right)\\}_{N\geq
1}$ is uniformly integrable; hence its expectation converges to the
expectation of $\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)$. ∎

######  Lemma C.3.

Suppose $x_{0}\neq y_{0}$ are Lebesgue points of $f_{1},\ldots,f_{M}$. There
exist homogeneous Poisson processes $\mathcal{P}_{\phi(x_{0})}$ and
$\mathcal{P}_{\phi(y_{0})}$ coupled with $\mathcal{P}_{N\phi_{N}}$ such that
$\mathcal{P}_{\phi(x_{0})}$ is independent of $\mathcal{P}_{\phi(y_{0})}$ and
for all $K>0$:

$\lim_{N\to\infty}\mathbb{P}\left[N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\cap
B(\mathbf{0},K)=\mathcal{P}_{\phi(x_{0})}\cap B(\mathbf{0},K)\right]\ =1,$
$\lim_{N\to\infty}\mathbb{P}\left[N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\cap
B(\mathbf{0},K)=\mathcal{P}_{\phi(y_{0})}\cap B(\mathbf{0},K)\right]\ =1.$

###### Proof.

Let $\mathcal{P}$ be a homogeneous Poisson processes of rate 1 on
$\mathbb{R}^{d}\times[0,\infty)$. Let $\mathcal{P}_{N\phi_{N}}$ be the
projection of the set

$\left\\{(x,t)\in\mathcal{P}:t\leq N\phi_{N}(x)\right\\}$

onto the $x$ space, i.e., $(x,t)\mapsto x$. Let $\mathcal{Q}$ be an
independent copy of $\mathcal{P}$. Let
$F_{x_{0}}:=\\{z\in\mathbb{R}^{d}:\|x_{0}-z\|<\|y_{0}-z\|\\}$ be the half-
space of points in $\mathbb{R}^{d}$ closer to $x_{0}$ than to $y_{0}$ and let
$F_{y_{0}}:=\\{z\in\mathbb{R}^{d}:\|y_{0}-z\|<\|x_{0}-z\|\\}$ be the half-
space of points in $\mathbb{R}^{d}$ closer to $y_{0}$ than to $x_{0}$. Let
$\mathcal{P}_{\phi(x_{0})}$ be the image of the point set

$\\{(x,t)\in\mathcal{P}\cap
F_{x_{0}}\times[0,N\phi(x_{0})]\\}\cup\\{(x,t)\in\mathcal{Q}\cap
F_{y_{0}}\times[0,N\phi(x_{0})]\\}$

under the mapping

$(x,t)\mapsto N^{1/d}(x-x_{0}).$

Note that $\mathcal{P}_{\phi(x_{0})}$ is a homogeneous Poisson process on
$\mathbb{R}^{d}$ with intensity $\phi(x_{0})$ for all $N$. Note that when $N$
is large, $B(x_{0},N^{-1/d}K)\subset F_{x_{0}}$. Hence, the number of points
in
$\left(N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\triangle\mathcal{P}_{\phi(x_{0})}\right)\cap
B(\mathbf{0},K)$ equals the number of points $(x,t)\in\mathcal{P}$ such that
$x\in B(x_{0},N^{-1/d}K)$ and $t$ is between $\phi_{N}(x)$ and $\phi(x_{0})$,
and follows a Poisson distribution with mean

$N\int_{B(x_{0},N^{-1/d}K)}|\phi_{N}(x)-\phi(x_{0})|{\rm d}x\leq
N\int_{B(x_{0},N^{-1/d}K)}(|\phi_{N}(x)-\phi(x)|+|\phi(x)-\phi(x_{0})|){\rm
d}x,$

which converges to 0 as $x_{0}$ is also a Lebesgue point of $\phi$.

Let $\mathcal{P}_{\phi(y_{0})}$ be the image of the point set

$\\{(x,t)\in\mathcal{P}\cap
F_{y_{0}}\times[0,N\phi(y_{0})]\\}\cup\\{(x,t)\in\mathcal{Q}\cap
F_{x_{0}}\times[0,N\phi(y_{0})]\\}$

under the mapping

$(x,t)\mapsto N^{1/d}(x-y_{0}).$

$\mathcal{P}_{\phi(x_{0})}$ and $\mathcal{P}_{\phi(y_{0})}$ are independent as
they are constructed from Poisson processes on disjoint regions of space. The
number of points in
$\left(N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\triangle\mathcal{P}_{\phi(y_{0})}\right)\cap
B(\mathbf{0},K)$ also converges in mean to 0 since it follows a Poisson
distribution with mean

$N\int_{B(y_{0},N^{-1/d}K)}|\phi_{N}(x)-\phi(y_{0})|{\rm d}x\leq
N\int_{B(y_{0},N^{-1/d}K)}(|\phi_{N}(x)-\phi(x)|+|\phi(x)-\phi(y_{0})|){\rm
d}x,$

which converges to 0 as $y_{0}$ is also a Lebesgue point of $\phi$. ∎

######  Lemma C.4.

Suppose $x_{0}\neq y_{0}$ are Lebesgue points of $f_{1},\ldots,f_{M}$ and
$\phi(x_{0}),\phi(y_{0})>0$. Let $\mathcal{P}_{\phi(x_{0})}$ and
$\mathcal{P}_{\phi(y_{0})}$ be independent homogeneous Poisson processes with
intensity $\phi(x_{0})$ and $\phi(y_{0})$ respectively. Suppose $\xi$ is
translation invariant, almost surely stabilizing on
$\mathcal{P}_{\phi(x_{0})}$ and $\mathcal{P}_{\phi(y_{0})}$. Suppose the
moment condition:

$\sup_{N\geq
1}\mathbb{E}\left[\left|\xi\left(x_{0},x_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\right)\xi\left(y_{0},y_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\right)\right|^{p}\right]<\infty$

is satisfied for some $p>1$. Then:

$\xi\left(x_{0},x_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\right)\xi\left(y_{0},y_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\right)\to\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)\xi\left(\mathbf{0},\mathcal{P}_{\phi(y_{0})}\right)$

in distribution and in expectation.

###### Proof.

Consider the coupling of $\mathcal{P}_{N\phi_{N}}$ and
$\mathcal{P}_{\phi(x_{0})},\mathcal{P}_{\phi(y_{0})}$ in the previous lemma.
The same argument in Lemma C.2 shows

$\displaystyle\mathbb{P}\left[\left|\xi\left(x_{0},x_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\right)-\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)\right|>\varepsilon\right]\to
0,$
$\displaystyle\mathbb{P}\left[\left|\xi\left(y_{0},y_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\right)-\xi\left(\mathbf{0},\mathcal{P}_{\phi(y_{0})}\right)\right|>\varepsilon\right]\to
0.$

Hence

$\xi\left(x_{0},x_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{y_{0}}-x_{0})\right)\xi\left(y_{0},y_{0}+N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x_{0}}-y_{0})\right)\overset{d}{\longrightarrow}\xi\left(\mathbf{0},\mathcal{P}_{\phi(x_{0})}\right)\xi\left(\mathbf{0},\mathcal{P}_{\phi(y_{0})}\right).$

Convergence in expectation again follows from uniform integrability. ∎

######  Lemma C.5 (Palm theory for Poisson processes [53, Theorem 1.6]).

Suppose $N>0$, $j\in\mathbb{N}$, and $h(\mathcal{Y},\mathcal{X})$ is a bounded
measurable function defined on all pairs of the form
$(\mathcal{Y},\mathcal{X})$ where $\mathcal{X}\subset\mathbb{R}^{d}$ is
finite, and $\mathcal{Y}\subset\mathcal{X}$, satisfying
$h(\mathcal{Y},\mathcal{X})=0$ when $\mathcal{Y}$ does not contain $j$
elements. Then

$\mathbb{E}\sum_{\mathcal{Y}\subset\mathcal{P}_{Nf}}h(\mathcal{Y},\mathcal{P}_{Nf})=\frac{N^{j}}{j!}\mathbb{E}h(x_{j}^{\prime},x_{j}^{\prime}\cup\mathcal{P}_{Nf}),$

where the sum on the left-hand side is over all subsets
$\mathcal{Y}\subset\mathcal{P}_{Nf}$, and on the right-hand side
$x_{j}^{\prime}$ is a set of $j$ i.i.d. observations from $f$, independent of
$\mathcal{P}_{Nf}$.

####  C.12.2 Proof of Theorem B.1

Write
$g_{N}(\Delta_{i}):=2\sum_{p=1}^{M}\frac{N_{p}}{N}K(\Delta_{i},p)-\sum_{p,q=1}^{M}\frac{N_{p}N_{q}}{N^{2}}K(p,q)$.
By standard U-statistics projection theory [21, Lemma D.4], there exists a
constant $C$ such that for $n\geq 2$,

$\mathbb{E}\left[\left(\frac{1}{\sqrt{n}(n-1)}\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}g_{N}(\Delta_{i})\right)^{2}\Big{|}n\right]\leq\frac{C}{n}.$

Since $n\sim{\rm Poisson}(N)$, $\mathbb{P}(n<2)$ decays exponentially as
$N\to\infty$, and recall that we set $\frac{1}{\sqrt{n}(n-1)}\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})$ to 0 when $n<2$. This implies

$\mathbb{E}\left[\left(\frac{1}{\sqrt{n}(n-1)}\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{n}}\sum_{i=1}^{n}g_{N}(\Delta_{i})\right)^{2}\right]\to
0,\ {\rm as\ }N\to\infty.$

Hence $\frac{1}{n-1}\sum_{i\neq
j}K(\Delta_{i},\Delta_{j})=\sum_{i=1}^{n}g_{N}(\Delta_{i})+o_{p}(\sqrt{n})$.
It suffices to derive the asymptotic distribution of

$H_{n}:=\frac{1}{\sqrt{N}}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}g_{N}(\Delta_{i}).$
(C.10)

Let $\mathcal{F}_{n}=\sigma(n,Z_{1},\ldots,Z_{n})$ be the $\sigma$-algebra
generated by the unlabelled data and the number of total observations. The
strategy of the proof goes as follows: First, we show that given
$\mathcal{F}_{n}$, $H_{n}$ — centered by its conditional mean — converges
conditionally to a normal distribution with constant variance, i.e.,
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\mid\mathcal{F}_{n}\overset{d}{\to}N(0,\kappa_{1}^{2})$.
Second, we show that
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{2}^{2})$.
Finally, a simple argument using the characteristic function yields that
$H_{n}-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{1}^{2}+\kappa_{2}^{2})$.
Observe that, for $t\in\mathbb{R}$,

$\displaystyle\mathbb{E}\left[e^{it\left(H_{n}-\mathbb{E}(H_{n})\right)}\right]$
$\displaystyle=\mathbb{E}\left[e^{it\left(H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\right)}e^{it\left(\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\right)}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[e^{it\left(H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\right)}\big{|}\mathcal{F}_{n}\right]e^{it\left(\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\right)}\right].$

Note that for any $t\in\mathbb{R}$, by continuous mapping,

$\mathbb{E}\left[e^{it\left(H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\right)}\big{|}\mathcal{F}_{n}\right]e^{it\left(\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\right)}\overset{d}{\longrightarrow}e^{-\frac{1}{2}\kappa_{1}^{2}t}\cdot
e^{it\cdot N(0,\kappa_{2}^{2})}.$

Hence, by choosing an almost sure representative and applying the dominated
convergence theorem, we have

$\mathbb{E}\left[\mathbb{E}\left[e^{it\left(H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\right)}\big{|}\mathcal{F}_{n}\right]e^{it\left(\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\right)}\right]\to\mathbb{E}\left[e^{-\frac{1}{2}\kappa_{1}^{2}t}\cdot
e^{it\cdot
N(0,\kappa_{2}^{2})}\right]=e^{-\frac{1}{2}(\kappa_{1}^{2}+\kappa_{2}^{2})t}.$

This implies
$H_{n}-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{1}^{2}+\kappa_{2}^{2})$.

Step 1. We wil first show that
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\mid\mathcal{F}_{n}\overset{d}{\to}N(0,\kappa_{1}^{2})$.

We first compute the limiting variance of $H_{n}$ given $\mathcal{F}_{n}$.
Observe that

$\displaystyle{\rm Var}(H_{n}|\mathcal{F}_{n})$
$\displaystyle=\frac{1}{N}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}^{2}}{\rm
Var}(K(\Delta_{i},\Delta_{j})|\mathcal{F}_{n})$ (C.11a)
$\displaystyle\quad+\frac{1}{N}\sum_{i\neq
j:(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{j}}{\rm
Var}(K(\Delta_{i},\Delta_{j})|\mathcal{F}_{n})$ (C.11b)
$\displaystyle\quad+\frac{1}{N}\sum_{\begin{subarray}{c}i,j,l{\rm\
distinct:}\\\
(Z_{i},Z_{j}),(Z_{i},Z_{l})\in\mathcal{E}(\mathcal{G}_{n})\end{subarray}}\frac{1}{d_{i}^{2}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{i},\Delta_{l})|\mathcal{F}_{n})$
(C.11c) $\displaystyle\quad+\frac{1}{N}\sum_{\begin{subarray}{c}i,j,l{\rm\
distinct:}\\\
(Z_{i},Z_{j}),(Z_{j},Z_{l})\in\mathcal{E}(\mathcal{G}_{n})\end{subarray}}\frac{1}{d_{i}d_{j}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{j},\Delta_{l})|\mathcal{F}_{n})$
(C.11d) $\displaystyle\quad+\frac{1}{N}\sum_{\begin{subarray}{c}i,j,l{\rm\
distinct:}\\\
(Z_{i},Z_{j}),(Z_{l},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})\end{subarray}}\frac{1}{d_{i}d_{l}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{l},\Delta_{i})|\mathcal{F}_{n})$
(C.11e) $\displaystyle\quad+\frac{1}{N}\sum_{\begin{subarray}{c}i,j,l{\rm\
distinct:}\\\
(Z_{i},Z_{j}),(Z_{l},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})\end{subarray}}\frac{1}{d_{i}d_{l}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{l},\Delta_{j})|\mathcal{F}_{n})$
(C.11f) $\displaystyle\quad-\frac{2}{N}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),g_{N}(\Delta_{i})|\mathcal{F}_{n})$ (C.11g)
$\displaystyle\quad-\frac{2}{N}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),g_{N}(\Delta_{j})|\mathcal{F}_{n})$ (C.11h)
$\displaystyle\quad+\frac{1}{N}\sum_{i=1}^{n}{\rm
Var}(g_{N}(\Delta_{i})|\mathcal{F}_{n}).$ (C.11i)

We will show the convergence of each term above. Let
$\phi_{N}(z)=\sum_{i=1}^{M}\frac{N_{i}}{N}f_{i}(z)$ be the marginal density of
$Z$, and $\phi(z):=\sum_{i=1}^{M}\pi_{i}f_{i}(z)$ be its limit. Then, for
$p\in\\{1,\ldots,M\\}$,

$\displaystyle\left|\frac{\frac{N_{p}}{N}f_{p}(z)}{\phi_{N}(z)}-\frac{\pi_{p}f_{p}(z)}{\phi(z)}\right|$
$\displaystyle=\left|\frac{\frac{N_{p}}{N}f_{p}(z)\phi(z)-\pi_{p}f_{p}(z)\phi_{N}(z)}{\phi(z)\phi_{N}(z)}\right|$
$\displaystyle=\left|\sum_{i=1}^{M}\frac{\frac{N_{p}}{N}f_{p}(z)\pi_{i}f_{i}(z)-\pi_{p}f_{p}(z)\frac{N_{i}}{N}f_{i}(z)}{\phi(z)\phi_{N}(z)}\right|$
$\displaystyle\leq\sum_{i=1}^{M}\frac{|\frac{N_{p}}{N}-\pi_{p}|f_{p}(z)\pi_{i}f_{i}(z)+\pi_{p}f_{p}(z)|\pi_{i}-\frac{N_{i}}{N}|f_{i}(z)}{\phi(z)\phi_{N}(z)}$
$\displaystyle=\sum_{i=1}^{M}\left(\frac{|\frac{N_{p}}{N}-\pi_{p}|}{\frac{N_{p}}{N}}\frac{\frac{N_{p}}{N}f_{p}(z)}{\phi_{N}(z)}\frac{\pi_{i}f_{i}(z)}{\phi(z)}+\frac{|\pi_{i}-\frac{N_{i}}{N}|}{\frac{N_{i}}{N}}\frac{\pi_{p}f_{p}(z)}{\phi(z)}\frac{\frac{N_{i}}{N}f_{i}(z)}{\phi_{N}(z)}\right)$
$\displaystyle\leq\sum_{i=1}^{M}\left(\frac{|\frac{N_{p}}{N}-\pi_{p}|}{\frac{N_{p}}{N}}+\frac{|\pi_{i}-\frac{N_{i}}{N}|}{\frac{N_{i}}{N}}\right)\to
0$

uniformly over $z\in\\{z:\phi(z)>0\\}$. Observe that

${\rm
Var}(K(\Delta_{i},\Delta_{j})|\mathcal{F}_{n})=\underbrace{\sum_{p,q=1}^{M}K^{2}(p,q)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{q}}{N}f_{q}(Z_{j})}{\phi_{N}(Z_{j})}-\left(\sum_{p,q=1}^{M}K(p,q)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{q}}{N}f_{q}(Z_{j})}{\phi_{N}(Z_{j})}\right)^{2}}_{{\rm
denoted\ by\ }h_{N}(Z_{i},Z_{j})}.$

By the uniform convergence of $\frac{\frac{N_{i}}{N}f_{i}(z)}{\phi_{N}(z)}$ to
$\frac{\pi_{i}f_{i}(z)}{\phi(z)}$, we know that $h_{N}$ converges uniformly to
$h$, a bounded function defined by:

$h(Z_{i},Z_{j}):=\sum_{p,q=1}^{M}K^{2}(p,q)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{q}f_{q}(Z_{j})}{\phi(Z_{j})}-\left(\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{q}f_{q}(Z_{j})}{\phi(Z_{j})}\right)^{2}.$

To show the convergence of the first term in the expansion of ${\rm
Var}(H_{n}|\mathcal{F}_{n})$, i.e.,

$\frac{1}{N}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}^{2}}\cdot
h_{N}(Z_{i},Z_{j})=:f(Z),$

where $Z:=(Z_{1},\ldots,Z_{n})$. We will show that there exists a constant
$\mu$ such that $\mathbb{E}f(Z)\to\mu$ and $\mathbb{E}f^{2}(Z)\to\mu^{2}$,
which implies $f(Z)\overset{L^{2}}{\longrightarrow}\mu$.

By Palm theory (Lemma C.5),

$\mathbb{E}f(Z)=\int\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h_{N}(z,y)\right]\phi_{N}(z){\rm
d}z.$ (C.12)

To show the convergence of the integrand, we can first replace $h_{N}$ by its
uniform limit $h$.

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h_{N}(z,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h(z,y)\right].$

Note that $(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))$
implies that $y$ is a neighbor of $z$, and Lemma C.1 together with the graph
being stabilizing on $\mathcal{P}_{\phi(z)}$ implies
$\sup_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\|y-z\|=O_{p}\left(N^{-1/d}\right)$.
If $z$ is a continuity point of $f_{1},\cdots,f_{M}$, then it is also a
continuity point of $h(z,\cdot)$. Hence we can further replace $y$ by $z$ to
get

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h(z,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{y:(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h(z,z)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\right]h(z,z)$
$\displaystyle=\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(z)}^{\mathbf{0}}))}\right]h(z,z)=\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(z,z).$

The last line follows from Lemma C.2 and $\mathcal{G}$ being translation and
scale invariant. Now, by dominated convergence,

$\mathbb{E}f(Z)\to\int\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(z,z)\phi(z){\rm
d}z.$

Now we consider $\mathbb{E}f^{2}(Z)$. Write
$\xi(x,\mathcal{X}):=\sum_{y:(x,y)\in\mathcal{E}(\mathcal{G}(\mathcal{X}\cup\\{x\\}))}\frac{1}{d_{x}^{2}}h_{N}(x,y)$.
Then
$f(Z)=\frac{1}{N}\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})$,
and

$\displaystyle\mathbb{E}f^{2}(Z)$
$\displaystyle=\frac{1}{N^{2}}\mathbb{E}\left[\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})\right]+\frac{2}{N^{2}}\mathbb{E}\left[\sum_{\\{x,y\\}\subset\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})\xi(y,\mathcal{P}_{N\phi_{N}})\right].$

Here $\sum_{\\{x,y\\}\subset\mathcal{P}_{N\phi_{N}}}$ means summing over all
$\frac{n(n-1)}{2}$ size-2 subsets of $\mathcal{P}_{N\phi_{N}}$. The first term
converges to 0 since $\xi$ is bounded and the number of points in
$\mathcal{P}_{N\phi_{N}}$ follows ${\rm Poisson}(N)$. The second term, by Palm
theory, equals

$\int\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]\phi_{N}(x)\phi_{N}(y){\rm
d}x{\rm d}y.$

By the uniform convergence of $h_{N}$ to $h$, assuming $x,y$ are continuity
point of $f_{1},\ldots,f_{M}$, the integrand has a limit given by

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\sum_{w:(x,w)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{x}^{2}}h_{N}(x,w)\sum_{w:(y,w)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{y}^{2}}h_{N}(y,w)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\sum_{w:(x,w)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{x}^{2}}h(x,x)\sum_{w:(y,w)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{y}^{2}}h(y,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d(x,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d(y,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\right]h(x,x)h(y,y)$
$\displaystyle=\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(x)}^{\mathbf{0}}))}\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(y)}^{\mathbf{0}}))}\right]h(x,x)h(y,y)=\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]^{2}h(x,x)h(y,y).$

The last line follows from Lemma C.4 and $\mathcal{G}$ being translation and
scale invariant. Now, by dominated convergence,

$\mathbb{E}f^{2}(Z)\to\int\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]^{2}h(x,x)h(y,y)\phi(x)\phi(y){\rm
d}x{\rm d}y.$

Hence $f(Z)$ converges in $L^{2}$ to
$\int\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(z,z)\phi(z){\rm
d}z$.

With the distribution of $\tilde{\Delta},\tilde{Z}$ defined in Section 2.1, if
$\tilde{\Delta},\tilde{\Delta}^{\prime}$ are independently drawn from the
conditional distribution of $\tilde{\Delta}\mid\tilde{Z}=z$, then $h(z,z)={\rm
Var}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]$. Hence the first
term in ${\rm Var}(H_{n}|\mathcal{F}_{n})$ converges to:

$\displaystyle\frac{1}{N}\sum_{i\neq
j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}^{2}}{\rm
Var}(K(\Delta_{i},\Delta_{j})|\mathcal{F}_{n})\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]\int{\rm
Var}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]\phi(z){\rm d}z.$

The same strategy applied to the second term (C.11b) yields:

$\displaystyle\quad\frac{1}{N}\sum_{i\neq
j:(Z_{i},Z_{j}),(Z_{j},Z_{i})\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{i}d_{j}}{\rm
Var}(K(\Delta_{i},\Delta_{j})|\mathcal{F}_{n})$
$\displaystyle\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{y:(y,\mathbf{0}),(\mathbf{0},y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]\int{\rm
Var}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]\phi(z){\rm d}z.$

The convergence of the third term (C.11c) can be handled similarly. Note that

$\displaystyle{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{i},\Delta_{l})|\mathcal{F}_{n})=\sum_{p,q,r=1}^{M}K(p,q)K(p,r)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{q}}{N}f_{q}(Z_{j})}{\phi_{N}(Z_{j})}\frac{\frac{N_{r}}{N}f_{r}(Z_{l})}{\phi_{N}(Z_{l})}$
$\displaystyle-\sum_{p,q=1}^{M}K(p,q)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{q}}{N}f_{q}(Z_{j})}{\phi_{N}(Z_{j})}\sum_{p,r=1}^{M}K(p,r)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{r}}{N}f_{r}(Z_{l})}{\phi_{N}(Z_{l})}.$

Denote the right-hand side of the above display by $h_{N}(Z_{i},Z_{j},Z_{l})$.
Then $h_{N}$ converges uniformly to $h$ defined as:

$\displaystyle
h(Z_{i},Z_{j},Z_{l})=\sum_{p,q,r=1}^{M}K(p,q)K(p,r)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{q}f_{q}(Z_{j})}{\phi(Z_{j})}\frac{\pi_{r}f_{r}(Z_{l})}{\phi(Z_{l})}$
$\displaystyle-\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{q}f_{q}(Z_{j})}{\phi(Z_{j})}\sum_{p,r=1}^{M}K(p,r)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{r}f_{r}(Z_{l})}{\phi(Z_{l})}.$

With $\xi(z,\mathcal{P}_{N\phi_{N}}):=\frac{1}{d_{x}^{2}}\sum_{x\neq
y:(z,x),(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}))}h_{N}(z,x,y)$,
(C.11c) can be written as

$f(Z)=\frac{1}{N}\sum_{z\in\mathcal{P}_{N\phi_{N}}}\xi(z,\mathcal{P}_{N\phi_{N}}).$

By Palm theory,

$\mathbb{E}f(Z)=\int\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{x\neq
y:(z,x),(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h_{N}(z,x,y)\right]\phi_{N}(z){\rm
d}z.$ (C.13)

Since $h_{N}$ converges uniformly to $h$, using the same argument as before,
if $z$ is a continuity point of $f_{1},\ldots,f_{M}$, then

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{x\neq
y:(z,x),(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h_{N}(z,x,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{x\neq
y:(z,x),(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h(z,x,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\sum_{x\neq
y:(z,x),(z,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}h(z,z,z)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{d(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))(d(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))-1)}{d^{2}(z,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))}\right]h(z,z,z)$
$\displaystyle=\mathbb{E}\left[1-\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(z)}^{\mathbf{0}}))}\right]h(z,z,z)=\mathbb{E}\left[1-\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(z,z,z).$

Hence by dominated convergence theorem applied to (C.13),

$\mathbb{E}f(Z)\to\int\mathbb{E}\left[1-\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(z,z,z)\phi(z){\rm
d}z.$

Note that

$\displaystyle\mathbb{E}f^{2}(Z)$
$\displaystyle=\frac{1}{N^{2}}\mathbb{E}\left[\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})\right]+\frac{2}{N^{2}}\mathbb{E}\left[\sum_{\\{x,y\\}\subset\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})\xi(y,\mathcal{P}_{N\phi_{N}})\right].$

The first term converges to 0 again as $\xi$ is bounded and the number of
points in $\mathcal{P}_{N\phi_{N}}$ follows ${\rm Poisson}(N)$. The second
term equals
$\int\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]\phi_{N}(x)\phi_{N}(y){\rm
d}x{\rm d}y$ by Palm theory. By the uniform convergence of $h_{N}$ to $h$,
assuming $x,y$ are continuity point of $f_{1},\ldots,f_{M}$, the integrand has
a limit given by

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\sum_{u\neq
v:(x,u),(x,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{x}^{2}}h_{N}(x,u,v)\sum_{u\neq
v:(y,u),(y,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{y}^{2}}h_{N}(y,u,v)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\sum_{u\neq
v:(x,u),(x,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{x}^{2}}h(x,u,v)\sum_{u\neq
v:(y,u),(y,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{y}^{2}}h(y,u,v)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\sum_{u\neq
v:(x,u),(x,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{x}^{2}}h(x,x,x)\sum_{u\neq
v:(y,u),(y,v)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{1}{d_{y}^{2}}h(y,y,y)\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{d(x,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))-1}{d(x,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\frac{d(y,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))-1}{d(y,\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{x,y}))}\right]h(x,x,x)h(y,y,y)$
$\displaystyle=\mathbb{E}\left[\frac{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(x)}^{\mathbf{0}}))-1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(x)}^{\mathbf{0}}))}\frac{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(y)}^{\mathbf{0}}))-1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{\phi(y)}^{\mathbf{0}}))}\right]h(x,x,x)h(y,y,y)$
$\displaystyle=\mathbb{E}\left[\frac{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))-1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]^{2}h(x,x,x)h(y,y,y).$

The last line again follows from Lemma C.4 and $\mathcal{G}$ being translation
and scale invariant. Now, by dominated convergence,

$\mathbb{E}f^{2}(Z)\to\int\mathbb{E}\left[\frac{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))-1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]^{2}h(x,x,x)h(y,y,y)\phi(x)\phi(y){\rm
d}x{\rm d}y.$

Hence $f(Z)$ converges in $L^{2}$ to
$\int\mathbb{E}\left[\frac{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))-1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]h(x,x,x)\phi(x){\rm
d}x$.

Therefore, (C.11c) converges in $L^{2}$ to

$\displaystyle\quad\frac{1}{N}\sum_{\begin{subarray}{c}i,j,l{\rm\
distinct:}\\\
(Z_{i},Z_{j}),(Z_{i},Z_{l})\in\mathcal{E}(\mathcal{G}_{n})\end{subarray}}\frac{1}{d_{i}^{2}}{\rm
Cov}(K(\Delta_{i},\Delta_{j}),K(\Delta_{i},\Delta_{l})|\mathcal{F}_{n})$
$\displaystyle\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[1-\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]\int
h(z,z,z)\phi(z){\rm d}z$
$\displaystyle=\mathbb{E}\left[1-\frac{1}{d(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\right]\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta},\tilde{\Delta}^{\prime\prime})|\tilde{Z}=z)\phi(z){\rm
d}z,$

where $\tilde{\Delta},\tilde{\Delta}^{\prime},\tilde{\Delta}^{\prime\prime}$
are independently drawn from the conditional distribution of
$\tilde{\Delta}\mid\tilde{Z}=z$.

The other terms in the right side of (C.11) converge similarly:

$\displaystyle\quad\eqref{eq:subeq4}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{x\neq
y:(x,\mathbf{0}),(\mathbf{0},y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(x,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta}^{\prime},\tilde{\Delta}^{\prime\prime})|\tilde{Z}=z)\phi(z){\rm
d}z,$
$\displaystyle\quad\eqref{eq:subeq5}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{x\neq
y:(\mathbf{0},x),(y,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta}^{\prime\prime},\tilde{\Delta})|\tilde{Z}=z)\phi(z){\rm
d}z,$
$\displaystyle\quad\eqref{eq:subeq6}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[\sum_{x\neq
y:(x,\mathbf{0}),(y,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(x,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(y,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta}^{\prime\prime},\tilde{\Delta}^{\prime})|\tilde{Z}=z)\phi(z){\rm
d}z,$ $\eqref{eq:subeq7}\overset{L^{2}}{\longrightarrow}-2\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),g(\tilde{\Delta})|\tilde{Z}=z)\phi(z){\rm
d}z,$ $\eqref{eq:subeq8}\overset{L^{2}}{\longrightarrow}-2\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),g(\tilde{\Delta}^{\prime})|\tilde{Z}=z)\phi(z){\rm
d}z,$ $\eqref{eq:subeq9}\overset{L^{2}}{\longrightarrow}\mathbb{E}\left[{\rm
Var}(g(\tilde{\Delta})|\tilde{Z})\right].$

where
$g(\Delta)=2\sum_{p=1}^{M}\pi_{p}K(\Delta,p)-\sum_{p,q=1}^{M}\pi_{p}\pi_{q}K(p,q)$
is the uniform limit in $\Delta$ of $g_{N}(\Delta)$.

Combining the results above, we have ${\rm
Var}(H_{n}|\mathcal{F}_{n})\overset{L^{2}}{\longrightarrow}\kappa_{1}^{2}$,
where

$\displaystyle\kappa_{1}^{2}$ $\displaystyle=(g_{1}+g_{3})\int{\rm
Var}[K(\tilde{\Delta},\tilde{\Delta}^{\prime})|\tilde{Z}=z]\phi(z){\rm d}z$
$\displaystyle\quad+(3-2g_{1}-2g_{3}+g_{2})\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta},\tilde{\Delta}^{\prime\prime})|\tilde{Z}=z)\phi(z){\rm
d}z$ $\displaystyle\quad-4\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),g(\tilde{\Delta})|\tilde{Z}=z)\phi(z){\rm
d}z\;+\;\mathbb{E}\left[{\rm Var}(g(\tilde{\Delta})|\tilde{Z})\right],$

where when simplifying the coefficient of $\int{\rm
Cov}(K(\tilde{\Delta},\tilde{\Delta}^{\prime}),K(\tilde{\Delta},\tilde{\Delta}^{\prime\prime})|\tilde{Z}=z)\phi(z){\rm
d}z$, we have used that

$\displaystyle\quad\mathbb{E}\left[\sum_{x\neq
y:(x,\mathbf{0}),(\mathbf{0},y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(x,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{x:(x,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)-1_{(\mathbf{0},x),(x,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}}{d\left(x,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)d\left(\mathbf{0},\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right],$

and

$\mathbb{E}\left[\sum_{x:(x,\mathbf{0})\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}}))}\frac{1}{d\left(x,\mathcal{G}(\mathcal{P}_{1}^{\mathbf{0}})\right)}\right]=1$

since it is the limit of the expectation of
$\frac{1}{n}\sum_{x\in\mathcal{X}_{n}}\sum_{y:(y,x)\in\mathcal{E}(\mathcal{G}_{n})}\frac{1}{d_{y}}\equiv
1$ [56, Lemma 3.2]. In particular, if the null hypothesis holds true (in which
case $\tilde{\Delta}\perp\\!\\!\\!\perp\tilde{Z}$), then $\kappa_{1}^{2}$
reduces exactly to the asymptotic variance of the numerator derived in Theorem
4, which is distribution-free.

If $\kappa_{1}^{2}=0$, then
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})\mid\mathcal{F}_{n}\overset{d}{\to}N(0,\kappa_{1}^{2})$
follows trivially from the convergence of the conditional variance of
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})$ to 0. We will assume
$\kappa_{1}^{2}>0$ and prove the convergence using Stein’s method based on
dependency graphs (see Theorem C.3 below). First note that
$H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})=\sum_{i=1}^{n}V_{i}$, where

$\displaystyle V_{i}$
$\displaystyle=\frac{1}{\sqrt{N}}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}g_{N}(\Delta_{i})$
$\displaystyle-\mathbb{E}\left[\frac{1}{\sqrt{N}}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}g_{N}(\Delta_{i})\Big{|}\mathcal{F}_{n}\right].$

Construct a graph $\mathcal{D}(\mathcal{G}_{n})$ on $\\{V_{1},\ldots,V_{n}\\}$
as follows: for $i\neq j$, there is an edge between $V_{i}$ and $V_{j}$ in
$\mathcal{D}(\mathcal{G}_{n})$ if and only if there is a path of length $\leq
2$ joining $Z_{i}$ and $Z_{j}$ in $\mathcal{G}_{n}$ (ignoring the direction of
edges in $\mathcal{G}_{n}$). Then for any pair of disjoint sets
$\Gamma_{1},\Gamma_{2}\subset\\{V_{1},\ldots,V_{n}\\}$ such that no edge in
$\mathcal{D}(\mathcal{G}_{n})$ has one endpoint in $\Gamma_{1}$ and the other
in $\Gamma_{2}$, $\\{V_{i}\\}_{i\in\Gamma_{1}}$ is independent of
$\\{V_{i}\\}_{i\in\Gamma_{2}}$ conditioned on $\mathcal{F}_{n}$. This implies
that $\mathcal{D}(\mathcal{G}_{n})$ is a dependency graph [18]. Suppose
$t_{n}$ is the maximum degree of $\mathcal{G}_{n}$. Then the maximal degree in
$\mathcal{D}(\mathcal{G}_{n})$ has an upper bound of $t_{n}^{2}$. We will use
the following CLT for dependency graph (with $p=3$).

######  Theorem C.3 (18, Theorem 2.7).

Suppose $\\{X_{i}\\}_{i\in\mathcal{V}}$ are random variables indexed by
vertices of a dependency graph, whose maximal degree is $D$. Set
$W=\sum_{i\in\mathcal{V}}X_{i}$. If ${\rm Var}(W)=1$, $\mathbb{E}X_{i}=0$, and
$\mathbb{E}|X_{i}|^{p}\leq\theta^{p}$ for $i\in\mathcal{V}$, then we have

$\sup_{z\in\mathbb{R}}|\mathbb{P}(W\leq z)-\Phi(z)|\leq
75D^{5(p-1)}|\mathcal{V}|\theta^{p}.$

Note that

$\displaystyle\quad\left|\frac{1}{\sqrt{N}}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}g_{N}(\Delta_{i})\right|^{3}$
$\displaystyle\leq\left(\left|\frac{1}{\sqrt{N}}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})\right|^{3}+\left|\frac{1}{\sqrt{N}}g_{N}(\Delta_{i})\right|^{3}\right)(1+1)(1+1)$
$\displaystyle\leq
4\left(\left|\frac{1}{\sqrt{N}}\|K\|_{\infty}\right|^{3}+\left|\frac{1}{\sqrt{N}}\|g_{N}\|_{\infty}\right|^{3}\right)\;\lesssim\;\frac{1}{N^{3/2}}.$

Hence

$\displaystyle\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\frac{H_{n}-\mathbb{E}(H_{n}|\mathcal{F}_{n})}{{\rm
Var}(H_{n}|\mathcal{F}_{n})}\leq
z\Big{|}\mathcal{F}_{n}\right)-\Phi(z)\right|$
$\displaystyle=\;\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\sum_{i=1}^{n}\frac{V_{i}}{{\rm
Var}(H_{n}|\mathcal{F}_{n})}\leq
z\Big{|}\mathcal{F}_{n}\right)-\Phi(z)\right|$ $\displaystyle\lesssim
75D^{5(3-1)}\cdot n\cdot\frac{1}{N^{3/2}({\rm
Var}(H_{n}|\mathcal{F}_{n}))^{3}}$
$\displaystyle\leq\frac{75t_{n}^{20}}{n^{1/2}({\rm
Var}(H_{n}|\mathcal{F}_{n}))^{3}}\cdot\frac{n^{3/2}}{N^{3/2}}\overset{p}{\to}0.$

Step 2. We will now show that
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{2}^{2})$.

Note that $\mathbb{E}(H_{n}|\mathcal{F}_{n})$ can be written as:

$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\frac{1}{d_{i}}\sum_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\sum_{p,q=1}^{M}K(p,q)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}\frac{\frac{N_{q}}{N}f_{q}(Z_{j})}{\phi_{N}(Z_{j})}-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\sum_{p=1}^{M}g_{N}(p)\frac{\frac{N_{p}}{N}f_{p}(Z_{i})}{\phi_{N}(Z_{i})}.$

Since $\sqrt{N}\left(\frac{N_{p}}{N}-\pi_{p}\right)\to 0$, $p=1,\ldots,M$, the
difference between $\frac{\frac{N_{p}}{N}f_{p}(z)}{\phi_{N}(z)}$ and
$\frac{\pi_{p}f_{p}(z)}{\phi(z)}$, and the difference between

$g_{N}(z):=2\sum_{p=1}^{M}\frac{N_{p}}{N}K(z,p)-\sum_{p,q=1}^{M}\frac{N_{p}N_{q}}{N^{2}}K(p,q)$

and

$g(z):=2\sum_{p=1}^{M}\pi_{p}K(z,p)-\sum_{p,q=1}^{M}\pi_{p}\pi_{q}K(p,q)$

are $o(\frac{1}{\sqrt{N}})$, uniformly in $z$. Hence the above statistic has
the same limit in probability as:

$\displaystyle\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\frac{1}{d_{i}}\sum_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})}\frac{\pi_{q}f_{q}(Z_{j})}{\phi(Z_{j})}-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(Z_{i})}{\phi(Z_{i})},$

We can write the above statistic as
$\frac{1}{\sqrt{N}}\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})$,
where

$\xi(x,\mathcal{P}_{N\phi_{N}}):=\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(y)}{\phi(y)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}.$
(C.14)

From the radius of stabilization for $\mathcal{G}$, we also have the same
radius of stabilization [55, 54] $R=R(N,x)$ for
$\xi(x,\mathcal{P}_{N\phi_{N}})$ for which the points outside $B(x,N^{-1/d}R)$
cannot impact the value of $\xi(x,\mathcal{P}_{N\phi_{N}})$, i.e.,

$\xi\left(x,[\mathcal{P}_{N\phi_{N}}\cap
B(x,N^{-1/d}R)]\cup\mathcal{X}\right)=\xi\left(x,\mathcal{P}_{N\phi_{N}}\cap
B(x,N^{-1/d}R)\right),$ (C.15)

for all finite $\mathcal{X}\subset A\backslash B(x,N^{-1/d}R)$.

The power-law stabilization of a function
$\xi:(x,\mathcal{X})\mapsto\mathbb{R}$, where $x\in\mathbb{R}^{d}$ and
$\mathcal{X}\subset\mathbb{R}^{d}$ is a finite set, can be defined in a
similar fashion as the power-law stabilization of a graph (see (B.3)):

######  Definition C.2 (Power-law stabilizing functions [55]).

Let $\xi(x,\mathcal{X})$ be a real-valued function defined for any
$x\in\mathbb{R}^{d}$ and any finite set $\mathcal{X}\subset\mathbb{R}^{d}$
such that $R\equiv R(N,x)$ is a radius of stabilization as defined in (C.15).
Let

$\tau(t):=\sup_{N\geq 1,x\in A}\mathbb{P}[R(N,x)>t],\qquad\mbox{for}\;\;t>0.$

Then, $\xi$ is said to be:

  1. 1.

power-law stabilizing of order $q$ with respect to $\phi_{N}$ if $\sup_{t\geq
1}t^{q}\tau(t)<\infty$,

  2. 2.

exponentially stabilizing with respect to $\phi_{N}$ if
$\limsup_{t\to\infty}t^{-1}\log\tau(t)<0$.

Let

$T:=\frac{1}{\sqrt{N}}\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})$

with $\xi(x,\mathcal{P}_{N\phi_{N}})$ defined in (C.14). The goal is to show
the variance of $T$ converges. By Palm theory,

$\mathbb{E}[T]=\sqrt{N}\int\mathbb{E}\xi(x,\mathcal{P}_{N\phi_{N}})\phi(x){\rm
d}x,$ $\displaystyle\mathbb{E}[T^{2}]$
$\displaystyle=\frac{1}{N}\mathbb{E}\left(\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})\right)+\frac{2}{N}\mathbb{E}\left[\sum_{\\{x,y\\}\subset\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})\xi(y,\mathcal{P}_{N\phi_{N}})\right]$
$\displaystyle=\int\mathbb{E}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})\phi_{N}(x){\rm
d}x+N\int\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]\phi_{N}(x)\phi_{N}(y){\rm
d}x{\rm d}y$

Hence,

$\displaystyle{\rm Var}[T]$
$\displaystyle=\int\mathbb{E}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})\phi_{N}(x){\rm
d}x$ (C.16)
$\displaystyle+N\int\left(\mathbb{E}\left[\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]-\mathbb{E}\xi(x,\mathcal{P}_{N\phi_{N}})\mathbb{E}\xi(y,\mathcal{P}_{N\phi_{N}})\right)\phi_{N}(x)\phi_{N}(y){\rm
d}x{\rm d}y.$

The convergence of the first term above is easy to see. Similar to previous
arguments, by the coupling lemma C.1 and the $\mathcal{P}_{\lambda}$
stabilization assumption,
$(x,y)\in\mathcal{E}(\mathcal{G}(\mathcal{P}_{N\phi_{N}}^{z}))$ implies that
$y$ is a neighbor of $x$, and lies within $O_{p}\left(N^{-1/d}\right)$
distance to $x$. If $x$ is a continuity point of $f_{q}$, for $q=1,\ldots,M$,
then together with boundedness of

$\xi^{2}(x,\mathcal{P}_{N\phi_{N}})=\left(\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(y)}{\phi(y)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2},$

we can replace $f_{q}(y)$ and $\phi(y)$ in the expression above by $f_{q}(x)$
and $\phi(x)$. Thus,

$\displaystyle\lim_{N\to\infty}\mathbb{E}\xi^{2}(x,\mathcal{P}_{N\phi_{N}})$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\left(\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(y)}{\phi(y)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2}\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\left(\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2}\right]$
$\displaystyle=\left(\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2}.$

To show the convergence of the second term in (C.16), we first perform a
change of variable $y=x+N^{-1/d}z$:

$\displaystyle\quad
N\int\left(\mathbb{E}\Big{[}\xi(x,\mathcal{P}_{N\phi_{N}}^{y})\xi(y,\mathcal{P}_{N\phi_{N}}^{x})\right]-\mathbb{E}\xi(x,\mathcal{P}_{N\phi_{N}})\mathbb{E}\xi(y,\mathcal{P}_{N\phi_{N}})\Big{)}\phi_{N}(x)\phi_{N}(y){\rm
d}x{\rm d}y=$ (C.17)
$\displaystyle\int\Big{(}\mathbb{E}\big{[}\underbrace{\xi(x,\mathcal{P}_{N\phi_{N}}^{x+N^{-\frac{1}{d}}z})}_{{\rm
denoted\ by\
}X}\underbrace{\xi(x+N^{-\frac{1}{d}}z,\mathcal{P}_{N\phi_{N}}^{x})}_{Z}\big{]}-\mathbb{E}\underbrace{\xi(x,\mathcal{P}_{N\phi_{N}})}_{X^{\prime}}\mathbb{E}\underbrace{\xi(x+N^{-\frac{1}{d}}z,\mathcal{P}_{N\phi_{N}})}_{Z^{\prime}}\Big{)}\phi_{N}(x)\phi_{N}(x+N^{-\frac{1}{d}}z){\rm
d}x{\rm d}z.$

When taking the limit as $N\to\infty$, we want to pass the limit inside the
integral. We will show that the integrand is small for large $|z|$, uniformly
in $N$. More specifically, there exists $C_{1}>0$ free of $N,x,z$ such that

$|\mathbb{E}[XZ]-\mathbb{E}[X^{\prime}]\mathbb{E}[Z^{\prime}]|\leq
C_{1}(|z|^{-d-1/C_{1}}\wedge 1),$ (C.18)

for all $N\geq 1,x\in\mathbb{R}^{d},z\in\mathbb{R}^{d}$, where
$X,Z,X^{\prime},Z^{\prime}$ are defined in (C.17). Let
$\tilde{X}:=X1_{R(x,N)\leq|z|/3}$,
$\tilde{Z}:=Z1_{R(x+N^{-1/d}z,N)\leq|z|/3}$. Then $\tilde{X}$ and $\tilde{Z}$
are independent because they are determined by the points of
$\mathcal{P}_{N\phi_{N}}$ in $B(x,N^{-\frac{1}{d}}|z|/3)$,
$B(x+N^{-\frac{1}{d}}z,N^{-\frac{1}{d}}|z|/3)$ respectively. So
$\mathbb{E}[\tilde{X}\tilde{Z}]=\mathbb{E}[\tilde{X}]\mathbb{E}[\tilde{Z}]$,
and

$\mathbb{E}[XZ]=\mathbb{E}[\tilde{X}]\mathbb{E}[\tilde{Z}]+\mathbb{E}[\tilde{X}(Z-\tilde{Z})]+\mathbb{E}[(X-\tilde{X})Z]$

while

$\mathbb{E}[X^{\prime}]\mathbb{E}[Z^{\prime}]=\mathbb{E}[\tilde{X}]\mathbb{E}[\tilde{Z}]+\mathbb{E}[\tilde{X}]\mathbb{E}[Z^{\prime}-\tilde{Z}]+\mathbb{E}[X^{\prime}-\tilde{X}]\mathbb{E}[Z^{\prime}].$

Note that the absolute values of
$X,Z,\tilde{X},\tilde{Z},X^{\prime},Z^{\prime}$ are all bounded by a constant
$C>0$ depending only on the kernel $K$, $\pi$, and
$\sup_{n}\\{\frac{t_{n}}{r_{n}}\\}$. Hence

$|\mathbb{E}[(X-\tilde{X})Z]|\leq
C\mathbb{E}|X-\tilde{X}|=C\mathbb{E}|X|1_{R(x,N)>|z|/3}\leq
C^{2}\mathbb{P}(R(x,N)>|z|/3).$

If we have power-law stabilization of order $q>d$, then there exists $C_{2}>0$
such that $C^{2}\mathbb{P}(R(x,N)>|z|/3)\leq C_{2}(|z|^{-d+1/C_{2}}\wedge 1)$.
Similarly, $\mathbb{E}[\tilde{X}(Z-\tilde{Z})]$,
$\mathbb{E}[\tilde{X}]\mathbb{E}[Z^{\prime}-\tilde{Z}]$, and
$\mathbb{E}[(X-\tilde{X})Z]$ are all controlled by
$C_{2}(|z|^{-d+1/C_{2}}\wedge 1)$. Hence (C.18) is proved. This implies that
the contribution of $|z|>K$ to the integral in (C.17) is sufficiently small
for large $K$, and hence we can pass the limit into the integral by dominated
convergence. It remains to compute the limit of the integrand:

$\Big{(}\mathbb{E}\big{[}\xi(x,\mathcal{P}_{N\phi_{N}}^{x+N^{-\frac{1}{d}}z})\xi(x+N^{-\frac{1}{d}}z,\mathcal{P}_{N\phi_{N}}^{x})\big{]}-\mathbb{E}\xi(x,\mathcal{P}_{N\phi_{N}})\mathbb{E}\xi(x+N^{-\frac{1}{d}}z,\mathcal{P}_{N\phi_{N}})\Big{)}\phi_{N}(x)\phi_{N}(x+N^{-\frac{1}{d}}z)$

as $N\to\infty$. We show each term above converges. Suppose $x$ is a
continuity point of $f_{p}$, $p=1,\ldots,M$. With the same coupling technique
used previously, we have:

$\displaystyle\lim_{N\to\infty}\mathbb{E}\xi(x,\mathcal{P}_{N\phi_{N}})$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(y)}{\phi(y)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right]$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right]$
$\displaystyle=\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}.$

Lemma C.1 shows that $N^{1/d}(\mathcal{P}_{N\phi_{N}}-x)$ can be locally
approximated by $\mathcal{P}_{\phi(x)}$, and thus
$N^{1/d}(\mathcal{P}_{N\phi_{N}}^{x+N^{-1/d}z}-x)$ can be approximated by
$\mathcal{P}_{\phi(x)}\cup\\{z\\}$. Since the graph is translation invariant
and stabilizing on $\mathcal{P}_{\phi(x)}^{\mathbf{0}}$, it is also
stabilizing on $\mathcal{P}_{\phi(x)}^{z}$ [56, Lemma 3.3]. Therefore,

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\xi(x+N^{-1/d}z,\mathcal{P}_{N\phi_{N}})$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\Bigg{[}\frac{1}{d_{x+N^{-1/d}z}}\sum_{y:(x+N^{-1/d}z,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x+N^{-1/d}z})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x+N^{-1/d}z)}{\phi(x+N^{-1/d}z)}\frac{\pi_{q}f_{q}(y)}{\phi(y)}$
$\displaystyle\quad-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x+N^{-1/d}z)}{\phi(x+N^{-1/d}z)}\Bigg{]}$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\left[\frac{1}{d_{x+N^{-1/d}z}}\sum_{y:(x+N^{-1/d}z,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x+N^{-1/d}z})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right]$
$\displaystyle=\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}.$

Further, as shown in the proof of Penrose et al. [56, Lemma 3.3], if a
geometric graph is translation invariant and stabilizing on
$\mathcal{P}_{1}^{\mathbf{0}}$, then it is also stabilizing on
$\mathcal{P}_{1}\cup\\{\mathbf{0},z\\}$. Hence,

$\displaystyle\quad\lim_{N\to\infty}\mathbb{E}\big{[}\xi(x,\mathcal{P}_{N\phi_{N}}^{x+N^{-1/d}z})\xi(x+N^{-1/d}z,\mathcal{P}_{N\phi_{N}}^{x})\big{]}$
$\displaystyle=\lim_{N\to\infty}\mathbb{E}\Bigg{[}\Big{(}\frac{1}{d_{x}}\sum_{y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x,x+N^{-1/d}z})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\Big{)}\times$
$\displaystyle\Big{(}\frac{1}{d_{x+N^{-1/d}z}}\sum_{y:(x+N^{-1/d}z,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}}^{x,x+N^{-1/d}z})}\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\Big{)}\Bigg{]}$
$\displaystyle=\left(\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2}.$

Hence the integrand in (C.17) converges pointwise to 0. Combining the above
convergence results we have:

${\rm
Var}[T]\to\int\left(\sum_{p,q=1}^{M}K(p,q)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\frac{\pi_{q}f_{q}(x)}{\phi(x)}-\sum_{p=1}^{M}g(p)\frac{\pi_{p}f_{p}(x)}{\phi(x)}\right)^{2}\phi(x){\rm
d}x=:\kappa_{2}^{2}.$

If $\kappa_{2}^{2}=0$, then $T-\mathbb{E}T\overset{d}{\to}N(0,\kappa_{2}^{2})$
is trivial. Suppose $\kappa_{2}^{2}>0$. We will use the following proposition.

######  Proposition C.1.

Suppose that $\xi(x,\mathcal{X})$ is a bounded function, defined for any
$x\in\mathbb{R}^{d}$ and any finite set $\mathcal{X}\subset\mathbb{R}^{d}$,
which is power-law stabilizing with respect to $\phi_{N}$ with order
$q_{0}>16d$ (see Definition C.2). Suppose $\phi(\cdot)$ is bounded and has
bounded support. Let
$T=\frac{1}{\sqrt{N}}\sum_{x\in\mathcal{P}_{N\phi_{N}}}\xi(x,\mathcal{P}_{N\phi_{N}})$.
Suppose ${\rm Var}(T)\to\sigma^{2}>0$. Then
$T-\mathbb{E}T\overset{d}{\to}N(0,\sigma^{2})$.

###### Proof.

We first cover ${\rm supp}(\phi)$ by $V=V(N)$ cubes of the form
$Q=\prod_{i=1}^{d}[j_{i}s_{N},(j_{i}+1)s_{N}]$, where
$s_{N}:=N^{-1/d}\rho_{N}$. Here $\rho_{N}$ will be taken as $N^{a}$ such that
$a<\frac{1}{8d}$ and $(aq_{0}-1)>1$. Since $\phi$ has bounded support,
$V(N)=O\left(\frac{1}{s_{N}^{d}}\right)=O\left(N\rho_{N}^{-d}\right)$. The
points in $\mathcal{P}_{N\phi_{N}}$ can be labelled as:

$X_{i,j},\ 1\leq i\leq V,\ 1\leq j\leq M_{i},\ M_{i}\sim{\rm
Poisson}\left(N\int_{Q_{i}}\phi_{N}(x){\rm d}x\right)$

where $Q_{1},\ldots,Q_{V}$ is an enumeration of the cubes that cover ${\rm
supp}(\phi)$ as described above. As $\phi$ is bounded, each $f_{i}\
(i=1,\ldots,M)$ is also bounded, and consequently $\|\phi_{N}\|_{\infty}$ is
bounded by a constant for all $N$, and $N\int_{Q_{i}}\phi_{N}(x){\rm d}x\leq
N\|\phi_{N}\|_{\infty}\int_{Q_{i}}{\rm
d}x=Ns_{N}^{d}\|\phi_{N}\|_{\infty}=\rho_{N}^{d}\|\phi_{N}\|_{\infty}$. If we
write $\xi_{i,j}:=\xi(X_{i,j},\mathcal{P}_{N\phi_{N}})$, then we have the
following bound on the $q$-th moment (i.e., for any random variable $X$,
$\|X\|_{q}:=(\mathbb{E}|X|^{q})^{1/q}$):

$\left\|\sum_{j=1}^{M_{i}}|\xi_{i,j}|\right\|_{q}\leq\|\xi\|_{\infty}\|M_{i}\|_{q}\leq\|\xi\|_{\infty}\left\|{\rm
Poisson}\left(\rho_{N}^{d}\|\phi_{N}\|_{\infty}\right)\right\|_{q}.$

Using the fact that $\|{\rm
Poisson}(\lambda)\|_{q}\leq\lambda\exp(q/(2\lambda))$ (see [1]), and
$\|\phi_{N}\|_{\infty}$ being bounded, we have:

$\left\|\sum_{j=1}^{M_{i}}|\xi_{i,j}|\right\|_{q}\lesssim\rho_{N}^{d}.$ (C.19)

Let $R_{i,j}$ denote the radius of stabilization of $\xi$ at $X_{i,j}$,
$E_{i,j}:=\\{R_{i,j}\leq\rho_{N}\\}$,
$E:=\cap_{i=1}^{V}\cap_{j=1}^{M_{i}}E_{i.j}$. Then by Palm theory [53, Theorem
1.6],

$\mathbb{P}[E^{c}]\leq\mathbb{E}\left[\sum_{i=1}^{V}\sum_{j=1}^{M_{i}}1_{E_{i,j}^{c}}\right]=N\int\mathbb{P}[R(x,N)>\rho_{N}]\phi_{N}(x){\rm
d}x\leq N\tau(\rho_{N}).$

Since $\xi$ is power law stabilizing of order $q_{0}$,
$\tau(t)\leq\frac{C}{t^{q_{0}}}$ for some $C>0$. Since $\rho_{N}=N^{a}$, we
have $\mathbb{P}[E^{c}]\leq\frac{C}{N^{aq_{0}-1}}$. Consider

$T^{\prime}:=\frac{1}{\sqrt{N}}\sum_{i=1}^{V(N)}\sum_{j=1}^{M_{i}}\xi(X_{i,j},\mathcal{P}_{N\phi_{N}})\cdot
1_{E_{i,j}}$

which equals $T$ on $E$, and is a sum of $V(N)$ “near independent” random
variables. We show that ${\rm Var}[T^{\prime}]$ is close to ${\rm Var}[T]$. By
Hölder’s inequality, for $q>2$:

$\displaystyle\|T-T^{\prime}\|_{2}$
$\displaystyle\leq\|T-T^{\prime}\|_{q}\mathbb{P}[E^{c}]^{1-\frac{2}{q}}$
$\displaystyle\leq\left(\|T\|_{q}+\|T^{\prime}\|_{q}\right)\left(\frac{C}{N^{aq_{0}-1}}\right)^{1-\frac{2}{q}}.$

To bound $\|T\|_{q}$ and $\|T^{\prime}\|_{q}$, we use (C.19) and
$V(N)=O\left(N\rho_{N}^{-d}\right)$ to get, for $q\geq 2$,

$\|T^{\prime}\|_{q}=\frac{1}{\sqrt{N}}\left\|\sum_{i=1}^{V(N)}\sum_{j=1}^{M_{i}}\xi(X_{i,j},\mathcal{P}_{N\phi_{N}})\cdot
1_{E_{i,j}}\right\|_{q}\lesssim\frac{1}{\sqrt{N}}\cdot
N\rho_{N}^{-d}\cdot\rho_{N}^{d}=\sqrt{N}.$

Similarly $\|T\|_{q}\lesssim\sqrt{N}$. Hence

$\|T-T^{\prime}\|_{2}\lesssim\sqrt{N}\cdot\left(\frac{C}{N^{aq_{0}-1}}\right)^{1-\frac{2}{q}}.$

Hence,

$\displaystyle|{\rm Var}(T)-{\rm Var}(T^{\prime})|$ $\displaystyle=|{\rm
Var}(T-T^{\prime})+2{\rm Cov}(T-T^{\prime},T^{\prime})|$
$\displaystyle\leq\|T-T^{\prime}\|_{2}^{2}+2\|T-T^{\prime}\|_{2}\|T^{\prime}\|_{2}$
$\displaystyle\lesssim
N\left(\frac{C}{N^{ap-1}}\right)^{2-\frac{4}{q}}+2\cdot\sqrt{N}\left(\frac{C}{N^{ap-1}}\right)^{1-\frac{2}{q}}\cdot\sqrt{N}.$

We can choose $q$ large enough such that $(aq_{0}-1)(1-2/q)-1>0$. Then we have
$|{\rm Var}(T)-{\rm Var}(T^{\prime})|\to 0$. In particular, ${\rm
Var}(T^{\prime})\to\sigma^{2}>0$. Also,
$|\mathbb{E}T-\mathbb{E}T^{\prime}|\leq\mathbb{E}|T-T^{\prime}|\leq\|T-T^{\prime}\|_{2}\to
0$.

Write $\left(T^{\prime}-\mathbb{E}[T^{\prime}]\right)/\sqrt{{\rm
Var}(T^{\prime})}=\sum_{i=1}^{V}S_{i}$, where

$S_{i}=\frac{1}{\sqrt{N{\rm
Var}(T^{\prime})}}\left(\sum_{j=1}^{M_{i}}\xi(X_{i,j},\mathcal{P}_{N\phi_{N}})\cdot
1_{E_{i,j}}-\mathbb{E}\left[\sum_{j=1}^{M_{i}}\xi(X_{i,j},\mathcal{P}_{N\phi_{N}})\cdot
1_{E_{i,j}}\right]\right).$

Construct a dependency graph on $\\{S_{1},\ldots,S_{V}\\}$ as follows: for
$i\neq j$, there is an edge between $S_{i}$ and $S_{j}$ if $d(Q_{i},Q_{j})\leq
2N^{-1/d}\rho_{N}$, where $d(Q_{i},Q_{j}):=\inf\\{|x-y|:x\in Q_{i},y\in
Q_{j}\\}$. By definition of the radius of stabilization $R(x,N)$, the value of
$S_{i}$ is determined by the restriction of $\mathcal{P}_{N\phi_{N}}$ to the
$N^{-1/d}\rho_{N}$-neighborhood of the cube $Q_{i}$, i.e.,
$\\{x\in\mathbb{R}^{d}:|x-y|\leq N^{-1/d}\rho_{N}\textrm{ for some }y\in
Q_{i}\\}$. By the independence property of the Poisson process, for any pair
of disjoint sets $\Gamma_{1},\Gamma_{2}\subset\\{S_{1},\ldots,S_{V}\\}$ such
that no edge has one endpoint in $\Gamma_{1}$ and the other in $\Gamma_{2}$,
$\Gamma_{1}$ is independent of $\Gamma_{2}$. Since the number of cubes in
$Q_{1},\ldots,Q_{V}$ that are at most $2N^{-1/d}\rho_{N}$ distant from a given
cube is bounded by $5^{d}$, it follows that the maximal degree of the
dependency graph is bounded by $D\leq 5^{d}$, a constant.

By (C.19) and ${\rm Var}(T^{\prime})\to\sigma^{2}>0$,

$\mathbb{E}[|S_{i}|^{3}]\lesssim\frac{1}{N^{3/2}}\rho_{N}^{3d}=\frac{1}{N^{3(0.5-ad)}}.$

Hence we can take $p=3$ and $\theta=\frac{C^{\prime}}{N^{0.5-ad}}$ for some
constant $C^{\prime}$ in Theorem C.3 to yield

$\displaystyle\sup_{z\in\mathbb{R}}\left|\mathbb{P}\left(\frac{T^{\prime}-\mathbb{E}[T^{\prime}]}{\sqrt{{\rm
Var}(T^{\prime})}}\leq z\right)-\Phi(z)\right|$ $\displaystyle\leq
75D^{5(p-1)}|\mathcal{V}|\theta^{p}$ $\displaystyle\lesssim
N\rho_{N}^{-d}\cdot\frac{1}{N^{3(0.5-ad)}}=\frac{1}{N^{0.5-4ad}}\to 0,$

as $N\to\infty$, which implies
$\frac{T^{\prime}-\mathbb{E}[T^{\prime}]}{\sqrt{{\rm
Var}(T^{\prime})}}\overset{d}{\to}N(0,1)$. Since $|{\rm Var}(T)-{\rm
Var}(T^{\prime})|\to 0$, $|\mathbb{E}T-\mathbb{E}T^{\prime}|\to 0$, and
$T,T^{\prime}$ coincide on the set $E$ whose probability converges to 1, we
have $\frac{T-\mathbb{E}[T]}{\sqrt{{\rm Var}(T)}}\overset{d}{\to}N(0,1)$. ∎

###  C.13 Proof of Theorem B.2

Note that Lemmas C.1 and C.3 still hold if $\mathcal{P}_{N\phi_{N}}$ is
replaced by $\mathcal{P}_{N\phi_{N}^{N}}$ and $\phi(x_{0})$ is replaced by
$f(x_{0})$, assuming $f_{1}^{N},\ldots,f_{M}^{N}$ are equicontinuous at
$x_{0}$ and $f_{i}^{N}(x_{0})\to f(x_{0})$ as $N\to\infty$. The slight change
in the proof is that (C.9) now becomes:

$N\int_{B(x_{0},N^{-1/d}K)}|\phi_{N}^{N}(x)-f(x_{0})|{\rm d}x\leq
N\int_{B(x_{0},N^{-1/d}K)}(|\phi_{N}^{N}(x)-\phi_{N}^{N}(x_{0})|+|\phi_{N}^{N}(x_{0})-f(x_{0})|){\rm
d}x\to 0,$

since $\phi_{N}^{N}$ is equicontinuous at $x_{0}$ and $\phi_{N}^{N}(x_{0})\to
f(x_{0})$. The rest of the proof follows almost verbatim from the proof of
Theorem B.1, by replacing every $f_{i}$ by $f_{i}^{N}$, every $\phi_{N}$ by
$\phi_{N}^{N}$, every $\phi$ by $\phi^{N}:=\sum_{i=1}^{M}\pi_{i}f_{i}^{N}$,
and applying an extra limit (as $N\to\infty$):

$\left\\{\begin{aligned} &\int|\phi^{N}(z)-f(z)|{\rm d}z\to 0,\\\
&\frac{\pi_{i}f_{i}^{N}(z)}{\phi^{N}(z)}\to\pi_{i},\qquad\mbox{for
a.e.}\;z\in{\rm supp}(f),\qquad\mbox{for }\;i=1,\ldots,M.\end{aligned}\right.$

The limiting variance $\kappa_{1}^{2}$ reduces to the numerator of the
asymptotic null variance in (8), which is distribution-free. The limiting
variance

$\kappa_{2}^{2}=\int\left(\sum_{p,q=1}^{M}K(p,q)\pi_{p}\pi_{q}-\sum_{p=1}^{M}g(p)\pi_{p}\right)^{2}f(x){\rm
d}x=0.$

Recall that
$g(\Delta)=2\sum_{i=1}^{M}\pi_{i}K(\Delta,i)-\sum_{i,j=1}^{M}\pi_{i}\pi_{j}K(i,j)$
was defined in Theorem B.1. Hence,
$\mathbb{E}(H_{n}|\mathcal{F}_{n})-\mathbb{E}(H_{n})\overset{d}{\to}N(0,\kappa_{2}^{2})\equiv\delta_{0}$
as the variance of $\mathbb{E}(H_{n}|\mathcal{F}_{n})$ converges to 0. Note
that we need $q_{0}>d$ here to establish the convergence of ${\rm
Var}\left(\mathbb{E}[H_{n}|\mathcal{F}_{n}]\right)$; see (C.18). Recall that
the power-law stabilization with order $q_{0}>16d$ was used for proving
Proposition C.1 which establishes the CLT assuming the variance converges to a
non-zero quantity, and is not needed in this proof.

###  C.14 Proof of Theorem B.3

Since

$\frac{n\tilde{H}_{n}-\mathbb{E}[n\tilde{H}_{n}]}{\sqrt{n}\kappa_{1,{\rm
null}}}\overset{d}{\to}N(0,1),$

the local power of the test in (B.7) is

$\Phi\left(z_{\alpha}+\lim\frac{\mathbb{E}[n\tilde{H}_{n}]}{\kappa_{1,{\rm
null}}\sqrt{n}}\right).$

As $\frac{n}{N}\overset{p}{\to}1$, we know that
$\frac{\mathbb{E}[n\tilde{H}_{n}]}{\sqrt{n}}$ has the same limit in
probability as the expectation of

$\frac{1}{\sqrt{N}}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})-\frac{1}{\sqrt{N}}\sum_{i=1}^{n}g_{N}(\Delta_{i}).$

Note that the expectation of the second term
$\frac{1}{\sqrt{N}}\sum_{i=1}^{n}g_{N}(\Delta_{i})$ equals the expectation of
the first term
$\frac{1}{\sqrt{N}}\sum\limits_{i=1}^{n}\frac{1}{d_{i}}\sum\limits_{j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}K(\Delta_{i},\Delta_{j})$
under the null, i.e., $\theta_{1}=\theta_{2}$. Hence we study the first term
as follows. With a $k$-NN graph and the discrete kernel, the expectation of
the first term can be written as

$\displaystyle\quad\frac{1}{k\sqrt{N}}\mathbb{E}\left[\sum_{i,j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}\left(1-I(\Delta_{i}\neq\Delta_{j})\right)\right]$
$\displaystyle=\sqrt{N}-\frac{1}{k\sqrt{N}}\mathbb{E}\left[\sum_{i,j:(Z_{i},Z_{j})\in\mathcal{E}(\mathcal{G}_{n})}I(\Delta_{i}\neq\Delta_{j})\right]$
$\displaystyle=\sqrt{N}-\frac{1}{k\sqrt{N}}\mathbb{E}\left[\sum_{x,y:(x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}^{N}})}\left(\frac{\frac{N_{1}}{N}f(x|\theta_{1})}{\phi_{N}^{N}(x)}\frac{\frac{N_{2}}{N}f(y|\theta_{2})}{\phi_{N}^{N}(y)}+\frac{\frac{N_{2}}{N}f(x|\theta_{2})}{\phi_{N}^{N}(x)}\frac{\frac{N_{1}}{N}f(y|\theta_{1})}{\phi_{N}^{N}(y)}\right)\right]$
$\displaystyle=\sqrt{N}-\frac{1}{k\sqrt{N}}N_{1}N_{2}\int_{S\times
S}(f(x|\theta_{1})f(y|\theta_{2})+f(x|\theta_{2})f(y|\theta_{1}))\rho_{K}^{\theta_{1},\theta_{2}}(x,y){\rm
d}x{\rm d}y,$

where
$\phi_{N}^{N}(x)=\frac{N_{1}}{N}f(x|\theta_{1})+\frac{N_{2}}{N}f(x|\theta_{2})$,
and
$\rho_{K}^{\theta_{1},\theta_{2}}(x,y)=\mathbb{P}\left((x,y)\in\mathcal{E}(\mathcal{P}_{N\phi_{N}^{N}}^{x,y})\right)$.
The last equality above follows from Palm theory. Define

$\mu_{N}^{(S)}(\theta_{1},\theta_{2}):=N^{2}\int_{S\times
S}(f(x|\theta_{1})f(y|\theta_{2})+f(x|\theta_{2})f(y|\theta_{1}))\rho_{K}^{\theta_{1},\theta_{2}}(x,y){\rm
d}x{\rm d}y.$

Then,

$\displaystyle\lim_{N\to\infty}\frac{\mathbb{E}[n\tilde{H}_{n}]}{\sqrt{n}}$
$\displaystyle=\lim_{N\to\infty}-\frac{N_{1}N_{2}}{kN^{2}}\frac{\mu_{N}^{(S)}(\theta_{1},\theta_{2})-\mu_{N}^{(S)}(\theta_{1},\theta_{1})}{\sqrt{N}}$
$\displaystyle=\lim_{N\to\infty}-\frac{N_{1}N_{2}}{kN^{2}}\left(\frac{\varepsilon_{N}^{\top}\nabla\mu_{N}^{(S)}(\theta_{1},\theta_{1})}{\sqrt{N}}+\frac{1}{2}\frac{\varepsilon_{N}^{\top}{\rm
H}\mu_{N}^{(S)}(\theta_{1},\theta_{1})\varepsilon_{N}}{\sqrt{N}}+\mathcal{R}_{N}\right),$

where the gradient and Hessian are taken with respect to $\theta_{2}$. The
limit of the gradient and Hessian follows the argument in Appendix E in the
supplementary file of Bhattacharya [11]. Note that there is a typo in the
limit of the Hessian term in Bhattacharya [11, Appendix E]. In fact, with
$\varepsilon_{N}=hN^{-1/4}$, we can show that

$\frac{\varepsilon_{N}^{\top}{\rm
H}\mu_{N}^{(S)}(\theta_{1},\theta_{1})\varepsilon_{N}}{\sqrt{N}}\to-2rk\mathbb{E}\left[\frac{h^{\top}\nabla_{\theta_{1}}f(X|\theta_{1})}{f(X|\theta_{1})}\right]^{2},$

where $r=2\pi_{1}\pi_{2}$. The coefficient above is 2 instead of $3/2$ in the
original version of Bhattacharya [11]. The rest of the proof follows from the
arguments at the beginning of Bhattacharya [11, Appendix B]. ∎

###  C.15 Proof of Theorem C.1

The proof is similar to that in Pham et al. [58]. But for completeness, we
provide the entire proof here. The proof applies the method of moments, an
idea that dates back to Pafnutii Lvovich Chebyshev (1821–1894) [27]. More
specifically, we will show that for all $r\in\mathbb{N}$,

$\mathbb{E}\left[\left(\frac{\sum a_{ij}B_{ij}}{{\rm Var}(\sum
a_{ij}B_{ij})}\right)^{r}\right]\to\mathbb{E}\left[Z^{r}\right],$

where $Z\sim N(0,1)$. This implies that for any subsequence that converges in
distribution, the limiting distribution has the same moments as $N(0,1)$, and
therefore must be $N(0,1)$ as these moments uniquely determine the
distribution $N(0,1)$ (see e.g., [12, Theorem 30.1]). Since weak convergence
to $N(0,1)$ holds for any subsequence, it further implies that the entire
sequence converges in distribution as weak convergence is metrizable.

In this proof, we will decompose $\left(\sum a_{ij}B_{ij}\right)^{r}$ into
“sums corresponding to different graphs” as in Bloemena [14]. The sums
corresponding to most of the graphs will be negligible when compared to
$\left[{\rm Var}(\sum a_{ij}B_{ij})\right]^{r}$. The remaining dominating
terms will lead to the moments of $N(0,1)$.

We first define equivalent graphs. In this proof, a graph will be considered
as a collection of edges (so it does not contain isolated vertices), and we
will consider the edges in a graph to be different, each having a label (from
$1,\ldots,r$) and a direction. Two graphs are equivalent if they can be mapped
to each other while keeping the direction and the labeling of the edges.
Multiple edges are allowed to exist between 2 vertices. For a graph $G$,
denote by $|G|$ the number of its vertices, and arbitrarily label the vertices
of $G$ as $1,\ldots,|G|$. Then an edge in $G$ can be given by $(\mu,\nu)$,
where $\mu,\nu\in\\{1,\ldots,|G|\\}$ are integers. Define

$\Sigma(a,G):={{\sum}^{\prime}_{i_{1},\ldots,i_{|G|}}}\prod_{(\mu,\nu)\in
G}a_{i_{\mu}i_{\nu}},$ (C.20)

where $\sum^{\prime}$ means the summation is over distinct
$i_{1},\ldots,i_{|G|}\in\\{1,\ldots,n\\}$.

Denote by $G(r)$ the set of all graphs with $r$ edges (equivalent graphs will
only be counted once). Then

$(\sum a_{ij})^{r}=\sum_{G\in G(r)}\Sigma(a,G).$ (C.21)

It can be seen from (C.20) that if $G$ has a self-loop, then $\Sigma(a,G)=0$
as $a_{ii}=b_{ii}=0$, $\forall i$. Hence, we can ignore graphs in $G(r)$ that
contain self-loops. To understand (C.21), for example, when $r=2$,

$\displaystyle\left(\sum a_{ij}\right)^{2}$
$\displaystyle={\sum}^{\prime}a_{ij}a_{kl}$
$\displaystyle+{\sum}^{\prime}a_{ij}a_{ik}+{\sum}^{\prime}a_{ij}a_{ki}+{\sum}^{\prime}a_{ij}a_{jk}+{\sum}^{\prime}a_{ij}a_{kj}$
$\displaystyle+{\sum}^{\prime}a_{ij}^{2}+{\sum}^{\prime}a_{ij}a_{ji}.$

The sum in the first line corresponds to a graph with two isolated edges. The
sums in the second line correspond to graphs with two edges sharing exactly
one vertex, but the edges can have different orientations. The sums in the
third line correspond to graphs with two edges sharing two vertices (which
forms a cycle), with two possible orientations.

Now,

$\displaystyle\mathbb{E}\left(\sum a_{ij}B_{ij}\right)^{r}$
$\displaystyle=\sum_{G\in
G(r)}\operatorname*{{\sum}^{\prime}}_{i_{1},\ldots,i_{|G|}}\mathbb{E}\prod_{(\mu,\nu)\in
G}a_{i_{\mu}i_{\nu}}B_{i_{\mu}i_{\nu}}$ (C.22) $\displaystyle=\sum_{G\in
G(r)}\mathbb{E}\left[\prod_{(\mu,\nu)\in
G}B_{i_{\mu}i_{\nu}}\right]\Sigma(a,G)$ $\displaystyle=\sum_{G\in
G(r)}\frac{1}{n(n-1)\cdots(n-|G|+1)}\Sigma(b,G)\Sigma(a,G).$

The second equality follows from the fact that each term $\prod_{(\mu,\nu)\in
G}B_{i_{\mu}i_{\nu}}$ has the same expectation, which is
$\frac{1}{n(n-1)\cdots(n-|G|+1)}\Sigma(b,G)$ as there are
$n(n-1)\cdots(n-|G|+1)$ many terms in $\Sigma(b,G)$.

We first estimate the order of $\Sigma(b,G)$. Suppose $G$ has an isolated
edge, say $(1,2)$. By definition,
$\Sigma(b,G)={{\sum}^{\prime}_{i_{1},\ldots,i_{|G|}}}\prod_{(\mu,\nu)\in
G}b_{i_{\mu}i_{\nu}}$ with the summation indices $i_{1},\ldots,i_{|G|}$
required to be distinct. If we relax the constraint for $i_{1},i_{2}$, i.e.,
$i_{1},i_{2}$ can freely take values from $1,\ldots,n$ while
$i_{3},\ldots,i_{|G|}$ are still required to be distinct, then the sum will be
0 as $\sum b_{ij}=0$. But then we need to subtract back the terms which were
originally not in $\Sigma(b,G)$ and came in because of relaxing the
constraints on $i_{1},i_{2}$. These terms are sum of a number of
$\Sigma(b,G^{\prime})$, where $G^{\prime}$ is the graph obtained from $G$ by
identifying vertex 1 or 2 or both with some vertex in $3,\ldots,|G|$. In
general, if $G$ has multiple isolated edges, $\Sigma(b,G)$ can be written as a
linear combination of $\Sigma(b,G^{\prime})$, where $G^{\prime}$ can be
obtained by sequentially identifying the vertices of isolated edges with other
vertices in the graph, and $G^{\prime}$ no longer has any isolated edge. Note
that

$|\Sigma(b,G^{\prime})|=\left|\operatorname*{{\sum}^{\prime}}_{i_{1},\ldots,i_{|G^{\prime}|}}\prod_{(\mu,\nu)\in
G^{\prime}}b_{i_{\mu}i_{\nu}}\right|\leq\prod_{(\mu,\nu)\in
G^{\prime}}\left(\sum_{i_{1},\ldots,i_{|G^{\prime}|}}|b_{i_{\mu}i_{\nu}}|^{r}\right)^{\frac{1}{r}},$

where we have used the fact that for $C_{ij}\geq 0$, $1\leq i\leq I$, $1\leq
j\leq J$,

$\prod_{i=1}^{I}\left(\sum_{j=1}^{J}C_{ij}^{r}\right)\geq\left(\sum_{j=1}^{J}\prod_{i=1}^{I}C_{ij}\right)^{r}.$

Since
$\sum_{i_{1},\ldots,i_{|G^{\prime}|}}|b_{i_{\mu}i_{\nu}}|^{r}=n^{|G^{\prime}|-2}\sum|b_{ij}|^{r}$
and by (B2), it is $O\left(n^{|G^{\prime}|-r}\|b\|^{r}\right)$, where $\|b\|$
is defined as $(\sum b_{ij}^{2})^{1/2}$, we have
$|\Sigma(b,G^{\prime})|=O\left(n^{|G^{\prime}|-r}\|b\|^{r}\right)$.

Let $s$ be the number of isolated edges in $G$. Note that if $(1,2)$ is an
isolated edge, identifying $1$ or $2$ with some vertex in $3,\ldots,|G|$
reduces the number of vertices by 1, and the number of isolated edges by at
most 2. Identifying both $1$ and $2$ with some vertices in $3,\ldots,|G|$
reduces the number of vertices by 2, and the number of isolated edges by at
most 3. Hence in either way, when the number of vertices is reduced by 1, the
number of isolated edges is reduced by at most 2. Since $G^{\prime}$ has no
isolated edge, $|G^{\prime}|\leq|G|-\lceil\frac{s}{2}\rceil$. Combining the
discussions above:

$\Sigma(b,G)=O\left(n^{|G|-\lceil\frac{s}{2}\rceil-r}\cdot\|b\|^{r}\right).$
(C.23)

Next, we estimate the order of $\Sigma(a,G)$. Without loss of generality, we
can suppose $\max|a_{ij}|=1$. Then, if $G$ is connected, by suppressing
redundant edges in $G$ to a spanning tree $T$, we have:

$\Sigma(|a|,G):=\operatorname*{{\sum}^{\prime}}_{i_{1},\ldots,i_{|G|}}\prod_{(\mu,\nu)\in
G}|a_{i_{\mu}i_{\nu}}|\leq\operatorname*{{\sum}^{\prime}}_{i_{1},\ldots,i_{|G|}}\prod_{(\mu,\nu)\in
T}|a_{i_{\mu}i_{\nu}}|.$

Here $|a|$ means taking the absolute value of each $a_{ij}$. Recall assumption
(A1): $\max_{i}\sum_{j}|a_{ij}|=O(\max|a_{ij}|\cdot t_{n})$. Hence by
sequentially summing over the index corresponding to a leaf,

$\Sigma(|a|,G)\leq\operatorname*{{\sum}^{\prime}}_{i_{1},\ldots,i_{|G|}}\prod_{(\mu,\nu)\in
T}|a_{i_{\mu}i_{\nu}}|\leq nt_{n}^{|G|-1}.$

If $G$ has $c$ connected components, we can relax the summation constraint
“$i_{1},\ldots,i_{|G|}$ are distinct” to “the indices within each component
are distinct”, which implies

$\Sigma(|a|,G)\leq n^{c}t_{n}^{|G|-c}\leq n^{c}t_{n}^{|G|}.$

If among these $c$ connected components, $s$ of them are isolated edges, then
similar to the previous argument, using $\sum a_{ij}=0$, $\Sigma(a,G)$ can be
written as a linear combination of $\Sigma(a,G^{\prime})$, where $G^{\prime}$
can be obtained by sequentially identifying the vertices of isolated edges
with other vertices in the graph, and $G^{\prime}$ no longer has any isolated
edge. Note that if $(1,2)$ is an isolated edge, identifying $1$ or $2$ with
some vertex in $3,\ldots,|G|$ reduces the number of connected components by 1
and the number of isolated edges by 1 or 2. Identifying both $1$ and $2$ with
some vertices in $3,\ldots,|G|$ may reduce the number of isolated edges by 3,
but in such a case the number of connected components is reduced by 2. Hence
in either case, when the number of connected components is reduced by 1, the
number of isolated edges is reduced by at most 2. Because $G^{\prime}$ no
longer has any isolated edge, it has at most $c-\lceil\frac{s}{2}\rceil$
connected components. Therefore,

$\Sigma(a,G)=O\left(n^{c-\lceil\frac{s}{2}\rceil}t_{n}^{|G|}\right).$

Recall assumption (A2): $\liminf\sum a_{ij}^{2}/(n\max a_{ij}^{2})>0$, which
implies

$\Sigma(a,G)=O\left(n^{c-\lceil\frac{s}{2}\rceil-\frac{r}{2}}t_{n}^{|G|}\cdot\|a\|^{r}\right).$
(C.24)

By considering the number of edges: $2(c-s)+s\leq r\Leftrightarrow
c\leq\frac{r}{2}+\frac{s}{2}$, we know

$\Sigma(a,G)\leq n^{-1/2}t_{n}^{|G|}\cdot\|a\|^{r}=o\left(\|a\|^{r}\right),$

unless every connected component has at most 2 edges and
$\lceil\frac{s}{2}\rceil=\frac{s}{2}$. Together with
$\Sigma(b,G)=O\left(n^{|G|-\lceil\frac{s}{2}\rceil-r}\cdot\|b\|^{r}\right)$
deducted previously, we have

$\frac{1}{n(n-1)\cdots(n-|G|+1)}\Sigma(b,G)\Sigma(a,G)=o\left(\frac{\|a\|^{r}\|b\|^{r}}{n^{r}}\right),$

unless $s=0$ and every connected component has at most 2 edges, which is then
equivalent to all connected components having exactly 2 edges.

We now claim that $w_{n}\gtrsim\frac{\|a\|\|b\|}{n}$. Hence in the expansion
(C.22) of $\mathbb{E}\left(\sum a_{ij}B_{ij}\right)^{r}$, the sum
corresponding to these graphs are negligible compared to $w_{n}^{r}$. To see
the claim, note that $\sum a_{ij}=0$ implies $\sum^{\prime}a_{ij}a_{ik}=\sum
a_{i+}^{2}-\sum a_{ij}^{2}$, and similarly $\sum^{\prime}b_{ij}b_{ik}=\sum
b_{i+}^{2}-\sum b_{ij}^{2}$. Hence,

$\displaystyle w_{n}^{2}$ $\displaystyle=\frac{4(n-2)}{n^{4}}\left(\sum
a_{i+}^{2}\right)\left(\sum b_{i+}^{2}\right)$ (C.25)
$\displaystyle\quad+\frac{2}{n^{2}}\left(\sum a_{ij}^{2}-\frac{2}{n}\sum
a_{i+}^{2}\right)\left(\sum b_{ij}^{2}-\frac{2}{n}\sum
b_{i+}^{2}\right)+\frac{4}{n^{3}}\left(\sum a_{ij}^{2}\right)\left(\sum
b_{ij}^{2}\right).$

From (A1), (A2), and (B1), we know that the second term in the right-hand side
of (C.25) is $\gtrsim\frac{\|a\|^{2}\|b\|^{2}}{n^{2}}$. Hence the third term
in (C.25) is negligible compared to the second term. Moreover, since ${\rm
Var}\left(\sum a_{ij}B_{ij}\right)$ has the following expression [58, Equation
(2.6)]:

$\displaystyle{\rm Var}\left(\sum a_{ij}B_{ij}\right)$
$\displaystyle=\frac{4}{(n-1)(n-2)^{2}}\left(\sum a_{i+}^{2}\right)\left(\sum
b_{i+}^{2}\right)$ $\displaystyle\quad+\frac{2}{n(n-3)}\left(\sum
a_{ij}^{2}-\frac{2n-5}{(n-2)^{2}}\sum a_{i+}^{2}\right)\left(\sum
b_{ij}^{2}-\frac{2n-5}{(n-2)^{2}}\sum b_{i+}^{2}\right),$

we have $w_{n}^{2}/{\rm Var}(\sum a_{ij}B_{ij})\to 1$.

Now consider the remaining graphs $G$ consisting of connected components
$G_{1},\ldots,G_{p}$, each of which has exactly 2 edges. Hence necessarily,
$r=2p$ is even in order that $\mathbb{E}\left[\left(\frac{\sum
a_{ij}B_{ij}}{{\rm Var}(\sum a_{ij}B_{ij})}\right)^{r}\right]$ has a non-zero
limit. Each $G_{q}$ has two possible configurations141414the configuration of
a graph is the blank graph obtained by disregarding the direction of all
edges, either having 2 vertices and 2 edges forming a loop, or having 3
vertices and 2 edges forming a tree. For the loop configuration,
$\Sigma(a,G_{q})=\sum^{\prime}a_{ij}^{2}=\sum a_{ij}^{2}$. For the tree
configuration, $\Sigma(a,G_{q})=\sum^{\prime}a_{ij}a_{ik}=\sum a_{ij}a_{ik}$.
Since $\Sigma(a,G)=\sum_{i_{1},\ldots,i_{|G|}}^{\prime}\prod_{(\mu,\nu)\in
G}a_{i_{\mu}i_{\nu}}$, if we relax the constraint that $i_{1},\ldots,i_{|G|}$
are distinct, then we get $\Sigma(a,G_{1}\cup\cdots\cup
G_{p})=\prod_{q=1}^{p}\Sigma(a,G_{q})+R_{1}$. The difference $R_{1}$ is the
sum of a number of $\Sigma(a,G^{\prime})$, where $G^{\prime}$ is obtained from
identifying vertices of $G_{1}\cup\cdots\cup G_{p}$. Since $G^{\prime}$ has at
most $p-1$ connected components, by (C.24),
$\Sigma(a,G^{\prime})=O\left(n^{p-1-\frac{r}{2}}t^{|G|}_{n}\cdot\|a\|^{r}\right)=O\left(\frac{t^{|G|}_{n}}{n}\cdot\|a\|^{r}\right)$.
Similarly, we can relax the constraint that $i_{1},\ldots,i_{|G|}$ are
distinct in $\Sigma(b,G_{1}\cup\cdots\cup G_{p})$, which leads to
$\Sigma(b,G)=\prod_{q=1}^{p}\Sigma(b,G_{q})+R_{2}$ and the difference $R_{2}$
is a sum of a number of $\Sigma(b,G^{\prime})$. Since $|G^{\prime}|\leq|G|-1$,
by (C.23), $\Sigma(b,G^{\prime})=O\left(n^{|G|-r-1}\cdot\|b\|^{r}\right)$.
Note that (C.23) and (C.24) also give
$\Sigma(a,G)=O\left(t^{|G|}_{n}\cdot\|a\|^{r}\right)$,
$\Sigma(b,G)=O\left(n^{|G|-r}\cdot\|b\|^{r}\right)$. Hence together with
$w_{n}\gtrsim\frac{\|a\|\cdot\|b\|}{n}$,

$\displaystyle\frac{1}{n(n-1)\cdots(n-|G|+1)}\Sigma(b,G)\Sigma(a,G)$
$\displaystyle=\frac{(n-|G|)!}{n!}\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})+o(w_{n}^{r})$
$\displaystyle=\frac{1}{n^{|G|}}\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})+o(w_{n}^{r}).$

The second equality in the above display follows as
$\frac{(n-|G|)!}{n!}-\frac{1}{n^{|G|}}=O\left(\frac{1}{n^{|G|+1}}\right)$ and
$\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})=O\left(n^{|G|-r}t^{|G|}_{n}\cdot\|a\|^{r}\|b\|^{r}\right)$.
Therefore, (C.22) reduces to

$\displaystyle\mathbb{E}\left(\frac{\sum a_{ij}B_{ij}}{w_{n}}\right)^{r}$
$\displaystyle=\frac{1}{w_{n}^{r}}\sum_{G\in\tilde{G}(r)}\frac{1}{n^{|G|}}\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})+o(1),$

where $\tilde{G}(r)$ contains the graphs that consists of connected components
$G_{1},\ldots,G_{p}$, each of which has exactly 2 edges.

Suppose $i$ of $G_{1},\ldots,G_{p}$ have a loop configuration and $p-i$ of
them have a tree configuration. Then $|G|=2i+3(p-i)$ and

$\displaystyle\frac{1}{n^{|G|}}\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})$
$\displaystyle=\frac{1}{n^{3p-i}}\left(\sum a_{ij}^{2}\right)^{i}\left(\sum
b_{ij}^{2}\right)^{i}\left({\sum}^{\prime}a_{ij}a_{ik}\right)^{p-i}\left({\sum}^{\prime}b_{ij}b_{ik}\right)^{p-i}$
$\displaystyle=\frac{1}{n^{3p-i}}\xi^{i}\eta^{p-i},$

where $\xi:=\left(\sum a_{ij}^{2}\right)\left(\sum b_{ij}^{2}\right)$ and
$\eta:=\left({\sum}^{\prime}a_{ij}a_{ik}\right)\left({\sum}^{\prime}b_{ij}b_{ik}\right)$.

The number of graphs having $i$ loops and $p-i$ trees is given by:

$\frac{\tbinom{2p}{2}\tbinom{2p-2}{2}\cdots\tbinom{2}{2}}{i!(p-i)!}\cdot
2^{i}\cdot 4^{p-i}=\frac{(2p)!2^{p-i}}{i!(p-i)!},$

where we first divide the $2p$ different edges into $i$ loops and $p-i$ trees,
and then determine the orientation of each loop and tree. Therefore,

$\displaystyle\sum_{G\in\tilde{G}(r)}\frac{1}{n^{|G|}}\prod_{q=1}^{p}\Sigma(a,G_{q})\Sigma(b,G_{q})$
$\displaystyle=\sum_{i=0}^{p}\frac{(2p)!2^{p-i}}{i!(p-i)!}\cdot\frac{1}{n^{3p-i}}\xi^{i}\eta^{p-i}$
$\displaystyle=\frac{(2p)!(2\eta)^{p}}{n^{3p}p!}\sum_{i=0}^{p}\frac{p!}{i!(p-i)!}\left(\frac{\xi}{2\eta
n}\right)^{i}$
$\displaystyle=\frac{(2p)!(2\eta)^{p}}{n^{3p}p!}\left(1+\frac{\xi}{2\eta
n}\right)^{p}$
$\displaystyle=\frac{(2p)!}{p!}\left(\frac{2\eta}{n^{3}}+\frac{\xi}{n^{4}}\right)^{p}$
$\displaystyle=\frac{(2p)!}{2^{p}p!}\left(\frac{4\left({\sum}^{\prime}a_{ij}a_{ik}\right)\left({\sum}^{\prime}b_{ij}b_{ik}\right)}{n^{3}}+\frac{2\left(\sum
a_{ij}^{2}\right)\left(\sum b_{ij}^{2}\right)}{n^{4}}\right)^{p}$
$\displaystyle=\frac{(2p)!}{2^{p}p!}w_{n}^{2p}.$

Note that $\frac{(2p)!}{2^{p}p!}$ is the $2p$-th moment of $N(0,1)$. Hence the
proof is completed.∎

##  Appendix D Further Simulations

In this section we provide further simulation studies to support the major
results provided in the main paper.

###  D.1 Validity of Theorem 3

Figure 3 shows the histogram of $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}$ from 20000 independent replications
constructed using the directed 1-NN graph and the discrete kernel. Here we
take $M=3$ and the three distributions are $P_{1}=P_{2}=P_{3}\equiv
N(\mathbf{0},I_{d})$ with equal sample sizes $n_{i}$ and $d=2$. The red curve
is the standard normal density function. It can be seen from the plots that
the empirical distribution is already close to the standard normal
distribution for $n_{i}=100$. The approximation gets even better when
$n_{i}=1000$, and when $n_{i}=10000$, the empirical distribution is almost
identical to the standard normal.

Figure 3: Histograms of $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}$ as $n_{1}=n_{2}=n_{3}$ vary. The red
curves are the standard normal density.

###  D.2 Validation of Theorem B.3

Here we empirically check the validity of Theorem B.3 on the detection
threshold of our method.

Figure 4: Empirical power of our method with $n_{1}$ samples from
$P_{1}=N\left(\mathbf{0},3^{2}\cdot I_{d}\right)$ and $n_{2}$ samples from
$P_{2}=N\left(\mathbf{0},(3+2n^{b})^{2}\cdot I_{d}\right)$. $(n_{1},n_{2})$ is
taken as either $(2\cdot 10^{4},10^{4})$ or $(10^{4},2\cdot 10^{4})$,
corresponding to $(\pi_{1},\pi_{2})=(\frac{2}{3},\frac{1}{3})$ or
$(\frac{1}{3},\frac{2}{3})$. The red vertical line shows the detection
threshold predicted by Theorem B.3.

Although Theorem B.3 was proved in the Poissonized setting, our simulation
results suggest that the conclusions also hold true in the usual non-
Poissonized regime. We consider $M=2$ with $n_{1}$ samples from
$P_{1}=N\left(\mathbf{0},3^{2}\cdot I_{d}\right)$ and $n_{2}$ samples from
$P_{2}=N\left(\mathbf{0},(3+2n^{b})^{2}\cdot I_{d}\right)$ for $b\in(-1,0)$,
$n=n_{1}+n_{2}$. Here $(n_{1},n_{2})$ is taken as either $(2\cdot
10^{4},10^{4})$ or $(10^{4},2\cdot 10^{4})$, corresponding to
$(\pi_{1},\pi_{2})=(\frac{2}{3},\frac{1}{3})$ or $(\frac{1}{3},\frac{2}{3})$.
Figure 4 shows the empirical power of the asymptotic test (7) averaged over
200 replications. The level of the test is set as $\alpha=0.05$ and $k$-NN
graphs are used with $k=1,2,3$. The red vertical line shows the detection
threshold predicted by Theorem B.3. When $d=5$, the detection threshold is at
$n^{-\frac{1}{4}}$. It can be seen that for either choice of
$(\pi_{1},\pi_{2})$, the power increases from 0.05 to 1 around
$b=-\frac{1}{4}$. For $d=10$ and $d=20$, the detection threshold depends on
the sign of $a_{k,\theta_{1}}$ (see (B.8)). When
$(\pi_{1},\pi_{2})=(\frac{2}{3},\frac{1}{3})$, $a_{k,\theta_{1}}>0$, and the
detection threshold is $n^{-\frac{1}{2}+\frac{2}{d}}$, so the power increases
from 0.05 to 1 around $b=-\frac{1}{2}+\frac{2}{d}$ (see the middle and right
plots in the first row of Figure 4). When $d$ is large, this threshold gets
closer and closer to the parametric threshold $-\frac{1}{2}$ (see the top
right plot in Figure 4). When $(\pi_{1},\pi_{2})=(\frac{1}{3},\frac{2}{3})$,
$a_{k,\theta_{1}}<0$, and the detection threshold is $n^{-\frac{2}{d}}$, so
there is a rapid increase of the power from 0 to 1 around $b=-\frac{2}{d}$. In
such a case, the limiting power for
$b\in(-\frac{1}{2}+\frac{2}{d},-\frac{2}{d})$ is 0 as predicted by Theorem
B.3, and this is also supported by the last two plots in the second row of
Figure 4, where the empirical power for
$b\in(-\frac{1}{2}+\frac{2}{d},-\frac{2}{d})$ is close to 0.

###  D.3 Choice of $k$ for the $k$-NN Graph

The choice of $k$ for the $k$-NN graph may depend on the task at hand. From
the following experiments we see that for testing the equality of the $M$
distributions, the empirical criterion $k=0.10n$, for samples with up to a few
hundred observations [57], seems to provide a good choice. However, for
estimating $\eta$ using its empirical version $\hat{\eta}$, a much smaller $k$
is recommended — often $k=1$ works best.

Choice of $k$ in testing: In Figure 5 we show the empirical power (over 1000
replications) of our test statistics, when the level is set at 0.05, using
different $k$-NN graphs.

For $M=3$, the three distributions we consider are:

  1. 1.

Normal location problem: $P_{1}=N(\mathbf{0},I_{d})$,
$P_{2}=N(0.1\cdot\mathbf{1},I_{d})$, $P_{3}=N(0.2\cdot\mathbf{1},I_{d})$.

  2. 2.

Normal scale problem: $P_{1}=N(\mathbf{0},I_{d})$,
$P_{2}=N(\mathbf{0},1.5\cdot I_{d})$, $N(\mathbf{0},2\cdot I_{d})$.

  3. 3.

Non-Gaussian problem with $t$-distribution: $P_{1}$ has each entry following
independent $t(1)$ with noncentrality parameter $\delta$; $P_{2}$ and $P_{3}$
have each entry following independent $t(1)$. The dimension $d$ is set to be
16.

For $M=5$, the distributions we consider are:

  1. 1.

Normal location problem: $P_{1}=N(\mathbf{0},I_{d})$,
$P_{2}=N(0.05\cdot\mathbf{1},I_{d})$, $P_{3}=N(0.1\cdot\mathbf{1},I_{d})$,
$P_{4}=N(0.15\cdot\mathbf{1},I_{d})$, $P_{5}=N(0.2\cdot\mathbf{1},I_{d})$.

  2. 2.

Normal scale problem: $P_{1}=N(\mathbf{0},I_{d})$,
$P_{2}=N(\mathbf{0},1.25\cdot I_{d})$, $P_{3}=N(\mathbf{0},1.5\cdot I_{d})$,
$P_{4}=N(\mathbf{0},1.75\cdot I_{d})$, $P_{5}=N(\mathbf{0},2\cdot I_{d})$.

  3. 3.

Non-Gaussian problem with $t$-distribution: $P_{1}$ has each entry following
independent $t(1)$ with noncentrality parameter $\delta$; $P_{i}$, $i=2,3,4,5$
have each entry following $t(1)$. The dimension is set to be 16.

Here $n_{i}$, for $i=1,\ldots,M$, are set to be equal. Different combinations
of $n$ and $d$ are considered in our simulation experiments. It can be seen
from Figure 5 that the empirical criterion $k=0.1n$ indeed provides reasonably
good performance: For the normal scale problem with discrete kernel $k=0.1n$
is close to the optimal $k$, while for the problem with $t$-distribution the
power continues to increase until $k\approx 0.2n$. For the normal location
problem and the normal scale problem with kernel $K_{2}$ (recall that
$K_{2}(1,1)=10$, $K_{2}(2,2)=K_{2}(3,3)=1$, $K_{2}(i,j)=0$,for $i\neq j$), the
power still increases as $k$ increases even after $k\approx 0.5n$ (see the
plots corresponding to $n=300$ in Figure 5).

Figure 5: Power of the tests using different $k$-NN graphs.

Choice of $k$ in estimating $\eta$: Although power increases as $k$ increases
until or even after $k=0.1n$, the estimate of $\eta$ may no longer be accurate
with a large $k$, since the $k$-NNs of a point may no longer be close to that
point. When estimating $\eta$, we face a bias-variance trade-off. When
$\eta=0$, $\hat{\eta}$ is unbiased, so only the variance needs to be
controlled. When $\eta$ is large, the variance is of the order $n^{-1}$, but
the bias can be of the order $n^{-1/d}$, a typical distance between a point
and its nearest neighbor, and is the dominating term when the dimension gets
large. This issue of the bias could be alleviated by choosing a small $k$ — as
in other $k$-NN applications, a smaller $k$ usually has less bias (see Figure
6).

Another observation is that, in both simulated and real data, as $k$
increases, the actual value of $\hat{\eta}$ often decreases. In the most
extreme case where $k=n-1$, $\hat{\eta}$ is exactly 0. For the above reasons,
we would suggest that when $\eta$ or the dimension $d$ are suitably large, $k$
should be chosen to be as small as possible; often $k=1$ would be the best
choice. However, when $\eta$ is close to 0 and the dimension $d$ is small, $k$
can be suitably increased to reduce the variance of $\hat{\eta}$.

Figure 6 shows the mean of $\hat{\eta}$ under different settings, where we
consider

  1. 1.

$\eta=1$:

    1. (a)

$M=2$, $P_{1}={\rm uniform}[0,1]^{d}$, $P_{2}={\rm uniform}[0,1]^{d}+e_{1}$,

    2. (b)

$M=3$, $P_{1}={\rm uniform}[0,1]^{d}$, $P_{2}={\rm uniform}[0,1]^{d}+e_{1}$,
$P_{3}={\rm uniform}[0,1]^{d}+e_{2}$,

  2. 2.

$\eta=0.5$:

    1. (a)

$M=2$, $P_{1}={\rm uniform}[0,1]^{d}$, $P_{2}={\rm
uniform}[0,1]^{d}+\frac{1}{2}e_{1}$,

    2. (b)

$M=3$, $P_{1}=P_{2}={\rm uniform}[0,1]^{d}$, $P_{3}={\rm
uniform}[0,1]^{d}+e_{1}$,

  3. 3.

$\eta=0.1$:

    1. (a)

$M=2$, $P_{1}={\rm uniform}[0,1]^{d}$, $P_{2}={\rm
uniform}[0,1]^{d}+0.1e_{1}$,

    2. (b)

$M=3$, $P_{1}=P_{2}={\rm uniform}[0,1]^{d}$, $P_{3}={\rm
uniform}[0,1]^{d}+0.2e_{1}$.

Here $e_{i}$ denotes an all-zero vector except a 1 at the $i$-th entry.

Figure 6: The mean of $\hat{\eta}$ for different values of $\eta$ as $k$ (in
the $k$-NN graph) varies (for different choices of $M,n$, and $d$), computed
using $10^{4}$ independent replications.

For $M=2$, $\pi_{1}=\pi_{2}=1/2$, and for $M=3$,
$\pi_{1}=\pi_{2}=\pi_{3}=1/3$. We consider
$(n,d)\in\\{60,300\\}\times\\{2,200\\}$. A larger $d$ indicates a higher noise
level. In each of our $10^{4}$ replications, we have $n_{i}=n\pi_{i}$ samples
from $P_{i}$, and the mean of $\hat{\eta}$ is reported.

We observe in Table 6 that with a discrete kernel, $\hat{\eta}$ tends to
underestimate $\eta$, and the bias increases as $k$ increases. Here $k=1$
gives the most accurate estimate in the sense that it has the least bias in
all scenarios.

The empirical mean squared error of $(\hat{\eta}-\eta)^{2}$ is given in Figure
7.

Figure 7: The mean squared error of $\hat{\eta}$ for different values of
$\eta$ as $k$ (in the $k$-NN graph) varies (for different choices of $M,n$,
and $d$), computed using $10^{4}$ independent replications.

It can be seen that when $\eta=1$, the mean squared error (MSE) strictly
increases as $k$ increases; when $\eta=0.5$, $k=1$ is also close to the
optimal $k$. When $\eta=0.1$, which is close to 0, the MSE first decreases and
then increases as we increase $k$. Moreover, when $\eta=0.1$, the MSE is small
for all choices of $k$ compared to the previous examples where $\eta=0.5$ or
$1$. When $d=20$, it is clear from Figures 6 and 7 that the bias dominates the
variance in its contribution to the MSE. When we get closer to the null
$\eta=0$ (where $\hat{\eta}$ is unbiased), the bias decreases to 0, so the MSE
also decreases dramatically. For the case $\eta=0.1$, one may slightly
increase $k$ to 4 or 5 to reduce the variance and obtain the optimal MSE.

###  D.4 Validation of Theorem 4

Here we empirically verify our theoretical results in Theorem 4, which states
that under ${\rm H}_{0}$, the (conditional) variance of $\hat{\eta}$ converges
to a distribution-free constant not depending on the common distribution, thus
yielding an asymptotic test using this distribution-free variance. We first
examine the convergence of $\tilde{g}_{i}$ to $g_{i}$.

Figure 8: The computed $\tilde{g}_{3}$ as the sample size $n$ varies. The red
line shows $g_{3}$.

Consider $n$ observations from three identical distributions
$P_{1}=P_{2}=P_{3}\equiv N(\mathbf{0},I_{d})$. Each observation has
probability $1/3$ to be sampled from $P_{i}$, $i=1,2,3$; so
$\pi_{1}=\pi_{2}=\pi_{3}=1/3$. With a directed 1-NN graph, $g_{1}=g_{2}=1$,
and $g_{3}$ depends on $d$. The exact value of $g_{3}$ can be obtained from
Henze [39]. When $d=1$, it can be seen from the left panel of Figure 8 that
$\tilde{g}_{3}$ converges to $g_{3}$ as $n$ increases. When $n=2000$, the
difference between $\tilde{g}_{3}$ and $g_{3}$ is only 0.006. However, the
convergence becomes slower as the dimension $d$ becomes higher, as is seen
from the right panel of Figure 8. When $d=5$, the difference between
$\tilde{g}_{3}$ and $g_{3}$ is 0.043 when $n=2000$, and it takes $n=3.2\times
10^{6}$ observations to bring the error down to 0.006.

Next, since both Theorems 3 and 4 provide a CLT for $\hat{\eta}$, we compare
the quality of the approximation given by the two theorems.

Figure 9: Histograms of normalized $\hat{\eta}$, either by the permutation
variance in Theorem 3 or by the limiting variance in Theorem 4. The red curves
show the standard normal densities.

We consider $n_{1}=n_{2}=n_{3}=150$, $P_{1}=P_{2}=P_{3}\equiv
N(\mathbf{0},I_{d})$ and vary $d\in\\{1,5,10\\}$, and use a directed 1-NN
graph to compute $\hat{\eta}$. Figure 9 shows the histograms of $\hat{\eta}$
normalized its limiting variance, i.e.,
$\frac{\sqrt{n}\hat{\eta}}{\sigma_{\mathcal{G},K,d,\pi}}$, along with the
histograms of $\hat{\eta}$ normalized by its permutation variance (as in
Theorem 3), i.e., $\frac{\hat{\eta}}{\sqrt{{\rm
Var}(\hat{\eta}|\mathcal{F}_{n})}}$. The red curves are the standard normal
densities. It can be seen that all the histograms are close to the standard
normal density, and using the permutation variance seems to yield a better
approximation than using the limiting variance. Therefore, we recommend using
Theorem 3 over Theorem 4 for testing equality of distributions in practice.

## Acknowledgement

The authors would like to thank Bhaswar Bhattacharya and Nabarun Deb for a
number of useful comments and references.

## References

  * Ahle [2021] Ahle, T. D. (2021).  Sharp and simple bounds for the raw moments of the binomial and poisson distributions.  arXiv preprint arXiv:2103.17027.
  * Arias-Castro and Pelletier [2016] Arias-Castro, E. and B. Pelletier (2016).  On the consistency of the crossmatch test.  J. Statist. Plann. Inference 171, 184–190.
  * Arya et al. [2019] Arya, S., D. Mount, S. E. Kemp, and G. Jefferis (2019).  RANN: Fast Nearest Neighbour Search (Wraps ANN Library) Using L2 Metric.  R package version 2.6.1.
  * Azadkia and Chatterjee [2021] Azadkia, M. and S. Chatterjee (2021).  A simple measure of conditional dependence.  Ann. Statist. 49(6), 3070–3102.
  * Baryshnikov and Yukich [2005] Baryshnikov, Y. and J. E. Yukich (2005).  Gaussian limits for random measures in geometric probability.  Ann. Appl. Probab. 15(1A), 213–253.
  * Bentley [1975] Bentley, J. L. (1975).  Multidimensional binary search trees used for associative searching.  Communications of the ACM 18(9), 509–517.
  * Beran [1977] Beran, R. (1977).  Minimum Hellinger distance estimates for parametric models.  Ann. Statist. 5(3), 445–463.
  * Berlinet and Thomas-Agnan [2004] Berlinet, A. and C. Thomas-Agnan (2004).  Reproducing Kernel Hilbert Spaces in Probability and Statistics.  Kluwer Academic Publishers, Boston, MA.
  * Berndt and Clifford [1994] Berndt, D. J. and J. Clifford (1994).  Using dynamic time warping to find patterns in time series.  In KDD Workshop, Volume 10, pp. 359–370.
  * Bhattacharya [2019] Bhattacharya, B. B. (2019).  A general asymptotic framework for distribution-free graph-based two-sample tests.  J. R. Stat. Soc. Ser. B. Stat. Methodol. 81(3), 575–602.
  * Bhattacharya [2020] Bhattacharya, B. B. (2020).  Asymptotic distribution and detection thresholds for two-sample tests based on geometric graphs.  Ann. Statist. 48(5), 2879–2903.
  * Billingsley [2012] Billingsley, P. (2012).  Probability and Measure.  Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ.
  * Biswas et al. [2014] Biswas, M., M. Mukhopadhyay, and A. K. Ghosh (2014).  A distribution-free two-sample run test applicable to high-dimensional data.  Biometrika 101(4), 913–926.
  * Bloemena [1964] Bloemena, A. R. (1964).  Sampling from a Graph, Volume 2 of Mathematical Centre Tracts.  Mathematisch Centrum, Amsterdam.