# Minimax estimation of the structure factor of spatial point processes

Gabriel Mastrilli

► To cite this version:

Gabriel Mastrilli. Minimax estimation of the structure factor of spatial point processes. 2025. hal-05520136

HAL Id: hal-05520136

https://hal.science/hal-05520136v1

Preprint submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Minimax estimation of the structure factor of spatial point processes

Gabriel Mastrilli¹,²

¹ Univ Rennes, Ensai, CNRS, CREST

² Inria, Paris, France.

## Abstract

We investigate the problem of estimating the structure factor, or spectra, of stationary spatial point processes. In the first part, we establish a minimax lower bound for this estimation problem, using an approach tailored to second-order properties of spatial point processes. Although not the main focus, this methodology also extends naturally to a minimax lower bound for the estimation of the pair correlation function of spatial point processes. In the second part, we construct a multitaper estimator that achieves the optimal rate of convergence in squared risk. Under a Brillinger-mixing condition, we further show that its expected chi-square-type concentration for a fixed number of tapers remains valid when this number diverges. Finally, we propose a data-driven procedure for selecting the number of tapers, supported by an oracle inequality, and we demonstrate the practical effectiveness of the method through numerical experiments.

MSC Classification: Primary, 62M15, 62M30, 62G05; Secondary, 62C20, 60G55.

Keywords: Cross-validation, Deviation inequality, Minimax, Multitaper, Pair-correlation function.

## 1 Introduction

Spectral methods for second-order inference in spatial point processes were introduced in the seminal works of Bartlett [5, 6]. These methods aim to estimate the spatial analogue of the spectral density of a time series, typically Bartlett's spectral measure, often referred to simply as the *spectra*, or its normalized version, the *structure factor*. While spectral density estimation for time series has been a central topic of statistical research for nearly a century (see, e.g., [62] and references therein), its spatial counterpart remains comparatively underdeveloped, despite its relevance in many applications. In particular, the structure factor plays a crucial role in statistical physics, where it underpins the analysis of crystalline and disordered materials [23, 48, 76]. In spatial statistics, following Bartlett's foundational contributions, the works [58, 57, 34, 59, 65, 66, 42] highlighted the practical value of spectral methods for analyzing spatial point patterns. However, several theoretical statistical questions had long remained only partially addressed. In the past few years, however, this line of research has received renewed and sustained attention, notably through a series of works including [64, 36, 31, 83, 24], which have significantly advanced the theoretical understanding of spectral estimators for spatial point processes. Despite these contributions, the question of optimal
rates of convergence for structure factor estimation remains unresolved. This paper addresses this gap by establishing minimax lower bounds and constructs an estimator that achieves these optimal rates.

Specifically, we establish that the minimax optimal rate of convergence for estimating a structure factor with Hölder regularity $\beta \in (0, \infty)$, based on a single realization of a stationary point process observed within a compact window $W \subset \mathbb{R}^d$, is of order $|W|^{-\beta/(2\beta+d)}$ (see Theorem 4.1). Equiva- lently, this rate can be expressed in terms of the expected number of observed points. Identifying $|W|$ with the sample size, this convergence rate formally matches the classical nonparametric rate for estimating Hölder functions [77], and coincides with the minimax rates known for spectral density estimation in time series [70, 9, 27, 1, 45]. As in these settings, the smoothness parameter $\beta$ has a direct interpretation in terms of second-order correlations: larger values of $\beta$, and hence faster rates, correspond to processes with more rapidly decaying pairwise correlations.

Nevertheless, although the overall structure of the proof of the minimax lower bound follows standard minimax arguments from nonparametric estimation theory, its implementation requires significant adaptation due to the different nature of the objects involved. In the time series setting, minimax lower bounds for spectral density estimation typically rely on Gaussian processes. Such Gaussian structure is not readily available for point processes. Moreover, while minimax lower bounds have been established for first-order properties of non-stationary Poisson point processes [67, 49, 68], these results yield trivial bounds when applied to second-order estimation in the stationary case, since the structure factor of any stationary Poisson is constant equal to one. To overcome this, and to leverage the analytical tractability of both Gaussian fields and Poisson point processes, we consider a suitable class of Cox point processes [56], where the random intensity field is given by $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$, with $(N(x))_{x \in \mathbb{R}^d}$ a centered stationary Gaussian field. For this class, we verify that standard minimax techniques apply (see [77, Chapter 2]). In addition, because this construction is specifically chosen to facilitate minimax arguments for second-order statistics of spatial point processes, it also leads to a minimax lower bound for the estimation of the pair correlation function (see Theorem 4.2). Although this result lies outside the main focus of the present work, we believe it is of independent interest for the spatial statistics community and worth reporting.

The second part of this work is devoted to the construction of an estimator of the structure factor that attains the optimal minimax rate of convergence. This estimator relies on the principle of multitapering [75], which we motivate and introduce in what follows. Earlier contributions to spectral inference have emphasized the benefits of refining the basic periodogram [78, 12, 5, 17]. A common refinement is to apply a taper to the periodogram to reduce bias. Nevertheless, the variance remains high. To mitigate this issue, one standard approach is to average multiple single-taper estimators. Two averaging schemes have been considered in the spatial statistics literature, inspired by analogous methods in time series analysis [62, Chapters 7–8]. The first scheme averages single-taper estimators evaluated at distinct frequencies near the target frequency, using a fixed taper function [64, 83, 24]. The second scheme, which we adopt in this work, averages tapered estimators evaluated at a fixed frequency across a family of orthogonal taper functions [64, 31]. This approach, called multitapering, is well-established in signal processing [75, 81, 44] and has recently proven effective for spectral analysis of spatial point processes [36, 64, 31, 53]. Importantly, it enables estimation of the structure factor at a fixed frequency without requiring information from nearby frequencies. Moreover, with a fixed number of tapers, it typically has a simple asymptotic $\chi^2$-type distribution[6, 75, 62, 64]; see equation (10).

To study the theoretical performance of multitaper estimators and prove that they can attain the minimax rate of convergence, we first establish a non-asymptotic upper bound on their squared risk (Theorem 5.1). This result makes explicit how the bias–variance trade-off is governed by the number of tapers and their spatial and spectral localizations. This insight motivates the use of
Hermite functions, which are simultaneously well localized in both domains and numerically easy to compute. We show that this choice leads to an estimator achieving the minimax-optimal convergence rate for structure factors with Hölder regularity $0 < \beta \le 2$. In addition, under Brillinger-mixing assumptions, we strengthen our squared risk guarantee by deriving a non-asymptotic high-probability bound. Specifically, we show that the usual $\chi^2$-type concentration of the multitaper estimator remains valid even when the number of tapers depends on the observation window (Theorem 5.2).

Finally, we propose a cross-validation procedure to select the number of tapers. This selection is local, allowing for different numbers of tapers across frequencies. The method adapts “leave-one-out” techniques, originally developed for bandwidth selection in kernel-smoothed spectral estimation [8, 41, 79, 24], to the multitaper framework for spatial point processes. We provide theoretical guarantees for this procedure by establishing an oracle inequality (see Theorem 6.1), and illustrate its practical effectiveness through numerical experiments, whose code is available at https://github.com/gabrielmastrilli/multitaper_spp.

The paper is organized as follows. Section 2 gathers our general notations. In Section 3, we introduce the notions that are specific to stationary point processes, in order to define the structure factor. Section 4 states minimax lower bounds for the problems of estimating the structure factor and the pair correlation function of a stationary point process. The multitaper estimator considered in this work is introduced in Section 5.1. In this section, we derive heuristic arguments supporting it and highlighting the necessity to consider a non-asymptotic framework to discuss the choice of its parameters. Section 5.2 studies its $L^2(\mathbb{P})$-risk and proves that for Hermite tapers, the multitaper estimator is optimal for the minimax criteria. Section 5.3 studies the concentration properties of the multitaper estimator and shows that the results of Section 5.2 remain valid with high probability. Finally, Section 6 presents the data-driven procedure for selecting the number of tapers and investigates its numerical performance. The proofs of the minimax lower bounds are provided in Section 7.

Appendices A, B, and C contain the proofs corresponding to Sections 5.2, 5.3, and 6, respectively. These proofs rely on technical lemmas presented in Appendix D concerning the cumulants of linear statistics of point processes. Appendix E concerns a structure factor estimator in which the unknown intensity is replaced by its empirical estimate. Appendix F gathers results on Hermite functions, while Appendix G collects technical tools on cumulants of random variables.

## 2 General notations

We consider the Euclidean space $\mathbb{R}^d$ of dimension $d \ge 1$. For $a, b \in \mathbb{C}^d$, $a \cdot b = \sum_{i=1}^d a_i \bar{b}_i$ denotes the Hermitian scalar product of $\mathbb{C}^d$ between $a$ and $b$. The associated Euclidean norm is denoted by $|a|$ and the supremum norm is $|a|_{\infty} := \max_{i=1,...,d} |a_i|$. For any reals $x$ and $y$, we use the notations $x \vee y := \max(x, y)$ and $x \wedge y := \min(x, y)$. For $x \in \mathbb{R}$, $\lfloor x \rfloor$ denotes the greatest integer less than $x$. For a discrete set $X$, $|X|$ denotes its cardinal.

For $1 \le p < \infty$, we denote by $L^p(\mathbb{R}^d)$ the space of measurable functions $f : \mathbb{R}^d \to \mathbb{C}$ such that $\int_{\mathbb{R}^d} |f(x)|^p dx < \infty$. The $L^p$-norm of $f \in L^p(\mathbb{R}^d)$ is defined by $\|f\|_p := \int_{\mathbb{R}^d} |f(x)|^p dx$. We denote by $L^\infty(\mathbb{R}^d)$ the space of functions $f : \mathbb{R}^d \to \mathbb{C}$ such that $\|f\|_\infty := \sup_{x \in \mathbb{R}^d} |f(x)| < \infty$. We adopt the following definition for the Fourier transform of a function $f \in L^1(\mathbb{R}^d)$:

$$ \forall k \in \mathbb{R}^d, \quad \mathcal{F}[f](k) := \frac{1}{(2\pi)^{d/2}} \int_{\mathbb{R}^d} f(x)e^{ik \cdot x} dx. $$

As usual, the Fourier transform is extended to $L^2(\mathbb{R}^d)$ functions thanks to the Plancherel theorem [29]: $\forall f_1, f_2 \in L^2(\mathbb{R}^d)$, $\langle f_1, f_2 \rangle = \langle \mathcal{F}[f_1], \mathcal{F}[f_2] \rangle$. The homogeneous Sobolev space $\dot{H}^s(\mathbb{R}^d)$ is the
space of functions $f \in L^2(\mathbb{R}^d)$ with $\|f\|_{\dot{H}^s}^2 := \int_{\mathbb{R}^d} |\mathcal{F}[f](k)|^2 |k|^{2s} dk < \infty$. We denote by $\delta_0$ the Dirac mass at zero.

We also denote, for $n \ge 1$, $[n] = \{1, \dots, n\}$ and $\Pi[n]$ the set of partitions of $[n]$. For a partition $\sigma \in \Pi[n]$, $|\sigma|$ denotes its number of blocks. For two partitions $\sigma$ and $\tau$ of $[n]$, where $n \ge 1$, their least upper bound (refer to Chapter 3 of [54] or Chapter 2 of [61]) is denoted $\sigma \vee \tau$. We define algebraically the cumulants. For random variables $X_1, \dots, X_m$, we denote their $m$-th cumulants as:

$$ \kappa(X_1, \dots, X_m) := \sum_{\sigma \in \Pi[n]} (-1)^{|\sigma|-1} (|\sigma|-1)! \prod_{b \in \sigma} \mathbb{E} \left[ \prod_{s \in b} X_s \right], \quad (1) $$

where $\prod_{b \in \sigma}$ denotes the product over all the blocks of the partition $\sigma$. Moreover, for a random variable $X$ having moments of order $m$, its $m$-th cumulants are defined as $\kappa_m(X) = \kappa(X, \dots, X)$, where $X$ appears $m$-times.

Finally, the total variation of a signed measure $\nu$ is defined by $|\nu| := \nu^+ + \nu^-$, where $\nu = \nu^+ - \nu^-$ is the Jordan decomposition of $\nu$ (see, e.g., [26], Theorem 5.6.1).

# 3 Notations and definitions specific to point processes

In the following, we introduce the key concepts and results concerning the structure factor in Section 3.1. Then cumulants of point processes are introduced in Section 3.2. For an in-depth treatment of point processes, we refer readers to the classic two-volume work [20, 21] or to the more recent, concise exposition in [4] and for fundamental understanding of the power spectra of point processes, we recommend [14].

## 3.1 Second order properties and structure factor of stationary point processes

The set of points configurations in $\mathbb{R}^d$ is defined as:

$$ \mathrm{Conf}(\mathbb{R}^d) := \{\phi \subset \mathbb{R}^d \mid \text{For all } K \text{ compact of } \mathbb{R}^d, \text{ then } |\phi \cap K| < \infty\}. $$

This set is endowed with the $\sigma$-algebra generated by the mapping $\phi \mapsto |\phi \cap K|$ for all compact sets $K$. A point process $\Phi$ is a random element of $\mathrm{Conf}(\mathbb{R}^d)$. A point process $\Phi$ is called stationary if for all $x \in \mathbb{R}^d$, $\Phi + x := \{y + x \mid y \in \Phi\}$ is equal in distribution to $\Phi$. As a consequence, the intensity measure $\rho^{(1)}$ of a stationary point process $\Phi$ (defined for any subset $A$ of $\mathbb{R}^d$ by $\rho^{(1)}(A) = \mathbb{E}(|\Phi \cap A|)$), is proportional to the Lebesgue measure on $\mathbb{R}^d$: $\rho^{(1)} = \lambda dx$. The scalar $\lambda \ge 0$ is called the intensity of the point process. The second order factorial moment measure $\rho^{(2)}$ of a point process $\Phi$ is a measure on $(\mathbb{R}^d)^2$ defined by

$$ \rho^{(2)}(A_1 \times A_2) = \mathbb{E} \left[ \sum_{x,y \in \Phi} \mathbf{1}_{x \in A_1, y \in A_2} \right], \quad (2) $$

for all $A_1, A_2$ subsets of $\mathbb{R}^d$. The symbol $\neq$ over the sum means that we consider only distinct points. To define the structure factor of a point process and ensure that it is a function and not a general measure, we introduce the following assumption.
**Assumption 3.1.** Throughout the paper, we tacitly assume that the point process $\Phi$ is stationary, with an intensity $\lambda > 0$, and that its second-order intensity measure $\rho^{(2)}$ is absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}^d \times \mathbb{R}^d$. This allows us to represent it as follows:

$$ \rho^{(2)}(dx, dy) = \lambda^2 g(x - y)dxdy, $$

where $g(x)$ is a function on $\mathbb{R}^d$. Additionally, we assume that $g - 1 \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d)$.

The function $g: \mathbb{R}^d \to \mathbb{R}$ is known as the pair-correlation function of $\Phi$. With this established, we can proceed to define the structure factor of $\Phi$ as

$$ \forall k \in \mathbb{R}^d, S(k) := 1 + \lambda \int_{\mathbb{R}^d} (g(x) - 1)e^{-ik \cdot x} dx. \quad (3) $$

The previous definition follows the convention from the physics literature [18, 35, 76]. Indeed, $S$ is a standardized version of $\lambda S$, which is the Bartlett spectral measure of $\Phi$ [6, 21], called spectra in related statistical works [64, 31, 83]. In the following, we assume $\lambda = 1$, so that the structure factor and the spectra are equal. This is not restrictive in practice, when $\lambda$ is unknown and has to be estimated, as shown in Appendix E.

By application of the Campbell formula [4] and of the Plancherel Theorem [29], we deduce the following representation of the covariance between two linear statistics of $\Phi$. For all $f_1, f_2 \in L^1(\mathbb{R}^d) \cap L^2(\mathbb{R}^d)$, we have

$$ \mathrm{Cov} \left[ \sum_{x \in \Phi} f_1(x), \sum_{x \in \Phi} f_2(x) \right] = \lambda \int_{\mathbb{R}^d} \mathcal{F}[f_1](k) \overline{\mathcal{F}[f_2]}(k) S(k) dk. \quad (4) $$

Moreover, standard results concerning the Fourier transform of $L^1(\mathbb{R}^d)$ functions [29] imply that, under Assumption 3.1, then $S$ is bounded and continuous.

## 3.2 Cumulants of stationary point processes

In this section we recall the definition of the factorial cumulant measures of a simple point process and their reduced counterpart in the stationary case; see [20, 21].

As a generalization of the second order factorial moment measure defined in (2), the factorial moment measures of order $m \ge 1$, $\rho^{(m)}$ of a simple point process $\Phi$ is defined for any bounded measurable sets $A_1, \dots, A_m$ in $\mathbb{R}^d$ by

$$ \rho^{(m)}(A_1 \times \dots \times A_m) = \mathbb{E} \left( \sum_{x_1, \dots, x_m \in \Phi} \mathbb{1}\{x_1 \in A_1, \dots, x_m \in A_m\} \right). $$

Then, the $m$-th order factorial cumulant moment measure $\gamma^{(m)}$ is defined by:

$$ \gamma^{(m)}(A_1 \times \dots \times A_m) = \sum_{\sigma \in \Pi[m]} (-1)^{|\sigma|-1} (|\sigma|-1)! \prod_{b \in \sigma} \rho^{|b|} \left( \prod_{s \in b} A_s \right). $$

Furthermore, when the point process $\Phi$ is stationary, we may consider for any $m \ge 2$ the reduced $m$-th order factorial cumulant moment measure $\gamma_{\mathrm{red}}^{(m)}$. This is a locally finite signed measure on $(\mathbb{R}^d)^{(m-1)}$ that satisfies

$$ \gamma^{(m)}(A_1 \times \dots \times A_m) = \int_{A_m} \gamma_{\mathrm{red}}^{(m)} \left( \prod_{i=1}^{m-1} (A_i - x) \right) dx, \quad (5) $$

where for $i = 1, \dots, m-1$, $(A_i - x)$ is the translation of the set $A_i$ by $x$.
# 4 Minimax lower bounds

In this section, Theorem 4.1 and Theorem 4.2 provide minimax lower bounds for the problems of estimating the structure factor and the pair correlation function of a stationary spatial point process, based on its observation over a compact window. They state that we recover the classical non parametric rate: the expected number of observed points at power $-\beta/(2\beta+d)$, when considering the problem of estimating a function of Hölder-regularity $\beta$, as defined below.

**Definition 1.** Let $\beta > 0$ and $L > 0$. We denote by $l$ the greatest integer strictly less than $\beta$. A function $f: \mathbb{R}^d \to \mathbb{R}$ belongs to the Hölder-class $\Theta(\beta, L)$, if $f$ is $l$-times differentiable and for all multi-indexes $i = (i_1, \dots, i_d) \in \mathbb{N}^d$ such that $i_1 + \dots + i_d = l$, we have

$$ \forall (x, y) \in \mathbb{R}^{2d}, |\partial_i^l f(x) - \partial_i^l f(y)| \le L|x - y|^{\beta-l}. $$

To state the minimax lower bounds, we introduce a notion of a family of regular compact sets. For example, $W = (\{Rw | w \in W_0\})_{R>0}$, where $W_0$ is a non empty compact set, is regular.

**Definition 2.** A family $W$ of compact sets of $\mathbb{R}^d$ is said regular if there exists $\eta > 0$ such that for all $W \in W$ there exists $R_W > 0$ such that $W \subset [-R_W, R_W]^d$ and $|[-R_W, R_W]^d| \le \eta|W|$.

The following theorem, whose proof is postponed to Section 7, establishes the minimax rate of convergence for estimating the structure factor of Hölder regularity $\beta$ of a point process $\Phi$, based on the observation of $\Phi \cap W$, where $W$ is a compact observation window. This result is stated in an expanding window regime, which is a standard asymptotic setting for structure factor estimation [64, 36, 83]. Finally, as we have assumed a unit intensity, the expected number of observed points is $E[|\Phi \cap W|] = |W|$, and the minimax rate is in fact expressed in terms of this expected count.

**Theorem 4.1.** Let $h: [0, \infty) \to [0, \infty)$ be an non-identically null increasing function with $h(0) = 0$. Consider $W$, a regular family of compact sets of $\mathbb{R}^d$. Let $(\beta, M, L) \in (0, \infty)^3$. For $S \in \Theta(\beta, L)$, we denote by $\mathcal{L}_{(S,M)}$ the set all laws of stationary point processes of unit intensity satisfying Assumption 3.1, having structure factor $S$ and such that $\max_{m \in \{2,3,4\}} |\gamma_{red}^{(m)}|(\mathbb{R}^{d(m-1)}) \le M$. Then,

$$ \inf_{W \in W} \inf_{\hat{S}} \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \sup_{k_0 \in \mathbb{R}^d} \sup_{\mathbb{E}} \left[ h \left( |W|^{\frac{\beta}{2\beta+d}} |\hat{S}(k_0; \Phi \cap W) - S(k_0)| \right) \right] > 0, \quad (6) $$

where $\sup_{S \in \Theta(\beta, L)}$ denotes the supremun over structure factors of regularity $\Theta(\beta, L)$ and $\inf_{\hat{S}}$ denotes the infimum over all the estimators $\hat{S}: \mathbb{R}^d \times \text{Conf}(\mathbb{R}^d) \to \mathbb{R}$.

Theorem 4.1 implies that the estimator of $S(0)$ constructed in [39] is minimax when $\beta \ge d^2/2$ (see Theorem A.3 in [39]). While this manuscript was being written, the preprint [24] introduced a kernel estimator of the structure factor that is minimax in dimensions $d \le 4$ for $\beta = 2$. In Section 5, we construct, via multitapering, a minimax estimator of $S$ for $\beta \in (0, 2]$.

We comment the obtained minimax rate of convergence. It is faster for large $\beta$, as expected, since larger value of $\beta$ are related to faster decay of correlations. Indeed, for a point process $\Phi$ having pair correlation function $g$, then if $x^\beta(g(x) - 1) \in L^1(\mathbb{R}^d)$, i.e. if $\Phi$ has second order structure that decorrelates quickly enough, then classical Fourier arguments imply that $S \in \Theta(\beta, L)$ for some $L \ge 0$ (see, e.g., [29]). To make it more concrete, for DDPs [74, 51] with stationary kernel $K \in L^2(\mathbb{R}^d \times \mathbb{R}^d)$ satisfying $|K(x, y)|^2 = k_0^2(x - y)$ for all $(x, y) \in \mathbb{R}^{2d}$, this condition can be rephrased $x^\beta k_0(x)^2 \in L^1(\mathbb{R}^d)$. For a Cox process [56] with directional intensity $(Z(x))_{x \in \mathbb{R}^d}$, where $Z(x)$ is a stationary random field, it reads, $x^\beta \text{Cov}(Z(x), Z(0)) \in L^1(\mathbb{R}^d)$.
As outlined in the introduction, the proof of Theorem 4.1 differs from classical minimax arguments for spectral estimation in time series, which rely on Gaussian processes [70, 9, 27, 1, 45]. Such constructions are not directly applicable to point processes, and existing minimax results for first-order estimation in non-stationary Poisson processes [67, 49, 68] yield trivial bounds for second-order spectra estimation for stationary spatial point processes.

To overcome this, we study a Cox process $\Phi_{\text{Cox}}$ with underlying intensity field $(1+\sin(N(x)))_{x \in \mathbb{R}^d}$, where $(N(x))_{x \in \mathbb{R}^d}$ is a centered stationary Gaussian field with covariance function $c$. The key idea is that as $\text{Var}[N(0)]$ tends to zero, the distribution of $\Phi_{\text{Cox}}$ converges to that of a stationary Poisson point process with intensity 1 (see Lemma 7.1). At the same time, by increasing the range of correlations of the Gaussian field $N$, we ensure that the second-order properties of $\Phi_{\text{Cox}}$ remain sufficiently distinct from those of a Poisson process. To this end, we set $c = \sigma^2 c_0(\cdot/\rho)$ as the covariance function of the Gaussian field, where $c_0$ is compactly supported, and let $\rho \to \infty$ and $\sigma \to 0$ simultaneously, with a proper scaling with respect to the volume $|W|$ of the observation window. This construction enables us to apply the $\chi^2$-divergence version of Theorem 2.2 from [77].

An important technical point in the proof, specific to point processes, concerns the control of higher-order cumulants. Since we let $\rho \to \infty$, meaning that the range of correlations in $\Phi_{\text{Cox}}$ diverges, it is not immediate that the total variation of the third and fourth order reduced factorial cumulant measures remains bounded. However, this control is crucial, as $L^2(\mathbb{P})$-risk upper bounds typically rely on the finiteness of these cumulants or related quantities. For Gaussian processes, this issue does not arise since joint cumulants of order higher than two vanish. In our setting, we prove that for $\Phi_{\text{Cox}}$, we have $\max_{m \in \{2,3,4\}} |\gamma_{\text{red}}^{(m)}|(\mathbb{R}^{d(m-1)}) \le M$, for some constant $M$, by leveraging the fact that when $\sigma \to 0$, the underlying intensity becomes nearly Gaussian: $1+\sin(N(x)) \simeq 1+N(x)$ (see Lemma 7.2).

To conclude this section, we note that the previous proof can be adapted with minor modifications to establish a minimax lower bound for the estimation of the pair correlation function of a point process. This result is stated in the next theorem, which is proved in Section 7. We also mention that the estimator introduced in [82] achieves the minimax convergence rate corresponding to the following theorem when $\beta \in \mathbb{N} \setminus \{0\}$ (see Lemma 3.2 in [82]).

**Theorem 4.2.** Let $(\beta, M, L) \in (0, \infty)^3$ and $h: [0, \infty) \to [0, \infty)$ be non-identically null, increasing and with $h(0) = 0$. Let $W$ be a regular family of compact sets of $\mathbb{R}^d$. For $g \in \Theta(\beta, M)$, we denote by $\mathcal{L}_{(g,M)}$ the set all laws of stationary point processes of unit intensity satisfying Assumption 3.1, having pair correlation function $g$ and such that $\max_{m \in \{2,3,4\}} |\gamma_{\text{red}}^{(m)}|(\mathbb{R}^{d(m-1)}) \le M$. Then,

$$ \inf_{W \in \mathcal{W}} \inf_{\hat{g}} \sup_{g \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(g,M)}} \sup_{x_0 \in \mathbb{R}^d} \mathbb{E} \left[ h \left( |W|^{\frac{\beta}{2\beta+d}} |\hat{g}(x_0; \Phi \cap W) - g(x_0)| \right) \right] > 0, \quad (7) $$

where $\sup_{g \in \Theta(\beta, L)}$ denotes the supremun over pair correlation functions of regularity $\Theta(\beta, L)$ and $\inf_{\hat{g}}$ denotes the infimum over all the estimators $\hat{g}: \mathbb{R}^d \times \text{Conf}(\mathbb{R}^d) \to \mathbb{R}$.

# 5 Minimax estimator of the structure factor

In this section, we construct an estimator based on multitapering that achieves the minimax rate of convergence stated in Theorem 4.1. Specifically, Section 5.1 introduces a multitaper estimator. Then, Section 5.2 analyses its squared risk. Finally, Section 5.3 establishes that the $L^2(\mathbb{P})$ minimax bound can be strengthened to a high-probability bound. As before, we assume for simplicity that the intensity is $\lambda = 1$. Appendix E extends the results to practical estimator, where the unknown intensity is replaced by the standard estimator $\tilde{\lambda} = |\Phi \cap W|/|W|$.
## 5.1 Multitaper estimator: introduction

The multitaper estimator of the structure factor of a spatial point process is defined as follows.

**Definition 3.** Let $\Phi$ be a stationary point process of $\mathbb{R}^d$, $W$ be a compact set of $\mathbb{R}^d$ and $I$ be a subset of $\mathbb{N}^d$. Let $k_0 \in \mathbb{R}^d$ be a fixed frequency. We consider a family of orthogonal tapers functions $(f_i)_{i \in I} \in L^2(\mathbb{R}^d)^{|I|}$. The multitaper estimator of the structure factor $S$ of $\Phi$ associated to this family is defined as

$$ \hat{S}(k_0) := \frac{1}{|I|} \sum_{i \in I} \hat{S}_i(k_0), \qquad (8) $$

where the single-taper estimator of $S$ is $\hat{S}_i(k_0) := |C_i(k_0)|^2$, with

$$ C_i(k_0) := T_i(k_0) - \mathbb{E}[T_i(k_0)], \quad T_i(k_0) := \sum_{x \in \Phi \cap W} e^{-ik_0 \cdot x} f_i(x), \quad \mathbb{E}[T_i(k_0)] = \int_W f_i(x)dx. \qquad (9) $$

**Remark 1.** Unlike previous works on taper-based estimators for spatial point processes [64, 83, 31], which consider tapers supported within the observation window $W$ and define $T_i(k_0)$ via sums over all points $x \in \Phi$, we consider here a more general setting where the tapers may have infinite support. However, in the following, we will require strong localization properties within the window to control the truncation error.

The standard framework for analyzing the statistical properties of multitaper estimators in spatial statistics is the expanding-window regime, where the observation window $W$ grows to fill $\mathbb{R}^d$. In this setting, the tapers are chosen so that they become increasingly concentrated at the origin in the Fourier domain, i.e., $\mathcal{F}[f_i] \to \delta_0$ as $W \to \mathbb{R}^d$ for every $i \in I$. When the number of tapers $|I|$ is kept fixed while $W$ expands, it is well established that multitaper estimators are asymptotically unbiased. Furthermore, under suitable mixing conditions, the linear statistics $(T_i(k_0))_{i \in I}$ (defined in (9)) become asymptotically Gaussian and pairwise uncorrelated. As a result, one obtains the following convergence in distribution, for $k_0 \neq 0$ (see [15, 57, 64, 31] for spatial point processes and [81, 44], Section 8 of [62] for similar time series results):

$$ \hat{S}(k_0) \xrightarrow{d} S(k_0) \frac{\chi^2(2|I|)}{2|I|}, \quad \text{as } W \to \mathbb{R}^d. \qquad (10) $$

This limit shows, however, that if $|I|$ remains fixed as $W$ grows, then $\hat{S}(k_0)$ is not consistent. This key observation, emphasized in the works cited above, compels us to let the number of tapers $|I|$ increase with the window size $W$. But this raises theoretical challenges that are critical for the implementation and the choice of the parameters of the multitaper estimators. For instance, increasing $|I|$ reduces variance by averaging, but at the same time it typically increases bias, since higher-order tapers have poorer spectral concentration. This naturally leads to a set of fundamental questions:

* How should $|I|$ scale with $W$ in order to preserve asymptotic unbiasedness?

* Does the convergence in distribution in (10) still hold when $|I|$ grows with $W$?

* More broadly, how does the rate of convergence depend on properties of the taper family, such as spatial support and spectral concentration?
The answers to the first and third questions are well established in the time series literature, typically under the assumption that the spectral density is twice differentiable (see, e.g., [75, 69, 62]). Recent works on spatial point processes [64, 31] have addressed analogous questions, emphasizing the asymptotic interplay between spatial and spectral localization of the tapers. Within our minimax multitaper construction for the structure factor, we obtain a non-asymptotic $L^2(\mathbb{P})$ risk bound (Theorem 5.1) that encapsulates this trade-off and connects back to these classical results. By contrast, the second question appears to be less thoroughly investigated in the existing time series literature (see, e.g., [44]). In our spatial point process setting, we address this question by establishing a deviation inequality under a (Brillinger) mixing condition, which is well suited to point processes.

## 5.2 Multitaper estimators: mean integrated squared risk

