# On the excursion area of perturbed Gaussian fields

Elena Di Bernardino, Anne Estrade, Maurizia Rossi

► To cite this version:

Elena Di Bernardino, Anne Estrade, Maurizia Rossi. On the excursion area of perturbed Gaussian fields. ESAIM: Probability and Statistics, 2020, 24, pp.252-274. 10.1051/ps/2020002. hal-02904016

HAL Id: hal-02904016

https://hal.science/hal-02904016

Submitted on 21 Jul 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# ON THE EXCURSION AREA OF PERTURBED GAUSSIAN FIELDS

ELENA DI BERNARDINO<fnref n="1" />*, ANNE ESTRADE<fnref n="2" /> AND MAURIZIA ROSSI<fnref n="3" />

**Abstract.** We investigate Lipschitz-Killing curvatures for excursion sets of random fields on $\mathbb{R}^2$ under a very specific perturbation, namely a small spatial-invariant random perturbation with zero mean. An expansion formula for mean curvatures is derived when the magnitude of the perturbation vanishes, which recovers the Gaussian Kinematic Formula at the limit by contiguity of the model. We develop an asymptotic study of the perturbed excursion area behaviour that leads to a quantitative non-Gaussian limit theorem, in Wasserstein distance, for fixed small perturbations and growing domain. When letting both the perturbation vanish and the domain grow, a standard Central Limit Theorem follows. Taking advantage of these results, we propose an estimator for the perturbation variance which turns out to be asymptotically normal and unbiased, allowing to make inference through sparse information on the field.

**Mathematics Subject Classification.** 60G60, 60F05, 60G15, 62M40, 62F12.

Received June 3, 2019. Accepted January 13, 2020.

## 1. INTRODUCTION

A wide range of phenomena can be seen as single realizations of a random field, for instance the Cosmic Microwave Background radiation (CMB) (see [24]), medical images of brain activity (see [34]) and of mammary tissue (see [11]), and many others. Their features can be investigated through geometrical functionals, among them the well-known class of Lipschitz-Killing (LK) curvatures of excursion sets (see e.g. [31, 33] for a precise definition and [16] for some applications in cosmology). From a theoretical point of view, probabilistic and statistical properties of the latter have been widely studied in the last decades. For instance, in the two-dimensional Euclidean setting, in [5, 7, 12], the length of the level sets (i.e., the perimeter of the excursion sets) is taken into account, in [15] the Euler-Poincaré characteristic, while several limit theorems are obtained for the excursion area in [10, 32]. See [22, 27] for higher dimensions. In this manuscript we focus on the two-dimensional setting, i.e., random fields defined on $\mathbb{R}^2$ endowed with the standard Euclidean metric.

In many cases, the LK curvatures are studied for Gaussian excursion sets via the Gaussian Kinematic Formula (see, e.g., [1, 8]). In this framework, a natural question is the following: *how do these geometric quantities change under small perturbations of the underlying field?* The present work gives an answer in the very specific case

<footnote>*Keywords and phrases:* LK curvatures, Gaussian fields, perturbed fields, quantitative limit theorems, sojourn times, sparse inference for random fields, spatial-invariant random perturbations.</footnote>

<footnote n="1">¹ Conservatoire National des Arts et Métiers, Paris, EA4629, 292 rue Saint-Martin, Paris Cedex 03, France.</footnote>

<footnote n="2">² MAP5 UMR CNRS 8145, Université Paris Descartes, 45 rue des Saints-Pères, 75270 Paris Cedex 06, France.</footnote>

<footnote n="3">³ Dipartimento di Matematica e Applicazioni, Università degli Studi di Milano-Bicocca Via Roberto Cozzi, 55, 20125 Milano MI, Italy.</footnote>

<footnote>* Corresponding author: elena.dibernardino@lecnam.net</footnote>
of an independent, additive, spatial-invariant noise with zero mean perturbing a stationary isotropic Gaussian field. Indeed, this model naturally arises when taking into account measurement errors that globally affect all the observations in physical experiments. As briefly anticipated above, LK curvatures have been very extensively exploited in the recent cosmological literature as a tool to probe non-Gaussianity and anisotropies in the CMB (see e.g., [13, 16]). Our setting could be viewed as the representation of a Gaussian field contaminated by superimposed point sources (i.e., galaxies and other astrophysical objects), and in this sense it could be used for point source detection or map validation in the framework of CMB data analysis. We remark that the perturbation of a Gaussian field obtained by adding either an independent Gaussian field or a function of the field itself can be fully treated through Gaussian techniques (see [6]). Some computations of expected values of LK curvatures of excursion sets in the latter setting have been given in the physical literature by e.g. [18, 26] in order to derive a promising method to constrain the primordial non-Gaussianity of the universe by temperature fluctuations in the CMB.

However, our aim is not to develop the theory of LK curvatures for excursion sets of general non-Gaussian fields (see for instance [3, 7, 23]) but to go beyond Gaussianity by introducing a small perturbation of the underlying Gaussian field. This perturbation clearly appears in the LK curvatures of excursion sets, allowing us to measure the discrepancy between the original and perturbed fields. Moreover, we are able to recover the classical Gaussian case by contiguity, i.e., when the perturbation vanishes.

Our model can be also seen as a random affine transformation of the initial excursion level. A deterministic and more challenging counterpart has been recently studied in [4] where a different geometrical functional is considered, namely, the number of connected components of excursion sets.

At last, considering excursion sets instead of the whole field is a sparse information that is commonly used by practitioners (see, for instance, Chap. 5 in [2]). Furthermore, it is equivalent to consider thresholded fields, which is a standard model in physics literature (see, e.g., [9, 29, 30]).

## 1.1. Main contributions

In this paper, we provide an expansion formula for the perturbed LK curvatures (see Prop. 2.5) where the contiguity property clearly appears for a vanishing perturbation. Visually, the perturbation is not evident by looking at the image of excursion sets (see Fig. 1) but its impact can be detected by a statistical and numerical analysis through the evaluation of their LK curvatures (see Fig. 3). Moreover, an asymptotically normal and unbiased estimator for the variance of the perturbation magnitude is proposed in Proposition 4.1. In order to get the Gaussian limiting behaviour of the latter estimator, we develop an asymptotic study of the second LK curvature, i.e. the area, of the perturbed excursion sets. We analyze both the case when the perturbation vanishes and the domain of observation grows to $\mathbb{R}^2$ (see Thm. 3.5) and the case of a fixed small perturbation and growing domain (see Thm. 3.2). The former is a standard CLT result, the latter is a quantitative limit theorem towards a non-Gaussian distribution, giving an upper bound for the convergence rate in Wasserstein distance. We deeply study the unusual non-Gaussian limiting law (see Thm. 3.3 and Figs. 4, 5 and 6).

An auxiliary result which is of some interest for its own is collected in Lemma A.1 where uniform rates (w.r.t. the level) of convergence for sojourn times of general Gaussian fields are proved. An argument similar to the one in the proof of Lemma A.1 allows to obtain uniform rates of convergence also for sojourn times of random hyperspherical harmonics, at the cost of getting worse rates than those found in [25].

## 1.2. Outline of the paper

Section 2 is devoted to the study of mean LK curvatures of our perturbed model. In particular, in Section 2.1 we introduce our setting and recall the notion of LK curvatures for excursion sets; in Section 2.2 we give unbiased estimators of the LK densities; in Section 2.3 we derive the asymptotic expansion for the mean curvatures as the perturbation vanishes (Prop. 2.5) providing some numerical evidence in Section 2.4 (see Figs. 2 and 3). In Section 3.1 we state and prove the quantitative limit theorem, in Wasserstein distance, for the excursion area of the perturbed model for fixed small perturbations and growing domain (Thm. 3.2). Theorem 3.3 characterizes
the unusual non-Gaussian limiting distribution whose numerical investigation leads to Figures 5 and 6. In Section 3.2 we state and prove the standard CLT for the excursion area for growing domain and disappearing perturbation (Thm. 3.5).

Taking advantage of the asymptotic studies for LK curvatures in Sections 2 and 3.1, in Proposition 4.1 we prove that the proposed estimator for the perturbation variance is unbiased and asymptotically normal. Its performance can be appreciated in Figure 7.

Finally, Appendix A collects the auxiliary result on uniform rates of convergence for sojourn times of Gaussian fields.

## 2. LK CURVATURES FOR THE CONSIDERED PERTURBED GAUSSIAN MODEL

### 2.1. Definitions and preliminary notions

In the present paper we consider the three additive functionals, called in the literature intrinsic volumes, Minkowski functionals or Lipschitz-Killing curvatures, $L_j$ for $j = 0, 1, 2$, defined on subsets of Borelians in $\mathbb{R}^2$. Roughly speaking, for A a Borelian set in $\mathbb{R}^2$, $L_0(A)$ stands for the Euler characteristic of A, $L_1(A)$ for the half perimeter of its boundary and $L_2(A)$ is equal to its area, i.e., the two-dimensional Lebesgue measure. Taking inspiration from the unidimensional framework, the $L_2$ functional is also called sojourn time, although no time is involved in this context.

**Notations.** All over the paper, $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{R}^2$ and $I_2$ the $2 \times 2$ identity matrix.

We will also denote by $| \cdot |$ the two-dimensional Lebesgue measure of any Borelian set in $\mathbb{R}^2$ and by $| \cdot |_1$ its one-dimensional Hausdorff measure. In particular, when $T$ is a bounded rectangle in $\mathbb{R}^2$ with non empty interior,

$$L_0(T) = 1, \quad L_1(T) = \frac{1}{2}|\partial T|_1, \quad L_2(T) = |T|,$$

where $\partial T$ stands for the boundary of $T$.

Let $T$ be a bounded closed rectangle in $\mathbb{R}^2$ with $0 \in \partial T$. In the following notation $T \nearrow \mathbb{R}^2$ stands for the limit along any sequence of bounded rectangles that grows to $\mathbb{R}^2$. For that, set $N > 0$ and define

$$T^{(N)} := \{Nt : t \in T\}$$

the image of a fixed rectangle $T$ by the dilatation $t \mapsto Nt$; then letting $T \nearrow \mathbb{R}^2$ is equivalent to $N \to \infty$, as soon as $0 \in \partial T$. Remark that $T^{(N)}$ is a Van Hove (VH)-growing sequence (see Def. 6 in [10]), i.e., $|\partial T^{(N)}|_1/|T^{(N)}| \to 0$ as $N \to \infty$. In the sequel, we sometimes drop the dependency in $N$ of the rectangle $T^{(N)}$ to soften notation.

We now define the main notions that we will deal with.

**Definition 2.1 (Considered Gaussian field).** Let $g$ be a Gaussian random field defined on $\mathbb{R}^2$ with $C^3$ trajectories. We assume that $g$ is

