\section{Introduction}
On one hand, point processes are a very important special form of stochastic processes. They have been the focus of numerous books, and are applied to model patterns in criminology, ecology and earth science, just to list a few concrete examples. On the other hand spectral/harmonic analysis of stochastic processes have led to making key insights into the nature of spatial and temporal processes. Despite the long development of spectral analysis in time series, hearkening back to the 1890's, their usage is far from standard in the spatial point process literature, in fact the lack of methodological development has virtually made application non-existent. This paper aims to develop the fundamentals of spectral analysis methodology for stationary point processes. Many of the basic choices for time series and random fields are not well developed for point processes, which we answer head on.

If one only observed a stochastic process once over a given domain, then additional assumptions are required to enable estimation, basically to enable averaging within rather than across ensembles. The standard assumption that enables averaging for estimation, is that of spatial homogeneity or stationarity. We shall also make this assumption throughout this article. Without this assumption any spectral decomposition needs additional care in interpretation, and as many questions remain outstanding even in this case, we confine our interest to stationary processes. In addition, it is possible to assume isotropy, so that dependence becomes characterised in one variable. We shall discuss isotropy, but not confine our interest to solely this special case.

The spectral decomposition of a random process on one hand gives a decomposition
of variance of a stochastic process into a scale--by--scale contribution of energy; on the other hand, it gives the representation of the process as a weighted average of complex exponentials so that given behaviour can be associated with a given physical scale (the inverse of its wavenumber or magnitude wavenumber). The latter notion of scale is key to many further developments in the representation of a particular stochastic process. For point processes the corresponding results are harder to derive as our observations do not take Euclidean values at fixed regular locations but are either absent or present at that given (but random) location. This means that the random location enters non-linearly in the Fourier transform, as the argument of the complex exponential, rather than as in a linear combination.

Spectral analysis of regularly sampled processes has already been studied extensively in the 1-D case~\citep{percival1993spectral}, but even as we study results into 2-D understanding the corresponding developments become more challenging, especially as non-parametric methods can become hard to interpret. For random fields~\citep{stein2012interpolation,Adler,rosenblatt2012gaussian,deb2017asymptotic,guillaumin2019efficient} spectral analysis does resemble the analysis of time series, but the additional challenges of embedded simplified structure (lower dimensional structure) in 2-D and higher D are harder to address, and the appropriate visualization remains an outstanding challenge, see for example~\citep{errico1985spectra,durran2017practical,osgood2019lectures}.
Also, there is some study of random fields sampled by 2-D point processes, as for example by~\citep{chou1992spectral}, in 1-D~\citep{lii1994spectral}.

Despite the potential promise of spectral analysis of point processes, after the initial flurry of development~\citep{bartlett1963spectral,bartlett1964spectral,mugglestone1996practical}, also summarized in~\cite[Ch.~8]{daley2003introduction}, some results followed for its practical utility ~\citep{mugglestone1990spectral,mugglestone1996practical,mugglestone1998detection}, but its promise was never realised into broad usage. Somehow spectral analysis was never converted into a well--used practical tool set for data analysis, and normally it is referred to only when modelling and analysis depend on it \citep{lavancier_determinantal_2015, Moller2007}. We suspect the reason why is, at least partly, that a number of practical choices necessary for constructing spectral summaries have never been substantiated by the proponents of spectral analysis from a theoretical or methodological point-of-view, and therefore remain somewhat mysterious and unclear.  Naturally in the limit of low wavenumbers spectral analysis of point processes must coincide what with spectral analysis of random fields to be self--consistent, but despite this fact, even at low wavenumbers, there are many unexpected features as we shall show. We shall assume that the estimation of the spectrum is the end point of the analysis of the point process.

We can simply contrast this state-of-the-art with the results available for random fields~\citep{rosenblatt2012gaussian,Adler,stein2012interpolation}.  The usage of spectral methods for point processes is limited in parts just because the Cooley--Tukey FFT algorithm~\citep{cooley1965algorithm} cannot be used when performing analysis when by definition, the spacing of observations is not regularly spaced.

The aim of this paper is to understand the properties of any direct spectral estimator defined from a spatial point process, assuming we observe a spatially stationary process. The spectral density captures all relevant second order features of the process, and so this is a fundamental characterization of the process in question. Taking inspiration from standard practice in signal processing, we shall focus on the class of bilinear spectral estimators~\citep{cohen1995time} that form a quadratic
either in the Direct Fourier Transform (DFT) of the data or the data directly, or a linearly modified version, bias-correcting or such as using tapering~\citep{walden2000unified}.

The well-known periodogram non-parametric spectral estimator is naturally a special (and important) case of a bilinear spectral estimator.
Two important tools from time series, of \underline{de-trending}~\citep{percival1993spectral,velasco2000whittle} and \underline{tapering}~\citep{thomson1982spectrum} will be transferred
to point processes, and they will have a unique form in this context. It goes without saying, that the point process setting is much more involved than that of time series, and we introduce appropriate separable and radial tapers, drawing on our expertise from signal processing.

To understand the spectral analysis of point processes therefore we will address the following key points:
\begin{enumerate}
\item How do we define spectral summaries for a point process? We recall existing frameworks for spectral analysis of a point process in Section~\ref{not}.
\item We discuss in  Section~\ref{spectralsum} how direct spectral summaries~\citep{percival1993spectral} are defined, both for general point processes and isotropic processes, drawing upon understanding from spectral analysis for random fields. This includes understanding how to taper a spatial point process.
\item To define method of moments estimators, we answer the question: What are the moments of the DFT?  This is a more complex question than that for random fields or time series,
    and is treated exhaustively in Section~\ref{spectralmom}.
\item To understand linear estimation
    we need both to determine the bias and variance/covariance; Section~\ref{spectralmom2} determines the second order properties of direct spectral estimators.
    \item Armed with the bias and variance, we can address non-parametric estimation.
    This question is treated in Section~\ref{sec:linear} and Section~\ref{sec:sims}.
    \item How does tapering of a point process work? This is further discussed in Section~\ref{sec:linear}.
    \item How should we compute the periodogram (the spectral ``energy'' of a process at a given wavenumber), and does it matter if we assume the process is isotropic or anisotropic?
    This question is further treated in Section~\ref{sec:isotropy}.
    \item How do we sample the spectrum in wavenumber?
    Both to achieve near uncorrelatedness, but also determining when the quality of the estimator becomes so poor that the periodogram contains no more information.
\end{enumerate}

All these points will be clarified in order to enable spectral estimation, and determine the mean square error of a spectral estimator.

Note that the Direct Fourier Transform (DFT) of the point process is \underline{not} zero--mean under the assumption of stationarity, unlike for random fields, as has been noted by~\citep{diggle1987nonparametric}, and discussed more abstractly in terms of random measures by~\citep{daley2003introduction}. We refer to work by~\citep{biscio2019general} and deduce Gaussianity of the DFT.
 We show that the DFT of a point process will not be
zero--mean, no matter its marginal distribution (compare with discussion in ~\citealt{diggle1987nonparametric} for small wavenumbers).  To be concrete, we shall determine the first and second moments of the DFT, which, combined with the Gaussian limiting results, will completely specify the distribution (under growing--domain asymptotics, compare with~\citet{stein1995fixed} for fixed-domain asymptotics).

The Fourier Transform (FT) will be (modulus) squared to yield the periodogram, which will have asymptotically a non-central chi-square distribution, but after de-trending, this will be a chi-square just like the modulus square of the DFT of a random field. This shows that we cannot simply assume that the properties of DFT of a random field will be directly inherited by a FT of a point process, and shows the challenges in using Diggle's~\citep{diggle1987nonparametric} isotropic estimator, as the problem of its sampling in wavenumber is much harder to treat.

Tapering~\citep{walden2000unified} is the second method that we shall explore more fulsomely. It is not clear how to taper in the point process context, or what taper to use, or indeed how to use them. Tapers are normally defined for regularly spaced data. We shall use more rarely applied continuous--space tapers, and show for what classes of processes, multi--tapering produces reduced variance estimators. We note that other choices of implementing tapering could have been made; and this will be discussed, such as tapers for irregularly spaced stochastic processes.
We shall finish up with an extensive simulation study to illustrate the performance of our proposed estimators.

Therefore, we aim to settle many basic outstanding issues in spectral analysis of point processes. We also connect to related work on sampling a continuously indexed random field, see~\citep{chou1992spectral} for example (already developed fully for 1-$D$ processes by~\citealt{lii1994spectral}). Our end result is the arrival at an appropriate smoothing method for the raw spectral estimator of a point process, and also suitable 1-D summaries of such a process, whether the process is isotropic or not. Section~\ref{sec:isotropy} discusses isotropic processes and isotropic summaries. This all combines to underpin our understanding of the spectral representation of point processes. Section~\ref{sec:sims} explores the implications of our results via simulations. We conclude the paper in Section~\ref{sec:disc} with a discussion of the derived results.

\section{Notation and Definitions}\label{not}

We assume that we observe a spatially homogeneous point process $X$
with intensity $\lambda=\rho^{(1)}>0$ and with $\rho^{(2)}(z)$ defined as the (second order) product
density of $X$. This means that $\rho^{(2)}(z)$ only has one argument ($z\in {\mathbb{R}}^d$, namely the spatial shift) rather than depending on two local spatial variables (say $x$ and $y$ rather than $z=x-y$). We assume that $X$ is a simple point process on ${\mathbb{R}}^d$, meaning no duplicate points are allowed. We take $d=2,3,\dots$, and so exclude the case of $d=1$, which is relatively well-studied, see references in e.g.~\citep{daley1971weakly,jowett1972prediction,vere1974elementary}. For a (Borel) set $A$ in ${\mathbb{R}}^d$, we define $|A|$ for the volume of $A$ and $N(A)$ for the number of points of $X\subset {\mathbb{R}}^d$ in $A$. We denote by $B$ any subset of $ {\mathbb{R}}^d$
where points are observed (or equivalently our observation domain).

In complete analogue with a random field we shall define the spectral density of $X$ as the Fourier transform of the complete covariance function of $X$. This is the standard approach, and was first proposed by~\cite{bartlett1964spectral}.
The complete covariance function of a stationary point process is therefore
\begin{equation}
\label{covy}
\gamma(z) \equiv \lambda\delta(z) + \rho^{(2)}(z) - \lambda^2,\quad z\in {\mathbb{R}^d},
\end{equation}
where $\delta(\cdot)$ is Dirac's delta function.
The spectral density function (sdf) ~\citep{bartlett1964spectral} of the point process $X$ is then defined as the Fourier transform of the complete covariance function:
\begin{equation}
\label{def:spectrum}
f(\k) \equiv \mathcal{F}[\gamma](\k)\ = \int_{{\mathbb{R}}^d}e^{-2\pi i\k\cdot z}\gamma(z)\,dz\  =\  \lambda + \int_{{\mathbb{R}^d}}e^{-2\pi i\k\cdot z}[\rho^{(2)}(z)-\lambda^2]dz, \qquad \k
\in {\mathbb{R}}^d.
\end{equation}
This representation is in direct analogy with the corresponding spectral decomposition of a random field or a time series.
The symbol $\mathcal{F}$ denotes the Fourier transform, with $i$ as the imaginary unit. We refer to the argument $\k$ as the ``wave number'' rather than frequency to knowledge its multi-dimensionality. This is a natural choice of a Fourier transform in analogy to the analysis of random fields~\citep{stein2012interpolation}. To avoid additional constants in the inverse Fourier transforms,
we parameterise the Fourier transform with wave number instead of the customary angular frequency $\xi = 2\pi \k$ used by \cite{bartlett1964spectral} and some others.