This section establishes an upper bound on the pointwise $L^2(\mathbb{P})$-risk of the multitaper estimator of the structure factor, as stated in Theorem 5.1. The result holds for any orthonormal family of tapers in $L^2(\mathbb{R}^d)$. It is further specialized in Corollary 5.1 to the case of Hermite tapers, yielding the minimax rate of convergence. Theorem 5.1 and Corollary 5.1 are proved in Appendix A, where the variance and the bias of the multitaper estimators are analyzed separately in Lemma A.1 and Lemma A.2.

**Theorem 5.1.** Let $\Phi$ be a stationary point process satisfying Assumption 3.1 with $\lambda = 1$. Let $W$ be a subset of $\mathbb{R}^d$, $I$ be a discrete subset of $\mathbb{N}^d$. We consider a family of orthonormal functions $(f_i)_{i \in I}$ of $L^2(\mathbb{R}^d)$. Assume that $S \in \Theta(\beta, L)$ with $L \in (0, \infty)$ and $\beta \in (0, 2]$ (see Definition 1). Then,

$$ \sup_{k_0 \in \mathbb{R}^d} \mathbb{E} \left[ \left| \hat{S}(k_0) - S(k_0) \right|^2 \right]^{\frac{1}{2}} \le \frac{\sqrt{2} \|S\|_{\infty}}{|I|^{\frac{1}{2}}} + \mathcal{F}\text{-loc} + B_4 + 6 \|S\|_{\infty} W\text{-loc}, \quad (11) $$

where

$$ \begin{aligned} \mathcal{F}\text{-loc} &= \sqrt{dL} \frac{1}{|I|} \sum_{i \in I} \|f_i\|_{H^{\beta/2}}^2, \quad W\text{-loc} := \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_{1\mathbb{R}^d \setminus W} \|_2^2 \right)^{\frac{1}{2}}, \\ B_4 &:= \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_4^2 \right) \left( 1 + 7|\gamma_{red}^{(2)}(\mathbb{R}^d)| + 6|\gamma_{red}^{(3)}(\mathbb{R}^{2d})| + |\gamma_{red}^{(4)}(\mathbb{R}^{3d})| \right)^{\frac{1}{2}}. \end{aligned} $$

We now discuss the quantities involved in the upper bound (11). The term $\|S\|_{\infty}/|I|^{1/2}$ follows from the orthonormality of the taper family $(f_i)_{i \in I}$, reflecting the asymptotic pairwise decorrelation of the linear statistics $(T_i(k_0))_{i \in I}$, defined in (9), as $W \to \mathbb{R}^d$ and $\mathcal{F}[f_i] \to \delta_0$ for all $i \in I$, at a fixed number of tapers $|I|$.

The remaining terms in (11), namely $\mathcal{F}\text{-loc}$, $W\text{-loc}$, and $B_4$, capture non-asymptotic effects. These quantities reflect, respectively, the lack of spectral concentration of the tapers around the origin, their imperfect spatial localization within the observation window $W$, and higher-order cumulants' contributions. We analyze them below.

1. $W\text{-loc}$: border effects. The term $W\text{-loc}$ arises because only the points of $\Phi$ within $W$ are observed. Consequently, when the taper functions $(f_i)_{i \in I}$ have too much mass outside of $W$, border effects occur. Accordingly, to ensure negligible border effects, one should consider tapers with good spatial localization inside $W$.
2. B$_4$: Fourth order points correlations. The term B$_4$ keeps track of the higher order correlations present in the underlying point process. This term is expected as the previous theorem deals with a squared risk of a second order statistics of stationary point processes. Note that this term is finite for Brillinger-mixing point processes (see Definition 5.1), including most DPPs [38, 11] and Cox processes [40, 84, 63, 43]. In addition, we observe a multiplicative term involving the average of the L$^4$(R$^d$)-norm of the tapers. This term is related to the spatial localization of the tapers functions: for well-localized taper functions, this average decreases to zero as the number of tapers |I| increases, providing a second argument in favor of well-localized tapers.

3. F-loc: Frequency-localization bias. The term F-loc represents the bias arising from the localization of the taper functions in the Fourier domain. Since this term typically diverges as |I| → ∞ at fixed W, it highlights the non-asymptotic trade-off between using many tapers, which reduces variance but increases bias, and using few tapers, which decreases bias but results in larger variance. In addition, its expression naturally motivates the use of an orthogonal family of tapers with small Sobolev norms, and hence with well-localized Fourier transforms. In this way, we recover, within a non-asymptotic framework, a classical requirement in taper theory [69, 62, 44, 64].

The previous discussion motivates the consideration of taper functions with good localization property in both the spatial and the Fourier domain. Accordingly, a balance must be found between good spatial localization, to reduce the bias W-loc caused by border effects and the variance due to higher-order point correlations B$_4$, and good localization in the Fourier domain, to reduce the bias F-loc. Common choices of families satisfying both requirements are the Slepian functions [73, 50, 72], sinusoidal tapers [69, 81] and the Hermite tapers [7, 28, 60, 22]. In this work, we focus on the latter family.

**Definition 4.** For $i = (i_1, \dots, i_d) \in \mathbb{N}^d$, the Hermite taper $\psi_i$ is defined as

$$ \forall x = (x_1, \dots, x_d) \in \mathbb{R}^d, \quad \psi_i(x) = e^{-\frac{1}{2}|x|^2} \prod_{l=1}^{d} H_{i_l}(x_l), \qquad (12) $$

where for $n \in \mathbb{N}$, $y \in \mathbb{R}$, $H_n(y) = \frac{(-1)^n}{(2^n n! \sqrt{\pi})^{\frac{1}{2}}} e^{y^2} \frac{d^n}{dy^n} e^{-y^2}$ are the Hermite polynomials.

The next corollary simplifies the terms appearing in Theorem 5.1 when considering scaled Hermite tapers. It states that, under appropriate choices of the parameters, the corresponding multitaper estimator is optimal regarding the minimax rate of Theorem 4.1, for structure factor of Hölder-regularity $\beta \le 2$ and for square observation windows.

**Corollary 5.1.** Let $\Phi$ be a stationary point process satisfying Assumption 3.1 with $\lambda = 1$. Assume that $S \in \Theta(\beta, L)$ with $L \in (0, \infty)$ and $\beta \in (0, 2]$ (see Definition 1). Let $R > 0$, $W = [-R, R]^d$ and $I = \{i \in \mathbb{N}^d | |i|_{\infty} < i_{\max}\}$ with $i_{\max} \ge 1$. We consider, for all $i \in I$, $f_i = r^{-\frac{d}{2}} \psi_i(\cdot/r)$, where $\psi_i$ is the i-th Hermite function introduced in Definition 4 and

$$ r = R/\sqrt{2i_{\max} + i_{\max}^{1/3+\theta}} \qquad (13) $$

where $\theta \in (0, 2/3)$. Suppose that

$$ |I| \left( |W|^{2\beta/(2\beta+d)} \left( \frac{\|S\|_{\infty}}{L} \right)^{\frac{d}{2\beta+d}} \right)^{-1} \in (c_2, c_3), \qquad (14) $$
for some constants $0 < c_2 \le c_3 < \infty$. Then, we have

$$ \sup_{k_0 \in \mathbb{R}^d} \mathbb{E} \left[ \left| \hat{S}(k_0) - S(k_0) \right|^2 \right]^{\frac{1}{2}} \le K \left( \frac{(L^d \|S\|_{\infty}^{\beta})^{\frac{1}{2\beta+d}}}{|W|^{\frac{\beta}{2\beta+d}}} + \sqrt{\frac{|\gamma_{red}^4| \log(|W|)^d}{|W|}} + e^{-c|W|^\alpha} \right), $$

where $|\gamma_{red}^4| := 1 + 7|\gamma_{red}^{(2)}|(\mathbb{R}^d) + 6|\gamma_{red}^{(3)}|(\mathbb{R}^{2d}) + |\gamma_{red}^{(4)}|(\mathbb{R}^{3d})$, $\alpha = 3\beta\theta/(d(2\beta+d))$ and $0 < c, K < \infty$ are constants that do not depend on $W$, $r$ and $I$.

In the previous result, condition (13) ensures exponentially negligible border effects W-loc by localizing all the Hermite tapers $(\psi_i(\cdot/r))_{i \in I}$ in the observation window $W$. Then, condition (14), concerning the number of tapers, balances the frequency bias term F-loc of order $L(|I|/|W|)^{\beta/d}$, and the asymptotic variance term $\|S\|_{\infty}|I|^{-1/2}$ (see equation (42)). It leads to the dominant term of order $(L^d \|S\|_{\infty}^{\beta})^{1/(2\beta+d)}/|W|^{\beta/(2\beta+d)}$. The term accounting for higher-order correlations, namely $|\gamma_{red}^4|^{1/2} \log(|W|)^{d/2}/|W|^{1/2}$, vanishes at the parametric rate (up to a logarithmic factor).

## 5.3 Multitaper estimators: non asymptotic fluctuations

In this section, we establish in Theorem 5.2 a high-probability counterpart of Theorem 5.1, which provides a justification of (10) when the number of tapers diverges as the observation window grows. Then, in Corollary 5.2, we specify this bound for Hermite tapers.

The proof of Theorem 5.2, deferred to Appendix B, relies on bounding the cumulants of the multitaper estimator $\hat{S}(k_0)$ and on showing that the orthogonality of the tapers, which at first sight affects only second-order moments, in fact propagates to higher-order cumulants. Then, the deviation inequality follows from the method of cumulants [10, 71, 25]. Accordingly, the results in this section rely on the following strong Brillinger mixing condition.

**Assumption 5.1.** A simple stationary point process $\Phi$ is said is $(c, \gamma, q)$-Brillinger-mixing if there exist $c \ge 1, \gamma \ge 0$ and $q \ge 1$, such that for all $m \ge 2$

$$ |\gamma_{red}^{(m)}|(\mathbb{R}^{d(m-1)}) \le c((m-1)!)^{1+\gamma}q^{m-1}. $$

Assumption 5.1 is standard in spatial statistic and allows for a quantitative analysis of the fluctuations of linear statistics of point processes; see, e.g., [15, 16, 40, 37, 38, 46]. When $\gamma = 0$, it can be seen as a sub-exponential condition (see, e.g., [80], Chapter 2). Note that, for a stationary $\alpha$-determinantal point process with stationary kernel $K \in L^1(\mathbb{R}^d)$, previous assumption holds with $\gamma = 0, c = 1$ and $q = \|K\|_1$ (refer to Theorem 2 of [38]). Assumption 5.1 is well-suited to our setting, since the multitaper estimator is built from linear statistics.

The following theorem provides a high-probability bound on the fluctuations of the multitaper estimator. It highlights three distinct concentration regimes, which are detailed in the discussion following the statement.

**Theorem 5.2.** Let $I$ be a discrete subset of $\mathbb{N}^d$, $W$ be a compact set of $\mathbb{R}^d$ and $(f_i)_{i \in I}$ be a family of orthonormal functions of $L^2(\mathbb{R}^d)$ that are in $L^1(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d)$. Assume that there exist $c_f \in (1, \infty)$, $\eta_1, \eta_2 \ge 0$ and $\rho > 0$ such that

$$ \forall s \ge 4, \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_s^2 \right)^{\frac{s}{2}} \le \frac{c_f \log(|I|^{\frac{1}{d}})}{\rho^{d(s/2-1)} |I|^{s\eta_1 + \eta_2}}, \quad (15) $$
$$ \forall s \ge 3, \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_s^2 \right)^{\frac{2}{s}} \le \frac{c_f}{\rho^{d(s/2-1)} |I|^{s\eta_1}}. \quad (16) $$

Let $\Phi$ be a stationary point process satisfying Assumption 3.1 with $\lambda = 1$. Moreover, suppose that $\Phi$ is $(c, \gamma, q)$-Brillinger-mixing (see Assumption 5.1) and that $S \in \Theta(\beta, L)$ with $L \in (0, \infty)$ and $\beta \in (0, 2]$. If $\beta > 1$, assume for all $(i, i') \in I^2$, $\mathcal{F}[f_i](k)\mathcal{F}[f_{i'}](k) = \mathcal{F}[f_i](-k)\mathcal{F}[f_{i'}](-k)$ for all $k \in \mathbb{R}^d$. Suppose that

$$ c\rho^d \ge (1+q) \frac{|I|^{\eta_2}}{\log(|I|^{\frac{1}{d}})} \vee 1. \quad (17) $$

Let $k_0 \in \mathbb{R}^d$. We denote

$$ \Sigma := \frac{S(k_0)}{|I|} + \sqrt{dL} \left( \frac{1}{|I|^2} \sum_{i \in I} \|f_i\|_{H^\beta}^2 \right)^{\frac{1}{2}} + 4\|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_{1\mathbb{R}^d \setminus W}^2 \right)^{\frac{1}{2}} \\ + \frac{\|S\|_{\infty}}{|I|^{\frac{1}{2}}} \wedge 2\|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|\mathcal{F}[f_i]\|_{1\setminus|\cdot| \ge |k_0|}^2 \right)^{\frac{1}{2}}, \quad (18) $$

$$ \Lambda := \left[ c_f \vee \frac{6c\rho^{\frac{3d}{2}} |I|^{3\eta_1}}{1+q} \left( \sup_{i \in I} \|f_i\|_{1\mathbb{R}^d \setminus W} \|3\|_{i \in I} \sup_{i \in I} \|f_i\|_{3}^{2} + \sup_{i \in I} \|\mathcal{F}[f_i]\|_{1\setminus|\cdot| \ge |k_0|/3} \|1\|_{i \in I} \right) \right], \quad (19) $$

$$ B_{\infty} := \frac{\Lambda^2 c^{1/2} (1+q)^{3/2} \log(|I|^{1/2})^{1/2}}{\rho^{d/2} |I|^{2\eta_1+\eta_2/2}} \wedge \frac{[c_f]^2 c^{2/3} (1+q)^{4/3}}{\rho^{d/3} |I|^{2\eta_1}}. $$

Then, for all $x \ge 0$, we have

$$ \mathbb{P}[|\hat{S}(k_0) - \mathbb{E}[\hat{S}(k_0)]| \ge x] \le C \max(a_1, a_2, a_3), \quad (20) $$

where $C < \infty$ is a numerical constant and

$$ a_1 := \exp \left( - \frac{1}{12} \frac{x^2 |I|}{S(k_0)^2 + xS(k_0)} \right), \quad a_2 := \exp \left( - \frac{x^2}{13^\gamma 39 \Sigma^2 + 1.5 x^2 (x\Sigma)^{-1/(1+\gamma)}} \right), $$

$$ a_3 := \exp \left( - \frac{x^2}{13^\gamma (71 B_\infty)^2 + 1.5 x^2 (2x B_\infty)^{-1/(2+\gamma)}} \right). $$

We begin by discussing the assumptions of the previous theorem. Assumptions (15) and (16) concern the localization properties of the tapers and are required to control high-order cumulants. The parameter $\rho$ in these conditions should be thought of as a scaling factor related to the localization of the tapers within the observation window $W$, analogous to the dilation factor $r$ in Corollary 5.1. Assumption (17) is of a technical nature. Additionally, the symmetry condition on the Fourier transform of the tapers is used to cancel a first-order term in the Taylor expansion of the structure factor near $k_0$.

The terms $a_1$, $a_2$, and $a_3$ in (20) correspond to different concentration types: $\chi^2$, $1+\gamma$ stretched exponential, and $2+\gamma$ stretched exponential, respectively, as detailed below.

1. $\chi^2$ concentration (fixed-$|I|$ regime). The term $a_1$ reflects the concentration of the multitaper estimator in the asymptotic regime where $|I|$ is fixed and $\mathcal{F}[f_i] \to \delta_0$ for all $i \in I$ as $W \to \mathbb{R}^d$ (see the discussion following Remark 1). Indeed, in this setting, $\hat{S}(k_0)$ is approximately
distributed as $S(k_0)\chi^2(2|I|)/(2|I|)$ (see (10)). Remarkably, Corollary 5.2 above shows that even when the number of tapers $|I|$ diverges with $W$, provided it scales appropriately, this term remains the leading contribution to the fluctuations of $\hat{S}(k_0)$, when $k_0$ is not too close to zero.

2. $1 + \gamma$ stretched exponential concentration (second order decorrelation error). The term $a_2$ accounts for the fact that the linear statistics $T_i(k_0)$ (see Definition 3) are only asymptotically uncorrelated. For $i_1, i_2 \in I$,

$$ \mathrm{Cov}[T_{i_1}(k_0), T_{i_2}(k_0)] = S(k_0)\mathbf{1}_{i_1=i_2} + \mathrm{Err}(i_1, i_2), $$

where the error term $\mathrm{Err}(i_1, i_2)$ captures residual non-asymptotic correlations. This non-asymptotic correlation leads to an $1 + \gamma$ stretched exponential deviation component. Whether $a_2$ dominates depends on the decay of

$$ \frac{\|S\|_{\infty}}{|I|^{1/2}} \wedge 2\|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|\mathcal{F}[f_i]\mathbf{1}_{|\cdot| \ge k_0}\|_2^2 \right)^{1/2}, $$

which typically depends on how small $k_0$ is.

3. $2 + \gamma$ stretched exponential concentration (high-order correlations). The term $a_3$ captures the contribution of higher-order cumulants of the statistics $(T_i(k_0))_{i \in I}$ (see Definition 3). For the discussion, assume that $\gamma = 0$. Under Brillinger-mixing, the total variation of the $m$-th order reduced factorial cumulant measure are bounded by $(m-1)!$ (up to constants). But the sequence $((m-1)!)_{m \ge 1}$ is the cumulants of a standard exponential random variable. Accordingly, the higher-order cumulants of the statistics $\hat{S}_i(k_0) = |T_i(k_0) - \mathbb{E}[T_i(k_0)]|^2$, for $i \in I$ behave, in the non-asymptotic regime, similarly to the higher-order cumulants of squares of centered exponential random variable variables. The contribution of $a_3$ is determined by $B_{\infty}$, which splits into two cases. To fix the idea, in the case of the Hermite tapers (see the next corollary), if $\log(|I|^{1/d})^{1/2}/(\rho^{d/2}|I|^{2\eta_1+\eta_2/2})$ dominates, then $B_{\infty}$ matches the fourth-order cumulant bound $B_4$ from Theorem 5.1. If instead $1/(\rho^{d/3}|I|^{2\eta_1})$ dominates (which could occur for small $k_0$ since in that case $\Lambda$ could be large), the higher-order joint cumulants contribute more significantly than the term $B_4$ appearing in the $L^2(\mathbb{P})$-risk.

The following corollary establishes that the multitaper estimator based on Hermite tapers concentrates at the minimax rate with $\chi^2$-type concentration, anticipated by (10), provided the target frequency $k_0$ is not too small. When $k_0$ is small, the estimator still achieves the minimax rate if $\beta \le d$, but with stretched-exponential concentrations. We detail these points in the discussion following the theorem.

**Corollary 5.2.** Let $\Phi$ be a stationary point process satisfying Assumption 3.1 with $\lambda = 1$. Moreover, suppose that $\Phi$ is $(c, \gamma, q)$-Brillinger-mixing (see Assumption 5.1) and that $S \in \Theta(\beta, L)$ with $L \in (0, \infty)$ and $\beta \in (0, 2]$ (see Definition 1). Let $R > 0$ and $W = [-R, R]^d$. Consider $I = \{i \in \mathbb{N}^d | |i|_{\infty} < i_{\max}\}$ with $i_{\max} \ge 1$. Let $\epsilon \in \{-1, 1\}$. If $\beta \ge 1$, we add the constraint $\psi_i(-\cdot) = \epsilon\psi_i$ in the indexes of $I$, where $\psi_i$ is the $i$-th Hermite function introduced in Definition 4. We consider, for all $i \in I$, $f_i = r^{-d/2}\psi_i(\cdot/r)$ where

$$ r = R/\sqrt{2i_{\max} + i_{\max}^{1/3+\theta}}, \qquad (21) $$

where $\theta \in (0, 2/3)$. Suppose that $|I| \ge 2$ and

$$ |I||W|^{-2\beta/(2\beta+d)} \in (c_2, c_3), \qquad (22) $$
for some constants $0 < c_2 \le c_3 < \infty$. Then, there exists $v_0, c_0, c_4 > 0$ and $K < \infty$ such that for $|W| \ge v_0$ and $x \ge 0$, we have with probability at least $1-x$:

$$ \sup_{|k_0| \ge |k_W^*|} \mathbb{P} \left[ \left| \hat{S}(k_0) - S(k_0) \right| \ge x + \frac{K}{|W|^{\frac{\beta}{2\beta+d}}} \right] \le K \left( \exp \left( -\frac{c_2^2}{12} \frac{(x|W|^{\frac{\beta}{2\beta+d}})^2}{S(k_0)^2 + xS(k_0)} \right) + \varepsilon_1 \right), \quad (23) $$

where $|k_W^*| := c_0|W|^{-\frac{1}{2\beta+d}}$ and

$$ \varepsilon_1 := \exp \left( -c_4(x|W|^{\frac{2\beta}{2\beta+d}})^2 \wedge (x|W|^{\frac{2\beta}{2\beta+d}})^{\frac{1}{1+\gamma}} \right) + \exp \left( -c_4 \left( \frac{x|W|^{\frac{1}{2}}}{\log(|W|)^{\frac{\beta}{2}}} \right)^2 \wedge \left( \frac{x|W|^{\frac{1}{2}}}{\log(|W|)^{\frac{\beta}{2}}} \right)^{\frac{1}{2+\gamma}} \right). $$

Moreover,

$$ \sup_{|k_0| \le |k_W^*|} \mathbb{P} \left[ \left| \hat{S}(k_0) - S(k_0) \right| \ge x + \frac{K}{|W|^{\frac{\beta}{2\beta+d}}} \right] \le K(A_1 + A_2 + A_3), \quad (24) $$

where

$$ A_1 := \exp \left( -\frac{c_2^2}{12} \frac{(x|W|^{\frac{\beta}{2\beta+d}})^2}{S(k_0)^2 + xS(k_0)} \right), \quad A_2 := \exp \left( -c_4(x|W|^{\frac{\beta}{2\beta+d}})^2 \wedge (x|W|^{\frac{\beta}{2\beta+d}})^{\frac{1}{1+\gamma}} \right), $$

$$ A_3 := \exp \left( -c_4(x|W|^{\frac{1}{3}})^2 \wedge (x|W|^{\frac{1}{3}})^{\frac{1}{2+\gamma}} \right). $$

We now discuss the assumptions of the theorem. Conditions (21) and (22) are identical to those of Corollary 5.1. Condition (21) localizes the tapers inside the observation window while Condition (22) balances bias and variance. The threshold $|k_W^*|$, corresponding to the minimal allowed frequency, separates the regime of large frequencies $|k_0| \ge |k_W^*|$, where $T_i(k_0)$ and $T_i(-k_0)$ (see Definition 3) asymptotically decorrelate, from the small frequency regime $|k_0| \le |k_W^*|$, where they do not.

For large frequencies $|k_0| \ge |k_W^*|$, the remainder term $\varepsilon_1$, involving the rates $|W|^{2\beta/(2\beta+d)}$ and $|W|^{1/2}/\log(|W|)^{d/2}$, becomes negligible as $|W| \to \infty$. For this situation, even if $|I|$ diverges with $W$, the bound matches the minimax rate, with the anticipated $\chi^2(2|I|)/(2|I|)$ concentration occurring at fixed $|I|$ (see the heuristic equation (10)). In contrast, for small frequencies $|k_0| \le |k_W^*|$, the terms $A_1$ and $A_2$ contribute to the bound. For small deviations $x$, $A_2$, which captures the non-asymptotic correlation between tapers, dominates $A_1$. The third term, $A_3$, related to higher-order correlations between single-taper estimators, is negligible compared to $A_1$ and $A_2$ if $\beta < d$, since in this case $1/3 > \beta/(2\beta+d)$. Accordingly, in that case the multitaper estimator still concentrates at the minimax rate, but with $1 + \gamma$ stretched exponential concentration. However, when $\beta = d$, $A_3$ dominates $A_1$ and $A_2$ for small $x$. But, the bound remains minimax with a $2 + \gamma$ stretched exponential concentration. Since $\beta$ is assumed to be smaller than two, the Hermite multitaper estimator always concentrates at the minimax rate in dimensions $d \ge 2$, under the conditions of the previous corollary.

# 6 Data-driven procedure for the choice of the number of tapers

To select the number of tapers in practice, we design a local cross-validation procedure. In Section 6.1, we describe the approach and provide theoretical justifications. Then, in Section 6.2, we evaluate its performances through numerical simulations. The implementation used for these experiments is available at https://github.com/gabrielmastrilli/multitaper_spp.
## 6.1 Description and theoretical guarantee

To select the number of tapers for an estimation of S at frequency $k_0$, we aim to choose $|I|$ to minimize the MSE (mean squared error) between $\hat{S}_I(k_0)$ and a pilot estimator $\hat{S}_p(k_0)$. The procedure relies on the following principles.

1. Rather than considering a single frequency $k_0$, we minimize the MSE between $\hat{S}_I(k_j)$ and $\hat{S}_p(\tilde{k}_j)$ over a collection of frequency pairs $(k_j, \tilde{k}_j)_{j \in [N]}$, for $N \ge 1$. These frequencies are chosen such that, for all $j \in [N]$, $S(k_j) = S(\tilde{k}_j)$ and, on the other hand, such that $S(k_j)$ equals $S(k_0)$ or is at least close to it. The pairs are moreover chosen sufficiently separated to ensure that:

* $\hat{S}_I(k_j)$ and $\hat{S}_p(\tilde{k}_j)$ are as decorrelated as possible, and

