Testing goodness of fit for point processes via topological data analysis

TL;DR

This paper develops new goodness-of-fit tests for point patterns using topological data analysis, specifically persistent Betti numbers, and demonstrates their effectiveness through simulations and a neuroscience application.

Contribution

It introduces a novel statistical framework based on persistent Betti numbers for assessing point process models, including theoretical conditions and practical testing procedures.

Findings

Persistent Betti numbers are asymptotically Gaussian for large observation windows.

The proposed tests outperform global envelope tests in power.

Application to neuroscience data demonstrates practical utility.

Abstract

We introduce tests for the goodness of fit of point patterns via methods from topological data analysis. More precisely, the persistent Betti numbers give rise to a bivariate functional summary statistic for observed point patterns that is asymptotically Gaussian in large observation windows. We analyze the power of tests derived from this statistic on simulated point patterns and compare its performance with global envelope tests. Finally, we apply the tests to a point pattern from an application context in neuroscience. As the main methodological contribution, we derive sufficient conditions for a functional central limit theorem on bounded persistent Betti numbers of point processes with exponential decay of correlations.

Tables

Table 1. Table 1. Rejection rates for the test statistics T 𝖢 subscript 𝑇 𝖢 T_{\mathsf{C}} and T 𝖫 subscript 𝑇 𝖫 T_{\mathsf{L}} under the null model and the alternatives.

	$\mathrm{𝖯𝗈𝗂}$	$\mathrm{𝖬𝖺𝗍𝖢}$	$\mathrm{𝖲𝗍𝗋}$
$T_{𝖢}$	5.1%	59.3%	60.7%
$T_{𝖫}$	4.8%	94.7%	71.4%

Table 2. Table 2. Rejection rates for envelope tests based on the L 𝐿 L function and cluster- and loop-based functional statistics.

	$\mathrm{𝖬𝖺𝗍𝖢}$	$\mathrm{𝖲𝗍𝗋}$
Ripley’s $L$	42.6%	20.5%
$\mathrm{𝖢𝗅𝗎𝗌𝗍𝖾𝗋}$	41.5%	26.3%
$\mathrm{𝖫𝗈𝗈𝗉}$	27.0%	32.2%

Equations

$$U_r(\mathcal{X})=\bigcup _{x\in \mathcal{X}}{B}_r(x).$$

$$R(x,y)=\inf \{r>0:{\mathcal{C}}_r(x)={\mathcal{C}}_r(y)\}.$$

$${\mathsf{PD}}^{M,q}(\mathcal{X})=\sum _{i\in I^{M,q}(\mathcal{X})}{\delta }_{(B_i^M,D_i^M)},$$

$${\beta }_{b,d}^{M,q}(\mathcal{X})={\mathsf{PD}}^{M,q}(\mathcal{X})([0,b]\times [d,r_{\mathsf{f}}])$$

$$\displaystyle\mathbb{E}\Big{[}\prod_{i\leq p}\mathcal{P}(A_{i})\Big{]}=\int_{A_{1}\times\cdots\times A_{p}}\rho^{(p)}(\boldsymbol{x}){\rm d}\boldsymbol{x}$$

$$\displaystyle\mathbb{E}\Big{[}\sum_{{(X_{1},\dots,X_{p})\in\mathcal{P}_{\neq}^{p}}}f(X_{1},\dots,X_{p};\mathcal{P})\Big{]}=\int_{\mathbb{R}^{2p}}\mathbb{E}_{\boldsymbol{x}}^{!}[f(\boldsymbol{x};\mathcal{P}\cup\boldsymbol{x})]\rho^{(p)}(\boldsymbol{x}){\rm d}\boldsymbol{x},$$

$$\underset{\begin{array}{c}l\le p\\ \mathbf{x}\in {\mathbb{R}}^{2l}\end{array}}{\sup }{\mathbb{E}}_{\mathbf{x}}^{!}[\mathcal{P}(W_1{)}^p]<\infty .$$

$$\langle f,\mathsf{PD}^{M,q}(\mathcal{P}_{n})\rangle=\int_{[0,r_{\mathsf{f}}]^{2}}f(b,d)\mathsf{PD}^{M,q}(\mathcal{P}_{n})({\rm d}b,{\rm d}d)=\sum_{i\in I^{M,q}(\mathcal{P}_{n})}f\big{(}B_{i}^{M},D_{i}^{M}\big{)}$$

$$\frac{⟨f,{\mathsf{PD}}^{M,q}({\mathcal{P}}_n)⟩-\mathbb{E}[⟨f,{\mathsf{PD}}^{M,q}({\mathcal{P}}_n)⟩]}{\sqrt{\mathsf{Var}(⟨f,{\mathsf{PD}}^{M,q}({\mathcal{P}}_n)⟩)}}$$

$$\mathbb{E}\big{[}\min_{i\in\{1,2\}}\mathbb{P}\big{(}\mathcal{P}_{r_{\mathsf{AC}}^{2}}\in E_{i}\,|\,\Lambda,\mathcal{P}\setminus W_{r_{\mathsf{AC}}^{2}}\big{)}\big{]}>0.$$

$$\big{\{}n^{-1/2}\big{(}\beta^{M,0}_{d}(\mathcal{P}_{n})-\mathbb{E}[\beta^{M,0}_{d}(\mathcal{P}_{n})]\big{)}\big{\}}_{d\leq r_{\mathsf{f}}}$$

$$\big{\{}n^{-1/2}\big{(}\beta^{M,1}_{b,d}(\mathcal{P}_{n})-\mathbb{E}[\beta^{M,1}_{b,d}(\mathcal{P}_{n})]\big{)}\big{\}}_{b,d\leq r_{\mathsf{f}}}$$

$${\mathsf{APF}}_r^{M,0}({\mathcal{P}}_n)=\sum _{i\in I^{M,0}({\mathcal{P}}_n)}{D}_i^M\text{mathbbmss}1\{D_i^M\le r\}$$

$${\mathsf{APF}}_r^{M,1}({\mathcal{P}}_n)=\sum _{i\in I^{M,1}({\mathcal{P}}_n)}(D_i^M-B_i^M)\text{mathbbmss}1\{B_i^M\le r\}.$$

$$\mathsf{dist}(\mathbf{x},{\mathbf{x}}^{\prime })=\underset{\begin{array}{c}i\le p\\ j\le q\end{array}}{\inf }|x_i-x_{p+j}|.$$

$$|{\rho }^{(p+q)}(\mathbf{x}\cup {\mathbf{x}}^{\prime })-{\rho }^{(p)}(\mathbf{x}){\rho }^{(q)}({\mathbf{x}}^{\prime })|\le (p+q{)}^{a(p+q)}\varphi (\mathsf{dist}(\mathbf{x},{\mathbf{x}}^{\prime }))$$

$${\rho }^{(j)}(u_1,\dots ,u_j)=\exp \left(j\mu +\frac{jc(0)}{2}\right)\prod _{1\le i<i^{\prime }\le j}\exp (c(u_i-u_{i^{\prime }})).$$

$${\mathbb{E}}_{\mathbf{x}}^{!}[\mathcal{P}(W_1{)}^p]=\mathbb{E}[{\mathcal{P}}_{\mathbf{x}}(W_1{)}^p]=\sum _{1\le j\le p}{\Delta }_{j,l,p}{\int }_{W_1^j}{\rho }_{\mathbf{x}}^{(j)}(u_1,\dots ,u_j)\mathrm{d}{u}_1\cdots \mathrm{d}{u}_j$$

$$\underset{\mathbf{x}\in {\mathbb{R}}^{2l}}{\sup }{\int }_{W_1^j}{\rho }_{\mathbf{x}}^{(j)}(u_1,\dots ,u_j)\mathrm{d}{u}_1\cdots \mathrm{d}{u}_j<\infty ,$$

$$\rho^{(j)}_{\boldsymbol{x}}(u_{1},\dots,u_{j})=\exp\Big{(}j\mu+\frac{jc(0)}{2}+\sum_{\begin{subarray}{c}1\leq i\leq j\\ 1\leq k\leq l\end{subarray}}c(u_{i},x_{k})\Big{)}\prod_{1\leq i<i^{\prime}\leq j}\exp(c(u_{i}-u_{i^{\prime}})).$$

$$f_{\mathbf{\Lambda }}(\varphi )=\exp (|W_{r_{\mathsf{AC}}^2}|-\mathbf{\Lambda }(W_{r_{\mathsf{AC}}^2}))\prod _{x\in \varphi }\exp (Y(x)),$$

$$\lambda (x)=\gamma \eta (B_R(x)).$$

$$\mathbb{E}^{!}_{\boldsymbol{x}}[\mathcal{P}(W_{1})^{p}]=\frac{1}{\mathbb{E}\big{[}\prod_{i\leq l}\lambda(x_{i})\big{]}}\cdot\mathbb{E}\big{[}\mathcal{P}(W_{1})^{p}\prod_{i\leq l}\lambda(x_{i})\big{]}.$$

$$\displaystyle\mathbb{E}\big{[}\prod_{i\leq l}\lambda(x_{i})\big{]}\geq\prod_{i\leq l}\mathbb{E}[\lambda(x_{i})]=(\gamma\pi R^{2})^{l},$$

$$\displaystyle\mathbb{E}\Big{[}\mathcal{P}(W_{1})^{p}\prod_{i\leq l}\lambda(x_{i})\Big{]}\leq\mathbb{E}\big{[}\mathcal{P}(W_{1})^{p(l+1)}\big{]}^{1/(l+1)}\mathbb{E}[\lambda(o)^{l+1}]^{l/(l+1)}.$$

$$\displaystyle\mathbb{E}\big{[}\mathcal{P}(W_{1})^{p(l+1)}\big{]}$$