- stationary, isotropic with $\mathbb{E}[g(0)] = 0$, $\mathrm{Var} g(0) = \sigma_g^2$ and covariance matrix of $\nabla g(0)$ equals to $\lambda I_2$ for some $\lambda > 0$, $\sigma_g > 0$,

- such that its covariance function $r(t) = \mathrm{Cov}(g(0), g(t))$ satisfies

$$|r(t)| = O(\|t\|^{-\alpha}), \quad \text{for some } \alpha > 2 \text{ as } \|t\| \to \infty. \qquad (2.1)$$

We will consider perturbations of the above Gaussian field prescribed by the following.

**Definition 2.2 (Perturbed Gaussian field).** Let $X$ be a random variable such that $\mathbb{E}[|X|^3] < +\infty$. Let $g$ be a Gaussian random field as in Definition 2.1, with $X$ independent of $g$. We consider the following perturbed field
$$f(t) = g(t) + \epsilon X, \ t \in \mathbb{R}^2, \text{ with } \epsilon > 0.$$

Let $u \in \mathbb{R}$ and $T$ a bounded rectangle in $\mathbb{R}^2$. For $h$ any real-valued stationary random field, we consider the excursion set within $T$ above level $u$:

$$\{t \in T : h(t) \geq u\} = T \cap E_h(u), \quad \text{where } E_h(u) := h^{-1}([u, +\infty)).$$

We now introduce the Lipschitz-Killing curvatures for the excursion set $E_h(u)$, $u \in \mathbb{R}$ (see [1, 7, 8] for more details).

**Definition 2.3** (LK curvatures of $E_h(u)$). Let $h$ be a real-valued stationary field that is almost surely of class $C^2$. Define the following Lipschitz-Killing curvatures for the excursion set $T \cap E_h(u)$, $u \in \mathbb{R}$, $T$ bounded rectangle in $\mathbb{R}^2$,

$$L_2(h, u, T) := |T \cap E_h(u)|$$

$$L_1(h, u, T) := \frac{|\partial(T \cap E_h(u))|_1}{2},$$

$$L_0(h, u, T) := \sharp \text{ connected components in } T \cap E_h(u) - \sharp \text{ holes in } T \cap E_h(u).$$

Furthermore, the normalized LK curvatures are given by

$$C_i^{\prime T}(h, u) := \frac{L_i(h, u, T)}{|T|}, \quad \text{for } i = 0, 1, 2,$$

and the associated LK densities, if they exist, are

$$C_i^*(h, u) := \lim_{T \nearrow \mathbb{R}^2} \mathbb{E}[C_i^{\prime T}(h, u)], \quad \text{for } i = 0, 1, 2.$$

The existence of the above limits for $h$ a Gaussian field as in Definition 2.1 is given in Proposition 2.3 in [8].

Figure 1 displays a realization of a Gaussian random field (first row) and of the associated perturbed one (second row) and two excursion sets for these fields for $u = 0$ (center) and $u = 1$ (right). We chose here a Student distributed centered random variable $X$ with $\nu = 5$ degrees of freedom, i.e., $X \sim t(\nu = 5)$, and $g$ a Bergmann-Fock Gaussian field prescribed by its covariance function $r(s) = \sigma_g^2 e^{-\kappa^2 \|s\|^2}$. In Figure 1 one can appreciate a visual similarity between these images and in particular in terms of their excursion sets. Then it could be difficult to evaluate the perturbation behind the considered Gaussian model by looking exclusively at Figure 1.

This motivates the necessity of a statistical and numerical processing in order to measure the impact of the perturbation. The goal of the next sections will be to study the LK curvatures of the perturbed field in order to both quantify the discrepancy between these black-and-white images and evaluate the robustness with respect to a small perturbation of the considered geometrical characteristics of the excursion sets.

## 2.2. Mean LK curvatures of excursion sets of perturbed Gaussian model

Let $g$ be as in Definition 2.1. The Gaussian kinematic formula provides the mean LK curvatures of excursion sets of $g$ within a rectangle $T$ (see, e.g., [1], Thm. 13.2.1 or [2], Thm. 4.3.1), for $u \in \mathbb{R}$,

$$\mathbb{E}[C_2^{\prime T}(g, u)] = \Psi\left(\frac{u}{\sigma_g}\right); \qquad (2.2)$$

$$\mathbb{E}[C_1^{\prime T}(g, u)] = \Psi\left(\frac{u}{\sigma_g}\right) \frac{|\partial T|_1}{2|T|} - \frac{\sqrt{2\pi\lambda}}{4\sigma_g} \Psi'\left(\frac{u}{\sigma_g}\right); \qquad (2.3)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Gaussian random field and its perturbed counter-part as in Definitions 2.1 and 2.2 with covariance $r(s) = \sigma_g^2 e^{-\kappa^2 \|s\|^2}$, for $\sigma_g = 2$, $\kappa = 100/2^{10}$ in a domain of size $2^{10} \times 2^{10}$ pixels, with $\epsilon = 1$ and $X \sim t(\nu = 5)$. First row: A realization of Gaussian random field $g$ (left) and the two associated excursion sets for $u = 0$ (center) and $u = 1$ (right). Second row: The associated realization of a perturbed Gaussian random field $f$ (left) and two excursion sets for $u = 0$ (center) and $u = 1$ (right).</figcaption></figure>

$$ \mathbb{E}[C_0^{T}(g, u)] = \Psi\left(\frac{u}{\sigma_g}\right) \frac{1}{|T|} - \sqrt{\frac{\lambda}{2\pi \sigma_g^2}} \Psi'\left(\frac{u}{\sigma_g}\right) \frac{|\partial T|_1}{2|T|} + \frac{\lambda}{2\pi \sigma_g^2} \Psi''\left(\frac{u}{\sigma_g}\right), \quad (2.4) $$

where $\lambda$ is the second spectral moment of $g$ and $\Psi(x) := \frac{1}{\sqrt{2\pi}} \int_x^{+\infty} e^{-t^2/2} dt$ the Gaussian tail distribution with zero mean and unit variance. Then the LK densities of the considered Gaussian field exist and are given by

$$ C_2^*(g, u) = \Psi\left(\frac{u}{\sigma_g}\right), \quad C_1^*(g, u) = -\frac{\sqrt{2\pi\lambda}}{4\sigma_g} \Psi'\left(\frac{u}{\sigma_g}\right), \quad C_0^*(g, u) = \frac{\lambda}{2\pi\sigma_g^2} \Psi''\left(\frac{u}{\sigma_g}\right). \quad (2.5) $$

We now focus on the modified field $f$. In the next result, we state the existence of the LK densities for the considered perturbed Gaussian field.
**Proposition 2.4.** Let $u \in \mathbb{R}$, $T$ a bounded rectangle in $\mathbb{R}^2$ and $f$ as in Definition 2.2. Then the LK densities $C_i^*(f, u)$, $i \in \{0, 1, 2\}$ exist and it holds that

$$
\begin{align*}
\mathbb{E}[C_2^{\prime T}(f, u)] &= C_2^*(f, u), \\
\mathbb{E}[C_1^{\prime T}(f, u)] &= C_1^*(f, u) + \frac{1}{2} C_2^*(f, u) \frac{|\partial T|_1}{|T|}, \\
\mathbb{E}[C_0^{\prime T}(f, u)] &= C_0^*(f, u) + \frac{1}{\pi} C_1^*(f, u) \frac{|\partial T|_1}{|T|} + C_2^*(f, u) \frac{1}{|T|}.
\end{align*}
$$

*Proof of Proposition 2.4.* Notice that for $i \in \{0, 1, 2\}$, $L_i(f, u, T) = L_i(g, u - \epsilon X, T)$ and hence,

$$
\mathbb{E}[C_i^{\prime T}(f, u)] = \mathbb{E}[\mathbb{E}[C_i^{\prime T}(g, u - \epsilon X)|X]]. \quad (2.6)
$$

Let us start with the area. Using the Gaussian kinematic formula (2.2) and the associated density $C_2^*(g, u)$ in (2.5), one yields

$$
\mathbb{E}[C_2^{\prime T}(f, u)] = \mathbb{E}[\mathbb{E}[C_2^{\prime T}(g, u - \epsilon X)|X]] = \mathbb{E}[\Psi((u - \epsilon X)/\sigma_g)].
$$

Letting $T \nearrow \mathbb{R}^2$, it proves that $C_2^*(f, u)$ exists and equals $\mathbb{E}[\Psi((u - \epsilon X)/\sigma_g)]$. Moreover, for any fixed rectangle $T$, $\mathbb{E}[C_2^{\prime T}(f, u)] = C_2^*(f, u)$.

Concerning the perimeter, using the Gaussian kinematic formula (2.3) and the associated density $C_1^*(g, u)$ in (2.5), one yields