* for $j \ne j'$, $(\hat{S}_I(k_j), \hat{S}_p(\tilde{k}_j))$ and $(\hat{S}_I(k_{j'}), \hat{S}_p(\tilde{k}_{j'}))$ are also as decorrelated as possible.

In the isotropic case, one may choose the pairs $(k_j, \tilde{k}_j)_{j \in [N]}$ as sufficiently well-spread frequencies with norm $|k_0|$.

2. To further reduce correlations, $\hat{S}_I$ and $\hat{S}_p$ are computed on two subsamples obtained by a 0.5-0.5 independent thinning of the original point pattern. This thinning is performed independently for each frequency pairs.

3. The pilot estimator $\hat{S}_p$ is chosen to have small bias, which is achieved by using a small number of tapers.

More rigorously, let $\Phi$ be a stationary point process. We denote by $N$ the number of frequency pairs $(k_j, \tilde{k}_j)_{j=1}^N$. Let $I \subset \mathbb{N}^d$ be a finite index set. For each $j \in [N]$, define $\Phi_1^j$ as an independent thinning of $\Phi$ with retention probability $1/2$. Let $\Phi_2^j := \Phi \setminus \Phi_1^j$ be its complement. The multitaper estimator $\hat{S}_I(k_j)$ is computed from $\Phi_1^j$ at frequency $k_j$ and using Hermite tapers $(r^{-d/2}\psi_i(\cdot/r))_{i \in I}$ for some dilation factor $r > 0$. Similarly, let $I_p \subset \mathbb{N}^d$ be another index set and define the pilot estimator $\hat{S}_p(\tilde{k}_j)$ computed from $\Phi_2^j$ at frequency $\tilde{k}_j$, using Hermite tapers $(r_p^{-d/2}\psi_i(\cdot/r_p))_{i \in I_p}$, for some $r_p > 0$. Let $\mathcal{I}$ be a finite set of candidate taper index sets. We define the cross-validation criterion:

$$ \hat{\mathcal{R}}(I) := \frac{1}{N} \sum_{j=1}^{N} \left( \hat{S}_{I}(k_j) - \hat{S}_{p}(\tilde{k}_{j}) \right)^2. \quad (25) $$

Finally, for the estimation of S at frequency $k_0$, we select the taper index set $I \in \mathcal{I}$ that minimizes $\hat{\mathcal{R}}(I)$.

**Remark 2.** The use of the families of thinned point processes $(\Phi_1^j)_{j \in [N]}$ and $(\Phi_2^j)_{j \in [N]}$ enhances the decorrelation between $\hat{S}_I(k_j)$ and $\hat{S}_p(\tilde{k}_j)$. For instance, in the Poisson case, for each $j \in [N]$, $\Phi_1^j$ and $\Phi_2^j$ are independent, and averaging over thinnings can be shown to be beneficial [52]. Even in the non-Poisson case, thinning has been observed to perform well for cross-validation [19]. However, thinning implies that both estimators $\hat{S}_I(k_j)$ and $\hat{S}_p(\tilde{k}_j)$ target the structure factor $S_{1/2} := 1/2 + S/2$ rather than S itself (see Lemma D.1). Nevertheless, due to the similarity between S and $S_{1/2}$, the number of tapers that minimizes the risk for estimating $S_{1/2}$ is also suitable for estimating S.

We provide with Theorem 6.1 a theoretical justification for the proposed cross-validation procedure. To state it, we introduce the notion of $k_0$-allowed frequency pairs.
**Definition 5.** Let $k_0 \in \mathbb{R}^d$ and $N \ge 1$. The family $(k_j, \tilde{k}_j)_{1 \le j \le N}$ of $\mathbb{R}^d$ is said $k_0$-allowed if $S(k_j) = S(\tilde{k}_j)$ for all $j \in [N]$ and if there exists $(c_1, c_2, c_3) \in (0, \infty)^3$ such that for all $j \in [N]$

1. $|k_j| \ge c_1|k_0|$, $|\tilde{k}_j| \ge c_1|k_0|$, $|k_j - \tilde{k}_j| \ge c_2|k_0|$ and $|k_j + \tilde{k}_j| \ge c_2|k_0|$,

2. for all $a > 0$

$$ \sum_{j'=1}^{N} \mathbb{1}_{|k_j - k_{j'}| \le a} + \mathbb{1}_{|\tilde{k}_j - k_{j'}| \le a} + \mathbb{1}_{|k_j - \tilde{k}_{j'}| \le a} + \mathbb{1}_{|\tilde{k}_j - \tilde{k}_{j'}| \le a} \le c_3 \left( N \left( \frac{a}{|k_0|} \right)^{d-1} + 1 \right). $$

The $k_0$-allowed assumption on the frequency pairs $(k_j, \tilde{k}_j)_{1 \le j \le N}$ ensures sufficient separation between frequencies, as stated in point 1. This guarantees the required decorrelation between single-taper estimators evaluated at distinct frequencies. Point 2 ensures that the frequencies are well spread and do not form clusters. For instance, in dimension 2, the $k_0$-allowed assumption is satisfied by choosing $\{(k_j, \tilde{k}_j)\}_{j \in [N]} = C(|k_0|)$, where

$$ C(q) := \{q(\cos(\frac{2\pi j}{N}), \sin(\frac{2\pi j}{N})), q(\sin(\frac{2\pi j}{N}), \cos(\frac{2\pi j}{N}))\}_{j \in [N]}, \quad N \ge 1. \quad (26) $$

With this choice, the symmetry condition $S(k_j) = S(\tilde{k}_j)$ for all $j \in [N]$ is satisfied whenever the structure factor is symmetric under coordinate exchange, i.e. $S(k_1, k_2) = S(k_2, k_1)$ for all $(k_1, k_2) \in \mathbb{R}^2$. In practice, it may be desirable to relax this strict symmetry requirement in order to increase the number of frequency pairs. A natural compromise is to allow pairs $(k_j, \tilde{k}_j)$ such that $S(k_j)$ is close to $S(k_0)$, though not necessarily equal. This enlargement promotes greater decorrelation and reduces the variance of the criteria $\hat{\mathcal{R}}(I)$. In particular, one may consider $\bigcup_{q \in Q} C(q)$, for a finite set $Q$ of positive numbers with $q$ close to $|k_0|$.

The following theorem, proved in Appendix C, shows that under suitable conditions, the taper index set selected by minimizing the empirical criterion $\hat{\mathcal{R}}$ achieves, asymptotically, near-optimal risk for estimating the structure factor of a 1/2-independent thinning of $\Phi$.

**Theorem 6.1.** Assume that $d \ge 2$. Let $\Phi$ be a stationary point process satisfying Assumption 3.1 with structure factor $S \in \Theta(\beta, L)$ where $\beta \le 2$ and $L \in (0, \infty)$. Assume that for all $2 \le m \le 8$: $|\gamma_{\text{red}}^{(m)}|(\mathbb{R}^{d(m-1)}) < \infty$. Let $k_0 \in \mathbb{R}^d$ be a target frequency. Consider $N \ge 1$ $k_0$-allowed frequencies $(k_j, \tilde{k}_j)_{1 \le j \le N}$ in $\mathbb{R}^d$. We denote $A_W = |k_0||W|^{1/(2\beta+d)}$. Suppose that $N \ge A_W^{2d(d-1)/(d+1)}$ and that $A_W \to \infty$ as $|W| \to \infty$. We consider

$$ \mathcal{I} := \left\{ I \mid c_5(|W|^{\frac{2\beta}{2\beta+d}}/A_W^{d(d-1)/(d+1)-\eta}) \vee \log(|W|)^{2d/(3\theta)+\eta} \le |I| \le c_6|W|^{\frac{2\beta}{2\beta+d}} \right\}, $$

where $0 < c_5 \le c_6 < \infty$, $\eta > 0$, $\theta \in (0, 2/3)$ and the tapers index sets $I$ are of the form $\{i \in \mathbb{N}^d | |i|_{\infty} < i_{\max}\}$ for some $i_{\max} \ge 1$. Assume that the dilation factor $r$ of the tapers $(r^{-d/2}\psi_i(\cdot/r))_{i \in I}$ corresponding to $I \in \mathcal{I}$ satisfies

$$ r = R/\sqrt{2i_{\max} + i_{\max}^{1/3+\theta}}. \quad (27) $$

Assume that the tapers index set $I_p$ of the pilot estimator belongs to $\mathcal{I}$ and satisfies $|I_p| = o(|W|^{2\beta/(2\beta+d)})$ as $|W| \to \infty$. Consider $\hat{I} := \arg \min_{I \in \mathcal{I}} \hat{\mathcal{R}}(I)$, where $\hat{\mathcal{R}}(I)$ is defined in (25). Then, as $|W| \to \infty$ we have

$$ \frac{1}{N} \sum_{j=1}^{N} \mathbb{E}[(\hat{S}_{\hat{I}}(k_j) - S_{\frac{1}{2}}(k_j))^2] \le \min_{I \in \mathcal{I}} \left\{ \frac{1}{N} \sum_{j=1}^{N} \mathbb{E}[(\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))^2] \right\} + o(|\mathcal{I}|^{\frac{1}{2}} |W|^{\frac{2\beta}{2\beta+d}}), $$

where $S_{\frac{1}{2}} = 1/2 + S/2$.
Before turning to the assumptions of the theorem, we first clarify the guarantee it provides. When the frequency pairs are chosen so that $S(k_j) = S(k_0)$, the result ensures that the cross-validation method selects asymptotically the optimal number of tapers in the sense that it minimizes the total MSE of the estimators $\hat{S}_I(k_j)$ of $S_{\frac{1}{2}}(k_0)$:

$$ \text{tMSE}(k_0) : I \in \mathcal{I} \mapsto \sum_{j=1}^{N} \mathbb{E}[(\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_0))^2]. \quad (28) $$

As discussed after Definition 5, it may be useful to consider a weaker requirement than $S(k_j) = S(k_0)$. For example, suppose that the target frequencies are taken from $\cup_{q \in Q} C(q)$, where $Q$ is a finite set of positive numbers and $C(q)$ is defined in (26). If, for each $q \in Q$, there exists a frequency $k^q$ such that $S(k) = S(k^q)$ for all $k \in C(q)$, then the theorem provides a guarantee for the sum, over $q$, of the total mean squared errors tMSE($k^q$) defined in (28). Provided that for each $q \in Q$, $S(k^q)$ does not differ too much from $S(k_0)$, this still provides a local guarantee.

Concerning the assumption of the previous theorem, the quantity $A_W$ already appeared in Theorem 5.2 as a threshold under which correlation between tapers associated to opposite frequencies becomes non negligible. The condition $A_W \to \infty$, together with the $k_0$-allowed assumption, ensures asymptotic decorrelation between frequency pairs by forcing the frequencies $(k_j, \tilde{k}_j)_{1 \le j \le N}$ to be sufficiently far from zero and well-spread. As a consequence, we expect the data-driven procedure to perform better at larger frequencies $k_0$. Regarding the taper index sets $I \in \mathcal{I}$, the lower bound on their cardinality $|I|$ prevents the selection of estimators with excessively high variance, while the upper bound excludes estimators with large bias. Lastly, the condition $|I_p| = o(|W|^{2\beta/(2\beta+d)})$ ensures that the pilot estimator $\hat{S}_p$ has negligible bias.

The remainder term $o(|\mathcal{I}|^{1/2}|W|^{2\beta/(2\beta+d)})$ in the previous theorem is negligible compared to the minimax rate established in Theorem 4.1, provided that the number of candidate taper index sets $|\mathcal{I}|$ remains bounded as $|W| \to \infty$. The factor $|\mathcal{I}|^{1/2}$ arises from the fact that the proof of Theorem 6.1 relies on controlling the second moment of a suitable transformation of the criterion $\hat{R}(I)$ (see Section C).

## 6.2 Numerical study

In this section we assess the practical performance of the data-driven criterion.

### 6.2.1 Description of the models

We consider the models listed in Table 1, all studied in dimension $d=2$ over the observation window $W = [-20, 20]^2$, with the parameters specified in the second column. These models were chosen because their structure factors exhibit distinct levels of smoothness and illustrate either clustering (Thomas, Matérn, and exponential LGCP) or repulsive behavior (Ginibre, Bessel DPP, and 1/2-perturbed lattice). For the parameters listed in Table 1, the Thomas cluster, Matérn cluster, and 1/2-perturbed lattice processes all have unit intensity, resulting in about 1600 points within the observation window. The exponential LGCP has intensity $\exp(0.5^2) \approx 1.28$, producing approximately 2000 points, whereas the Ginibre process (intensity $1/\pi$) and the Bessel DPP (intensity 0.3) each yield about 500 points.

The Thomas and Matérn cluster point processes are Neyman-Scott models with intensity field $\Lambda(x) = \alpha \sum_{y \in PP_\mu} p(y-x)$, where $\alpha > 0$ and $PP_\mu$ is a stationary Poisson point process of intensity $\mu$ [21]. For the Thomas, $p$ is the density of a zero-mean Gaussian distribution with variance $\sigma^2$. For the Matérn, $p$ is the density of a uniform random variable over $B(0, \rho)$. The Ginibre and the Bessel
<table><caption>Table 1: Models considered and corresponding structure factors.</caption><thead><tr><th>Model</th><th>Parameter(s)</th><th>Structure factor</th></tr></thead><tbody><tr><td>Thomas cluster</td><td>(α, σ², μ) = (5, 0.5², 0.2)</td><td>1 + 5 exp(-|k|^2/4)</td></tr><tr><td>Matérn cluster</td><td>(α, ρ, μ) = (5, 1.5, 0.2)</td><td>1 + 5(J₁(1.5r|k|)/(1.5r|k|))²</td></tr><tr><td>Exponential LGCP</td><td>(α, σ², μ) = (1, 0.5², 0)</td><td>1 + e<sup>0.25</sup> ∫<sub>R<sup>2</sup></sub>(e<sup>0.5e<sup>-|x|</sup></sup> - 1)e<sup>-ik·x</sup>dx</td></tr><tr><td>Ginibre</td><td>—</td><td>1 - exp(-|k|^2/4)</td></tr><tr><td>Bessel DPP</td><td>(ρ, α) = (0.3, 1)</td><td>1 - 0.3 1<sub>|k|≤4</sub> [2 arccos(<sup>|k|</sup>/<sub>4</sub>) - <sup>|k|</sup>/<sub>2</sub>√(1 - (<sup>|k|</sup>/<sub>4</sub>)<sup>2</sup>)]</td></tr><tr><td>1/2-perturbed lattice</td><td>—</td><td>1 - sinc(<sup>k₁</sup>/<sub>2</sub>)<sup>2</sup> sinc(<sup>k₂</sup>/<sub>2</sub>)<sup>2</sup> exp(-√|k/2|)</td></tr></tbody></table>

DPP are determinantal point processes. The kernel of the Ginibre is $(x, y) \mapsto e^{-x\bar{y}-|x|^2/2-|y|^2/2}/\pi$, where $\mathbb{R}^2$ is identified with $\mathbb{C}$ [30], and the one of the Bessel DPP is $(x, y) \mapsto 2\rho J_1(2|x-y|/\alpha)/(2|x-y|/\alpha)$, where $\rho > 0$, $\alpha > 0$ satisfy $\rho < (\pi\alpha^2)^{-1}$ and $J_1$ is the Bessel function of the first kind of order one. Finally, the 1/2-perturbed lattice [47] is defined as the point process $\{x+V_x+U_x+U \mid x \in \mathbb{Z}^2\}$, where $(U, (U_x)_{x \in \mathbb{Z}^2})$ are i.i.d. uniformly distributed on $[-1/2, 1/2]^2$ and the $(V_x)_{x \in \mathbb{Z}^2}$ are i.i.d. random variables with a symmetric 1/2-stable distribution, whose characteristic function is given by $k \mapsto \exp(-\sqrt{|k/2|}/2)$ for $k \in \mathbb{R}^2$.

### 6.2.2 Choice of the hyper-parameters of the cross-validation

To implement the multitaper estimator, we use Hermite tapers and substitute the standard estimator $\tilde{\lambda} = |\Phi \cap W|/|W|$ for the intensity (see Appendix E). We consider a family of nodes

$$ \{k_n = q_n(\cos(2\pi/3), \sin(2\pi/3)) \mid q_n \in \{0, 0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5\}\}. $$

For the estimation of $S$ at these nodes, we use the data-driven procedure introduced in Section 6 to select the number of tapers, as described below. Then for the estimation of $S$ at an arbitrary frequency $k_0$, we choose the number of tapers by interpolating the selected number of tapers between the nodes of closest norm.

For the cross-validation procedure at a given node $k_n$, we choose:

$$ \mathcal{I} := \{I_{i_{\max}} \mid i_{\max} \le 25\}, \quad \text{where } I_{i_{\max}} := \{(i_1, i_2) \in \mathbb{N}^2 \mid i_1 < i_{\max}, i_2 < i_{\max}\}. $$

For each taper set $I_{i_{\max}}$, we take $r = R/\sqrt{2i_{\max} + i_{\max}^{2/3}}$ for the dilation factor. The tapers set for the pilot estimator is fixed as $I_p = \{(i_1, i_2) \in \mathbb{N}^2 \mid i_1 < 8, i_2 < 8\}$ for the Thomas/Matérn cluster and the exponential LGCP. For the Ginibre, the Bessel DPP and the 1/2-perturbed lattices, where $S$ takes smaller values near the origin, we instead use a smaller taper set $I_p = \{(i_1, i_2) \in \mathbb{N}^2 \mid i_1 < 2, i_2 < 2\}$. We consider the frequency pairs $\{(k_{j_1,j_2}, \tilde{k}_{j_1,j_2})\}$ given by

$$ \{(q_{j_1}(\cos(\frac{2\pi j_2}{n_2}), \sin(\frac{2\pi j_1}{n_2})), q_{j_1}(\sin(\frac{2\pi j_2}{n_2}), \cos(\frac{2\pi j_1}{n_2})))\}_{j_1,j_2 \in [n_1] \times [n_2]}, $$

where $n_1 \in \{4, 5\}$ and $n_2$ is chosen such that for all $j_2 \in \{0, \dots, n_2-1\}$, $|k_{j_1,j_2} - k_{j_1,j_2+1}| = 0.02$ if $|k_n| \le 1$ and $|k_{j_1,j_2} - k_{j_1,j_2+1}| = 0.1$ if $|k_n| \ge 1$. Moreover, for $k_n \ne 0$, we take $q_{j_1} \in |k_n| + \{-0.2, -0.1, 0, 0.1, 0.2\}$, and for $k_n = 0$, we remove the frequency zero from this last set and consider $q_{j_1} \in \{-0.2, -0.1, 0.1, 0.2\}$.
### 6.2.3 Numerical results

Figure 1 presents a qualitative study of the performance of the multitaper estimator, using the local cross-validation procedure to select the number of tapers (see the red curves). The blue curves correspond to the kernel estimator considered in [83], using the cross-validation approach introduced in [24].

The use of the local cross-validation method for selecting the number of tapers is advantageous, since the optimal number of tapers can vary widely across frequencies and cannot be reliably inferred by merely inspecting the curves for different taper numbers. For example, in simulations not shown here (but available upon request), the optimal number of tapers for the Ginibre point process varies from lower than 10 at low frequencies to about 400 at high frequencies. For all models presented in Figure 1, at frequencies that are not too small, the data-driven procedure (corresponding to the red curves) consistently selects an appropriate number of tapers, as predicted by Theorem 6.1. In the more challenging low-frequency regime, the procedure remains effective for repulsive point processes (Ginibre, Bessel DPP and 1/2-perturbed lattice), but it encounters more difficulty for clustered point processes (Thomas, Matérn cluster and exponential LGCP). Finally, as shown by the repulsive models, the estimator achieves higher accuracy when the structure factor $S$ exhibits greater Hölder regularity $\beta$, in agreement with theoretical predictions.

Regarding the comparison with the kernel estimator of $S$ considered in [83] and [24] (see the blue curves in Figure 1), we observe that, for clustering point processes, the multitaper estimator with a locally selected number of tapers exhibits lower bias at low frequencies than the kernel estimator. However, it fluctuates more, which, in the case of the Thomas and Matérn cluster point processes, leads to larger mean squared errors. For the LGCP, the reduction in bias at low frequencies compensates for the increase in fluctuations, resulting in smaller mean squared errors overall. At larger frequencies for these three clustering models, the multitaper and kernel estimators perform comparably. For the other repulsive models considered, however, the multitaper method with a locally data-driven selection of the number of tapers outperforms the kernel estimator having a globally selected bandwidth, across both low and high frequencies, highlighting its advantages in these scenarios.

## 7 Proofs of the minimax lower bounds

In this section, we prove Theorems 4.1 and 4.2. The first key ingredient is Lemma 7.1, which provides an upper bound on the $\chi^2$ divergence between an arcsin-Gaussian Cox point process and a unit-intensity stationary Poisson.

**Definition 6.** Let $P$ and $Q$ be two probability measures. The $\chi^2$ divergence between $P$ and $Q$ is defined by:

$$ \chi^2(P, Q) := \int \left( \frac{dP}{dQ} \right)^2 dQ - 1 $$

if $P$ is absolutely continuous with respect to $Q$ and $+\infty$ otherwise.

Considering the $\chi^2$ divergence is particularly useful in our context, as its quadratic structure enables a symmetrization argument. It allows for reducing the proof of Lemma 7.1 to bounding the exponential moment of an integral of a random field.

**Lemma 7.1.** Let $\Phi_{PPP}$ be a Poisson point process of intensity 1. Let $\sigma \in (0, 1)$, $\rho > 0$ and $(N(x))_{x \in \mathbb{R}^d}$ be a stationary centered Gaussian field with covariance function $c(x) = \sigma^2 c_0(x/\rho)$ where $c_0$ has support included in $B(0, 1/2)$ and $0 \le c_0 \le 1$. We consider a Cox point process $\Phi_{Cox}$ with
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 1: Each model is simulated 500 times in the observation window $[-20, 20]^2$. The Hölder regularity $\beta$ of the corresponding structure factors are indicated between parenthesis. The dashed black curve represents the theoretical structure factor. The red color corresponds to the multitaper method, with the number of tapers locally selected via the cross-validation procedure of Section 6. The blue color corresponds to the kernel estimator [83], with the cross-validation method introduced in [24]. For both methods, the solid line represents the average of the structure factor estimates, while the shaded region corresponds to their 0.05–0.95 empirical quantiles.</figcaption></figure>
intensity field $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$ and an observation window $W = [-R, R]^d$ with $R > 0$. Assume that $8^d \sigma^4 \rho^{2d} < 1$ and $R/\rho \in \mathbb{N}$. Then

$$ \chi^2 := \chi^2(\Phi_{Cox} \cap W, \Phi_{PPP} \cap W) \le (1 - 8^d \sigma^4 \rho^{2d})^{-|W|/(2^d \rho^d)} - 1 $$

*Proof.* The density of $\Phi_{Cox} \cap W$ with respect to $\Phi_{PPP} \cap W$ is given by

$$ \mathbb{E}_N \left[ \exp \left( - \int_W \sin(N(x)) dx \right) \prod_{y \in \Phi_{PPP} \cap W} (1 + \sin(N(y))) \right], $$

where the expectation $\mathbb{E}_N$ is with respect to the Gaussian field $(N(x))_{x \in \mathbb{R}^d}$ [55]. Accordingly, the $\chi^2$ divergence writes

$$ \chi^2 = \mathbb{E}_P \left[ \mathbb{E}_N \left[ \exp \left( - \int_W \sin(N(x)) dx \right) \prod_{y \in \Phi_{PPP} \cap W} (1 + \sin(N(y))) \right]^2 \right] - 1, $$

where $\mathbb{E}_P$ denotes the expectation with respect to the Poisson point process $\Phi_{PPP} \cap W$. We first simplify the nested expectations by introducing $(\tilde{N}(x))_{x \in \mathbb{R}^d}$ an independent copy of $(N(x))_{x \in \mathbb{R}^d}$. Using the Fubini's theorem and the expression for the generating function of the Poisson point process, we obtain

$$ \begin{aligned} & \mathbb{E}_P \left[ \mathbb{E}_N \left[ \exp \left( - \int_W \sin(N(x)) dx \right) \prod_{y \in \Phi_{PPP} \cap W} (1 + \sin(N(y))) \right]^2 \right] \\ &= \mathbb{E}_P \left[ \mathbb{E}_{N, \tilde{N}} \left[ e^{-\int_W (\sin(N(x)) + \sin(\tilde{N}(x))) dx} \prod_{y \in \Phi_{PPP} \cap W} (1 + \sin(N(y)))(1 + \sin(\tilde{N}(y))) \right] \right] \\ &= \mathbb{E}_{N, \tilde{N}} \left[ e^{-\int_W (\sin(N(x)) + \sin(\tilde{N}(x))) dx} e^{\int_W ((1+\sin(N(x)))(1+\sin(\tilde{N}(x))) - 1) dx} \right] \\ &= \mathbb{E}_{N, \tilde{N}} \left[ \exp \left( \int_W \sin(N(x)) \sin(\tilde{N}(x)) dx \right) \right]. \end{aligned} $$

The next step is to decompose $W$ as an union of families of cubes such that the cubes of each family are well-enough separated to exploit the decorrelation of the fields $N$ and $\tilde{N}$. Denote, for $\tau \in \mathbb{Z}^d$ and $\varepsilon \in \{0, 1\}^d$, $C_{\varepsilon, \tau} := \prod_{i=1}^d [(2\tau_i + \varepsilon_i)\rho, (2\tau_i + 1 + \varepsilon_i)\rho]$. We decompose $W$ as

$$ W = \bigsqcup_{\varepsilon \in \{-1, 1\}^d} \bigsqcup_{\tau \in \mathbb{Z}^d} C_{\varepsilon, \tau} \cap W $$

Note that the union is finite. Hölder's inequality yields the following bound

$$ \begin{aligned} \chi^2 &\le \prod_{\varepsilon \in \{-1, 1\}^d} \mathbb{E}_{N, \tilde{N}} \left[ \exp \left( 2^d \int_{\bigsqcup_{\tau \in \mathbb{Z}^d} C_{\varepsilon, \tau} \cap W} \sin(N(x)) \sin(\tilde{N}(x)) dx \right) \right]^{2-d} - 1 \\ &= \prod_{\varepsilon \in \{-1, 1\}^d} \mathbb{E}_{N, \tilde{N}} \left[ \exp \left( 2^d \sum_{\tau \in \mathbb{Z}^d; C_{\varepsilon, \tau} \cap W \ne \emptyset} X_{\varepsilon, \tau} \right) \right]^{2-d} - 1, \end{aligned} $$
with $X_{\varepsilon, \tau} = \int_{C_{\varepsilon, \tau}} \sin(N(x)) \sin(\tilde{N}(x)) dx$. Recall that the covariance function $c$ of $N$ and $\tilde{N}$ has support included in $B(0, \rho/2)$. For $x_{\tau} \in C_{\varepsilon, \tau}$ and $x_{\tau'} \in C_{\varepsilon, \tau'}$ with $\tau \neq \tau'$, we have $|x_{\tau} - x_{\tau'}| \geq \rho$, so $c(x_{\tau} - x_{\tau'}) = 0$. The random integrals $(X_{\varepsilon, \tau})_{\tau \in \mathbb{Z}^d}$ are thus independent. Accordingly,

$$ \chi^2 \leq \prod_{\varepsilon \in \{-1, 1\}^d} \prod_{\tau \in \mathbb{Z}^d, C_{\varepsilon, \tau} \cap W \neq \emptyset} \underbrace{\mathbb{E}_{N, \tilde{N}} \left[ \exp \left( 2^d \int_{C_{\varepsilon, \tau} \cap W} \sin(N(x)) \sin(\tilde{N}(x)) dx \right) \right]}_{=: A_{\varepsilon, \tau}} \quad 2^{-d} - 1. $$

The previous argument have controlled the correlations of the intensity field $1 + \sin(N)$. The remaining of the proof concerns its amplitude. Expanding the exponential function in power series and using the Fubini's theorem, we get

$$ \begin{aligned} A_{\varepsilon, \tau} &= \sum_{n \ge 0} \frac{(2^d)^n}{n!} \mathbb{E} \left[ \left( \int_{C_{\varepsilon, \tau} \cap W} \sin(N(x)) \sin(\tilde{N}(x)) dx \right)^n \right] \\ &= \sum_{n \ge 0} \frac{(2^d)^n}{n!} \int_{(C_{\varepsilon, \tau} \cap W)^n} \mathbb{E} \left[ \prod_{i=1}^n \sin(N(x_i)) \right]^2 dx_1 \dots dx_n \\ &= \sum_{n \ge 0} \frac{4^{dn}}{(2n)!} \int_{(C_{\varepsilon, \tau} \cap W)^{2n}} \mathbb{E} \left[ \prod_{i=1}^{2n} \sin(N(x_i)) \right]^2 dx_1 \dots dx_{2n}. \end{aligned} $$

In the last line, we used $\mathbb{E}[N(x)] = 0$. Using $|\sin(x)| \le |x|$ for all $x \in \mathbb{R}$, the Hölder-inequality and the stationarity of $N$, we obtain

$$ \begin{aligned} A &\le \sum_{n \ge 0} \frac{4^{dn}}{(2n)!} \int_{(C_{j,\tau} \cap W)^{2n}} \prod_{i=1}^{2n} \mathbb{E} [\sin(N(x_i))^{2n}]^{1/n} dx_1 \dots dx_{2n} \\ &\le \sum_{n \ge 0} \frac{4^{dn} |C_{j,\tau} \cap W|^{2n}}{(2n)!} \mathbb{E} [N(0)^{2n}]^2 \\ &\le \sum_{n \ge 0} \frac{(4^d \rho^d)^{2n}}{(2n)!} \left( \sigma^{2n} \frac{(2n)!}{n!2^n} \right)^2 \\ &\le \sum_{n \ge 0} (8^d \sigma^4 \rho^{2d})^n = (1 - 8^d \sigma^4 \rho^{2d})^{-1}, \end{aligned} $$

since $\binom{2n}{n} \le 2^{2n}$ and $8^d \sigma^4 \rho^{2d} < 1$. This yields

$$ \chi^2(\Phi_{Cox} \cap W, \Phi_{PPP} \cap W) \le (1 - 8^d \sigma^4 \rho^{2d})^{-2-d} \sum_{\varepsilon \in \{-1, 1\}^d} \sum_{\tau \in \mathbb{Z}^d} \mathbb{1}_{C_{\varepsilon, \tau} \cap W \neq \emptyset}. $$

To conclude, we use $R/\rho \in \mathbb{N}$ to get $\sum_{\varepsilon \in \{-1, 1\}^d} \sum_{\tau \in \mathbb{Z}^d} \mathbb{1}_{C_{\varepsilon, \tau} \cap W \neq \emptyset} = |W|/\rho^d$. $\square$

The next crucial result is Lemma 7.2, which addresses the behavior of the total variation of the reduced cumulant measure of the Cox processes of Lemma 7.1. This result is important, as it ensures that the constraint $\max_{m \in \{2, 3, 4\}} |\gamma_{red}^{(m)}| (\mathbb{R}^{d(m-1)}) \le M$ is satisfied in Theorem 4.1. It proof consists in transposing the property of the vanishing cumulants of order larger than three for Gaussian fields to these specific Cox point processes. This lemma motivates our choice of considering an underlying intensity field of the form $(1 + \sin(\sigma^2 N_0(x)))_{x \in \mathbb{R}^d}$, where $(N_0(x))_{x \in \mathbb{R}^d}$ is a unit-variance Gaussian field. Using $(1 + \sigma^2 \sin(N_0(x)))_{x \in \mathbb{R}^d}$ would have simplified the proof of the previous lemma, but it would not allow us to control higher-order cumulants in the low-regularity setting $\beta < \frac{d}{2}$.
**Lemma 7.2.** Let $\sigma \in (0, 1)$ and $\rho > 0$. Let $(N(x))_{x \in \mathbb{R}^d}$ be a stationary centered Gaussian field with covariance function $c(x) = \sigma^2 c_0(x/\rho)$ where $c_0$ has support included in $B(0, 1/2)$ and $0 \le c_0 \le 1$. We consider a Cox point process $\Phi_{Cox}$ with intensity field $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$. Then, the reduced factorial cumulants measures of $\Phi_{Cox}$ satisfies

$$|\gamma_{\text{red}}^{(3)}(\mathbb{R}^{2d})| \le K\sigma^4\rho^{2d}, \quad |\gamma_{\text{red}}^{(4)}(\mathbb{R}^{3d})| \le K\sigma^6\rho^{3d},$$

where $K < \infty$ is a numerical constant.

*Proof.* Let $m \ge 2$. The $m$-th order cumulant measure of $\Phi_{Cox}$ is given by [32]

$$\gamma^{(m)}(x_1, \dots, x_m) = \kappa_m (1 + \sin(N(x_1)), \dots, 1 + \sin(N(x_m))).$$

Using the multi-linearity of the cumulants, we get:

$$\gamma^{(m)}(x_1, \dots, x_m) = \kappa_m (\sin(N(x_1)), \dots, \sin(N(x_m))).$$

Therefore, $\gamma_{\text{red}}^{(m)}(y_1, \dots, y_{m-1}) = \kappa_m (\sin(N(y_1)), \dots, \sin(N(y_{m-1})), \sin(N(0)))$. The heuristic for the remainder of the proof is as follows. Assume that $m = 3$. Then, as $\sigma \to 0$, we have

$$\gamma_{\text{red}}^{(3)}(y_1, y_2) = \kappa_3 (\sin(N(y_1)), \sin(N(y_2)), \sin(N(0))) \\ \approx \kappa_3 \left( 1 + N(y_1) + \frac{1}{2}N(y_1)^2, 1 + N(y_2) + \frac{1}{2}N(y_2)^2, 1 + N(0) + \frac{1}{2}N(0)^2 \right).$$

Then, using the multi-linearity of the cumulants, we obtain a sum of terms, whose main order with respect to $\sigma \to 0$ is: $\sigma^{1+1+1}\kappa_3(N_1, N_2, N_3) + \sigma^{1+1+2}\kappa_3(N_1, N_2, N_3^2)$, where $(N_1, N_2, N_3)$ is Gaussian. Accordingly, the first order term $\kappa_3(N_1, N_2, N_3)$ vanishes, and the scaling $\sigma^4$ appears. Below, we make the previous heuristic rigorous.

We only consider the case $m = 4$, the case $m = 3$ being similar. Accordingly, in the following, we fix $m = 4$ and denote $y_4 = 0$. Note that if $[4] = L \sqcup L'$ with two disjoint non-empty sets $L$ and $L'$ such that the family of random variables $(\sin(N(y_l)))_{l \in L}$ and $(\sin(N(y_l)))_{l \in L'}$ are independent, then $\gamma_{\text{red}}^{(4)}(y_1, \dots, y_{m-1}) \equiv 0$ (see, e.g., [61], Section 3.1. Property (ii)). Since the Gaussian field $(N(x))_{x \in \mathbb{R}^d}$ has covariance function having support included in $B(0, \rho/2)$, it implies that if $\gamma_{\text{red}}^{(m)}(y_1, \dots, y_{m-1}) \ne 0$ then there exists a cycle $\pi : [4] \mapsto [4]$ such that for all $s \in [d]$ and $j \in [4]$, then $|y_j \cdot e_s - y_{\pi(j)} \cdot e_s| \le \rho$. Accordingly,

$$\begin{aligned} |\gamma_{\text{red}}^{(4)}(\mathbb{R}^{3d})| & \le \sum_{\text{cycle } \pi: [4] \mapsto [4]} \int_{\mathbb{R}^{3d}} \prod_{i=1}^{d} \prod_{j=1}^{4} \mathbb{1}_{|y_j \cdot e_i - y_{\pi(j)} \cdot e_i| \le \rho_i} dy_1 dy_2 dy_2 \sup_{z \in \mathbb{R}^{3d}} |\gamma_{\text{red}}^{(4)}(z)| \\ & = (4-1)! \left(\prod_{i=1}^{d} 2\rho_i\right)^3 \sup_{z \in \mathbb{R}^{3d}} |\gamma_{\text{red}}^{(4)}(z)|. \end{aligned} \quad (29)$$

It remains to bound the supremum of $h$. Note that $N = \sigma N_0$ in distribution, where $N_0$ is a Gaussian field with covariance function $c_0$. Moreover, Rolle's theorem implies that for all $j \in [4]$ there exists a random variable $Y_j \in \mathbb{R}$, such that $\sin(N(y_j)) = \sigma N_0(y_j) - \sigma^3 N_0(y_j)^3 \cos(Y_j)/6$. Using the multi-linearity of the cumulants we get

$$\gamma_{\text{red}}^{(4)}(y_1, y_2, y_3) = \sum_{P \subset [4]} \sigma^{|P^c|+3|P|} \left(-\frac{1}{6}\right)^{|P|} \kappa_4(N_0(y_l), \cos(Y_{l'})N_0(y_{l'})^3; l \in P^c, l' \in P),$$
where the sum runs over all the subsets $P$ of $[4]$. Since the fourth order joint cumulants of a Gaussian field are zero, the sums runs in fact over the non empty subsets $P$. Hence, $|P^c|+3|P| = 4+2|P| \ge 6$. Using $\sigma \le 1$, we have

$$|\gamma_{\text{red}}^{(4)}(y_1, y_2, y_3)| \le \sigma^6 \sum_{P \subset [4]} |\kappa_4(N_0(y_l), \cos(Y_{l'}) N_0(y_{l'})^3; l \in P^c, l' \in P)| \quad (30)$$

We denote $Z_l = N_0(x_l)^3$ is $l \in P$ and $Z_l = N_0(x_l)$ if $l \in P^c$. By Hölder-inequality, we have

$$|\kappa_4(N_0(y_l), \cos(Y_{l'}) N_0(y_{l'})^3; l \in P^c, l' \in P)| \le \sum_{\pi \subset \Pi[4]} (|\pi| - 1)! \prod_{b \in \pi} \prod_{l \in b} \mathbb{E} [|Z_l|^{|b|}]^{1/|b|},$$

which is a numerical constant. Recalling (29) and (30), this yields $|\gamma_{\text{red}}^{(4)}|(\mathbb{R}^{3d}) \le K\sigma^6\rho^{3d}$. With similar arguments, we obtain $|\gamma_{\text{red}}^{(3)}|(\mathbb{R}^{2d}) \le K\sigma^4\rho^{2d}$. $\square$

The next lemma derives the structure factor and the pair correlation function of the arcsin-Cox processes considered in this section.

**Lemma 7.3.** Let $(N(x))_{x \in \mathbb{R}^d}$ be a stationary centered Gaussian field with bounded covariance function $c$ having a compact support. We consider a Cox point process $\Phi_{\text{Cox}}$ with intensity field $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$. Then, the pair correlation function $g_c$ of $\Phi_{\text{Cox}}$ is

$$\forall x \in \mathbb{R}^d, g_c(x) = 1 + \frac{e^{-c(0)}}{2} (e^{c(x)} - e^{-c(x)})$$

Moreover, the structure factor $S_c$ of $\Phi_{\text{Cox}}$ is

$$\forall k \in \mathbb{R}^d, S_c(k) = 1 + \frac{e^{-c(0)}}{2} \int_{\mathbb{R}^d} (e^{c(x)} - e^{-c(x)}) e^{-ik \cdot x} dx.$$

*Proof.* The intensity of $\Phi_{\text{Cox}}$ is $\mathbb{E}[1 + \sin(N(0))] = 1$. Moreover, its pair correlation function $g_c$ is [32]:

$$g_c(x) = \mathbb{E}[(1 + \sin(N(x)))(1 + \sin(N(0)))] \\ = 1 + 2\mathbb{E}[\sin(N(0))] + \mathbb{E}[\sin(N(x))\sin(N(0))] = 1 + \mathbb{E}[\sin(N(x))\sin(N(0))].$$

The Moivre's formula yields

$$\begin{aligned} \mathbb{E}[\sin(N(x)) \sin(N(0))] &= -\frac{1}{2} \left( \mathbb{E}[e^{i(N(x)+N(0))}] - \mathbb{E}[e^{i(N(x)-N(0))}] \right) \\ &= \frac{e^{-c(0)}}{2} (e^{c(x)} - e^{-c(x)}). \end{aligned}$$

Since $c$ is compactly supported, $g_c - 1 \in L^1(\mathbb{R}^d)$. This concludes the proof, recalling (3). $\square$

We are now in position to prove Theorem 4.1.

*Proof of Theorem 4.1.* Let $W \in \mathcal{W}$. We prove the lower bound when estimating $S$ at frequency zero.

$$r_W := \inf_{\hat{S}} \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \sup_{k_0 \in \mathbb{R}^d} \mathbb{E} \left[ h \left( |\Phi \cap W|^{\frac{\beta}{2\beta+d}} \left| \hat{S}(k_0; \Phi \cap W) - S(k_0) \right| \right) \right] \\ \ge \inf_{\hat{S}} \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \mathbb{E} \left[ h \left( |\Phi \cap W|^{\frac{\beta}{2\beta+d}} \left| \hat{S}(0; \Phi \cap W) - S(0) \right| \right) \right].$$
The newt step is to reduce the proof to the square window setting. Since $\mathcal{W}$ is regular, there exists $\eta > 0$ such that for all $W \in \mathcal{W}$ there exists $R_W$ such that $W \subset [-R_W, R_W]^d$ and $|[-R_W, R_W]^d| \le \eta|W|$. Let $\hat{S} : \mathbb{R}^d \times \text{Conf}(\mathbb{R}^d) \to \mathbb{R}$ be any estimator. Using the fact that $h$ is increasing and $|[-R_W, R_W]^d| \le \eta|W|$, we have

$$ \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \mathbb{E} \left[ h(|W|^{\frac{\beta}{2\beta+d}} |\hat{S}(0; \Phi \cap W) - S(0)|) \right] \\ \geq \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \mathbb{E} \left[ h_{\eta}(|[-R_W, R_W]^d|^{\frac{\beta}{2\beta+d}} |\hat{S}_W(0; \Phi \cap [-R_W, R_W]^d) - S(0)|) \right], $$

where $h_{\eta} := h(\cdot/\eta^{\frac{2\beta}{2\beta+d}})$ and $\hat{S}_W(\cdot; X) := \hat{S}(\cdot; X \cap W)$. Note that in the last line, we used that $W \cap [-R_W, R_W]^d = W$. Accordingly,

$$ r_W \geq \inf_{\hat{S}} \sup_{S \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(S,M)}} \mathbb{E} \left[ h_{\eta} \left( \left( |[-R_W, R_W]^d|^{\frac{\beta}{2\beta+d}} - \hat{S}(0; \Phi \cap [-R_W, R_W]^d) - S(0) \right) \right) \right]. \quad (31) $$

To derive a lower bound for this quantity, we follow the method presented in Chapter 2 of [77], and specifically apply Theorem 2.2 therein.

To alleviate the following computations, we use the notation $R := R_W$. Let $\Phi_{PPP}$ be a Poisson point process of intensity 1. Let $\phi : \mathbb{R}^d \to \mathbb{R}^+$ be a non identically null, infinitely differentiable function with support included in $B(0, 1/4)$. Then, $c_0 = \phi * \phi / \|\phi * \phi\|_{\infty}$ is positive definite and has support included in $B(0, 1/2)$. Moreover $0 \le c_0 \le 1$. Let $\varepsilon \in (0, 1)$ with $1/\varepsilon \in \mathbb{N}$, that will be chosen small enough during the proof. Let

$$ \rho := \varepsilon R / \left[ \frac{R}{R^{d/(d+2\beta)}} \right], \quad \sigma^2 = (R^d)^{-(d+\beta)/(2\beta+d)}. $$

We consider a stationary centered Gaussian field $(N(x))_{x \in \mathbb{R}^d}$ with covariance function $c(x) = \sigma^2 c_0(x/\rho)$ and a Cox point process $\Phi_{Cox}$ with intensity field $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$. In order to apply Theorem 2.2 of [77] we prove that:

1. The structure factor $S_c$ of $\Phi_{Cox}$ is $\Theta(\beta, L)$,

2. $|S_c(0) - 1| \ge A[[-R, R]^d]^{-\beta/(2\beta+d)}$ for some constant $A > 0$ that does not depend on $R$,

3. $\chi^2(\Phi_{PPP} \cap [-R, R]^d, \Phi_{Cox} \cap [-R, R]^d) \le \alpha$ where $\alpha < \infty$ does also not depend on $R$,

4. The total variations of the reduced factorial cumulant measures of $\Phi_{Cox}$ satisfies

$$ \max_{m \in \{2,3,4\}} |\gamma_{red}^{(m)}|(\mathbb{R}^{d(m-1)}) \le M. $$

For the moment, assume that all previous points are satisfied. Then, using (31), we obtain

$$ r_W \geq \inf_{\hat{S}} \max_{j=0,1} \mathbb{E}_{\Phi_j \cap [-R, R]^d} \left[ h_{\eta} \left( \left( |[-R, R]^d|^{\frac{\beta}{2\beta+d}} - \hat{S}(0; \Phi_j \cap [-R, R]^d) - S(0) \right) \right) \right], $$

where $\Phi_0 := \Phi_{PPP}$ and $\Phi_1 := \Phi_{Cox}$. Note that the expectations are with respect to the distributions of $\Phi_j \cap [-R, R]^d$. Then, employing the method of Chapter 2, Section 2 of [77], we obtain $r_W \ge c_1 \inf_{\Psi} \max_{j=0,1} \mathbb{P}_{\Phi_j \cap [-R, R]^d}[\Psi \ne j]$, where the infimum is over all tests $\Psi : \text{Conf}(\mathbb{R}^d) \to \{0, 1\}$, and $c_1$ is a constant that does not depend on $R$. Using Theorem 2.2 of [77], we obtain $r_W \ge c_1 \exp(-\alpha)/2 > 0$, which does not depend on $W$. This yields (6).
Now, we prove that points 1, 2, 3 and 4 are satisfied. We consider the point 1. According to Lemma 7.3, the structure factor $S_c$ of $\Phi_{Cox}$ is

$$ 
\begin{aligned}
S_c(k) &= 1 + \frac{e^{-\sigma^2 c_0(0)}}{2} \int_{\mathbb{R}^d} (e^{\sigma^2 c_0(x/\rho)} - e^{-\sigma^2 c_0(x/\rho)}) e^{-ik \cdot x} dx \\
&= 1 + \rho^d \frac{e^{-\sigma^2 c_0(0)}}{2} \int_{B(0,1)} (e^{\sigma^2 c_0(x)} - e^{-\sigma^2 c_0(x)}) e^{-i\rho k \cdot x} dx.
\end{aligned} 
\quad (32) $$

Since $c_0$ has compact support, the Lebesgue dominated convergence theorem ensures that $S_c$ is infinitely differentiable. Let $l \ge 0$ be the smallest integer that is strictly smaller than $\beta$. In particular, $S_c$ is $l$-times differentiable. Moreover, for all $n = (n_1, \dots, n_d) \in \mathbb{N}^d$ such that $n_1 + \dots + n_d = l$, we have for $(k, k') \in \mathbb{R}^{2d}$:

$$ 
\begin{aligned}
|\partial_n^l S_c(k) - \partial_n^l S_c(k')| &\le \rho^{d+l} \int_{B(0,1)} \frac{1}{2} (e^{\sigma^2 c_0(x)} - e^{-\sigma^2 c_0(x)}) |x|^l |e^{i\rho k \cdot x} - e^{i\rho k' \cdot x}| dx \\
&\le e\rho^{d+l}\sigma^2 \int_{B(0,1)} c_0(x)|x|^l 2^{1-(\beta-l)} \rho^{\beta-l} |k-k'|^{\beta-l} |x|^{\beta-l} dx \\
&\le 2e\rho^{d+\beta}\sigma^2 \left( \int_{B(0,1)} c_0(x)|x|^{\beta} dx \right) |k-k'|^{\beta-l}.
\end{aligned} 
$$

The choices of $\rho$ and $\sigma^2$ ensure:

$$ \rho^{d+\beta}\sigma^2 \le \varepsilon^{d+\beta} \frac{R^{d+\beta}}{(R/R^{d/(2\beta+d)})^{d+\beta}} (R^d)^{-(d+\beta)/(2\beta+d)} = \varepsilon^{d+\beta} \quad (33) $$

By reducing $\varepsilon$ if necessary, we have $S_c \in \Theta(\beta, L)$. Concerning point 2, with (32), we get

$$ 
\begin{aligned}
|S_c(0) - 1| &= \rho^d \frac{e^{-\sigma^2 c_0(0)}}{2} \int_{B(0,1)} (e^{\sigma^2 c_0(x)} - e^{-\sigma^2 c_0(x)}) dx \\
&\ge \rho^d \sigma^2 e^{-c_0(0)} \int_{B(0,1)} e^{\sigma^2 c_0(x)} c_0(x) dx \ge \rho^d \sigma^2 e^{-c_0(0)} \|c_0\|_1.
\end{aligned} 
$$

Note that

$$ \rho^d \sigma^2 \ge \varepsilon^d \frac{R^d}{(R/R^{d/(2\beta+d)} + 1)^{d+\beta}} (R^d)^{-\frac{d+\beta}{2\beta+d}} \ge 2^{-\beta-d} \varepsilon^d (R^d)^{-\frac{\beta}{2\beta+d}}. $$

Subsequently, point 2 holds with $A = 2^{-\beta-d}\varepsilon^d e^{-c_0(0)} \|c_0\|_1$. To verify point 3, we use Lemma 7.1. If $\varepsilon = 1/n_0$ with $n_0 \in \mathbb{N}$ and $n_0 \ge 5$, then its assumptions are satisfied. Indeed, with computation similar to (33), we get $8^d \sigma^4 \rho^{2d} \le 8^d \varepsilon^{2d} (R^d)^{-2\beta/(2\beta+d)} \le 8^d \varepsilon^{2d} < 1$. Hence Lemma 7.1 ensures that

$$ 
\begin{aligned}
\chi^2(\Phi_{PPP} \cap [-R, R]^d, \Phi_{Cox} \cap [-R, R]^d) &\le (1 - 8^d \sigma^4 \rho^{2d})^{-\frac{R^d}{\rho^d}} \\
&\le (1 - 8^d \varepsilon^{2d} (R^d)^{-\frac{2\beta}{2\beta+d}})^{-\frac{(R^d)^{2\beta/(2\beta+d)}}{\varepsilon^d}} \le \sup_{x \in [1, \infty)} (1 - 2d 8^d \varepsilon^{2d}/x)^{-x/\varepsilon^d} := \alpha.
\end{aligned} 
$$
Accordingly, point (3) holds. It remains to verify point 4. Firstly, using sinh($x$) $\le xe^x$ for $x \ge 0$, we have

$$ \begin{aligned} |\gamma_{\text{red}}^{(2)}|(\mathbb{R}^d) &= \frac{e^{-\sigma^2 c_0(0)} \rho^d}{2} \int_{B(0,1)} (e^{\sigma^2 c_0(x)} - e^{-\sigma^2 c_0(x)}) dx \\ &\le \rho^d \sigma^2 e \int_{B(0,1)} c_0(x) dx \le \varepsilon^{2d} e \int_{B(0,1)} c_0(x) dx. \end{aligned} $$

Moreover, according to Lemma 7.2, $|\gamma_{\text{red}}^{(3)}|(\mathbb{R}^{2d}) \le K\sigma^4\rho^{2d}$, $|\gamma_{\text{red}}^{(4)}|(\mathbb{R}^{3d}) \le K\sigma^6\rho^{3d}$, for some numerical constant $K < \infty$. Since,

$$ \sigma^6 \rho^{3d} \le \varepsilon^{3d} (R^d)^{3d/(2\beta+d)} R^{-3(d+\beta)/(2\beta+d)} \le \varepsilon^{3d}, \quad \sigma^4 \rho^{2d} \le \varepsilon^{2d}, $$

then point 4 holds, if $\varepsilon$ is small enough. Theorem 2.2. of [77] concludes the proof. $\square$

**Proof of Theorem 4.2.** Let $W \in W$. Arguing as in the beginning of the proof of Theorem 4.1, we obtain

$$ \begin{aligned} r &:= \inf_{\hat{g}} \sup_{g \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(g,M)}} \sup_{x_0 \in \mathbb{R}^d} \mathbb{E} \left[ h \left( |W|^{\frac{\beta}{2\beta+d}} |\hat{g}(x_0; \Phi \cap W) - g(x_0)| \right) \right] \\ &\ge \inf_{\hat{g}} \sup_{g \in \Theta(\beta, L)} \sup_{\Phi \sim \mathcal{L}_{(g,M)}} \mathbb{E} \left[ h_{\eta} \left( \left( |-R, R|^d \right)^{\frac{\beta}{2\beta+d}} |\hat{g}(0; \Phi \cap [-R, R]^d) - g(0)| \right) \right], \end{aligned} $$

where $\eta$ does not depend on $W$, $h_{\eta} := h(\cdot/\eta^{\beta/(2\beta+d)})$ and $R = R_W > 0$. As in the proof of Theorem 4.1, the result will follow from the method of Chapter 2 of [77].

Let $\Phi_{PPP}$ be a Poisson point process of intensity 1. Let $\phi: \mathbb{R}^d \to \mathbb{R}^+$ be a non identically null, infinitely differentiable function with support included in $B(0, 1/4)$. Then, $c_0 = \phi * \phi / \|\phi * \phi\|_{\infty}$ is positive definite and has support included in $B(0, 1/2)$. Moreover $0 \le c_0 \le 1$. Let $\varepsilon \in (0, 1)$, that will be chosen small enough during in the proof. Consider

$$ \rho := 1 / \left[ R^{d/(d+2\beta)} \right], \quad \sigma^2 = \varepsilon R^{-d\beta/(2\beta+d)}. $$

We consider a stationary centered Gaussian field $(N(x))_{x \in \mathbb{R}^d}$ with covariance function $c(x) = \sigma^2 c_0(x/\rho)$ and a Cox point process $\Phi_{Cox}$ with intensity field $(1 + \sin(N(x)))_{x \in \mathbb{R}^d}$. In order to apply Theorem 2.2 of [77] we prove that:

1. The pair correlation function $g_c$ of $\Phi_{Cox}$ is $\Theta(\beta, L)$,

2. $|g_c(0) - 1| \ge A[|-R, R]^d|^{-\beta/(2\beta+d)}$ for some constant $A > 0$ that does not depend on $R$,

3. $\chi^2(\Phi_{PPP} \cap [-R, R]^d, \Phi_{Cox} \cap [-R, R]^d) \le \alpha$ where $\alpha < \infty$ does also not depend on $R$,

4. The total variations of the reduced factorial cumulant measures of $\Phi_{Cox}$ satisfies

$$ \max_{m \in \{2,3,4\}} |\gamma_{\text{red}}^{(m)}|(\mathbb{R}^{d(m-1)}) \le M. $$

We check point 1. Recall from Lemma 7.3, that

$$ \forall x \in \mathbb{R}^d, g_c(x) = 1 + \frac{e^{-\sigma^2 c_0(0)}}{2} \left( e^{\sigma^2 c_0(x/\rho)} - e^{-\sigma^2 c_0(x/\rho)} \right). $$
Let $n = (n_1, \dots, n_d) \in \mathbb{N}^d$ and $l = n_1 + \dots + n_d$. According to Faà-di-Bruno formula,

$$ \partial_n^l (g_c(x) - 1) = \frac{e^{-\sigma^2 c_0(0)}}{2} \sum_{\pi \in [l]} \left( e^{\sigma^2 c_0(\frac{x}{\rho})} \prod_{b \in \pi} \frac{\sigma^2}{\rho^{|b|}} \partial_b c_0(\frac{x}{\rho}) + (-1)^{|\pi|} e^{-\sigma^2 c_0(\frac{x}{\rho})} \prod_{b \in \pi} \frac{\sigma^2}{\rho^{|b|}} \partial_b c_0(\frac{x}{\rho}) \right), $$

where in the partition $\pi$, the numbers $1, \dots, n_1$ correspond to a differentiation with respect to $x_1$, and $n_1 + 1, \dots, n_1 + n_2$ to $x_2$ and so on. This yields, using $\sigma \le 1$,

$$ \|\partial_n^l (g_c - 1)\|_{\infty} \le e^{\frac{\sigma^2}{\rho l}} \sum_{\pi \in [m]} \prod_{b \in \pi} \|\partial_b c_0\|_{\infty}. $$

Let $l \ge 0$ be the smallest integer that is strictly smaller than $\beta$. Then, for $(x, y) \in \mathbb{R}^{2d}$,

$$ \begin{aligned} |\partial_n^l g_c(x) - \partial_n^l g_c(x)| &\le (2\|\partial_n^l (g_c - 1)\|_{\infty})^{1-(\beta-l)} (\|\nabla \partial_n^l (g_c - 1)\|_{\infty} |x-y|)^{\beta-l} \\ &\le C_1 \left(\frac{\sigma^2}{\rho^l}\right)^{1-(\beta-l)} \left(\frac{\sigma^2}{\rho^{l+1}}\right)^{\beta-l} |x-y|^{\beta-l} = C_1 \frac{\sigma^2}{\rho^{\beta}} |x-y|^{\beta-l}, \end{aligned} $$

for some constant $C_1 < \infty$ that depends only on $c_0$. The choices of $\sigma^2$ and $\rho$ ensure that $c_1\sigma^2/\rho^\beta \le L$ for $\varepsilon$ small enough. This yields point 1. Concerning point 2, we have, using $\sinh(x) \ge x$ for $x \ge 0$,

$$ g_c(0) - 1 = \frac{e^{-\sigma^2 c_0(0)}}{2} \left( e^{\sigma^2 c_0(0)} - e^{-\sigma^2 c_0(0)} \right) \ge \frac{e^{-\varepsilon}}{2} \sigma^2 c_0(0) = A |[-R, R]^d|^{-\beta/(2\beta+d)}, $$

with $A = \exp(-\varepsilon)\varepsilon^2 c_0(0)/2 > 0$. We deduce point 3 from Lemma 7.1:

$$ \begin{aligned} \chi^2(\Phi_{PPP} \cap [-R, R]^d, \Phi_{Cox} \cap [-R, R]^d) &\le (1 - \varepsilon^2 8^d R^{-(d+\beta)/(2\beta+d)})^{\varepsilon R^{(d+\beta)/(2\beta+d)}} \\ &\le \sup_{x \ge 1} (1 - \varepsilon^2 8^d/x)^{\varepsilon x} =: \alpha, \end{aligned} $$

provided that $\varepsilon^2 8^d < 1$. Finally, since $\sigma$ and $\rho$ are smaller than 1, Lemma 7.2 ensures that point 4 is satisfied. This concludes the proof, using Theorem 2.2. of [77]. □

## Acknowledgment

The author sincerely thanks Frédéric Lavancier for thoroughly reviewing the first versions of the manuscript and for the enriching discussions, as well as Bartłomiej Błaszczyszyn for the stimulating discussions surrounding this work.

## References

[1] ANEVSKI, D., AND SOULIER, P. Monotone spectral density estimation. Ann. Statist. 39, 1 (2011), 418–438.

[2] APTEKAREV, A., DEHESA, J., SÁNCHEZ-MORENO, P., AND TULYAKOV, D. Asymptotics of Lp-norms of hermite polynomials and rényi entropy of rydberg oscillator states. Contemp. Math 578 (2012), 19–29.
[3] ASKEY, R., AND WAINGER, S. Mean convergence of expansions in laguerre and hermite series. *American Journal of Mathematics* 87, 3 (1965), 695-708.

[4] BACCELLI, F., BLASZCZYSZYN, B., AND KARRAY, M. *Random Measures, Point Processes, and Stochastic Geometry*. hal-02460214, Inria, 2020.

[5] BARTLETT, M. S. The spectral analysis of point processes. *Journal of the Royal Statistical Society Series B: Statistical Methodology* 25, 2 (1963), 264-281.

[6] BARTLETT, M. S. The spectral analysis of two-dimensional point processes. *Biometrika* 51, 3/4 (1964), 299-311.

[7] BAYRAM, M., AND BARANIUK, R. G. Multiple window time-frequency analysis. In *Proceedings of third international symposium on time-frequency and time-scale analysis (tfts-96)* (1996), IEEE, pp. 173-176.

[8] BELTRATO, K. I., AND BLOOMFIELD, P. Determining the bandwidth of a kernel spectrum estimate. *Journal of time series analysis* 8, 1 (1987), 21-38.

[9] BENTKUS, R. Rate of uniform convergence of statistical estimators of spectral density in spaces of differentiable functions. *Lithuanian Mathematical Journal* 25, 3 (1985), 209-219.

[10] BENTKUS, R., AND RUDZKIS, R. On exponential estimates of the distribution of random variables. *Liet. Mat. Rinkinys* 20, 1 (1980), 15-30.

[11] BISCIO, C. A. N., AND LAVANCIER, F. Brillinger mixing of determinantal point processes and statistical applications. *Electronic Journal of Statistics* 10, 1 (2016), 582-607.

[12] BLACKMAN, R., TUKEY, J., AND TEICHMANN, T. The measurement of power spectra. *Physics Today* 13, 2 (1960), 52-54.

[13] BONAN, S. S., AND CLARK, D. S. Estimates of the hermite and the freud polynomials. *Journal of Approximation Theory* 63, 2 (1990), 210-224.

[14] BRÉMAUD, P. Fourier analysis of stochastic processes. In *Fourier analysis and stochastic processes*. Springer, 2014, pp. 119-179.

[15] BRILLINGER, D. R. The spectral analysis of stationary interval functions. In *Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Theory of Statistics* (1972), vol. 6, University of California Press, pp. 483-514.

[16] BRILLINGER, D. R. Uses of cumulants in wavelet analysis. In *Advanced Signal Processing: Algorithms, Architectures, and Implementations V* (1994), vol. 2296, SPIE, pp. 2-18.

[17] BRILLINGER, D. R. John w. tukey's work on time series and spectrum analysis. *The Annals of Statistics* 30, 6 (2002), 1595-1618.

[18] COWLEY, J. M. *Electron diffraction techniques*, vol. 2. Oxford University Press, 1992.

[19] CRONIE, O., MORADI, M., AND BISCIO, C. A. A cross-validation-based statistical theory for point processes. *Biometrika* 111, 2 (2024), 625-641.

[20] DALEY, D. J., AND VERE-JONES, D. *An introduction to the theory of point processes: volume II: general theory and structure*. Springer Science & Business Media, 2007.
[21] DALEY, D. J., VERE-JONES, D., ET AL. *An introduction to the theory of point processes: volume I: elementary theory and methods*. Springer, 2003.

[22] DAUBECHIES, I. Time-frequency localization operators: a geometric phase space approach. *IEEE transactions on information theory* 34, 4 (2002), 605–612.

[23] DAWSON, B. The significance of accurate structure factors. *Foundations of Crystallography* 25, 1 (1969), 12–29.

[24] DING, Q.-W., YANG, J., AND SHIN, J. Pseudo-spectra of multivariate inhomogeneous spatial point processes. *arXiv preprint arXiv:2502.09948* (2025).

[25] DÖRING, H., JANSEN, S., AND SCHUBERT, K. The method of cumulants for the normal approximation. *Probability Surveys* 19 (2022), 185–270.

[26] DUDLEY, R. M. *Real analysis and probability*, vol. 79. Cambridge University Press, 2002.

[27] EFROMOVICH, S. Data-driven efficient estimation of the spectral density. *Journal of the American Statistical Association* 93, 442 (1998), 762–769.

[28] FLANDRIN, P. Maximum signal energy concentration in a time-frequency domain. In *ICASSP-88., International Conference on Acoustics, Speech, and Signal Processing* (1988), IEEE Computer Society, pp. 2176–2177.

[29] FOLLAND, G. B. *Fourier analysis and its applications*, vol. 4. American Mathematical Soc., 2009.

[30] GINIBRE, J. Statistical ensembles of complex, quaternion, and real matrices. *Journal of Mathematical Physics* 6, 3 (1965), 440–449.

[31] GRAINGER, J. P., RAJALA, T. A., MURRELL, D. J., AND OLHEDE, S. C. Spectral estimation for spatial point processes and random fields. *arXiv preprint arXiv:2312.10176* (2023).

[32] GRANDELL, J. *Doubly stochastic Poisson processes*, vol. 529. Springer, 2006.

[33] GRONWALL, T. The gamma function in the integral calculus. *Annals of mathematics* 20, 2 (1918), 35–124.

[34] HALLIDAY, D., ROSENBERG, J., AMJAD, A., BREEZE, P., CONWAY, B., AND FARMER, S. A framework for the analysis of mixed time series/point process data—theory and application to the study of physiological tremor, single motor unit discharges and electromyograms. *Progress in biophysics and molecular biology* 64, 2-3 (1995), 237–278.

[35] HANSEN, J.-P., AND MCDONALD, I. R. *Theory of simple liquids: with applications to soft matter*. Academic press, 2013.

[36] HAWAT, D., GAUTIER, G., BARDENET, R., AND LACHIÈZE-REY, R. On estimating the structure factor of a point process, with applications to hyperuniformity. *Statistics and Computing* 61, 33 (2023).

[37] HEINRICH, L. Asymptotic gaussianity of some estimators for reduced factorial moment measures and product densities of stationary poisson cluster processes. *Statistics* 19, 1 (1988), 77–86.
[38] HEINRICH, L. On the strong Brillinger-mixing property of α-determinantal point processes and some applications. *Applications of Mathematics* 61 (2016), 443–461.

[39] HEINRICH, L., AND PROKEŠOVÁ, M. On estimating the asymptotic variance of stationary point processes. *Methodology and Computing in Applied Probability* 12, 3 (2010), 451–471.

[40] HEINRICH, L., AND SCHMIDT, V. Normal convergence of multidimensional shot noise and rates of this convergence. *Advances in Applied Probability* 17, 4 (1985), 709–730.

[41] HURVICH, C. M. Data-driven choice of a spectrum estimate: extending the applicability of cross-validation methods. *Journal of the American Statistical Association* 80, 392 (1985), 933–940.

[42] JARRAH, A., NIANGA, J.-M., IOST, A., GUILLEMOT, G., AND NAJJAR, D. On the detection of corrosion pit interactions using two-dimensional spectral analysis. *Corrosion science* 52, 2 (2010), 303–313.

[43] JOVANOVIĆ, S., HERTZ, J., AND ROTTER, S. Cumulants of hawkes point processes. *Physical Review E* 91, 4 (2015), 042802.

[44] KARNIK, S., ROMBERG, J., AND DAVENPORT, M. A. Thomson's multitaper method revisited. *IEEE Transactions on Information Theory* 68, 7 (2022), 4864–4891.

[45] KARTSIOUKAS, R., STOEV, S., AND HSING, T. Spectral density estimation of function-valued spatial processes. *arXiv preprint arXiv:2302.02247* (2023).

[46] KHABOU, M., AND TORRISI, G. L. Gaussian approximation and moderate deviations of poisson shot noises with application to compound generalized hawkes processes. *Advances in Applied Probability* (2024), 1–41.

[47] KIM, J., AND TORQUATO, S. Effect of imperfections on the hyperuniformity of many-body systems. *Physical Review B* 97, 5 (2018), 054105.

[48] KITTEL, C., AND MCEUEN, P. *Introduction to solid state physics*. John Wiley & Sons, 2018.

[49] KUTOYANTS, Y. A. *Statistical inference for spatial Poisson processes*, vol. 134. Springer Science & Business Media, 2012.

[50] LANDAU, H. J., AND POLLAK, H. O. Prolate spheroidal wave functions, fourier analysis and uncertainty—ii. *Bell System Technical Journal* 40, 1 (1961), 65–84.

[51] LAVANCIER, F., MÖLLER, J., AND RUBAK, E. Determinantal point process models and statistical inference. *Journal of the Royal Statistical Society Series B: Statistical Methodology* 77, 4 (2015), 853–877.

[52] LIN, X., AND KANG, H. Efficient, cross-fitting estimation of semiparametric spatial point processes. *arXiv preprint arXiv:2410.04359* (2024).

[53] MASTRILLI, G., BŁASZCZYSZYN, B., AND LAVANCIER, F. Estimating the hyperuniformity exponent of point processes. *arXiv preprint arXiv:2407.16797* (2024).

[54] McCULLAGH, P. *Tensor methods in statistics: Monographs on statistics and applied probability*. Chapman and Hall/CRC, 2018.
[55] MOLLER, J., AND WAAGEPETERSEN, R. P. *Statistical inference and simulation for spatial point processes*. CRC press, 2003.

[56] MØLLER, J., AND WAAGEPETERSEN, R. P. Modern statistics for spatial point processes. *Scandinavian Journal of Statistics* 34, 4 (2007), 643–684.

[57] MUGGLESTONE, M. A., AND RENSHAW, E. The exploratory analysis of bivariate spatial point patterns using cross-spectra. *Environmetrics* 7, 4 (1996), 361–377.

[58] MUGGLESTONE, M. A., AND RENSHAW, E. A practical guide to the spectral analysis of spatial point processes. *Computational Statistics & Data Analysis* 21, 1 (1996), 43–65.

[59] MUGGLESTONE, M. A., AND RENSHAW, E. Spectral tests of randomness for spatial point patterns. *Environmental and Ecological Statistics* 8, 3 (2001), 237–251.

[60] PARKS, T. W., AND SHENOY, R. G. Time-frequency concentrated basis functions. In *International Conference on Acoustics, Speech, and Signal Processing* (1990), IEEE, pp. 2459–2462.

[61] PECCATI, G., AND TAQQU, M. S. *Wiener Chaos: Moments, Cumulants and Diagrams: A survey with computer implementation*, vol. 1. Springer Science & Business Media, 2011.

[62] PERCIVAL, D. B., AND WALDEN, A. T. *Spectral analysis for univariate time series*, vol. 51. Cambridge University Press, 2020.

[63] PROKEŠOVÁ, M., AND JENSEN, E. B. V. Asymptotic palm likelihood theory for stationary point processes. *Annals of the Institute of Statistical Mathematics* 65, 2 (2013), 387–412.

[64] RAJALA, T. A., OLHEDE, S. C., GRAINGER, J. P., AND MURRELL, D. J. What is the fourier transform of a spatial point process? *IEEE Transactions on Information Theory* (2023).

[65] RENSHAW, E. Spectral techniques in spatial analysis. *Forest Ecology and Management* 94, 1-3 (1997), 165–174.

[66] RENSHAW, E. Two-dimensional spectral analysis for marked point processes. *Biometrical Journal: Journal of Mathematical Methods in Biosciences* 44, 6 (2002), 718–745.

[67] REYNAUD-BOURET, P. Adaptive estimation of the intensity of inhomogeneous poisson processes via concentration inequalities. *Probability Theory and Related Fields* 126, 1 (2003), 103–153.

[68] REYNAUD-BOURET, P., AND RIVOIRARD, V. Near optimal thresholding estimation of a poisson intensity on the real line. *Electronic Journal of Statistics* 4 (2010), 172–238.

[69] RIEDEL, K. S., AND SIDORENKO, A. Minimum bias multiple taper spectral estimation. *IEEE Transactions on Signal Processing* 43, 1 (1995), 188–195.

[70] SAMAROV, A. Lower bound on the risk for spectral density estimates. *Problemy Peredachi Informatsii* 13, 1 (1977), 67–72.

[71] SAULIS, L., AND STATULEVICIUS, V. *Limit theorems for large deviations*, vol. 73. Springer Science & Business Media, 2012.

[72] SIMONS, F. J., AND WANG, D. V. Spatiospectral concentration in the cartesian plane. *GEM-International Journal on Geomathematics* 2, 1 (2011), 1–36.
[73] SLEPIAN, D., AND POLLAK, H. O. Prolate spheroidal wave functions, fourier analysis and uncertainty—i. *Bell System Technical Journal* 40, 1 (1961), 43–63.

[74] SOSHNIKOV, A. Determinantal random point fields. *Russian Mathematical Surveys* 55, 5 (2000), 923.

[75] THOMSON, D. J. Spectrum estimation and harmonic analysis. *Proceedings of the IEEE* 70, 9 (1982), 1055.

[76] TORQUATO, S. Hyperuniform states of matter. *Physics Reports* 745 (2018), 1–95.

[77] TSYBAKOV, A. B. *Introduction à l'estimation non paramétrique*, vol. 41. Springer Science & Business Media, 2003.

[78] TUCKEY, J. The sampling theory of power spectrum estimates. In *Symposium on Applications of Autocorrelation Analysis to Physical Problems. US Office of Naval Research* (1950), pp. 47–67.

[79] VELASCO, C. Local cross-validation for spectrum bandwidth choice. *Journal of Time Series Analysis* 21, 3 (2000), 329–361.

[80] VERSHYNIN, R. *High-dimensional probability: An introduction with applications in data science*, vol. 47. Cambridge university press, 2018.

[81] WALDEN, A. T. A unified view of multitaper multivariate spectral estimation. *Biometrika* 87, 4 (2000), 767–788.

[82] XU, G., ZHAO, C., JALILIAN, A., WAAGEPETERSEN, R., ZHANG, J., AND GUAN, Y. Non-parametric estimation of the pair correlation function of replicated inhomogeneous point processes. *Electronic Journal of Statistics* 14, 2 (2020).

[83] YANG, J., AND GUAN, Y. Fourier analysis of spatial point processes. *arXiv preprint arXiv:2401.06403* (2024).

[84] ZHU, L., YANG, J., JUN, M., AND COOK, S. On minimum contrast method for multivariate spatial point processes. *Electronic Journal of Statistics* 19, 1 (2025), 1889–1941.

# A Proofs related to Section 5.2

## A.1 Variance study

The next lemma provides an upper bound on the variance of the multitaper estimators. Its main argument relies on Bessel’s inequality:

$$ \sum_{i \in I} |\langle f_i, h \rangle|^2 \le \|h\|_2^2, $$

for any orthonormal family of function $(f_i)_{i \in I}$ of $L^2(\mathbb{R}^d)$ and $h \in L^2(\mathbb{R}^d)$. This allows us to exploit the orthogonality of the tapers in a non-asymptotic regime. This argument will be used extensively throughout the proofs.
**Lemma A.1.** Let $\Phi$ satisfying Assumption 3.1 with intensity $\lambda = 1$. Let $W$ be a subset of $\mathbb{R}^d$, $I$ be a discrete subset of $\mathbb{N}^d$ and let $k_0 \in \mathbb{R}^d$. We consider a family of orthonormal functions $(f_i)_{i \in I} \in L^2(\mathbb{R}^d)^{|I|}$. Then,

$$ \mathrm{Var}[\hat{S}(k_0)] \le \frac{2\|S\|_{\infty}^2}{|I|} \left( 1 + \left( \sum_{i \in I} \|f_i\|_{1^{\mathbb{R}^d} \setminus W}^2 \right)^{\frac{1}{2}} \right)^2 + \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_4^2 \right)^2 |\gamma_{\mathrm{red}}^4|, \quad (34) $$

where $|\gamma_{\mathrm{red}}^4| := \lambda + 7|\gamma_{\mathrm{red}}^{(2)}|(\mathbb{R}^d) + 6|\gamma_{\mathrm{red}}^{(3)}|(\mathbb{R}^{2d}) + |\gamma_{\mathrm{red}}^{(4)}|(\mathbb{R}^{3d})$.

*Proof.* Recall that $\hat{S}(k_0) = |I|^{-1} \sum_{i \in I} |C_i(k_0)|^2$ where $(C_i(k_0))_{i \in I}$ is defined in (9). Using the bilinearity of the covariance, we get

$$ \mathrm{Var}[\hat{S}(k_0)] = \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} \kappa_2(|C_{i_1}(k_0)|^2, |C_{i_2}(k_0)|^2). $$

According to Corollary G.1, we have

$$ \kappa_2(|C_{i_1}(k_0)|^2, |C_{i_2}(k_0)|^2) = \sum_{\substack{\sigma \in \Pi[4], \\ \sigma \lor \tau = \mathbf{1}_4}} \prod_{b \in \sigma} \kappa(D_{i_j}(k_0); j \in b), $$

where $\tau = 12|34, \lor$ denotes the least-upper-bound between partitions (see Section G) and $D_{i_1}(k_0) = C_{i_1}(k_0), D_{i_2}(k_0) = \overline{C_{i_1}(k_0)}, D_{i_3}(k_0) = C_{i_2}(k_0), D_{i_4}(k_0) = \overline{C_{i_2}(k_0)}$. Note that we are considering cumulants of centered random variables, so the sum runs over the partitions with blocks of size larger than 2. Hence:

$$ \begin{aligned} \kappa_2(|C_{i_1}(k_0)|^2, |C_{i_2}(k_0)|^2) &= \sum_{\substack{\sigma \in \Pi_2[4], \\ \sigma \lor \tau = \mathbf{1}_4}} \prod_{b \in \sigma} \kappa(E_{i_j}(k_0); j \in b) \\ &= \kappa(T_{i_1}(k_0), \overline{T_{i_1}(k_0)}, T_{i_2}(k_0), \overline{T_{i_2}(k_0)}) + \sum_{\substack{\sigma \in \Pi_2[4], \\ \sigma \lor \tau = \mathbf{1}_4, |\sigma| = 2}} \prod_{b \in \sigma} \mathrm{Cov}[E_{i_{b_1}}(k_0), E_{i_{b_2}}(k_0)], \end{aligned} \quad (35) $$

where $\Pi_2[4]$ denotes the sets of partitions of [4] having only blocks of size larger than 2 and $E_{i_1}(k_0) = T_{i_1}(k_0), E_{i_2}(k_0) = \overline{T_{i_1}(k_0)}, E_{i_3}(k_0) = T_{i_2}(k_0), E_{i_4}(k_0) = \overline{T_{i_2}(k_0)}$. We bound the cumulant corresponding to the block of order 4 using Lemma D.2:

$$ |\kappa(T_{i_1}(k_0), \overline{T_{i_1}(k_0)}, T_{i_2}(k_0), \overline{T_{i_2}(k_0)})| \le \|f_{i_1}\|_4^2 \|f_{i_2}\|_4^2 \sum_{\pi \in \Pi[4]} |\gamma_{\mathrm{red}}^{(|\pi|)}(\mathbb{R}^{d(|\pi|-1)})|. \quad (36) $$

Subsequently,

$$ \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} |\kappa(T_{i_1}(k_0), \overline{T_{i_1}(k_0)}, T_{i_2}(k_0), \overline{T_{i_2}(k_0)})| \le \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_4^2 \right)^2 |\gamma_{\mathrm{red}}^4|. \quad (37) $$

We now focus on the covariance terms. Note that $\tau = 12|34$, so the only two partitions $\sigma \in \Pi_2[4]$ such that $\sigma \lor \tau = \mathbf{1}_4$ are $13|24$ and $14|23$. Accordingly,

$$ \begin{aligned} \sum_{\substack{\sigma \in \Pi[4], \\ \sigma \lor \tau = \mathbf{1}_4, |\sigma| = 2}} \prod_{b \in \sigma} \mathrm{Cov}[E_{i_{b_1}}(k_0), E_{i_{b_2}}(k_0)] &= \mathrm{Cov}[E_{i_1}(k_0), E_{i_3}(k_0)] \mathrm{Cov}[E_{i_2}(k_0), E_{i_4}(k_0)] \\ &\quad + \mathrm{Cov}[E_{i_1}(k_0), E_{i_4}(k_0)] \mathrm{Cov}[E_{i_2}(k_0), E_{i_3}(k_0)] \\ &= |\mathrm{Cov}[T_{i_1}(k_0), T_{i_2}(k_0)]|^2 + |\mathrm{Cov}[T_{i_1}(k_0), T_{i_2}(-k_0)]|^2. \end{aligned} $$
According to Lemma D.3, we have

$$ \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} \left( |\mathrm{Cov}[T_{i_1}(k_0), T_{i_2}(k_0)]|^2 + |\mathrm{Cov}[T_{i_1}(k_0), T_{i_2}(-k_0)]|^2 \right) \\ \leq \frac{2 \|S\|_{\infty}^2}{|I|} \left( 1 + \left( \sum_{i \in I} \|f_i \mathbf{1}_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{1/2} \right)^2, \quad (38) $$

which concludes the proof of this lemma.

<figure><img src="image_8.png" /></figure>

## A.2 Bias study

Next lemma studies the bias for the multitaper estimator. Note that the factor $\sqrt{d}$ appears only when $\beta > 1$ and is due to the assumption that $S \in \Theta(\beta, L)$, which concerns the coefficient of the gradient of $S$ and not their $\ell^2([d])$ norm.

**Lemma A.2.** Let $\Phi$ satisfying Assumption 3.1 with intensity $\lambda = 1$. Let $W$ be a subset of $\mathbb{R}^d$, $I$ be a discrete subset of $\mathbb{N}^d$, $k_0 \in \mathbb{R}^d$ and let $(f_i)_{i \in I}$ be a family of functions with unit $L^2(\mathbb{R}^d)$ norm. Assume that $S \in \Theta(\beta, L)$ with $(\beta, L) \in (0, \infty)^2$ (see Definition 1). Assume that $\beta \le 2$. Then,

$$ \left| \mathbb{E} \left[ \hat{S}(k_0) \right] - S(k_0) \right| \le \frac{L \sqrt{d}}{|I|} \sum_{i \in I} \|f_i\|_{\hat{H}^{\beta/2}}^2 + \frac{4 \|S\|_{\infty}}{|I|} \sum_{i \in I} \|f_i \mathbf{1}_{\mathbb{R}^d \setminus W}\|_2. \quad (39) $$

*Proof.* Assume first that $\beta \le 1$. Using the fact that $\|f_i\|_2 = 1$ for all $i \in I$, we get

$$ \left| \mathbb{E} \left[ \hat{S}(k_0) \right] - S(k_0) \right| \le A + B, $$

where

$$ A := \frac{1}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} \left| \mathcal{F}[f_i \mathbf{1}_W](k) \right|^2 - \left| \mathcal{F}[f_i](k) \right|^2 \left( S(k_0 + k) + S(k_0) \right) dk, \quad (40) $$