$$\displaystyle\leq\mathbb{E}\Big{[}\prod_{x\in\eta\cap(W_{1}\oplus B_{R}(o))}e^{p(l+1)\#\Phi_{x}}\Big{]}$$

$$\displaystyle=\exp\big{(}\gamma|W_{1}\oplus B_{R}(o)|(\mathbb{E}[e^{p(l+1)\#\Phi_{0}}]-1)\big{)}.$$

$$f_{\eta }(\varphi )=\gamma \exp (|W_{r_{\mathsf{AC}}^2}|-\mathbf{\Lambda }(W_{r_{\mathsf{AC}}^2}))\prod _{x\in \varphi }\eta (B_R(x))$$

$$\mathcal{E}=\{W_{r_{\mathsf{AC}}^2}\subset \eta \oplus B_R(o)\},$$

$$\underset{i\in \{1,2\}}{\min }\mathbb{P}({\mathcal{P}}_{r_{\mathsf{AC}}^2}\in E_i|\eta )1_{\mathcal{E}}(\eta )>0.$$

Full text

1. Introduction

Topological data analysis (TDA) provides insights into a variety of datasets by capturing their most salient properties via refined topological features. Since the mathematical field of topology specializes in describing invariants of objects independently of the choice of a precise metric, these features are robust against small perturbations or different embeddings of the object [11, 12]. Among the most classical topological invariants are the Betti numbers. Loosely speaking, they capture the number of $k$-dimensional holes of the investigated structure. TDA refines this idea substantially by constructing filtrations and tracing when topological features appear and disappear. In point pattern analysis, simplicial complexes are built so that they are topologically equivalent to a union of balls with the same radius and centered at the data points, see the first three panels of Figure 1. As the radius increases, a sequence of simplicial complexes is then defined. Examples of such complexes are the basic Čech complex or the more elaborate $\alpha$-complex, which is based on the Delaunay triangulation, see [18]. In that framework, 1-dimensional features correspond to loops in the simplicial complexes while 0-dimensional features correspond to connected components. When moving up in the filtration, additional edges appear and at some point create new loops. On the other hand, more and more triangles also appear, thereby causing completely filled loops to disappear. Usually, the filtration is indexed by time, and we refer to the appearance and disappearance of features as births and deaths. We refer the reader to [18] for a detailed presentation of these concepts. The persistence diagram visualizes the time points when the features are born and die, see the bottom-right panel in Figure 1. Persistent Betti numbers count the number of events in upper-left blocks of the persistence diagram and are also illustrated in the figure.

In this paper, we leverage persistent Betti numbers to derive goodness-of-fit tests for planar point processes. In this setting, the abstract general definition of persistent Betti numbers gives way to a clear geometric intuition induced by a picture of growing disks centered at the points of the pattern and all having radius $r$, corresponding to the index of the filtration. Features of dimension 0 correspond to connected components in the union of balls, interpreted as point clusters, whereas boundaries of the complement set can be considered as the loops forming the 1-dimensional features. Since the notion of clusters in the sense of connected components lies at the heart of persistent Betti numbers in degree 0, they become highly attractive as a tool to detect clustering in point patterns. Our tests are based on a novel functional central limit theorem (CLT) for the persistent Betti numbers in large domains. The present work embeds into two active streams of current research.

First, now that TDA has become widely adopted, the community is vigorously working towards putting the approach on a firm statistical foundation paving the way for hypothesis testing. On the one hand, this encompasses large-sample Monte Carlo tests when working on a fixed domain [6, 9, 13]. Although these tests are highly flexible, the test statistics under the null hypothesis must be re-computed each time when testing observations in a different window. In large domains, this becomes time-consuming. On the other hand, there has been substantial progress towards establishing CLTs in large domains for functionals related to persistent Betti numbers [39, 40, 29, 35, 25]. However, these results are restricted to the null hypothesis of complete spatial randomness – i.e., the Poisson point process – and establish asymptotic Gaussianity on a multivariate, but not on a functional level. Our proof of a functional CLT is based on recently developed stabilization techniques for point processes with exponential decay of correlations [8]. As explained in the final section of [10], the main technical step towards a functional CLT are bounds on the cumulants.

Second, the introduction of global rank envelope tests has lead to a novel surge of research activity in goodness-of-fit tests for point processes [34]. One of the reasons for their popularity is that they rely on functional summary statistics rather than scalar quantities. Thus, they reveal a substantially more fine-grained picture of the underlying point pattern. In the overwhelming majority of cases, variants of the $K$-function are used as a functional summary statistic, thereby essentially capturing the relative density of point pairs at different distances. Here, the persistent Betti numbers offer an opportunity to augment the basic second-order information by more refined characteristics of the data. Still, even for classical summary statistics, rigorous limit theorems in large domains remain scarce. For instance, a functional central limit theorem of the estimated $K$-function is proven in detail only for the Poisson point process in [22] and an extension to $\alpha$-determinantal point processes is outlined in [23].

The rest of the manuscript is organized as follows. First, in Section 2, we introduce the concepts of $M$-bounded persistence diagrams and $M$-bounded persistent Betti numbers. Next, in Section 3, we state the two main results of the paper, a CLT for the $M$-bounded persistence diagram and a functional CLT for the $M$-bounded persistent Betti numbers. In Section 4, we provide specific examples of point processes satisfying the conditions of the main results. Sections 5 and 6 explore TDA-based tests for simulated and real datasets, respectively. Finally, Section 7 summarizes the findings and points to possible avenues of future research. The proofs of the main results are deferred to Sections 8 and 9 of the appendix.

2. $M$-bounded persistent Betti numbers

For a locally finite point set $\mathcal{X}\subset\mathbb{R}^{2}$, persistent Betti numbers provide refined measures for the amount of clusters and voids on varying length scales. More precisely, we let

[TABLE]

denote the union of closed disks of radius $r\geq 0$ centered at points in $\mathcal{X}$. A 0-dimensional topological feature is a connected component of this union, corresponding to a cluster of points in $\mathcal{X}$, while a 1-dimensional feature can be thought of as a bounded connected component of the background space, often identified with its boundary loop, and describes a vacant area in the plane. As the disks grow, new features arise and vanish; we say that they are born and die again. The persistent Betti numbers quantify this evolution of clusters and loops. Henceforth, we consider the persistence diagram only until a fixed deterministic radius $r_{\mathsf{f}}\geq 0$.

As $r$ approaches the critical radius for continuum percolation, long-range phenomena emerge [31]. Thus, determining whether two points are connected could require exploring large regions in space. While useful quantitative bounds on cluster sizes are known for Poisson point processes [1], for more general classes of point processes the picture remains opaque and research is currently at a very early stage [27, 7]. Recently, a central limit theorem for persistent Betti numbers has been established in the Poisson setting [29, 25], but for general point processes the long-range interactions pose a formidable obstacle towards proving a fully-fledged functional CLT.

From a more practical point of view, these long-range dependencies are of less concern. Although large features can carry interesting information, we expect that spatially bounded topological features already provide a versatile tool for the statistical analysis of both simulated point patterns and real datasets, even when focusing only on features of a bounded size. For that purpose, we concentrate on features whose spatial diameter does not exceed a large deterministic threshold $M$.

To define these $M$-bounded features, we introduce the Gilbert graph $G_{r}(\mathcal{X})$ on the vertex set $\mathcal{X}$. The Gilbert graph $G_{r}(\mathcal{X})$ has for vertices the points in $\mathcal{X}$ and two points are connected by an edge if the distance between them is at most $2r$ or, equivalently, if the two disks of radius $r$ centered at the points intersect.

2.1. $M$-bounded clusters

The 0-dimensional $M$-bounded features alive at time $r>0$ are the connected components of $G_{r}(\mathcal{X})$ with diameter at most $M$. Starting at $r=0$, all points belong to separate connected components that merge into larger clusters when $r$ increases. We thus say that all components are born at time 0.

To define the death time of a component, let $\mathcal{C}_{r}(x)$ denote the connected component of $x\in\mathcal{X}$ in $G_{r}(\mathcal{X})$. The components of $x,y\in\mathcal{X}$ meet at time

[TABLE]

Then, the death time of $x\in\mathcal{X}$ is the smallest $R(x,y)$ such that the spatial diameter of $\mathcal{C}_{r}(x)$ exceeds $M$ or such that $P_{x}$ is lexicographically larger than $P_{y}$, where $P_{x},P_{y}$ are the points of $\mathcal{C}_{r}(x)\cap\mathcal{X}$ and $\mathcal{C}_{r}(y)\cap\mathcal{X}$ whose associated disks meet at time $R(x,y)$. This ordering determines which component dies when two of them meet.

2.2. $M$-bounded loops

Next, we introduce 1-dimensional features. At time $r>0$, these correspond to holes, i.e., bounded connected components in the vacant phase $V_{r}(\mathcal{X})=\mathbb{R}^{2}\setminus U_{r}(\mathcal{X})$. In contrast to the clusters, there are no holes at time [math], so that both birth and death times must be specified. Moreover, it needs to be defined how holes are related for different radii $r$.

The death time of a hole $H_{s}$ in $V_{s}(\mathcal{X})$ is the first time $r>s$ when the hole is completely covered by disks, i.e., $H_{s}\subseteq U_{r}(\mathcal{X})$. We identify a hole $H$ with the point $p(H)$ that is covered last. Thus, holes $H_{s}$ in $V_{s}(\mathcal{X})$ and $H_{r}$ in $V_{r}(\mathcal{X})$, are identified if $p(H_{r})=p(H_{s})$.

New holes in $V_{r}(\mathcal{X})$ can only appear when two balls merge, which corresponds to including a new edge in $G_{r}(\mathcal{X})$. When a new hole is formed, it can happen in two ways: either a finite component is separated from the infinite component, or an existing hole is split in two. In both cases, we define the size of the newly created hole(s) as follows: Let $x_{1},\dots,x_{k}\in\mathcal{X}$ be the points in $\mathcal{X}$ such that the disks of radius $r$ around the points intersect the boundary of the hole $H$ in $V_{r}(\mathcal{X})$. Then, the size of $H$ is the diameter of the set $\{x_{1},\dots,x_{k}\}$. The size remains unchanged until the next time the hole is split into smaller pieces. Then the size is recomputed for both new holes. This definition ensures that the size decreases when the balls grow and only changes when a new edge is added to $G_{r}(\mathcal{X})$.

The birth time of a hole $H$ is the minimal $s$ such that there is a hole $H_{s}$ in $V_{s}(\mathcal{X})$ with $p(H)=p(H_{s})$ and size less than $M$. By an $M$-bounded loop, we mean a loop with size lower than $M$.

2.3. The persistence diagram

We now adapt the definition of the persistence diagram in [25] to only include $M$-bounded features. That is, we define the $q$th $M$-bounded persistence diagram, $q\in\{0,1\}$, as the empirical measure

[TABLE]

where $I^{M,q}(\mathcal{X})$ is an index set over all $M$-bounded $q$-dimensional features that die before time $r_{\mathsf{f}}$ and $B_{i}^{M},D_{i}^{M}$ are the birth and death times of the $i$th feature. Then, the $q$th $M$-bounded persistent Betti numbers

[TABLE]

are the number of $M$-bounded features born before time $b\geq 0$ and dead after time $d\leq r_{\mathsf{f}}$. When $q=0$, all features are born at time 0, so that only death times are relevant. Hence, we write $\beta^{M,0}_{d}$ instead of the more verbose $\beta^{M,0}_{b,d}$.

3. Main results

Henceforth, $\mathcal{P}$ denotes a simple stationary point process in $\mathbb{R}^{2}$. We think of $\mathcal{P}$ as a random variable taking values in the space of locally finite subsets $\mathcal{N}$ of $\mathbb{R}^{2}$ endowed with the smallest $\sigma$-algebra $\mathfrak{N}$ such that the number of points in any given Borel set becomes measurable. Throughout the manuscript, we assume that the factorial moment measures exist and are absolutely continuous. In particular, writing $\boldsymbol{x}=(x_{1},\dots,x_{p})\in\mathbb{R}^{2p}$, the $p$th factorial moment density $\rho^{(p)}$ is determined via the identity

[TABLE]

for any pairwise disjoint bounded Borel sets $A_{1},\dots,A_{p}\subset\mathbb{R}^{2}$, where $\mathcal{P}(A_{i})$ denotes the number of points of $\mathcal{P}$ in $A_{i}$. Moreover, as we rely on the framework of [8], we also require that $\mathcal{P}$ exhibits exponential decay of correlations. Loosely speaking this expresses an approximate factorization of the factorial moment densities and is made precise in Section 4 below. Many of the most prominent examples of point processes appearing in spatial statistics exhibit exponential decay of correlations [8, Section 2.2].

Our first main result is a CLT for the persistence diagram built on the restriction $\mathcal{P}_{n}=\mathcal{P}\cap W_{n}$ of the point process $\mathcal{P}$ to a large observation window $W_{n}=[-\sqrt{n}/2,\sqrt{n}/2]^{2}$. With a slight abuse of notation, we write $\mathcal{P}\cup\boldsymbol{x}=\mathcal{P}\cup\{x_{1},\dots,x_{p}\}$. To prove the CLT, we impose an additional condition concerning moments under the reduced $p$-point Palm distribution $\mathbb{P}_{\boldsymbol{x}}^{!}$. We recall that this distribution is determined via

[TABLE]

for any bounded measurable $f:\mathbb{R}^{2p}\times\mathcal{N}\to\mathbb{R}$, where $\mathcal{P}_{\neq}^{p}$ denotes $p$-tuples of pairwise distinct points in $\mathcal{P}$. Then, we impose the following moment condition.

• (M)

For every $p\geq 1$

[TABLE]

To state the CLT for the persistence diagram precisely, we let

[TABLE]

denote the integral of a bounded measurable function $f:\,[0,r_{\mathsf{f}}]^{2}\to\mathbb{R}$ with respect to the measure $\mathsf{PD}^{M,q}(\mathcal{P}_{n})$.

Theorem 3.1** (CLT for persistence diagrams).**

Let $M>0$, $q\in\{0,1\}$ and $f:\,[0,r_{\mathsf{f}}]^{2}\to\mathbb{R}$ be a bounded measurable function. Assume that $\mathcal{P}$ exhibits exponential decay of correlations and satisfies condition (M). Furthermore, assume that $\liminf_{n\to\infty}\mathsf{Var}(\langle f,\mathsf{PD}^{M,q}(\mathcal{P}_{n})\rangle)n^{-\nu}=\infty$ for some $\nu>0$. Then,

[TABLE]

converges in distribution to a standard normal random variable as $n\to\infty$.

In order to derive a functional CLT for the persistent Betti numbers, we add a further constraint on $\mathcal{P}$, which is needed to establish a lower bound on the variance via a conditioning argument in the vein of [38, Lemma 4.3]. For this purpose, we consider a random measure $\Lambda$, which is jointly stationary with $\mathcal{P}$ and which we think of as capturing additional useful information on the dependence structure of $\mathcal{P}$. For instance, if $\mathcal{P}$ is a Cox point process, we choose $\Lambda$ to be the random intensity measure. If $\mathcal{P}$ is a Poisson cluster process, then $\Lambda$ would describe the cluster centers. If the dependence structure is exceptionally simple, it is also possible to take $\Lambda=0$. The idea of using additional information is motivated from conditioning on the spatially refined information coming from the clan-of-ancestors construction in Gibbsian point processes [38].

The point process $\mathcal{P}$ is conditionally $m$-dependent if $\mathcal{P}\cap A$ and $\mathcal{P}\cap A^{\prime}$ are conditionally independent given $\sigma(\Lambda,\mathcal{P}\cap{A^{\prime\prime}})$ for any bounded Borel sets $A,A^{\prime},A^{\prime\prime}\subset\mathbb{R}^{2}$ such that the distance between $A$ and $A^{\prime}$ is larger than some $m>0$. Here, $\sigma(\Lambda,\mathcal{P}\cap{A^{\prime\prime}})$ denote the $\sigma$-algebra generated by $\Lambda$ and $\mathcal{P}\cap{A^{\prime\prime}}$.

Finally, we impose an absolute continuity-type assumption on the Poisson point process in a fixed box with respect to $\mathcal{P}$ when conditioned on $\Lambda$ and the outside points. More precisely, we demand that there exists $r_{\mathsf{AC}}>6M$ with the following property, where $\mathcal{Q}$ denotes a homogeneous Poisson point process in the window $W_{r_{\mathsf{AC}}^{2}}$.

1. (AC)

Let $E_{1},E_{2}\in\mathfrak{N}$ be such that $\min_{i\in\{1,2\}}\mathbb{P}(\mathcal{Q}\in E_{i})>0$. Then,

[TABLE]

Although (AC) appears technical, Section 4 illustrates that it is tractable for many commonly used point processes.

Since the persistent Betti numbers exhibit jumps at the birth- and death times of features, we work in the Skorokhod topology [5, Section 14].

Theorem 3.2** (Functional CLT for persistent Betti numbers).**

Let $M>0$ and $\mathcal{P}$ be a conditionally $m$-dependent point process with exponential decay of correlations and satisfying conditions (M) and (AC). Then, the following convergence statements hold true.

1. q=0.

The one-dimensional process

[TABLE]

converges weakly in Skorokhod topology to a centered Gaussian process. 2. q=1.

The two-dimensional process

[TABLE]

converges weakly in Skorokhod topology to a centered Gaussian process.

Additionally, [8, Theorem 1.12] implies convergence of the rescaled variances. While Theorem 3.1 is an adaptation of [8], Theorem 3.2 is much more delicate. As an application of Theorem 3.2, we obtain a functional CLT for the following two characteristics, which are modified variants of the accumulated persistence function from [6]:

[TABLE]

and

[TABLE]

Corollary 3.3** (Functional CLT for the APF).**

Let $M>0$ and $\mathcal{P}$ be as in Theorem 3.2. Then, both \big{\{}n^{-1/2}(\mathsf{APF}^{M,0}_{r}(\mathcal{P}_{n})-\mathbb{E}[\mathsf{APF}^{M,0}_{r}(\mathcal{P}_{n})])\big{\}}_{r\leq r_{\mathsf{f}}} and \big{\{}n^{-1/2}(\mathsf{APF}^{M,1}_{r}(\mathcal{P}_{n})-\mathbb{E}[\mathsf{APF}^{M,1}_{r}(\mathcal{P}_{n})])\big{\}}_{r\leq r_{\mathsf{f}}} converge to centered Gaussian processes.

4. Examples of point processes

In this section, we give examples of point processes which satisfy the assumptions of our main theorems. More precisely, we show that log-Gaussian Cox processes with compactly supported covariance functions and Matérn cluster processes both satisfy the conditions of Theorems 3.1 and 3.2. We also show that the Ginibre point process satisfies the conditions of Theorem 3.1.

Conversely, we do not expect that hard-core point processes satisfy the functional central limit theorem in the generality of Theorem 3.2. Indeed, hard-core conditions put a strict lower bound on the death time of clusters and the birth time of loops. However, we believe that suitable repulsive point processes, where the hard-core conditions only need to be imposed with a certain probability can be embedded in the framework of Theorem 3.2.

We first recall the definition of exponential decay of correlations from [8]. To this end, we define the separation distance between $\boldsymbol{x}=\{x_{1},\dots,x_{p}\}\subset\mathbb{R}^{2}$ and $\boldsymbol{x}^{\prime}=\{x_{p+1},\dots,x_{p+q}\}\subset\mathbb{R}^{2}$ as in [8, Formula (1.3)] via

[TABLE]

Definition 4.1**.**

Let $\mathcal{P}$ be a stationary point process in $\mathbb{R}^{2}$, such that the $k$-point correlation function $\rho^{(k)}$ exists for all $k\geq 1$. Then, $\mathcal{P}$ exhibits exponential decay of correlations if there exist $a<1$, $\phi:[0,\infty)\to[0,\infty)$ such that

1. (1)

$\lim_{t\to\infty}t^{n}\phi(t)=0$ for all $n\geq 1$, 2. (2)

$\liminf_{t\to\infty}\log\phi(t)/t^{b}<0$ for some $b>0$, 3. (3)

[TABLE]

for any $\boldsymbol{x}=\{x_{1},\dots,x_{p}\},\boldsymbol{x}^{\prime}=\{x_{p+1},\dots,x_{p+q}\}\subset\mathbb{R}^{2}$.

4.1. Log-Gaussian Cox process

Let $Y=\{Y(x)\}_{x\in\mathbb{R}^{2}}$ be a stationary Gaussian process with mean $\mu\in\mathbb{R}$ and covariance function $c(x,x^{\prime})=c(x-x^{\prime})$. Then, the random measure on $\mathbb{R}^{2}$ defined as $\mathbf{\Lambda}(B)=\int_{B}\exp(Y(x)){\rm d}x$, for any Borel subset $B\subset\mathbb{R}^{2}$ has moments of any order. Let $\mathcal{P}$ be a Cox process with random intensity measure $\mathbf{\Lambda}$, referred to as a Log-Gaussian Cox process. By [15, Equation (7)], the factorial moment densities of $\mathcal{P}$ are given by

[TABLE]

To apply Theorems 3.1 and 3.2, we assume that $c$ is bounded and of compact support, which ensures that $\mathcal{P}$ exhibits exponential decay of correlation.

We show below that condition $\mathbf{(M)}$ is satisfied. Let $\boldsymbol{x}=(x_{1},\dots,x_{l})\in\mathbb{R}^{2l}$. According to [16, Theorem 1], the Log-Gaussian Cox process $\mathcal{P}$ under the reduced Palm version is also a Log-Gaussian Cox process $\mathcal{P}_{\boldsymbol{x}}$ with underlying Gaussian process $Y_{\boldsymbol{x}}(x)=Y(x)+\sum_{i\leq l}c(x,x_{i})$. According to [17, Equation (5.4.5)],

[TABLE]

for suitable coefficients $\Delta_{j,l,p}\in\mathbb{R}$, where $\rho^{(j)}_{\boldsymbol{x}}(u_{1},\dots,u_{j})$ denotes the $j$th factorial moment density with respect to $\mathcal{P}_{\boldsymbol{x}}$. Therefore, it is enough to prove that

[TABLE]

for all $j,l\geq 1$. Now, Equation (8) in [15] gives that

[TABLE]

where the right-hand side is bounded as $\mu$ and $c$ are bounded independently of $\boldsymbol{x}$. This verifies condition (M).

Since conditionally on $\mathbf{\Lambda}$, the point process $\mathcal{P}$ is a Poisson point process, the conditional $m$-dependence property holds with $\Lambda=\mathbf{\Lambda}$.

It remains to verify condition $\mathbf{(AC)}$. By [33, Equation (6.2)], conditionally on $\Lambda$, the distribution of the point process $\mathcal{P}_{r_{\mathsf{AC}}^{2}}$ admits the density with respect to a homogeneous Poisson point process $\mathcal{Q}$ with intensity 1 in $W_{r_{\mathsf{AC}}^{2}}$ given by

[TABLE]

where $\phi\in\mathfrak{N}$. In particular, $f_{\mathbf{\Lambda}}(\phi)$ is strictly positive for all $\phi$. Therefore, if $E_{1},E_{2}$ are two events such that $\min_{i\in\{1,2\}}\mathbb{P}(\mathcal{Q}\in E_{i})>0$, then $\mathbb{P}(\mathcal{P}_{r_{\mathsf{AC}}^{2}}\in E_{i}\,|\,\Lambda=\mathbf{\Lambda})>0$. This verifies condition $\mathbf{(AC)}$.

4.2. Matérn cluster process

Let $\eta$ be a homogeneous Poisson point process in $\mathbb{R}^{2}$ with intensity $\gamma>0$. Given a realization of $\eta$, we define a family of independent point processes $(\Phi_{x})_{x\in\eta}$, where $\Phi_{x}$, $x\in\eta$, is a homogeneous Poisson point process with intensity 1 in the disk $B_{R}(x)$ of radius $R>0$ centered at $x\in\mathbb{R}^{2}$. The point process $\mathcal{P}=\bigcup_{x\in\eta}\Phi_{x}$ is referred to as a Matérn cluster process. Since $\mathcal{P}$ is $2R$-dependent, it exhibits exponential decay of correlations.

Next, we verify condition (M). For this purpose, we deduce from [16, Section 5.3.2] that a Matérn cluster process is a Cox process whose random intensity measure $\mathbf{\Lambda}$ has as density the random field $(\lambda(x))_{x\in\mathbb{R}^{2}}$ given by

[TABLE]

Now, let $\boldsymbol{x}=(x_{1},\dots,x_{l})\in\mathbb{R}^{2l}$ and $p\geq 1$ be fixed. From [16, Equations (19) and (20)] we obtain that

[TABLE]

Since $\eta(B_{R}(x))$ is increasing in $\eta$ for every $x\in\mathbb{R}^{2}$ in the sense of [30], the Harris-FKG inequality [30, Theorem 20.4] gives that

[TABLE]

where we used that $\lambda(x_{i})=\eta(B_{R}(x_{i}))$ is a Poisson random variable with parameter $\pi R^{2}$. In order to bound \mathbb{E}\Big{[}\mathcal{P}(W_{1})^{p}\prod_{i\leq l}\lambda(x_{i})\Big{]}, we first apply the Hölder inequality and stationarity, to arrive at

[TABLE]

First, $\mathbb{E}[\lambda(o)^{l+1}]=\mathbb{E}[\eta(B_{R}(o))^{l+1}]$ is finite since $\eta$ is a Poisson point process. For the remaining part, we note that $\mathcal{P}(W_{1})\leq\sum_{y\in\eta\cap(W_{1}\oplus B_{R}(o))}\#\Phi_{x}$, where $W_{1}\oplus B_{R}(o)=\{x+y:\,x\in W_{1},\,y\in B_{R}(o)\}$ denotes the Minkowski sum. Hence,

[TABLE]

where $\Phi_{0}$ is a homogeneous Poisson point process of intensity 1 in the disk $B_{R}(o)$ [30, Theorem 3.9]. Again, since $\#\Phi_{0}$ is a Poisson random variable with parameter $\pi R^{2}$, the latter expression is finite. Taking the supremum over all $\boldsymbol{x}$ and all $l\leq p$, this verifies condition (M). The point process $\mathcal{P}$ is also conditionally $m$-dependent, by taking $m=2R$ and $\Lambda=\eta$.

It remains to prove $\mathbf{(AC)}$. By [33, Equation (6.2)], conditional on $\Lambda=\eta$, the distribution of $\mathcal{P}_{r_{\mathsf{AC}}^{2}}$ admits the density

[TABLE]

with respect to the distribution of a homogeneous Poisson point process. Now, consider the event on the event

[TABLE]

the density $f_{\eta}$ is positive. Therefore, if $E_{1},E_{2}$ are such that $\min_{i\in\{1,2\}}\mathbb{P}(\mathcal{Q}\in E_{i})>0$, then almost surely

[TABLE]

Since $\mathcal{E}$ occurs with positive probability, this proves condition (AC).

4.3. Ginibre point process

The Ginibre point process is a determinantal point process with kernel

[TABLE]

with $z_{1},z_{2}\in\mathbb{C}$. As mentioned in [8, p. 19], this point process exhibits exponential decay of correlation. According to [21, Theorem 2], for $\boldsymbol{x}=(x_{1},\dots,x_{l})\in\mathbb{R}^{2l}$ we have $\mathbb{E}^{!}_{\boldsymbol{x}}\left[\mathcal{P}(W_{1})^{p}\right]\leq\mathbb{E}{\left[\mathcal{P}(W_{1})^{p}\right]}$, where the right-hand side is finite by [26, Lemma 4.2.6]. Hence, we obtain an upper bound for $\mathbb{E}^{!}_{\boldsymbol{x}}\left[\mathcal{P}(W_{1})^{p}\right]$, which is independent of $\boldsymbol{x}$, thereby verifying condition (M).

5. Simulation study

We elucidate in a simulation study, how cluster- and loop-based test statistics derived from Theorem 3.2 can detect deviations from complete spatial randomness. The simulations are carried out on top of the R-packages spatstat and TDA [20, 2].

For the entire simulation study, the null model $\mathsf{Poi}(2)$ is a Poisson point process with intensity $2$ in a $10\times 10$ observation window. Moreover, we fix $M=\sqrt{2}\cdot 10$ so large that it encompasses the entire sampling window and therefore suppress its appearance in the notation. Although the proof of Theorem 3.2 relies on the $M$-boundedness, the simulation study illustrates that it is not critical to impose this condition when testing hypotheses on common point patterns.

5.1. Deviation tests

As a first step, we derive scalar cluster- and loop-based test statistics.

5.1.1. Definition of test statistics

As a test statistic based on clusters, we use the integral over the number of cluster deaths in a time interval $[0,r_{\mathsf{C}}]$ with $r_{\mathsf{C}}\leq r_{\mathsf{f}}$, i.e.,

[TABLE]

After subtracting the mean, this test statistic becomes reminiscent of the classical Cramér-von-Mises statistic except that we do not consider squared deviations. Although squaring would make it easier to detect two-sided deviations, it would also require knowledge of quantiles of the square integral of a centered Gaussian process. Albeit possible, this incurs substantial computational expenses. Our simpler alternative has the appeal that as an integral of a Gaussian process, $T_{\mathsf{C}}$ is asymptotically normal and therefore characterized by its mean and variance.

As a test statistic based on loops, we use the accumulated persistence function, which aggregates the life times of all loops with birth times in a time interval $[0,r_{\mathsf{L}}]$ with $r_{\mathsf{L}}\leq r_{\mathsf{f}}$, i.e.,

[TABLE]

By Corollary 3.3, after centering and rescaling, the statistic $T_{\mathsf{L}}$ converges in the large-volume limit to a normal random variable.

The statistics $T_{\mathsf{C}}$ and $T_{\mathsf{L}}$ are specific possibilities to define scalar characteristics from the persistence diagram. Depending on the application context other choices, such as $\mathsf{APF}^{0}$ instead of $T_{\mathsf{C}}$ could be useful. However, in the simulation study below we found the weighting by life times of clusters to be detrimental.

5.1.2. Exploratory Analysis

As alternatives to the Poisson null hypothesis, we consider the attractive Matérn cluster and the repulsive Strauss processes. More precisely, the Matérn cluster process $\mathsf{MatC}(2,0.1,1)$ features a Poisson parent process with intensity 2 and generates a $\mathsf{Poi}(1)$ number of offspring uniformly in a disk of radius $0.1$ around each parent. The Strauss process $\mathsf{Str}(4.5,0.1,0.35)$ has interaction parameter $0.1$ and interaction radius $0.35$. The intensity parameter $4.5$ was tuned so as to match approximately the intensity of the null model. Figure 2 shows realizations of the null model and the alternatives.

In a first step, in Figure 3, we plot the persistence diagrams of samples from the null model and of the alternatives.

From the cluster-based diagrams, it becomes apparent that in comparison to the null model, in the Matérn cluster process, features can die also at rather late times, whereas this happens very rarely in the Strauss process. When analyzing loops, we see that loops with long life times can appear earlier in the null model than in the Matérn cluster process. Conversely, while some loops with substantial life time emerge at later times in the null model, there are very few such cases in the Strauss model.

5.1.3. Mean and variance under the null model

Now, we determine the mean and variance of $T_{\mathsf{C}}$ and $T_{\mathsf{L}}$ under the null model with $r_{\mathsf{f}}=1.5$. For this purpose, we compute the number of cluster deaths and accumulated loop life times for 10,000 independent draws of the null model.

Comparing the mean curves for the number of cluster deaths in the null model with those of the alternatives matches up nicely with the intuition about attraction and repulsion. For late times, they all approach a common value, namely the expected number of points in the observation window. However, Figure 4 shows that for the Matérn model, the slope is far steeper for early times, caused by merging of components of points within a cluster. In contrast, for the Strauss process the increase is at first much less pronounced than in the Poisson model, thereby reflecting the repulsive nature of the Gibbs potential.

For the loops, a radically different picture emerges. Here, the curve for the Strauss process lies above the accumulated loop life times of the null model. The Strauss model spawns substantially more loops than the Poisson model, although most of them live for a shorter period. Still, taken together these competing effects lead to a net increase of the accumulated loop life times in the Strauss model.

5.1.4. Type I and II errors

By Theorem 3.2, the statistics $T_{\mathsf{C}}$ and $T_{\mathsf{L}}$ are asymptotically normal, so that knowing the mean and variance allows us to construct a deviation test whose nominal confidence level is asymptotically exact. For the loops, we can choose the entire relevant time range, so that $r_{\mathsf{L}}=0.5$. For the cluster features, this choice would be unreasonable, as for late times, we simply obtain the number of points in the observation window, which is not discriminative. Hence, we set $r_{\mathsf{C}}=0.1$. We stress that in situations with no a priori knowledge of a good choice of $r_{\mathsf{C}}$, the test power can degrade substantially.

To analyze the type I and II errors, we draw 1,000 realizations from the null model and from the alternatives, respectively. Table 1 shows the rejection rates of this test setup. Under the null model the rejection rates are close to the nominal 5%-level, thereby illustrating that already for moderately large point patterns the approximation by the Gaussian limit is accurate. Using the mean and variance from the null model, we now compute the test powers for the alternatives. Already $T_{\mathsf{C}}$ leads to a test power of approximately $60\%$ for both alternatives. When considering $T_{\mathsf{L}}$, we obtain a type I error rate of 4.8%, so that the confidence level is kept. Moreover, the power analysis reveals that in the present simulation set-up, $T_{\mathsf{L}}$ is better in detecting deviations from the null hypothesis than $T_{\mathsf{C}}$.

5.2. Envelope Tests

Leveraging Theorem 3.2 shows that the deviation statistics $T_{\mathsf{C}}$ and $T_{\mathsf{L}}$ are asymptotically normal. Using a simulation-based estimate for the asymptotic mean and variance under the null model allowed us to construct a deviation test whose confidence level is asymptotically precise. A caveat of the above analysis is that the magnitude of clustering and repulsion is strong and clearly visible in the samples.

Recently, global envelope tests have gained widespread popularity, because they are both powerful and provide graphical insights as to why a null hypothesis is rejected [34]. The global envelope tests are fundamentally Monte Carlo-based tests and therefore do not relate directly to the large-volume CLT. However, they also rely on a functional summary statistics as input. Most of the applications in spatial statistics use a distance-based second-order functional such as Ripley’s $L$-function. In this section, we compare such classical choices with cluster- and loop-based statistics.

5.2.1. Alternatives

Since envelope tests excel at detecting subtle changes from the null model, we consider now a new parameter set-up to compare the $L$-function with the cluster- and loop-based statistics. Here, both the Matérn cluster as well as the Strauss process are substantially more similar to the Poisson point process. Hence, for the alternatives, we use again Matérn cluster and Strauss processes, but choose different parameters.

We found that the cluster- and loop-based statistics were particularly powerful in situations involving small interaction radii. Hence, as alternatives we choose the $\mathsf{MatC}(20,0.1,0.1)$ process and the $\mathsf{Str}(2.1,0.1,0.1)$ process, see Figure 5. The interaction parameter of the Strauss process was again tuned to match approximately the intensity of the null model.

5.2.2. Power analysis

To analyze the power of the envelope test, we generate $s=4,999$ realizations of the null model and 1,000 realizations of the alternatives. Then, we perform the global envelope test from [34] with three functional summary statistics. The first is Ripley’s $L$-function [33, Definition 4.6]. Second, we consider the number of cluster deaths as illustrated in Figure 4. Third, for the loops, we use a two-dimensional functional statistics derived from the persistent Betti numbers $\{\beta^{*,1}_{b,l}\}_{b,l}$ associated with life times $l$ rather than death times $d$ in order to expand the support of the statistic to the entire first quadrant. More precisely, $\beta^{*,1}_{b,l}$ counts the number of loops born before time $b$ and with life time at least $l$.

The rejection rates from Table 2 illustrate that for the alternatives described above, the cluster-based test gives a similar test power as the $L$-function-based test for the Matérn cluster process and a substantially higher test power for the Strauss process. Moreover, the loop-based test works even better in the Strauss case, but performs substantially worse for the Matérn alternative.

6. Analysis of the minicolumn dataset

In this section, we explore to what extent the deviation tests from Section 5 provide insights when dealing with real data. For this purpose, we analyze the minicolumn dataset provided by scientists at the Centre for Stochastic Geometry and Advanced Bioimaging.

As it should serve only to illustrate the application of Theorem 3.2, the present analysis is very limited in scope, and we refer to [14] for a far more encompassing study. For instance, that work considers two datasets and investigates 3D data together with marks for the directions attached to the neurons.

6.1. Exploratory analysis

The minicolumn dataset consists of 634 points emerging as two-dimensional projections of a three-dimensional point pattern of neurons. As neurons are believed to arrange in vertical columns, the projections are expected to exhibit clustering, see [32, 37]. The projections are taken along $z$-axis, since neuroscientists expect an arrangement in vertical columns. A visual inspection of the point pattern in Figure 6 supports this hypothesis.

As a first step, we explore whether the purported clustering already manifests in the persistence diagram. Comparing the loop-based persistence diagram of the minicolumn data with the persistence diagram of a homogeneous Poisson point process in Figure 7 shows that loops with substantial life times tend to be born later in the minicolumn model. This suggests clustering since loops formed by points within a cluster typically disappear rapidly.

Now, we explore whether the impressions from the persistence diagrams are reflected in the summary statistics from Section 5. When comparing in Figure 8 (left) the number of cluster death at different points in time, we note that until time 35, the curve for the observed data runs a bit above the curve for the null model. This provides already a first indication towards clustering. Next, we proceed to the loop-based features. As shown in Figure 8 (right), the curve for the observed pattern runs substantially below the one of the null model. This reflects a property that we have seen already in the persistence diagram: clusters with substantial life time tend to be born earlier in the null model, thereby leading to a steeper increase of the accumulated life times.

6.2. Test for complete spatial randomness

Under the impression of the previous visualizations, we now test the minicolumn pattern against the null model. As in Section 5, we deduce from Theorem 3.2 that the statistics are asymptotically normal under the null model, so that we only need to determine means and variances.

A subtle issue concerns the choice of the integration interval. The simplest option would be to take the whole intervals shown in Figure 8. For instance, for the loop-based features, this means $r_{\mathsf{L}}=r_{\mathsf{f}}=120$. However, for the cluster-based features the choice of the interval is less clear, since taking the whole interval is not discriminatory. The experiences from the simulation study indicate that the test is most powerful for early death times. Therefore, we choose $r_{\mathsf{C}}=10$.

With these choices, both the cluster-based and the loop-based test reject the null-hypothesis at the 5% level, since the corresponding $p$-values are $1.7\%$ and $1.2\%$. However, the tests are sensitive to the choice of the integration bound. This is especially true for the cluster-based test, where going from $r_{\mathsf{C}}=10$ to $r_{\mathsf{C}}=15$ results in a $p$-value of $5.2\%$, so that the null hypothesis is no longer rejected. On the other hand, the loop-based test still yields a $p$-value of $1.9\%$ $r_{\mathsf{L}}=75$. However, when reducing even further to $r_{\mathsf{L}}=50$, then the $p$-value increases sharply to $34.1\%$, so that the null-hypothesis is no longer rejected.

7. Discussion

In this paper, we elucidated how to apply tools from TDA to derive goodness-of-fit tests for planar point patterns. For this purpose, we derived sufficient conditions for a large-domain functional CLT for the $M$-bounded persistent Betti numbers on point processes exhibiting exponential decay of correlations. Following the framework developed in [8], the main difficulty arose from a detailed analysis of geometric configurations when bounding higher-order cumulants.

A simulation study revealed that the asymptotic Gaussianity is already accurate for patterns consisting of a few hundred data points. Additionally, as functional summary statistics, the persistent Betti numbers can also be used in the context of global envelope tests. Here, our finding is that TDA-based statistics can provide helpful additional information for point patterns with small interaction radii.

Finally, we applied the TDA-based tests on a point pattern from a neuroscientific dataset. As conjectured from the application context, the functional summary statistics indicate a clustering of points and the tests reject the Poisson null-model. However, the analysis also reveals a sensitivity to the range of birth times considered in the statistics.

In future work, we plan to extend the present analysis to dimensions larger than 2. On a technical level, the definition of higher-dimensional features requires a deeper understanding of persistent homology groups. Additionally, when thinking of broader application scenarios, a further step is to extend the testing framework from mere point patterns to random closed sets involving a richer geometric structure.

Acknowledgments

We thank J. Møller and A. D. Christoffersen for valuable discussions on the minicolumn dataset and for providing references. We also thank J. Yukich for inspiring discussions and helpful remarks. Finally, we thank the Centre for Stochastic Geometry and Advanced Bioimaging for collecting and sharing the data. CB, CH and AS are supported by The Danish Council for Independent Research — Natural Sciences, grant DFF – 7014-00074 Statistics for point processes in space and beyond, and by the Centre for Stochastic Geometry and Advanced Bioimaging, funded by grant 8721 from the Villum Foundation. NS was partially supported by the French ANR grant ASPAG (ANR-17-CE40- 0017).

8. Proof of Theorem 3.1

The main tool to prove Theorem 3.1 is the general CLT of [8, Theorem 1.14]. To be in that framework, we need to express the quantity \langle f,\mathsf{PD}^{M,q}(\mathcal{P}_{n})\rangle=\sum_{i\in I^{M,q}(\mathcal{P}_{n})}f\big{(}B_{i}^{M},D_{i}^{M}\big{)} in the form $\sum_{x\in\mathcal{P}_{n}}\xi(x,\mathcal{P}_{n})$ for a suitable score function $\xi(x,\mathcal{P}_{n})$.

In other words, we need to transform the indexing over features into an indexing over the points of the point process $\mathcal{P}_{n}$. We achieve this goal by assigning to each feature a point $x\in\mathcal{P}_{n}$ that either kills or gives birth to this feature, depending on whether $q=0$ or $q=1$.

First, the death of a cluster at time $r>0$ is always caused by the merging of two points $x,x^{\prime}\in\mathcal{P}_{n}$ at distance $2r$. Indeed, when the size of a component has a jump, this can only appear by attaching to another component. If $\mathcal{C}_{r}(x)$ dies by this merging, we say that $x^{\prime}$ kills $\mathcal{C}_{r}(x)$. This ensures that if two components both die when they merge, their deaths are caused by different points.

Similarly, if $q=1$, then the birth of a hole at time $r>0$ is caused by two points $x,x^{\prime}\in\mathcal{P}_{n}$ at distance $2r$ whose connection creates a new hole. If only one feature is born at time $r$, we choose the lexicographic minimum of $x$ and $x^{\prime}$ and say that it gives birth to this hole. However, if a large hole is split into two $M$-bounded holes, it can happen that two holes $H,H^{\prime}$ are born at the same time. In this situation, we assign one hole to each of $x$ and $x^{\prime}$. Hence, we define the score functions as

[TABLE]

Definition (8) translates the desired CLT for $\langle f,\mathsf{PD}^{M,q}(\mathcal{P}_{n})\rangle$ into the framework of [8, Theorem 1.14]. It remains to verify the conditions stated therein.

Proof of Theorem 3.1.

According to [8, Theorem 1.14] we have to verify that the pair $(\mathcal{P},\xi_{q})$ belongs to class (A2) (see Definition [8, Definition 1.7]) and that the $p$th moment condition [8, Equation (1.19)] holds for every $p$.

Belonging to class (A2) involves itself three conditions. The first is exponential decay of correlations, one of our standing assumptions on the point process $\mathcal{P}$. The second asks for an exponentially decaying radius of stabilization. Since we work with $M$-bounded features, this radius is finite. Finally, we need to verify the power-growth condition [8, Equation (1.18)] stating that for $W_{r}(x)=x+[-\sqrt{r}/2,\sqrt{r}/2]^{2}$, the upper bound

[TABLE]

holds for every $r>0$, locally finite $\mathcal{X}\subset\mathbb{R}^{2}$ and $x\in\mathcal{X}$. To achieve this goal, we note that in the worst case $x$ can be responsible for the death of all other points of $\mathcal{X}$. Similarly, it can give birth to at most $\mathcal{X}(W_{r}(x))-1$ holes. Hence,

[TABLE]

Finally, we verify the $p$th moment condition [8, Equation (1.19)]. That is, we prove that for every $p>0$ there exists $M_{p}>0$ such that

[TABLE]

We explain in detail how this is achieved if $q=0$, noting that the case $q=1$ can be deduced after minor modifications. If $x\in\mathcal{P}$ is responsible for the death of a component at time $r$, then there exists $x^{\prime}\in\mathcal{P}_{n}$ at distance $2r$ from $x$. Since each ball grows for time at most $r_{\mathsf{f}}$, we see that

[TABLE]

and an application of condition (M) concludes the proof. ∎

9. Proofs of Theorem 3.2 and Corollary 3.3

In the following, we assume $q=1$, since the proofs for $q=0$ are similar but easier. Hence, to simplify notation, we write $\beta_{b,d}(\mathcal{P}_{n})$ for $\beta^{M,1}_{b,d}(\mathcal{P}_{n})$.

Proof of Corollary 3.3.

Note that if $(X(s))_{s\leq r_{\mathsf{f}}}$ is a Gaussian process, then the process $(\int_{0}^{r}X(s){\rm d}s)_{r\leq r_{\mathsf{f}}}$ is also Gaussian. As mentioned above, the plan is to start from Theorem 3.2 and then apply the continuous mapping theorem [28, Theorem 4.27]. To this end, we show that $\{\mathsf{APF}^{M,1}_{r}(\mathcal{P}_{n})\}_{r\leq r_{\mathsf{f}}}$ is a continuous functional of the persistent Betti numbers $\{\beta_{b,d}(\mathcal{P}_{n})\}_{b,d\leq r_{\mathsf{f}}}$. We assert that

[TABLE]

The remainder of the proof proceeds in two steps. First, we verify identity (10). Second, we show that the right-hand side is continuous in $\beta$ with respect to the Skorokhod topology.

To prove identity (10), linearity allows us to reduce the claim to the case where the persistence diagram consists of a single $\delta$-measure at a point ${(B_{0},D_{0})}$ for some $D_{0}>B_{0}>0$. If $B_{0}>r$, then both sides vanish. If $B_{0}\leq r$, then $\beta_{b,0}=\mathbbmss{1}\{b\geq B_{0}\}$ and $\beta_{r,t}=\mathbbmss{1}\{t\leq D_{0}\}$, so that the right-hand side of (10) gives the asserted

[TABLE]

Let $\beta\in D([0,r_{\mathsf{f}}]^{2},\mathbb{R})$, where $D([0,r_{\mathsf{f}}]^{2},\mathbb{R})$ is the Skorokhod space of càdlàg functions from $[0,r_{\mathsf{f}}]^{2}$ to $\mathbb{R}$. For any $r\geq 0$ put

[TABLE]

According to (10), it is sufficient to prove that the function $\Phi_{r}:\,D([0,r_{\mathsf{f}}]^{2},\mathbb{R})\to D([0,r_{\mathsf{f}}],\mathbb{R})$, $\beta\mapsto(\Phi_{r}(\beta))_{r\leq r_{\mathsf{f}}}$ is continuous with respect to the Skorokhod topology. We prove this for the first integral. The arguments for the second are similar. Let $\beta^{\prime}:\,[0,r_{\mathsf{f}}]^{2}\to\mathbb{R}$ be càdlàg and $\lambda:[0,r_{\mathsf{f}}]\to[0,r_{\mathsf{f}}]$ be an increasing continuous bijection. Then,

[TABLE]

If $\beta^{\prime}$ approaches $\beta$ in the Skorokhod metric, then by definition of this metric, we can choose $\lambda$ such that the first two expressions become arbitrarily small. Moreover, since $\beta$ itself is càdlàg, it follows that also the third expression tends to 0. ∎

The proof of Theorem 3.2 decomposes into two steps: lower and upper variance bounds and an upper bound on fourth-order cumulants. In what follows, we write

[TABLE]

for the increment of $\beta_{b,d}$ in the block $E=(b_{-},b_{+}]\times(d_{-},d_{+}]$ with $b_{-}<b_{+}$ and $d_{-}<d_{+}$. Notice that this is minus the measure $\mathsf{PD}^{M,q}(\mathcal{P}_{n})$ from (2) evaluated at the block $E$. Moreover, $\beta(E,\mathcal{P}_{n})$ is the number of holes with birth time before $b_{+}$ and death time between $d_{-}$ and $d_{+}$ minus the number of holes with birth time before $b_{-}$ and death time between $d_{-}$ and $d_{+}$. Following [4], two blocks $E,E^{\prime}\subset[0,r_{\mathsf{f}}]^{2}$ are neighboring if they share a common side.

Proposition 9.1** (Variance lower bound).**

Let $\mathcal{P}$ be a conditionally $m$-dependent point process with exponential decay of correlations. Moreover, let $a_{1},\dots,a_{k}\neq 0$ and $E_{1},\dots,E_{k}\subset[0,r_{\mathsf{f}}]^{2}$ be pairwise disjoint blocks such that each $E_{i}$ contains some $(b,d)\in[0,r_{\mathsf{f}}]^{2}$ with $d>b$. Then,

[TABLE]

Proposition 9.2** (Variance upper bound).**

Let $\mathcal{P}$ be a conditionally $m$-dependent point process with exponential decay of correlations. Then, there exist $n_{0}\geq 1$ and $\varepsilon_{0},C_{0}>0$ such that

[TABLE]

holds for all $n\geq n_{0}$ and blocks $E\subset[0,r_{\mathsf{f}}]^{2}$.

Now, the $k$th cumulant $c^{k}$ of $k\geq 1$ real random variables $Y_{1},\dots,Y_{k}$ equals

[TABLE]

provided that all appearing moments are well-defined [36, Proposition 3.2.1]. Here, the sum ranges over all partitions $\{T_{1},\dots,T_{p}\}$ of the set $\{1,\dots,k\}$.

Proposition 9.3** (Cumulant bound).**

Let $\mathcal{P}$ be a conditionally $m$-dependent point process with exponential decay of correlations satisfying conditions (AC) and (M). Then, there exist $n_{0}^{\prime}\geq 1$ and $\varepsilon_{0}^{\prime},C_{0}^{\prime}>0$ such that

[TABLE]

holds for all $n\geq n_{0}^{\prime}$ and neighboring blocks $E,E^{\prime}\subset[0,r_{\mathsf{f}}]^{2}$.

We postpone the proofs of Propositions 9.1–9.3 to Sections 9.1–9.3, respectively. To deduce Theorem 3.2 from these two central auxiliary results, we write

[TABLE]

for the centered persistent Betti numbers.

Proof of Theorem 3.2.

Let $a_{1}^{\prime},\dots,a_{k^{\prime}}^{\prime}\neq 0$ and $(b_{1},d_{1}),\dots,(b_{k^{\prime}},d_{k^{\prime}})\in[0,r_{\mathsf{f}}]^{2}$ be pairwise distinct, and put

[TABLE]

Then, after suitable regrouping of terms, we can express $X_{n}$ in the form

[TABLE]

as in Proposition 9.1. Now, combining Proposition 9.2 with Theorem 3.1 and the variance asymptotics [8, Theorem 1.12] shows that the centered and rescaled random variable $n^{-1/2}(X_{n}-\mathbb{E}[X_{n}])$ converges in distribution to a Gaussian. Hence, the Cramér-Wold device yields convergence of the finite-dimensional distributions of $n^{-1/2}\overline{\beta}_{b,d}(\mathcal{P}_{n})$.

Next, [36, Proposition 3.2.1] gives the general cumulant identity

[TABLE]

for centered random variables $X,Y$. Hence, by Propositions 9.2 and 9.3,

[TABLE]

for some $\varepsilon_{0}^{\prime\prime}>0$. In particular, the process \big{\{}n^{-1/2}\overline{\beta}_{b,d}(\mathcal{P}_{n})\big{\}}_{b,d\leq r_{\mathsf{f}}} is tight in Skorokhod topology [24, Lemma 3]. In this context, we note that condition (8.4) of [24, Lemma 3] follows from the variance upper bound derived in Proposition 9.2 and that similar as in (2.18) [22], we have replaced the equality in (8.5) of [24, Lemma 3] by an inequality. Combining this property with the convergence of finite-dimensional distributions yields the asserted weak convergence. ∎

9.1. Proof of Proposition 9.1

To show the variance lower bound, we adapt a conditioning argument that has already been successfully applied in the setting of Gibbsian point processes [38]. More precisely, we subdivide the window $W_{n}$ into blocks of a fixed size and use the law of conditional variance to obtain a lower bound in the order of the number of blocks.

Associate with the $j$th feature $H_{j}$ in $\mathsf{PD}^{M,1}(\mathcal{P}_{n})$ a center point $y_{j}\in W_{n}$, for instance by taking the point $p(H_{j})$ as defined in Section 2.2. Then,

[TABLE]

defines a signed measure of total mass $\nu_{n}(W_{n})=\sum_{i\leq k}a_{i}\beta(E_{i},\mathcal{P}_{n})$.

In the vein of [38], the key towards proving a lower bound on the variance is the following non-degeneracy property, where $r_{\mathsf{AC}}$ is introduced in Section 3.

Lemma 9.4** (Non-degeneracy).**

It holds that

[TABLE]

Before proving Lemma 9.4, we explain how it implies Proposition 9.1. In essence, the proof follows along the lines of [38, Lemma 4.3]. Nevertheless, since the details of the conditioning argument differ a bit from the corresponding picture for Gibbs processes, we explain how to adapt the main steps from [38, Lemma 4.3] in the present setting.

Proof of Proposition 9.1.

The idea of proof is to consider a family of well-separated blocks in $W_{n}$. Then, we leverage the conditional $m$-dependence of the point process and the $M$-boundedness of the features to decompose the variance of their contributions as the sum of the variances. More precisely, we apply the assumption of conditional $m$-dependence with the conditioning set

[TABLE]

chosen as the complement of the union of well-separated blocks of side length $\rho=m\vee r_{\mathsf{AC}}$. Then, the law of total variance yields the lower bound

[TABLE]

Moreover, since $\rho>M$ the statistics $\nu_{n}((A^{\prime\prime})^{-})$ in the smaller domain

[TABLE]

is measurable with respect to $\mathcal{P}\cap{A^{\prime\prime}}$. We obtain that

[TABLE]

because $\nu_{n}((A^{\prime\prime})^{-})$ is $\mathcal{P}\cap A^{\prime\prime}$ measurable. Thanks to the conditional $m$-dependence, we have

[TABLE]

Now, the number of $6\rho$-blocks contained in $W_{n}$ is of order $n$, and we conclude by noting that Lemma 9.4 and $\rho>r_{\mathsf{AC}}$ imply that each of the contributions is bounded away from 0. ∎

To verify non-degeneracy, we rely on the techniques introduced in [38]. In particular, we make use of [38, Lemma 2.3], which we restate below to render the presentation self-contained.

Lemma 9.5**.**

Let $Y$ be a real random variable and $A_{1},A_{2}$ be Borel sets of $\mathbb{R}$. Then,

[TABLE]

Proof of Lemma 9.4.

Write

[TABLE]

for the events that there are no points in $W_{r_{\mathsf{AC}}^{2}/9}$ and $W_{r_{\mathsf{AC}}^{2}}\setminus W_{r_{\mathsf{AC}}^{2}/9}$, respectively. Next, let

[TABLE]

denote the event that all but the first of the considered persistent Betti numbers vanish. Now, let $I_{0}$ denote the indices of all features that are entirely contained in $\mathbb{R}^{2}\setminus W_{r_{\mathsf{AC}}^{2}}$ and put

[TABLE]

Then, by Lemma 9.5 with $A_{1}=[a_{1},\infty)$ and $A_{2}=\{0\}$,

[TABLE]

and it remains to show that the right-hand side is non-zero.

Since $E_{1},\dots,E_{k}$ are pairwise disjoint and contain points above the diagonal, [25, Example 1.8] shows that under the homogeneous Poisson point process the event $F^{\prime}\cap F_{2}$ has positive probability. Also $F^{\prime}\cap F_{1}$ is of positive probability. Hence, an application of condition (AC) concludes the proof. ∎

9.2. Proof of Proposition 9.2

For a block $E=(b_{-},b_{+}]\times(d_{-},d_{+}]\subset[0,r_{\mathsf{f}}]^{2}$, we let $\xi_{E}$ denote the score function associated with $\beta(E,\mathcal{P}_{n})$. That is,

[TABLE]

is the number of holes born by $x$ with birth and death times in $E$. Note that if $x$ gives birth to the $i$th hole, then it gets in contact with another point at time $B^{M}_{i}\in(b_{-},b_{+}]$. In particular, $\mathcal{P}$ contains a point in the annulus $A_{2b_{-},2b_{+}}(x)=B_{2b_{+}}(x)\setminus B_{2b_{-}}(x)$.

Moreover, if the $i$th hole dies at time $D^{M}_{i}\in(d_{-},d_{+}]$, then a previously vacant component is covered completely, which is caused by three disks centered at points in $\mathcal{P}$ meeting at a single point in the plane. The three center points of the disks must form a triangle with no obtuse angle. Otherwise, two of the disks would meet for the first time in the interior of the third and hence no connected component in the background was covered by the merging. This could be interpreted as a feature that is born and dies at the same time, but we chose to exclude such features in our definition of 1-features.

Henceforth, let $B^{\pm}_{d}(x,y)\subset\mathbb{R}^{2}$ denote the two disks of radius $d>0$ whose boundary passes through $x,y\in\mathbb{R}^{2}$. If $|x-y|/2>d$, we let $B^{\pm}_{d}(x,y)$ be empty. The points in $B^{+}_{d}(x,y)\cup B^{-}_{d}(x,y)$ are exactly the points $z$ such that the time when the boundaries of the three disks around $x$, $y$, and $z$ meet in one point is at most $d$. For $d_{+}>d_{-}\geq 0$, we let

[TABLE]

where $a=|x-y|/2$. This set consists of all points $z$ such that the boundaries of the three disks around $x$, $y$ and $z$ meet at time $r$ with $d_{-}<r\leq d_{+}$. Some $z\in D_{d_{-},d_{+}}(x,y)$ may still form a triangle having an obtuse angle with $x$ and $y$, that is, the disks around $x$, $y$, and $z$ already met earlier in an interior point of one of the disks. However, all $z$ that can cause the death of a hole in $E$ together with $x$ and $y$ must be contained in $D_{d_{-},d_{+}}(x,y)$.

Now,

[TABLE]

where $E_{x}$ denotes the event that for some $\mathcal{P}^{\prime}\subset\mathcal{P}$ with $x\in\mathcal{P}^{\prime}$ the event $E_{x,\mathsf{b}}(\mathcal{P}^{\prime})\cap E_{x,\mathsf{d}}(\mathcal{P}^{\prime})$ occurs, where

[TABLE]

Here, we say that $y_{1},y_{2},y_{3}\in\mathcal{P}^{\prime}$ kill the hole $H$ if the disks around the points meet for the first time at $p(H)$. In particular, any three points can kill at most one hole.

Similarly, for a block $E^{\prime}=(b_{-}^{\prime},b_{+}^{\prime}]\times(d_{-}^{\prime},d_{+}^{\prime}]\subset[0,r_{\mathsf{f}}]^{2}$,

[TABLE]

where we let $E_{x,x^{\prime}}^{\prime\prime}$ denote the event that for some $\mathcal{P}^{\prime}\subset\mathcal{P}$ with $x,x^{\prime}\in\mathcal{P}^{\prime}$ the event

[TABLE]

occurs.

Using this notation, the proof of the variance upper bound is now based on the following pivotal geometric moment bound. In the following, $\mathbb{P}_{\boldsymbol{x}}$ denotes the unreduced Palm measure characterized via $\mathbb{E}_{\boldsymbol{x}}[f(\mathcal{P})]=\mathbb{E}_{\boldsymbol{x}}^{!}[f(\mathcal{P}\cup\boldsymbol{x})]$ for any non-negative measurable $f:\,\mathcal{N}\to[0,\infty)$. We recall from (3) that $\rho^{(p)}$ denotes the $p$th factorial moment density. In the following, we adhere to the convention $\int_{B^{0}}f(x){\rm d}z=f(x)$.

Lemma 9.6** (Moment bound).**

Let $\mathcal{P}$ be a stationary point process having fast decay of correlations and satisfying condition (M). Let $p\geq 0$ and $K_{0}>0$. Then, there exist $\varepsilon>0$ and $C_{\mathsf{g}}>0$ such that for all $n>0$ and any ball $B\subset\mathbb{R}^{2}$ of radius $K>K_{0}$,

1. (1)

[TABLE]

holds for all blocks $E\subset[0,r_{\mathsf{f}}]^{2}$. 2. (2)

[TABLE]

holds for all neighboring blocks $E,E^{\prime}\subset[0,r_{\mathsf{f}}]^{2}$, and 3. (3)

[TABLE]

holds for all neighboring blocks $E,E^{\prime}\subset[0,r_{\mathsf{f}}]^{2}$.

The proof of Lemma 9.7 relies on a delicate geometric analysis that we defer to Section 9.4. We now prove Proposition 9.2. As in [8, Equation (1.6)], for $\boldsymbol{x}=(x_{1},\dots,x_{p})\in\mathbb{R}^{2p}$ and $k_{1},\dots,k_{p}\geq 0$, we introduce the mixed $\xi_{E}$-moments

[TABLE]

In the rest of the manuscript, we freely use that exponential decay of correlations implies boundedness of the factorial moment densities [8, Inequality (1.11)].

Proof of Proposition 9.2.

To lighten notation, we write $\xi$ instead of $\xi_{E}$. To give the paper more pleasant to read, we have not attempted to optimize the exponents occurring in the course of this proof. Proceeding as in [8, Equation (4.1)], the refined Campbell-Mecke formula [8, Equation (1.9)] gives that \mathsf{Var}\big{(}\beta(E,\mathcal{P}_{n})\big{)} equals

[TABLE]

We derive bounds for the two summands separately.

By stationarity, (11) and Hölder’s inequality, the first expression is at most

[TABLE]

Hence, Lemma 9.6(1) with $p=0$ yields the asserted upper bound.

To deal with the double integral in (13), we recall that $\xi$ is a local score function and that $\mathcal{P}$ exhibits exponential decay of correlations. Hence, as in [8, Equation (3.26)], the factorial moment measure expansion shows that

[TABLE]

for some $c>0$. In particular, choosing a cut-off $K=|E|^{-1/128}$, we see that

[TABLE]

holds for a suitable $C>0$ and it suffices to derive an upper bound for

[TABLE]

For the second summand, we can argue similarly as in (14), so that it remains to bound the integral involving $m_{n}^{(1,1)}(x,y)$. Here, we set $z=y-x$, note that $\rho^{(2)}(x,y)=\rho^{(2)}(o,z)$ and combine (11) with Hölder’s inequality to arrive at

[TABLE]

We bound $\mathbb{E}_{o,z}[\mathcal{P}(B_{M}(z))^{16}]$ thanks to condition (M). Finally, by Jensen’s inequality applied to the uniform distribution on $B_{K}(o)$,

[TABLE]

so that applying Lemma 9.6(1) with $p=1$ shows that the right-hand side is of order at most $(|E|^{3/4}|B_{K}(o)|)^{7/8}=|E|^{7\cdot 47/512}$, thereby concluding the proof. ∎

9.3. Proof of Proposition 9.3

To prove Proposition 9.3, we take up the idea suggested in [24, Theorem 8] and [10, Theorem 8.1] and express $c^{4}$ in terms of cumulant measures induced by the functional of interest. A slight technical nuisance in the present setting comes from dealing with a product of two different functionals – one associated with the block $E$ and the other with $E^{\prime}$ – whereas the semi-cluster measure machinery from [8, Section 4.3] relies on a single score function. However, this artificial difficulty can be overcome by formally attaching $\{1,2\}$-valued marks to $\mathcal{P}_{n}$. Taking up the notation from [19], we let $\breve{\mathbb{R}}^{2}=\mathbb{R}^{2}\times\{1,2\}$ and $\breve{\mathcal{P}_{n}}=\mathcal{P}_{n}\times\{1,2\}$ denote the correspondingly marked space and point process. Writing $E^{\prime\prime}=(E,E^{\prime})$, we define an augmented score function $\xi_{E^{\prime\prime}}$, where points with mark 1 are evaluated with the first score function and points with mark 2 are evaluated with the second score function. In other words,

[TABLE]

We take the concise proof of Proposition 9.2 as a blueprint for the strategy of the more involved setting laid out in Proposition 9.3. In particular, we need to address two main steps: bounds for mixed moments and a reduction of the integral to the diagonal.

In order to reduce to the diagonal, we decompose the cumulant measure into semi-cluster measures as in [3, Section 5.1] and [19, Section 3.2]. For the convenience of the reader, we reproduce the basic definitions. First, the $k$th moment measure $M^{k}(\mu_{n})$ is given as

[TABLE]

where $\boldsymbol{f}=f_{1}\otimes\cdots\otimes f_{k}$ is non-negative and measurable with each $f_{i}$ defined on $\breve{\mathbb{R}}^{2}$, and

[TABLE]

denotes the empirical measure associated with $\xi_{E^{\prime\prime}}$ and $\breve{\mathcal{P}_{n}}$. In terms of mixed $\xi$-moments, with $\breve{\boldsymbol{x}}_{T_{i}}$ the projection of $\breve{\boldsymbol{x}}$ to the coordinates in $T_{i}$, we write

[TABLE]

where ${\rm d}\breve{\boldsymbol{x}}_{T_{i}}$ are the singular differentials determined via

[TABLE]

where $f:\breve{\mathbb{R}}^{2|T_{i}|}\to[0,\infty)$ is any non-negative measurable function [19, Section 3.1]. As in (12), for $k_{1},\dots,k_{p}\geq 0$, the mixed $\xi_{E^{\prime\prime}}$-moments are given as

[TABLE]

for every $\breve{\boldsymbol{x}}=((x_{1},\tau_{1}),\dots,(x_{k},\tau_{k}))\in\breve{\mathbb{R}}^{2k}$.

Similarly, the $k$th cumulant measure $c^{k}_{n}=c^{k}(\mu_{n})$ equals

[TABLE]

so that

[TABLE]

where

[TABLE]

denotes the moment measure with coordinates in $T_{i}$.

Next, the space $\breve{W}_{n}^{4}$ decomposes into a union of subsets according to which coordinate is most distant from the diagonal [19, Lemma 3.1]. More precisely, write

[TABLE]

for the maximal separation of $\breve{\boldsymbol{x}}_{S}$ and $\breve{\boldsymbol{x}}_{T}$, where $\mathsf{dist}(\breve{\boldsymbol{x}}_{S},\breve{\boldsymbol{x}}_{T})=\mathsf{dist}(\boldsymbol{x}_{S},\boldsymbol{x}_{T})$. Then, put

[TABLE]

Here, the marks are ignored for the diagonal $\Delta\subset\breve{W}_{n}^{4}$. We also put $W^{(1,2)}_{n}=(W_{n}\times\{1\})^{2}\times(W_{n}\times\{2\})^{2}$.

Lemma 9.7** (Off-diagonal bounds).**

Let $S,T$ denote a non-trivial partition of $\{1,2,3,4\}$. Then, there exist $n_{S,T}\geq 1$ and $\varepsilon_{S,T},C_{S,T}>0$ such that

[TABLE]

holds for all $n\geq n_{S,T}$ and neighboring blocks $E,E^{\prime}\subset[0,r_{\mathsf{f}}]^{2}$.

Before proving Lemma 9.7, we elucidate how to deduce Proposition 9.3.

Proof of Proposition 9.3.

First, integration over the cumulant measure decomposes into a diagonal and an off-diagonal part [19, Equation (3.28)]. That is,

[TABLE]

where $\boldsymbol{f}=\mathbbmss{1}_{W^{(1,2)}_{n}}$ is the indicator function of the domain $W^{(1,2)}_{n}$ and the sum is over all non-trivial partitions $S,T$. By Lemma 9.7, the off-diagonal contributions in this decomposition are bounded above by $\sum_{S,T}C_{S,T}|E|^{1/2+\varepsilon_{S,T}}|E^{\prime}|^{1/2+\varepsilon_{S,T}}$.

Next, when integrating over the diagonal, we leverage that in the decomposition (17), only $p=1$ contributes [19, Lemma 3.1]. Hence,

[TABLE]

so that applying Lemma 9.6(2) with $p=1$ and noting the convention preceding that result concludes the proof. ∎

To prove Lemma 9.7, we decompose the cumulant measures into semi-cluster measures [3, Lemma 5.1]. More precisely, as in [3, 19], any two disjoint non-empty subsets $S^{\prime},T^{\prime}\preceq\{1,2,3,4\}$, induce a cluster measure

[TABLE]

Now, $c^{4}_{n}$ decomposes into semi-cluster measures

[TABLE]

where the sum runs over all partitions such that $S^{\prime}$ and $T^{\prime}$ are non-empty subsets of $S$ and $T$, respectively [19, Lemma 3.2].

Equipped with these ingredients, we now prove Lemma 9.7. Since the basic structure of the proof parallels that of Proposition 9.2, we only provide details for the steps that are substantially different.

Proof of Lemma 9.7.

Putting $D_{K}=\{\breve{\boldsymbol{x}}\in W^{(1,2)}_{n}\cap\sigma(S,T):\,D(\breve{\boldsymbol{x}})>K\}$ for $K\geq 1$, we first derive an upper bound for

[TABLE]

For this purpose, we decompose the moment measures ${\rm d}M^{S^{\prime}\cup T^{\prime}}$, ${\rm d}M^{S^{\prime}}$ and ${\rm d}M^{T^{\prime}}$ according to (16). Hence, we need bounds for the absolute value of differences of mixed $\xi$-moments of the form

[TABLE]

where $\{S^{\prime\prime}_{1},\dots,S^{\prime\prime}_{p^{\prime\prime}}\}$ and $\{T^{\prime\prime}_{1},\dots,T^{\prime\prime}_{r^{\prime\prime}}\}$ are partitions of $S^{\prime}$ and $T^{\prime}$, respectively. Since we are working on the set $\sigma(S,T)$, as in the proof of Proposition 9.2, the fast decay of $\xi$-correlations bounds (19) by $c\phi(D(\breve{\boldsymbol{x}}_{S^{\prime}\cup T^{\prime}})/2)$ for a suitable $c>0$.

Next, as in [19, Section 3.1] the singular differentials occurring in the expansion (16) of the moment measure $M^{k}$ can be grouped into a single object. More precisely, we write $\tilde{\rm d}\breve{\boldsymbol{x}}$ for the measure that equals ${\rm d}\breve{\boldsymbol{x}}_{T_{1}}\cdots{\rm d}\breve{\boldsymbol{x}}_{T_{p}}$ on the subset of $\breve{\mathbb{R}}^{2k}$ consisting of all $\breve{\boldsymbol{x}}=(\breve{x}_{1},\dots,\breve{x}_{k})$ such that $\breve{x}_{i}=\breve{x}_{j}$ if $i,j\in T_{r}$ for some $r\leq p$ and $\breve{x}_{i}\neq\breve{x}_{j}$ otherwise.

In the setting of the present proof, we note that the bounds on the mixed moments from (19) only involve coordinates with indices in the set $S^{\prime}\cup T^{\prime}$. Hence, we need to consider also singular differentials only with respect to these coordinates, i.e., integrate with respect to $\tilde{\rm d}\breve{\boldsymbol{x}}_{S^{\prime}\cup T^{\prime}}$. In particular, we arrive at the bound

[TABLE]

Now, setting $K=|E|^{-\varepsilon/128}|E^{\prime}|^{-\varepsilon/128}$, the exponential decay assumption on the function $\phi$ gives control on one integral over the window, while the integrals with respect to the remaining variables are controlled by the volume of balls. Then, a repeated application of Hölder’s inequality provides suitable bounds on the moment measures such that

[TABLE]

holds for some $C>0$. Hence, it suffices to provide upper bounds for

[TABLE]

where $\{T^{\prime}_{1},\dots,T^{\prime}_{p^{\prime}}\}$ is an arbitrary partition of $\{1,2,3,4\}$. We explain how to proceed for $p^{\prime}=1$, noting that for $p^{\prime}>1$ the arguments are similar but easier.

We claim that for some $C^{\prime}>0$,

[TABLE]

To prove this claim, decompose $M^{\{1,2,3,4\}}$ according to (16) and let $\{T^{\prime\prime}_{1},\dots,T^{\prime\prime}_{p^{\prime\prime}}\}$ be an arbitrary partition of $\{1,2,3,4\}$. As in the proof of Proposition 9.2, a repeated use of Hölder’s inequality shows that on $W^{(1,2)}_{n}$, the mixed moments of the form

[TABLE]

are bounded above by c^{\prime}\big{(}\mathbb{P}_{\boldsymbol{x}}(E_{x_{1},x_{i}}^{\prime\prime})\rho^{(p^{\prime\prime})}(\boldsymbol{x})\big{)}^{1-\varepsilon} for a suitable $c^{\prime}>0$ and some $i\leq 4$. At this point, we may proceed similarly as in (15) by invoking Lemmas 9.6(2) and 9.6(3). As an illustration consider the setting where $p^{\prime\prime}=4$ and $i=2$. Then, we set $z^{\prime}=x_{2}-x_{1}$, $z_{3}=x_{3}-x_{1}$ and $z_{4}=x_{4}-x_{1}$. We combine Jensen’s inequality with Lemma 9.6(3) to show that

[TABLE]

Hence, inserting the definition of $K$ concludes the proof. ∎

9.4. Proof of Lemma 9.6

We now turn to the proof of Lemma 9.6. The proof is based on the following four lemmas that are used to bound the probability with which certain point configurations occur. Throughout we use the notation

[TABLE]

The proofs make use of the inequalities

[TABLE]

where $C_{0}>0$ is some constant. Moreover, we repeatedly use that the volume of an annulus is given by

[TABLE]

Lemma 9.8**.**

Let $x,y\in\mathbb{R}^{2}$ and $a=|x-y|/2$. There is a constant $C>0$ such that for all $0\leq a\leq d_{+}\leq r_{\mathsf{f}}$,

[TABLE]

Proof.

Recall that

[TABLE]

The line through $x$ and $y$ cuts the disk $B_{d}^{+}(x,y)$ into two parts. The area of the larger part is given by

[TABLE]

$D_{d_{-},d_{+}}(x,y)$ is the union of two such sets of radius $d_{+}$ from which we remove two sets of the same type with radius $d_{-}\vee a$ from the interior. This yields the formula for the area.

The inequality follows from

[TABLE]

and, using (22) and (23),

[TABLE]

∎

Lemma 9.9**.**

Let $0\leq b_{-}<b_{+}\leq r_{\mathsf{f}}$ and $0\leq d_{-}<d_{+}\leq r_{\mathsf{f}}$ and let $B_{M}$ be a disk of radius $M$. Then, there is a constant $C>0$ such that

[TABLE]

Proof.

Integration with respect to $y_{3}$ yields:

[TABLE]

When $\delta_{b}\leq\delta_{d}$, the claim follows directly from Lemma 9.8. Otherwise, letting $a=|y_{1}-y_{2}|/2$, we split the integral in two terms according to whether $a<d_{-}$ or $a\geq d_{-}$. Applying Lemma 9.8 yields the bound

[TABLE]

To bound (24), we apply the mean value theorem and perform the integration to obtain the bound

[TABLE]

To bound (25), we bound the integrand using Lemma 9.8 and note that

[TABLE]

This proves the claim when $\delta_{d}\leq\delta_{b}$. ∎

Lemma 9.10**.**

Let $B_{M}$ be a disk of radius $M>0$. There is a constant $C>0$ such that for all $b_{-},b_{+},b_{-}^{\prime},b_{+}^{\prime},d_{-},d_{+}\in[0,r_{\mathsf{f}}]$ with $d_{-}<d_{+}$ and either $b_{-}<b_{+}=b_{-}^{\prime}<b_{+}^{\prime}$ or $b_{-}=b_{-}^{\prime}$ and $b_{+}=b_{+}^{\prime}$,

[TABLE]

Proof.

We may assume $d_{+}>3\delta_{d}\vee 8\sqrt{r_{\mathsf{f}}(\delta_{b}+\delta_{b^{\prime}})}$. Indeed, if $d_{+}\leq 3\delta_{d}$, we can show the claim by first integrating with respect to $y_{1}$, then using that by Lemma 9.9,

[TABLE]

and finally integrating with respect to $x_{2}$ and $x_{3}$ to provide a factor $|B_{M}|\delta_{b}\delta_{b^{\prime}}$. If $d_{+}\leq 8\sqrt{r_{\mathsf{f}}(\delta_{b}+\delta_{b^{\prime}})}$, we first integrate with respect to $x_{1}$, which yields the area of $A_{2b_{-},2b_{+}}(x_{2})\cap A_{2b_{-}^{\prime},2b_{+}^{\prime}}(x_{3})$. This is bounded by $C_{2}\delta_{b}\wedge\delta_{b^{\prime}}$, and by Lemma 9.9 the remaining integral is bounded by

[TABLE]

Let $a=|x_{2}-x_{3}|/2$. We write the integral as a sum of three terms corresponding to whether I: $a<d_{+}/4$, II: $d_{+}/4\leq a<b_{-}\wedge b_{-}^{\prime}$, or III: $b_{-}\wedge b_{-}^{\prime}\leq a\leq b_{+}\vee b_{+}^{\prime}$.

Term I: We first integrate with respect to $y_{1}$. Since

[TABLE]

the mean value theorem applied to the formula in Lemma 9.8 implies that $|D_{d_{-},d_{+}}(x_{2},x_{3})|\leq C_{4}\delta_{d}$. We then integrate with respect to $x_{2}$ and $x_{3}$ to obtain the bound $C_{5}|B_{M}|\delta_{b}\delta_{b^{\prime}}\delta_{d}$.

Term II: When $d_{+}/4\leq a\leq b_{-}\wedge b_{-}^{\prime}$, we first integrate with respect to $x_{1}$ to obtain the area of $A_{2b_{-},2b_{+}}(x_{2})\cap A_{2b_{-}^{\prime},2b_{+}^{\prime}}(x_{3})$. To bound term II, we need to explicitly compute this area. For this, we first compute the area $A_{a}(b_{1},b_{2})$ of the intersection $B_{2b_{1}}(x_{2})\cap B_{2b_{2}}(x_{3})$ where $b_{1},b_{2}\in\{b_{+},b_{-},b_{+}^{\prime},b_{-}^{\prime}\}$. By the assumption on $d_{+}$,

[TABLE]

This ensures that the line containing the two points where the boundaries of the disks $B_{2b_{1}}(x_{2})$ and $B_{2b_{2}}(x_{3})$ meet separates $x_{2}$ and $x_{3}$. The area of $B_{2b_{1}}(x_{2})\cap B_{2b_{2}}(x_{3})$ is

[TABLE]

The area of $A_{2b_{-},2b_{+}}(x_{2})\cap A_{2b_{-}^{\prime},2b_{+}^{\prime}}(x_{3})$ is given by

[TABLE]

It is a straightforward computation to see that $\tfrac{\partial^{2}}{\partial b_{1}\partial b_{2}}\mathcal{A}_{a}(b_{1},b_{2})$ is uniformly bounded by $C_{6}/d_{+}^{2}$ on the set of $a,b_{1},b_{2}\leq r_{\mathsf{f}}$ satisfying (26) and $d_{+}/4\leq a\leq b_{1}\wedge b_{2}$. In particular, (26) guarantees that

[TABLE]

such that $\arccos$ and $x\mapsto\sqrt{1-x^{2}}$ have bounded derivatives for the relevant values of $x$. It follows that (27) is bounded by $C_{7}\delta_{b}\delta_{b^{\prime}}/d_{+}^{2}$. The remaining integral is of order $|B_{M}|d_{+}^{2}\delta_{d}$ by Lemma 9.9, which yields the appropriate bound.

Term III: In this case, we first integrate with respect to $x_{1}$ providing a factor $\delta_{b}\wedge\delta_{b^{\prime}}$. The remaining integral is bounded using Lemma 9.9. ∎

The fourth lemma allows us to analyze which point configurations can cause the birth and death of $M$-bounded features. To state it, we recall the $\alpha$-complex associated with a locally finite point set $\mathcal{X}\subseteq\mathbb{R}^{2}$, see e.g. [18, Sec. III.4] for details. It is built from the Delaunay triangulation, which is a triangulation of the plane with vertex set $\mathcal{X}$. For $r>0$, $\alpha_{r}(\mathcal{X})$ is the union of all edges in the Delaunay triangulation with length at most $2r$ and all triangles such that the three balls of radius $r$ centered at its vertices cover the triangle. Then $\alpha_{r}(\mathcal{X})\subseteq U_{r}(\mathcal{X})$ and the inclusion is a homotopy equivalence, i.e. it preserves the topology.

Lemma 9.11**.**

Let $\mathcal{X}\subseteq\mathbb{R}^{2}$ be locally finite.

• (i)

Each connected component of $\mathbb{R}^{2}\backslash\alpha_{r}(\mathcal{X})$ contains at most one $M$-bounded connected component of $\mathbb{R}^{2}\backslash U_{r}(\mathcal{X})$.

• (ii)

If an $M$-bounded loop is born at time $b$ because two balls centered at $x_{1},x_{2}$ meet, then there is an edge of length $2b$ joining $x_{1},x_{2}$ in the $\alpha$-complex.

• (iii)

If an $M$-bounded feature dies at time $d$ because exactly three balls centered at points $y_{1},y_{2},y_{3}$ meet, then $y_{1},y_{2},y_{3}$ form a triangle with no obtuse angle in the $\alpha$-complex.

Proof.

The analogous statements hold for unbounded loops by the homotopy equivalence between the $\alpha$-complex and the union of balls. (i) follows because any $M$-bounded loop is also an unbounded loop. An $M$-bounded feature is either born the same way as the corresponding unbounded component or when two balls meet to split off a component. In both cases, some unbounded loop is born by the merging, and hence an edge is added to the $\alpha$-complex. This shows (ii). When an $M$-bounded loop dies, so does the corresponding unbounded loop, hence (iii) is clear. ∎

We are now ready to prove Lemma 9.6.

Proof of Lemma 9.6.

Proof of (1). Stationarity and Equation (4) yield

[TABLE]

In the following, we let $\boldsymbol{y}=(y_{1},y_{2},y_{3})$, and

[TABLE]

for simplicity. By definition of $E_{x}$, (LABEL:disint) is bounded by

[TABLE]

Here, we have used that $g(x_{1},x_{2},\boldsymbol{y})$ is symmetric in $x_{1}$ and $x_{2}$ and in $y_{1}$, $y_{2}$, and $y_{3}$. Applying (4) again, we may bound the last term in (LABEL:disjointSum) by

[TABLE]

since $b_{+}\leq M$. The remaining terms are treated similarly. Now choose a covering $B+x_{1}\subseteq\bigcup_{i\leq\ell}W_{1}^{(i)}$, where each $W_{1}^{(i)}$ is a translation of $W_{1}$ and such that $\ell\leq C_{1}|B|$ for some $C_{1}$ independent of $K$ (for instance using that $B_{K}\subseteq W_{4\lceil K\rceil^{2}}$). Then, by the moment condition (M) for $\boldsymbol{x}=(x_{1},\dots,x_{k})$,

[TABLE]

We apply this in (30) together with Lemma 9.9. Since each $\rho^{(k)}$ is bounded according to the assumption of fast decay of correlations, we obtain the bound $C_{4}|B|^{p+1}|E|^{1/2+\varepsilon}$.

**Proof of (2). ** In the following, we use the notation

[TABLE]

Note that since the blocks $E$ and $E^{\prime}$ are neighboring, the features in $E$ and $E^{\prime}$ are different. Putting $\boldsymbol{x}=(x_{1},x_{2},x_{3})$, we now expand as in (LABEL:disint)

[TABLE]

The condition $x_{2}\neq x_{3}$ comes from the fact that $x_{1}$ can give birth to at most one feature when connecting to another point, and since $E$ and $E^{\prime}$ are neighboring, $x_{2}$ and $x_{3}$ correspond to different features. Similarly, $\boldsymbol{y}^{\prime}\neq\boldsymbol{y}$ comes from the fact that a triangle can kill at most one feature.

The event $A$ excludes certain point configurations that are not possible. If the triangles formed by $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime}$ share an edge, and the vertices of this edge coincide with $x_{2}$ and $x_{3}$, then $|x_{2}-x_{3}|>2(b_{+}\vee b_{+}^{\prime})$ is not allowed. Indeed, it follows from Lemma 9.11 that the triangles correspond to the same feature in the $\alpha$-complex until $x_{2}$ and $x_{3}$ are joined. Thus, this must happen before both triangles are born, that is, at the latest at time $b_{+}\vee b_{+}^{\prime}$. Moreover, if the two triangles share an edge, then the two points in $\boldsymbol{y},\boldsymbol{y}^{\prime}$ not lying on this edge cannot be equal to $x_{1}$ and $x_{2}$ or to $x_{1}$ and $x_{3}$, as this would lead to crossing edges in the $\alpha$-complex by Lemma 9.11 (since the triangles formed by $\boldsymbol{y},\boldsymbol{y}^{\prime}$ cannot have any obtuse angles).

We now write the sum in (LABEL:9pointsum) as a sum where each term is a sum over $\mathcal{P}_{\neq}^{k}$, $4\leq k\leq 9$, as in (LABEL:disjointSum). Each such term comes from grouping $\boldsymbol{x},\boldsymbol{y},\boldsymbol{y}^{\prime}$ into sets of equal points. Consider for illustration the term corresponding to the situation $x_{2}=y_{1}^{\prime},x_{1}=y_{2}=y_{2}^{\prime},x_{3}=y_{3}=y_{3}^{\prime}$. The sum is handled as in the proof of Lemma 9.6(1) by applying (4) and bounding the involved Palm means. For this special point configuration, it is sufficient to bound $\mathbbmss{1}_{A}$ by 1.

[TABLE]

Now, we apply the Hölder inequality with $\frac{1}{q_{1}}+\frac{1}{q_{2}}=1$ to obtain the bound

[TABLE]

In the first integral, we first integrate with respect to $x_{2}$ and then apply Lemma 9.9, while in the second integral we first integrate with respect to $y_{1}$ and use the bound in Lemma 9.8 and then apply Lemma 9.9 again. Next we use that $E$ and $E^{\prime}$ are neighboring blocks so that either $\delta_{b}=\delta_{b^{\prime}}$ or $\delta_{d}=\delta_{d^{\prime}}$.

When $\delta_{b}=\delta_{b^{\prime}}$, we get the bound

[TABLE]

so we take $1/q_{1}>1/4$ and $1/q_{2}>2/3$.

When $\delta_{d}=\delta_{d^{\prime}}$, we use Lemma 9.9 to get the bound

[TABLE]

so we take $1/q_{1}>1/2$ and $1/q_{2}>1/3$.

For a general term, note that there are at least four different points among $\boldsymbol{y},\boldsymbol{y}^{\prime}$, so one of them, say $y_{1}$, cannot be equal to any of $\boldsymbol{x}$. We consider two cases:

• I

$y_{1}$ is not among $y_{1}^{\prime},y_{2}^{\prime},y_{3}^{\prime}$.

• II

$y_{1}=y_{1}^{\prime}$, $y_{2}=y_{2}^{\prime}$, and $y_{3}=x_{2}$ and $y_{3}^{\prime}=x_{3}$.

Since we no longer keep track of which edge kills which triangle, all possible point configurations allowed by $A$ fall into one of the above cases after possibly renaming the variables.

In particular, if $y_{1}=y_{1}^{\prime}$ and the points $y_{2},y_{3},y_{2}^{\prime},y_{3}^{\prime}$ are all different, one of them cannot be any of $x_{1},x_{2},x_{3}$, and we could have taken this as $y_{1}$ and be in Case I. If $y_{1}=y_{1}^{\prime}$, $y_{2}=y_{2}^{\prime}$ and, say, $y_{3}$ is not any of $x_{1}$, $x_{2}$, $x_{3}$, we could have chosen $y_{3}$ as $y_{1}$ and be in Case I.

We further divide the Case I configurations allowed by $A$ into the following two sub-cases that have to be treated separately:

• Ia.

$x_{3}$ is not any of $y_{2},y_{3}$.

• Ib.

$x_{2}=y_{2}=y_{2}^{\prime}$, $x_{3}=y_{3}=y_{3}^{\prime}$, $|x_{2}-x_{3}|/2\leq b_{+}\vee b_{+}^{\prime}$.

Again, after renaming the variables, we are always in one of the two sub-cases.

Case Ia: We apply the Hölder inequality to

[TABLE]

The first factor is integrated with respect to $x_{3}$ and the remaining integral is bounded using Lemma 9.9. The second factor is first integrated wrt. $y_{1}$, the result is bounded using Lemma 9.8, and the remaining integral is bounded using Lemma 9.9. The rest of the argument proceeds as in the special case treated above.

Case Ib: The claim follows by applying the Hölder inequality to (35) and arguing as in Case Ia using Lemma 9.10 to bound the first integral.

Case II: We apply the Hölder inequality exactly as in (35) and argue as in Case Ia, except that the second integral is first integrated with respect to $y_{3}$ rather than $y_{1}$.

Proof of (3). As in (LABEL:disint), we find

[TABLE]

The set $\tilde{A}$ consists of tuples of points $(x_{1},x_{2},x_{3},x_{4},\boldsymbol{y},\boldsymbol{y}^{\prime})\in\mathbb{R}^{20}$ and, similar to $A$, it excludes certain configurations of the points $(x_{1},x_{2},x_{3},x_{4},\boldsymbol{y},\boldsymbol{y}^{\prime})$ that are not allowed by Lemma 9.11. If the triangles formed by $\boldsymbol{y}$ and $\boldsymbol{y}^{\prime}$ share an edge, then the length of this edge must be at most $2(b_{+}\vee b_{+}^{\prime})$. Moreover, if the two triangles share an edge, then the two points in $\boldsymbol{y},\boldsymbol{y}^{\prime}$ not lying on this edge cannot be equal to $x_{1}$ and $x_{3}$ or to $x_{2}$ and $x_{4}$.

The contribution from the cases where two of the points $x_{1},x_{2},x_{2}^{\prime},z^{\prime}$ are identical is bounded by

[TABLE]

which is handled exactly as in the proof of Lemma 9.6(2). Thus, it remains to treat the terms where $x_{1},x_{2},x_{2}^{\prime},z^{\prime}$ are all different. Therefore, if we put $\boldsymbol{x}=(x_{1},x_{2},x_{3},x_{4})$, we must bound

[TABLE]

The rest of the proof proceeds as the proof of Lemma 9.6(2) by suitable applications of the Hölder inequality. We divide into two cases according to whether all points in $\boldsymbol{y}$, $\boldsymbol{y}^{\prime}$ are one of $\boldsymbol{x}$ or not. After renaming the variables, we may assume

• I

$y_{1}=y_{1}^{\prime}=x_{1}$, $y_{2}=y_{2}^{\prime}=x_{2}$, $y_{3}=x_{3}$, and $y_{3}^{\prime}=x_{4}$, or

• II

$y_{1}$ is not any of $\boldsymbol{x}$.

Notice that in Case I we exclude the case $y_{1}=y_{1}^{\prime}=x_{1}$, $y_{2}=y_{2}^{\prime}=x_{3}$, $y_{3}=x_{2}$, and $y_{3}^{\prime}=x_{4}$ because it was excluded by definition of $\tilde{A}$. After renaming variables, Case II is divided into

• IIa

$y_{1}$ is not any of $\boldsymbol{x}$ or $\boldsymbol{y}^{\prime}$, and $x_{1}$ is not any of $y_{2},y_{3}$.

• IIb

$y_{1}=y_{1}^{\prime}$ and $y_{1}$ is not any of $\boldsymbol{x}$, $y_{2}=y_{2}^{\prime}\neq x_{3}$, $y_{3}=x_{1}$.

• IIc

$y_{1}=y_{1}^{\prime}$, $y_{2}=x_{2}$, $y_{3}=x_{4}$, $y_{2}^{\prime}=x_{1}$, $y_{3}^{\prime}=x_{3}$.

• IId

$y_{1}=y_{1}^{\prime}$, $y_{2}=x_{1}$, $y_{3}=x_{2}$, $y_{2}^{\prime}=x_{3}$, $y_{3}^{\prime}=x_{4}$.

In Case IIa, $y_{1}$ is not one of $\boldsymbol{y}^{\prime}$, while in Case IIb, IIc, and IId it is. Case IIb corresponds to the situation in which the triangles formed by $\boldsymbol{y},\boldsymbol{y}^{\prime}$ share an edge, while in Case IIc and IId they share only one vertex. In Case IIc, each triangle contains one of the edges joining $x_{1}$ to $x_{3}$ and $x_{2}$ to $x_{4}$, while in Case IId they do not.

Case I: When $\delta_{b}=\delta_{b^{\prime}}$, we first write

[TABLE]

We then apply the Hölder inequality. Integrating first with respect to $x_{4}$ and then $y_{2}$ in (36) and integrating with respect to $y_{3}$ first in (37) yields a bound of order

[TABLE]

This is the same as (33) since $\delta_{b}=\delta_{b^{\prime}}$. When $\delta_{d}=\delta_{d^{\prime}}$, we replace $\mathbbmss{1}_{D_{d_{-},d_{+}}(y_{1},y_{2})}(y_{3})$ by $\mathbbmss{1}_{D_{d_{-}^{\prime},d_{+}^{\prime}}(y_{1}^{\prime},y_{2}^{\prime})}(y_{3}^{\prime})$ in (36), to obtain a bound of order

[TABLE]

which reduces to the same form as (34).

Case IIa: We apply the Hölder inequality to (36)–(37) and integrate first with respect to $x_{1}$ and then $y_{1}$ in (36) and with respect to $y_{1}$ first in (37). The remaining argument proceeds as in the proof of Lemma 9.6(2) Ia.

Case IIb: We apply the Hölder inequality to (36)–(37) and integrate first with respect to $x_{3}$ and then $y_{1}$ in (36) and with respect to $y_{3}$ first in (37) and argue as in the proof of Lemma 9.6(2) Ia.

Case IIc: In (36), we first integrate with respect to $x_{1}$. In (37), we first integrate with respect to $y_{2}^{\prime}$ and $y_{3}^{\prime}$ to obtain a factor $\delta_{d^{\prime}}$. Then we integrate with respect to $y_{1}$ and $x_{2}$ and apply Lemma 9.9 to obtain a factor $\delta_{d}^{1/2}\delta_{b^{\prime}}$. The resulting bounds are stricter than (33) and (34).

Case IId: Here we integrate (36) with respect to $x_{3}$ first and then $y_{1}$ while (37) is integrated first with respect to $y_{2}$ and then $y_{1}$.

In all cases treated above, a minor difference to (LABEL:holderStep) is that the integration domains are slightly more complicated due to the indicator $\mathbbmss{1}_{B+x_{1}}(x_{2})$. However, it contributes at most a factor $C_{7}|B|$ to the bound, and this cancels when we divide by $|B|^{p+2}$.

∎