$$
\begin{align*}
\mathbb{E}[C_1^{\prime T}(f, u)] &= \mathbb{E}[\mathbb{E}[C_1^{\prime T}(g, u - \epsilon X)|X]] \\
&= -\frac{\sqrt{2\pi\lambda}}{4\sigma_g} \mathbb{E}[\Psi'((u - \epsilon X)/\sigma_g)] + \mathbb{E}[\Psi((u - \epsilon X)/\sigma_g)] \frac{|\partial T|_1}{2|T|}.
\end{align*}
$$

Letting $T \nearrow \mathbb{R}^2$, it proves that $C_1^*(f, u)$ exists and equals $-\frac{\sqrt{2\pi\lambda}}{4\sigma_g} \mathbb{E}[\Psi'((u - \epsilon X)/\sigma_g)]$. Moreover, for any fixed rectangle $T$,

$$
\mathbb{E}[C_1^{\prime T}(f, u)] = C_1^*(f, u) + \frac{1}{2} C_2^*(f, u) \frac{|\partial T|_1}{|T|},
$$

as announced in the proposition.

A similar procedure with Formula (2.4) and associated density $C_0^*(g, u)$ in (2.5) yields the desired formula for $\mathbb{E}[C_0^{\prime T}(f, u)]$ in the proposition. $\square$

Notice that Proposition 2.4 proves that $f$ in Definition 2.2 is a “standard random field” in the sense of Definition 2.1 in [8]. For this reason, the same unbiased estimates can be applied also in this setting.

Actually, as products of (2.2)–(2.3)–(2.4)–(2.5) and Proposition 2.4 respectively, one can build the following unbiased estimators of $C_i^*(h, u)$, $i = 0, 1, 2$ for $h$ equal $g$ or $f$.

$$
\hat{C}_{0,T}(h, u) = C_0^{\prime T}(h, u) - \frac{|\partial T|_1}{\pi |T|} C_1^{\prime T}(h, u) + \left( \frac{1}{2\pi} \left( \frac{|\partial T|_1}{|T|} \right)^2 - \frac{1}{|T|} \right) C_2^{\prime T}(h, u), \quad (2.7)
$$

$$
\hat{C}_{1,T}(h, u) = C_1^{\prime T}(h, u) - \frac{|\partial T|_1}{2|T|} C_2^{\prime T}(h, u), \qquad (2.8)
$$

$$
\hat{C}_{2,T}(h, u) = C_2^{\prime T}(h, u). \tag{2.9}
$$
In the forthcoming section, we will provide a numerical illustration of the performance of the unbiased estimators $\hat{C}_{i,T}(h, u)$ above, for both $h = g$ and $h = f$.

## 2.3. Asymptotics of the mean LK curvatures for small perturbations

As we consider our perturbation as an additive spatial-invariant noise, from now on we will assume that $\mathbb{E}[X] = 0$. In this particular setting, we now provide Taylor expansions of the LK curvatures for the perturbed Gaussian model for small $\epsilon$.

**Proposition 2.5** (LK curvatures for the perturbed Gaussian model). Let $f(t) = g(t) + \epsilon X$, $t \in \mathbb{R}^2$ as in Definition 2.2 and $\mathbb{E}[X] = 0$. Then, for small $\epsilon$, it holds that

$$ \mathbb{E}[C_0^{I/T}(f, u)] = C_0^*(g, u) \left( 1 + \frac{\epsilon^2 \mathbb{E}[X^2]}{2\sigma_g^2} \left( H_2 \left( \frac{u}{\sigma_g} \right) - 2 \right) \right) \quad (2.10) \\ + \frac{1}{\pi} C_1^*(g, u) \left( 1 + \frac{\epsilon^2 \mathbb{E}[X^2]}{2\sigma_g^2} H_2 \left( \frac{u}{\sigma_g} \right) \right) \frac{|\partial T|_1}{|T|} \\ + \left( C_2^*(g, u) + \epsilon^2 \mathbb{E}[X^2] \frac{\pi}{\lambda} C_0^*(g, u) \right) \frac{1}{|T|} + O \left( \epsilon^3 \left( 1 + \frac{|\partial T|_1}{2|T|} + \frac{1}{|T|} \right) \right), $$

$$ \mathbb{E}[C_1^{I/T}(f, u)] = C_1^*(g, u) + C_2^*(g, u) \frac{|\partial T|_1}{2|T|} \quad (2.11) \\ + \epsilon^2 \mathbb{E}[X^2] \left( \frac{C_1^*(g, u)}{2\sigma_g^2} H_2 \left( \frac{u}{\sigma_g} \right) + C_0^*(g, u) \frac{\pi}{\lambda} \frac{|\partial T|_1}{2|T|} \right) + O \left( \epsilon^3 \left( 1 + \frac{|\partial T|_1}{2|T|} \right) \right), $$

$$ \mathbb{E}[C_2^{I/T}(f, u)] = C_2^*(g, u) + \epsilon^2 \mathbb{E}[X^2] \frac{\pi}{\lambda} C_0^*(g, u) + O(\epsilon^3), \quad (2.12) $$

where $H_2(y) = y^2 - 1$, for $y \in \mathbb{R}$ (i.e., the second Hermite polynomial) and the constants involved in the $O$-notation only depend on $g$ and $X$.

*Proof of Proposition 2.5.* Let $a := \frac{\lambda}{\sigma_g^2}$. In the following we will use that $\Psi'(x) = -\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$, $\Psi''(x) = x\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$, $\Psi'''(x) = \Psi'(x)H_2(x)$ and $\Psi''''(x) = \Psi''(x)H_2(x) + 2x\Psi'(x) = \Psi''(x)(H_2(x) - 2)$. From (2.6), Taylor developing the Gaussian tail distribution $\Psi$ and bearing in mind that $X$ is a centered random variable we have

$$ \begin{aligned} \mathbb{E}[C_2^{I/T}(f, u)] &= \mathbb{E}[\mathbb{E}[C_2^{I/T}(g, u - \epsilon X)|X]] = \mathbb{E}[\Psi((u - \epsilon X)/\sigma_g)] \\ &= \Psi(u/\sigma_g) + \frac{\epsilon^2 \Psi''(u/\sigma_g)}{2\sigma_g^2} \mathbb{E}[X^2] + O\left(\frac{\epsilon^3}{\sigma_g^3} \mathbb{E}[|X|^3]\right), \end{aligned} \quad (2.13) $$

where the constant involved in the $O$-notation is universal, in particular it does not depend on $u$ since the third derivative of $\Psi$ is bounded. One can rewrite (2.13), by using the kinematic formula for LK densities $C_i^*(g, u)$ of the Gaussian field $g$ in (2.5). Hence the result in (2.12). Analogously,

$$ \begin{aligned} \mathbb{E}[C_1^{I/T}(f, u)] &= \mathbb{E}[\mathbb{E}[C_1^{I/T}(g, u - \epsilon X)|X]] \\ &= -\frac{\sqrt{a}}{4} \sqrt{2\pi} \mathbb{E}[\Psi'((u - \epsilon X)/\sigma_g)] + \mathbb{E}[\Psi((u - \epsilon X)/\sigma_g)] \frac{|\partial T|_1}{2|T|} \\ &= -\frac{\sqrt{a}}{4} \sqrt{2\pi} \left( \Psi'(u/\sigma_g) + \frac{\epsilon^2 \Psi'''(u/\sigma_g)}{2\sigma_g^2} \mathbb{E}[X^2] + O(\epsilon^3) \right) \\ &\quad + \left( C_2^*(g, u) + \epsilon^2 C_0^*(g, u) \frac{\pi \mathbb{E}[X^2]}{\lambda} + O(\epsilon^3) \right) \frac{|\partial T|_1}{2|T|}. \end{aligned} $$
Then by using the Gaussian LK densities $C_i^*(g, u)$ in (2.5), we get equation (2.11). Finally,