$$ B = \frac{1}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} |\mathcal{F}[f_i](k)|^2 |S(k_0 + k) - S(k_0)| dk. $$

Since $S$ is $\beta$-Hölder with $\beta \le 1$, we have

$$ B \le \frac{1}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} |\mathcal{F}[f_i](k)|^2 L |k|^\beta dk = \frac{L}{|I|} \sum_{i \in I} \|f_i\|_{\hat{H}^{\beta/2}}^2. $$

Using the Cauchy-Schwarz inequality and the Plancherel theorem we obtain

$$ \begin{aligned} A &\le \frac{2 \|S\|_{\infty}}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} |\mathcal{F}[f_i \mathbf{1}_{\mathbb{R}^d \setminus W}](k)|^2 \|\mathcal{F}[f_i \mathbf{1}_W](k)\| + |\mathcal{F}[f_i](k)| \|f_i\|_{\hat{H}^{\beta/2}}^2 dk \\ &\le \frac{2 \|S\|_{\infty}}{|I|} \sum_{i \in I} \|\mathcal{F}[f_i \mathbf{1}_{\mathbb{R}^d \setminus W}]\|_2 \|\mathcal{F}[f_i \mathbf{1}_W]\|_2 + |\mathcal{F}[f_i]| \|f_i\|_{\hat{H}^{\beta/2}}^2. \end{aligned} $$
which concludes the proof when $\beta \le 1$. Assume now that $\beta \in (1, 2]$. Hence $S$ is differentiable and using the fact that $k \mapsto |\mathcal{F}[f_i 1_W](k)|^2 \nabla S(k_0) \cdot k$ is even, we have