For a time series or random field there are a number of spectral results, ranging from the spectral representation theorem (I)~\cite[pp.130]{Percival1993SAPA} and to Bochner's (Herglotz) theorem (II)~\citep{brockwell1991time,Percival1993SAPA}.
Both these results  are not exactly available in the point process setting, but for a time series $X_t$ we can note them:
\begin{align}
\label{Xt}
X_t&\overset{(I)}{=}\mu+\int_{-\frac{1}{2}}^{\frac{1}{2}} dZ(\omega)e^{2i \pi \omega t},\quad {\mathrm{\bf Var}}\{dZ(\omega)\}\overset{(II)}{=}f(\omega)d\omega
\\
\gamma(t)&\overset{(III)}{=}\int_{-\frac{1}{2}}^{\frac{1}{2}}f(\omega)e^{2\pi i \omega t}\,d\omega,\quad f(\omega)\overset{(IV)}{=}\sum_\tau \gamma(\tau) e^{-2i \pi \omega \cdot \tau}.
\label{gamt}
\end{align}
Equation~\eqref{Xt} decomposes $X_t$ into random contributions associated with each angular frequency $\omega$. The covariance in Eqn~\eqref{gamt} is also decomposed into weighted frequency contribution. Both of these observations are important for the interpretation of the Fourier representation of $X_t$. ``Important'' contributions correspond to $f(\omega)$ being considerably larger relative to other  $f(\omega)$ as in that scenario ${\mathrm{\bf Var}}\{dZ(\omega)\}$ is bigger than ${\mathrm{\bf Var}}\{dZ(\omega')\}$ and therefore $|dZ(\omega)|$ likelier to be larger than $|dZ(\omega')|$.
Spectral representations of point processes are also discussed in~\citep[Ch.~8]{daley2003introduction}. We can yet again decompose the (complete) covariance
in terms of a spectral representation like in (II), and it takes the form in terms of a $d$-dimensional Fourier transform of
\begin{equation}
f(w)=\int e^{-2\pi i w\cdot z}\gamma(z)\;dz=\lambda+\int e^{-2\pi i w\cdot z} \left\{\rho^{(2)}(z)-\lambda^2\right\}\;dz.
\label{weightav}
\end{equation}
We refer to $f(w)$ as the \underline{spectrum} of $X$. We note that the Fourier transform can be inverted to yield the equality of (an analogy of III):
\begin{equation}
    \rho^{(2)}(z)=\lambda^2+\int_{{\mathbb{R}}^d} f(w) e^{2\pi i w\cdot z}\,dz.
    \label{eqn:rhobyscale}
\end{equation}

We assume the spectrum is the primary object of interest in our study of point patterns, as the covariance $\gamma(z)$ can be fully determined from it, and the covariance specifies the point process. We see directly from~\eqref{weightav}
that $\rho^{(2)}(z)-\lambda^2$ and {\em not} the complete covariance $\gamma(z)$ is playing the role of a time series covariance. This is as when transforming $\gamma(z)$ rather than $\rho^{(2)}(z)-\lambda^2$ we otherwise get a constant boost of $\lambda$ that spreads across all frequencies. This does not indicate to us which frequencies are ``more'' important unlike for a random field. By important we here mean contribute to the process more than any other wavenumber.

We have to exercise some caution when interpreting the Fourier transform as a bijective transform. Yes, it should contain the same information about scale, but the meaning of the word scale will be different. In time or space the notion is associated with the support of $\gamma(z)$ or $\rho^{(2)}(z)$. In frequency being supported over a wavenumber $w$ means variation over scales $1/\|w\|$ is present in the correlation function, or to represent the variability in $\gamma(z)$ we need frequencies $\|w\|$. The function $\rho^{(2)}(z)-\lambda^2$ is often approached in terms of what scales it is non-zero at, but the variation in the function can be associated with long or short scales. So for a covariance function, we now have two notions of what it means to possess scales $\|w\|$; either $\gamma(1/\|w\|\mathbf{u})$ is non-zero for unit vector $\mathbf{u}$ or $\gamma(z)$ is variable (changing) and the scale of the change is $\|w\|$. Imagine observing a sinusoid with period $k$. The function will hit unity at a regular set of $z$ values, and so it will be supported at those $z$ values, but also its Fourier transform will be supported at $1/k$.

\begin{figure}[!h]
	\centering
	\includegraphics[width=.9\linewidth]{fig-pcf_sdf_pairs_1x2-poissonEq1.pdf}
	\caption{Pair correlation $pcf = \rho^{(2)}/\lambda^2$ (pcf, left) and corresponding scaled spectral density function $f/\lambda$ (sdf/$\lambda$, right) for two stationary and isotropic non-Poisson point process models, two variants each, with known pcfs. For the Thomas processes the sdf is given in the text. For the regular Matern II process the sdf was numerically approximated using the Hankell transform (cf. Section 7). See details of the processes and the variants in Section \ref{sec:sims}.}
	\label{fig:fig-1}
\end{figure}

This is clear from Figure~\ref{fig:fig-1} where covariance is present at smaller or medium scales smoothly for the clustered Thomas processes, or repelled for the Mat\'{e}rn hard-core processes. As the Mat\'{e}rn process's complete covariance function is discontinuous, its Fourier transform is supported over all scales (due to the Heisenberg-Gabor uncertainty principle~\citep{cohen1995time}). This can be deduced from~\eqref{eqn:rhobyscale}. To reproduce the discontinuity we need all scales.

We see from~\eqref{weightav} that not all frequencies or wavenumbers are given equal weighting. First all wavenumbers are given an equal weighting $\lambda$ and then the Fourier transform of $\rho^{(2)}(z)-\lambda^2$ determines the wavenumbers that are up--weighted (added) or down--weighted (subtracted) relative to the overall level of $\lambda$. We then observe what wavenumbers are important to the representation of the point process, which gives more information, conveniently decomposed on a scale--by--scale manner.

For a random field or a stochastic process in $d$--dimensions a few cartoon characteristics are important. For a discretely sampled process in $\mathbb{Z}^d$ the simplest random process, white noise, is constant across wavenumbers yielding a spectrum  that takes the form of $\sigma^2$ on $[-\frac{\pi}{\Delta},\frac{\pi}{\Delta}]^d$, where $\Delta>0$ is the sampling period, and zero otherwise. For a point process
the choice of definition of the spectral density does not imply a decay because of the inclusion of the term $\lambda \delta(z)$ in space. However once we remove this term, we expect a decay of the remaining spectrum $f(w)-\lambda$ as   $|w|\rightarrow \infty$. Otherwise
~\eqref{weightav} would contain additional singularities.

Furthermore at $w=0$ we retain
\begin{equation}
f(0)=\lambda+\int \left\{\rho^{(2)}(z)-\lambda^2\right\}\;dz.
\label{weightav2}
\end{equation}
The second term in this expansion is not required to be zero, and for the Thomas process, it is not zero. For a Thomas process the spectrum is~\citep[p~55]{mugglestone1996practical}
\begin{equation}
f(\omega)=\lambda_0+\lambda_0'e^{-\frac{1}{2}a^2\omega\cdot\omega}.
\end{equation}
It is clear from this expression that we cannot arrive at a constant contribution due to the exponential.

As we shall be studying the moments of a point process, it is convenient to restate
Campbell's theorem \citep[Sec. 4.3.3]{SKM}, and we shall use this result multiple times. The  theorem applies to any measurable function $h:R^{nd}\mapsto R^+$ with $n=1,2,...$, and states that
\begin{equation}
\label{eqn:Campbell}
{\mathbf{E}} \sum_{x_1,...,x_n\in X}^{\neq}h(x_1,...,x_n)=\int_{{\mathbb{R}}^{nd}}h(x_1,...,x_n)\rho^{(n)}(x_1,...,x_n)dx_1...dx_n,
\end{equation}
where the summation is over distinct point tuples. Note that if the point pattern $X$ is stationary then $ \rho^{(1)}(x)=\lambda$ is constant for any $x$, and $\rho^{(n)}(x_1,...,x_n)=\rho^{(n)}(o, x_2-x_1,...,x_n-x_1) \equiv \bar{\rho}^{(n)}(x_2-x_1,...,x_n-x_1)$.
We shall use the $\bar{\rho}^{(n)}$ notation to explicitly avoid confusion when parameterizing in terms of the shift $z_1=x_2-x_1,\dots, z_{n-1}=x_n-x_1$ (and the stationary parameterization) rather than the basic coordinates.

We define a new parameterisation to capture how the process's $f(w)$ differs from that of a Poisson process. We define the $n$th deviation as
\begin{equation}
\label{eqn:excess}
    \widetilde{\rho}^{(n)}(z)=\frac{\bar{\rho}^{(n)}(z)-\lambda^n}{\lambda^n},\quad n=2,3,\dots .
\end{equation}
For a Poisson process these $n=2,3,\dots$ deviations from Poissonianity
will all be identically zero. By understanding deviations from Poisson behaviour, we get greater understanding of the underlying process under study.
For completeness we also define the Fourier transform of the deviations:
\[ \widetilde{f}^{(n)}(w)={\cal F}\left\{\widetilde{\rho}^{(n)}(z)\right\},\]
where if $n=2$ we suppress the superscript.
With this definition we find that
\begin{equation}
\label{eqn:fourierexcess}
    \widetilde{f}(w)=\frac{f(w)-\lambda}{\lambda^2},\quad w \in{\mathbb{R}}^d,\quad f(w)=\lambda+\lambda^2\widetilde{f}(w).
\end{equation}
Thus in a sense, $\widetilde{f}(w)$ encapsulates the deviation of the function from a constant spectral density of $\lambda$ via the term $\lambda^2\widetilde{f}(w)$.

\citet[p.~337]{daley2003introduction} discusses some differences in spectral representation of time series versus that of point processes. We would argue that decay of the spectrum is still reasonable to assume once $\lambda$ has been subtracted (so the decay of $\widetilde{f}(w)$ is reasonable to assume).  Having established the theoretical spectral description of point processes, we now turn to their sampling properties.

\section{Direct Spectral Summaries of the Data}\label{spectralsum}

We shall re-examine spectral estimation of point processes, and clarify a number of points that are currently unclear.
We start from Bartlett's periodogram estimator based on observing a point pattern in region $B\subset {\mathbb{R}}^d$ which can be written as:
\begin{align}
\label{bartlett2}
I_0(\k) &:=\hat{\lambda}  + |B|^{-1}\sum_{x\neq y,\in X\cap B}e^{-2\pi i\k\cdot (x-y)} ,\quad \k\in  {\mathbb{R}^d},
\end{align}
where we have set $\hat{\lambda} =|B|^{-1}|X\cap B|$. This definition uses all the data available to us, $\{x:\;x\in X\cap B\}$, but is simply {\em a possible choice} amongst all possible direct spectral estimators.  In general a direct spectral estimator is one which is \underline{bilinear} in the DFT of the data, see for example the discussion in~\citep[Ch 5--6]{Percival1993SAPA}. We recall a bilinear form is a function
$D({\mathbf{x}},{\mathbf{y}})$ satisfying $D({\mathbf{x}}_1+{\mathbf{x}}_2,{\mathbf{y}})=D({\mathbf{x}}_1,{\mathbf{y}})+D({\mathbf{x}}_2,{\mathbf{y}})$, and equivalently in the second argument.

Bilinear forms have been thoroughly discussed in signal processing for the analysis of time series
~\citep{loughlin1993bilinear}. As a point process consists of locations the best we can hope for in terms of bilinearity will be in terms of \underline{sesquilinearity} of the Fourier transform as the point locations will appear in the argument of the complex exponential. Sesquilinearity of $D({\mathbf{x}},{\mathbf{y}})$ simply generalizes a bilinear form to the Hermitian symmetry, additionally requiring $D({\mathbf{x}},{\mathbf{y}})=D({\mathbf{y}}^\ast,{\mathbf{x}})$ as well as $D({\mathbf{x}},{\mathbf{y}})$ satisfying $D({\mathbf{x}}_1+{\mathbf{x}}_2,{\mathbf{y}})=D({\mathbf{x}}_1,{\mathbf{y}})+D({\mathbf{x}}_2,{\mathbf{y}})$.

In time series the bilinear form has been chosen to ensure that spectral estimators are real-valued, and often non-negative though some bilinear estimators are not,  see e.g. the Kirkwood-Rihaczek spectrum~\citep{hindberg2006kernel}, and Guyon's spectral estimator~\citep{guyon1982parameter}.
Spectral estimators in time series are bilinear in the observed real-valued  process $X$
and sesquilinear in its Fourier transform $J_h(\k)$. We shall still require that the form of the estimator is sesquilinear in the DFT of the point process, but the estimators will not be bilinear in $X$, the point pattern, as this will not be possible.
We define the tapered DFT of $X$ of a point pattern for  a specified general (square) integrable function $h(x)$ (referred to as a `data taper' by ~\citealt{percival1993spectral}) with Fourier transform $H(\k)$, to be
\begin{equation}
J_h(\k):=\sum_{x \in X\cap B} h(x) e^{-2\pi i\k\cdot x},\quad \k\in {\mathbb{R}}^d,
\label{taper1}
\end{equation}
for any unit norm $h$, e.g. $\|h\|_2=1$.

Because $x$ appears in the argument of the complex exponential, we cannot define estimator bilinear in the data. In practice one has to decide what taper function to use. In 1D for point processes tapering has been used~\citep{cohen2014multi,taleb2018multiresolution} with continuous wavelets, but also multitapering has been used on interpolated data~\citep{das2020multitaper}. Here, we do not interpolate, do not implement localized analysis, and do not implement analysis in 1D. Higher dimension is more complex to interpolate, as it is not orderly.

We note that in general $h(x) \exp{(-2\pi i\k\cdot x)}$ is identically distributed but not independent for many choices of point processes $x$, making the choice of a Central Limit Theorem (CLT) a bit more complex. For a large class of processes one such CLT is given by \cite{biscio2019general}. If we compare the quantity in \eqref{taper1} with~\citep{biscio2019general}, then we see that $q=1$ and $p=2$, in their notation.
We do not have to worry about $e^{-2\pi i\k\cdot x}$ being a function of several of the points, and this is why $q=1$. Requirement ${\cal H}1$ is therefore trivially not a problem as an assumption, ${\cal H}2$ is about the scaling of the observation area, as we assume the sides scale the same not a problem, ${\cal H}3$ (higher order centred powers) are limited through the usage of Campbell's theorem and as $q=1$ we do not have to worry about pathological scaling between components in ${\cal H}4$.
Therefore, citing~\citep{biscio2019general}, we can deduce from their Theorem 1 that as $|B|$ diverges, and recover the second moments by direct calculations in Section~\ref{spectralmom}:
\begin{align}
    \label{CLT}
    J_h(\k)&\overset{L}{\rightarrow}N_C\left(\lambda H(\k),f(w),r(w)\right),\\
r(w)&=\lambda \int_{\mathbb{R}^d}
H(w'-2w) H(w')dw'+\int_{\mathbb{R}^d} U(w,z) e^{-2\pi i w\cdot  z}\{\rho(z)-\lambda^2\}dxdz.
\end{align}
Here $|B|\rightarrow \infty$ in equal speed in all directions of $B$ and $B$ is the cuboid, and we define the tapered related term $H(w)$ in~\eqref{eqn:Hdef}. We deduce the moments from
Proposition 5.1 combined with Theorem 5.1.
We can view $J_h(\k)$ as a complex-valued scalar or a real-valued two-vector. $N_C(\mu,\Sigma,C)$ is the general complex normal~\citep{schreier2003second} and its arguments are its mean $\mu$, it covariance $\Sigma$, as well as its complimentary covariance $C$. It should be contrasted with the complex propert normal $N_C(\mu,\Sigma),$ that has zero relation.

Based on $J_h(\k)$ the spectral estimator becomes its modulus square or
\begin{equation}
I_h(\k) :=\left|J_h(\k)\right|^2=\left| \sum_{x \in X\cap B} h(x) e^{-2\pi i\k\cdot x} \right|^2=\sum_{x\in X\cap B}h(x)h^*(x) + \sum_{x,y\in X\cap B}^{\neq}h(x)h^*(y)e^{-2\pi i\k\cdot (x-y)}.
\label{periodogram2}
\end{equation}
If we take $h(z)=h_0(z)=\frac{1}{\sqrt{|B|}}1(z\in B)$ then we recover Bartlett's peridogram in~\eqref{bartlett2} from~\eqref{periodogram2}.
The point of using a general function $h(z)$
is that abruptly ending the inclusion of points when we leave the region $B$ causes ripples in wavenumber, and therefore a worse estimation of the spectrum, as is also the case of time series when not using a taper~\citep{percival1993spectral}.
We have deduced an asymptotic distribution for $J_h(\k)$
using~\citep{biscio2019general} but this result cannot be applied to $I_h(\k)$ even if it superficially may seem to be of the correct form.

What can we then say about $I_h(\k)$?
Using the continuous mapping theorem~\citep{dasgupta2008asymptotic} we can deduce from~\eqref{CLT} that if we consider only arguments $\k$ so that the complementary covariance is negligible then
\begin{equation}
    J_h(\k)\overset{L}{\rightarrow}N(\lambda H(\k),f(w))\Rightarrow \quad
     |J_h(\k)|^2\overset{L}{\rightarrow}\frac{f(w)}{2}\chi^2_2(\lambda^2 |H(\k)|^2),
\end{equation}
where the parameters of the non--central $\chi^2_k(m^2)$
($k$ denoting the degrees of freedom, and $m^2$ the non-centrality parameter). We give $H(w)$ in~\eqref{eqn:Hdef}.
We can also use the uniform integrability of the variable $ J_h(\k)$, to get the asymptotic moments  of $I_0(\k)$ from these limits.

Note that the sum in~\eqref{periodogram2} cannot be evaluated in a computationally efficient manner (compare to~\citep{cooley1965algorithm}), as the locations $\{x\}$ are \underline{not} regularly spaced, which is unfortunate, but unavoidable. Finally~\eqref{periodogram2} is  bilinear in
$J_h(w)$ but it is not bilinear in $X$, unlike the commensurate expressions for the DFT of time series and random fields.

For completeness we also note the isotropic estimator of
\citep[Eqn.~3.3]{diggle1987nonparametric} in 2D, namely
\begin{align}
\label{diggle}
\bar{I}_0(\|\k\|) &:=\hat{\lambda}+\frac{1}{|B|}\sum_{x,y\in X\cap B}^{\neq}
\mathcal{J}_0\left(2\pi\|\k \| \|x-y\|\right),
\end{align}
where $\mathcal{J}_0(x)$ denotes the Bessel function of order 0, specified in Section \ref{sec:isotropy}. The $d$-dimensional extension is given later in equation~\eqref{diggleD}.
$\bar{I}_0$ is less clearly bilinear in the data, but if we start from the estimator of~\eqref{bartlett2} and average it over orientations analytically, then we arrive at this form. We therefore abuse our terminology to also refer to it as a bilinear estimator.
There is one additional correction made by~\citep{diggle1987nonparametric}. As we shall see, for low wave numbers there is a bias inherent in~\eqref{diggle}. To address this problem~\citep[Eqn.~3.4]{diggle1987nonparametric} suggests taking for some choice of lower bound $t_0$
\begin{equation}\label{diggle4}
I_D(\|\k\|)=\left\{ \begin{array}{lcr}
\bar{I}_0(t_0) & {\mathrm{if}} & \|\k\|\leq t_0\\
\bar{I}_0(\|\k\|) & {\mathrm{if}} &  \|\k\| > t_0
\end{array}\right. .
\end{equation}
\cite{diggle1987nonparametric} suggest taking $t_0$ so that $\bar{I}_0(t_0)$ is at a minimum, and also suggest smoothing the estimated
${I}_D$, which we clarify in the linear estimation section. The authors also discuss iterated methods of bias correction in their Section~5, and correcting biased estimators of the correlation. They discuss difficulties of preserving short range features in the spectrum.

Note also that we have modified the estimator in~\eqref{diggle} to
divide by
$|B|$ rather then the observed number of points in the region, $N(B)$, as do for example~\citep{diggle1987nonparametric}, as it is much preferable not to have a random denominator.  We shall discuss the usage of $I_0(\cdot)$ versus $I_D(\cdot)$ further on in Section~\ref{sec:isotropy}, and the implications of when the process is isotropic or not.

As can be surmised from~\eqref{periodogram2} the expectation of the general estimator ${\mathbf{E}}\left\{ I_h(\k)\right\}$ is a convolution between the Fourier transform of the window and the true spectrum, and if $h(z)$ does not go to zero nicely at the boundary of $B$ then the periodogram becomes quite complicated in terms of its expectation.
Looking at~\eqref{periodogram2} also raises a second problem, namely that there is a constant contribution
$|B|^{-1}|X\cap B|$ which does not give any wave number specific information, and is correlated between wave numbers. This term is new to time series and random fields, if not new to the expectation of the periodogram of randomly sampled stochastic processes, see~\citep{lii1994spectral,chou1992spectral}.

We can generalize the family of estimators from Eqn.~\eqref{periodogram2} to take the form of
\begin{align}
\label{Ih}
I_h(\k) &:= \left| \sum_{x \in X\cap B} h(x) e^{-2\pi i\k\cdot x} \right|^2 =|J_h(\k)|^2,\quad \k\in  {\mathbb{R}^d}.
\end{align}
We study the {\em family of estimators} $I_h(\k)$. We also note that~\eqref{Ih} is a direct spectral estimator~\cite[p.~207]{percival1993spectral}. Multidimensional tapers are available to us~\citep{hanssen1997multidimensional,simons2011spatiospectral}, and can be precomputed. Our choice of tapering corresponds to using continuous tapers evaluated at random locations.

Should we use no data taper, or equivalently use the classical taper $h_0$ from (8), then we introduce the notation of writing $I_0(w)$ for the periodogram for clarity.
We then also define
\begin{equation}
\label{eqn:Hdef}
H=\mathcal{F}[h],\quad H_0=\mathcal{F}[h_0].
\end{equation}

The expectation of any member of that family takes the form of:
\begin{equation}
\label{expectation}
{\mathbf{E}}\left\{ I_h(\k) \right\}= \int_{{\mathbb{R}^d}} |H(\k'-\k)|^2 f(\k')d\k' + \lambda^2 |H(\k)|^2.
\end{equation}
This can be derived straightforwardly by direct computation from~\eqref{Ih} by taking expectations with Campbell's formula and using the convolution theorem.

Should we chose to taper isotropically in 2-D-space then instead we take with isotropic taper $h_I(\cdot)$
\begin{equation}
\label{diggle2}
    \bar{I}_{h}(\|\k\|) :=\hat{\lambda}+\sum_{x,y\in X\cap B}^{\neq}
h_I\left(\|x-y\| \right) \mathcal{J}_0\left(2\pi \|\k \| \cdot \|x-y\|\right).
\end{equation}
Note that this cannot necessarily be constructed. Say $J_h()$ with an isotropic taper would be made with the largest spherical domain and a compact support could be constructed, but what does that imply for $I_h()$?
Radial tapers can be determined as described in~\citep{simons2011spatiospectral}, whether in 2-D continuous space or discrete space.
Equation~\eqref{diggle2} can be of course be extended into higher dimensions. In general we would have in dimension $d=1,2,3,4,\dots$
\begin{equation}
\label{diggleD}
    \bar{I}_{h}(\|\k\|) :=\hat{\lambda}+\sum_{x,y\in X\cap B}^{\neq}
h_I\left(\|x-y\| \right)  \mathcal{J}_{d/2-1}\left(\|\k \| \|x-y\|\right) \|x-y\|^{-(d/2-1)}.
\end{equation}
We can see directly that
\begin{align}
\nonumber
    {\mathbf{E}}\{ \bar{I}_{h}(\|\k\|)\}
    &=\lambda+\int_{B}\int_{B\cap B_{-z}}
    h_I\left(\|z\| \right)\ \mathcal{J}_{d/2-1}\left(2\pi\|\k \| \|z\|\right)\ \|z\|^{d/2-1}\ \tilde{\rho}^{(2)}(\|z\|) \ dz dx.
\end{align}
As the change of variable to polar coordinates means we will multiply by $\|z\|^{d-1}$, and this will correspond to a $d/2-1$ Hankel transform of the taper $h(\cdot)$ multiplied by the product density. For more details on the Hankel transform see~\citep{grafakos2013fourier,cressie2015statistics}.

We can also use more than one taper, and will use a sequence of orthogonal functions
$\{h_k\}$~\citep{percival1993spectral}, that will be used to get a new estimator. We will assume they are all of unit energy and are all (pairwise) orthogonal. This means we will require $||h_k||_2=1$ and
\begin{equation}
\int_{{\mathbb{R}}^d} h_k(z) h_l(z)\,dz=\delta_{kl}.
\end{equation}
Tapering has been viewed with variable approbation. Criticism comes especially for throwing away data~\citep{fougere1977solution,yuen1979comments}.
We shall study the properties of direct spectral estimators of the form of Eqn~\eqref{Ih}. This will help us to characterise the second order properties of the process $X$.

Just as a brief note: it is not obvious to know how to taper. Most of time series analysis is based on discrete regular sampling for instance. We chose to use continuous tapers, rather than interpolating the points to a regular grid.
This leaves the spheroidal wavefunctions (continuous but hard to compute), as well as the cosine tapers of~\cite{riedel1995minimum}. These correspond to separable taper choices; non-separable choices with be discussed in Section \ref{sec:isotropy}.

\section{Expectation of Bilinear Spectral Estimators}\label{spectralmom}
We shall start from the simplest spectral estimator, namely the Bartlett periodogram as specified by~\eqref{bartlett2}.
We shall now be quite concrete and understand some common cases of spatial data, and their sampling, and choose as early special cases Cartesian product domains.
We shall focus on the cuboid sampling domain:
\[{B}_\square=[0,\ell_1]\times\dots \times [0,\ell_d].\]

\begin{Lemma}\label{lemma1}
Assume that $X$ is a homogeneous point process with intensity $\lambda$ and twice differentiable spectrum $f(\k)$ at all values of $\k\neq 0$. Assume $X$ is observed in a cuboid domain $B_\square$ with a centroid at 0, which is growing in every dimension. Then the expected value of the periodogram $I_0(w)$ satisfies the equation
\begin{equation}
\label{EIk}
{\mathbf{E}} I_0(\k)= {f}(\k)+ |B|^{-1}\lambda^2 T(B, \k)+o(1),\quad \k\neq 0,
\end{equation}
for
\[T(B,\k)=\prod_{j=1}^d \frac{\sin^2(\pi \k_j\ell_j)}{(\pi \k_j)^2}.\]
\end{Lemma}
\begin{proof}
See Appendix~\ref{sec:Pf1}.
\end{proof}
This is not what we would expect as a large sample expectation of the spectrum, given our experience for time series and random fields in that the term $|B|^{-1}\lambda^2 T(B, \k)$ has been added. To get there, we need to understand the DFT
$J_h(w)$ better.

\begin{Proposition}\label{DFTprop}
Assume that $X$ is a stationary point process with intensity $\lambda$ and with spectrum $f(w)$, and that $h$ is a unit energy taper supported in the domain $B_\square$ itself only. Then
the first moments of the direct spectral estimator $J_h(w)$ are given by:
\begin{align}
\label{EDFT}
{\mathbf{E}}\left\{ J_h(w)\right\}&={\mathbf{E}} \sum_{x\in X\cap B}h(x)\exp(-2\pi i w\cdot x)\\
\nonumber
&=\lambda \int_{\mathbb{R}^d}h(x)\exp(-2\pi i w\cdot  x)\;dx=
\lambda H(w),\quad w\in {\mathbb{R}^d}\\
\nonumber
{\mathrm{\bf Var}}\left\{ J_h(w)\right\}&=
\lambda  + \lambda^2\int_{\mathbb{R}^d}
\left| H(k'-w)\right|^2\widetilde{f}(k')\;dk'+\lambda^2 \left|H(w)\right|^2-\lambda^2 |H(w)|^2\\
\label{varDFT}
&=\lambda  + \lambda^2\int_{\mathbb{R}^d}
\left| H(k'-w)\right|^2\widetilde{f}(k')\;dk'\\
{\mathrm{Rel}}\left\{ J_h(w)\right\}&=
\lambda \int_{\mathbb{R}^d}
H(k'-2w) H(k')dk'+\int_{\mathbb{R}^d} U(w,z) e^{-2\pi i w\cdot z}\{\rho(z)-\lambda^2\}dxdz,
\end{align}
with $U(w,z)=\int h(x)h(z+x)e^{-2\pi i w\cdot (2x)}\;dx$.
\end{Proposition}
\begin{proof}
See Appendix~\ref{DFTproppf} and notice the definition of
~\eqref{eqn:Hdef}.
\end{proof}
In the above Proposition the term ${\mathrm{Rel}}\left\{\cdot\right\}$ denotes ``relation'', also known as the ``complimentary covariance'', see e.g.~\citep{miller1974complex}.
These moments will have implications for the distributional properties of
$I_h(w)=\left| J_h(w)\right|^2$.
 As we consider $\k=0$ then the variance and the relation are the same. For $w$ small if $J_h(w)$
is Gaussian but the mean is non-zero and the moments are different than from the moments of a zero-mean time series and likely to lead to a non-centred $\chi^2$ distribution.

As the DFT is the sufficient statistic (asymptotically) of a stationary process, in a sense, determining its first two moments are pivotal. Perhaps most important is~\eqref{EDFT} because of the non-central nature of the distribution. Understanding such distrubutions, e.g. the non-central $\chi^2$ of the square of a Rician distribution, is not uncommon in medical imaging~\citep{olhede2011nonparametric}. Inference is often not facilitated subsequently to making this observation--instead we shall propose to remove the mean, removing the non-centrality parameter.

We thus define a bias-corrected periodogram (see~\citealt[p.~292]{daley2003introduction}, where the signed measure is used to make the equivalent definition), or
\begin{equation}
\label{biascorrect}
\widetilde{I}_h(w)=\left| J_h(w)-\lambda H(w)\right|^2=
\left| \widetilde{J}_h(w)\right|^2.
\end{equation}
Thus $J_h(w)-\lambda H(w)$ is often called the ``signed measure''. This is the DFT of the process $X^0$ dual to the mean-corrected random measure $N^0(dx)=N(dx) - \lambda dx$ where $N$ is the dual of $X$. As we have subtracted $H(w)$ off the discrete Fourier transform, it is no longer strictly positive, but once we take the modulus square we are guaranteed to arrive at a real-valued positive quantity.

There is one wavenumber which is problematic, namely $w=0$. For the periodogram we have $h_0(x)=\mathbf{1}_B(x)/\sqrt{|B|}$. Thus we have
\begin{equation}
J_0(w)=\frac{1}{\sqrt{|B|}}\sum_{x\in X\cap B}e^{-2\pi i w\cdot x},\quad
J_0(0)=\frac{1}{\sqrt{|B|}} N(B)=\sqrt{\|B\|}\widehat{\lambda}.
\end{equation}

This removes the problem of a non-zero mean of the DFT, as long as we assume that we know the intensity $\lambda$, and so are in the position to remove this effect. Assuming knowledge of this quantity is not a major hurdle, as it can be estimated consistently for largish areas ($|B|\rightarrow \infty$), by just counting the number of points and dividing by the area.
For completeness we also define the signed measure DFT as
\begin{equation}
\widetilde{J}_h(w)=J_h(w)-\hat{\lambda} H(w),\quad
.
\label{signed}
\end{equation}
We note directly from~\eqref{CLT}, that
\[\widetilde{J}_h(w)\rightarrow N_C\{0,f(w),r(w)\}.\]
Also it follows
\begin{equation}
\widetilde{J}_0(0)=\sqrt{|B|}\widehat{\lambda}-\widehat{\lambda}H_0(0)=\sqrt{|B|}\widehat{\lambda}-\sqrt{|B|}\widehat{\lambda}=0,
\end{equation}
trivially. Thus we cannot estimate the DFT at wave number zero.
For a time series analysis when calculating DFTs we subtract off the sample mean, and then also the mean corrected DFT is zero at wave number zero. In a time series setting the periodogram at frequency zero is often not plotted.

We note that two ways of bias subtraction for the periodogram is possible, namely $\widetilde{I}(w)$ as described by~\eqref{biascorrect}
and additionally we could have defined
\begin{equation}
\label{eqn:u2}
\breve{I}(w;h)= \sum_{x\in X\cap B}\sum_{y\in X\cap B} h(x)h(y)\exp\left\{-2\pi i w\cdot (x-y)\right\}-\widehat{\lambda^2}\left|H(w)\right|^2.
\end{equation}
The advantage of the first estimator is that it is naturally non-negative and removes the error before we take the modulus square. This also should have a positive impact on the variance. In addition,  $\widehat{\lambda}^2$ can be a biased estimator for $\lambda^2$.

If we are looking at a more sophisticated estimator than the separable estimator of~\eqref{eqn:u2} or~\eqref{biascorrect} see e.g.~\citep[Ch.~6]{cohen1995time}, then we have
\begin{equation}
\label{Ig}
    I_g(w)=\sum_{x\in X\cap B} \sum_{y\in X\cap B} g(x,y) e^{-2\pi iw\cdot \left(x-y\right)},
\end{equation}
where $g(x,y)$ is a {\em kernel} which may have a specified support, that has to be selected. For example we could define $g(x,y)=h(x)h^\ast(y)$ with $h(x)$ a multi-dimensional data taper ~\citep{hanssen1997multidimensional}, or we could make a non-separable choice of $g(x,y)=h(\|x-y\|)$. Depending on the choice of
$g(x,y)$ the quantity $I_g(w)$ may satisfy a number of desirable characteristics. For example positivity, asymptotic unbiasedness as an estimator, and  computational efficiency.

The quantity $\widetilde{I}_h(w)$ will be referred to as the ``debiased periodogram'', (and the subscript unity will be removed) and is preferable to use to the normal periodogram as it removes a number of unpalatable effects.

\begin{Lemma}\label{lemdebias}
Assume that $X$ is a homogeneous point process with a spectral density $f(w)$. Then the bias-corrected tapered periodogram
has a first moment given by:
\begin{align}
\nonumber
{\mathbf{E}} \left\{ \widetilde{I}_h(w)\right\}
&={\mathbf{E}} I_h(\k)-2\lambda\Re\left\{J_h(w)H^\ast(w) \right\}
+\lambda^2 \left|H(w)\right|^2\\
&=\int_{{\mathbb{R}^d}} |H(\k'-\k)|^2 f(\k')d\k'.
\label{Ih2}
\end{align}
\end{Lemma}

\begin{proof}
See Appendix~\ref{pflemdebias}.
\end{proof}
Thus as along as $H(w)$ is getting more concentrated in wave number (e.g. $H(w)\rightarrow \delta(w)$), this is asymptotically in $|B|$ an unbiased estimator of $f(w)$.
Using the signed measure DFT of~\eqref{signed} to study $X$ is more convenient, as this lets us avoid the contribution that affects the low wave numbers.

We now look to determine further properties of $I_h(\k)$, or properties of
$\widetilde{I}_h(\k)$. To produce an estimator that is smoothed we need to study the covariance and variance of $\widetilde{I}_h(\k_1)$ and
$\widetilde{I}_h(\k_2)$.
We first look at ${\mathbf{E}} \left\{ \widetilde{I}_h(\k)\right\}$. We find that
its form for large spatial regions is specified by the following theorem.

Let us write down what grid we get large sample unbiasesness for, customarily called the \emph{Fourier grid}, and this is useful when the observational domain $B$ is a box.
\begin{Definition}
The Fourier wave number grid for a point process observed on a cuboid domain $B_\square(\bm{l})$ corresponds to the points
\begin{equation}
{\cal K}(\bm{l}):={\cal K}(B_\square(\bm{l}))=\{\k_n,\;
\k_n=\begin{pmatrix} p_{n1}/l_1,\dots, p_{nd}/l_d\end{pmatrix},\quad p_{nj}\in{\mathbb{Z}}\}.
\end{equation}
\end{Definition}
Note that the physical units of the wave number elements is per unit length, such as $m^{-1}$. Referring to Lemma~\ref{lemma1} we see that using this grid both removes any asymptotic bias at low wave numbers, but also gives us a sampling of the wave numbers that approximately leads to independent periodogram ordinates, rather like in the random field case. Note that the centering of the box $B_\square$ at zero, e.g. requiring the box to take the form
\[B_\square=\begin{pmatrix} -\frac{\ell_1}{2} & \frac{\ell_1}{2}\end{pmatrix}\times \dots \times
\begin{pmatrix} -\frac{\ell_d}{2} & \frac{\ell_d}{2}\end{pmatrix},
\]
is required for the stated result. If the observational domain is not centred, then additional phase factors enter the expression. We can compare these results to wave number grid we chose to evaluate a standard spectral estimator of a time series, the \underline{Fourier frequencies}~\citep[p.~197--198]{percival1993spectral}. Their basic importance follows because we can expect the direct Fourier transform to be Gaussian and so uncorrelated implies independence.

A second feature of time series is the notion of the Nyquist wave number~\citep[p.~98]{percival1993spectral}. This does not exist for point processes.
It may seem counter intuitive that there is no upper limit to the wave numbers we can estimate. When analysing a process that has been sampled in space, such as a random field, we expect to see aliasing. Aliasing is when variation that is happening very rapidly is confused with slower cycles, because when we have sampled more sparsely, rapid variation cannot be resolved. Because the distances can be any real-valued value, the Nyquist wave number does not exist.

\begin{Theorem}{Large--Domain Expectation of the tapered periodogram.}\label{thm1}
Assume that $X$ is a stationary point process with intensity $\lambda$ and twice differentiable spectrum $f(\k)$ at all values of $\k\neq 0$, and that $h$ is a unit energy taper (e.g. $\|h\|=1$) supported in the cuboid domain $B_\square$ only. Assume $X$ is observed in $B_\square$, which is growing in every dimension, that is $\min \ell_j\rightarrow \infty$. Then the expected value of the tapered periodogram $I_h(w)$ satisfies the equation
\begin{equation}
\label{Prop2expr}
{\mathbf{E}} I_h(w)= f(w)+\lambda^2|H(w)|^2+o(1),\quad w\neq 0,
\end{equation}
where $H(w)$ is the Fourier transform of $h(x)$, as defined by Eqn.~\eqref{eqn:Hdef}.
\end{Theorem}
\begin{proof}
See Appendix~\ref{sec:Pf2}.
\end{proof}
In fact, asymptotics are well--understood both for time series and random fields. Our understanding of asymptotics for spectral estimation of point processes falls back to~\cite{daley2003introduction} or \cite{diggle1987nonparametric}. In the latter reference \cite{diggle1987nonparametric} refers back to Chemistry theory by Hansen and McDonald, discussing the properties of gasses.

\begin{Corollary}\label{prop1}
Assume that $X$ is a stationary point process with intensity $\lambda$ and twice differentiable spectrum $f(\k)$ at all values of $\k\neq 0$, and that $h$ is a unit energy taper supported in the domain $B$ only. Assume $B$ is growing in every dimension. Then the expected value of the tapered periodogram $I_h(\k)$ satisfies the equation
\begin{equation}
\label{Prop2expr3}
{\mathbf{E}} I_h(\k)= f(\k)+\lambda^2|H(w)|^2+o(1),\quad \k\neq 0,
\end{equation}
where $H(\k)$ is the Fourier transform of $h(x)$.
\end{Corollary}

\begin{proof}
The proof can be found in Appendix~\ref{pfthm2}.
\end{proof}
This result is not unexpected and generalizes Theorem~\ref{thm1}. The main difference here is that we consider a general domain $B$ rather than a cuboid domain $B_\square$.
This theorem allows us to determine the expectation of $\widetilde{I}_h(w)$, in general, as is made clear in the following corollary.

\begin{Corollary}{Bias--corrected periodogram.}\label{cor1}
Assume that $X$ is a stationary point process with intensity $\lambda$ and twice differentiable spectrum $f(\k)$ at all values of $\k\neq 0$, and that $h$ is a unit energy taper supported in the cuboid domain $B_\square$ only. Assume $X$ is observed in $B_\square$, which is growing in every dimension. Then the expected value of the tapered periodogram $\widetilde{I}_h(\k)$ satisfies the equation
\begin{equation}
\label{Prop2expr2}
{\mathbf{E}} \widetilde{I}_h(\k)= f(\k)+o(1),\quad \k\neq 0.
\end{equation}
\end{Corollary}
\begin{proof}
See Appendix~\ref{proofcor1}.
\end{proof}

This corollary shows the importance of removing the spectral bias before squaring--otherwise it gets hard to isolate and remove the effect at zero wave number--which is exacerbated at higher intensities (growing $\lambda$), as is clear from~\eqref{Prop2expr}.
So the estimator $\widetilde{I}_h(\k)$ is debiased relative to ${\mathbf{E}}\left\{  I_h(w)\right\}$ (as we shall see from Theorem~\ref{prop1}), but is still guaranteed to be non-negative. This can be compared to removing $\widehat{\mu}$ from a time series before calculating the periodogram.

This gives us insight into the act of debiasing. In a sense our understanding of the spectral nature of point processes is informed by the spectral representation of randomly sampled processes~\citep{chou1992spectral}. We can continue to explore
the differences and similarities between spatial point processes and random fields.

Finally, to enable smoothing, define $P\ge 1$ estimates of the spectrum via
\begin{equation}
I_p^t(\k) =\sum_{x,y\in X \cap B}h_p(x) h^\ast_p(y)e^{-i2\pi \k\cdot (x-y)},\quad p=1,\dots P,
\label{tapered}
\end{equation}
where we assume
\[\int_{{\mathbb{R}^d}} h_p(x) h_l^\ast(x)dx = \delta_{pl}\]
where $\delta_{pl}=1$ if $p=l$ and $0$ otherwise.
These estimates can also be arrived at starting from~\eqref{taper1} and $J_h(\k)$. As we shall use $P$ tapers linearly, and $P^2$ for the quadratic estimators, we shall use the subscript `0' to refer to the untapered periodogram, and $p=1,\dots,P$ to refer to the subsequent tapers.

Tapering is very common for stochastic processes. Initially the idea was hotly contested, but its utility is now firmly established both for time series and random fields~\citep{dahlhaus1987edge,walden2000unified}. The idea of tapering is to ameliorate edge effects that leads to the asymptotically leading bias term. This is important in 1-D~\citep{thomson1982spectrum,walden2000unified} but of greater importance in 2-D and higher~\citep{dahlhaus1987edge}. In 2-D and higher the edge effects become more pronounced, and indeed asymptotically dominant~\citep{kent1996spectral,mardia1984maximum}. This paper also uses multitapering to stabilize the variance of any estimator, as described for time series in~\citep{walden2000unified}. We shall describe the necessary steps to perform linear estimation of the spectrum of a 2-D and higher dimension point process.

\begin{figure}[!ht]
	\centering
	\includegraphics[width=1\linewidth]{fig-4x3-example-data-periodograms-nonsquare-truncated.pdf}
	\caption{\emph{Top row}: Four example data patterns in a rectangle of area 2e4 units, exhibiting regularity, complete randomness, small scale clustering, and clustering with a left-right axis directionality. \emph{Second row}: Their debiased periodograms excluding $\k=0$. \emph{Third row}: Their multitapered debiased periodograms, $3\times 3$ sine-tapers. \emph{Bottom row}: Their rotationally averaged periodograms, including the case without debiasing. The 2D summaries use a $\k$-grid with a step size $0.005$ in both directions. The rotational averaging was estimated at $|w|=0.0050, 0.0075, ..., 0.2500$ using a box kernel having a radius 1.5 times the $\k$-grid step size 0.005.}
	\label{fig:dataexample}
\end{figure}

Let us set theory aside for a moment and have a look at some spectral estimates: Figure \ref{fig:dataexample} demonstrates the debiased periodogram, the debiased multitapered periodogram with sine-tapers and a rotationally averaged 1D summary of the peridogram (Sections \ref{sec:linear}\&\ref{sec:isotropy}) for four point patterns exhibiting different structural behaviour. The aspect ratio of the plots are kept the same only to have an orderly figure: Due to debiasing we can estimate the periodogram on a wavenumber grid of our own choosing, and not just on the Fourier grid (which in these finite examples would be different for different tapers, cf. Section \ref{spectralmom2}). The  wavenumber-grid in this illustration is a regular grid with a step size 0.005 in both dimensions.

We see that information is present at wavenumbers up to $|w|\approx 0.2$. The tapering smooths the periodogram, as expected. There is a dark well near the origin for the regular pattern, and a bright hump for the clustered pattern, and both features transfer to the rotationally averaged curves (compare with Figure \ref{fig:fig-1}). The anisotropy of the fourth pattern is hinted by the anisotropy visible in the 2D periodograms as elongation of the "hump" in the second dimension. The elongation is perpendicular to the elongations of the clusters in the data because wavenumbers relate inversely to spatial units. The rotationally averaged periodogram was also computed from the (not shown) non-debiased periodogram for comparison. There is a very prominent bias near $|w|=0$ for the non-debiased version, which hopefully illustrates why the proposed debiasing step is relevant when applying the method in practice on freely chosen wavenumber grids.

\section{Variance and Covariance of Bilinear Spectral Estimators}
\label{spectralmom2}

Before determining how to do bilinear estimation of the spectral density let us determine the second order properties of spectral estimators. Core to our understanding of smoothing will be the variance and covariance of the tapered functions. This will be established in the following theorem. For brevity define
\begin{equation}
    \phi_\k(x)=\exp\{-2\pi i \k\cdot x \},\quad \k,x\in{\mathbb{R}}^d.
\end{equation}
With this notation we have
\begin{equation}
I_p^t(\k) =\sum_{x,y\in X\cap B}h_{p}(x)h_{p}^\ast(y)\phi_\k(x-y),\quad p=1,\dots P.
\label{tapered3}
\end{equation}
Above we note that we have the same value of $p$ across $x$ and $y$ in~\eqref{tapered3} as this corresponds to a modulus square.

\begin{Theorem}\label{Propcovar}
Let $I^t_p(k)$ denote the (tapered) periodogram given by
~\eqref{tapered}. The covariance between the (tapered)  periodogram and itself across two wavenumbers $\k_1,\k_2\in{\mathbb{R}}^d$ then takes the form of
\begin{align*}
{\mathrm{\bf Cov}}& \{I^t_p(\k_1) , I^t_l(\k_2)\}
=\int_{B^4}
 \rho^{(4)}(x,y,z,v)
 h_p(x) h_p(y)h_l(z) h_l(v)\phi_{\k_1}(x-y) \phi_{\k_2}(z-v)\;dxdydvdz\\
 \nonumber
 &+\int_{B^3} \rho^{(3)}(x,y,v)\left[h_p(x) h_p(y)h_l(x) h_l(v)\phi_{\k_1}(x-y)\left\{ \phi_{\k_2}(x-v)+  \phi_{-\k_2}(x-v)\right\}\right.\\
&+h_p(x) h_p(y)h_l(y) h_l(v)\phi_{\k_1}(x-y)\left\{ \phi_{\k_2}(y-v)+  \phi_{-\k_2}(y-v)\right\}\\
&\left.+\phi_{\k_2}(x-y)h_p^2(v)h_l(x)h_l(y)+\phi_{\k_1}(x-y) h_l^2(v)h_p(x)h_p(y)\right]\;
dxdydv\\
&+\int_{B^2}\rho^{(2)}(x,y)\left[ h_p(x) h_p(y)h_l(x) h_l(y)\phi_{\k_1}(x-y)\left\{
\phi_{\k_2}(x-y)+\phi_{-\k_2}(x-y)\right\}\right.\\
&+\left\{h_p^2(x)h_l(x)h_l(y)+ h_l^2(x)h_p(x)h_p(y)\right\}
\left.\left\{\phi_{\k_1}(x-y)+\phi_{\k_2}(x-y)\right\}+h_p^2(x)h_l^2(y)\right]\,dxdy\\
 &+\lambda
 \int_{B} h_p^2(x)h_l^2(x)\,dx-\lambda^2-\lambda \int_{B^2} h_p(x)h_p(y)h_l(x)h_l(y)(\phi_{\k_1}(x-y)+\phi_{\k_2}(x-y))\rho^{(2)}(x,y)dxdy\\
 &-\left\{  \int_{B^2} h_p(x)h_p(y)\phi_{\k_1}(x-y)\rho^{(2)}(x,y)dxdy \right\}\cdot \left\{  \int_{B^2} h_l(x)h_l(y)\phi_{\k_2}(x-y)\rho^{(2)}(x,y)dxdy \right\}.
\end{align*}
\end{Theorem}
\begin{proof}
The proof is in Appendix~\ref{sec:Pf3}.
\end{proof}
Note that we can in fact write:
\[\widetilde{I}_h(w)=I_h(w)-2\lambda\Re\left\{ J_h(w)H^\ast(w)\right\}
+\lambda^2 \left|H(w)\right|^2.\]

The expression derived in Theorem~\ref{Propcovar}  is quite general, and only requires the assumptions of homogeneity in space, but it is hard to understand, especially as it is hard to order the terms in size without putting very restrictive assumptions in place.
We can still consider the untapered case by setting $p=l$ and taking $h(x)=|B|^{-1/2}\mathbf{1}_B(x)$, to arrive at an expression.
Let us study this general correlation in the instance of the Poisson process where we know that
\begin{equation}
\label{assumpPoisson}
\rho^{(n)}(x_1,\dots, x_n)=\lambda^n,
\end{equation}
and consider the case of $\k_1=\k_2=\k$.

For completeness we also define the spectral bandwidth $b_h$ by
\begin{align}
    b^2_h&=\int_{{\mathbb{R}}^d}\|\k\|^2 \left|H(\k) \right|^2d\k.
\end{align}

\begin{Proposition}\label{Propcovar2}
Let $I^t_p(\k)$ denote the tapered periodogram given by
~\eqref{tapered} for a Poisson process with intensity $\lambda$. The covariance between the tapered periodogram with one taper and the tapered periodogram with another taper then takes the form for $\|k\|>b_h$  of
\begin{align*}
{\mathrm{\bf Cov}}& \{I^t_p(\k) , I^t_l(\k)\}=\lambda
 \int_{B} h_p^2(x)h_l^2(x)\,dx+\lambda^2\delta_{pl}+o(1).
\end{align*}
The form for $\|\k\|<b_h$ is much more complex, and is provided in Appendix~\ref{sec:Pf4}.
\end{Proposition}
\begin{proof}
The proof can be found in Appendix~\ref{sec:Pf4}.
\end{proof}
\begin{Remark}
At this stage we have simplified our assumptions to the Poisson process which is quite disappointing. The reason why we have chosen to do so is clear from Theorem~\ref{Propcovar}. If we assume the process is Poisson then a lot of terms cancel. The correct treatment of non-Poisson processes would correspond to transforming the terms into a wavenumber representation, doing a Taylor series because of the concentration of the tapers, and then implementing the same calculation after the Taylor series.
\end{Remark}
\begin{Proposition}\label{Propcovar3}
Let $I^t_p(\k)$ denote the tapered periodogram given by
~\eqref{tapered} for a Poisson process with intensity $\lambda$. The covariance between the (tapered) periodogram and itself at different wave numbers takes the form of
\begin{align*}
{\mathrm{\bf Cov}}& \{I^t_p(\k_1) , I^t_p(\k_2)\}=\lambda \|h_p\|_4^4=o(1).
\end{align*}
\end{Proposition}
\begin{proof}
The proof can be found in Appendix~\ref{sec:Pf5}.
\end{proof}
Fix $0<\epsilon<1$. This is so that we can approximately find an uncorrelated grid. Recall that $h(x)$ is unit energy in space or wave number.

\begin{Definition}
The tapered Fourier wave number grid for a point process observed on a cuboid domain $B_\square(\bm{l})$ corresponds to the points
\begin{equation}
{\cal K}(B_\square(\bm{l}))=\{\k_n,\;
\k_n=\begin{pmatrix} p_{n1}\tau(\epsilon) /l_1,\dots, p_{nd}\tau(\epsilon)/l_d\end{pmatrix},\quad p_{nj}\in{\mathbb{Z}}\},
\end{equation}
where $\tau$ is defined from the equation
\[\min_{\tau(\epsilon)}\left|\int H_p(\tau(\epsilon)-\k') H_l^\ast(\k') \;d\k'\right|<\epsilon.\]
\end{Definition}

These results hint at producing general methodology; we can for the Poisson process figure out how to do non-parametric estimation. This is however not enough. We need to determine the variance of the periodogram in more general. The following proposition establishes its form.

\begin{Theorem}\label{Propvar1}
Assume $X$ is a stationary process. The variance of the (tapered) periodogram then takes the form of
\begin{equation*}
{\mathrm{\bf Var}}\{I_p^t(\k)\}=\lambda^4\left\{  \widetilde{f}^{(4)}({\k},\k,\k)-\widetilde{f}^{(2)}(\k)(1+2\widetilde{f}^{(2)}(\k)\right\}+\lambda^2\widetilde{f}(\k)+\lambda^2
  +\lambda\|h\|_4^4+o(1).
\end{equation*}
\end{Theorem}
\begin{proof}
The proof can be found in Appendix~\ref{sec:Pf6}.
\end{proof}
With these moments we can now consider the problem of estimation.
We might feel an increasing degree of unease as we have succumbed to the usual fallacy of only deriving the most general results for a Poisson process. However as the observational domain becomes larger we would expect $|H_p(w)|$ to concentrate and so the lack of constancy will become less of an issue.

\section{Linear Smoothing}\label{sec:linear}
Having derived the first and second order properties of the periodogram,
the next task is then estimation of $f(\k)$, as the primal quantity of interest in studying the point process $X$. If we fully characterise $f(\k)$ then we have fully characterised $\gamma(z)$, as the Fourier transform is a bijection. Just like for a time series or random field, a key question concerns the correct method for quadratic characterisation and estimation.

For the process itself ($X$), then
its first moment is encapsulated by $\lambda$ for a stationary (or homogeneous) processes. We recall the spectrum uniquely maps to the complete covariance function $\gamma(z)$.  We want to start from the random variables $\{I_h(w)\}$ (or $\{\widetilde{I}_h(w)\}$) and then non-parametrically estimate the spectral density directly from those quantities. As long as we seek to estimate $f(\k)$ at points $\k$ where the function $f(\k)$ is continuous,
averaging $\{I_0(\k)\}$ (or $\{\widetilde{I}_h(w)\}$) locally seems like a sensible strategy.
\citep{diggle1987nonparametric} proposed smoothing $\widetilde{I}_D(\k)$ as defined in~\eqref{diggle4}. They just note that they apply a weighted average smoother with user tuned weights. This will not determine how to sample wavenumbers, or remove correlation, and we also have the alternative form of bias removal to them before smoothing.

We shall only treat growing domain scenarios ($|B|\rightarrow \infty$
(one can contrast growing domain observations and increasing density of observations), see e.g.~\cite{zhang2005towards} or \cite{stein1995fixed}. Of course when faced with a real point pattern we cannot tell if asymptotic methods offer a good approximation to the finite sample distribution.

For sufficiently high wave numbers the constant effect in the spectral density with expectation $\lambda$ will dominate both in the mean and variance of $I_0(w)$. Thus for higher wave numbers in $w$ the periodogram will eventually become perfectly correlated over wavenumbers, and will be converging to $\widehat{\lambda}$. For a growing domain we can show that
$\hat{\lambda}$ converges in probability to $\lambda$.

Our understanding of the correlation of higher wave-numbers
should be contrasted with the discussion of~\cite{mugglestone1996exploratory}. The latter authors
remark that there is an upper limit on the number of independent periodogram ordinates, see either of ~\cite{mugglestone2001spectral,mugglestone1996practical}, but do not justify that statement (which is based on random field analysis rather than point processes). It is not clear from these manuscripts exactly how this upper limit is arrived at, as~\cite{mugglestone2001spectral}
it is based on heuristic (or one could say practical) arguments on the nominal degrees of freedom present in the data, in the random field limit.

We have already described how to de--bias $J_h(w)$ into
$\widetilde{J}_h(w)$, which is a properly zero--mean quantity.
We clearly do not need to do this as long as we stick to the Fourier grid.
Mean correction has also been discussed by ~\cite[p.332]{daley2003introduction} to achieve analogous results to the spectral representation of time series.
This mean correction has the advantage of letting us determine the properties of $J_h(w)$ directly, and as we chose to evaluate Bartlett's rather than Diggle's spectral estimator, we have the complex--valued quantity available to us, and can do the mean correction while still retaining a real-valued positive spectral estimate. Mean correction for the periodogram after squaring is more precarious, as we are not guaranteed a positive estimate.

Having understood the periodogram point-wise in terms of its expectation, to understand how to smooth it, we also need to make sense of its variance/covariance structure.  Smoothing requires an explicit understanding of local degrees of freedom.
For this purpose we shall use tapered estimates, see also~\citep{walden2000unified}. We refer to the previous section and~\eqref{tapered} for the basic taper definitions.

We have two possible estimators, namely the bias--corrected multi-taper estimator of~\citep{walden2000unified} as well as
a multi-dimensional kernel density estimator~\citep{duong2005cross}. They are both linear in the estimated spectrum. We shall use the periodogram rather than another estimator for the kernel density estimation, as it is easier to keep track of bias and correlation.

We therefore define the multitaper estimator~\citep{walden2000unified} (and refer to~\citep{hanssen1997multidimensional,liu1992multiple})
\begin{equation}
    \overline{I}_P^t(w)=\frac{1}{P}\sum_{p=1}^P
    \widetilde{I}_{h_p}(w),\quad w\in{\mathbb{R}}^d,
\end{equation}
where $ \widetilde{I}_{h_p}(w)$ is defined in~\eqref{biascorrect}. Choosing to use continuous tapers in~\eqref{biascorrect} is not the only option, and this definition allows for nonseparable tapers~\citep{simons2011spatiospectral}, as it does not loop over Cartesian product tapers. We could have calculated tapers for our spatial samples, see~\citep{wingham1992reconstruction,bronez1988spectral}, or we could have interpolated the points~\citep{das2020multitaper}. Neither solution is palatable, and so we have stuck with continuous space tapers evaluated at random locations. This is a deliberate choice, as otherwise our method would become opaque to analysis.

If we use tapers then to avoid bias, it is necessary to debias $\widetilde{I}_{h_p}(w)$. If we stick to the Fourier grid, then the bias is zero and not necessary to remove. Thus starting from multi-dimensional kernel density estimators~\citep{duong2005cross} we instead chose a window function $W(x)$
and  matrix bandwith
$\bm{\Omega}$. We then define
\begin{equation}
\label{smoothingI}
    \overline{I}_h(w;\bm{\Omega})=\sum_{\k'\in {\cal K}} |\bm{\Omega}|^{-1}W\left(\bm{\Omega}^{-1}\left(\k'-\k \right) \right)\widetilde{I}_h\left(\k'\right),
\end{equation}
where $\widetilde{I}_h\left(\k\right)$ is defined by~\eqref{biascorrect}. The `$h$' just denotes the taper, and the tilde captures that we are using a bias-corrected periodogram. A useful modification is to consider one-dimensional window $W$ and $\bm\Omega=\Omega$ and consider the averaging over directions, thus reducing the statistic to magnitude $t=\|\k\|$ only,
\begin{equation}
\label{angularaveragedI}
    \overline{I}_h(t;\Omega)=\sum_{\k'\in {\cal K}} {\Omega}^{-1}W\left({\Omega}^{-1}\left(||\k'||- t \right) \right)\widetilde{I}_h\left(\k'\right),
\end{equation}
which we refer to as the rotationally averaged periodogram. \cite{mugglestone1996practical} call a version of it the $R$-spectrum. Alternatively, one could average also over the magnitudes to arrive at a summary function in direction. In 2D the convenient argument would be the polar angle; \cite{mugglestone1996practical} call this version the $\theta$-spectrum.

Basically both the multi-tapering and the local averaging smooth any direct spectral estimate (or bilinear estimator). Both of them have a smoothing parameter, namely the number of tapers $P$ or the smoothing matrix bandwidth $\bm{\Omega}$. Let us determine how we select them, and discuss how this might be done in theory as well as practice. We shall focus on 2D as this is the best explored dimensionality after 1D.

\subsection{Bandwidth selection in 2D}
We have already discussed the tapering in~\eqref{tapered3}. Some care has to be taken over the indices we loop over, especially the spatial index and the taper index.
Eqn~\eqref{tapered3} has to contain the same index over the taper as it is a modulus square. On the other hand in~\eqref{multitaper} we loop over the taper indices to combine the 1-d Cartesian tapers to make 2-d Cartesian tapers, e.g. we would use $00$ then $01$, $10$ as well as $11$. The two types of indices should not be confused.

With this, we are now ready to turn our theoretical understanding of
$I_p^t(k)$, the tapered periodogram, into practical methodology.  Starting from the multiple estimates of the spectrum we define the multitaper estimates
for $x\in B\in{\mathbb{R}}^2$
\begin{equation}
\label{multitaper}
\overline{I}_{P_1P_2}^t(\k)=\frac{1}{P_1P_2}\sum_{p_1=1}^{P_1}
\sum_{p_2=1}^{P_2} \widetilde{I}_{h_{p_1p_2}}(\k),\quad \k\in{\mathbb{R}}^2.
\end{equation}
Note that when looking at 1D tapers $p$ ranges from $1$ to $P_1$ but when going higher dimension we need to loop over each dimension in turn. The higher the dimension, the more
complex are our options. In 2D we shall simply explore Cartesian products of 1D tapers.

The number of tapers $P_1$ and $P_2$ in each of the two dimensions in~\eqref{multitaper} needs to be selected to balance variance and bias. For intuition we shall first understand the mean square error of a local averaging procedure as that defined in~\eqref{smoothingI}. We shall thus be inspired by regular $d$-dimensional smoothing to determine how to smooth. Starting from standard understanding available for $d$-dimensional density estimation~\citep{duong2005cross,wand1993comparison} we use the kernel function $W(k)$ with matrix-bandwidth $\bm{\Omega}$. Scaling the kernel function using the bandwidth matrix we obtain
\begin{equation}
W_{\bm{\Omega}}(\k)=|\bm{\Omega}|^{-1}\cdot W\left( \bm{\Omega}^{-1}\k\right),\quad \k\in {\mathbb{R}}^2.
\end{equation}
The bandwidth function satisfies
\begin{equation}
\label{normalization}
\sum_{w\in {\cal K}(B_\square({\bm{l}}))} W\left(w\right)=1.
\end{equation}
Note that this is not the same as the continuous version of the taper integrating to unity, see also the Section discusses the tapers.

With this definition we get the kernel density estimator
\begin{equation}
\label{multitaper2}
\overline{I}_h(\k;\mathbf{H})=\sum_{\k'\in {\cal K}(B_\square(\bm{l}))} \widetilde{I}_h(\k')W_{\bm{\Omega}}(\k-\k'),\quad \k\in{\mathbb{R}}^d.
\end{equation}
Note that~\eqref{multitaper} is not a special case of
~\eqref{multitaper2}. As each quadratic estimate
$I_p^t(\k)$ is arrived at by linear operations on the DFT,
and then implementing a modulus square, not all direct spectral estimates $I_p^t(\k)$ can be arrived at from smoothing the periodogram $I_0(\k)$.

\begin{Lemma}\label{biaslemma}
Assume $\widetilde{f}(\k)$ for all $\k\neq 0$ is twice differentiable in $\k$, and denote the mixed second derivative by $\frac{\partial^2 \widetilde{f}}{\partial w_i \partial w_j}=\widetilde{f}_{ij}(\k)$. Then the absolute bias of~\eqref{multitaper2} as an estimator of $f(\k)$ is given for $\k\neq 0$ by
\begin{align}
\label{biasexp}
\left|{\mathrm{Bias}}_0\left\{ \overline{I}_h(\k;\mathbf{H})\right\}\right|
&\leq \frac{\lambda^2}{2}\max_i|\widetilde{f}_{ii}(\k)|\tau^{2}_h,\\
\nonumber
\tau^{2}_h&=\sum_{\k'\in {\cal K}(B_\square(\bm{l}))}
(\k')\cdot \k' W_h(\k').
\end{align}
\end{Lemma}
\begin{proof}
See Appendix~\ref{biasproof}.
\end{proof}
The next step is the variance and the mean square error of our smoothed spectral estimator.
\begin{Lemma}\label{variancelemma}
Assume $\widetilde{f}(\k)$ for all $\k\neq 0$ is twice differentiable. Then the variance of~\eqref{multitaper2} as an estimator of $f(\k)$ is given for $\k\neq 0$ by
\begin{equation}
{\mathrm{\bf Var}}\{\overline{I}_h(\k;\bm{\Omega})\}=\frac{f(\k_n)-\lambda+\lambda^2+\lambda\|h\|^4_4}{\sigma_k \ell}+\lambda\|h\|^4_4.
\end{equation}
\end{Lemma}
\begin{proof}
See Appendix~\ref{varproof}.
\end{proof}

We put these results together and arrive at the optimal bandwidth.
\begin{Theorem}\label{mselemma}
We assume that $X$ is a point process with spectrum $f(\k)$ that is twice differentiable for all $\k\neq 0$.
Assume we smooth the periodogram using a box-car window with width $\sigma$, and take $\tilde{\sigma}=\sigma \ell h$ as the bandwidth measures in terms of the length $\ell$. In this case the mean square error takes the form of
\begin{align*}
{\mathrm{MSE}}\left\{ \widetilde{I}(\k_n)\right\}&\leq
\frac{\lambda^4}{4}\max_i|\widetilde{f}_{ii}(\k')|^2\left\{d(\tilde{\sigma}^2 )\frac{1}{12}\right\}^2+\frac{f(\k_n)-\lambda+\lambda^2+O(|B|^{-1})}{\tilde{\sigma}}+O(|B|^{-1}).
\end{align*}
This mean square error is minimized at
\begin{equation}
\label{sigma1}
\tilde{\sigma}=\left(12^2\frac{\lambda^2}{d^2\lambda^4\max_i|\widetilde{f}_{ii}(\k')|^2}\right)^{1/5}.
\end{equation}
\end{Theorem}
\begin{proof}
See Appendix~\ref{mseproof}.
\end{proof}
Naturally for a real data set, the parameters in this expression must be estimated, and the second derivative depends on the quantity we want to estimate. Given the theorem shows that an optimal smoothing bandwidth exists, it makes sense to pick one, and we propose, as we shall demonstrate in simulation studies, that visually a suitable degree of smoothing can be arrived at in examples. Thus the main point of this section is that an optimal wavenumber can be found, and it depends on the intensity of the point process as well as the roughness of the spectrum.

\subsection{Multitaper number}
This is all very well to find the bandwidth for the locally averaged estimator specified by
\eqref{multitaper2}, but we also want it for the multi--taper estimator in \eqref{multitaper}. The two most popular sequences of tapers to use are the prolate spheroidal wave sequences~\citep{landau1962prolate} and the sinusoidal functions~\citep{riedel1995minimum}. The former is used because it is maximally concentrated in wave number whilst being compact in time/space; the second because it produces the least bias formally.

Note that our choice is not trivial; some authors chose to interpolate the point process to a regularly spaced process, and then taper. We could chose to design taper to match the observed locations, e.g.~\cite{bronez1988spectral}, but the analysis of that choice would be complicated.

We shall use $P^d$ tapers in $d$ dimensions, and in this case we can derive the same expression for the bias as in~\eqref{biasexp} but we now need to use the properties of the discrete prolate spheroidal eigensequences--the number of well-concentrated tapers with a bandwidth $\sigma$ is
$P=2\ell \sigma$ for a cuboid data analysis window~\citep[p~369]{percival1993spectral}.  We then get
the same bias expression for each tapered periodogram as for $k>\sigma$ of
\begin{equation}
\label{expfrequency2}
{\mathbf{E}}\left\{I_p^t(\k;h)\right\}=\lambda  + \lambda^2\int_{\mathbb{R}^d}
\left| H_p(\k'-\k)\right|^2\widetilde{f}(\k')\;d\k',\quad p=1,\dots, P^d.
\end{equation}
We do this by mapping the numbering of each dimension into a single enumeration.
Having determined the expectation of each term we find that
\begin{align*}
{\mathbf{E}}\overline{I}_P(\k)&=\frac{1}{P^d}\sum_{p=1}^{P^d} {\mathbf{E}} I_p^t(\k)\\
&=\lambda+\lambda^2\int_{\mathbb{R}^d}
\frac{1}{P^d}\sum_{p=1}^{P^d}\left| H_p(\k'-\k)\right|^2\widetilde{f}(\k')\;d\k',
\end{align*}
this defining the spectral window of the taper by
\begin{equation*}
{\cal H}_P(\k)=\frac{1}{P^d}\sum_{p=1}^{P^d}\left| H_p(\k)\right|^2,\;p=1,\dots, P^d.
\end{equation*}
Using spheroidal wave functions, assuming the observational spatial domain is a Cartesian product of intervals of length $\ell$, we find that~\cite[p~342]{percival1993spectral}
with $P=2\ell \sigma$  then
\begin{equation}\label{av:window}\frac{1}{P^d}\sum_{p=1}^{P^d}\left| H_p(\k'-\k)\right|^2\approx
\left\{\begin{array}{lcr}
1 & {\mathrm{if}} & \forall \;|\k_j|<\sigma\\
0& {\mathrm{if}} & \exists \;j:\; \;|\k_j|<\sigma
\end{array}\right. .\end{equation}
This follows directly from the choice of direct product domains, and the use of the spheroidal wave functions.

We again expand the spectral density function to
\begin{equation}
\widetilde{f}(\k')=\widetilde{f}(\k)+\nabla \widetilde{f}(\k)\cdot (\k'-\k)+\frac{1}{2}(\k'-\k)\cdot \widetilde{\bm{H}}(\k'')(\k'-\k).
\end{equation}
We may as before note that
\begin{align*}
{\mathbf{E}}\overline{I}_P(\k)&=\lambda  + \lambda^2\int_{\mathbb{R}^d}
{\cal H}_P(\k'-\k)\left[\widetilde{f}(\k)+\nabla \widetilde{f}(\k)\cdot (\k'-\k)+\frac{1}{2}(\k'-\k)^T\widetilde{\bm{H}}(\k'')(\k'-\k)\right]\;d\k'\\
&=\lambda  + \lambda^2\widetilde{f}(\k)+\frac{\lambda^2}{2}\int_{\mathbb{R}^d}
{\cal H}_P(\k'-\k)(\k'-\k)^T\widetilde{\bm{H}}(\k'')(\k'-\k)\;
d\k'.
\end{align*}
We may argue about the bias in a manner analogous to before. We find:
\begin{align*}
|{\mathrm{bias}}\left\{\overline{I}_P(\k)\right\}|
&\leq \frac{\lambda^2}{2}\int_{\mathbb{R}^d}
{\cal H}_P(\k'-\k)\left|(\k'-\k)^T\widetilde{\bm{H}}(\k'')
(\k'-\k)\right|\;d\k'
\\
&\leq \frac{\lambda^2}{2}\max_{k,i}|\widetilde{f}_{ii}(\k')|^2 \int_{\mathbb{R}^d}
{\cal H}_P(\k'-\k)\sum_{i=1}^d \left(\k'-\k \right)_{ii}^2 \;d\k'\\
&=\frac{\lambda^2}{2}\max_{\k',i}|\widetilde{f}_{ii}(\k')|^2 \cdot d \cdot \frac{1}{12}\sigma^2.
\end{align*}
We now use that $P=2\ell \sigma$. For this reason our bias remains the same as before:
\begin{align*}
|{\mathrm{bias}}\left\{\overline{I}_P(\k)\right\}|
&\leq \frac{\lambda^2}{2}\max_{\k',i}|\widetilde{f}_{ii}(\k')|^2 \cdot d \cdot \frac{1}{12}(P/\{2\ell\})^2.
\end{align*}
The variance of the estimator comes from individual contributions from covariance expressions:
\begin{align*}
{\mathrm{\bf Var}} \{\overline{I}_P(\k)\}&=\frac{1}{P^d}\sum_{p=1}^{P^d}
\frac{1}{P^d}\sum_{q=1}^{P^d} {\mathrm{\bf Cov}}\{I_p^t(\k),I_q^t(\k)\}.
\end{align*}
To approximate the covariance we observe Theorem~\ref{Propvar1} which gives the form of the variance, where we ignore the higher order terms.

We then note
\begin{align}
{\mathrm{\bf Cov}}\{I_p^t(\k),I_q^t(\k)\}&\asymp \delta_{pq}\left\{
\lambda^2(1+\tilde{f}(\k))+\lambda {\cal O}(\|h\|_4^4)
\right\}+\lambda{\cal O}(\|h\|_4^4).
\end{align}
Because of the scaling we get
\begin{equation}
\int_B h_p^4(x)\,dx={\cal O}\left(\frac{1}{|B|} \right),
\end{equation}
and so putting these together we get that
\begin{align*}
{\mathrm{\bf Var}} \{\overline{I}_P(\k)\}&=\frac{1}{P^{2d}}\sum_{p=1}^{P^d}
\sum_{q=1}^{P^d} \left\{  \delta_{pq}\left\{
\lambda^2(1+\tilde{f}(\k))+\lambda {\cal O}(\frac{1}{|B|} )
\right\}+\lambda{\cal O}(\frac{1}{|B|} )\right\}\\
&=\frac{\lambda^2(1+\tilde{f}(\k))+ {\cal O}(\frac{1}{|B|} )}{P^{d}}+{\cal O}\left(\frac{1}{|B|} \right).
\end{align*}
We can then write
\begin{align*}
{\mathrm{MSE}}\{ \overline{I}_P(\k)\}&=
\left\{\frac{\lambda^2}{2}\max_{\k',i}|\widetilde{f}_{ii}(\k')|^2 \cdot d \cdot \frac{1}{12}(P/\{2\ell\})^2\right\}^2+\frac{1}{P^d}\left\{\lambda^2(1+\tilde{f}(\k))+{\cal O}\left(\frac{1}{|B|} \right)\right\}+{\cal O}\left(\frac{1}{|B|} \right)\\
&=\frac{\lambda^4}{4}\max_{i,\k'}|\widetilde{f}_{ii}(\k')|^2d^2P^4/\{16 \ell^4\}+\frac{1}{P^d}\left\{\lambda^2(1+\tilde{f}(\k))\right\}+{\cal O}\left(\frac{1}{|B|} \right).
\end{align*}
We can yet again minimize in $P$ and so use calculus to do so
\begin{align}
\nonumber
\frac{\partial}{\partial P}{\mathrm{MSE}}\{ \overline{I}_P(\k)\}&=\frac{\lambda^4}{16}d^2\max_{\k',i}|\widetilde{f}_{ii}(\k')|^2P^3/\ell^4-\frac{d}{P^{d+1}}\left\{\lambda^2(1+\tilde{f}(\k))+{\cal O}\left(\frac{1}{|B|} \right)\right\}\\
\Rightarrow P^{d+4}&=\frac{16\left\{\lambda^2(1+\tilde{f}(\k))+{\cal O}\left(\frac{1}{|B|} \right)\right\}\ell^4}{\lambda^4d^2\max_{\k',i}|\widetilde{f}_{ii}(\k')|^2}.
\label{Popt}
\end{align}
Again we see how the optimal taper number scales. There is a clear problem with the ``curse of dimensionality'', as the range of $P$ decays with $d$.
This indeed gives a very similar answer to our previous smoothing, and so we can estimate how many taper we need to use using the previous parameters. Again, the choice of taper necessitates knowledge of the second derivative of what we are trying to estimate--if we truly knew the second derivative we would not need to implement that estimation. Instead what the proof here shows is the scaling in the fundamental parameters of the process, and that a balance between variance and bias can be made. Furthermore this is not trivially done at the boundary of using $P^d$ or 1 taper.

In practice we do not have the unknowns in equations~\eqref{sigma1} and \eqref{Popt}.
But the advantage of determining their value is that we know that an optimum exists: we shall therefore empirically look for it using techniques applied in time series.
An approach in time series is to pre-select the wave number resolution $2W$, see for example~\citep[p.~369]{Percival1993SAPA}. The authors note that if the basic resolution is $W$ then we select $K=2NW\Delta t-1$.
In fact for a time series sampled at $\Delta t$ steps apart we call $2NW\Delta t$ the \underline{Shannon number}~\citep[p.~335]{Percival1993SAPA}. If $2W$ is the bandwidth in wave number of the data-taper concerntation operator, then $k=0,...,K-1$ where $K<2NW\Delta t$ will yield a set of tapers $\{h_k(\cdot)\}$ that are optimally concentrated to the wave number band $[-W\Delta t,W\Delta t]$.
So one method of selecting the number of data tapers is to just select the wave number resolution $w$, and then take the value of $K$ that matches the Shannon number.

\subsection{Plotting Spectral Estimates}
For regularly spaced time series it is very clear what frequencies/wave numbers we plot. We stick to the Fourier wave numbers, and plot up to the Nyquist wave number. But for a point process no such upper limit exists, as is also the case for randomly sampled stochastic processes, see~\citep{lii1994spectral}.

What we will find at larger wave numbers is that the baseline contribution of $\lambda$ dominates the spectrum. Unlike
time series that are regularly sampled in time, and it is natural to assume that the spectrum of the process is decaying in wave number, we can only assume the spectral deviations $\widetilde{f}(w)$ are decaying with $\|w\|$.

How do we then determine what is the highest wave number we can do spectral estimation? To determine this number we must guess at the decay of $\widetilde{f}(w)$. We shall assume that
\begin{equation}
\widetilde{f}(w)=C_0\|w\|^{-\alpha},\quad    \|w\|>\|w'\|
\end{equation}
If we knew $\|w'\|$ then our approach would be simple. We would take
\begin{align}
\nonumber
I(w_k)&=\hat{\lambda}+\hat{\lambda}^2C_0\|w_k\|^{-\alpha},\quad
k=1,\dots, K\\
\nonumber
Y_k&=I(w_k)-\hat{\lambda}=\hat{\lambda}^2C_0\|w_k\|^{-\alpha},\quad
k=1,\dots, K\\
Z_k&=\log|I(w_k)-\hat{\lambda}|=2\log \hat{\lambda}+\log C_0-\alpha \log \|w_k\|.
\end{align}
We can then assign $\|w_1\|$ as the wave number taking the value $\frac{1}{2}I(\frac{1}{\ell_1},\frac{1}{\ell_2})$ up to that taking the value $2\widehat{\lambda}$ using the Fourier grid and then solve for $(C_0,\alpha)$ using least squares, this giving us $(\hat{C}_0,\hat{\alpha})$.
With these estimates we solve
\[\hat{C}_0\|w_h\|^{-\hat{\alpha}}=\widehat{\lambda},\]
and use this as the highest wave number, $w_h$.

\section{Isotropy versus Anisotropy}\label{sec:isotropy}

We now consider the special case of isotropic point processes for which the complete covariance function depends only on distance, $\gamma(z)= \gamma_I(||z||)$ \citep[see also][p.~310--311]{daley2003introduction}. Such radial functions \citep{grafakos2013fourier} are special as we do not require any orientational specificity in our spectral representation. The isotropy of $\gamma(z)$ transfers through the Fourier transform to the sdf so $f(w)=f_I(||w||)$. We shall refer to the treatment of isotropic random fields when implementing analysis~\citep{Ponomarenko2006spectral}. An alternative is discussed by~\citep{errico1985spectra} and~\citep{durran2017practical}, for random fields.

The orientation invariance leads to dimensional reduction of the multidimensional Fourier transform. Three basic concepts are useful for the isotropic analysis: Recall the Bessel function of order $m>-1/2$
$$\mathcal{J}_m(t) = \frac{t^m}{2^m\pi^{1/2}\Gamma(m+1/2)}\int_0^\pi \sin^{2m}(u)~\cos[t\ \cos(u)]\;du\qquad t \ge 0.$$
It is connected to the $d$-dimensional Fourier transform via \citep[Sec 3]{vembu1961fourier}
$$
\int_{S^{d-1}}e^{-2\pi l i v\cdot u} du=2\pi~l^{-(d/2-1)}\mathcal{J}_{d/2-1}(2\pi l)\quad \forall v\in S^{d-1}, l\ge 0,
$$
we can also note relationships between radial and Cartesian Fourier representations, see for example ~\cite{grafakos2013fourier}, or~\cite{estrada2014radial}.

Recall also the (1-dimensional) Hankel transform of order $m$ of a function $g:[0,\infty]\mapsto{\mathbb{R}}$
$$\mathcal{H}_m[g](l) = \int_0^\infty g(r)\mathcal{J}_m(lr)r~dr\qquad l \ge 0.$$
Then if we assume $X$ is a stationary and isotropic process with a radial complete covariance function $\gamma_I$, the sdf at wavenumber $w$ with  $t=\|w\|$ takes the form
\begin{align}
\nonumber
	\label{eq:isotropicsdf}
    f_I(t=\|w\|):=f(w)&= \lambda + \int_{{\mathbb{R}^d}}\left\{\rho^{(2)}(z)-\lambda^2\right\}e^{-2\pi i w\cdot z}\;dz\\
        &= \lambda + \int_{0}^\infty\left\{\rho_I^{(2)}(r)-\lambda^2\right\} r^{d-1} \int_{S^{d-1}}e^{-2\pi ri~ w\cdot u}\;du\;dr\nonumber\\
    	&= \lambda +  2\pi\int_0^\infty \left\{\rho_I^{(2)}(r)-\lambda^2\right\} \left(\frac{r}{t}\right)^{d/2-1}\mathcal{J}_{d/2-1}(2\pi tr) r dr \nonumber\\
    	&\equiv \lambda +  2\pi~\mathcal{H}_{d/2-1}\left[\left\{\rho_I^{(2)}(r)-\lambda^2\right\} \left(\frac{r}{t}\right)^{d/2-1}\right](2\pi t),
\end{align}
i.e. the sdf is related to the Hankel transform of the complete covariance function.

And finally recall the radially averaged set covariance function,
$$\bar{\nu}_B(r) := \frac{1}{|S^{d-1}|}\int_{S^{d-1}}|B\cap B_{ru}|du , \quad r\ge 0$$
where $\nu_B(z) := |B\cap B_z|$ is the set covariance function of the set $B$ \citep[Appendix B.3]{IPSS2008}.

To gain insight into how we might estimate the radial sdf besides numerically rotation-averaging the $d$-dimensional periodogram using Eq.~\eqref{angularaveragedI}, let us analytically rotation-average Bartlett's periodogram (Eq.~\ref{Ih}): with $t\ge 0$
\begin{align}
    \overline{I}_0(t)&:=\frac{1}{|S^{d-1}|}\int_{S^{d-1}} {I}_0(tu)\, du\\
    \nonumber
    &=\frac{1}{|B||S^{d-1}|}\sum_{x\in X\cap B} \sum_{y\in X\cap B} \int_{S^{d-1}} e^{-2\pi i (t u)\cdot (x-y)}\, du\\
    \nonumber
    &= \hat\lambda + \frac{1}{|B||S^{d-1}|}\sum_{x,y\in X\cap B}^{\neq} \int_{S^{d-1}} e^{-2\pi t i (x-y)\cdot u}\ du\\
\label{eq:rawRA}
    &= \hat\lambda + \frac{2\pi}{|B||S^{d-1}| t^{d/2-1}}\sum_{x,y\in X\cap B}^{\neq} \mathcal{J}_{d/2-1}(2\pi t\|x-y\|)~\|x-y\|^{-(d/2-1)}
\end{align}
which in the planar case simplifies to
\begin{align}
\nonumber
&\stackrel{d=2}{=}~ \hat\lambda + \frac{1}{|B|}\sum_{x,y\in X\cap B}^{\neq} \mathcal{J}_{0}(2\pi t\|x-y\|).
\end{align}

This explains the motivation behind the isotropic estimator discussed by \cite{bartlett1964spectral} and \cite{diggle1987nonparametric}. Them authors prefer to normalize by $N(B)$ rather than by $|B|$, corresponding to further dividing \eqref{eq:rawRA} by $\hat\lambda$.

Like the classical periodogram, the isotropic "shortcut" estimator of Eq.\eqref{eq:rawRA} is highly biased near $0$, as has been observed by \cite{diggle1987nonparametric}. The biases are described in the following proposition.

\begin{Proposition}{Expectation of the isotropic Estimator, $I$.}
	\label{BesselE0}
	With the definition of the isotropic estimator given by \eqref{eq:rawRA},
	\begin{align}
	\nonumber
	{\mathbf{E}}\left\{ \bar{I}_0(t)\right\} =&
	\lambda + 2\pi~\mathcal{H}_{d/2-1}\left[\left\{\rho_I^{(2)}(r)-\lambda^2\right\} \left(\frac{r}{t}\right)^{d/2-1}~\frac{\bar{\nu}_B(r)}{|B|}\right](2\pi t) +
	2\pi\lambda^2~\mathcal{H}_{d/2-1}\left[\left(\frac{r}{t}\right)^{d/2-1}~\frac{\bar{\nu}_B(r)}{|B|}\right](2\pi t)
	\end{align}
\end{Proposition}
\begin{proof}
	Follows by applying the Campbell's theorem and adding and subtracting $$2\pi\lambda^2~\mathcal{H}_{d/2-1}\left[\left(\frac{r}{t}\right)^{d/2-1}~\frac{\bar{\nu}_B(r)}{|B|}\right](2\pi t) = \lambda^2\frac{1}{|S^{d-1}||B|}\int_{S^{d-1}}T(B, tu) du.$$
\end{proof}
If we compare the formulas for the isotropic sdf in Eq.\eqref{eq:isotropicsdf} and the expectation in Proposition \ref{BesselE0}, we recognise two sources of bias, just like with the periodogram: A convolution with a function $\bar{\nu}_B/|B|$, coming from the finite observation window like in the periodogram but this time radially averaged, and a centering bias term that is a radially averaged version of the $d$-dimensional bias. Diggle's truncation-to-local-minimum construction, see~\eqref{diggle4}, tries to correct for the biases, but among other things, it will not work well for a clustered process as the spectrum near zero will be underestimated.

To work towards a tapered estimator, consider rotation averaging the debiased, tapered estimator of \eqref{biascorrect}, in expectation,
\begin{align}
\label{orientav}
&{\mathbf{E}}\frac{1}{|S^{d-1}|}\int_{S^{d-1}}
\widetilde{I}_h(tu)\, du
\\
=&{\lambda}+\frac{2\pi}{|S^{d-1}|t^{d/2-1}}{\mathbf{E}}\left[ \sum_{x,y\in X\cap B}^{\neq} h(x) h(y) \mathcal{J}_{d/2-1}\left(t \| x-y\| \right) \| x-y\|^{-(d/2-1)}\right] -\frac{\lambda^2}{|S^{d-1}|}\int_{S^{d-1}}\left| H(t u)
\right|^2\,du. \nonumber
\end{align}

The problem is that a radial function $h_I$ might not exist for which $h_I(||x-y||)=h(x)h(y)$. But if we consider data-tapering the observed differences $\{x-y:x,y\in X\cap B\}$ instead of in $B$, we can device a taper $h_I$ in $B\oplus B = \bigcup_{x\in B}\{B+x\}$, thus we propose the following isotropic estimator:
\begin{align}
\label{eq:debiasedIsotropic}
\bar I_h(t) := \hat\lambda + \frac{2\pi}{|S^{d-1}| t^{d/2-1}|B|}\sum_{x,y\in X\cap B}^{\neq} h_I(x-y)\mathcal{J}_{d/2-1}(2\pi t\|x-y\|)~\|x-y\|^{-(d/2-1)} - \widehat{\lambda^2}Bias_I(t)
\end{align}
where the bias correction term consists of the rotation average of both the taper and the set covariance of $B$,
\begin{eqnarray*}
    Bias_I(t)&=& 2\pi\mathcal{H}_{d/2-1}\left[\left(\frac{r}{t}\right)^{d/2-1}\frac{\bar{p}_{B,h}}{|B|}\right](2\pi t), \quad \text{where}\\
    \bar{p}_{B,h}(r) &=& |S^{d-1}|^{-1}\int_{S^{d-1}}|B\cap B_{ru}|h_I(ru)du,
\end{eqnarray*}
and depends only on $t=|w|$. We warn the reader that the estimator of~\eqref{eq:debiasedIsotropic} is not necessarily positive, but like previous authors ~\citep{lii1994spectral,bandyopadhyay2009asymptotic}, we sacrifice this to remove large part of the bias. Additional issue arises from estimation of the required $\lambda^2$, as the standard estimator $\hat{\lambda}^2$ is biased with a term depending on $\lambda/|B|$ and a term depending on the unknown second order product density of the underlying process. In our examples we use the slightly less biased $\widehat{\lambda^2}=N(B)[N(B)-1]/|B|^2$ which removes the first order bias. On the positive side, do note that the debiased isotropic estimator without particular tapering (i.e. $h_I\equiv 1$) can be estimated without any tuning parameters.

We may now ask what radial taper to use. Tapers are normally designed for processes observed in discrete time or space. We have discussed separable filters, and now note that the choice of filter needs to match the observational domain. For example it is difficult to match a square observational domain with a radial filter, unless one chooses a compact filter inscribed by the box. There are therefore three considerations 1) demanding the taper be separable in space, 2) isotropic in space, or 3) exactly compact in space. Discrete and isotropic filters have been been determined numerically by~\cite{simons2011spatiospectral}. This will not be an option for us as we need to evaluate the tapers at random locations.

Slepian and co-authors have studied the design of tapers in arbitrary dimensions~\citep{slepian1964prolate}. Often solutions are in terms of prolate spheroidal wavefunctions. Their discrete analogue is the prolate spheroidal wave sequences, which are similar to the set of Hermite function~\citep{xu1999multiple}. However, the isotropic estimator can be seen as an estimator for the Fourier transform of the (non-stationary) first moment of the difference process $Z=\{z=x-y:x,y\in X\cap B\}$, with the connection $\rho^{(2)}(z) = \lambda_Z(z)/|B\cap B_z|$.

Assuming the observation window is a cuboid $B=\prod_{j=1}^d[0,\ell_j]$, then the difference vector observation window is $B\oplus B=\prod_{j=1}^d[-\ell_j,\ell_j]$, and we can create the $d$-dimensional taper related to the Hermite functions as the product of the squared exponentials
\begin{equation}
    h_I(z) = \prod_{j=1}^d \exp\left[-az_j^2/(2\ell_j)^2\right], \quad z\in{\mathbb{R}^d}
\end{equation}
We see that if $\ell_j=\ell$ for all $j$, then $h_I$ becomes radial. This is related to circular harmonic decompositions, see~\citep{jacovitti2000multiresolution}. The authors of~\citep{xu1999multiple} recommend using a scaling factor which we have adjusted for a sampling region of $z_j\in (-\ell_j,\ell_j)$, but which they match to the spheroidal wave functions. We propose setting $a=25$, as this numerically seems to not down--weight too much data, or have a too large jump near the border of $B\oplus B$.

In a sense our above construction is a ``fix'' as it only works for the 0th Hermite function, and using more than one taper will even if $\ell_j=\ell$ for all $j$ not lead to isotropic functions. This is not unsurprising as the sampling domain will leave an inprint, and capturing this set of information requires using more than the first taper. Also, if we have a spatial sampling that is highly elongated rectangle rather than a square then another solution needs to be chosen.

Finally we may ask ourselves, what happens if we have a point process $X$ which is anisotropic but we still evaluate~\eqref{diggle} in 2D for the pattern. Because it incorporates angular averaging we expect the RHS of the expression to only depend on $\|w\|$ as it does.

\begin{Proposition}{Expectation of Diggle's 2D Estimator II.}\label{BesselE}
The large area expectation of Diggle's estimator given in~\eqref{diggle} when $B$ is a cuboidal box when the spectrum is not necessarily isotropic takes the form in terms of $\widetilde{f}(w)$ as defined in~\eqref{eqn:fourierexcess} of
\begin{align}
    {\mathbf{E}}\{ I_D(w)\}&=\lambda+\lambda^2\int_{{\mathbb{R}}^d}
    \widetilde{f}(v)\frac{\delta(\|w\|-\|v\|)}{\|w\|}dv
+\int_{{\mathbb{R}}^d}\frac{\delta(\|w\|-\|v\|)}{\|w\|}
|B|^{-1}T(B,w-v)dv+o(1),
\end{align}
where $T(B,w)$ is defined in Lemma~\ref{lemma1}.
\end{Proposition}
\begin{proof}
The proof is provided in Appendix~\ref{proofBesselE}.
\end{proof}

\section{Simulation studies}\label{sec:sims}

To study the behaviour of the estimators we have discussed, and to verify some of the debiasing results, we conducted a simulation study.

\subsection{Details for executing the simulations study}

We simulated three stationary models in 2D to represent the archetypal spatial patterns of regularity, complete randomness and clustering. Throughout the simulations the intensity, i.e. expected number of points per unit area, was kept constant at $\lambda\equiv 0.01$. The processes were simulated in increasing square observation windows $B_n=[-\ell_n/2, \ell_n/2]^2$ such that on average we observed $\lambda|B_n| = \lambda \ell_{n}^2 = n = 25, 50, 100, 200, 400$ or $800$ points. For example, $\ell_{n=100}=100$ and $\ell_{n=400}=200$ spatial units. We chose the regime both to explore the often encountered small data scenarios and to ascertain asymptotic behaviour. The models below are tuned so that they produce patterns with visually distinct structures at $\ell_{n=100}$.

Complete randomness i.e. the Poisson process has no adjustable parameters after $\lambda$ was fixed. Spatial regularity is represented by the Mat\'ern type II process \citep[Section 6.5.2]{IPSS2008} in two variants, one with hard-core radius $2$ (variant "r2") and one with radius $5$ (variant "r5") which correspond to about 40\% and 90\% of the maximum allowed radius for the process given $\lambda$, respectively. The clustering behaviour is represented by the Thomas process \citep[Section 6.3.2]{IPSS2008} with cluster intensity $\kappa$ and Gaussian dispersal kernel standard deviation $\sigma$ and again in two variants, one with many small or tight clusters ($\kappa = 0.6\lambda, \sigma=2$, variant "MS") and one with a few large clusters ($\kappa=0.3\lambda, \sigma=6$, variant "FL"). The per-cluster expected point count $\mu$ got then fixed via the property $\mu\kappa=\lambda$ to $1.67$ and $3.33$, respectively. The theoretical pair correlation functions and spectral density functions of the models are depicted in Figure \ref{fig:fig-1}. Example patterns of each model and observation window combination are given in Supplement Figure \ref{figS:patterns}.

For the estimation of the 2D spectra we fixed the wavenumber grid to $w\in[-0.3, 0.3]^2$ with 101 equidistant steps in each dimension, giving the step~size~0.006 in both dimensions. This scale was chosen to cover the interesting range of non-constant values for all models (cf. Figure \ref{fig:fig-1}). To reduce the effect of high wave number noise only the values on the sub-grid $w\in[-0.2,0.2]^2$ were considered when integrating over $w$ for the quality metrics discussed below. When a rotationally averaged curve was to be computed using Eq. \eqref{angularaveragedI}, the averaging was done over the magnitude-grid $|w|=0.003, 0.006,...,0.300$ using a box kernel and if not otherwise stated, radius of $1.25\cdot 0.006$. The radial Hermite taper parameter was fixed to $a=25$.

We summarised the quality of each estimator $I_.=I_.(w;{\mathbf{x}})$ under all combinations of a model $M$ and observation window $B_n$, say $M_n$, with an integrated summary. First, we estimated per-wavenumber variance $V(w;M_n)=\ \text{Var}\left[I_.(w;{\mathbf{x}}) | {\mathbf{x}} \sim  M_n\right]$, bias $Bias(w; M_n)=\ {\mathbf{E}} \left[I_.(w;{\mathbf{x}}) - f_{M}(w)|{\mathbf{x}}\sim M_n\right]$ and mean square error $MSE(w;M_n)\ =\ V(w;M_n) + Bias^2(w;M_n)$. The $f_{M}$ stand for the theoretical sdf of the model $M$. Then these were summarised further to $iVar(M_n)\ =\sum_w V(w;M_n)$, the integrated squared bias $iBias^2(M_n)=\sum_w Bias^2(w;M_n)$ and $iMSE(M_n) = iVar(M_n) + iBias^2(M_n)$. Smaller $iVar, iBias^2$ and $iMSE$ indicate better quality. The quantities were estimated from 1000 simulations of every $M_n$.

The estimator for a periodogram is given in equation \eqref{biascorrect} with the debiasing term included. Bartlett's periodogram ("periodogram") has the taper $h_0(x)=1_{B_n}(x)/\sqrt{|B_n|}$ for window $B_n$. For the multitapering ("mt $P$"), defined in equation \eqref{tapered} with $P^2$ tapers each having a parameter $p~=~(p_1, p_2)~\in~\{1,...,P\}^2$, we used the orthogonal sine-tapers  $h_p(x)=\mathbf{1}_{B_n}(x)\prod_{j=1}^2\sin[ \pi p_j (x_j+l_n/2)/l_n]$. For the kernel smoothed estimators ("smoothed $m$"), we first compute the Bartlett's periodogram and then, given the estimate $\hat I_0$ on the wavenumber gird, we convolute it with a $m\times m$ discrete template having approximately Gaussian weights.

All computations were done using the R-software \citep[v.3.6.3]{Rsoftware}. Simulations were done with the help the R-package \texttt{spatstat} \citep[v.1.64-1]{spatstat}. Each of the discussed estimators were programmed into an R-package \texttt{ppspectral}, available on request from the first author.

\subsection{The effect of debiasing}

As we saw in Figure \ref{fig:dataexample}, the debiasing term has a large effect on the quality of the estimates. Table \ref{tab:bias} provides the fractions of bias removed by the debiasing term for the Bartlett's periodogram and $P=3$ multitapered periodogram. Nearly all bias is removed, with the only notable exception being the 65\% bias left in the multitapered periodogram when data is a (arguably very) small sample of Thomas FL variant.

\begin{table}[ht]
	\centering
	\begin{tabular}{llrrrrrr}
		\hline
		Model & Estimator & $n=$25 & 50 & 100 & 200 & 400 & 800 \\
		\hline
		MaternII r5 & periodogram & 1.00 & 1.00 & 1.00 & 0.98 & 0.99 & 1.00 \\
		MaternII r5 & mt P=3 & 0.99 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
		MaternII r2 & periodogram & 1.00 & 1.00 & 1.00 & 0.99 & 0.99 & 1.00 \\
		MaternII r2 & mt P=3 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
		Poisson & periodogram & 1.00 & 1.00 & 1.00 & 0.98 & 0.99 & 1.00 \\
		Poisson & mt P=3 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 & 1.00 \\
		Thomas FL & periodogram & 0.93 & 0.98 & 0.99 & 0.97 & 0.99 & 1.00 \\
		Thomas FL & mt P=3 & 0.35 & 0.86 & 0.98 & 1.00 & 1.00 & 1.00 \\
		Thomas MS & periodogram & 0.98 & 0.99 & 1.00 & 0.97 & 0.99 & 1.00 \\
		Thomas MS & mt P=3 & 0.89 & 0.98 & 1.00 & 1.00 & 1.00 & 1.00 \\
		\hline
	\end{tabular}

\caption{The fraction of bias removed by the proposed bias correction. Maximum is 1.00.}
\label{tab:bias}
\end{table}

We illustrate the form of the biases in Figure \ref{fig:bias2d} for two model variants. The centering bias is always positive, as expected from equation \eqref{expectation}. Bartlett's periodogram bias (cf. Lemma 5.1) is concentrated along the axis where one of the sinc-functions is constant 1 and near the origin where the second sinc function grows to 1. The multitapered periodogram has a square-shape around the origin.

\begin{figure}[!h]
	\centering
	\includegraphics[width=.8\linewidth]{fig-2x5-bias2d.pdf}
	\caption{The effect of bias removal for two of the models described in the text. Theoretical sdf on the left, and then towards the right, estimated mean of: Bartlett's periodogram, biased and debiased; Multitapered $P=3$, biased and debiased. Estimated with sample sizes $n=800$.}
	\label{fig:bias2d}
\end{figure}

\subsection{Overall quality of the estimators}

Figure \ref{fig:Quality} summarizes the estimated integrated variance, squared bias and MSE of the estimators with various parameters. Only the debiased estimators are included, and the values are given in $\log_{10}$-scale and relative to the baseline given by Bartlett's periodogram.

Multitaper with $P=1$ (i.e. only one taper) does not differ in iMSE from the periodogram as the benefit of averaging does not apply. With more than one taper the variance is reduced, e.g. $P=3$ (9 tapers) it goes down by 90\%, and the bias stays at baseline. But adding more and more tapers ($P=6$, 36 tapers) starts to introduce bias especially for low sample sizes. The bias is more prominent for the clustered cases, and based on the earlier discussion this is due to oversmoothing at small wavenumbers. Both the variance reduction and the oversmoothing happens with the post-hoc smoothing as well.

From the results we can confirm that some smoothing is beneficial but again a balance must be struck to avoid oversmoothing (cf. Section \ref{sec:linear}). The regular process is more resilient to oversmoothing with reasonable sample sizes.

\begin{figure}[!h]
	\centering
	\includegraphics[width=1\linewidth]{fig_3x5_overall_quality_relative.pdf}
	\caption{Relative integrated variance (top row), squared integrated bias (middle), and integrate mean square error (bottom) for different models and estimators of spectral density at various levels of sample size $n$. The values are relative to debiased Bartlett's periodogram, which has the value 0 in the presented $\log_{10}$-scale. For example, a value of -1 means the error is of order 0.1 relative to the periodogram.}
	\label{fig:Quality}
\end{figure}

\subsection{Rotational averaging and isotropic estimation}

Isotropic estimation, either by averaging radially or by the isotropic shortcuts discussed in Section 7, is an important dimension reduction step for easier visual inspection of the periodogram (particularly for 3D data). The kernel-based rotation averaging depends on a bandwidth parameter, and we studied how this affects the quality metrics. We set the bandwidth radius to $BF\cdot 0.006$ and varied the bandwidth factor $BF$ from 1 to 10. Figure \ref{fig:rotave-bandwidth} provides the $iBias^2$ and $iMSE$ for the debiased Bartlett's periodogram when comparing the rotationally averaged estimates to the true 1D isotropic sdfs (cf. Figure \ref{fig:fig-1}).

The Poisson process has a constant sdf so oversmoothing will not be penalised. For the clustered processes even a small bandwidth introduces bias, but this is typically more than offset by the reduction in variance. The error for the few large clusters variant is more sensitive to changes in the bandwidth than the many small clusters. This is because the effective range of wavenumbers where the sdf exhibits structure is more concentrated, and thus easier to smoothen out, than the few large cluster variant (cf. Figure \ref{fig:fig-1}).

The strongly regular variant "r5" is similar to the clustered cases in that an optimal bandwidth is easy clearly present. For the less regular "r2" variant the error is like that of the Poisson process, possibly due to the target sdf being a weak oscillation around a constant and this being hard  to detect. If we look the integrated squared bias we however do see that with large enough data oversmoothing is detectable.

We note that since multitapered periodograms and locally averaged periodograms are already smoothed, their rotational averaged estimates exhibit smoothing bias at smaller bandwidth factors (results not worth showing here).

\begin{figure}[h]
	\centering
	\includegraphics[width=1\linewidth]{fig-2x5-rotave-bandwidth-ibias-imse}
	\caption{Integrated squared bias (top) and integrated MSE of the rotationally averaged Bartlett's periodogram when compared to the true 1D isotropic sdf. The rotational averaging bandwidth was $BF\cdot 0.006$, where the bandwidth factor $BF$ is on the x-axis. The models are described in the text.}
	\label{fig:rotave-bandwidth}
\end{figure}

To compare the radial averaging to direct isotropic estimation, we estimated the debiased isotropic periodogram and the squared-exponential tapered isotropic periodogram with $a=25$. Table \ref{tab:isotropic} provides the relative quality scores for the isotropic and the tapered isotropic periodograms relative to radially averaged periodogram with the best bandwidth for each model and data size, cf. Figure \ref{fig:rotave-bandwidth}. The values are medians over the sample sizes $n\ge 50$.

In iMSE sense the radially averaged periodogram performs the best, chiefly because iMSE is of order iVar and the radial averaging greatly reduces the variance. However, in terms of bias the isotropic estimator provides more accuracy, except for Poisson when no oversmoothing is possible. Adding a smooth taper to the isotropic estimator further reduces the bias, but can increase the variance, likely due to some data near the edges is being filtered out. The bias occurs near $|w|=0$, mostly when $|w| < 1/(\sqrt{2}l_n)$ (illustration in Supplements Figure \ref{figS:rotave-v-isotropic}).

\begin{table}[ht]
	\centering
	\begin{tabular}{l|rrr|rrr}
		\hline
		& \multicolumn{3}{c}{No taper} &  \multicolumn{3}{c}{Squared-exponential taper} \\
				\hline
	    Model    & iMSE & iBias$^2$ & iVar  & iMSE & iBias$^2$ & iVar\\
	\hline
MaternII r5 & 6.12 & 0.21 & 8.57 & 3.51 & 0.29 & 4.90 \\
  MaternII r2 & 22.40 & 3.52 & 23.07 & 15.59 & 1.09 & 16.12 \\
  Poisson & 15.72 & 10.53 & 15.89 & 13.76 & 5.75 & 13.77 \\
  Thomas FL & 1.22 & 0.04 & 1.81 & 1.32 & 0.14 & 1.97 \\
  Thomas MS & 4.34 & 0.01 & 6.50 & 3.58 & 0.04 & 5.35 \\
\hline
	\end{tabular}
\caption{Scores $iMSE,iBias$ and $iVar$, relative to radially averaged debiased periodogram with optimised smoothing bandwidth.}
\label{tab:isotropic}
\end{table}

\section{Discussion}\label{sec:disc}
The Fourier transform of a point process has been treated by numerous authors~\citep{mugglestone1996exploratory,bartlett1964spectral,bartlett1963spectral, diggle1987nonparametric}.
Despite this fact we have seen limited methodological development, leaving methodological understanding lagging much behind that for time series and random fields.
There are a number of barriers for using and understanding the DFT of a point process. The first is computational: unlike the FFT the act of calculating the DFT is not ${\cal O}(N\log N)$, where $N$ are the number of observed points.

Second, there has been an often said, that Fourier analysis of point processes are not very different from that of random fields.  However, the results we have developed here, show that the spectral moments are different between random fields and point processes, e.g. compare to~\citep{stein2012interpolation}. Following on from this realisation we have found a way of avoiding the treatment of low wave numbers differently.
An important aspect of the spectral moments is the impact of the spatial domain we observe the point process in: the geometry of this object determines what wave numbers we should evaluate the Fourier transform at. Unlike for a regularly sampled stochastic process like a random field, the point process is biased at low wave numbers. This is because, as we demonstrate, we need to de--mean the spectral measure.
By understanding the moments of the Fourier transform we show that we may remove this bias completely. Our realisation in this issue is that the sampling domain determines both bias and second order structure for a point process and needs to be studied directly. This has implications for our sampling of the wave number grid, where once we have debiased estimates can be formed at any wave numbers.

Third, we discuss various forms of variance reduction of a spectral variance estimator (the periodogram). We propose either a kernel smoothing of the debiased periodogram, or multitapering combined with a debiasing step. Our studies agree with the theory that unless the second derivative of the spectrum is zero, there is an optimal smoothing width. We also discuss constructing one-dimensional summaries from the two-dimensional spectrum, and tapering in many dimensions.

We have here taken the first steps to fully developing the methodology for spectral analysis of point processes, by calculating the moments of spectral summaries, deducing Gaussianity of the DFT, and fully determining how the wave number grid is influenced by the geometry of the region we observe the point process in. We imagine that authors starting from results such as those provided by this paper will seek to fill in the space remaining in between technical understanding of point processes~\citep{daley2003introduction} to their practical application~\citep{lavancier_determinantal_2015,Moller2007}.

\section{Acknowledgements}
This work was supported in part by
by the UK Engineering and Physical Sciences Research Council
under Mathematical Sciences Leadership Fellowship EP/I005250/1, Developing Leaders Award EP/L001519/1, and Award EP/N007336/1;  and the European Research Council under Grant CoG
2015-682172NETS within the Seventh European Union Framework Program.

\bibliographystyle{plainnat}