$$
\begin{align*}
\mathbb{E}[C_0^{T'}(f, u)] &= \mathbb{E}[\mathbb{E}[C_0^{T'}(g, u - \epsilon X)|X]] \\
&= \frac{a}{2\pi} \mathbb{E}[\Psi''((u - \epsilon X)/\sigma_g)] - \sqrt{\frac{a}{2\pi}} \mathbb{E}[\Psi'((u - \epsilon X)/\sigma_g)] \frac{|\partial T|_1}{2|T|} \\
&\quad + \mathbb{E}[\Psi(u - \epsilon X)/\sigma_g] \frac{1}{|T|} \\
&= \frac{a}{2\pi} \left( \Psi''(u/\sigma_g) + \frac{\Psi'''(u/\sigma_g)}{2} \frac{\epsilon^2}{\sigma_g^2} \mathbb{E}[X^2] + O(\epsilon^3) \right) \\
&\quad - \sqrt{\frac{a}{2\pi}} \left( \Psi'(u/\sigma_g) + \frac{\epsilon^2}{\sigma_g^2} \frac{\Psi'''(u/\sigma_g)}{2} \mathbb{E}[X^2] + O(\epsilon^3) \right) \frac{|\partial T|_1}{2|T|} \\
&\quad + \left( \Psi(u/\sigma_g) + \frac{\epsilon^2}{\sigma_g^2} \frac{\Psi''(u/\sigma_g)}{2} \mathbb{E}[X^2] + O(\epsilon^3) \right) \frac{1}{|T|}.
\end{align*}
$$

As before, by using (2.5), we get

$$
\begin{align*}
\mathbb{E}[C_0^{T'}(f, u)] &= C_0^*(g, u) \left( 1 + \epsilon^2 (H_2\left(\frac{u}{\sigma_g}\right) - 2) \frac{\mathbb{E}[X^2]}{2\sigma_g^2} \right) \\
&\quad + \frac{1}{\pi} C_1^*(g, u) \left( 1 + \epsilon^2 \frac{\mathbb{E}[X^2]}{2\sigma_g^2} H_2\left(\frac{u}{\sigma_g}\right) \right) \frac{|\partial T|_1}{|T|} \\
&\quad + \left( C_2^*(g, u) + \epsilon^2 \mathbb{E}[X^2] \frac{\pi}{\lambda} C_0^*(g, u) \right) \frac{1}{|T|} + O\left( \epsilon^3 \left( 1 + \frac{|\partial T|_1}{2|T|} + \frac{1}{|T|} \right) \right).
\end{align*}
$$

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

In the next result, we give close formulae for the LK densities for the considered perturbed Gaussian field.

**Corollary 2.6 (LK densities for the considered perturbed Gaussian model).** Under assumptions of Proposition 2.5 and using the same notations, it holds that

$$
C_0^*(f, u) = C_0^*(g, u) \left( 1 + \frac{\epsilon^2 \mathbb{E}[X^2]}{2\sigma_g^2} \left( H_2\left(\frac{u}{\sigma_g}\right) - 2 \right) \right) + O(\epsilon^3),
$$

$$
C_1^*(f, u) = C_1^*(g, u) \left( 1 + \frac{\epsilon^2 \mathbb{E}[X^2]}{2\sigma_g^2} H_2\left(\frac{u}{\sigma_g}\right) \right) + O(\epsilon^3),
$$

$$
C_2^*(f, u) = C_2^*(g, u) + \epsilon^2 \mathbb{E}[X^2] \frac{\pi}{\lambda} C_0^*(g, u) + O(\epsilon^3).
$$

The proof of Corollary 2.6 is based on the property of the VH-growing sequence of rectangles $T$ on $\mathbb{R}^2$.

**2.4. Numerical illustrations**

An illustration for the finite sample performance of the proposed three unbiased estimators $\hat{C}_{0,T}(g, u)$, $\hat{C}_{1,T}(g, u)$ and $\hat{C}_{2,T}(g, u)$ obtained by adapting equations (2.7)–(2.9) to $g$, is given in Figure 2. In this case, $g$ has
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Gaussian random field as in Definition 2.1 with covariance $r(s) = \sigma_g^2 e^{-\kappa^2 \|s\|^2}$, for $\sigma_g = 2$, $\kappa = 100/2^{10}$ in a domain of size $2^{10} \times 2^{10}$ pixels. Theoretical $u \mapsto C_0^*(g, u)$ (left panel), $C_1^*(g, u)$ (center panel) and $C_2^*(g, u)$ (right panel) in (2.5) are drawn in black lines. We display with red stars the averaged values on $M = 100$ sample simulations of $\hat{C}_{0,T}(g, u)$ (left panel), $\hat{C}_{1,T}(g, u)$ (center panel) $\hat{C}_{2,T}(g, u)$ (right panel), obtained by adapting to $g$ estimators in (2.7)–(2.9), as a function of the level $u$. We also present $u \mapsto C_0^{\prime T}(g, u)$ and $C_1^{\prime T}(g, u)$ in green dashed lines (left and center panels). These samples have been obtained with Matlab using circulant embedding matrix.</figcaption></figure>

a covariance function $r(s) = \sigma_g^2 e^{-\kappa^2 \|s\|^2}$. Analogously, a good statistical performance of $\hat{C}_{0,T}(f, u)$, $\hat{C}_{1,T}(f, u)$ and $\hat{C}_{2,T}(f, u)$ can be observed in Figure 3. Here, we observe small empirical intervals (in red vertical lines) showing an accurate precision in the estimations.

The quantities $C_0^{\prime T}$, $C_1^{\prime T}$ and $C_2^{\prime T}$ in (2.7)–(2.9) are computed with the Matlab functions `bweuler`, `bwperim` and `bwarea`, respectively. When it is required to specify the connectivity, we average between the 4 and 8 connectivity. Since $C_1^*$ is defined as the average half perimeter, we divide by 2 the output derived from `bwperim`.

The good quality of the estimates $\hat{C}_{i,T}$ (red stars) for $i = 0, 1, 2$ of LK densities $C_i^*$ (plain black lines) is globally evident in all panels of Figures 2 and 3 in particular if we compare with biased $C_i^{\prime T}$ (green dashed lines). The impact of this bias correction is clearly appreciable for small levels $u$ in center panels of Figures 2 and 3.

However, for values of $u$ near zero $\hat{C}_{1,T}$ overestimates $C_1^*$ in the three center panels of Figures 2 and 3. This problem, independent to the statistical procedure, is generated by a well known numerical phenomenon in the Matlab function `bwperim`, due to the discretized representation from a pixelated image of the continuous boundary of the excursion set (see, e.g., the discussion in [8]).

The functions $u \mapsto C_i^*(g, u)$, for $i = 0, 1, 2$ in blue dashed lines in Figure 3 could be used as reference values to visually appreciate the discrepancy between the considered geometrical characteristics of the excursion sets of the Gaussian model and the perturbed one (black plain lines for $C_i^*(f, u)$, red stars for $\hat{C}_{i,T}(f, u)$, for $i = 0, 1, 2$).

Conversely to Figure 1 where the quantification of the perturbation was hard to get, by providing this statistical and numerical analysis based on the LK curvatures, we are now able to precisely measure the impact of the perturbation. Furthermore, the contiguity of the Gaussian model $g$ with respect to the perturbed one $f$ can be observed in their LK curvatures when the magnitude of the perturbation decreases, i.e., $\epsilon \to 0$ (see Fig. 3: first row with $\epsilon = 0.4$, second row with $\epsilon = 1$).
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Perturbed Gaussian random field as in Definition 2.2 with covariance $r(s) = \sigma_g^2 e^{-\kappa^2 \|s\|^2}$, for $\sigma_g = 2$, $\kappa = 100/2^{10}$ in a domain of size $2^{10} \times 2^{10}$ pixels, with $X \sim t(\nu = 5)$ and $\epsilon = 0.4$ (first row); $\epsilon = 1$ (second row). Theoretical $u \mapsto C_0^*(f, u)$ (left panel), $C_1^*(f, u)$ (centered panel) and $C_2^*(f, u)$ (right panel) in Corollary 2.6 are drawn in black plain lines and relative $u \mapsto C_0^*(g, u)$, $C_1^*(g, u)$ and $C_2^*(g, u)$ in blue dashed lines. We present $u \mapsto C_0^{\prime T}(f, u)$ and $C_1^{\prime T}(f, u)$ in green dotted lines (left and center panels). We also display the averaged values on $M = 100$ sample simulations of $\hat{C}_{0,T}(f, u)$ (left panel), $\hat{C}_{1,T}(f, u)$ (centered panel) $\hat{C}_{2,T}(f, u)$ (right panel) in (2.7)–(2.9), as a function of the level $u$ by using red stars. The empirical intervals associated to the estimation of $\hat{C}_{i,T}(f, u)$, for $i = 0, 1, 2$ are given by using red vertical lines. These samples have been obtained with Matlab using circulant embedding matrix.</figcaption></figure>

### 3. ASYMPTOTICS FOR THE EXCURSION AREA OF PERTURBED GAUSSIAN FIELDS

Recall that $C_2^{\prime T}(f, u) = C_2^{\prime T}(g, u - \epsilon X)$ and that $\mathbb{E}[C_2^{\prime T}(g, u)] = C_2^*(g, u) = \Psi(u/\sigma_g)$. We are interested in the asymptotic distribution as $T \uparrow \mathbb{R}^2$ of

$$Y_T^\epsilon(u) := |T|^{1/2} (C_2^{\prime T}(f, u) - \mathbb{E}[C_2^{\prime T}(f, u)]). \quad (3.1)$$
Considering the unperturbed case, by Theorem 3 in [10] for instance, we know that for a Gaussian field $g$ as in Definition 2.1 and for any $u \in \mathbb{R}$, the following convergence in distribution holds,

$$|T|^{1/2} (C_2^{\prime T}(g, u) - C_2^*(g, u)) \xrightarrow[T \nearrow \mathbb{R}^2]{} \mathcal{N}(0, v(u)), \qquad (3.2)$$

with

$$v(u) = \frac{1}{2\pi} \int_{\mathbb{R}^2} \int_0^{\rho(t)} \frac{1}{\sqrt{1-x^2}} \exp \left\{ -\frac{u^2}{\sigma_g^2(1+x)} \right\} dx dt \qquad (3.3)$$

and $\rho(t) := \text{corr}(g(0), g(t)) = r(t)/\sigma_g^2$. The interested reader is also referred to [22, 27].

Actually, we are able to state a more powerful result. It is given in the next lemma where the convergence in (3.2) is proved to be uniform with respect to level $u$. In order to formulate our result, let us introduce the Wasserstein distance between random variables $Z_1$ and $Z_2$:

$$d_W(Z_1, Z_2) = \sup_{h \in \mathcal{H}} |\mathbb{E}[h(Z_1)] - \mathbb{E}[h(Z_2)]|,$$

where $\mathcal{H}$ denotes the set of Lipschitz functions whose Lipschitz constant is $\le 1$.

**Lemma 3.1.** Let $g$ be a Gaussian field as in Definition 2.1. Then,

$$d_W(|T|^{1/2} (C_2^{\prime T}(g, u) - C_2^*(g, u)), \mathcal{N}(0, v(u))) = O((\log |T|)^{-1/12}),$$

where $v(u)$ is defined in (3.3) and the $O$-constant does not depend on the level $u$.

*Proof of Lemma 3.1.* We apply Lemma A.1 with $X = g/\sigma_g$. Indeed, the condition (2.1) implies assumption (A.1). Hence, with the notations that are in force, in particular (A.2), conclusion of Lemma A.1 can be rewritten as

$$d_W(|T|^{1/2} (C_2^{\prime T}(g, u) - C_2^*(g, u)), \mathcal{N}(0, \sigma^2(u/\sigma_g))) = O((\log |T|)^{-1/12}), \qquad (3.4)$$

Note that obviously, (3.4) yields the same CLT as (3.2) for the excursion area as $T \nearrow \mathbb{R}^2$, and hence the variance $\sigma^2(u/\sigma_g)$ given by (A.2) equals variance $v(u)$ given by (3.3). $\square$

Let us come back to the study of the asymptotics of $Y_T^\epsilon$, defined by (3.1). We will use the next decomposition

$$\begin{aligned} Y_T^\epsilon(u) &= |T|^{1/2} \left(C_2^{\prime T}(f, u) - \Psi\left(\frac{u - \epsilon X}{\sigma_g}\right)\right) + |T|^{1/2} \left(\Psi\left(\frac{u - \epsilon X}{\sigma_g}\right) - \mathbb{E}\left[\Psi\left(\frac{u - \epsilon X}{\sigma_g}\right)\right]\right) \\ &=: Z_T^\epsilon(u) + R_T^\epsilon(u). \end{aligned} \qquad (3.5)$$

### 3.1. Asymptotics for fixed small $\epsilon$ and $T \nearrow \mathbb{R}^2$

In this section, we introduce a non-Gaussian random variable that we denote by $\Theta_\epsilon(u)$. We firstly provide an upper-bound for the Wasserstein distance between $Z_T^\epsilon(u)$ in (3.5) and $\Theta_\epsilon(u)$. Secondly, we describe the form of the density of $\Theta_\epsilon(u)$ by providing a Taylor expansion for small $\epsilon > 0$.

**Theorem 3.2** (Quantitative asymptotics for $Z_T^\epsilon(u)$). Let $f(t) = g(t) + \epsilon X$, $t \in \mathbb{R}^2$ as in Definition 2.2 with $\mathbb{E}[X] = 0$. For any fixed $\epsilon > 0$ and $u \in \mathbb{R}$, we consider $\Theta_\epsilon(u)$ a random variable whose conditional distribution
given {$X = x$} is centered Gaussian with variance $v(u - \epsilon x)$, $v(\cdot)$ being defined by (3.3). Then, as $T \nearrow \mathbb{R}^2$, it
holds that

$$
d_W(Z_T^\epsilon(u), \Theta_\epsilon(u)) = O\left( (\log|T|)^{-1/12} \right),
$$

where the constant involved in the O-notation depends neither on $\epsilon$ nor on $u$.

Proof of Theorem 3.2. By the definition of Wasserstein distance, we have

$$
\begin{align*}
d_W(Z_T^\epsilon(u), \Theta_\epsilon(u)) &= \sup_{h \in \mathcal{H}} |\mathbb{E}[h(Z_T^\epsilon(u)) - h(\Theta_\epsilon(u))]| \\
&\le \mathbb{E} \left[ \sup_{h \in \mathcal{H}} |\mathbb{E}[h(Z_T^\epsilon(u)) - h(\Theta_\epsilon(u))|X]| \right].
\end{align*}
$$

The latter supremum is equal to the Wasserstein distance between $Z_T^\epsilon(u)$ and $\Theta_\epsilon(u)$ with respect to the conditional expectation given $X$.

Actually, conditionally to {$X = x$}, $Z_T^\epsilon(u)$ equals $|T|^{1/2} (C_2^{1/T}(g, u - \epsilon x) - C_2^*(g, u - \epsilon x))$, since $C_2^*(g, \cdot) := \Psi(\frac{\cdot}{\sigma_g})$ and $\Theta_\epsilon(u)$ is $\mathcal{N}(0, v(u - \epsilon x))$ distributed.

Hence, applying Lemma 3.1 yields $\sup_{h \in \mathcal{H}} |\mathbb{E}[h(Z_T^\epsilon(u)) - h(\Theta_\epsilon(u))|X]| = O((\log|T|)^{-1/12})$, where the O-
constant does not depend on $u$ nor on $\epsilon$ and $X$. Lebesgue dominated convergence theorem allows us to
conclude. $\square$

We now focus on the random variable $\Theta_{\epsilon}(u)$ that has been introduced in Theorem 3.2. Let us quote that it is non-Gaussian, yielding an unusual non-Gaussian limit of $Z_{T}^{\epsilon}(u)$ as $T \nearrow \mathbb{R}^{2}$. In the next theorem, we provide the density distribution function of $\Theta_{\epsilon}(u)$ and a corresponding Taylor expansion for small $\epsilon > 0$.

**Theorem 3.3.** Under the same assumptions as Theorem 3.2, it holds that, for fixed $\epsilon > 0$,

$$
Z_T^\epsilon(u) = |T|^{1/2} \left( C_2^{1/T} (f, u) - \Psi \left( \frac{u - \epsilon X}{\sigma_g} \right) \right) \xrightarrow{d} \Theta_\epsilon(u), \quad (3.6)
$$

where $\Theta_{\epsilon}(u)$'s probability density function is given by

$$
h_{\epsilon}: y \mapsto \mathbb{E}[\phi(v(u - \epsilon X), y)], \quad y \in \mathbb{R}, \tag{3.7}
$$

where $\phi(v, \cdot)$ stands for the p.d.f. of $\mathcal{N}(0, v)$ and $v(\cdot)$ is given by (3.3). Furthermore, $h_\epsilon$ can be expanded for small $\epsilon > 0$, as

$$
\begin{equation}
\begin{split}
h_{\epsilon}(y) = f_{\text{BEP}}^{\delta=0}(y) &+ \frac{\epsilon^2 \mathbb{E}[X^2]}{2} \left[ \frac{3}{4} \frac{v'(u)^2}{v(u)^2} \left( f_{\text{BEP}}^{\delta=0}(y) - 2f_{\text{BEP}}^{\delta=2}(y) + f_{\text{BEP}}^{\delta=4}(y) \right) \right. \\
& \left. \qquad + \frac{1}{2} \frac{v''(u)}{v(u)} \left( -f_{\text{BEP}}^{\delta=0}(y) + f_{\text{BEP}}^{\delta=2}(y) \right) \right] + O(\epsilon^3),
\end{split}
\tag{3.8}
\end{equation}
$$

where

$$
f_{\text{BEP}}^{\delta}(y) = \left| \frac{y}{\sqrt{2v(u)}} \right|^{\delta} e^{-\frac{y^2}{2v(u)}} \left( \sqrt{2v(u)} \Gamma\left(\frac{\delta+1}{2}\right) \right)^{-1}, \quad y \in \mathbb{R}. \quad (3.9)
$$

The notation $f_{\text{BEP}}$ in (3.9) stands for *Bimodal Exponential Power* density function as detailed in the discussion below.
**Proof of Theorem 3.3.** The convergence in (3.6) is a direct consequence of Theorem 3.2. In order to get the probability density function of $\Theta_{\epsilon}(u)$ in (3.7), it is enough to compute $\mathbb{E}[\xi(\Theta_{\epsilon}(u))]$ for any bounded positive function $\xi$ as follows,

$$ \begin{aligned} \mathbb{E}[\xi(\Theta_{\epsilon}(u))] &= \mathbb{E}[\mathbb{E}[\xi(\Theta_{\epsilon}(u))|X]] = \mathbb{E}\left[\int_{\mathbb{R}} \xi(y) \phi(v(u - \epsilon X), y) \, dy\right] \\ &= \int_{\mathbb{R}} \xi(y) \mathbb{E}[\phi(v(u - \epsilon X), y)] \, dy, \end{aligned} $$

where Fubini-Tonelli theorem has been used for the last equality.

To get the approximation of $h_{\epsilon}$ in (3.8), we recall the following result that can be proved with Taylor expansion and easy algebra.

**Lemma 3.4.** For any function $\varphi$ in $C^3(\mathbb{R})$ with bounded derivatives up to order two and any random variable $\eta$ with finite third moment,

$$ \mathbb{E}[\varphi(\eta)] = \varphi(\mathbb{E}\eta) + \frac{1}{2}\varphi''(\mathbb{E}\eta) \mathrm{Var}\eta + O(\mathbb{E}[|\eta - \mathbb{E}\eta|^3]), $$

where the constant in O-notation depends on $\mathrm{Var}\eta$ and on the bounds of derivatives of $\varphi$.

Applying Lemma 3.4 with $\eta = \epsilon X$ and $\varphi(\cdot) = \phi(v(u - \cdot), y)$ for fixed $u, y$ and $\epsilon$, and bearing in mind that $\mathbb{E}[X] = 0$, one gets

$$ \mathbb{E}[\phi(v(u - \epsilon X), y)] = \phi(v(u), y) + \frac{1}{2}\varphi''(0) \epsilon^2 \mathbb{E}[X^2] + O(\epsilon^3 \mathbb{E}[|X|^3]), $$

where

$$ \varphi''(0) = \partial_{vv}^2 \phi(v(u), y) v'(u)^2 + \partial_v \phi(v(u), y) v''(u). $$

Since

$$ \partial_v \phi(v, y) = \frac{\sqrt{\pi v} e^{-\frac{y^2}{2v}} (y^2 - v)}{2\pi \sqrt{2} v^3} \quad \text{and} \quad \partial_{vv}^2 \phi(v, y) = \frac{\sqrt{\pi v} (y^4 + 3v^2 - 6vy^2) e^{-\frac{y^2}{2v}}}{4v^5 \pi \sqrt{2}}, $$

the proof is complete.

**Discussion on $h_{\epsilon}$ density.**

Since the density in (3.8) plays a crucial role in our asymptotics and as its non-Gaussian shape was not previously studied in the literature, in the following we propose an analysis of the truncated version of $h_{\epsilon}(y)$ (i.e., we neglect the $O(\epsilon^3)$ term),

$$ \tilde{h}_{\epsilon}(y) = f_{\text{BEP}}^{\delta=0}(y)(1 + \gamma_1 - \gamma_2) + f_{\text{BEP}}^{\delta=2}(y)(\gamma_2 - 2\gamma_1) + f_{\text{BEP}}^{\delta=4}(y)\gamma_1, \quad (3.10) $$

where $f_{\text{BEP}}^{\delta}(y)$ as in (3.9) and $\gamma_1 := \frac{3}{8} \epsilon^2 \mathbb{E}[X^2] \frac{v'(u)^2}{v(u)^2}$, $\gamma_2 := \frac{1}{4} \epsilon^2 \mathbb{E}[X^2] \frac{v''(u)}{v(u)}$.

Firstly, one can remark that coefficients $\gamma_1$ and $\gamma_2$ of the linear combination $\tilde{h}_{\epsilon}$ depend on the variance function $v(u)$ in (3.3) and on its first and second derivatives. For the nodal set with $u = 0$ one can easily evaluate $v(0) = (2\pi)^{-1} \int_{\mathbb{R}^2} \arcsin(\rho(t)) dt$. An illustration of theoretical $u \mapsto v(u)$ in (3.3), $u \mapsto v'(u)$ and $u \mapsto v''(u)$ can be found in Figure 4 (left panel).
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Preliminary studies to build $\tilde{h}_{\epsilon}$ as in (3.10), on BEP functions, $u \mapsto v(u)$ and its derivatives. Left panel: $u \mapsto v(u)$ in (3.3) (dotted black line), $u \mapsto v'(u)$ (dashed red line) and $u \mapsto v''(u)$ (crossed blue line), for several levels $u$. Black point represents the value $v(0)$. Center and right panels: theoretical Bimodal Exponential Power density functions for $u = 1.5$ (center panel) and $u = 3$ (right panel). We display $f_{\text{BEP}}^{\delta=0}$ (black line) $f_{\text{BEP}}^{\delta=2}$ (red line) $f_{\text{BEP}}^{\delta=4}$ (blue line), with $f_{\text{BEP}}^{\delta}$ as in (3.9) with $v(u)$ as in (3.3). The considered correlation function is $\rho(t) = e^{-\kappa^2 \|t\|^2}$, with $\sigma_g = 1$ and $\kappa = 100/2^{10}$.</figcaption></figure>

Furthermore, notice that function $f_{\text{BEP}}^{\delta}(y)$ in (3.9) is a particular case of the Bimodal Exponential Power density function, i.e., $f_{\text{BEP}}(y) = (\alpha |\frac{y-\mu}{\zeta}|^{\delta} e^{-|\frac{y-\mu}{\zeta}|^{\alpha}})(2\zeta\Gamma(\frac{\delta+1}{\alpha}))^{-1}$, for $y \in \mathbb{R}$ (see [17]) with fixed values of parameters $\alpha = 2, \mu = 0$ and $\zeta = \sqrt{2v(u)}$ and varying $\delta$. Obviously, for $y \in \mathbb{R}$, $f_{\text{BEP}}^{\delta=0}(y) = \phi(v(u), y)$, i.e., the Gaussian density with zero mean and variance $v(u)$. An illustration of the behaviour of these multimodal densities for two different values of $u$ is given in Figure 4 ($u = 1.5$ center panel, $u = 3$ right panel).

Theoretical resulting $\tilde{h}_{\epsilon}$ functions in (3.10), built by using $v(u)$, $v'(u)$, $v''(u)$ and $f_{\text{BEP}}^{\delta}$ functions studied above, are displayed below in Section 3.3 (see Figs. 5 and 6).

## 3.2. Asymptotics for $\epsilon \to 0$ and $T \nearrow \mathbb{R}^2$

Let $T^{(N)} = NT$, as introduced in Section 2.1. In the following we prove that $Y_{T^{(N)}}^{\epsilon_N}$ given by (3.1) satisfies a classical Central Limit Theorem as soon as $\epsilon_N$ goes to 0 sufficiently fast, for $N \to \infty$.

**Theorem 3.5.** Let $f(t) = g(t) + \epsilon X$, $t \in \mathbb{R}^2$ as in Definition 2.2 with $\mathbb{E}[X] = 0$ and $\epsilon_N$ be such that

$$ \lim_{N \to \infty} N \epsilon_N = 0. \qquad (3.11) $$

Then it holds that,

$$ Y_{T^{(N)}}^{\epsilon_N}(u) = |T^{(N)}|^{1/2} \left( C_2^{T^{(N)}}(f, u) - \mathbb{E}[C_2^{T^{(N)}}(f, u)] \right) \xrightarrow{d} \mathcal{N}(0, v(u)), $$

with $v(u)$ given by (3.3).

**Proof of Theorem 3.5.** We start by writing $Y_T^{\epsilon}(u) = (Y_T^{\epsilon}(u) - \Theta_{\epsilon}(u)) + \Theta_{\epsilon}(u)$.

On the one hand, by triangular inequality we have

$$ d_W(Y_T^{\epsilon}(u), \Theta_{\epsilon}(u)) \le d_W(Y_T^{\epsilon}(u), Z_T^{\epsilon}(u)) + d_W(Z_T^{\epsilon}(u), \Theta_{\epsilon}(u)) \qquad (3.12) $$
From (3.5), we have $d_W(Y_T^\epsilon(u), Z_T^\epsilon(u)) \le \sqrt{\mathbb{E}[R_T^\epsilon(u)^2]}$. Then, since

$$ \mathbb{E}[R_T^\epsilon(u)^2] = |T| \mathbb{E} \left[ \left( \Psi \left( \frac{u - \epsilon X}{\sigma_g} \right) - \mathbb{E} \left[ \Psi \left( \frac{u - \epsilon X}{\sigma_g} \right) \right] \right)^2 \right], $$

by performing the Taylor expansion of $\Psi(\frac{u-\epsilon X}{\sigma_g})$ at order 1 and using $|\Psi'| \le C$, we get

$$ \mathbb{E}[R_T^\epsilon(u)^2] = |T| \left( \Psi' \left( \frac{u}{\sigma_g} \right)^2 \frac{\epsilon^2 \mathbb{E}[X^2]}{\sigma_g^2} + O(\epsilon^3) \right) = \epsilon^2 |T| \left( \frac{C^2}{\sigma_g^2} \mathbb{E}[X^2] + O(\epsilon) \right), $$

where the constant involved in the O-notation does not depend on $T$. Then, from condition in (3.11), the first term on the r.h.s. of (3.12) with $\epsilon = \epsilon_N$ and $T = T^{(N)}$ goes to 0 as $N$ goes to infinity.

Concerning the second term, Theorem 3.2 yields $\kappa_1 (\log |T|)^{-1/12}$ as upper bound, where $\kappa_1$ does not depend on $\epsilon$. Therefore, the second term on the r.h.s. of (3.12) goes to 0 as $T \nearrow \mathbb{R}^2$ uniformly with respect to $\epsilon$ (see Thm. 3.3).

Finally, thanks to the Wasserstein distance in (3.12) that goes to 0, we get that $Y_{T^{(N)}}^{\epsilon_N}(u) - \Theta_{\epsilon_N}(u)$ converges to 0 in distribution.

On the other hand, $\Theta_{\epsilon_N}(u) \xrightarrow[N \to \infty]{} \mathcal{N}(0, v(u))$ since $h_\epsilon(y) \to \phi(v(u), y)$ as $\epsilon \to 0$. At last, Slutsky theorem allows us to conclude. $\square$

**Remark 3.6.** A careful inspection<fnref n="1" /> of the proof of Theorem 3.5 shows that under the assumption $1 + \rho(t) \ne 0$ for all $t \ne 0$ we have

$$ d_W(Y_{T^{(N)}}^{\epsilon_N}(u), \mathcal{N}(0, v(u))) \xrightarrow[N \to \infty]{} 0. $$

Actually,

$$ d_W(Y_{T^{(N)}}^{\epsilon_N}(u), \mathcal{N}(0, v(u))) \le d_W(Y_{T^{(N)}}^{\epsilon_N}(u), \Theta_{\epsilon_N}(u)) + d_W(\Theta_{\epsilon_N}(u), \mathcal{N}(0, v(u))). \quad (3.13) $$

The first term on the r.h.s. of (3.13) has been already shown to vanish as $N \to +\infty$ in the proof of Theorem 3.5. To deal with the second term, note first that for every $\epsilon > 0$

$$ \Theta_{\epsilon}(u) = \sqrt{v(u - \epsilon X)} U, $$

where $U$ is a standard Gaussian independent of $X$ and the equality is in law. Then

$$ \begin{aligned} d_W(\Theta_{\epsilon}(u), \mathcal{N}(0, v(u))) &= d_W(\sqrt{v(u - \epsilon X)}U, \sqrt{v(u)}U) \\ &\le \sqrt{\frac{2}{\pi v(u)}} \mathbb{E}[|v(u - \epsilon X) - v(u)|] \\ &\le C \left( \frac{1 + |u|}{\sqrt{v(u)}} \right) \epsilon, \end{aligned} \quad (3.14) $$

<footnote n="1">¹We thank an anonymous referee for this insight.</footnote>
<table><caption>TABLE 1. Parameters setting associated to Figure 5. Here $|T| = 1024^2$, $\sigma_g = 1$ and $\rho(t) = e^{-\kappa^2\|t\|^2}$, for $\kappa = 100/2^{10}$, i.e., $\lambda = 0.019$.</caption><thead><tr><th>u</th><th>X</th><th>ε</th><th>ε := ε<sup>2</sup> E[X<sup>2</sup>]</th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>Figure 5</th></tr></thead><tbody><tr><td rowspan="3">1.5</td><td rowspan="3">Skellam<sub>μ<sub>1</sub>=μ<sub>2</sub>=1</sub></td><td>0.5</td><td>0.5</td><td>0.979</td><td>0.485</td><td>left panel</td></tr><tr><td>0.3</td><td>0.18</td><td>0.352</td><td>0.175</td><td>center panel</td></tr><tr><td>0.1</td><td>0.02</td><td>0.039</td><td>0.019</td><td>right panel</td></tr><tr><td rowspan="3">3</td><td rowspan="3">Skellam<sub>μ<sub>1</sub>=μ<sub>2</sub>=1</sub></td><td>0.5</td><td>0.5</td><td>2.818</td><td>1.774</td><td>left panel</td></tr><tr><td>0.3</td><td>0.18</td><td>1.015</td><td>0.639</td><td>center panel</td></tr><tr><td>0.1</td><td>0.02</td><td>0.113</td><td>0.071</td><td>right panel</td></tr></tbody></table>

for some positive constant $C$ only depending on $E[X^2]$ and $\sigma_g^2$. Indeed, if $1 + \rho(t) \neq 0$ for all $t \neq 0$, from (3.3) we have

$$|v'(u)| \leq \frac{2}{\pi} |u| \int_{\mathbb{R}^2} \frac{|\rho(t)|}{\sqrt{1 + \rho(t)}(\sqrt{1 + \rho(t)} + \sqrt{1 - \rho(t)})} dt \leq C|u|,$$

for some positive constant $C$ depending only on $\sigma_g^2$. In particular, the inequality in (3.14) ensures the second term on the r.h.s. of (3.13) to go to zero as $N \to +\infty$.

### 3.3. Numerical illustrations

All over this section, $\sigma_g$ is assumed to be equal to 1. In the following, by using histograms we compare the empirical density of the random variable $Z_T^\epsilon(u) := |T|^{1/2}(C_2^{T/T}(f, u) - \Psi(u - \epsilon X))$ versus the truncated probability density function of $\Theta_\epsilon$, i.e., $\tilde{h}_\epsilon$ given in (3.10). Each histogram is built by reproducing 300 Monte-Carlo independent simulations in a large domain such that $|T| = 1024^2$.

#### Case 1: X is Skellam distributed.

Firstly, we consider the case where $X$ follows a discrete Skellam probability distribution which is the difference of two independent Poisson-distributed random variables with respective expected values $\mu_1$ and $\mu_2$. In particular we consider $\mu_1 = \mu_2$ to have $X$ centered. We choose the parameters setting gathered in Table 1. Obtained results are shown Figure 5 for $u = 1.5$ (first row) and for $u = 3$ (second row). Furthermore, necessary preliminary studies to built $\tilde{h}_\epsilon$ as in (3.10), on BEP functions, $u \mapsto v(u)$ and its derivatives are given in Section 3.1.

#### Case 2: X is t-distributed.

We now consider the case where $X$ follows a t-distribution and the parameters are those in Table 2. Obtained results are shown Figure 6 for $u = 1.5$ (first row) and for $u = 3$ (second row). Preliminary studies of BEP functions, $u \mapsto v(u)$ and its derivatives are identical to those in Section 3.1.

The multimodal behaviour of $\tilde{h}_\epsilon$ in (3.10) is clearly visible in Figures 5 and 6. Furthermore in the numerical studies above one can appreciate the contiguity property of the proposed model for $\epsilon \to 0$. Indeed since theoretically $h_\epsilon(y) \to \phi(v(u), y)$ as $\epsilon \to 0$, in Figures 5 and 6 the unimodal Gaussian behaviour appears when the perturbation magnitude decreases ($\epsilon = 0.5$ in first columns of Figures 5 and 6, $\epsilon = 0.3$ in second columns and $\epsilon = 0.1$ in third ones). Finally the choice of the level $u$ plays an important role in term of magnitude of obtained histograms (see the y-axis scale in Figs. 5 and 6). This behaviour was already visible in the theoretical $\tilde{h}_\epsilon$ function (see center and right panels of Fig. 4 for $u = 1.5$ and $u = 3$ respectively).
<figure><img src="image_7.png" /><figcaption>FIGURE 5. Histogram for the study of density of $Z_T^\epsilon$ when X is Skellam distributed, for $u = 1.5$ (first row) and $u = 3$ (second row), based on 300 Monte-Carlo independent simulations. The chosen parameters setting is gathered in Table 1. In particular we chose $\epsilon = 0.5$ (first column), $\epsilon = 0.3$ (second column) and $\epsilon = 0.1$ (third column). Necessary preliminary studies to build $\tilde{h}_\epsilon$ as in (3.10), on BEP functions, $u \mapsto v(u)$ and its derivatives are given in Figure 4. Resulting theoretical $\tilde{h}_\epsilon$ density is drawn by using red plain line.</figcaption></figure>

<table><caption>TABLE 2. Parameters setting associated to Figure 6. Here $|T| = 1024^2$, $\sigma_g = 1$ and $\rho(t) = e^{-\kappa^2\|t\|^2}$, for $\kappa = 100/2^{10}$, i.e., $\lambda = 0.019$.</caption><thead><tr><th>u</th><th>X</th><th>ε</th><th>ε := ε<sup>2</sup> E[X<sup>2</sup>]</th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>Figure 6</th></tr></thead><tbody><tr><td rowspan="3">1.5</td><td rowspan="3">t<sub>ν=5</sub></td><td>0.5</td><td>0.417</td><td>0.816</td><td>0.404</td><td>left panel</td></tr><tr><td>0.3</td><td>0.150</td><td>0.294</td><td>0.145</td><td>center panel</td><td rowspan="2">first row</td></tr><tr><td>0.1</td><td>0.017</td><td>0.033</td><td>0.016</td><td>right panel</td></tr><tr><td rowspan="3">3</td><td rowspan="3">t<sub>ν=5</sub></td><td>0.5</td><td>0.417</td><td>2.349</td><td>1.478</td><td>left panel</td></tr><tr><td>0.3</td><td>0.150</td><td>0.846</td><td>0.532</td><td>center panel</td><td rowspan="2">second row</td></tr><tr><td>0.1</td><td>0.017</td><td>0.094</td><td>0.059</td><td>right panel</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Histogram for the study of density of $Z_T^\epsilon$ when $X$ is $t$-distributed, for $u = 1.5$ (first row) and $u = 3$ (second row), based on 300 Monte-Carlo independent simulations. The chosen parameters setting is gathered in Table 2. In particular we chose $\epsilon = 0.5$ (first column), $\epsilon = 0.3$ (second column) and $\epsilon = 0.1$ (third column). Necessary preliminary studies to build $\tilde{h}_\epsilon$ as in (3.10), on BEP functions, $u \mapsto v(u)$ and its derivatives are given in Figure 4. Resulting theoretical $\tilde{h}_\epsilon$ density is drawn by using red plain line.</figcaption></figure>

# 4. INFERENCE FOR PERTURBATION

## 4.1. Unbiased estimator of the perturbation

In this section we will focus on the restrictive case $\sigma_g = 1$. If the variance of the considered field is unknown a preliminary estimation of $\sigma_g^2$ is necessary. In the case where only the sparse observations given by the excursion sets are available, this statistical issue is investigated in [14].

Let $u \neq 0$ being fixed. We introduce $\epsilon := \epsilon^2 \mathbb{E}[X^2]$. Since $\mathbb{E}[X] = 0$, it is clear that $\epsilon$ quantifies the variability around zero of the considered perturbation and it can be useful to measure the discrepancy between the observed excursion set $T \cap E_f(u)$ and the associated Gaussian one.

By using (2.12) and then (2.5), we can rewrite

$$ \mathbb{E}[C_2^{T}(f, u)] = C_2^*(g, u) + \epsilon \frac{\pi}{\lambda} C_0^*(g, u) + O(\epsilon^{3/2}) = \Psi(u) + \epsilon \frac{u}{2\sqrt{2\pi}} e^{-\frac{u^2}{2}} + O(\epsilon^{3/2}). $$
It appears then clearly that $\epsilon$ has the same order of magnitude than

$$ \epsilon_u := \frac{2\sqrt{2\pi} e^{\frac{u^2}{2}}}{u} (C_2^*(f, u) - \Psi(u)). \quad (4.1) $$

This means that $\epsilon$ in (4.1) can be estimated by using the LK curvature of order 2, i.e., the area of the excursion set at a (chosen) level $u$. Then, $\epsilon_u$ is completely empirically accessible by using this sparse observation because it does not depend on the (unknown) second spectral moment $\lambda$ of the Gaussian field. In Proposition 4.1 below, we present a consistent estimator based on the observation $T \cap E_f(u)$ for the perturbation error $\epsilon_u$.

**Proposition 4.1.** Let $f(t) = g(t) + \epsilon X$, $t \in \mathbb{R}^2$ as in Definition 2.2 with $\mathbb{E}[X] = 0$. Let $u \neq 0$ being fixed. Let consider the empirical counterpart of $\epsilon_u$ in (4.1), i.e.,

$$ \hat{\epsilon}_{u,T} := \frac{2\sqrt{2\pi}e^{\frac{u^2}{2}}}{u} (\hat{C}_{2,T}(f, u) - \Psi(u)), \quad (4.2) $$

with $\hat{C}_{2,T}(f, u)$ as in (2.9). Then, it holds that

(i) $\hat{\epsilon}_{u,T}$ is an unbiased estimator for $\epsilon_u$,

(ii) $|T^{(N)}|^{1/2}(\hat{\epsilon}_{u,T^{(N)}} - \epsilon_u) \xrightarrow{d} \mathcal{N}(0, \sigma_{\epsilon_u}^2)$, with $\sigma_{\epsilon_u}^2 = 8\pi \frac{e^{u^2}}{u^2} v(u)$ for $v(u)$ as in (3.3) with $\sigma_g = 1$.

*Proof of Proposition 4.1.* Since $\hat{C}_{2,T}(f, u)$ is an unbiased estimator of $C_2^*(f, u)$, one can easily see that $\mathbb{E}[\hat{\epsilon}_u - \epsilon_u] = 0$. Furthermore, using the fact that $|T^{(N)}|^{1/2}(\hat{\epsilon}_u - \epsilon_u) = \frac{\sqrt{8\pi}e^{\frac{u^2}{2}}}{u} Y_{T^{(N)}}(u)$, from Theorem 3.5 we get the result. $\square$

**Remark 4.2.** If $u = 0$, assuming that $\mathbb{E}[X^3] \neq 0$ and the fourth moment of $X$ is finite, then by Taylor developing the function $\mathbb{E}[C_2^{\prime T}(f, 0)]$ up to the order 3 (see Prop. 2.5) we easily get an unbiased and asymptotically normal estimator for $\epsilon^3 \mathbb{E}[X^3]$, similar to the r.h.s. of (4.2).

## 4.2. Numerical illustrations

In this section we provide an illustration of the inference procedure for the perturbation $\epsilon := \epsilon^2 \mathbb{E}[X^2]$ proposed in Section 4.1 above. The considered perturbed model and the associated parameters are gathered in Table 3. By using this framework, in Figure 7 one can appreciate the finite sample performance of the inference procedure proposed in Section 4.1 above, for several values of perturbation $\epsilon$ and for two levels $u$ ($u = 1.5$ in center panel and $u = 3$ in right one).

Unsurprisingly, we remark that the variability of the estimation is related to the choice of level $u$. The asymptotic standard deviation function $u \mapsto \sigma_{\epsilon_u}|T^{(N)}|^{-1/2}$ in the left panel of Figure 7 allows us to identify some choices of levels $u$ where the variance is minimum. Indeed, for large values of $|u|$, less observations are available than for intermediate values of $|u|$. This aspect can be appreciated by observing the larger confidence intervals in the case $u = 3$. For $u = 0$, the variance $\sigma_{\epsilon_u}^2$ diverges (see the left panel of Fig. 7) implying that this inference procedure will be not robust for $u \approx 0$ (see Rem. 4.2).

## 5. DISCUSSION

In this paper we study how geometric quantities change under small perturbations of the underlying Gaussian field. We focus on the particular case of an independent, additive, spatial-invariant noise with zero mean perturbing a stationary isotropic Gaussian field. This perturbed model, $f := g + \epsilon X$, can be seen as a small super-imposed noise to a Gaussian field. In this setting, the zero mean assumption of the noise is usually required in the literature. However similar Taylor expansion formulae as those in Proposition 2.5 and Corollary 2.6 are still possible with an extra linear $\epsilon$ term.
<table><caption>TABLE 3. Parameters setting associated to Figure 7. Furthermore we consider $|T| = 1024^2$, $\sigma_g = 1$ and $\rho(t) = e^{-\kappa^2\|t\|^2}$, for $\kappa = 100/2^{10}$, i.e., $\lambda = 0.019$.</caption><thead><tr><th>Level u</th><th>X</th><th>Chosen ε</th><th>ε := ε<sup>2</sup> E[X<sup>2</sup>]</th><th>average of estimated ε̂<sub>u</sub></th><th>Figure 7<br>on 100 Monte-Carlo Simulations</th></tr></thead><tbody><tr><td rowspan="5">1.5</td><td rowspan="5">Skellam<br>μ<sub>1</sub> = μ<sub>2</sub> = 1</td><td>0.1</td><td>0.02</td><td>0.023</td><td rowspan="5">first panel</td></tr><tr><td>0.2</td><td>0.08</td><td>0.085</td></tr><tr><td>0.3</td><td>0.18</td><td>0.182</td></tr><tr><td>0.4</td><td>0.32</td><td>0.324</td></tr><tr><td>0.5</td><td>0.50</td><td>0.492</td></tr><tr><td rowspan="5">3</td><td rowspan="5">Skellam<br>μ<sub>1</sub> = μ<sub>2</sub> = 1</td><td>0.1</td><td>0.02</td><td>0.033</td><td rowspan="5">second panel</td></tr><tr><td>0.2</td><td>0.08</td><td>0.072</td></tr><tr><td>0.3</td><td>0.18</td><td>0.158</td></tr><tr><td>0.4</td><td>0.32</td><td>0.345</td></tr><tr><td>0.5</td><td>0.50</td><td>0.510</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>FIGURE 7. Left panel: theoretical $u \mapsto \sigma_{\varepsilon_u} |T^{(N)}|^{-1/2}$ in Proposition 4.1. Center and right panels: Theoretical $\varepsilon := \varepsilon^2 \mathbb{E}[X^2]$ (black circles) and obtained average of $\hat{\varepsilon}_u$ on 100 Monte-Carlo independent simulations (red stars). Theoretical confidence intervals at level 0.95 prescribed by Proposition 4.1 are also displayed. The chosen parameters setting is gathered in Table 3. Here $u = 1.5$ (center panel) and $u = 3$ (right panel).</figcaption></figure>

Notice that we work here under the assumption $\mu = \mathbb{E}(g(0)) = 0$. Furthermore, the statistical inference procedure for perturbation in Section 4 is developed in the restrictive unit-variance case, $\sigma_g = \text{Var}(g(0)) = 1$. If the mean and the variance of the considered underlying Gaussian field are unknown a preliminary estimation is then necessary. In the case where only the sparse observations obtained from one or two excursion sets of the field are available, this important statistical issue has been recently investigated in [14].

The asymptotic analysis for the excursion area presented in this paper is based on a quantitative CLT in Wasserstein distance for $C_2(g, u)$, where in particular the bound does not depend on the level $u$ (see Lem. (3.1)). This result is a refinement of Theorem 2 in [28]. Concerning the other Lipschitz-Killing curvatures on $\mathbb{R}^2$, boundary length and Euler-Poincaré characteristic, it should be possible to get similar results, i.e. quantitative
CLT in Wasserstein distance with uniform bounds, carefully exploiting their chaotic decomposition. However, this interesting issue goes far beyond the purpose of the present work.

Finally, notice that the mean of LK curvatures in Proposition 2.5 can be derived also in the case of an additive spatially variant perturbation, i.e., $f(t) = g(t) + \epsilon X(t)$, for $t \in \mathbb{R}^2$ and $\epsilon > 0$, with $X$ a stationary random field with finite third moment and independent of $g$. The proof comes down in the same way and the results are completely analogous to those in equations (2.10), (2.11) and (2.12). However, the asymptotics results obtained in Section 3 will become more challenging in that case. Indeed, even in the classical case of excursion area of Gaussian fields, to the best of our knowledge we are not aware of any (quantitative) central limit theorem in the case of a non-constant level. This could represent an interesting point to investigate in a future work. For sake of completeness, the interested reader is referred to [21] where CLT results are obtained for the curve-crossings number of a stationary Gaussian process ($d = 1$) according to the form of the moving curve (periodic or linear).

## APPENDIX A. UNIFORM RATES OF CONVERGENCE FOR CLTS OF SOJOURN TIMES OF STATIONARY GAUSSIAN FIELDS

The following lemma is a refinement of a result in [28]. Therefore, we keep the notations introduced therein.

**Lemma A.1.** Let $\{X(t), t \in \mathbb{R}^d\}$ be a stationary centered Gaussian field with unit variance and covariance function $\rho \in L^1(\mathbb{R}^d)$. For $T > 0$ and $u \in \mathbb{R}$, we define $S_T(u)$ to be

$$S_T(u) := \int_{[0,T]^d} 1_{(X(t) \ge u)} dt$$

the excursion volume above level $u$. Under the hypothesis that

$$\int_{\mathbb{R}^d \setminus [-a, a]^d} |\rho(t)| dt = O(1/\log a), \quad a \to +\infty, \qquad (A.1)$$

we have, as $T \to +\infty$,

$$d_W \left( \frac{S_T(u) - T^d \Psi(u)}{\sqrt{T^d}}, \mathcal{N}(0, \sigma^2(u)) \right) = O \left( 1/(\log T)^{1/12} \right),$$

where the constant involved in the $O$-notation only depends on the field $\{X(t), t \in \mathbb{R}^d\}$ and where

$$0 < \sigma^2(u) := \sum_{n=1}^{+\infty} \frac{\phi^2(u) H_{n-1}^2(u)}{n!} \int_{\mathbb{R}^d} \rho^n(t) dt < +\infty, \qquad (A.2)$$

$\phi$ being the density function of the standard Gaussian and $\Psi$ the tail of its distribution.

**Remark A.2.** Actually, Theorem 2 in [28] ensures that, as $|T| \to +\infty$,

$$d_W \left( \frac{S_T(u) - T^d \Psi(u)}{\sqrt{T^d}}, \mathcal{N}(0, \sigma^2(u)) \right) = O \left( 1/(\log T)^{1/4} \right), \qquad (A.3)$$

where the constant involved in the $O$-notation depends on the field and the level. By adapting the proof, we provide a uniform rate of convergence w.r.t. the level $u$ in (A.3).
Proof of Lemma A.1. We will use the following estimate see e.g. ([19], (30) and [20], Prop. 3): there exists $K > 0$ such that, for every $n \in \mathbb{N}$ and $x \in \mathbb{R}$, we have

$$ \exp(-x^2/4)|H_n(x)| \le K \sqrt{n!} n^{-1/12}. \qquad (A.4) $$

We can write

$$ d_W \left( \frac{S_T(u) - T^d \Psi(u)}{\sqrt{T^d}}, \mathcal{N}(0, \sigma^2(u)) \right) \le d_W \left( \frac{S_T(u) - T^d \Psi(u)}{\sqrt{T^d}}, \frac{S_{T,N_T}(u) - \mathbb{E}[S_{T,N_T}(u)]}{\sqrt{T^d}} \right) \\ + d_W \left( \frac{S_{T,N_T}(u) - \mathbb{E}[S_{T,N_T}(u)]}{\sqrt{T^d}}, \mathcal{N}(0, \sigma_{N_T}^2(u)) \right) \\ + d_W(\mathcal{N}(0, \sigma^2(u)), \mathcal{N}(0, \sigma_{N_T}^2(u))) =: d_1 + d_2 + d_3, $$

where $S_{T,N_T}$ is the truncation of $S_T$ at position $N_T$ in the Wiener chaos expansion ($N_T$ will be chosen later on). For $d_1$ we have, due to (A.4),

$$ \begin{aligned} d_1 &\le \sqrt{\sum_{n=N_T+1}^{+\infty} \frac{\phi^2(u) H_{n-1}^2(u)}{n! T^d} \int_{[-T,T]^d} \rho^n(t) \prod_{j=1}^d (T-|t_j|) dt} \\ &\le K \sqrt{\phi(u)} \sqrt{\int_{\mathbb{R}^d} |\rho(t)| dt} \sqrt{\sum_{n=N_T+1}^{+\infty} \frac{1}{n(n-1)^{1/6}}} \\ &= O(N_T^{-1/12}), \end{aligned} \qquad (A.5) $$

where the constant involved in the $O$-notation only depends on the field.

Note that due to (A.4) we can give upper bounds for $d_2$ and $d_3$ (being inspired by the proof of Thm. 2 in [28]) independently of $u$,

$$ d_2 = O(3^{N_T}/\sqrt{T^d}), \qquad d_3 = O_{T \to +\infty}(N_T^{-1/12} + T^{-1/4} + (\log T)^{-1/2}). \qquad (A.6) $$

Summing up the bounds for $d_1, d_2, d_3$ in (A.5) and (A.6), and choosing $N_T = (\log T)/4$ we have

$$ d_1 + d_2 + d_3 = O_{T \to +\infty}((\log T)^{-1/12}), $$

that concludes the proof.

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

Acknowledgements. The authors would like to thank Domenico Marinucci for insightful comments. The research of MR has been supported by the Fondation Sciences Mathématiques de Paris, and is currently supported by the ANR project *Unirandom* (ANR-17-CE40-0008) and by the PRA 2018 49 project at the University of Pisa. The research of EDB and AE has been partially supported by the ANR project *MISTIC* (ANR-19-CE40-0005). The authors sincerely express their gratitude to two anonymous referees for their valuable comments and remarks which improve the quality of the present work.

REFERENCES

[1] R.J. Adler and J.E. Taylor, Random fields and geometry. Springer Monographs in Mathematics. Springer, New York (2007).

[2] R.J. Adler and J.E. Taylor, Topological complexity of smooth random functions. *Lectures from the 39th Probability Summer School held in Saint-Flour, 2009, École d'Été de Probabilités de Saint-Flour. [Saint-Flour Probability Summer School].* Vol. 2019 of Lecture Notes in Mathematics. Springer, Heidelberg (2011).
[3] R.J. Adler, G. Samorodnitsky and J.E. Taylor, Excursion sets of three classes of stable random fields. *Adv. Appl. Prob.* **42** (2010) 293–318.

[4] D. Beliaev, M. McAuley and S. Muirhead, On the number of excursion sets of planar Gaussian fields. Preprint to 1807.10209 (2018).

[5] C. Berzin, Estimation of Local Anisotropy Based on Level Sets. Preprint to 1801.03760 (2018).

[6] T.H. Beeman, A.M. Turner and V. Vitelli, Stochastic geometry and topology of non-Gaussian fields. *Proc. Natl. Acad. Sci.* **109** (2012) 19943–19948.

[7] H. Biermé and A. Desolneux, On the perimeter of excursion sets of shot noise random fields. *Ann. Probab.* **44** (2016) 521–543.

[8] H. Biermé, E. Di Bernardino, C. Duval and A. Estrade, Lipschitz-Killing curvatures of excursion sets for two-dimensional random fields. *Electr. J. Stat.* **13** (2019) 536–581.

[9] F. Bron and D. Jeulin, Modelling a food microstructure by random sets. *Image Anal. Stereol.* **23** (2011).

[10] A. Bulinski, E. Spodarev and F. Timmermann, Central limit theorems for the excursion set volumes of weakly dependent random fields. *Bernoulli* **18** (2012) 100–118.

[11] A.E. Burgess, Mammographic structure: Data preparation and spatial statistics analysis. *Medical Imaging'99. International Society for Optics and Photonics* (1999) 642–653.

[12] E.M. Cabaña, Affine processes: a test of isotropy based on level sets. *SIAM J. Appl. Math.* **47** (1987) 886–891.

[13] B. Casaponsa, B. Crill, L. Colombo, L. Danese, J. Bock, A. Catalano, A. Bonaldi, S. Basak, L. Bonavera, A. Coulais et al., Planck 2015 results: XVI. Isotropy and statistics of the CMB. *Astron. Astrophys.* **594** (2016) A16.

[14] E. Di Bernardino and C. Duval, Statistics for Gaussian Random Fields with Unknown Location and Scale using Lipschitz-Killing Curvatures. Preprint hal-02317747 (2019).

[15] A. Estrade and J.R. León, A central limit theorem for the Euler characteristic of a Gaussian excursion set. *Ann. Probab.* **44** (2016) 3849–3878.

[16] Y. Fantaye, D. Marinucci, F. Hansen and D. Maino, Applications of the Gaussian kinematic formula to CMB data analysis. *Phys. Rev. D* **91** (2015) 063501.

[17] M.Y. Hassan and R.H. Hijazi, A bimodal exponential power distribution. *Pak. J. Stat.* **26** (2010) 379–396.

[18] C. Hikage and T. Matsubara, Limits on second-order non-gaussianity from Minkowski functionals of WMAP data. *Mon. Not. R. Astron. Soc.* **425** (2012) 2187–2196.

[19] E. Hille, A class of reciprocal functions. *Ann. Math.* **27** (1926) 427–464.

[20] P. Imkeller, V. Pérez-Abreu and J. Vives, Chaos expansions of double intersection local time of Brownian motion in $R^d$ and renormalization. *Stoch. Process. Appl.* **56** (1995) 1–34.

[21] M. Kratz and J. León, Level curves crossings and applications for Gaussian models. *Extremes* **13** (2010) 315–351.

[22] M. Kratz and S. Vadlamani, Central limit theorem for Lipschitz-Killing curvatures of excursion sets of Gaussian random fields. *J. Theor. Probab.* **31** (2018) 1729–1758.

[23] R. Lachièze-Rey, Normal convergence of nonlocalised geometric functionals and shot-noise excursions. *Ann. Appl. Probab.* **29** (2019) 2613–2653.

[24] D. Marinucci and G. Peccati, Random fields on the sphere. Representation, limit theorems and cosmological applications. In Vol. 389 of *London Mathematical Society Lecture Note Series*. Cambridge University Press, Cambridge (2011).

[25] D. Marinucci and M. Rossi, Stein-Malliavin approximations for nonlinear functionals of random eigenfunctions on $S^d$. *J. Funct. Anal.* **268** (2015) 2379–2420.

[26] T. Matsubara, Analytic minkowski functionals of the cosmic microwave background: second-order non-gaussianity with bispectrum and trispectrum. *Phys. Rev. D* **81** (2010) 083505.

[27] D. Müller, A central limit theorem for Lipschitz-Killing curvatures of Gaussian excursions. *J. Math. Anal. Appl.* **452** (2017) 1040–1081.

[28] V.-H. Pham, On the rate of convergence for central limit theorems of sojourn times of Gaussian fields. *Stoch. Process. Appl.* **123** (2013) 2158–2174.

[29] A. Roberts and M. Teubner, Transport properties of heterogeneous materials derived from Gaussian random fields: bounds and simulation. *Phys. Rev. E* **51** (1995) 4141–4154.

[30] A. Roberts and S. Torquato, Chord-distribution functions of three-dimensional random media: approximate first-passage times of Gaussian processes. *Phys. Rev. E* **59** (1999) 4953–4963.

[31] R. Schneider and W. Weil, Stochastic and integral geometry. Probability and its Applications. Springer-Verlag, Berlin (2008).

[32] E. Spodarev, Limit theorems for excursion sets of stationary random fields. In *Modern stochastics and applications*, vol. 90 of Springer Optim. Appl. Springer, Cham (2014) 221–241.

[33] C. Thäle, 50 years sets with positive reach – a survey. *Surv. Math. Appl.* **3** (2008) 123–165.

[34] K.J. Worsley, The geometry of random images. *Chance* **1** (1997) 27–40.