$$ \left| \mathbb{E} \left[ \hat{S}(k_0) \right] - S(k_0) \right| = \left| \frac{1}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} |\mathcal{F}[f_i 1_W](k)|^2 (S(k_0 + k) - S(k_0) - \nabla S(k_0) \cdot k) dk \right| \le A + B', $$

where $A$ is defined in (40) and

$$ B' := \frac{1}{|I|} \sum_{i \in I} \int_{\mathbb{R}^d} |\mathcal{F}[f_i](k)|^2 |S(k_0 + k) - S(k_0) - \nabla S(k_0) \cdot k| dk. $$

Since $S$ is $\beta$-Hölder with $\beta \ge 1$, we have

$$ \begin{aligned} |S(k_0 + k) - S(k_0) - \nabla S(k_0) \cdot k| &= \left| \int_0^1 (\nabla S(k_0 + tk) - \nabla S(k_0)) \cdot k dt \right| \\ &\le \int_0^1 \sqrt{dL} t^{\beta-1} |k|^{\beta-1} |k| dt \le L \sqrt{d} |k|^{\beta}. \end{aligned} \quad (41) $$

Then, concluding as in the proof of the case $\beta \le 1$, we obtain the result when $\beta \in (1, 2]$. $\square$

## A.3 Proofs of Theorem 5.1 and Corollary 5.1

In this section, gathering Lemmas A.1 and A.2, we prove Theorem 5.1. Then, we specify Theorem 5.1 for the Hermite tapers with Corollary 5.1.

*Proof of Theorem 5.1.* The bias variance decomposition writes

$$ s := \mathbb{E} \left[ \left| \hat{S}(k_0) - S(k_0) \right|^2 \right] = \mathrm{Var}[\hat{S}(k_0)] + \left| \mathbb{E} \left[ \hat{S}(k_0) \right] - S(k_0) \right|^2. $$

According to Lemma A.1 and Lemma A.2, we have, by triangular inequality

$$ \begin{aligned} s &\le \left( \frac{2 \|S\|_{\infty}^2}{|I|} \left( 1 + \left( \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}} \right)^2 + \left( \frac{1}{|I|} \sum_{i \in I} \|f_i\|_4^2 \right)^2 |\gamma_{\mathrm{red}}^4|(\mathbb{R}^{3d}) \right. \\ & \qquad \left. + \left( \frac{\sqrt{dL}}{|I|^{\beta}} \sum_{i \in I} \|f_i\|_{\hat{H}^{\beta/2}}^2 + \frac{4 \|S\|_{\infty}}{|I|} \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^2 \right)^{\frac{1}{2}} \\ &\le \frac{\sqrt{2} \|S\|_{\infty}}{|I|^{\frac{1}{2}}} + \frac{1}{|I|} \sum_{i \in I} \|f_i\|_4^2 |\gamma_{\mathrm{red}}^4|(\mathbb{R}^{3d})^{\frac{1}{2}} + \frac{\sqrt{dL}}{|I|} \sum_{i \in I} \|f_i\|_{\hat{H}^{\beta/2}}^2 \\ & \qquad + \sqrt{2} \|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}} + \frac{4 \|S\|_{\infty}}{|I|} \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2. \end{aligned} $$

We conclude using $\sqrt{2} + 4 \le 6$ and the Cauchy-Schwarz inequality:

$$ \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2 \le |I|^{\frac{1}{2}} \left( \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}}. $$

$\square$
**Proof of Corollary 5.1.** The change of variable $x \leftrightarrow x/r$ and Lemma F.1 applied with $\rho = R/r$ give $W\text{-loc} \le Ce^{-c'(R/r)^{3\theta}} \le Ce^{-c'|I|^{\frac{3\theta}{2d}}}$, where $0 < c', C < \infty$ are constants. In the second inequality we used (13). With (14), we obtain $W\text{-loc} \le Ce^{-c|W|^{\frac{3\beta\theta}{d(2\beta+d)}}}$, for constants that do not depend on $r, I$ and $W$. Moreover, Lemma F.2 and Lemma F.3 yield

$$ \mathcal{F}\text{-loc} \le \frac{L\sqrt{d}|I|^{\beta/(2d)}}{r^{\beta}}, \quad \frac{1}{|I|} \sum_{i \in I} \|\psi_i\|_4^2 \le c_{\psi}^d \frac{\log(|I|^{\frac{1}{d}})^{\frac{d}{2}}}{r^{d/2}|I|^{1/4}}, $$

where $c_{\psi} \in (0, \infty)$ is a numerical constant. Condition (13) implies that $r \ge R/(\sqrt{2}|I|^{1/(2d)})$. Subsequently, we get, with Theorem 5.1

$$ \mathbb{E} \left[ \left| \hat{S}(k_0) - S(k_0) \right|^2 \right]^{\frac{1}{2}} \le \frac{\sqrt{2} \|S\|_{\infty}}{|I|^{\frac{1}{2}}} + \frac{L \sqrt{d}}{2^{\frac{\beta}{2}}} \left( \frac{|I|}{|W|} \right)^{\frac{\beta}{d}} + \left( \frac{c_{\psi}}{\sqrt{2}} \right)^d |\gamma_{\text{red}}^4|^{\frac{1}{2}} \frac{\log(|I|^{\frac{1}{d}})^{\frac{d}{2}}}{|W|^{\frac{1}{2}}} + 6C \|S\|_{\infty} e^{-c|W|^{\frac{3\beta\theta}{d(2\beta+d)}}} \quad (42) $$

Using condition (14) we conclude the proof with elementary computations. $\square$

## B Proof related to Section 5.3

The results of Section 5.3 rely on the control the cumulants of the multitaper estimator. Lemma B.1, which constitutes a first step for a workable bound on the cumulants of $\hat{S}(k_0)$, is the core of the proof of Theorem 5.2. It relies on the Isserlis-Wick formula of Corollary G.1 that allows to express the joint cumulants of square of random variables as a sum of products of joint cumulants of these random variables. The key point of this lemma is to provide a bound on the cumulants of the multitaper estimator where the tapers indexes are summed (see $\Sigma$ defined in (18)). To this end, we make essential use of the combinatorial Lemma B.2, which crucially exploits the connectivity condition of the partitions involved.

**Lemma B.1.** Let $\Phi$ satisfying Assumption 3.1 with intensity $\lambda = 1$ and such that $S \in \Theta(\beta, L)$ with $\beta \in (0, 2]$ and $L < \infty$. Let $W$ be a subset of $\mathbb{R}^d$ and $I$ be a discrete subset of $\mathbb{N}^d$. Assume that $(f_i)_{i \in I}$ is family of orthonormal functions that are $L^1(\mathbb{R}^d) \cap L^\infty(\mathbb{R}^d)$. If $\beta > 1$, assume for all $(i, i') \in I^2$, $\mathcal{F}[f_i](k)\mathcal{F}[f_{i'}](k) = \mathcal{F}[f_i](-k)\mathcal{F}[f_{i'}](-k)$ for all $k \in \mathbb{R}^d$. Let $k_0 \in \mathbb{R}^d$. Then, for all $m \ge 2$,

$$ |\kappa_m(\hat{S}(k_0))| \le 2^{m-1}(m-1)! \frac{S(k_0)^m}{|I|^{m-1}} + \sum_{\substack{\sigma \in \Pi[2m], \\ \sigma \lor \tau = 1_{2m}, |\sigma| \ge 2}} \prod_{b \in \sigma} K_{|b|}, \quad (43) $$

where, for $\Sigma$ defined in (18):

$$ K_1 = 0, K_2 = \Sigma, \quad \forall s \ge 3, K_s = \left( \frac{1}{|I|} \sum_{i \in I} a_s(i)^2 \right)^{\frac{s}{2}} \sum_{\pi \in \Pi[s]} |\gamma_{\text{red}}^{(|\pi|)}|(\mathbb{R}^d(|\pi|-1)). \quad (44) $$

Finally, $a_s(i) := \|f_i\|_s$ for $s \ge 4$ and

$$ a_3(i) := \left( 3 \sup_{i \in I} \|f_i\|_{\mathbb{R}^d \setminus W} \sup_{i \in I} \|f_i\|_3^2 + 3 \sup_{i \in I} \|\mathcal{F}[f_i]\|_{1-|\cdot| \ge |k_0|/3} \right)^{\frac{1}{3}} \wedge \|f_i\|_3. $$
Proof. Using the m-linearity of the cumulants we get

$$ \kappa_m(\hat{S}(k_0)) = \frac{1}{|I|^m} \sum_{(i_1, \dots, i_m) \in I^m} \kappa(|C_{i_1}(k_0)|^2, \dots, |C_{i_m}(k_0)|^2). $$

Recall that $(C_i(k_0))_{i \in I}$ are defined in (9). Then, using Corollary G.1 and arguing as in the proof of Lemma A.1, we obtain

$$ \kappa(|C_{i_1}(k_0)|^2, \dots, |C_{i_m}(k_0)|^2) = \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \vee \tau = 1_{2m}}} \prod_{b \in \sigma} \kappa(E_{i_j}(k_0); j \in b), \quad (45) $$

where $\Pi_2[2m]$ denotes the sets of partitions of $[2m]$ having only blocks of size larger than 2, $E_{i_{2p-1}}(k_0) = T_{i_p}(k_0)$ and $E_{i_{2p}}(k_0) = T_{i_p}(k_0)$ for $p \in [m]$, and $\tau = 12|\dots|(2m-1)2m$. We bound the cumulants corresponding to the blocks of size larger than 3 using Lemma D.2:

$$ |\kappa(E_{i_j}(k_0); j \in b)| \le \prod_{j \in b} \|f_{i_{\lfloor j/2 \rfloor}}\|_{|b|} \sum_{\pi \in \Pi[|b|]} |\gamma_{\text{red}}^{(\lfloor \pi \rfloor)}| (\mathbb{R}^{d(|\pi|-1)}). $$

Moreover, when $|b| = 3$, Lemma D.5 yields

$$ \begin{aligned} |\kappa(E_{i_j}(k_0); j \in b)| & \le 3 \sum_{\pi \in \Pi[3]} |\gamma_{\text{red}}^{(\lfloor \pi \rfloor)}| (\mathbb{R}^{d(|\pi|-1)}) \left( \sup_{i \in I} \|f_i\mathbf{1}_{\mathbb{R}^d \setminus W}\|_3 \sup_{i \in I} \|f_i\|_3^2 + \sup_{i \in I} \|\mathcal{F}[f_i]\mathbf{1}_{|\cdot| \ge |k_0|/3}\|_1 \right). \end{aligned} $$

Indeed, the sum of frequencies involved in the functions $f_i$ is of the form $k_0 \sum_{s=1}^3 \epsilon_s$, with $(\epsilon_s)_{s \in [3]} \in \{-1, 1\}^3$. Hence $|k_0 \sum_{s=1}^3 \epsilon_s| \ge |k_0|$. Accordingly,

$$ \begin{aligned} |\kappa_m(|C_{i_1}(k_0, r)|^2, \dots, |C_{i_m}(k_0)|^2)| & \le \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \vee \tau = 1_{2m}}} \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \prod_{\substack{j \in b \\ |b| \ge 3}} a_{|b|(i_{\lfloor j/2 \rfloor})} \sum_{\pi \in \Pi[|b|]} |\gamma_{\text{red}}^{(\lfloor \pi \rfloor)}| (\mathbb{R}^{d(|\pi|-1)}) \prod_{\substack{b \in \sigma \\ |b|=2}} |\kappa(E_{i_j}(k_0); j \in b)|. \end{aligned} $$

We have thus obtained:

$$ |\kappa_m(\hat{S}(k_0))| \le \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \vee \tau = 1_{2m}}} \left( \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \sum_{\pi \in \Pi[|b|]} |\gamma_{\text{red}}^{(\lfloor \pi \rfloor)}| (\mathbb{R}^{d(|\pi|-1)}) \right) A_{\sigma}, $$

with

$$ A_{\sigma} := \frac{1}{|I|^m} \sum_{(i_1, \dots, i_m) \in I^m} \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \prod_{\substack{j \in b \\ |b| \ge 3}} a_{|b|(i_{\lfloor j/2 \rfloor})} \prod_{\substack{b \in \sigma \\ |b|=2}} |\kappa(E_{i_j}(k_0); j \in b)|. $$

To continue the proof, we introduce some notations. Denote by $\mathcal{B}_2(\sigma)$ the blocks of size 2 of $\sigma$. Let $\epsilon_p = 1$ if $p$ is even and $-1$ otherwise. Moreover, denote

$$ \mathcal{B}_2^+(\sigma) := \{b = (b_1, b_2) \in \mathcal{B}_2(\sigma) | \epsilon_{b_1} \epsilon_{b_2} = +1\}, $$
and $\mathcal{B}_2^-(\sigma) := \mathcal{B}_2(\sigma) \setminus \mathcal{B}_2^+(\sigma)$. Finally, for $b = (b_1, b_2) \in \mathcal{B}_2(\sigma)$, let $\delta_b := \mathbb{1}_{i_{[b_1/2]}=i_{[b_2/2]}}$. We introduce these sets since we anticipate, for $b = (b_1, b_2) \in \mathcal{B}_2^+(\sigma)$,

$$ \begin{aligned} \kappa(E_{i_j}(k_0); j \in b) &= \int_{\mathbb{R}^d} \mathcal{F}[f_{i_{[\frac{b_1}{2}]}} \mathbb{1}_W](k + \varepsilon_{b_1} k_0) \overline{\mathcal{F}[f_{i_{[\frac{b_2}{2}]}} \mathbb{1}_W]}(k + \varepsilon_{b_2} k_0) S(k) dk \\ &\approx S(k_0) \langle f_{i_{[\frac{b_1}{2}]}}, f_{i_{[\frac{b_2}{2}]}} \rangle = S(k_0) \delta_b, \end{aligned} $$

whereas for $b = (b_1, b_2) \in \mathcal{B}_2^-(\sigma)$, with $b_1$ even and $b_2$ odd,

$$ \kappa(E_{i_j}(k_0); j \in b) \approx S(k_0) \langle \mathcal{F}[f_{i_{[\frac{b_1}{2}]}}](\cdot + k_0), \mathcal{F}[f_{i_{[\frac{b_2}{2}]}}](\cdot - k_0) \rangle. $$

To take into account this previous heuristic, we denote, for $b \in \mathcal{B}_2^+(\sigma)$,

$$ \Delta_b := \kappa(E_{i_j}(k_0); j \in b) - S(k_0)\delta_b. $$

In the term $A_\sigma$, we expand the product over the blocks of size two of $\sigma$ as follows:

$$ \begin{aligned} \prod_{\substack{b \in \sigma \\ |b|=2}} |\kappa(E_{i_j}(k_0); j \in b)| &= \prod_{b \in \mathcal{B}_2^-(\sigma)} |\kappa(E_{i_j}(k_0); j \in b)| \prod_{b \in \mathcal{B}_2^+(\sigma)} (S(k_0)\delta_b + \Delta_b) \\ &= \sum_{P \in \mathcal{P}(\mathcal{B}_2^+(\sigma))} \prod_{b \in \mathcal{B}_2^-(\sigma)} |\kappa(E_{i_j}(k_0); j \in b)| \prod_{b \in P} S(k)\delta_b \prod_{b \in P^c} \Delta_b, \end{aligned} $$

where $\mathcal{P}(\mathcal{B}_2^+(\sigma))$ denote the subsets of $\mathcal{B}_2^+(\sigma)$. This yields:

$$ A_{\sigma} = \frac{1}{|I|^m} \sum_{P \in \mathcal{P}(\mathcal{B}_2^+(\sigma))} S(k_0)^{|P|} B_{\sigma,P}, $$

where

$$ B_{\sigma,P} := \sum_{(i_1, \dots, i_m) \in I^m} \prod_{\substack{b \in \sigma \\ |b| \ge 3}} a_{|b|}(i_{[j/2]}) \prod_{b \in \mathcal{B}_2^-(\sigma)} |\kappa(E_{i_j}(k_0); j \in b)| \prod_{b \in P} \delta_b \prod_{b \in P^c} \Delta_b, $$

where we recall that $\delta_b$ and $\Delta_b$ depend on indexes $i_1, \dots, i_m$. We distinguish several cases. Assume that $\sigma$ contains only blocks of size 2 and that $|P| = m$. In that case, according to Lemma B.2, $B_{\sigma,P} = |I|$. Otherwise, Lemma B.2 yields:

$$ B_{\sigma,P} \le \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \left( \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}} \left( \sum_{(i,i') \in I^2} |a_2^-(i,i')|^2 \right)^{\frac{|\mathcal{B}_2^-(\sigma)|}{2}} \left( \sum_{(i,i') \in I^2} |a_2^+(i,i')|^2 \right)^{\frac{|P^c|}{2}}, $$

where

$$ a_2^-(i, i') := \int_{\mathbb{R}^d} \mathcal{F}[f_i \mathbb{1}_W](k + k_0) \overline{\mathcal{F}[f_{i'} \mathbb{1}_W]}(k - k_0) S(k) dk, $$

$$ a_2^+(i, i') := \int_{\mathbb{R}^d} \mathcal{F}[f_i \mathbb{1}_W](k + k_0) \overline{\mathcal{F}[f_{i'} \mathbb{1}_W]}(k + k_0) S(k) dk - S(k_0) \mathbb{1}_{i=i'}. $$
Note that, using $S(\cdot) = S(-\cdot)$ and changes of variable, we have $|a_2^\pm(i, i')| = |a_2^\pm(i', i)|$. If $\sigma$ contains at least one block of size larger than three, we have, using the binomial theorem:

$$
\begin{aligned}
A_{\sigma} &\le \frac{1}{|I|^m} \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \left( \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}} \left( \sum_{(i,i') \in I^2} |a_2^-(i, i')|^2 \right)^{\frac{|B_2^-(\sigma)|}{2}} \left( S(k_0) + \left( \sum_{(i,i') \in I^2} |a_2^+(i, i')|^2 \right)^{\frac{1}{2}} \right)^{|B_2^+(\sigma)|} \\
&\le \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \left( \frac{1}{|I|} \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}} \prod_{b \in B_2(\sigma)} \Sigma_0,
\end{aligned}
$$

where

$$
\Sigma_0 := \frac{S(k_0)}{|I|} + \frac{1}{|I|} \left( \sum_{(i,i') \in I^2} |a_2^+(i,i')|^2 \right)^{\frac{1}{2}} + \frac{1}{|I|} \left( \sum_{(i,i') \in I^2} |a_2^-(i,i')|^2 \right)^{\frac{1}{2}}.
$$

Indeed, for the term $|I|^m$, we use the fact that $2m = 2|\mathcal{B}_2(\sigma)| + \sum_{b \in \sigma, |b| \ge 3} |b|$. If $\sigma$ contains only blocks of size 2, we obtain:

$$
\begin{aligned}
A_{\sigma} &\le \frac{S(k_0)^m}{|I|^{m-1}} + \frac{1}{|I|^m} \left( \sum_{(i,i') \in I^2} |a_2^-(i,i')|^2 \right)^{\frac{|B_2^-(\sigma)|}{2}} \sum_{P \in \mathcal{P}(B_2^+(\sigma)), |P| < m} S(k_0)^{|P|} \left( \sum_{(i,i') \in I^2} |a_2^+(i,i')|^2 \right)^{\frac{|P^c|}{2}} \\
&\le \frac{S(k_0)^m}{|I|^{m-1}} + \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \left( \frac{1}{|I|} \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}} \prod_{b \in B_2(\sigma)} \Sigma_0.
\end{aligned}
$$

Gathering both cases and using the fact that $|\{\sigma \in \Pi[2m] | |\sigma| = m, \sigma \vee \tau = \mathbf{1}_{2m}\}| = 2^{m-1}(m-1)!$ (exercise 3.5 page 92 of [54]), we get

$$
|\kappa_m(\hat{S}(k_0))| \le 2^{m-1}(m-1)! \frac{S(k_0)^m}{|I|^{m-1}} \\
\qquad + \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \vee \tau = \mathbf{1}_{2m}}} \left( \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \sum_{\pi \in \Pi([b])} |\gamma_{\text{red}}^{(|\pi|)}(\mathbb{R}^d(\pi|-1))| \right) \left( \frac{1}{|I|} \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}} \prod_{b \in B_2(\sigma)} \Sigma_0,
$$

where in the second term we use the fact that $1_{|\sigma|=m} + 1_{|\sigma|<m} = 1$. The end of the proof consists in bounding $\Sigma_0$ by $\Sigma$ (see (18)). According to Lemma D.4,

$$
\frac{1}{|I|} \left( \sum_{(i,i') \in I^2} |a_2^-(i,i')|^2 \right)^{\frac{1}{2}} \le 2 \|S\|_{\infty} \left( \left( \frac{1}{|I|} \sum_{i \in I} \|f_i \mathbf{1}_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}} + \left( \frac{1}{|I|} \sum_{i \in I} \|\mathcal{F}[f_i] \mathbf{1}_{|\cdot| \ge k_0} \|_2^2 \right)^{\frac{1}{2}} \right).
$$

We also have, according to Lemma D.3.

$$
\frac{1}{|I|} \left( \sum_{(i, i') \in I^2} |a_2^-(i, i')|^2 \right)^{\frac{1}{2}} \le \frac{\|S\|_{\infty}}{|I|^{\frac{1}{2}}} + \|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|f_i \mathbf{1}_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}}.
$$
For the sum related to $a_2^+$, we obtain, using the triangular inequality, the Cauchy-Schwarz inequality, the Plancherel theorem and the fact that $(f_i)_{i \in I}$ is orthogonal:

$$|a_2^+(i, i')| \le \|S\|_{\infty} (\|f_i 1_{\mathbb{R}^d \setminus W}\|_2 + \|f_{i'} 1_{\mathbb{R}^d \setminus W}\|_2) + |\langle \mathcal{F}[f_i], \mathcal{F}[f_{i'}](S(\cdot - k_0) - S(k_0)) \rangle|.$$

Assume that $\beta \le 1$. Then, using again the fact that $(f_i)_{i \in I}$ is orthogonal we have:

$$\frac{1}{|I|} \left( \sum_{(i, i') \in I^2} |a_2^+(i, i')|^2 \right)^{\frac{1}{2}} \le L \left( \frac{1}{|I|^2} \sum_{i \in I} \|f_i\|_{H^\beta}^2 \right)^{\frac{1}{2}} + 2 \|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}}.$$

Assume that $\beta \in (1, 2]$. In that case $\mathcal{F}[f_i]\mathcal{F}[f_{i'}]$ is even. Hence $\langle \mathcal{F}[f_i], \mathcal{F}[f_{i'}]\nabla S(k_0) \cdot k \rangle = 0$. Accordingly, arguing as in the proof of Lemma A.2, we obtain

$$\frac{1}{|I|} \left( \sum_{(i, i') \in I^2} |a_2^+(i, i')|^2 \right)^{\frac{1}{2}} \le \sqrt{dL} \left( \frac{1}{|I|^2} \sum_{i \in I} \|f_i\|_{H^\beta}^2 \right)^{\frac{1}{2}} + 2 \|S\|_{\infty} \left( \frac{1}{|I|} \sum_{i \in I} \|f_i 1_{\mathbb{R}^d \setminus W}\|_2^2 \right)^{\frac{1}{2}}.$$

The previous bounds yields $\Sigma_0 \le \Sigma$, which concludes the proof. $\square$

The next lemma exploits the connectivity condition of partitions arising in the computation of joint cumulants of the squares of the linear statistics associated with the tapers. This condition allows for combinatorial simplifications due to the orthogonality of the tapers. The key point is to provide a bound involving the sums over indices $i, i' \in I$ of the sequences $(a_2^b(i, i'))_{b \in \mathcal{B}_2(\sigma)}$. Indeed, $a_2^b(i, i')$ typically corresponds to a scalar product between the taper $f_i$ and a function depending on $b$. In this context, Bessel's inequality yields a sharper bound on $\sum_{(i, i') \in I^2} |a_2^b(i, i')|^2$ than the trivial bound $\sum_{(i, i') \in I^2} |a_2^b(i, i')|^2 \le |I|^2 \sup_{(i, i')} |a_2^b(i, i')|^2$, for which the use of the combinatorial structure would not have been necessary. However, this refinement is essential to show that the non-asymptotic correlation term $a_2$ in Theorem 5.2 is negligible at sufficiently large frequencies.

**Lemma B.2.** Let $I$ be a discrete subset of $\mathbb{N}^d$. Let $m \ge 2$ and $\sigma$ be a partition of $[2m]$ having block of size at least 2. Assume that $\sigma \vee \tau = \mathbf{1}_m$, where $\tau = 12|\dots|(2m-1)2m$. Denote by $\mathcal{B}_2(\sigma)$ the blocks of size 2 of $\sigma$. For a partition $b \in \mathcal{B}_2(\sigma)$, we denote $b_1 < b_2$ such that $b = (b_1, b_2)$. Let $s \le m$. Consider $P$ a subset of cardinal $s$ of $\mathcal{B}_2(\sigma)$. Finally consider non-negative numbers $(a_n(i))_{n \ge 3, i \in I}$ and $(a_2^b(i, i'))_{b \in \mathcal{B}_2(\sigma), (i, i') \in I^2}$. Denote:

$$B_{\sigma,P} := \sum_{(i_1, \dots, i_m) \in I^m} \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \prod_{j \in b} a_{|b|(i_{\lceil b/2 \rceil})} \prod_{b \in \mathcal{B}_2(\sigma) \setminus P} a_2^b(i_{\lceil b_1/2 \rceil}, i_{\lceil b_2/2 \rceil}) \prod_{b \in P} 1_{i_{\lceil b_1/2 \rceil} = i_{\lceil b_2/2 \rceil}}.$$

Then, $B_{\sigma,P} = |I|$ if $|\sigma| = m = s$. Otherwise, we have

$$B_{\sigma,P} \le \prod_{\substack{b \in \sigma \\ |b| \ge 3}} \left( \sum_{i \in I} |a_{|b|(i)}|^2 \right)^{\frac{|b|}{2}} \prod_{b \in \mathcal{B}_2(\sigma) \setminus P} \left( \sum_{(i, i') \in I^2} |a_2^b(i, i')|^2 \right)^{\frac{1}{2}}. \quad (46)$$

*Proof.* For the proof, we introduce graphs to account for the connectivity condition $\sigma \vee \tau = \mathbf{1}_{[2m]}$. Let $G = G(\sigma)$ be an undirected graph whose vertices are the union of the blocks of $\sigma$ and the blocks of $\tau$. An edge is drawn between two vertices $b$ and $b'$ if $b \cap b' \ne \emptyset$. The condition $\sigma \vee \tau = \mathbf{1}_{[2m]}$
is equivalent to G being connected, that is, having a single connected component (see Section 3.2 of [54]).

For a vertex $v = (v_1, \dots, v_{|v|})$ in $G$, we call $v_1, \dots, v_{|v|}$ the labels of $v$. We now define $G_2$ as a graph having the same edges as $G$, but where each label $n$ appearing in the blocks of $G$ is replaced by $\lfloor n/2 \rfloor$. Since the vertices of $G$ corresponding to blocks of the partition $\tau$ yield in $G_2$ vertices where the two labels are equal, we remove such vertices from $G_2$.

Note that the labels of $G_2$ belong to $\{1, \dots, m\}$. Moreover, for any vertex $v$ of $G_2$ with only two labels $v_1$ and $v_2$, we necessarily have $v_1 \neq v_2$. To see this, assume the contrary. If $v_1 = v_2$, then $v$ comes from a vertex of $G$ of the form $(2n - 1, 2n)$ for some $n \in [m]$, corresponding to a block of $\sigma$. In that case, the block $(2n - 1, 2n)$ is only connected to the vertex of $G$ corresponding to the block $(2n - 1, 2n)$ of $\tau$, and conversely. Therefore, $G$ would not be connected, which contradicts the assumption.

The graph $G_2$ satisfies two key properties: $G_2$ is connected (since $G$ is connected) and each number $n \in [m]$ appears exactly twice in the labels of $G_2$.

For the following, we introduce additional notations. Consider a graph $G_0$ whose labels belong to $[m]$. Denote by $\mathcal{V}_2(G_0)$ the set of vertices of $G_0$ with exactly two labels, and by $\mathcal{V}_3(G_0)$ the set of vertices of $G_0$ with more than three labels. We also denote by $P_\mathcal{V}(G_0)$ the (possibly empty) set of vertices of $G_0$ corresponding to the blocks of $\sigma$ belonging to $P$.

In the following, products of the form $\prod_{j \in v}$, where $v \in \mathcal{V}_2(G_0) \cup \mathcal{V}_3(G_0)$, iterate over the labels of the vertex $v$. We denote by $\mathcal{L}(G_0)$ the set of distinct labels of $G_0$ (without repetition). Finally, for a vertex $v$ of $G_0$ with exactly two labels, we write $v = (v_1, v_2)$ with $v_1 < v_2$. With the previous notations, we can rewrite $B_{\sigma,P}$ as

$$B_{\sigma,P} := \sum_{i_l \in I; l \in \mathcal{L}(G_2)} \prod_{v \in \mathcal{V}_3(G_2)} \prod_{j \in v} a_{|v|}(i_j) \prod_{v \in \mathcal{V}_2(G_2) \setminus P_\mathcal{V}(G_2)} a_2^v(i_{v_1}, i_{v_2}) \prod_{v \in P_\mathcal{V}(G_2)} \mathbf{1}_{i_{v_1} = i_{v_2}}.$$

Consider the case where $\sigma$ contains at least one block of size larger than three, i.e. $|\sigma| < m$.

We proceed inductively. If $P_\mathcal{V}$ is empty, we do nothing. Otherwise, we take into account one constraint $\mathbf{1}_{i_{v_1} = i_{v_2}}$ in the sum (for instance, the one such that $(v_1, v_2)$ is the smallest for the lexicographical order). Specifically, we remove the vertex $v = (v_1, v_2)$ in $G_2$ and change the value of the label $v_1$ appearing in the other vertex $v'$ of $G_2$ by $v_2$. We obtain a new graph $G'_2$. Recall that every label appears only twice in $G_2$, so $v_1$ appears in only one block distinct of $(v_1, v_2)$. The case $v' = (v_2, v_1)$ is impossible since $\sigma$ contains at least one block of size larger than three. Indeed, otherwise $v, v'$ would be a connected component of $G_2$, distinct from $G_2$, which contradicts the fact that $G_2$ is connected. Hence, $v$ was connected to two (and only two) distinct vertices $v'$ and $v''$. We add in $G'_2$ an edge between $v'$ and $v''$.

With that procedure, $G'_2$ remains connected. Moreover, $G'_2$ has one fewer block than $G_2$ and contains the labels of $G_2$, except for $v_1$. Since we removed and added once the label $v_2$, every label in $G'_2$ appears exactly twice. Finally, we have $|P_\mathcal{V}(G'_2)| = |P_\mathcal{V}(G_2)| - 1$. The previous discussion yields

$$B_{\sigma,P} = \sum_{i_l \in I, l \in \mathcal{L}(G'_2)} \prod_{v \in \mathcal{V}_3(G'_2)} \prod_{j \in v} a_{|v|}(i_j) \prod_{v \in \mathcal{V}_2(G'_2) \setminus P_\mathcal{V}(G'_2)} a_2^v(i_{v_1}, i_{v_2}) \prod_{v \in P_\mathcal{V}(G'_2)} \mathbf{1}_{i_{v_1} = i_{v_2}}.$$

Note that the number of terms in the product over $\mathcal{V}_2(G'_2) \setminus P_\mathcal{V}(G'_2)$ is unchanged and still equals $|B_2(\sigma)| - s$. Since $G'_2$ satisfies the same properties as $G_2$ that we used previously, we can iterate this procedure to obtain

$$B_{\sigma,P} = \sum_{i_l \in I, l \in \mathcal{L}(G_2^s)} \prod_{v \in \mathcal{V}_3(G_2^s)} \prod_{j \in v} a_{|v|}(i_j) \prod_{v \in \mathcal{V}_2(G_2^s) \setminus P_\mathcal{V}(G_2^s)} a_2^v(i_{v_1}, i_{v_2}),$$
where $G_2^s$ is a graph having $m-s$ labels, each appearing exactly twice. To upper bound the last quantity, we establish the following intermediate result. Let $N \ge 2$ and $(c^n(i, i'))_{n \le N, (i, i') \in I^2}$ a sequence of non-negative numbers. Then

$$ \sum_{(i_1, \dots, i_N) \in I^N} \prod_{n=1}^{N} c^n(i_n, i_{n+1}) \le \prod_{n=1}^{N} \left( \sum_{(i, i') \in I^2} c^n(i, i')^2 \right)^{\frac{1}{2}}, \quad (47) $$

with the convention $i_{N+1} = i_1$. We prove it by induction. For $N=2$, it follows from the Cauchy-Schwarz inequality:

$$ \sum_{i_1, i_2 \in I^N} c^1(i_1, i_2) c^2(i_2, i_1) \le \sum_{i_1 \in I} \left( \sum_{i' \in I} c^1(i_1, i')^2 \right)^{\frac{1}{2}} \left( \sum_{i \in I} c^2(i, i_1)^2 \right)^{\frac{1}{2}} \\ \le \left( \sum_{(i, i') \in I^2} c^1(i, i')^2 \right)^{\frac{1}{2}} \left( \sum_{(i, i') \in I^2} c^2(i, i')^2 \right)^{\frac{1}{2}}. $$

Fix $N \ge 2$. Assume that the inequality (47) is true when the product of its left-hand side is over $N$ terms. Then, by induction we have

$$ \sum_{(i_1, \dots, i_{N+1}) \in I^{N+1}} \prod_{n=1}^{N+1} c^n(i_n, i_{n+1}) = \sum_{(i_1, \dots, i_N) \in I^N} \prod_{n=1}^{N-1} c^n(i_n, i_{n+1}) \left( \sum_{i_{N+1} \in I} c^N(i_N, i_{N+1}) c^{N+1}(i_{N+1}, i_1) \right) \\ \le \prod_{n=1}^{N-1} \left( \sum_{(i, i') \in I^2} c^n(i, i')^2 \right)^{\frac{1}{2}} \left( \sum_{(i, i') \in I^2} \left( \sum_{i_{N+1} \in I} c^N(i, i_{N+1}) c^{N+1}(i_{N+1}, i') \right)^2 \right)^{\frac{1}{2}} \\ \le \prod_{n=1}^{N+1} \left( \sum_{(i, i') \in I^2} c^n(i, i')^2 \right)^{\frac{1}{2}}, $$

which proves (47).

We now have the tools to upper bound $B_{\sigma,P}$. We introduce

$$ \mathcal{L}_3 := \{ l \in \mathcal{L}(G_2^s) | l \text{ appears only in vertices having more than 3 labels} \}. $$

Then, $B_{\sigma,P} = B_{\sigma,P}^3 B_{\sigma,P}^{23}$ where

$$ B_{\sigma,P}^3 := \sum_{i_l; l \in \mathcal{L}_3} \prod_{v \in \mathcal{V}_3(G_2^s)} \prod_{j \in v} a_{|v|(i_j)}^{1_j \in \mathcal{L}_3}, $$

$$ B_{\sigma,P}^{23} := \sum_{i_l; l \in \mathcal{L} \setminus \mathcal{L}_3} \prod_{v \in \mathcal{V}_3(G_2^s)} \prod_{j \in v} a_{|v|(i_j)}^{1_j \in \mathcal{L} \setminus \mathcal{L}_3} \prod_{v \in \mathcal{V}_2(G_2^s) \setminus \mathcal{P}_v(G_2^s)} a_2^{v}(i_{v_1}, i_{v_2}), $$

We now bound $B_{\sigma,P}^3$. For each label $l \in \mathcal{L}_3$, denote by $v(l)$ and $v'(l)$ the (possibly identical) vertices of size larger than three in which $l$ appears. Using the Cauchy-Schwarz inequality, we get:

$$ B_{\sigma,P}^3 = \prod_{l \in \mathcal{L}_3} \sum_{i \in I} a_{|v(l)|}(i) a_{|v'(l)|}(i) \le \prod_{l \in \mathcal{L}_3} \left( \sum_{i \in I} |a_{|v(l)|}(i)|^2 \right)^{\frac{1}{2}} \left( \sum_{i \in I} |a_{|v'(l)|}(i)|^2 \right)^{\frac{1}{2}}. $$
Next, consider the following product:

$$p := \prod_{v \in \mathcal{V}_3(G_2^s)} \prod_{j \in v} a_{|v|}(i_j)^{1_j} \in \mathcal{L} \setminus \mathcal{L}_3.$$

The number of terms in $p$ is equal to $2m - 2|\mathcal{B}_2(\sigma)| - 2|\mathcal{L}_3|$. Since this number is even, we can group the terms of $p$ two by two. Moreover, because each label appearing in $p$ appears exactly once, we can perform this grouping such that each pair contains distinct labels.

Using the fact that each label of $G_2^s$ appears exactly twice, we can rewrite $B_{\sigma,P}$ in the form of the left-hand side of (47). Applying this result yields:

$$B_{\sigma,P}^{23} \le \prod_{l \in \mathcal{L}_{23}} \left( \sum_{i \in I} |a_{|v_3(l)|}(i)|^2 \right)^{\frac{1}{2}} \prod_{v \in \mathcal{V}_2(G_2^s) \setminus P_v(G_2^s)} \left( \sum_{(i,i') \in I^2} a_2^v(i, i')^2 \right)^{\frac{1}{2}}, \quad (48)$$

where

$$\mathcal{L}_{23} := \{l \in \mathcal{L}(G_2^s) \mid l \text{ appears only in one vertex of size 2}\},$$

and for each $l \in \mathcal{L}_{23}$, $v_3(l)$ denotes the vertex of size larger than three in which it appears. Note that, to obtain (48), we used, after applying (47), the fact that for vertices of size larger than three:

$$\sum_{(i',i) \in I^2} (a_n(i)^2 a_{n'}(i'))^2 = \left(\sum_{i \in I} a_n(i)^2\right) \left(\sum_{i' \in I} a_{n'}(i')^2\right), \quad \text{for } n, n' \ge 3,$$

together with the fact that the grouping of terms in $p$ is performed so that each pair contains distinct labels.

Since the effect on $\mathcal{V}_2(G_2) \setminus \mathcal{P}(G_2)$ and $\mathcal{V}_3(G_2)$ of the procedure to go from $G_2$ to $G_2^s$ was only to change their labels, we have

$$\prod_{v \in \mathcal{V}_2(G_2^s) \setminus P_v(G_2^s)} \left( \sum_{(i, i') \in I^2} a_2^v(i, i')^2 \right)^{\frac{1}{2}} = \prod_{b \in \mathcal{B}_2(\sigma)} \left( \sum_{(i, i') \in I^2} a_2^b(i, i')^2 \right)^{\frac{1}{2}},$$

and

$$\prod_{l \in \mathcal{L}_3} \left( \sum_{i \in I} |a_{|v(l)|}(i)|^2 \right)^{\frac{1}{2}} \left( \sum_{i \in I} |a_{|v'(l)|}(i)|^2 \right)^{\frac{1}{2}} \prod_{l \in \mathcal{L}_{23}} \left( \sum_{i \in I} |a_{|v_3(l)|}(i)|^2 \right)^{\frac{1}{2}} \\ = \prod_{b \in \sigma; \ j \in b} \left( \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{1}{2}} = \prod_{b \in \sigma; \ |b| \ge 3} \left( \sum_{i \in I} a_{|b|}(i)^2 \right)^{\frac{|b|}{2}}.$$

It yields (46) when $|\sigma| < m$.

Now, assume that $\sigma$ contains only blocks of size two, i.e. $|\sigma| = m$. Then,

$$B_{\sigma,P} = \sum_{i_l \in I; l \in \mathcal{L}(G_2)} \prod_{v \in \mathcal{V}_2(G_2) \setminus P_v(G_2)} a_2^v(i_{v_1}, i_{v_2}) \prod_{v \in P_v(G_2)} 1_{i_{v_1} = i_{v_2}}.$$

Using similar connectivity arguments as in the case $|\sigma| < m$, we reduce by $s$ the dimension $|\mathcal{L}(G_2)|$ of the sum over the indices $i_l$ for $l \in \mathcal{L}(G_2)$, by introducing a sequence of graphs $G_2^2, \dots, G_2^s$, provided
that $G_2^2, \dots, G_2^{s-1}$ do not reduce to two vertices. This situation cannot occur if $s \le m-1$, since only one vertex is removed at each step and we start with a graph having $m$ vertices. Accordingly, if $s \le m-1$, using (47), we obtain (46), as in the case $|\sigma| = m$. Finally, when $m = s$, the term $B_{\sigma,P}$ reduces to

$$B_{\sigma,P} = \sum_{i_l \in I, l \in \mathcal{L}(G_2)} \prod_{v \in G_2} \mathbb{1}_{i_{v_1} = i_{v_2}}.$$

Reducing the dimension of the sum by $s-1$, using the previously mentioned method, yields:

$$B_{\sigma,P} = \sum_{i_1 \in I, i_2 \in I} \mathbb{1}_{i_1=i_2} \mathbb{1}_{i_2=i_1} = |I|.$$

This concludes the proof of the lemma.

<figure><img src="image_9.png" /></figure>

The next lemma shows that the Brillinger-mixing assumption on the total variation of reduced factorial cumulant measures extends to the total variation of the cumulant measures of linear statistics. This result is specific to point processes, where the discrete nature of the linear statistics implies that their cumulants of order $m \ge 2$ involve factorial cumulants of orders less than or equal to $m$.

**Lemma B.3.** Let $s \ge 1$ and $q, \gamma > 0$. Then

$$\sum_{\pi \in \Pi[s]} q^{\lvert \pi \rvert - 1} ((\lvert \pi \rvert - 1)!)^{1+\gamma} \le 2^s (1+q)^{s-1} ((s-1)!)^{1+\gamma}.$$

*Proof.* We denote by $S(s,l)$ the Stirling number of the second kind. We have

$$\sum_{\pi \in \Pi[s]} q^{\lvert \pi \rvert - 1} ((\lvert \pi \rvert - 1)!)^{1+\gamma} \le ((s-1)!)^{\gamma} \sum_{l=1}^{s} (l-1)! q^{l-1} S(s,l).$$

According to equation (4.5) of [71], $S(s,l) = \sum_{n_1+\dots+n_l=s} \frac{s!}{n_1!\dots n_l!l!}$, where the sum is over the decompositions of the integer $s$ into positive integers. Hence,

$$\sum_{\pi \in \Pi[s]} q^{\lvert \pi \rvert - 1} ((\lvert \pi \rvert - 1)!)^{1+\gamma} \le ((s-1)!)^{1+\gamma} s \sum_{l=1}^{s} \frac{1}{l} q^{l-1} \sum_{n_1+\dots+n_l=s} \frac{1}{n_1! \dots n_l!} \\ \le ((s-1)!)^{1+\gamma} s \sum_{l=1}^{s} \frac{1}{l} q^{l-1} \sum_{n_1+\dots+n_l=s} 1 = ((s-1)!)^{1+\gamma} s \sum_{l=1}^{s} \frac{1}{l} q^{l-1} \binom{s-1}{l-1},$$

where in the second line we used equation (4.2) of [71]. To conclude, we use

$$s \sum_{l=1}^{s} \frac{1}{l} q^{l-1} \binom{s-1}{l-1} = \sum_{l=1}^{s} q^{l-1} \binom{s}{l} \le 2^s (1+q)^{s-1}.$$

<figure><img src="image_10.png" /></figure>

The following lemma controls the cumulants appearing on the right-hand side of (43), by introducing suitable random variables: a Gaussian distribution to handle non-asymptotic correlations between tapers in second-order cumulants, and a Gamma distribution to account for higher-order cumulants. Its key idea is to separate second-order contributions from higher-order ones.
**Lemma B.4.** Let $\Phi$ satisfying Assumption 3.1 with intensity $\lambda = 1$, $I$ be a discrete subset of $\mathbb{N}$ and $(f_i)_{i \in I}$ be a family of $L^2(\mathbb{R}^d)$. Assume that $\Phi$ satisfies Assumption 5.1 with parameters $(c, q, \gamma) \in (0, \infty)^3$. Assume (17). We denote

$$ \mathcal{A} := \sum_{\substack{\sigma \in \Pi[2m], \\ \sigma \vee \tau = \mathbf{1}_{2m}, |\sigma| \ge 2}} \prod_{b \in \sigma} K_{|b|}, \quad (49) $$

where $(K_s)_{s \ge 0}$ is defined in (44). If (15) is satisfied with $c_f \in (1, \infty)$, $\eta_1, \eta_2 \ge 0$ and $\rho > 0$, then

$$ \mathcal{A} \le \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \left( \frac{m!}{2\sqrt{2\pi}} (4\Sigma)^m + \frac{(m!)^2}{2\sqrt{2\pi}} \left( 32 \frac{\Lambda^2 c^{1/2} (1+q)^{3/2} \log(|I|^{\frac{1}{d}})^{\frac{1}{2}}}{\rho^{d/2} |I|^{2\eta_1+\eta_2/2}} \right)^m \right), \quad (50) $$

where $\Lambda$ is defined in (19). If (15) is satisfied with $c_f \in (1, \infty)$, $\eta_1, \eta_2 \ge 0$ and $\rho > 0$, then

$$ \mathcal{A} \le \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \left( \frac{m!}{2\sqrt{2\pi}} (4\Sigma)^m + \frac{(m!)^2}{2\sqrt{2\pi}} \left( 32 \frac{[c_f]^2 c^{2/3} (1+q)^{4/3}}{\rho^{d/3} |I|^{2\eta_1}} \right)^m \right). \quad (51) $$

*Proof.* Using the fact that $\Phi$ satisfies Assumption 5.1, Lemma B.3, and the inequality $k_1! \dots k_n! \le (k_1 + \dots + k_n)!$ for all integers $k_1, \dots, k_n$, we obtain

$$ \prod_{|b| \ge 3} \sum_{\pi \in \Pi[|b|]} |\gamma_{\text{red}}^{(|\pi|)}| (\mathbb{R}^{d(|\pi|-1)}) \le \prod_{|b| \ge 3} c 2^{|b|} (1+q)^{|b|-1} ((|b|-1)!)^{1+\gamma} \\ \le ((2m)!)^{\gamma} \prod_{|b| \ge 3} c 2^{|b|} (1+q)^{|b|-1} (|b|-1)! $$

We define $\kappa_1 = 0$, $\kappa_2 = \Sigma$ and

$$ \kappa_3 = c 2^4 (1+q)^2 \left( 3 \left( \sup_{i \in I} \|f_i \mathbf{1}_{\mathbb{R}^d \setminus W} \|_3 \sup_{i \in I} \|f_i\|_3^2 + \sup_{i \in I} \|\mathcal{F}[f_i] \mathbf{1}_{|\cdot| \ge k_0/3} \|_1 \right) \wedge \frac{c_f}{\rho^{d/2} |I|^{3\eta_1}} \right), \\ \forall s \ge 4, \quad \kappa_s = c_f \frac{c \rho^d \log(|I|^{\frac{1}{d}})}{(1+q)|I|^{\eta_2}} \left( 2 \frac{1+q}{\rho^{\frac{d}{2}} |I|^{\eta_1}} \right)^s (s-1)! $$

Then, using Assumptions (15) and (16), we get

$$ \mathcal{A} \le \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \vee \tau = \mathbf{1}_{2m}}} \prod_{b \in \sigma} \kappa_{|b|}, \quad (52) $$

where in the last line we used the fact that for $m \ge 2$, $(2m)! \le 4^m (m!)^2 / \sqrt{2\pi}$, obtained by Stirling approximation [33].

We prove (50). To do so, we consider

$$ X = \Sigma^{\frac{1}{2}} \mathcal{N} + 2 \frac{(c(1+q)^3)^{1/4} \log(|I|^{\frac{1}{d}})^{\frac{1}{4}}}{\rho^{d/4} |I|^{\eta_1+\eta_2/4}} \Gamma, $$

where $\mathcal{N}$ is a standard Gaussian random variable and $\Gamma$ is an Gamma random variable with parameters $(\Lambda, 1)$. Moreover, $\mathcal{N}$ and $\mathcal{E}$ are independents. We also denote $\bar{X} = X - \mathbb{E}[X]$. We introduced $\bar{X}$ because $\kappa_s \le \kappa_s(\bar{X})$, for all $s \ge 1$. Indeed, $\kappa_1 = \kappa_1(\bar{X}) = 0$,

$$ \kappa_2 = \Sigma \le \Sigma + \Lambda \left( 2 \frac{(c(1+q)^3)^{1/4} \log(|I|^{\frac{1}{d}})^{\frac{1}{4}}}{\rho^{d/4} |I|^{\eta_1+\eta_2/4}} \right)^2 = \kappa_2(\bar{X}). $$
For $s = 3$, we have

$$ \kappa_3 = \kappa_3 \left( \frac{\rho^{d/4} |I|^{\eta_1+\eta_2/4}}{2(c(1+q)^3)^{1/4} \log\left(|I|^{1/4}\right)} \right)^3 \frac{1}{2} \left( 2 \frac{(c(1+q)^3)^{1/4} \log\left(|I|^{1/4}\right)}{\rho^{d/4} |I|^{\eta_1+\eta_2/4}} \right)^3 2! \le \kappa_3(\bar{X}), $$

where in the second line we used $c\rho^d \log(|I|^{1/4}) \ge (1+q)|I|^{\eta_2}$. Finally, for $s \ge 4$, using again $c\rho^d \log(|I|^{1/4}) \ge (1+q)|I|^{\eta_2}$, we have

$$ \kappa_s = c_f \frac{c \rho^d \log(|I|^{1/4})}{(1+q)|I|^{\eta_2}} \left( 2 \frac{1+q}{\rho^{d/2} |I|^{\eta_1}} \right)^s (s-1)! \le c_f \left( \frac{c \rho^d \log(|I|^{1/4})}{(1+q)|I|^{\eta_2}} \right)^{s/4} \left( 2 \frac{1+q}{\rho^{d/2} |I|^{\eta_1}} \right)^s (s-1)! \le \kappa_s(\bar{X}). $$

Hence, with equation 3.2.6 of [61], we obtain

$$ A \le \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \sum_{\substack{\sigma \in \Pi_2[2m], \\ \sigma \nmid \tau = \mathbf{1}_{2m}}} \prod_{b \in \sigma} \kappa_s(\bar{X}) \le \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \mathbb{E}[\bar{X}^{2m}]. \quad (53) $$

Let $a = \Sigma^{1/2}$ and $b = 2(c(1+q)^3)^{1/4} \log(|I|^{1/4})^{1/4} / (\rho^{d/4} |I|^{\eta_1+\eta_2/4})$. We have

$$ \begin{aligned} \mathbb{E}[\bar{X}^{2m}] &= \mathbb{E} \left[ (a\mathcal{N} + b(\Gamma - \mathbb{E}[\Gamma]))^{2m} \right] = \sum_{j=0}^{2m} \binom{2m}{j} a^j b^{2m-j} \mathbb{E}[\mathcal{N}^j] \mathbb{E}[(\Gamma - \mathbb{E}[\Gamma])^{2m-j}] \\ &= \sum_{j=0}^{m} \binom{2m}{2j} a^{2j} b^{2m-2j} \mathbb{E}[\mathcal{N}^{2j}] \mathbb{E}[(\Gamma - \mathbb{E}[\Gamma])^{2m-2j}]. \end{aligned} $$

Let $(\mathcal{E}_n)_{n \in [\Lambda]}$ be i.i.d. standard exponential random variables. Then, for even $j \ge 0$, we have

$$ \begin{aligned} \mathbb{E}[(\Gamma - \mathbb{E}[\Gamma])^j]^{1/j} &= \mathbb{E} \left[ \left( \sum_{n=1}^{\Lambda} \mathcal{E}_n - \mathbb{E}[\mathcal{E}_n] \right)^j \right]^{1/j} \le \sum_{n=1}^{\Lambda} \mathbb{E} \left[ (\mathcal{E}_n - \mathbb{E}[\mathcal{E}_n])^j \right]^{1/j} \\ &= \Lambda \left( j! \sum_{n=0}^{j} \frac{(-1)^n}{n!} \right)^{1/j} \le \Lambda (j!)^{1/j}. \end{aligned} $$

Using the expression of the even moments of the Gaussian distribution, we obtain

$$ \mathbb{E}[\bar{X}^{2m}] \le \sum_{j=0}^{m} \binom{2m}{2j} a^{2j} b^{2m-2j} \frac{(2j)!}{2^j j!} \Lambda^{2m-2j} (2m-2j)! = (2m)!(b\Lambda)^{2m} \sum_{j=0}^{m} \frac{1}{j!} d^j, $$

where $d = a^2/(2b^2\Lambda^2)$. Using the binomial theorem, we get

$$ \sum_{j=0}^{m} \frac{1}{j!} d^j = \frac{1}{m!} \int_{0}^{\infty} (d+t)^m e^{-t} dt \le \frac{2^{m-1}}{m!} (d^m + m!). $$

Accordingly,

$$ \begin{aligned} \mathbb{E}[\bar{X}^{2m}] &\le \frac{(2m)!}{m!} 2^{m-1} \frac{a^{2m}}{2^m} + (2m)!(b\Lambda)^{2m} 2^{m-1} \\ &\le \frac{m!}{2\sqrt{2\pi}} (4\Sigma)^m + \frac{(m!)^2}{2\sqrt{2\pi}} \left( 32 \frac{\Lambda^2 c^{1/2} (1+q)^{3/2} \log(|I|^{1/4})^{1/2}}{\rho^{d/2} |I|^{2\eta_1+\eta_2/2}} \right)^m, \end{aligned} $$
which concludes the proof of (50), recalling (53). The proof of (51) is done in a similar way, by considering:

$$X = \Sigma^{\frac{1}{2}} \mathcal{N} + 2 \frac{(c(1+q)^2)^{1/3}}{\rho^{d/6} |I|^{\eta_1}} \Gamma,$$

where $\mathcal{N}$ is a standard Gaussian random variable and $\Gamma$ is an Gamma random variable with parameters $(\lceil c_f \rceil, 1)$. $\square$

To prove Theorem 5.2, we combine the previous lemmas with Theorem G.1.

**Proof of Theorem 5.2.** According to lemma B.1 and Lemma B.4, we have

$$|\kappa_m(\hat{S}(k_0))| \le A_1 + A_2 + A_3 \wedge A_4,$$

where

$$A_1 := 2^{m-1}(m-1)! \frac{S(k_0)^m}{|I|^{m-1}}, \quad A_2 := \left(\frac{4^m m!}{\sqrt{2\pi}}\right)^\gamma \frac{m!}{2\sqrt{2\pi}} (4\Sigma)^m,$$

$$A_3 := \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \frac{(m!)^2}{2\sqrt{2\pi}} \left( 32 \frac{\Lambda^2 c^{1/2} (1+q)^{3/2} \log(|I|^{1/2})^{1/2}}{\rho^{d/2} |I|^{2\eta_1+\eta_2/2}} \right)^m,$$

$$A_4 := \left( \frac{4^m m!}{\sqrt{2\pi}} \right)^\gamma \frac{(m!)^2}{2\sqrt{2\pi}} \left( 32 \frac{[c_f]^2 c^{2/3} (1+q)^{4/3}}{\rho^{d/3} |I|^{2\eta_1}} \right)^m.$$

We have, since $m \ge 2$,

$$A_1 \le \frac{m!}{2} \left( \frac{2S(k_0)}{|I|} \right)^{m-2} \frac{2S(k_0)^2}{|I|}.$$

Then, after computations, we get

$$A_2 \le \left( \frac{m!}{2} \right)^{1+\gamma} \frac{16}{\sqrt{2\pi}} \left( \frac{32}{\sqrt{2\pi}} \right)^{\gamma} \Sigma^2 (4^{1+\gamma} \Sigma)^{m-2}.$$

Finally,

$$A_3 \wedge A_4 \le \left( \frac{m!}{2} \right)^{2+\gamma} \frac{2^{11}}{\sqrt{2\pi}} \left( \frac{32}{\sqrt{2\pi}} \right)^{\gamma} B_{\infty}^2 (32B_{\infty})^{m-2}.$$

Theorem G.1 and numerical computations conclude the proof. $\square$

Using the results of Section F, we specialize Theorem 5.2 to Hermite tapers.

**Proof of Corollary 5.2.** In the following, $K$ denotes a finite constant that may change from line to line. We have:

$$|\hat{S}(k_0) - S(k_0)| \le |\hat{S}(k_0) - \mathbb{E}[\hat{S}(k_0)]| + |S(k_0) - \mathbb{E}[\hat{S}(k_0)]|.$$

According to Lemma A.2 and arguing as in the proof of Corollary 5.1, the bias term is bounded as:

$$|S(k_0) - \mathbb{E}[\hat{S}(k_0)]| \le \frac{K}{|W|^{\beta/(2\beta+d)}}.$$
We consider the variance term $|\hat{S}(k_0) - \mathbb{E}[\hat{S}(k_0)]|$. According to Lemma F.3, Assumptions (15) and (16) are satisfied with $\rho = r$, $\eta_1 = 1/12$ and $\eta_2 = 1/6$. To apply Theorem 5.2, we check that condition (17) is satisfied.

$$ \frac{(1+q)|I|^{\eta_2}}{c\rho^d \log(|I|^{\frac{1}{d}})} \le \frac{\sqrt{2(1+q)}|I|^{\frac{2}{3}}}{c\log(2^{\frac{1}{d}})|W|} \le \frac{\sqrt{2(1+q)}c_3^{\frac{2}{3}}}{c\log(2^{\frac{1}{d}})}|W|^{\frac{2}{3}\frac{2\beta}{2\beta+d}-1} = \frac{\sqrt{2(1+q)}c_3^{\frac{2}{3}}}{c\log(2^{\frac{1}{d}})}|W|^{-\frac{2\beta+3d}{6\beta+3d}} \le 1, $$

provided that $|W|$ is large enough, say $|W| \ge v_0$ for some $v_0 > 0$. Moreover, with similar computations, we obtain

$$ \frac{1+q}{c\rho^d} \le \sqrt{2}\frac{1+q}{c}|W|^{-\frac{d+\beta}{2\beta+d}} \le 1, $$

by increasing $v_0$ if necessary. Accordingly, one can apply Theorem 5.2. It remains to bound $\Sigma$ and $B_\infty$. With the arguments of the proof of Corollary 5.1 (see Section A.3), we obtain

$$ \begin{aligned} \Sigma &\le \frac{S(k_0)}{|I|} + K \left( \frac{L}{|I|^{\frac{1}{2}}} \left( \frac{|I|}{|W|} \right)^{\frac{\beta}{d}} + \|S\|_{\infty} e^{-c'|I|^{\frac{3\theta}{2d}}} \right) + 2 \|S\|_{\infty} \sup_{i \in I} \|\psi_i \mathbf{1}_{|\cdot| \ge k_0}\|_2 \\ &\le K \frac{L + \|S\|_{\infty}}{|W|^{2\beta/(2\beta+d)}} + K \|S\|_{\infty} \left( \sup_{i \in I} \|\psi_i \mathbf{1}_{|\cdot| \ge k_0}\|_2 \wedge \frac{1}{|W|^{\beta/(2\beta+d)}} \right). \end{aligned} $$

Note that when $\beta \ge 1$, we use the fact that $|I| \ge c_d\{|i \in \mathbb{N}^d | |i|_\infty < i_{\max}\}$, for a constant $c_d > 0$ depending only on the dimension $d$. Assume that $|k_0| \ge c_3^{\frac{1}{d}} 2|W|^{-\frac{1}{2\beta+d}}$. Then, using assumptions (21) and (22), we prove that $r|k_0| \ge \sqrt{2}|I|^{1/(2d)} \ge \sqrt{i_{\max} + i_{\max}^{1/3+\theta}}$. So, we can apply Lemma F.1, to get

$$ \sup_{i \in I} \|\mathcal{F}[f_i]\mathbf{1}_{|\cdot| \ge k_0}\|_2 \le K e^{-c'(r|k_0|)^{3\theta}} \le K e^{-c'(\sqrt{2}|I|^{\frac{1}{2d}})^{3\theta}} \le \frac{K}{|W|^{2\beta/(2\beta+d)}}. $$

Accordingly, if $|k_0| \ge c_3^{\frac{1}{d}} 2|W|^{-\frac{1}{2\beta+d}}$, then $\Sigma \le K(L + \|S\|_{\infty})/|W|^{2\beta/(2\beta+d)}$. Without any assumption on $k_0$, we have $\Sigma \le K(L + \|S\|_{\infty})/|W|^{\beta/(2\beta+d)}$. To upper bound $B_\infty$, we start with the case $|k_0| \ge c_3^{\frac{1}{d}} 2|W|^{-\frac{1}{2\beta+d}}$. Using again Lemma F.1 and assumptions (21), (22), we prove that

$$ \sup_{i \in I} \|\psi_i \mathbf{1}_{\mathbb{R}^d \setminus W}\|_3 \sup_{i \in I} \|\psi_i\|_3^2 + \sup_{i \in I} \|\psi_i \mathbf{1}_{|\cdot| \ge |k_0|/3}\|_1 \le K e^{-K_1 |W|^\alpha}, $$

for positive constants $\alpha > 0$ and $K_1 > 0$. Subsequently, with assumptions (21) and (22), we obtain, for some $\alpha' > 0$,

$$ \Lambda \le K \left[ c_f \vee e^{-K_1 |W|^\alpha} |W|^{\alpha'} \right] \le K. $$

This yields, using $2\eta_1 + \eta_2/2 = 1/4$,

$$ B_{\infty} \le K \frac{c^{1/2}(1+q)^{3/2} \log(|W|)^{\frac{d}{2}}}{\rho^{d/2}|I|^{1/4}} \le K \frac{c^{1/2}(1+q)^{3/2} \log(|W|)^{\frac{d}{2}}}{|W|^{\frac{1}{2}}}. $$

Finally, without any assumption on $k_0$, we get

$$ B_{\infty} \le K \frac{c_3^{\frac{2}{3}}(1+q)^{\frac{4}{3}}}{r^{d/3}|I|^{\frac{1}{6}}} \le K \frac{c_3^{\frac{2}{3}}(1+q)^{\frac{4}{3}}}{|W|^{\frac{1}{3}}}. $$

This concludes the proof of this lemma.

□
C Proofs related to Section 6

To prove Theorem 6.1, we introduce the following convenient transformation of the data-driven criterion. With the notations of Section 6 and 5, we consider, for all $I \subset \mathbb{N}^d$:

$$ \Delta_I := \hat{\mathcal{R}}(I) - \frac{1}{N} \sum_{j=1}^{N} (\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))^2 - \frac{1}{N} \sum_{j=1}^{N} \hat{S}_p(\tilde{k}_j)^2 - 2\hat{S}_p(\tilde{k}_j)S_{\frac{1}{2}}(\tilde{k}_j). \quad (54) $$

For the following proofs, it is convenient to rewrite $\Delta_I$ as:

$$ \Delta_I = \frac{2}{N} \sum_{j=1}^{N} (\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))(S_{\frac{1}{2}}(\tilde{k}_j) - \hat{S}_p(\tilde{k}_j)), \quad (55) $$

provided that $S(k_j) = S(\tilde{k}_j)$ for all $j \in [N]$. Theorem 6.1 is a direct consequence of the following lemma.

**Lemma C.1.** Suppose the assumptions of Theorem 6.1 hold. Then,

$$ \mathbb{E} \left[ \max_{I \in \mathcal{I}} |\Delta_I| \right] = o \left( |\mathcal{I}|^{\frac{1}{2}} |W|^{-\frac{2\beta}{2\beta+d}} \right), \text{ as } |W| \to \infty. $$

The proof of the previous lemma is the main objective of this section. Before proceeding with it, we first show how to deduce Theorem 6.1.

**Proof of Theorem 6.1.** By definition, we have, for all $I \in \mathcal{I}$, $\hat{\mathcal{R}}(\hat{I}) \le \hat{\mathcal{R}}(I)$. This yields:

$$ \begin{aligned} \frac{1}{N} \sum_{j=1}^{N} (\hat{S}_{\hat{I}}(k_j) - S_{\frac{1}{2}}(k_j))^2 &\le \frac{1}{N} \sum_{j=1}^{N} (\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))^2 + \Delta_I - \Delta_{\hat{I}} \\ &\le \frac{1}{N} \sum_{j=1}^{N} (\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))^2 + 2 \max_{I' \in \mathcal{I}} |\Delta_{I'}|. \end{aligned} $$

Taking expectations in the previous inequality and applying Lemma C.1 completes the proof. □

To prove Lemma C.1, we upper bound both the expectation and the variance of $\Delta_I$. We introduce the next notation to highlight the dominant terms in the forthcoming bounds.

**Definition 7.** In this section, $\mathcal{E}(|I|, |I_p|)$ denotes a non negative exponentially small quantity such that

$$ \mathcal{E}(|I|, |I_p|, \theta) \le K \left( e^{-|I_p|^{3\theta/(2d)}/K} + e^{-|I|^{3\theta/(2d)}/K} \right), $$

for a constant $K \in (0, \infty)$ and where $\theta \in (0, 2/3)$ is the parameter of equation (27). When the context is clear, we use the abbreviation $\mathcal{E} = \mathcal{E}(|I|, |I_p|, \theta)$.

The next four lemmas provide upper bounds on the terms appearing in $\mathbb{E}[\Delta_I]$ and $\text{Var}[\Delta_I]$. The following computations rely on the two general underlying ideas, summarized as follows:

1. In general, cumulants of order $m \ge 4$ of linear statistics associated with a family of tapers are controlled by the product of the $L^m(\mathbb{R}^d)$-norms of the tapers (see Lemma D.2). However, when one of the Hermite tapers $\psi_i$ has an index $i$ appearing only once in a cumulant, and when we sum the square of this cumulant over $i \in I$, we instead use Lemma D.6 to exploit the orthogonality of the tapers.
2. For cumulants of order two or three, two cases are distinguished. If the cumulant involves tapers index $i \in I$ and/or $i_p \in I_p$ with $|I|$ and $|I_p|$ small enough we exploit frequency decorrelation via Lemmas D.4 and D.5. Otherwise, we take advantage of the orthogonality between tapers using Lemmas D.3 and D.6.

**Lemma C.2.** Suppose that the assumptions of Theorem 6.1 hold. Then,

$$ \forall I \in \mathcal{I}, \quad |\mathbb{E}[\Delta_I]| \le K \left( \left( \frac{|I||I_p|}{|W|^2} \right)^{\frac{\beta}{d}} + \frac{\log(|W|)^d}{|W|} + \frac{1}{|W||k_0|^d} + \mathcal{E} \right), $$

where $K < \infty$ does not depend on $I, r, W, N$ and $|k|$.

*Proof.* According to (55),

$$ \Delta_I = \frac{2}{N} \sum_{j=1}^{N} (\hat{S}_I(k_j) - S_{\frac{1}{2}}(k_j))(S_{\frac{1}{2}}(\tilde{k}_j) - \hat{S}_p(\tilde{k}_j)). $$

For $s \in \{1, 2\}$, denote by $C_i^s$ the quantity of equation (9) computed with $\Phi_s$. Hence,

$$ \mathbb{E}[\Delta_I] = \frac{2}{N} \sum_{j=1}^{N} \frac{1}{|I||I_p|} \sum_{(i_1, i_2) \in I \times I_p} \mathbb{E}[(|C_{i_1}^1(k_j)|^2 - S_{\frac{1}{2}}(k_j))(S_{\frac{1}{2}}(\tilde{k}_j) - |C_{i_2}^2(\tilde{k}_j)|^2)]. $$

Using Corollary G.1, we get

$$ |\mathbb{E}[\Delta_I]| \le \frac{2}{N} \sum_{j=1}^{N} \frac{1}{|I||I_p|} \sum_{(i_1, i_2) \in I \times I_p} (A + B + C + D), $$

where

$$ A := |\mathbb{E}[|C_{i_1}^1(k_j)|^2] - S_{\frac{1}{2}}(k_j)|\mathbb{E}[|C_{i_2}^2(\tilde{k}_j)|^2] - S_{\frac{1}{2}}(\tilde{k}_j)|, $$

$$ B := |\mathbb{E}[C_{i_1}^1(k_j)C_{i_2}^2(\tilde{k}_j)]|^2, \quad C := |\mathbb{E}[C_{i_1}^1(k_j)\overline{C_{i_2}^2(\tilde{k}_j)}]|^2, $$

$$ D := \kappa(C_{i_1}^1(k_j), \overline{C_{i_1}^1(k_j)}, C_{i_2}^2(\tilde{k}_j), \overline{C_{i_2}^2(\tilde{k}_j)}). $$

During the proof, we denote by $K$ a finite constant that may change from line to line. We start by bounding the sum concerning $A$. According to Lemmas A.2, F.2 and F.1, we have:

$$ \frac{1}{|I||I_p|} \sum_{(i_1, i_2) \in I \times I_p} A \le \prod_{(\rho, J) \in \{(r, I), (r_p, I_p)\}} \left( \frac{\sqrt{dL}}{2|J|\rho^\beta} \sum_{i \in J} \|\psi_i\|_{\dot{H}^{\frac{\beta}{2}}}^2 + \frac{4}{|J|} \sum_{i \in J} \|\psi_i \mathbb{1}_{\mathbb{R}^d \setminus W/\rho}\|_2 \right) \\ \le K \prod_{(\rho, J) \in \{(r, I), (r_p, I_p)\}} \left( \left( \frac{|J|^{\frac{1}{2d}}}{\rho} \right)^{\beta} + e^{-c'(R/\rho)^{3\theta}} \right). $$

In the second line we used assumption (27). Using again this assumption and the fact that $|J| \le |W|$, we obtain

$$ \frac{1}{|I||I_p|} \sum_{(i_1, i_2) \in I \times I_p} A \le K \prod_{(\rho, J) \in \{(r, I), (r_p, I_p)\}} \left( \left( \frac{|J|}{|W|} \right)^{\frac{\beta}{d}} + \mathcal{E} \right) \le K \left( \frac{(|I||I_p|)^{\frac{\beta}{d}}}{|W|^2} + \mathcal{E} \right). \quad (56) $$
Now, we consider B and C. The idea is to exploit the exponential decorrelation between linear statistics associated with tapers evaluated at distinct frequencies when the number of tapers is small. When this condition does not hold, we instead rely on the variance reduction resulting from the orthogonality of the tapers. Assume that $\rho c_2|k_0| \ge \sqrt{2}|J|^{\frac{1}{2d}}$ for $(\rho, J) \in \{(r, I), (r_p, I_p)\}$. According to Lemma D.4, we have

$$B^{\frac{1}{2}} \le K \left( \sum_{s=1}^{2} \| \rho_s^{-d/2} \psi_{i_s}(\cdot/\rho_s) \mathbf{1}_{\mathbb{R}^d \setminus [-R, R]^d} \|_2 + \| \rho_s^{-d/2} \mathcal{F}[\psi_{i_s}(\cdot/\rho_s)] \mathbf{1}_{|\cdot| \ge |k_j - \tilde{k}_j|/2} \|_2 \right) \\ = K \left( \sum_{s=1}^{2} \| \psi_{i_s} \mathbf{1}_{\mathbb{R}^d \setminus [-R/\rho_s, R/\rho_s]^d} \|_2 + \| \mathcal{F}[\psi_{i_s}] \mathbf{1}_{|\cdot| \ge \rho_s |k_j - \tilde{k}_j|/2} \|_2 \right),$$

where $\rho_1 = r$ and $\rho_2 = r_p$. The $k_0$-allowed assumption ensures that $|k_j - \tilde{k}_j| \ge c_2|k_0|$. Using Lemma F.1, we get, for some $c' > 0$,

$$\frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} B \le K \sum_{\substack{(\rho, J) \in \\ \{(r, I), (r_p, I_p)\}}} e^{-c'(R/\rho)^{3\theta}} + e^{-c'(c_2\rho|k_0|)^{3\theta}}.$$

Moreover, using assumption (27), we obtain $(\rho|k_0|)^2 \ge \sqrt{2}|J|^{\frac{1}{2d}}\rho|k_0| \ge \sqrt{2}|W|^{\frac{1}{d}}|k_0|$. This yields,

$$\frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} B \le K \left( \frac{1}{|W||k|^d} + \mathcal{E} \right), \qquad (57)$$

when $\rho c_2|k_0| \ge \sqrt{2}|J|^{\frac{1}{2d}}$ for $(\rho, J) \in \{(r, I), (r_p, I_p)\}$

Assume now that $\rho c_2|k_0| \le \sqrt{2}|J|^{\frac{1}{2d}}$ for $(\rho, J) = (r, I)$ or $(\rho, J) = (r_p, I_p)$. According to Lemma D.3 and Lemma F.1, we have

$$\frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} B \le \frac{K}{|I| \vee |I_p|} + \frac{K}{|I|} \sum_{i_1 \in I} \|r^{-\frac{d}{2}} \psi_{i_1} \left(\frac{\cdot}{r}\right) \mathbf{1}_{\mathbb{R}^d \setminus W} \|_2^2 + \frac{K}{|I_p|} \sum_{i_2 \in I_p} \|r_p^{-\frac{d}{2}} \psi_{i_2} \left(\frac{\cdot}{r_p}\right) \mathbf{1}_{\mathbb{R}^d \setminus W} \|_2^2 \\ \le K \left( \frac{1}{|I| \vee |I_p|} + \mathcal{E} \right).$$

If $r|k_0|c_2 \le \sqrt{2}|I|^{\frac{1}{2d}}$, we have, using (27) $\sqrt{2}|I|^{\frac{1}{2d}} \ge c_2|k_0||W|^{\frac{1}{d}}/|I|^{\frac{1}{2d}}$. Subsequently, $|I| \vee |I_p| \ge |I| \ge (c_2/\sqrt{2})^d |W||k_0|^d$. If $r_p|k_0|c_2 \le \sqrt{2}|I_p|^{\frac{1}{2d}}$, we have the same lower bound for $|I| \vee |I_p|$. Accordingly, the bound (57) is always valid.

Since we also have $|k_j + k_j| \ge c_2|k_0|$, the technique for bounding B can be applied to obtain

$$\frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} C \le K \left( \frac{1}{|W||k_0|^d} + \mathcal{E} \right). \qquad (58)$$

It remains to consider D. According to Lemma D.2

$$|D| \le K \|r^{-\frac{d}{2}} \psi_{i_1} \left(\frac{\cdot}{r}\right) \|_4^2 \|r_p^{-\frac{d}{2}} \psi_{i_2} \left(\frac{\cdot}{r_p}\right) \|_4^2 = K (rr_p)^{-\frac{d}{2}} \|\psi_{i_1}\|_4^2 \|\psi_{i_2}\|_4^2.$$
Using Lemma F.3 and assumption (27), we obtain

$$ \frac{1}{|I||I_p|} \sum_{(i_1, i_2) \in I \times I_p} |D| \le K \frac{(\log(|I|^{\frac{1}{d}}) \log(|I_p|^{\frac{1}{d}}))^{\frac{d}{2}}}{r^{\frac{d}{2}} |I|^{1/4} r_p^{\frac{d}{2}} |I_p|^{1/4}} \le K \frac{(\log(|W|)^d}{|W|}. \quad (59) $$

Gathering the bounds (56), (57), (58) and (59), we obtain the result. □

To upper bound the variance of $\Delta_I$, we proceed in three steps with Lemmas C.3, C.4 and C.5.

**Lemma C.3.** Suppose the assumptions of Theorem 6.1 hold and denote

$$ C_1 := \kappa(\hat{S}_I(k_{j_1}), \hat{S}_I(k_{j_2}))\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_0(\tilde{k}_{j_2})), $$

$$ C_2 := \kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2})). $$

Then, we have, for $s \in \{1, 2\}$

$$ \frac{1}{N^2} \sum_{\substack{1 \le j_1 \le N \\ 1 \le j_2 \le N}} |C_s| \le K \left( \frac{1}{|I||I_p|} \left( \left( \frac{(|I_p| \vee |I|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}} |k_0|} \right)^{\zeta} + \frac{1}{N} \right) + \frac{\log(|W|)^{2d}}{|W|^2} + \mathcal{E}(|I|, |I_p|) \right), \quad (60) $$

where $K < \infty$ is a constant and $\zeta = d - 1$.

*Proof.* During the proof, $K$ denotes a constant that may change from line to line. We start by proving (60) for $s = 1$. By bi-linearity, and using Corollary G.1, we obtain:

$$ \kappa(\hat{S}_I(k_{j_1}), \hat{S}_I(k_{j_2})) = \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} \kappa(C_{i_1}^1(k_{j_1}), \overline{C_{i_1}^1(k_{j_1})}, C_{i_2}^1(k_{j_2}), \overline{C_{i_2}^1(k_{j_2})}) \\ + |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^1(k_{j_2}))|^2 + |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^1(-k_{j_2}))|^2, \quad (61) $$

where for $s \in \{1, 2\}$, $C_i^s$ denotes the quantity of equation (9) computed with $\Phi_s$. Using Lemma D.2, Lemma F.3 and assumption (27), we get

$$ \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} |\kappa(C_{i_1}^1(k_{j_1}), \overline{C_{i_1}^1(k_{j_1})}, C_{i_2}^1(k_{j_2}), \overline{C_{i_2}^1(k_{j_2})})| \le K \frac{\log(|I|^{\frac{1}{d}})^d}{r^d |I|^{\frac{1}{2}}} \le K \frac{\log(|W|)^d}{|W|}. \quad (62) $$

To bound the other terms, we will take benefit of the decorrelation between distinct frequencies. We start with the one involving $|\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^1(k_{j_2}))|^2$. Assume that $|k_{j_1} - k_{j_2}|r \ge \sqrt{2}|I|^{\frac{1}{2d}}$. According to Lemma D.4, Lemma F.1 and condition (27), we have

$$ \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^1(k_{j_2}))|^2 \le K \mathcal{E}, $$

On the opposite, if $|k_{j_1} - k_{j_2}|r \le \sqrt{2}|I|^{\frac{1}{2d}}$, then using Lemma D.3 and Lemma F.1, we obtain

$$ \frac{1}{|I|^2} \sum_{(i_1, i_2) \in I^2} |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^1(k_{j_2}))|^2 \le K \left( \frac{1}{|I|} + \mathcal{E} \right). $$
We obtain the same upper bound for the term of (61) involving $|\kappa(C_{i_1}^1(k_{j_1}), C_{i_1}^1(-k_{j_2}))|^2$. Accordingly, we have

$$|\kappa(\hat{S}_I(k_{j_1}), \hat{S}_I(k_{j_2}))| \le K \left( \frac{\log(|W|)^d}{|W|} + \left( \mathbf{1}_{|k_{j_1}-k_{j_2}|r \le \sqrt{2}|I|^{\frac{1}{2d}}} + \mathbf{1}_{|k_{j_1}+k_{j_2}|r \le \sqrt{2}|I|^{\frac{1}{2d}}} \right) \frac{1}{|I|} + \mathcal{E} \right).$$

The same line of reasoning also applies to bound $\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))$, yielding

$$|\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))| \le K \left( \frac{\log(|W|)^d}{|W|} + \left( \mathbf{1}_{|\tilde{k}_{j_1}-\tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}} + \mathbf{1}_{|\tilde{k}_{j_1}+\tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}} \right) \frac{1}{|I_p|} + \mathcal{E} \right).$$

Subsequently, we get

$$\frac{1}{N^2} \sum_{1 \le j_1, j_2 \le N} |\kappa(\hat{S}_I(k_{j_1}), \hat{S}_I(k_{j_2}))\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))| \le K \left( \frac{\log(|W|)^{2d}}{|W|^2} + \delta + \mathcal{E} \right),$$

with

$$\begin{aligned}
\delta &= \frac{\log(|W|)^d}{|W|} \frac{1}{|I|N^2} \sum_{1 \le j_1, j_2 \le N} \left( \mathbf{1}_{|k_{j_1}-k_{j_2}|r \le \sqrt{2}|I|^{\frac{1}{2d}}} + \mathbf{1}_{|k_{j_1}+k_{j_2}|r \le \sqrt{2}|I|^{\frac{1}{2d}}} \right) \\
&\quad + \frac{\log(|W|)^d}{|W|} \frac{1}{|I_p|N^2} \sum_{1 \le j_1, j_2 \le N} \left( \mathbf{1}_{|\tilde{k}_{j_1}-\tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}} + \mathbf{1}_{|\tilde{k}_{j_1}+\tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}} \right) \\
&\quad + \frac{1}{|I||I_p|} \frac{1}{N^2} \sum_{1 \le j_1, j_2 \le N} \left( \mathbf{1}_{|k_{j_1}-k_{j_2}| \le \sqrt{2} \frac{|I|^{\frac{1}{2d}}}{r}} \wedge \frac{|I_p|^{\frac{1}{2d}}}{r_p} + \mathbf{1}_{|k_{j_1}+k_{j_2}|r_p \le \sqrt{2} \frac{|I|^{\frac{1}{2d}}}{r}} \wedge \frac{|I_p|^{\frac{1}{2d}}}{r_p} \right) \\
&=: \delta_1 + \delta_2 + \delta_3.
\end{aligned}$$

With the $k_0$-allowed assumption, we obtain

$$\delta_1 \le \frac{\log(|W|)^d}{|W|} \frac{c_3}{N^2|I|} \left( N^2 \left( \frac{|I|^{\frac{1}{2d}}}{r|k_0|} \right)^{\zeta} + N \right) \le K \frac{\log(|W|)^d}{|W||I|} \left( \left( \frac{|I|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right).$$

Similar computations yield

$$\delta_2 \le K \frac{\log(|W|)^d}{|W||I_p|} \left( \left( \frac{|I_p|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right), \quad \delta_3 \le K \frac{1}{|I||I_p|} \left( \left( \frac{(|I_p| \wedge |I|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right).$$

Using the fact that $|I| \vee |I_p| \le c_5|W|^{2\beta/(2\beta+d)}$, for upper bounding $\delta_1$ and $\delta_2$, we obtain

$$\delta \le \frac{K}{|I||I_p|} \left( \left( \frac{(|I_p| \vee |I|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right).$$

This concludes the proof of (60) for $s=1$. The proof of (60) for $s=2$ follows the same lines, but the details differs. By bi-linearity, and using Corollary G.1, we obtain:

$$\begin{aligned}
\kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_2})) = \frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} \kappa(C_{i_1}^1(k_{j_1}), \overline{C_{i_1}^1(k_{j_1})}, C_{i_2}^2(\tilde{k}_{j_2}), \overline{C_{i_2}^2(\tilde{k}_{j_2})}) \\
+ |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^2(\tilde{k}_{j_2}))|^2 + |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^2(-\tilde{k}_{j_2}))|^2.
\end{aligned} \quad (63)$$
Using Lemma D.2, Lemma F.3 and assumption (27), we obtain, as in (62),

$$ \frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} |\kappa(C_{i_1}^1(k_{j_1}), \overline{C_{i_1}^1(k_{j_1})}, C_{i_2}^2(\tilde{k}_{j_2}), \overline{C_{i_2}^2(\tilde{k}_{j_2})})| \le K \frac{\log(|W|)^d}{|W|}. $$

Assume that $|k_{j_1} - \tilde{k}_{j_2}| \ge \sqrt{2}(|I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p)$. Using Lemma D.4 and Lemma F.1 we have

$$ \frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^2(\tilde{k}_{j_2}))|^2 \le K\varepsilon. \quad (64) $$

Suppose that $|k_{j_1} - \tilde{k}_{j_2}| \le \sqrt{2}(|I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p)$. Using Lemma D.3 and Lemma F.1, we obtain

$$ \frac{1}{|I||I_p|} \sum_{\substack{i_1 \in I \\ i_2 \in I_p}} |\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^2(\tilde{k}_{j_2}))|^2 \le K \left( \frac{1}{|I| \lor |I_p|} + \varepsilon \right). \quad (65) $$

We obtain a similar bound for the term involving $|\kappa(C_{i_1}^1(k_{j_1}), C_{i_2}^2(\tilde{k}_{j_2}))|^2$, yielding

$$ \begin{aligned} |\kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))| &\le K \left( \frac{\log(|W|)^d}{|W|} + \varepsilon \right) \\ &\quad + K \left( \mathbf{1}_{|k_{j_1} - \tilde{k}_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} + \mathbf{1}_{|k_{j_1} + \tilde{k}_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} \right) \frac{1}{|I| \lor |I_p|}. \end{aligned} $$

By exchanging $j_1$ and $j_2$, we have

$$ \begin{aligned} |\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}))| &\le K \left( \frac{\log(|W|)^d}{|W|} + \varepsilon \right) \\ &\quad + K \left( \mathbf{1}_{|\tilde{k}_{j_1} - k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} + \mathbf{1}_{|\tilde{k}_{j_1} + k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} \right) \frac{1}{|I| \lor |I_p|}. \end{aligned} $$

The previous inequalities gives

$$ \begin{aligned} \frac{1}{N^2} \sum_{1 \le j_1, j_2 \le N} |\kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_2}))\kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}))| &\le K \left( \left( \frac{\log(|W|)^d}{|W|} + \varepsilon \right)^2 + \delta \right) \\ &\le K \left( \frac{\log(|W|)^{2d}}{|W|^2} + \varepsilon + \delta \right), \end{aligned} $$

where

$$ \begin{aligned} \delta \le \frac{\log(|W|)^d}{|W|N^2|I| \lor |I_p|} \left( \begin{aligned} & \sum_{1 \le j_1, j_2 \le N} \mathbf{1}_{|\tilde{k}_{j_1} - k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} + \mathbf{1}_{|\tilde{k}_{j_1} + k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} \\ & + \mathbf{1}_{|\tilde{k}_{j_2} - k_{j_1}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} + \mathbf{1}_{|\tilde{k}_{j_2} + k_{j_1}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} \end{aligned} \right) \\ + \frac{1}{N^2(|I| \lor |I_p|)^2} \sum_{1 \le j_1, j_2 \le N} \mathbf{1}_{|\tilde{k}_{j_1} - k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} + \mathbf{1}_{|\tilde{k}_{j_1} + k_{j_2}| \le \sqrt{2} |I|^{\frac{1}{2d}}/r \lor |I_p|^{\frac{1}{2d}}/r_p} \end{aligned} $$
Using the $k_0$-allowed assumption, condition (27) and $|I| \le c_5|W|^{2\beta/(2\beta+d)}$, we obtain

$$ 
\begin{aligned}
|\delta| &\le K \left( \frac{\log(|W|)^d}{|W||I| \vee |I_p|} + \frac{1}{(|I| \vee |I_p|)^2} \right) \left( \left( \frac{(|I| \vee |I_p|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right) \\
&\le \frac{K}{|I||I_p|} \left( \left( \frac{(|I_p| \vee |I|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right).
\end{aligned}
$$

This ends the proof for the case $s = 2$.

<figure><img src="image_11.png" /></figure>

**Lemma C.4.** Suppose the assumptions of Theorem 6.1 hold. Denote

$$ B_1 := \mathbb{E}[\hat{S}_I(k_{j_1}) - S_{\frac{1}{2}}(k_{j_1})] \kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}), \hat{S}_0(\tilde{k}_{j_2})), $$

$$ B_2 := \mathbb{E}[\hat{S}_0(\tilde{k}_{j_1}) - S_{\frac{1}{2}}(\tilde{k}_{j_1})] \kappa(\hat{S}_I(k_{j_1}), \hat{S}_I(k_{j_2}), \hat{S}_0(\tilde{k}_{j_2})), $$

$$ B_3 := \mathbb{E}[\hat{S}_I(k_{j_2}) - S_{\frac{1}{2}}(k_{j_2})] \kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_1}), \hat{S}_0(\tilde{k}_{j_2})), $$

$$ B_4 = \mathbb{E}[\hat{S}_0(\tilde{k}_{j_2}) - S_{\frac{1}{2}}(\tilde{k}_{j_2})] \kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2})). $$

Then, for all $s \in \{1, 4\}$,

$$ \frac{1}{N^2} \sum_{1 \le j_1, j_2 \le N} |B_s| \le K (b(I, I_p) + b(I_p, I) + b + \mathcal{E}(|I|, |I_p|)), \quad (66) $$

for a constant $K \in (0, \infty)$ and where, for $\zeta = d - 1$,

$$ b(J_1, J_2) := \left(\frac{|J_1|}{|W|}\right)^{\frac{\beta}{d}} \left( \frac{1}{|J_1||J_2|^{\frac{1}{2}}} \left( \left( \frac{|J_2|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right) + \frac{\log(|W|)^d}{|W||J_1|} \right), $$

$$ b := \frac{\log(|W|)^d}{|W|^2|k_0|^d} + \frac{\log(|W|)^d}{|W|^{\frac{5}{3}}} + \frac{1}{|k_0|^{2d}|W|^2}. $$

*Proof.* We only need to prove the upper bound for $|B_1|$ and $|B_2|$, since by exchanging $j_1$ and $j_2$ we have

$$ \sum_{1 \le j_1, j_2 \le N} |B_3| = \sum_{1 \le j_1, j_2 \le N} |B_1|, \quad \sum_{1 \le j_1, j_2 \le N} |B_4| = \sum_{1 \le j_1, j_2 \le N} |B_2|. $$

We start with $B_1$. Using Lemmas A.2 and F.2 together with assumption (27) we have

$$ |\mathbb{E}[\hat{S}_I(k_{j_1}) - S_{\frac{1}{2}}(k_{j_1})]| \le K \frac{|I|^{\frac{\beta}{2d}}}{r^{\beta}} \le K \left( \frac{|I|}{|W|} \right)^{\frac{\beta}{d}}. $$

For $s \in \{1, 2\}$, we denote by $C_i^s$ the quantity of equation (9) computed with $\Phi_s$. By multi-linearity, we have

$$ \kappa(\hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}), \hat{S}_0(\tilde{k}_{j_2})) = \frac{1}{|I| |I_p|^2} \sum_{i_2 \in I, (i_1, i_3) \in I_p} \kappa(C_{i_1}^2(\tilde{k}_{j_1}) \overline{C_{i_1}^2(\tilde{k}_{j_1})}, C_{i_2}^1(k_{j_2}) \overline{C_{i_2}^1(k_{j_2})}, C_{i_3}^2(\tilde{k}_{j_2}) \overline{C_{i_3}^2(\tilde{k}_{j_2})}). $$
We denote

$$A_1 := C_{i_1}^2(\tilde{k}_{j_1}), \quad A_2 := \overline{C_{i_1}^2(\tilde{k}_{j_1})}, \quad A_3 := C_{i_2}^1(k_{j_2}), \quad A_4 := \overline{C_{i_2}^1(k_{j_2})}, \quad A_5 := C_{i_3}^2(\tilde{k}_{j_2}), \quad A_6 := \overline{C_{i_3}^2(\tilde{k}_{j_2})}.$$

Using Corollary G.1, we expand $\kappa(A_1A_2, A_3A_4, A_5A_6)$:

$$\kappa(A_1A_2, A_3A_4, A_5A_6) = \sum_{\substack{\sigma \in \Pi[6], \\ \sigma \lor \tau_6 = \overline{1}_6}} \prod_{b \in \sigma} \kappa(A_s; s \in b),$$

where $\tau_6 = 12|34|56$ and $\overline{1}_6 = 123456$. Since the random variables considered are centered, the sum runs only over the partitions having blocks of size at least 2. In the following, we use the notation

$$\kappa_{\pi} = \prod_{b \in \pi} \kappa(A_s; s \in b).$$

In the sum over $\sigma$, 5 categories of partitions appear, denoted by **I**, **II**, **III**, **IV** and **V**. To fix ideas, the partitions involved are as follows (refer to [54], Table 1):

**I:** 123456

**II:** 13|2456, 14|2456, 15|2346, 16|2345, 23|1456, 24|1356, 25|1346, 26|1345, 35|1246, 36|1245, 45|1236, 46|1235

**III:** 123|456, 124|356, 125|346, 126|345, 234|156, 134|256

**IV:** 135|246, 235|146, 145|236, 136|245

**V:** 13|25|46, 13|26|45, 14|25|36, 14|26|35, 15|23|46, 15|24|36, 16|23|45, 16|24|35

We bound the cumulant of the partition of category **I**. According to Lemma D.2,

$$|\kappa_{123456}| \le K \|r_p^{-\frac{d}{2}} \psi_{i_1} \left(\frac{\cdot}{r_p}\right)\|_6^2 \|r^{-\frac{d}{2}} \psi_{i_2} \left(\frac{\cdot}{r}\right)\|_6^2 \|r_p^{-\frac{d}{2}} \psi_{i_3} \left(\frac{\cdot}{r_p}\right)\|_6^2 = K r_p^{-\frac{4d}{3}} r^{-\frac{2d}{3}} \|\psi_{i_1}\|_6^2 \|\psi_{i_2}\|_6^2 \|\psi_{i_3}\|_6^2.$$

Then, using Lemma F.3, we get

$$\begin{aligned} \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{123456}| &\le K \frac{\log(|W|^d)}{r_p^{\frac{4d}{3}} |I_p|^{\frac{4}{9}} r^{\frac{2d}{3}} |I|^{\frac{2}{9}}} = K \frac{\log(|W|^d) |I_p|^{\frac{2}{9}} |I|^{\frac{1}{9}}}{r_p^{\frac{4d}{3}} |I_p|^{\frac{2}{3}} r^{\frac{2d}{3}} |I|^{\frac{1}{3}}} \\ &\le K \frac{\log(|W|^d) |W|^{1/3}}{|W|^2} = K \frac{\log(|W|^d)}{|W|^{\frac{5}{3}}}, \end{aligned}$$

where in the last line we used assumption (27).

For the partitions of class **II**, two types of terms appear, depending on either the block of size 2 contains a taper of index $i_2 \in I$ or not. In the first case, the corresponding partitions are given by

$$\mathbf{II}_a := \{13|2456, 14|2456, 23|1456, 24|1356, 35|1246, 36|1245, 45|1236, 46|1235\}.$$
We bound the term 35|1246 and the computations are similar for the others. By Lemma D.2, we have

$$ |\kappa_{1246}| \le K \|r_p^{-\frac{d}{2}} \psi_{i_1} \left(\frac{\cdot}{r_p}\right) \|_4^2 \|r^{-\frac{d}{2}} \psi_{i_2} \left(\frac{\cdot}{r}\right) \|_4 \|r_p^{-\frac{d}{2}} \psi_{i_3} \left(\frac{\cdot}{r_p}\right) \|_4. \quad (67) $$

By Cauchy-Schwarz inequality, we have

$$ \begin{aligned} & \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{35|1246}| \\ & \le \frac{1}{|I_p| r_p^{\frac{d}{2}}} \sum_{i_1 \in I_p} \|\psi_{i_1}\|_4^2 \sqrt{\frac{1}{|I||I_p|} \sum_{\substack{i_2 \in I \\ i_3 \in I_p}} \frac{\|\psi_{i_2}\|_4^2 \|\psi_{i_3}\|_4^2}{(rr_p)^{\frac{d}{2}}}} \sqrt{\frac{1}{|I||I_p|} \sum_{\substack{i_2 \in I \\ i_3 \in I_p}} |\kappa(C_{i_2}^1(k_{j_2})), C_{i_3}^2(\tilde{k}_{j_2})|^2} \\ & \le \frac{K \log(|W|)^d}{|W|} \sqrt{\frac{1}{|I||I_p|} \sum_{\substack{i_2 \in I \\ i_3 \in I_p}} |\kappa(C_{i_2}^1(k_{j_2})), C_{i_3}^2(\tilde{k}_{j_2})|^2}, \end{aligned} $$

where in the last line, we used Lemma F.3 and assumption (27). Using the fact that the family $(r^{-\frac{d}{2}} \psi_i(\cdot/r) e^{-ik_{j_2} \cdot})_{i \in \mathbb{N}^d}$ is orthogonal, we have by Lemma D.1

$$ \frac{1}{|I||I_p|} \sum_{\substack{i_2 \in I \\ i_3 \in I_p}} |\kappa(C_{i_2}^1(k_{j_2})), C_{i_3}^2(\tilde{k}_{j_2})|^2 \le K \left( \frac{1}{|I|} + \varepsilon \right). $$

Accordingly,

$$ \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{35|1246}| \le K \left( \frac{\log(|W|)^d}{|W||I|} + \varepsilon \right). $$

By performing identical computations for the other partitions of category $\mathbf{II}_a$, we obtain

$$ \sum_{\sigma \in \Pi_a} |\kappa_{\sigma}| \le K \left( \frac{\log(|W|)^d}{|W||I|} + \varepsilon \right). $$

We now consider the case where the block of size 2 does not contain tapers of index $i_2$, i.e. the partitions

$$ \mathbf{II}_b := \{15|2346, 16|2345, 25|1346, 26|1345\}. $$

We detail the bound of the term 15|2346. Using Lemma D.2, Cauchy-Schwarz inequality and
Lemma F.3, we have

$$
\begin{aligned}
2_b &:= \frac{1}{|I| ||I_p||^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{15|2346}| \\
&\le K \frac{1}{|I| r^{\frac{d}{2}}} \sum_{i_2 \in I} \|\psi_{i_2}\|_4^2 \sqrt{\frac{1}{|I_p|^2} \sum_{\substack{i_1 \in I_p \\ i_3 \in I_p}} \frac{\|\psi_{i_2}\|_4^2 \|\psi_{i_3}\|_4^2}{r_p^d}} \sqrt{\frac{1}{|I_p|^2} \sum_{\substack{i_1 \in I_p \\ i_3 \in I_p}} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))|^2} \\
&\le \frac{K \log(|W|)^d}{|W|} \sqrt{\frac{1}{|I_p|^2} \sum_{\substack{i_1 \in I_p \\ i_3 \in I_p}} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))|^2}.
\end{aligned}
$$

Using (64) and (65), we get

$$
\frac{1}{|I_p|^2} \sum_{\substack{i_1 \in I_p \\ i_3 \in I_p}} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))|^2 \le \frac{K}{|I_p|} \mathbf{1}_{|\tilde{k}_{j_1} - \tilde{k}_{j_2}|_{rp} \le \sqrt{2}|I_p|^{\frac{1}{2d}}} + K\mathcal{E}.
$$

This yields, using the $k_0$-allowed assumption:

$$
\begin{aligned}
\frac{1}{N^2} \sum_{1 \le j_1, j_2 \le N} |2b| &\le \frac{K \log(|W|)^d}{|W||I_p|^{\frac{1}{2}}} \frac{N}{N^2} \left( N \left( \frac{|I_p|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + 1 \right) + K\mathcal{E} \\
&\le K \left( \frac{1}{|I||I_p|^{\frac{1}{2}}} \left( \left( \frac{|I_p|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right) \right).
\end{aligned}
$$

Now, we consider the partitions of categories III and IV such that $i_2$ appears once in each blocks of size 3, i.e.

$$
\{123|456, 124|356, 126|345, 135|246, 235|146, 145|236, 136|245\}.
$$

We detail the computations for the partition 123|456, that can be directly used for bounding the
other terms. Using the Cauchy-Schwarz inequality, we have

$$
\sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{123|456}| \le \sum_{(i_1, i_3) \in I_p} \sqrt{\sum_{i_2 \in I} \left| \kappa(C_{i_1}^2(\tilde{k}_{j_1}), \overline{C_{i_1}^2(\tilde{k}_{j_1})}, C_{i_2}^1(k_{j_2})) \right|^2} \sqrt{\sum_{i_2 \in I} \left| \kappa(\overline{C_{i_2}^1(k_{j_2})}, C_{i_3}^2(\tilde{k}_{j_2}), \overline{C_{i_3}^2(\tilde{k}_{j_2})}) \right|^2}.
$$

Since $(r^{-\frac{d}{2}} \psi_i(\cdot/r) e^{-ik_{j2}})_{i \in \mathbb{N}^d}$ is orthogonal, we can apply Lemma D.6 to get

$$
\sum_{i_2 \in I} \left| \kappa(C_{i_1}^2(\tilde{k}_{j_1}), \overline{C_{i_1}^2(\tilde{k}_{j_1})}, C_{i_2}^1(k_{j_2})) \right|^2 \leq K \|r_p^{-\frac{d}{2}} \psi_{i_2} \left( \frac{\cdot}{r_p} \right) \|_4^4 + K \delta,
$$
with

$$ \delta = |I| \|r_p^{-\frac{d}{2}} \psi_{i_2} \left(\frac{\cdot}{r_p}\right) 1_{\mathbb{R}^d \setminus W}\|_3 + \sum_{i \in I} \|r^{-\frac{d}{2}} \psi_i \left(\frac{\cdot}{r}\right) 1_{\mathbb{R}^d \setminus W}\|_3. $$

Using Lemma F.1 and assumption (27), we get, after computations, $|\delta| \le K|I|\mathcal{E}$. Subsequently,

$$ \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{123|456}| \le K \sum_{(i_1, i_3) \in I_p} \left( \frac{\|\psi_{i_1}\|_4^2 \|\psi_{i_3}\|_4^2}{r_p^d} + \mathcal{E} \left( \frac{\|\psi_{i_1}\|_4^2}{r_p^{\frac{d}{2}}} + \frac{\|\psi_{i_3}\|_4^2}{r_p^{\frac{d}{2}}} + \mathcal{E} \right) \right). $$

Using Lemma F.3, we obtain

$$ \frac{1}{|I| |I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{123|456}| \le K \left( \frac{\log(|W|^d)}{|W||I|} + \mathcal{E} \right). $$

For the partitions of category III such that $i_2 \in I$ appears in only one block of size 3, i.e.

$\{125|346, 126|345, 234|156, 134|256\},$

we use the decorrelation between frequencies in the cumulants of size 3 where only the indexes belonging to $I_p$ are present. We detail the computations on 125|346 and they can be used to bound the other terms. Assume that $|I_p| \le \sqrt{2}|W||k_0|^d$. According to Lemma D.5, we have

$$ \begin{aligned} |\kappa_{125}| & \le K \sup_{i \in I_p} \left( \|r_p^{-\frac{d}{2}} \psi_i \left(\frac{\cdot}{r_p}\right) 1_{\mathbb{R}^d \setminus W}\|_3 r_p^{-\frac{d}{6}} + \|r_p^{-\frac{d}{2}} \mathcal{F}[\psi_i \left(\frac{\cdot}{r_p}\right)] 1_{|\cdot| \ge \overline{k_{j_2}}/3} \|_1 \right) \\ & = K \sup_{i \in I_p} \left( \|\psi_i 1_{\mathbb{R}^d \setminus \overline{r_p}}\|_3 r_p^{-\frac{d}{3}} + \|\psi_i 1_{|\cdot| \ge c_1 r_p |k_0|/3} \|_1 \right), \end{aligned} $$

where we used the fact that $\psi_i$ is bounded by 1 and of unit $L^2(\mathbb{R}^d)$ norm in the first line. Since $|I_p| \le \sqrt{2}|W||k_0|^d$, we can apply Lemma F.1 to get

$$ |\kappa_{125}| \le K \left( e^{-c_6 |k_0||W|^{\frac{1}{d}}} + \mathcal{E} \right), $$

for some constant $c_6 > 0$. Using Lemma D.5 again, we prove that $|\kappa(A_3, A_4, A_6)|$ is bounded by a constant. Hence,

$$ \frac{1}{|I| |I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{125|346}| \le K \left( e^{-c_6 |k_0||W|^{\frac{1}{d}}} + \mathcal{E} \right) \le K \left( \frac{1}{|k_0|^{2d} |W|^2} + \mathcal{E} \right). $$

On the opposite, assume that $|I_p| > \sqrt{2}|W||k_0|^d$. The Cauchy-Schwarz inequality yields

$$ \begin{aligned} & \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{125|346}| \\ & \le \sum_{\substack{i_2 \in I \\ i_1 \in I_p}} \sqrt{\sum_{i_3 \in I_p} |\kappa(C_{i_1}^1(\widetilde{k_{j_1}}), \overline{C_{i_1}^2(\widetilde{k_{j_1}})}, C_{i_3}^2(\widetilde{k_{j_2}}))|^2} \sqrt{\sum_{i_3 \in I_p} |\kappa(C_{i_1}^1(k_{j_2}), \overline{C_{i_1}^1(k_{j_2})}, C_{i_3}^2(\widetilde{k_{j_2}}))|^2}. \end{aligned} $$
Using Lemma D.6 and Lemma F.2, we obtain after computations

$$ \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{125|346}| \le K \left( \frac{\log(|W|)^d}{|W||I_p|} + \mathcal{E} \right) \le K \left( \frac{\log(|W|)^d}{|W|^2|k_0|^d} + \mathcal{E} \right), $$

since $|I_p| > \sqrt{2}|W||k_0|^d$. Subsequently, we have proved that in both cases

$$ \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{125|346}| \le K \left( \frac{\log(|W|)^d}{|W|^2|k_0|^d} + \frac{1}{|k_0|^{2d}|W|^2} + \mathcal{E} \right), $$

In the study of $B_1$, it remains to consider the partitions of category V. All these terms are bounded by the following steps, that we details on 13|25|46. Assume that $|\tilde{k}_{j_1} - \tilde{k}_{j_2}|r_p \ge \sqrt{2}|I_p|^{\frac{1}{2d}}$. Then, using Lemmas D.4 and F.1, we get $|\kappa(A_2, A_5)| \le K\mathcal{E}$. This leads to

$$ \frac{1}{|I||I_p|^2} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{13|25|46}| \le K\mathcal{E}, $$

when $|\tilde{k}_{j_1} - \tilde{k}_{j_2}|r_p \ge \sqrt{2}|I_p|^{\frac{1}{2d}}$. Assume now that $|\tilde{k}_{j_1} - \tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}$. Using the Cauchy-Schwarz inequality we obtain

$$ \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{13|25|46}| \le \sum_{\substack{i_1, i_2 \in I_p \\ (i_1, i_3) \in I_p}} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))| \sqrt{\sum_{i_2 \in I} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}, C_{i_2}^1(k_{j_2})))|^2} \sqrt{\sum_{i_2 \in I} |\kappa(C_{i_3}^2(\tilde{k}_{j_2}, C_{i_2}^1(k_{j_2})))|^2}. $$

Then, using twice Lemma D.3, we obtain

$$ \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{13|25|46}| \le K \sum_{\substack{i_1, i_2 \in I_p}} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))| \\ \le K|I_p| \sqrt{\sum_{i_1, i_2 \in I_p} |\kappa(C_{i_1}^2(\tilde{k}_{j_1}), C_{i_3}^2(\tilde{k}_{j_2}))|^2} \le K|I_p|^{\frac{3}{2}}. $$

Accordingly, using the $k_0$-allowed assumption, we get

$$ \frac{1}{|I||I_p|^2 N^2} \sum_{1 \le j_1, j_2 \le N} \sum_{\substack{i_2 \in I \\ (i_1, i_3) \in I_p}} |\kappa_{13|25|46}| \le \frac{K}{N^2 |I||I_p|^{\frac{1}{2}}} \sum_{1 \le j_1, j_2 \le N} 1_{|\tilde{k}_{j_1} - \tilde{k}_{j_2}|r_p \le \sqrt{2}|I_p|^{\frac{1}{2d}}} \\ \le K \left( \frac{1}{|I||I_p|^{\frac{1}{2}}} \left( \left( \left( \frac{|I_p|^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right) \right) \right). $$

This concludes the study of $B_1$. We obtain the corresponding bound for $B_2$ by exchanging $I_p$ and $I$, and the bounds for $B_3$ and $B_4$ follow by exchanging $j_1$ and $j_2$. $\square$
**Lemma C.5.** Suppose the assumptions of Theorem 6.1 hold. Then,

$$ \frac{1}{N^2} \sum_{\substack{1 \le j_1 \le N \\ 1 \le j_2 \le N}} |\kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}), \hat{S}_0(\tilde{k}_{j_2}))| \le K \left( \frac{1}{|I||I_p|} \left( \left( \frac{(|I| \vee |I_p|)^{\frac{1}{d}}}{|W|^{\frac{1}{d}}|k_0|} \right)^{\zeta} + \frac{1}{N} \right) + \frac{\log(|W|)^d}{|W|^{\frac{5}{3}}} + \frac{1}{|k_0|^{2d}|W|^2} + \frac{\log(|W|)^d}{|W|^2|k_0|^d} \right) + \varepsilon, $$

for a constant $K \in (0, \infty)$ and $\zeta = d - 1$.

*Proof.* We proceed as in the previous lemma. By multi-linearity and using Corollary G.1, we have

$$ \kappa(\hat{S}_I(k_{j_1}), \hat{S}_0(\tilde{k}_{j_1}), \hat{S}_I(k_{j_2}), \hat{S}_0(\tilde{k}_{j_2})) = \frac{1}{|I|^2|I_p|^2} \sum_{(i_1, i_3) \in I^2} \sum_{\sigma \in \Pi[8], (i_2, i_4) \in I_p^2} \prod_{b \in \sigma} \kappa(A_s; s \in b), $$

where

$$ \begin{aligned} A_1 &:= C_{i_1}^1(k_{j_1}), & A_2 &:= C_{i_1}^1(-k_{j_1}), & A_3 &:= C_{i_2}^2(\tilde{k}_{j_1}), & A_4 &:= C_{i_2}^2(-\tilde{k}_{j_1}), \\ A_5 &:= C_{i_3}^1(k_{j_2}), & A_6 &:= C_{i_3}^1(-k_{j_2}), & A_7 &:= C_{i_4}^2(\tilde{k}_{j_2}), & A_8 &:= C_{i_4}^2(-\tilde{k}_{j_2}). \end{aligned} $$

In the sum over the partitions, there are 12 categories that appear. To fix the idea, one representative of each of them is given by (refer to [54], Table 1):

I : 12345678 (1), II : 12345768 (24), III : 12345678 (24), IV : 13578246 (32),
V : 12354678 (24), VI : 13572468 (8), VII : 12354768 (96), VIII : 13572468 (48)
IX : 12345768 (96), X : 12456728 (48), XI : 13725846 (96), XII : 13254768 (48).

and the number of partitions inside each category is written between the parenthesis. In the following, we use the notation

$$ \kappa_C := \frac{1}{|I|^2|I_p|^2} \sum_{(i_1, i_3) \in I^2} \sum_{\sigma \in C} \prod_{b \in \sigma} |\kappa(A_s; s \in b)|, $$

where C is a category of partitions. We start by bounding the partitions containing a block of size larger than 6. For the partitions of category I, we have, according to Lemma D.2,

$$ \begin{aligned} |\kappa(A_i; i \in [8])| &\le K \|r^{-\frac{d}{2}} \psi_{i_1}(\cdot)\|_8^2 \|r_p^{-\frac{d}{2}} \psi_{i_2}(\cdot)\|_8^2 \|r_p^{-\frac{3d}{2}} \psi_{i_3}(\cdot)\|_8^2 \|r_p^{-\frac{d}{2}} \psi_{i_4}(\cdot)\|_8^2 \\ &= r^{-\frac{3d}{2}} \|\psi_{i_1}\|_8^2 \|\psi_{i_3}\|_8^2 r_p^{-\frac{3d}{2}} \|\psi_{i_2}\|_8^2 \|\psi_{i_4}\|_8^2. \end{aligned} $$

Using Lemma F.3 and assumption (27), we obtain

$$ \kappa_I \le K \frac{\log(|W|)^d}{r^{\frac{3d}{2}} |I|^{\frac{5}{12}} r_p^{\frac{3d}{2}} |I_p|^{\frac{5}{12}}} \le \frac{K \log(|W|)^d (|I||I_p|)^{\frac{1}{3}}}{|W|^3} \le \frac{K \log(|W|)^d}{|W|^{\frac{7}{3}}}. $$