# Existence, properties, and parametric inference for possibly hyperuniform Gibbs perturbed lattices

Jean-François Coeurjolly, Christopher Renaud-Chan

► To cite this version:

Jean-François Coeurjolly, Christopher Renaud-Chan. Existence, properties, and parametric inference for possibly hyperuniform Gibbs perturbed lattices. 2026. hal-05526996

HAL Id: hal-05526996

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

Preprint submitted on 25 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
EXISTENCE, PROPERTIES, AND PARAMETRIC INFERENCE FOR
POSSIBLY HYPERUNIFORM GIBBS PERTURBED LATTICES

BY JEAN-FRANÇOIS COEURJOLLY <fnref n="1" />,ª AND CHRISTOPHER RENAUD CHAN<fnref n="1" />,ᵇ

This work lies at the intersection of Gibbs models and hyperuniform point processes. Classical Gibbs models, whether defined on lattices or in continuous space, provide flexible tools to describe interacting particle systems but are generally not hyperuniform. Conversely, known hyperuniform models such as the Ginibre process or perturbed lattices lack flexibility and typically cannot enforce physically relevant constraints such as hard-core interactions. We introduce a new class of models, termed Gibbs perturbed lattice models, which preserve a lattice structure while allowing interactions through a Hamiltonian defined on the perturbed particle locations. We establish existence results for the associated Gibbs measures, derive DLR-type equilibrium equations, and show that some models in this class exhibit hyperuniformity. Finally, we propose statistical inference methods based on the Takacs-Fiksel type approach and prove their asymptotic properties.

# 1. Introduction.

**Gibbs models.** Gibbs processes form a particular class of stochastic models designed to describe interacting geometric objects (e.g., points). Such processes are typically defined either on a lattice, say $S = \mathbb{Z}^d$, or on a continuous space, say $S = \mathbb{R}^d$. General references include [18, 19, 20, 31]. More precisely, when $S = \mathbb{Z}^d$, one speaks of a *Gibbs lattice system*, denoted by $\mathbf{X} = \{X_i, i \in \mathbb{Z}^d\}$. The locations $i$ are fixed, and the model takes values in $E^{\mathbb{Z}^d}$, where $E$ denotes the state space. The model is specified through conditional distributions and a Hamiltonian depending on the variables $X_i$. Well-known examples include the Ising model, where $E = \{-1, 1\}$, the Potts model, where $E = \{a_0, \dots, a_{K-1}\}$ with $K \ge 2$ (see, e.g., [17]), and Besag autoregressive-type models [4]. When $S = \mathbb{R}^d$, one speaks of a *Gibbs point process*, denoted by $\mathbf{X} = \{x_1, \dots, x_n, \dots\}$ with $x_i \in \mathbb{R}^d$. For standard Gibbs point processes, both the locations and the number of points are random, and the model can again be described through conditional densities with respect to the Poisson measure, see, for example, [9]. Interactions between points are encoded via an energy function, or Hamiltonian. When this Hamiltonian depends on pairwise distances, one obtains pairwise interaction point processes, including the classical Strauss model [32] and the Lennard-Jones model (see, e.g., [31]). Models with higher-order interactions include, for instance, the Quer-mass model [24], which encompasses the area-interaction and Widom-Rowlinson models as special cases, and the Delaunay-Voronoi interaction models [13]. Gibbs lattice systems and Gibbs point processes constitute two fundamentally different classes of models, each with its own characteristics, mathematical tools, and theoretical results. In particular, the existence and phase transition analysis for infinite-volume Gibbs point processes relies on techniques that differ from those used in the lattice setting [9, 10, 19].

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, LJK, 38000 Grenoble, France, ajean-francois.coeurjolly@univ-grenoble-alpes.fr;</footnote>

<footnote>ᵇchristopher.renaud-chan@univ-grenoble-alpes.fr</footnote>

<footnote>MSC2020 subject classifications: Primary 60G55; Secondary 62M30.</footnote>

<footnote>Keywords and phrases: Lattice systems; Gibbs processes; Hyperuniformity; Dobrushin-Lanford-Ruelle equations; Simulation; Takacs-Fiksel method.</footnote>
**Hyperuniformity.** Hyperuniformity is a relatively recent concept in probability theory. Introduced in the physics literature, most notably by [36] (see also the surveys by [34, 35]). It aims to describe systems of rigid particles exhibiting an unusual form of large-scale order. Roughly speaking, a point process is said to be hyperuniform if $\text{Var}[N(B(0, r))]/|B(0, r)|$, that is the normalized variance of the number of points falling in a Euclidean ball with radius $r$, tends to 0 as $r \to \infty$. This condition expresses that the fluctuations of the counting variable grow more slowly than those of a homogeneous Poisson point process. We refer the reader to [26] for a recent and comprehensive probabilistic and statistical overview of hyperuniform point processes. Hyperuniformity is not simply related to repulsiveness. For instance, Matérn hard-core point processes are known not to be hyperuniform (see, e.g., [25]). Rather, hyperuniformity is a global and subtle property.

Several classes of point processes are known to be hyperuniform, including the Sine-$\beta$ process [37], the Ginibre process, as particular examples of (continuous) determinantal point processes (see, e.g., [27]), and the zero set of Gaussian Analytic Functions [21]. Another important class consists of perturbed lattices, defined as $\Gamma = \{i+X_i, i \in L\}$, where $L$ is a lattice and $(X_i)_i$ is a random field indexed by $L$. When the random variables $X_i$, called *moves*, are i.i.d. with common distribution $Q$, [12] provide conditions on $Q$ and on the dimension $d$ that ensure hyperuniformity of $\Gamma$. Although the Ginibre and perturbed lattice models are simple to define and interpret, they suffer from important limitations. In particular, it is not possible to impose a hard-core condition of the form $|i + X_i - j - X_j| \ge \delta > 0$ for $i \ne j$. Such a constraint is, however, of significant interest in physics [5]. One might conjecture that this issue could be addressed through a dependent subsampling of a hyperuniform point process; however, such a procedure destroys the hyperuniformity property [26].

**Gibbs perturbed lattice model.** At this point, one may wonder whether continuous Gibbs point processes, which are natural and flexible models able to produce repulsive point patterns, can be hyperuniform. The answer is in fact negative, as shown recently by [11]. The authors prove that most Gibbs point processes defined on $\mathbb{R}^d$ are not hyperuniform and fail to identify any hyperuniform example within this class. In this paper, we introduce a new class of models lying at the intersection of perturbed lattices and continuous Gibbs point processes. On a finite sub-lattice, $\Delta \subset L$, and for ease of exposition, a Gibbs perturbed lattice model is defined through a density with respect to $Q^{\otimes \Delta}$. It is worth emphasizing that our framework fundamentally differs from both Gibbs lattice systems and continuous Gibbs point processes. Indeed, the resulting point process remains of the form $\Gamma = \{i + X_i, i \in L\}$. Each lattice site corresponds to a unique point of $\Gamma$, so that $\Gamma$ can also be interpreted as a lattice random field. The key difference, however, is that the Hamiltonian does not depend on the values of the random field itself, but instead on the locations of the perturbed points, that is, on the elements of $\Gamma$.

**Contributions.** Given the specific structure of our model, important theoretical questions, such as the existence of infinite-volume Gibbs measures and the formulation of Dobrushin-Lanford-Ruelle (DLR) equations [9], cannot be directly addressed using existing results for classical Gibbs lattice systems or continuous Gibbs point processes (see, e.g., [19, 20, 31]). We tackle these issues in the present paper. We provide sufficient conditions on the move distribution $Q$ and on the Hamiltonian of the particle system that ensure the existence of stationary Gibbs perturbed lattice models, both in finite and infinite volumes. Building on a general result obtained by [12], we then show that some models within this class are hyperuniform. In addition, we introduce a broad family of equilibrium equations (see Section 3.4), referred to as first-order, second-order, and variational DLR equations, which can be viewed as analogues of the corresponding first-order, second-order, and variational Georgii-Nguyen-Zessin (GNZ) equations [2, 19, 29].
Having at our disposal a broad new class of parametric models, we then address the problem of estimating a parameter $\theta \in \mathbb{R}^p$ ($p \ge 1$) from a single realization of a Gibbs perturbed lattice model. We consider two observation settings, summarized in Figure 1: either the pairs $(i, i + X_i)$ are observed, or only the perturbed locations $(i + X_i)$. Due to the presence of an intractable normalizing constant in the model definition, maximum likelihood estimation is impractical, and likelihood-based approximations are computationally expensive, as is also the case for standard Gibbs lattice systems and continuous Gibbs point processes, see e.g. [18]. We therefore propose alternative inference methods based on pseudo-likelihood [3, 4, 23] and variational approaches (see, e.g., [1, 2]). We show that both approaches can be formulated as particular cases of the Takacs-Fiksel method, originally introduced for continuous Gibbs point processes via Georgii-Nguyen-Zessin (GNZ) equations [6, 15, 16, 33]. In our framework, the Takacs-Fiksel estimator is defined by minimizing a contrast function based on the $L^2$ distance between empirical versions of DLR equations evaluated on $L \ge p$ test functions. It is worth emphasizing that the pseudo-likelihood and Takacs-Fiksel contrast functions differ substantially from their counterparts for Gibbs lattice systems and continuous Gibbs point processes, as they require nonstandard and careful border corrections involving both $i$ and $i+x$ for $x \in \mathbb{R}^d$. We work within an increasing-domain asymptotic framework and establish, using in particular novel ergodic results for empirical versions of DLR equations, that the resulting Takacs-Fiksel estimator of $\theta$ is consistent. Moreover, we derive its rate of convergence.

**Outline.** The remainder of the paper is organized as follows. Section 2 reviews Gibbs lattice systems and Gibbs point processes and introduces the main assumptions imposed on the move distribution $\mathcal{Q}$ and on the Hamiltonian of the particle system considered throughout the paper. Section 3 is devoted to the definition of a Gibbs perturbed lattice models, the existence of the associated Gibbs measures, their main properties, and the derivation of DLR-type equations. Section 4 addresses parametric statistical inference for the Gibbs perturbed lattice model. We present several inference methodologies and establish asymptotic results within the increasing-domain framework. Section 5 provides illustrative examples of the proposed model and statistical methodology. Finally, the proofs of the main results are deferred to A-F.

## 2. Background and notation.

### 2.1. Lattice, perturbed lattices and point processes.
Let $\mathcal{L}$ be a full rank lattice, subgroup of $\mathbb{R}^d$ generated by linear combinations with integer coefficients of the vectors $(a_1, \dots, a_d)$ that forms a basis of $\mathbb{R}^d$. The fundamental domain associated to a lattice $\mathcal{L}$ is the set

$$ \mathcal{D} := \left\{ \sum_{i=1}^{d} t_i a_i, (t_1, \dots, t_d) \in [0, 1)^d \right\}. $$

The fundamental domain is the polyhedra that is used to pave $\mathbb{R}^d$ according to the lattice. We denote by $\delta := \min\{|i|, i \in \mathcal{L}\}$ the size parameter of the lattice. The basic example is $\mathcal{L} = \mathbb{Z}^d$ and $\mathcal{D} = [0, 1)^d$. In this paper, the lattice $\mathcal{L}$ is viewed as the origin of the particles. We consider $\mathbb{R}^d$ to be the mark space, also called spin space or state space, which describes the displacement of a single particle from their lattice positions. We denote by $\mathcal{Q}$ a probability measure on $\mathbb{R}^d$, which serves as the reference distribution for this single displacement particle. For each $i \in \mathcal{L}$, we denote by $x_i \in \mathbb{R}^d$ the displacement (or perturbation) of the particle at site $i$. The lattice system configuration space is denoted by $\Omega := (\mathbb{R}^d)^{\mathcal{L}}$, and is endowed with the product $\sigma$-algebra, and a configuration is denoted by $\mathbf{x} := (x_i)_{i \in \mathcal{L}}$. For any subset $\Lambda \subset \mathcal{L}$, the configuration on the sub-lattice is denoted by $\Omega_{\Lambda} := (\mathbb{R}^d)^{\Lambda}$, and the projection of a configuration $\mathbf{x} \in \Omega$ into $\Omega_{\Lambda}$ is denoted by $\mathbf{x}_{\Lambda} := (x_i)_{i \in \Lambda}$. A random field on the lattice
$\mathcal{L}$, denoted by $\mathbf{X}$, is therefore a random variable on $\Omega$. In the following, we let $U \sim \mathcal{U}_D$ be a random variable uniformly distributed on $\mathcal{D}$. For $B \subseteq \mathbb{R}^d$, we let $\mathbb{1}_B : \mathbb{R}^d \to \{0, 1\}$ denote the indicator function on $B$. Finally, for any vector $z \in \mathbb{R}^k$, $k \ge 1$ or any bounded set $B \subset \mathbb{R}^d$, the notation $|z|$ and $|B| = \lambda^d(B)$ respectively stand for the Euclidean norm of $z$ and the volume (the $d$-dimensional Lebesgue measure $\lambda^d$) of $B$.

A perturbed lattice, denoted by $\Gamma$, is a random variable from $\Omega \times \mathbb{R}^d$ to the set of all countable subsets of $\mathbb{R}^d$ such that $\Gamma(\mathbf{X}, U) := \{i + X_i + U, i \in \mathcal{L}\}$, where we require that particles arising from a perturbed lattice are indistinguishable. Furthermore, for any subset $\Lambda \subset \mathcal{L}$, the perturbed sub-lattice $\Gamma_\Lambda$ is the random variable $\Gamma_\Lambda(\mathbf{X}_\Lambda, U) := \{i + X_i + U, i \in \Lambda\}$. We denote by $\mathcal{L}_U := \{i + U, i \in \mathcal{L}\}$ the shifted lattice. This shifted lattice is by construction invariant under any translation on $\mathbb{R}^d$. Therefore, $U$ stationarizes the underlying lattice from which we build the perturbed lattice. From now on, unless explicitly specified and when there is no ambiguity, the notation $P, E[\cdot], \mathrm{Var}[\cdot], \mathrm{Cov}[\cdot, \cdot]$, standing for probability, expectation, variance and covariance, is used without indicating the underlying probability measure(s). A point configuration in $\mathbb{R}^d$ is in general denoted by $\gamma$ and is a locally finite set in $\mathbb{R}^d$, i.e. for any bounded Borel set $W$ of $\mathbb{R}^d$, $\#\gamma \cap W$ is finite. The set of point configurations is denoted by $\mathcal{C}$ and the set of finite point configurations is denoted by $\mathcal{C}_f$. We endow $\mathcal{C}$ with the $\sigma$-algebra generated by the counting function $N_W$ for all Borel sets $W$ in $\mathbb{R}^d$. In the following, for any $W \subset \mathbb{R}^d$ and $\gamma \in \mathcal{C}$ we denote by $\gamma_W := \gamma \cap W$, the restriction of $\gamma$ to $W$. A point process is a random variable on $\mathcal{C}$. For $\gamma \in \mathcal{C}$ and $v \in \mathbb{R}^d$, we denote by $\gamma + v$ the translation of all points of $\gamma$ by $v$. A point process $\Gamma$ is said to be stationary when its distribution is translation invariant, i.e. the distributions of $\Gamma + v$ and $\Gamma$ are identical for any $v \in \mathbb{R}^d$.

**2.2. Hamiltonian of a system and move distribution.** The Gibbs perturbed lattice model depends mainly on two ingredients: the Hamiltonian $H$ which measures the interaction between particles of the system and the mark (or move) distribution $\mathcal{Q}$. In this section, we provide definitions, examples and the main assumptions on $H$ and $\mathcal{Q}$.

**Hamiltonian or energy functional.** We consider models of interaction that are common in the canonical and grand canonical Gibbs point process setting. The Hamiltonian, also called the energy functional, is a measurable function, $H : \mathcal{C}_f \to \mathbb{R} \cup \{\infty\}$. Let $W \subset \mathbb{R}^d$ be a bounded domain, the local energy on $W$ of a configuration $\zeta \in \mathcal{C}$ is defined as

$$ (2.1) \qquad H_W(\gamma) := \lim_{n \to \infty} \left( H(\gamma_{[-n,n]^d}) - H(\gamma_{[-n,n]^d \setminus W}) \right) $$

if the limit exists, where we use the convention $\infty - \infty = 0$. Another quantity of interest is the local energy of a point which measures the contribution of a particle at a certain position to the energy of the whole configuration. Let $\zeta \in \mathcal{C}$ and $x \in \mathbb{R}^d$, the local energy is a measurable function from $\mathbb{R}^d \times \mathcal{C}$ to $\mathbb{R} \cup \{\infty\}$ given by

$$ (2.2) \qquad h(x, \gamma) := \lim_{n \to \infty} \left( H(\{x\} \cup \gamma_{[-n,n]^d}) - H(\gamma_{[-n,n]^d}) \right). $$

We impose a few definitions, restrictions on $H$. These assumptions are in particular required to ensure the existence of a Gibbs perturbed lattice in the infinite volume, see Theorem 1.

(H1) $H$ is stable, i.e. there exists $A \ge 0$ such that for all $\gamma \in \mathcal{C}_f$, $H(\gamma) \ge -AN(\gamma)$.

(H2) $H$ is non-degenerate, i.e. there exists $\varepsilon > 0$ and $B \ge 0$ such that for all $i \in \mathcal{L}$, $|x_i| \le \varepsilon$, and for any finite $\Lambda \subset \mathcal{L}$, $H(\{i + x_i\}_{i \in \Lambda}) \le B|\Lambda|$.

(H3) $H$ is hereditary, i.e. for any $\gamma \in \mathcal{C}_f$ such that $H(\gamma) = \infty$, then, for any $x \in \mathbb{R}^d$, $H(\gamma \cup \{x\}) = \infty$.

(H4) $H$ is invariant under translation, i.e. for all $v \in \mathbb{R}^d$ and $\gamma \in \mathcal{C}_f$, $H(\gamma) = H(\gamma + v)$.
(H5.1) H has a finite range, i.e. there exists $R > 0$ such that for any $W \subset \mathbb{R}^d$ bounded, $H_W(\gamma) = H_W(\gamma_{W \oplus B(0,R)})$, where $W \oplus B(0, R) := \{x + y, x \in W, y \in B(0, R)\}$.

(H5.2) There exists $\phi : \mathbb{R}^+ \to \mathbb{R} \cup \{\infty\}$ monotone at infinity and $\epsilon > 0$ such that

$$ (2.3) \qquad H(\gamma) = \sum_{\{x,y\} \subset \gamma} \phi(|x-y|), \quad \text{and} \quad \int_{\mathbb{R}^+ \setminus [0, \epsilon]} r^{d-1} \phi(r) dr < \infty. $$

(H5.3) The local energy function $h$ is uniformly bounded.

Let us discuss these assumptions through the presentation of a few examples of standard Hamiltonians used in the literature [24, 28, 31]. The first three ones are pairwise interaction point processes, while the last one acts on geometric features of a random set. All these examples satisfy (at least) (H1)-(H4) (with some restriction on the parameter space for some of them).

1. **Strauss interaction with possible hard-core interaction.** Let $A \in \mathbb{R}$ and $0 \le r < R$. For any $\gamma \in C_f$, the Hamiltonian is given by

$$ (2.4) \qquad H(\gamma) = \sum_{\{x,y\} \in \gamma} A \mathbb{1}_{[r,R]}(|x-y|) + (\infty) \mathbb{1}_{[0,r]}(|x-y|) $$

This model satisfies (H5.1) (and thus (H5.2)). Note in particular that (H1) is satisfied if $r=0$ and $A \ge 0$ or if $r > 0$ and that (H2) is satisfied for $\epsilon < \delta - 2r$.

2. **Lennard-Jones interaction.** Let $A > 0$, $B \in \mathbb{R}$, $R > 0$ and $m > n > d$. For any $\gamma \in C_f$, the Hamiltonian is given by

$$ (2.5) \qquad H(\gamma) = \sum_{\{x,y\} \in \gamma} A \left( \frac{R}{|x-y|} \right)^m - B \left( \frac{R}{|x-y|} \right)^n. $$

This infinite range model satisfies (H5.2). Furthermore, for $\Lambda \subset \mathcal{L}$, $\epsilon < \frac{\delta}{2}$ and $|x_i| \le \epsilon$, by integral-series comparison there is $c > 0$ such that

$$ (2.6) \qquad H(\{i + x_i\}_{i \in \Lambda}) \le c|\Lambda| \int_{\delta-2\epsilon}^{\infty} r^{d-1} \left| A \left( \frac{R}{r} \right)^m - B \left( \frac{R}{r} \right)^n \right| dr, $$

and thus the interaction satisfies (H2).

3. **Quermass interaction.** Let $\theta_i \in \mathbb{R}$ for $i \in [0, d]$ and $R > 0$. For any $\gamma \in C_f$, the Hamiltonian is given by

$$ (2.7) \qquad H(\gamma) = \sum_{i=0}^{d} \theta_i M_i^d(L_R(\gamma)) $$

where $M_i^d$ is the $i$-th Minkowski functional, in particular $M_d^d$ is the Lebesgue measure, $M_{d-1}^d$ is the $d-1$-Hausdorff measure of the surface and $M_0^d$ is the Euler-Poincaré characteristic, and $L_R(\gamma)$ is the halo of a configuration, i.e.

$$ L_R(\gamma) = \bigcup_{x \in \gamma} B(x, R). $$

When $\theta_i = 0$ for $i = 0, \dots, d-1$ and $\theta_d \in \mathbb{R}$, the Quermass interaction corresponds to the one component Widom-Rowlinson model and it satisfies (H5.3). When $d=2$, the Quermass interaction verifies (H1), as proven in [24]. More precisely, they show that there is $c > 0$ such that $|M_0^2(L_R(\gamma))| \le cN(\gamma)$. However, in higher dimensions ($d \ge 3$), the
stability of this Hamiltonian remains unclear; see [24] for further discussion. Furthermore, in dimension 2, for $\Lambda \subset \mathcal{L}$, $\varepsilon > 0$ and $|x_i| \le \varepsilon$ we have

$$ (2.8) \qquad H(\{i + x_i\}_{i \in \Lambda}) \le |\Lambda| (c + |M_1^2(B(0, R))| + |M_2^2(B(0, R))|) $$

and thus the interaction verifies (H2).

**Move (or mark) distribution.** The Gibbs perturbed lattice model, precisely described in the next section, also depends on a move distribution $\mathbb{Q}$. Again to ensure probabilistic properties of the resulting Gibbs perturbed lattice, we require some of the assumptions from the following list.

(Q0) $\mathbb{Q}(B(0, \varepsilon)) > 0$ for the real number $\varepsilon$ given in (H2).

(Q1)[m]. There exists $c > 0$ such that $\mathbb{E}(e^{c|\tilde{X}|^m}) < \infty$ where $X \sim \mathbb{Q}$.

(Q2) The set $\mathbb{Q} = \text{supp}(\mathbb{Q})$ is bounded.

(Q3) $\mathbb{E}(|X|^d) < \infty$ for $X \sim \mathbb{Q}$.

Obviously, if $\mathbb{Q}$ has a continuous density on $\mathbb{Q}$, $\mathbb{Q}(B(0, \varepsilon)) > 0$ is satisfied for any $\varepsilon$ such that $B(0, \varepsilon) \cap \mathbb{Q} \ne \emptyset$. Second, (Q2) implies (Q1)[m] (for any $m > 0$) which itself implies (Q3). Three main examples are considered in this paper: (a) The uniform distribution on $\mathbb{Q}$ which satisfies (Q2); (b) the isotropic Gaussian distribution with variance $\sigma^2$, which satisfies (Q1)[m] (for $m = 1, 2$); (c) The exponential distribution on the positive orthant, i.e. the distribution with density $\lambda_1(x) = \exp(-q(x))$ with $q(x) = \theta_1 \sum_i x_i + c(\theta_1)$ for any $x \in \mathbb{Q} = \{x \in \mathbb{R}^d : x_i \ge 0, i = 1, \dots, d\}$. For this model $p_1 = 1$, $\theta_1 \in \mathbb{R}^+$ and (Q1)[1] is satisfied.

### 3. Existence and properties of Gibbs perturbed lattices.

**3.1. Gibbs measure on a finite sub-lattice.** Before we introduce and define the Gibbs perturbed lattice, we define a Gibbs measure on a finite sub-lattice with spin space $\mathbb{R}^d$ and a Hamiltonian indexed by the lattice and given by

$$ (3.1) \qquad \tilde{H}(\mathbf{x}_{\Lambda}) = H(\Gamma_{\Lambda}(\mathbf{x}_{\Lambda}, 0)) $$

for any finite $\Lambda \subset \mathcal{L}$ and any $\mathbf{x}_{\Lambda} \in (\mathbb{R}^d)^{\Lambda}$.

**DEFINITION 3.1.** For any configuration $\mathbf{x} \in (\mathbb{R}^d)^{\mathcal{L}}$ and any finite subset $\Delta \subset \mathcal{L}$, we define the local energy of $\mathbf{x}$ moved from $\Delta$ as

$$ (3.2) \qquad \tilde{H}_{\Delta}(\mathbf{x}) = \lim_{n \to \infty} \left( \tilde{H}(\mathbf{x}_{\Lambda_n}) - \tilde{H}(\mathbf{x}_{\Lambda_n \setminus \Delta}) \right) $$

where $\Lambda_n = [-n, n]^d \cap \mathcal{L}$ and with the convention $\infty - \infty = 0$.

It is worth pointing out that if $H$ has a finite range $R$, i.e. satisfies (H5.1), $\tilde{H}$ has a finite range if and only if $\mathbb{Q}$ is bounded, i.e. satisfies (Q2). In such a case, the finite range of $\tilde{H}$ is $\tilde{R} = R + \text{diam} \mathbb{Q}$ and for any finite $\Delta \subset \mathcal{L}$ we have

$$ (3.3) \qquad \tilde{H}_{\Delta}(\mathbf{x}) = \tilde{H}_{\Delta}(\mathbf{x}_{\Delta \oplus B(0, \tilde{R}) \cap \mathcal{L}}). $$

**DEFINITION 3.2.** Let $\Lambda \subset \mathcal{L}$ finite, the distribution of the Gibbs measure on the sub-lattice $\Lambda$ is given by

$$ (3.4) \qquad \mathbb{P}_{\Lambda}(d\mathbf{x}_{\Lambda}) = \frac{1}{Z_{\Lambda}} e^{-\tilde{H}(\mathbf{x}_{\Lambda})} \mathbb{Q}^{\otimes \Lambda}(d\mathbf{x}_{\Lambda}) $$
where

$$Z_{\Lambda} = \int e^{-\tilde{H}_{\Delta}(\mathbf{x}_{\Delta})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}_{\Delta})$$

is the normalizing constant, also called the partition function.

Under the assumptions (H1)-(H2) we have $(\log \mathbb{Q}(B(0, \varepsilon)) - B)|\Lambda| \le \log Z_{\Lambda} \le A|\Lambda|)$, which ensures that the Gibbs measure on the sub-lattice is well-defined. The next result states that it also satisfies the Dobrushin-Lanford-Ruelle (DLR) equations, i.e. provides the conditional distribution of the lattice for any finite subset $\Delta \subset \Lambda$ given the configuration outside of $\Delta$.

**PROPOSITION 3.3 (DLR Equations).** Let $\Delta \subset \Lambda \subset \mathcal{L}$ be finite, we have for $\mathbb{P}_{\Lambda}$-a.e $\mathbf{x}_{\Lambda \setminus \Delta}$,

$$(3.5) \qquad \mathbb{P}_{\Lambda}(d\mathbf{x}'_{\Delta}|\mathbf{x}_{\Lambda \setminus \Delta}) = \frac{1}{Z_{\Delta}(\mathbf{x}_{\Lambda \setminus \Delta})} e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Lambda \setminus \Delta})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}'_{\Delta})$$

where $Z_{\Delta}(\mathbf{x}_{\Lambda \setminus \Delta})$ is the normalizing constant given by

$$Z_{\Delta}(\mathbf{x}_{\Lambda \setminus \Delta}) = \int e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Lambda \setminus \Delta})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}'_{\Delta}).$$

**PROOF.** By definition of $H_{\Delta}$ we have for any $\mathbf{x}_{\Lambda}$, $\tilde{H}_{\Delta}(\mathbf{x}_{\Lambda}) = \tilde{H}(\mathbf{x}_{\Lambda}) - \tilde{H}(\mathbf{x}_{\Lambda \setminus \Delta})$ and thus

$$\mathbb{P}_{\Lambda}(d\mathbf{x}_{\Lambda}) = \frac{1}{Z_{\Lambda}} e^{-\tilde{H}_{\Delta}(\mathbf{x}_{\Lambda}) - \tilde{H}(\mathbf{x}_{\Lambda \setminus \Delta})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}_{\Lambda}) \mathbb{Q}^{\otimes \Lambda \setminus \Delta}(d\mathbf{x}_{\Lambda \setminus \Delta})$$

whereby we deduce that the density of the conditional probability $\mathbb{P}_{\Lambda}(d\mathbf{x}'_{\Delta}|\mathbf{x}_{\Lambda \setminus \Delta})$ with respect to $\mathbb{Q}^{\otimes \Delta}$ is: $\mathbf{x}'_{\Delta} \to e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Lambda \setminus \Delta})}$ and $Z_{\Delta}(\mathbf{x}_{\Lambda \setminus \Delta})$ is necessarily the normalizing constant. $\square$

**3.2. Infinite Gibbs perturbed lattice.** In the previous section, Gibbs measures on finite sub-lattices are shown to satisfy the DLR equations. These equations remain key ingredients in the context of an infinite lattice.

**DEFINITION 3.4.** A probability measure $\mathbb{P}$ on $(\mathbb{R}^d)^{\mathcal{L}}$ invariant under translation on $\mathcal{L}$ is a Gibbs measure for the Hamiltonian $H$ if it satisfies the DLR equations for any finite $\Delta \subset \mathcal{L}$, i.e. for $\mathbb{P}$-a.e. $\mathbf{x}_{\Delta^c}$,

$$(3.6) \qquad \mathbb{P}(d\mathbf{x}'_{\Delta}|\mathbf{x}_{\Delta^c}) = \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}'_{\Delta})$$

where

$$Z_{\Delta}(\mathbf{x}_{\Delta^c}) = \int e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} \mathbb{Q}^{\otimes \Delta}(d\mathbf{x}'_{\Delta}).$$

We denote by $\mathcal{G}_{\mathcal{L}}(H, \mathbb{Q})$ the set of all Gibbs measures.

It is worth pointing out that if $\mathcal{G}_{\mathcal{L}}(H, \mathbb{Q})$ is non-empty, then any Gibbs measure $\mathbb{P}$ on the lattice is by definition invariant under translation on $\mathcal{L}$. Therefore, if $\mathbf{X} = \{X_i, i \in \Lambda\} \sim \mathbb{P}$, then $(X_i)_{i \in \mathcal{L}}$ are identically distributed.
**DEFINITION 3.5.** Let $\mathbb{P} \in \mathcal{G}_L(H, \mathbb{Q})$, the Gibbs perturbed lattice is defined as $\mathbf{\Gamma} := \mathbf{\Gamma}(\mathbf{X}, U)$ where $\mathbf{X} \sim \mathbb{P}$, $U \sim U_D$ and $\mathbf{X}$ is independent of $U$. We denote by $P_L(H, \mathbb{Q})$ the set of stationary Gibbs perturbed lattices.

We now present our main result and provide assumptions on the point process Hamiltonian $H$ and on the mark distribution $\mathbb{Q}$ ensuring the existence of a solution to DLR equations.

**THEOREM 1.** Assume (H1)-(H4). Assume also either (H5.1)-(Q1)[1] or (H5.2)-(Q2), then $\mathcal{G}_L(H, \mathbb{Q}) \neq \emptyset$, $P_L(H, \mathbb{Q}) \neq \emptyset$ and any Gibbs perturbed lattice $\mathbf{\Gamma}$ is a stationary point process on $\mathbb{R}^d$. Moreover, for any $u \in \mathbb{R}^d$, the distribution of $\mathbf{\Gamma}$ given $U = u$, denoted by $P_u$, is invariant under any translation on $\mathcal{L}$.

The proof of Theorem 1 relies on specific-entropy techniques developed in [20] to establish the existence of a limiting process. We adapt parts of this approach to show that the limit satisfies the DLR equations. The details are provided in A–B. It is worth mentioning that, by construction, the set $P_L(H, \mathbb{Q})$ is a Choquet simplex and any Gibbs perturbed lattice is a convex combination of extremal ergodic measures. This property follows from the fact that $\mathcal{G}_L(H, \mathbb{Q})$ is a Choquet simplex; see Chapter 14 of [20] for further details. Thus, $\mathbb{P}$ and $P_u$ can be decomposed as mixtures of ergodic measures.

**3.3. Hyperuniformity of Gibbs perturbed lattices.** A point process $\mathbf{\Gamma}$ is said to be hyperuniform if, as $R \to \infty$,

$$ (3.7) \qquad \frac{\mathrm{Var}[\#\mathbf{\Gamma} \cap B(0, R)]}{|B(0, R)|} \to 0 $$

where $B(z, r)$ denotes the Euclidean ball centered at $z \in \mathbb{R}^d$ with radius $r > 0$. Hyperuniformity of perturbed lattices has been investigated in [12]. Their Theorem 1, in particular, establishes the hyperuniformity of any point process obtained from a perturbed lattice as soon as $\mathbb{E}(|X_0|^d) < \infty$, where $X_0$ is the mark from a typical location on the lattice, say 0 by stationarity. It is worth emphasizing that the global shift $U$ is necessary to establish the hyperuniformity of the perturbed lattice. The next result investigates this condition for Gibbs perturbed lattices.

**PROPOSITION 3.6.** We have the following two statements.

(i) Assume (H1)-(H4). We also assume either (H5.1)-(Q1)[d] or (H5.2)-(Q2), then there exists $\mathbb{P} \in \mathcal{G}_L(H, \mathbb{Q})$ and $\mathbf{X} \sim \mathbb{P}$ such that $\mathbb{E}(|X_0|^d) < \infty$. This implies that when $d \le 2$, there exists a hyperuniform point process $\mathbf{\Gamma} \in P_L(H, \mathbb{Q})$.

(ii) Assume (H1)-(H4) and (H5.3)-(Q3), then for any $\mathbb{P} \in \mathcal{G}_L(H, \mathbb{Q})$ and $\mathbf{X} \sim \mathbb{P}$, we have $\mathbb{E}(|X_0|^d) < \infty$. This implies that when $d \le 2$, any $\mathbf{\Gamma} \in P_L(H, \mathbb{Q})$ is hyperuniform.

The key ingredient for the proof of Proposition 3.6, proposed in C, relies on the variational principle that requires the specific entropy to be finite; see (A.3). What distinguishes (i) from (ii) is precisely the fact that specific entropies are not necessarily finite for every Gibbs measure.

**3.4. DLR type equations.** The Gibbs perturbed lattice satisfies several equations that are derived from the DLR equations (3.6). These equations can be derived not only for the expectation with respect to the Gibbs measure but also for the conditional expectation given $U = u \in D$, denoted by $\mathbb{E}_u[\cdot]$ in the rest of the paper.
We start with a presentation of the DLR equation at a single site $i \in \mathcal{L}_U$ by applying equation (3.6) with $\Delta = \{i\}$. Let $f: C \to \mathbb{R}$ be a bounded measurable function. Then, we have

$$ (3.8) \qquad \mathbb{E}_u[f(\mathbf{\Gamma})] = \mathbb{E}_u \left[ \frac{1}{Z_i(\mathbf{\Gamma}_{i^c})} \int f(\mathbf{\Gamma}_{i^c} \cup \{i+x+U\})e^{-h(i+x+U, \mathbf{\Gamma}_{i^c})} \mathcal{Q}(dx) \right] $$

where for $i \in \mathcal{L}_U = \{i+U, i \in \mathcal{L}\}$ and $\mathcal{L}_U$ is the shifted lattice,

$$ (3.9) \qquad \mathbf{\Gamma}_{i^c} = \{j + X_j, j \in \mathcal{L}_U, j \neq i\}. $$

In the following, we assume that $\mathcal{Q}$ is absolutely continuous with respect to the Lebesgue measure and we denote by $\lambda_1(x) = e^{-q(x)}$ its density with support $\mathcal{Q}$. We also denote, for any $i \in \mathcal{L}_U$, $x \in \mathcal{Q}$ and $\gamma \in \mathcal{C}$, $\lambda_2(i+x, \gamma) = e^{-h(i+x, \gamma)}$, the two measurable functions

$$ (3.10) \qquad \lambda(i, x, \gamma) := \lambda_1(x)\lambda_2(i+x, \gamma) \quad \text{and} \quad \Lambda(i, x, \gamma) = \frac{\lambda(i, x, \gamma)}{Z_i(\gamma)} = \frac{\lambda(i, x, \gamma)}{\int_{\mathcal{Q}} \lambda(i, x, \gamma)dx}. $$

These functions respectively act as unnormalised and normalised Papangelou conditional intensities. We are now in a position to state the main equations used in Section 4 to derive parametric inference methods. The next three results are proved in D.

**PROPOSITION 3.7.** Assume the conditions of Theorem 1 are satisfied and thus let $\mathbf{\Gamma} \sim \mathbb{P} \in P_{\mathcal{L}}(H, \mathcal{Q}) \neq \emptyset$. For $u \in \mathcal{D}$, let $f: \mathcal{L}_u \times \mathbb{R}^d \times C \to \mathbb{R}$ be any measurable function.

(i) Assume that

$$ \mathbb{E}_u \sum_{i \in \mathcal{L}_U} |f(i, X_i, \mathbf{\Gamma}_{i^c})| < \infty \quad \text{and} \quad \mathbb{E}_u \int_{\mathcal{Q}} |f(i, x, \mathbf{\Gamma}_{i^c})| \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx < \infty, $$

then we have

$$ (3.11) \qquad \mathbb{E}_u \sum_{i \in \mathcal{L}_U} f(i, X_i, \mathbf{\Gamma}_{i^c}) = \mathbb{E}_u \sum_{i \in \mathcal{L}_U} \int_{\mathcal{Q}} f(i, x, \mathbf{\Gamma}_{i^c}) \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx. $$

(ii) Assume that $f$ is such that $f(i, x, \gamma) = f(0, x, \gamma - i)$ for any $i \in \mathcal{L}_u$, $x \in \mathcal{Q}$, $\gamma \in \mathcal{C}$ and such that

$$ \mathbb{E}_u \int_{\mathcal{Q}} |f(0, x, \mathbf{\Gamma}_{0^c})| \Lambda(0, x, \mathbf{\Gamma}_{0^c}) dx < \infty, $$

then, for any bounded domain $W \subset \mathbb{R}^d$, we have

$$ (3.12) \qquad \mathbb{E}_u \sum_{i \in \mathcal{L}_U \cap W} f(i, X_i, \mathbf{\Gamma}_{i^c}) = \mathbb{E}_u \sum_{i \in \mathcal{L}_U \cap W} \int_{\mathcal{Q}} f(i, x, \mathbf{\Gamma}_{i^c}) \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx. $$

Starting from the DLR equations on two sites, we can obtain a similar second-order GNZ-type equation [19, 20, 29].

**PROPOSITION 3.8.** Assume the conditions of Theorem 1 are satisfied and thus let $\mathbf{\Gamma} \sim \mathbb{P} \in P_{\mathcal{L}}(H, \mathcal{Q}) \neq \emptyset$. For any $u \in \mathcal{D}$ and any measurable function $f: \mathcal{L}_u^2 \times (\mathbb{R}^d)^2 \times C \to \mathbb{R}$ such that

$$ \mathbb{E}_u \sum_{\{i,j\} \in \mathcal{L}_U} |f(i, j, X_i, X_j, \mathbf{\Gamma}_{\{i,j\}^c})| < \infty $$

$$ \mathbb{E}_u \int_{\mathcal{Q}} \int_{\mathcal{Q}} |f(i, j, x, y, \mathbf{\Gamma}_{\{i,j\}^c})| \Lambda(i, x, \mathbf{\Gamma}_{\{i,j\}^c}) \Lambda(j, y, \mathbf{\Gamma}_{\{i,j\}^c} \cup \{i+x\}) dx dy < \infty. $$
Then, we have

$$
\begin{align*}
\mathbb{E}_u \sum_{\{i,j\} \in \mathcal{L}_U} f(i, j, X_i, X_j, \mathbf{\Gamma}_{\{i,j\}^c}) &= \\
&\mathbb{E}_u \sum_{\{i,j\} \in \mathcal{L}_U} \iint_{\mathbb{Q}^2} f(i, j, x, y, \mathbf{\Gamma}_{\{i,j\}^c}) \Lambda(i, x, \mathbf{\Gamma}_{\{i,j\}^c}) \Lambda(j, y, \mathbf{\Gamma}_{\{i,j\}^c} \cup \{i+x\}) dy dx.
\end{align*}
$$

To complete these DLR-type equations, we now propose a variational equation, similarly to the one obtained by [2]. This requires a few additional definitions and notation. Let $\tilde{f}: \mathbb{R}^d \times \mathcal{C} \to \mathbb{R}$ be a measurable function and $x \in \mathbb{R}^d$. We say that $\tilde{f}$ is differentiable at $x$, if for every $\gamma \in \mathcal{C}$ the function: $x \mapsto \tilde{f}(x, \gamma)$ is differentiable at $x$. Furthermore, $\tilde{f}$ is $\lambda^d$-a.e. differentiable, if $\tilde{f}$ is differentiable for $\lambda^d$-a.e. $x \in \mathbb{R}^d$. We denote by $\nabla \tilde{f}(x, \gamma)$ the corresponding gradient. By convention, if $\tilde{f}(y) = \infty$ for $y \in \mathcal{V}$, where $\mathcal{V}$ is an open neighborhood of $x$, we have $\nabla \tilde{f}(x, \gamma) = 0$. Let $\tilde{f}, h: \mathbb{R}^d \times \mathcal{C} \to \mathbb{R}$ be measurable functions. We say that $\tilde{f}$ is regularizing with respect to $h$ if for every $\gamma \in \mathcal{C}$ the function

$x \mapsto \tilde{f}(x, \gamma) e^{-h(x, \gamma)}$

is absolutely continuous. As a consequence, if $\tilde{f}$ and $h$ are $\lambda^d$-a.e. differentiable and $\tilde{f}$ is regularizing with respect to $h$, then we have

$$
(3.13) \quad \nabla \left( \tilde{f} e^{-h} \right) (x, \gamma) = \left( \nabla \tilde{f}(x, \gamma) - \tilde{f}(x, \gamma) \nabla h(x, \gamma) \right) e^{-h(x, \gamma)}.
$$

Finally, for $u \in \mathcal{D}$, we denote by $\mathcal{F}_R(H, \mathbb{Q}, \mathcal{L}_u)$ the class of functions that are $\lambda^d$-a.e. differentiable, regularizing with respect to $h+q$, and such that for any $x \in \partial \mathbb{Q}$, $i \in \mathcal{L}_u$ and $\gamma \in \mathcal{C}$, we have $\tilde{f}(x+i, \gamma) = 0$ and $\tilde{f}$ vanishes at infinity.

**PROPOSITION 3.9.** Assume the conditions of Theorem **1** are satisfied and thus let $\Gamma \sim \mathbb{P} \in P_{\mathcal{L}}(H, \mathbb{Q}) \neq \emptyset$. For $u \in \mathcal{D}$ and $\tilde{f} \in \mathcal{F}_R(H, \mathbb{Q}, \mathcal{L}_u)$ such that

$$
(3.14) \quad \mathbb{E}_u \sum_{i \in \mathcal{L}_U} \left\{ \left| \nabla \tilde{f}(i+X_i, \mathbf{\Gamma}) \right| + \left| \tilde{f}(i+X_i, \mathbf{\Gamma}) \right| \left| \nabla h(i+X_i, \mathbf{\Gamma}_{i^c}) + \nabla q(X_i) \right| \right\} < \infty,
$$

then we have the variational equation

$$
(3.15) \quad \mathbb{E}_u \sum_{i \in \mathcal{L}_U} \left\{ \nabla \tilde{f}(i+X_i, \mathbf{\Gamma}) - \tilde{f}(i+X_i, \mathbf{\Gamma}) (\nabla h(i+X_i, \mathbf{\Gamma}_{i^c}) + \nabla q(X_i)) \right\} = 0.
$$

**REMARK 3.10.** It is worth mentioning that Proposition 3.9, and in particular the variational equation (3.15), can be seen as a specific case of Proposition 3.7 by taking $f(i, x, \gamma) = \nabla \tilde{f}(i+x, \gamma) - \tilde{f}(i+x, \gamma)(\nabla h(i+x, \gamma) + \nabla q(x_i))$.

4. Parametric inference of Gibbs perturbed lattices.

4.1. Notation and parametric model. In this section, we assume that assumptions (**H1**)-(**H4**) and (**H5.1**)-(**Q1**)[1] are satisfied and thus, by Theorem **1**, there exists $\mathbf{\Gamma} \sim \mathbb{P} \in P_{\mathcal{L}}(H, \mathbb{Q})$. Assuming the finite range property (**H5.1**) eases in many ways the statistical methodology and its asymptotic properties. We leave for future research the situation where the Hamiltonian satisfies (**H5.2**), as is the case, for example, for the Lennard-Jones potential, as considered by [7, 8] for continuous Gibbs point processes.

We recall the notation $\mathbf{\Gamma} = \{i + X_i, i \in \mathcal{L}_U\}$, where $X_i$ are realizations of the mark distribution $\mathbb{Q}$, assumed to have a density $\mathbb{Q}$ with known support $\mathbb{Q}$ (possibly $\mathbb{R}^d$). Finally, we also remind that the notation $\mathbf{\Gamma}_{i^c}$, given by (3.9), stands for the set of perturbed points except the one coming from $i$. For any $u \in \mathcal{D}$ and any function $f: \mathcal{L}_u \times \mathbb{Q} \times \mathcal{C} \to \mathbb{R}$, we consider the following properties and definitions:
(TI) $f$ is said to be invariant under translation if for any $i \in \mathcal{L}_u$, $x \in \mathcal{Q}$ and $\gamma \in \mathcal{C}$, $f(i, x, \gamma) = f(0, x, \gamma - i)$.

(FR) $f$ is said to have a finite range $R > 0$ if for any $i \in \mathcal{L}_u$, $x \in \mathcal{Q}$ and $\gamma \in \mathcal{C}$, $f(i, x, \gamma) = f(i, x, \gamma_{B(i+x,R)})$.

In particular, the local energy function $\tilde{h}(i, x, \gamma) = h(i + x, \gamma)$ as well as $\lambda_2(i + x, \gamma)$ satisfy (TI)-(FR).

This section focuses on parametric modelling and estimation: for any $u \in \mathcal{D}$, $i \in \mathcal{L}_u$, $x \in \mathcal{Q}$, $\gamma \in \mathcal{C}$, we let $\lambda(i, x, \gamma; \theta)$ (resp. $\Lambda(i, x, \gamma; \theta)$) denote a parametric version of (3.10), where $\theta \in \Theta \subseteq \mathbb{R}^p$ for some $p \ge 1$. Hence, we consider parametric forms of $q$ and the local energy function $h$. To ease statistical methodologies and asymptotic results, it is standard to focus on exponential family models [23]. So, on the one hand, we assume that $-\log \lambda_1(x) = q(x)$ is of the form $\theta_1^\top S_1(x) + c(\theta_1)$ where $\theta_1 \in \mathbb{R}^{p_1}$ for some $p_1 \ge 1$ and $c(\theta_1)$ is the log-partition constant. Here are a few examples:

* Uniform distribution on $\mathcal{Q}$ (with $|\mathcal{Q}| < \infty$). In this case, there is no parameter to estimate and we take the convention $p_1 = 0$.

* Centered Gaussian density on $\mathcal{Q} = \mathbb{R}^d$: $p_1 = 1$, $\theta_1 \in \mathbb{R}^+$ and $S_1(x) = \|x\|^2$ and $c(\theta_1) = -d/2 \log(\pi\theta_1)$.

* Exponential distribution on $\mathcal{Q} = \{x \in \mathbb{R}^d : x_i \ge 0, i = 1, \dots, d\}$: $p_1 = 1$, $\theta_1 \in \mathbb{R}^+$, $S_1(x) = \sum_{i=1}^d x_i$ and $c(\theta_1) = d \log(\theta_1)$.

On the other hand, we assume that the local energy function $h = -\log \lambda_2$ can be written as $-\log \lambda_2(i + x, \gamma) = h(i + x, \gamma) = \theta_2^\top S_2(i + x, \gamma)$ where $\theta_2 \in \mathbb{R}^{p_2}$, $p_2 \ge 1$. By assumption, $\lambda_2$, $h$ and $S_2$ satisfy (TI)-(FR).

As a first example, let us cite the well-known Strauss hard-core model for which $p_2 = 1$ and

$$S_2(y, \gamma) = \sum_{y' \in \gamma} \{1_{[r,R]}(|y' - y|) + (\infty) \times 1_{[0,r]}(|y' - y|)\}$$

for any $y \in \mathbb{R}^d$ and known real parameters $r, R$ such that $0 \le r < R$. The exponential family models also include for instance piecewise versions of the Strauss hard-core model, the Geyer saturation model and the Widom-Rowlinson area-interaction [14, 22, 28]. Note that all these models have finite range. Moreover, they also all satisfy some local stability property which means that for any $y \in \mathbb{R}^d$ and $\gamma \in \mathcal{C}$, $e^{-h(y, \gamma)} \le \kappa$ for some $\kappa > 0$ (independent of $y, \gamma$).

As a summary let $\theta = (\theta_1^\top, \theta_2^\top)^\top \in \mathbb{R}^p$ and for any $x \in \mathcal{Q}$, $S(i, x, \gamma) := (S_1(x)^\top, S_2(i+x, \gamma)^\top)^\top$ with $p = p_1 + p_2$. All in all, we consider the following exponential parametric model

$$\begin{align}
-\log \lambda(i, x, \gamma; \theta) &= -\log \lambda_1(x; \theta) - \log \lambda_2(i + x, \gamma; \theta) = h(i + x, \gamma) + q(x) \tag{4.1} \\
&= \theta_1^\top S_1(x) + \theta_2^\top S_2(i + x, \gamma) \nonumber
\end{align}$$

Finally, we denote by $\theta^*$ the true parameter vector to be estimated, i.e. we assume that $\Gamma \sim \mathbb{P} = \mathbb{P}_{\theta^*}$.

4.2. Statistical context and methodologies. In this paper, we consider two different statistical settings, illustrated by Figure 1. Let $W_n$ be a bounded domain of $\mathbb{R}^d$ which plays the role of the observation domain.

* **Framework 1.** We observe the shifted lattice $\mathcal{L}_U$ and all perturbed points falling in $W_n$, i.e. we observe a realization of

$$\mathcal{D}_{1,n} = \{(i, X_i), \text{ with } i \in \mathcal{L}_U \cap W_n \text{ and } i + X_i \in W_n\}.$$
* **Framework 2.** We observe the shifted lattice $\mathcal{L}_U$, and only the perturbed points falling in $W_n$, i.e. we observe a realization of

$$ (4.3) \qquad \mathcal{D}_{2,n} = \{i + X_i \in W_n \text{ for } i \in \mathcal{L}_U\}. $$

<figure><img src="image_2.png" /><figcaption>Fig 1: Illustration of the two statistical frameworks considered in this paper. In (a)-(b) the inside grey rectangle corresponds to $W_n$, crosses to points $i \in \mathcal{L}_U$ and dots to $i + x_i$, for $i \in \mathcal{L}_U$. For (a) Framework 1, we observe crosses (i) and dots (i + x_i) that fall in $W_n$. For (b) Framework 2, we observe only dots falling in $W_n$.</figcaption></figure>

The rest of this section discusses strategies to estimate $\theta^*$ based on a single observation of either $(i, x_i)_i$ or $(i + x_i)_i$ on a bounded domain $W_n \subset \mathbb{R}^d$. The domain $W_n$ plays the role of the observation domain and grows to $\mathbb{R}^d$ as $n \to \infty$. We also let $W_n \ominus r$ for some $r \ge 0$ be the domain $W_n$ eroded by $r$. The standard approach would consist in using the likelihood method. However, due to the unknown normalizing constant depending on the parameter, it is well-known that the likelihood method is already computationally expensive even for estimating an unmarked Gibbs model on a lattice [18]. Many alternative methods have been proposed in the literature for Gibbs models on the lattice or in the continuous case, which are in particular appropriate for Framework 1 [18]. Unlike standard Gibbs lattice models, it is worth noting that this model has a continuous mark distribution and has the particularity that the Hamiltonian depends only on $i + X_i$, $i \in \mathcal{L}_U$ and $X_i \in \mathcal{Q}$. The most popular alternative method is the (log-)pseudolikelihood. To account for edge-effects problems, we have to introduce a specific border correction factor, $b_n$, described below. The border-corrected (log)-pseudolikelihood we consider is given by

$$ (4.4) \qquad \text{LPL}_n(\Gamma; \theta) = \sum_{i \in \mathcal{L}_U} \log \Lambda_n(i, X_i, \Gamma_{i^c}; \theta), $$

where for any $i \in \mathcal{L}_U$, $x \in \mathcal{Q}$ and $\gamma \in \mathcal{C}$

$$ (4.5) \qquad \Lambda_n(i, x, \gamma; \theta) = \frac{b_n(i, x) e^{-\theta^\top S(i, x, \gamma)}}{\int_{\mathcal{Q}} b_n(i, y) e^{-\theta^\top S(i, y, \gamma)} dy} $$

where the term $b_n$ is called a border correction factor. For technical reasons, mainly explained in Theorem 2(ii), we define this border correction factor as

$$ (4.6) \qquad b_n(i, x) = \mathbb{1}(i \in W_n \ominus (m_n + R), i + x \in W_n \ominus R) $$
where $(m_n)_n$ is a sequence of real numbers, specified later in assumption (M1). It is aimed to tend to $\infty$ such that $|W_n \ominus (m_n+R)| \sim |W_n|$ as $n \to \infty$. The interest of the eroded domain $W_n \ominus R$ is also revealed. By the finite range assumption, for any $i \in \mathcal{L}_u$ such that $i+z \in W_n \ominus R$ for $z = X_i, x \in \mathcal{Q}$, $\Lambda(i,z, \Gamma_{ic}; \theta) = \Lambda(i,z, \Gamma_{ic} \cap W_n; \theta)$ can always be evaluated with data from Framework 1. The gradient (with respect to $\theta$), denoted by $\text{LPL}_n^{(1)}(\Gamma; u, \theta)$ of $\text{LPL}_n(\Gamma; u, \theta)$ for exponential family models is given by

$$ (4.7) \qquad \text{LPL}_n^{(1)}(\Gamma; \theta) = \sum_{i \in \mathcal{L}_U} \left\{ -b_n(i, X_i)S(i, X_i, \Gamma_{ic}) + \int_{\mathcal{Q}} b_n(i, x)S(i, x, \Gamma_{ic})\Lambda_n(i, x, \Gamma_{ic}; \theta)dx \right\}. $$

At that point, it is worth observing that, under appropriate assumptions, not detailed for now, and replacing $\Lambda_n$ by $\Lambda$ in the previous equation, Proposition 3.7 can be applied and would show that (4.7) is centered when $\theta = \theta^*$.

Let us now consider the situation of Framework 2. Since we lose the links between $i$ and $i+X_i$, it is straightforwardly seen that some elements of (4.7) cannot be evaluated: (a) for $i+X_i \in \Gamma$, $q(X_i)$ and consequently $S(i, X_i, \Gamma_{ic})$ are not observed; (b) for any $x \in \mathcal{Q}$, $S(i, x, \Gamma_{ic})$ is also incomputable since $\Gamma_{ic}$ is not observed. In this situation we did not come up with a general estimating equation similar to (4.7) that would overcome problems mentioned in (a)-(b). However, for some particular cases, the variational equation (3.15), allows us to partially solve the problem. Let $\tilde{f} = (\tilde{f}_1, \dots, \tilde{f}_L)$, for some $L \ge p$, be a set of functions satisfying the assumptions of Proposition 3.9, and consider the empirical version of (3.15) given, for the whole set of test functions $\tilde{f}$ and for exponential family models, by

$$ (4.8) \qquad \text{VARE}_n(\tilde{f}; \theta) = \sum_{\substack{i \in \mathcal{L}_U \\ i+X_i \in W_n \ominus R}} \left\{ \nabla \tilde{f}(i+X_i; \Gamma_{ic}) - \tilde{f}(i+X_i; \Gamma_{ic}) (\nabla S_1(X_i)^\top \theta_1 + \nabla S_2(i, X_i, \Gamma_{ic})^\top \theta_2) \right\}. $$

If one considers Framework 2, the main problem of (3.15) lies in $\nabla S_1$, which is not observable. However, we can exhibit two particular examples, for which this term either cancels out or is equal to one:

(i) The distribution $\mathbb{Q}$ is the uniform distribution on $\mathcal{Q}$, then in this case $p_1 = \nabla S_1 = 0$ and we propose to take $L = p$ and consider $\tilde{f}(i+x, \gamma) = \nabla S_2(i+x, \gamma)\psi(i+x, \gamma)$ with $\psi: \mathbb{R}^d \to \mathbb{R}$.

(ii) The density $\lambda_1(x) = e^{-q(x)}$ corresponds to the exponential distribution on the positive orthant described in Section 2.2. For this specific distribution, $\nabla S_1 = 1$ and we propose to take $L = p$ and consider $\tilde{f}(i+x, \gamma) = (1 + \nabla S_2(i+x, \gamma))\psi(i+x, \gamma)$ with $\psi: \mathbb{R}^d \to \mathbb{R}$.

The variational equation has also the advantage for exponential models to be linear in $\theta$ which yields an explicit and simple estimator. Remark that for both cases (i)-(ii), one can always find a smooth function $\psi$ such that given $U = u$, $\tilde{f} \in \mathcal{F}_R(H, \mathbb{Q}, \mathcal{L}_u)$. Hence, Proposition 3.9 can be applied which also means that given $U = u$ (3.15) is an unbiased estimating equation, i.e. $E_u[\text{VARE}_n(\tilde{f}; \theta)] = 0$ when $\theta = \theta^*$.

From a practical point of view, methodologies summarized by (4.7) and (4.8) appear to be different since the first one (resp. second one) is more appropriate for Framework 1 (resp. Framework 2) and that edge effects are treated differently. Nevertheless, from a theoretical point of view, by Remark 3.10, estimates obtained from these equations follow from the same DLR-type equation (3.6). Therefore to encompass both methods from a theoretical point of
view, we propose to consider the Takacs-Fiksel [6, 15, 33] methodology for which (4.7) and (4.8) are particular cases.

Let $L \ge p$ and for any $u \in \mathcal{D}$, let $\mathbf{f}: \mathcal{L}_u \times \mathcal{Q} \times \mathcal{C} \times \mathbb{R}^p \to \mathbb{R}^L$. Then, we construct the least squares criterion, called Takacs-Fiksel criterion

$$ (4.9) \qquad \text{TF}_n(\mathbf{f}, \mathbf{\Gamma}; \theta) = \sum_{l=1}^{L} \text{DLR}_n(f_l, \mathbf{\Gamma}; \theta)^2 \quad \text{with} \\ (4.10) \qquad \qquad \text{DLR}_n(f_l, \mathbf{\Gamma}; \theta) = \sum_{i \in \mathcal{L}_u} \left\{ b_n(i, X_i) f_l(i, X_i, \mathbf{\Gamma}_{i^c}; \theta) \right. \\ \qquad \qquad \qquad \left. - \int_{\mathcal{Q}} b_n(i, x) f_l(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda_n(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx \right\} $$

The Takacs-Fiksel estimator is defined as $\hat{\theta} = \argmin_{\theta} \text{TF}_n(\mathbf{f}, \mathbf{\Gamma}; \theta)$. Again, $\hat{\theta}$ corresponds to the pseudo-likelihood estimator (4.4) when $f(i, x, \gamma; \theta) = S(i, x, \gamma)$ and to the variational estimator (4.8) for the two particular cases (i)-(ii) presented above when $f(i, x, \gamma) = \nabla \hat{f}(i + x, \gamma) - \hat{f}(i + x, \gamma)(\nabla h(i + x, \gamma) + \nabla q(x_i))$ and $\hat{f}(i + x, \gamma) = (z + \nabla S_2(i, x, \gamma))\psi(i + x, \gamma)$ (with $z$ the $d$-dimensional 0 or 1 vector). We intrinsically assume in the following that in the latter case, $\psi$ is chosen such that given $U = u, \hat{f} \in \mathcal{F}_R(H, \mathcal{Q}, \mathcal{L}_u)$.

For any function or stochastic function of $\theta$, $\ell(\theta)$, we denote by $\ell^{(1)}(\theta)$ and $\ell^{(2)}(\theta)$ its gradient and Hessian respectively. Minimizing (4.9) is equivalent to finding the zeros of its gradient

$$ \text{TF}_n^{(1)}(\mathbf{f}, \mathbf{\Gamma}; \theta) = 2 \sum_{l=1}^{p} \text{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta) \times \text{DLR}_n(f_l, \mathbf{\Gamma}; \theta) $$

### 4.3. Asymptotic results.
In this section, we investigate the properties of the Takacs-Fiksel estimator minimizing (4.9). We let (M1)-(M5) denote the set of assumptions on the model and test functions $f_j$:

(M1) We assume that the perturbed Gibbs lattice model satisfies (H1)-(H4), (H5.1)-(Q1)[1], is parameterized by (4.1) where $\theta \in \Theta$, an open bounded convex set of $\mathbb{R}^p$, with interior containing $\theta^*$ the true parameter vector to be estimated, and is observed in the sequence of increasing cubes $(W_n)_n$ with volume $|W_n| \to \infty$ as $n \to \infty$. We assume that the sequence of real numbers $(m_n)$ defined in (4.6) is such that $m_n = \beta \log |W_n|$ for $\beta > 1/2$ such that $\sum_{n \ge 1} |W_n|^{-\beta} < \infty$.

(M2) Let $L \ge p$. For any $l = 1, \dots, L$, $u \in \mathcal{D}$, $i \in \mathcal{L}_u$, $x \in \mathbb{R}^d$ and $\gamma \in \mathcal{C}$, the functions $f_l(i, x, \gamma; \cdot)$ and $\Lambda(i, x, \gamma; \cdot)$ are twice continuously differentiable on $\Theta$. For $k = 0, 1, 2$, we assume that each component $f_l^{(k)}(i, x, \gamma; \theta)$ can be decomposed as $\phi_1(x; \theta)\phi_2(i+x, \gamma; \theta)$, $\phi_1, \phi_2$ being different for different components of $f_l^{(k)}$, such that $\phi_2$ satisfies (TI)-(FR).

(M3) For any $j, k = 1, \dots, p$, $l = 1, \dots, L$ and $\theta \in \Theta$, (A)[f] is satisfied for $f = f_l$, $(f_l^{(1)})_j$ and $(f_l^{(2)})_{jk}$, and that (B)[f] is satisfied for $f = f_l$, $(f_l^{(1)})_j\Lambda$, $f_l\Lambda_j^{(1)}/\Lambda \times$, $(f_l^{(2)})_{jk}$, $(f_l^{(1)})_j(f_l^{(1)})_k$, $(f_l^{(1)})_j(\Lambda^{(1)})_j/\Lambda$ and $f_l(\Lambda^{(2)})_{jk}/\Lambda$. The assumptions (A)[f]-(B)[f] correspond to:
(A)[f] For any $u \in \mathbb{R}^d$ and $\theta \in \Theta$

$$ \sup_{x \in \mathcal{Q}} \mathbb{E}_u \left( |f(0, x, \mathbf{\Gamma}_{0^c}; \theta)| \lambda_2(0, x, \mathbf{\Gamma}_{0^c}; \theta) \right) < \infty. $$
(B)[f] For any $u \in \mathbb{R}^d$ and $\theta \in \Theta$

$$ \sup_{x \in \mathcal{Q}} \mathbb{E}_u \left( |f(0, x, \mathbf{\Gamma}_{0^c}; \theta)| \frac{\Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta)}{\lambda_1(x; \theta^*)} \right) < \infty $$

$$ \sup_{n \ge 1} \sup_{x,y \in \mathcal{Q}} \mathbb{E}_u \left( |f(0, x, \mathbf{\Gamma}_{0^c}; \theta)| \frac{\Lambda_n(0, x, \mathbf{\Gamma}_{0^c}; \theta) \Lambda(0, y, \mathbf{\Gamma}_{0^c}; \theta)}{\lambda_1(x; \theta^*) \lambda_1(y; \theta^*)} \right) < \infty. $$

(M4) We assume that for any $u \in \mathcal{D}$, the quantity $\mathcal{E}_u$ defined by (F.5) is finite.

(M5) $\lambda_{\min}(\mathcal{I}(\theta^*)) > 0$ where $\mathcal{I}(\theta^*) = \sum_{l=1}^L \mathcal{I}_l(\theta^*) \mathcal{I}_l(\theta^*)^\top$ with

$$ (4.11) \qquad \mathcal{I}_l(\theta^*) = \mathbb{E}_u \left[ \int_{\mathcal{Q}} f_l(0, x, \mathbf{\Gamma}_{0^c}; \theta^*) \Lambda^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta^*) dx \right] $$

and where $\lambda_{\min}(M)$, for a symmetric squared matrix $M$, stands for the smallest eigenvalue of $M$.

Assumption (M1) ensures in particular, from Theorem 1, the existence of at least one perturbed lattice stationary Gibbs measure. The assumption on $W_n$ is for instance satisfied for the simple sequence of cubes $W_n = [-n^{1/d}/2, n^{1/d}/2]^d$ with volume $|W_n| = n$ and $\beta > 1$. Domains other than cubes could be considered. Finally, a look at the proof of Theorem 2 in particular shows that the sequence $m_n$ is actually related to the tail of the distribution $\mathcal{Q}$. More precisely it corresponds to an upper-bound of $-\log \mathbb{P}(|X| > m_n)$ for $X \sim \mathcal{Q}$. Note that if $\mathcal{Q}$ is compactly supported, such that $\mathbb{P}(|X| \ge K) = 0$, one may replace $m_n$ by $K$ in the methodology. Assumption (M2) requires specific form of the test functions $f_j$ and regularity conditions. For exponential models, such an assumption is satisfied for test functions leading to (4.4) and (4.8). Assumption (M3) is essentially used in Theorem 2 to prove ergodic results for empirical versions of DLR-type equations. Assumption (M4) is used to control the conditional variance of $\overline{\text{DLR}}_n(f_l; \theta^*)$. Finally, Assumption (M5) corresponds to an identifiability condition.

We start with an ergodic theorem for DLR-type equations. Before that, for $A_n$ a random variable, $a_n$ a sequence of positive real numbers and $u \in \mathcal{D}$, we write $A_n = O_{\mathbb{P}_u}(a_n)$ if the following property is satisfied: for any $\varepsilon = \varepsilon(u)$, there exists some finite $M = M(u) > 0$ such that for $n$ large enough

$$ (4.12) \qquad \mathbb{P}(|A_n| \ge M a_n | U = u) \le \varepsilon. $$

It is worth mentioning that if $\mathbb{E}_u |A_n|^p < \infty$ for some $p > 0$, by conditional Markov's inequality, we have that $A_n = O_{\mathbb{P}_u}((\mathbb{E}_u |A_n|^p)^{1/p})$.

**THEOREM 2.** Assume (M1) and let $f$ be a function satisfying (M2) and (M3). We define the random variables

$$ A_{W_n}(f, \mathbf{\Gamma}; \theta) = \sum_{i \in \mathcal{L}_U} b_n(i, X_i) f(i, X_i, \mathbf{\Gamma}_{i^c}; \theta) $$

$$ B_{W_n}(f, \mathbf{\Gamma}; \theta) = \sum_{i \in \mathcal{L}_U} \int_{\mathcal{Q}} b_n(i, x) f(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda_n(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx $$

$$ \tilde{B}_{W_n}(f, \mathbf{\Gamma}; \theta) = \sum_{i \in \mathcal{L}_U} \int_{\mathcal{Q}} b_n(i, x) f(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx $$

$$ \text{DLR}_n(f, \mathbf{\Gamma}; \theta) = A_{W_n}(f, \mathbf{\Gamma}; \theta) - B_{W_n}(f, \mathbf{\Gamma}; \theta) $$

$$ \overline{\text{DLR}}_n(f, \mathbf{\Gamma}; \theta) = A_{W_n}(f, \mathbf{\Gamma}; \theta) - \tilde{B}_{W_n}(f, \mathbf{\Gamma}; \theta). $$
Then, as $n \to \infty$, given $U = u$, we have the following two statements.

(i)

$$ (4.13) \qquad |W_n|^{-1} A_{W_n}(f, \Gamma; \theta) \xrightarrow{a.s.} \mathbb{E}_u \int_Q f(0, x, \Gamma_{0c}; \theta) \Lambda(0, x, \Gamma_{0c}; \theta^*) dx $$

$$ (4.14) \qquad |W_n|^{-1} B_{W_n}(f, \Gamma; \theta) \xrightarrow{a.s.} \mathbb{E}_u \int_Q f(0, x, \Gamma_{0c}; \theta) \Lambda(0, x, \Gamma_{0c}; \theta) dx $$

$$ (4.15) \qquad |W_n|^{-1} \mathrm{DLR}_n(f, \Gamma; \theta) \xrightarrow{a.s.} \mathbb{E}_u \int_Q f(0, x, \Gamma_{0c}; \theta) \{ \Lambda(0, x, \Gamma_{0c}; \theta^*) - \Lambda(0, x, \Gamma_{0c}; \theta) \} dx $$

$$ (4.16) \qquad |W_n|^{-1} (\mathrm{DLR}_n(f, \Gamma; \theta) - \widetilde{\mathrm{DLR}}_n(f, \Gamma; \theta)) \xrightarrow{a.s.} 0. $$

(ii)

$$ (4.17) \qquad \mathrm{DLR}_n(f, \Gamma; \theta) - \widetilde{\mathrm{DLR}}_n(f, \Gamma; \theta) = O_{\mathbb{P}_u} \left( |W_n|^{1-\beta} \right) $$

where $\beta$ is defined in (M1).

The proof of this result is postponed to E. The next result establishes the consistency of the Takacs-Fiksel estimator.

**THEOREM 3.** Under the assumptions (M1)-(M5), we have the following two statements obtained as $n \to \infty$.

(i) Let $v_n = \sup_l \mathrm{Var}_u[\widetilde{\mathrm{DLR}}_n(f_l, \Gamma; \theta^*)]$, where $\mathrm{Var}_u[\widetilde{\mathrm{DLR}}_n(f_l, \Gamma; \theta^*)]$ stands for the conditional variance of $\widetilde{\mathrm{DLR}}_n(f_l, \Gamma; \theta^*)$ given $U = u$. Then, $v_n = O(|W_n|)$.

(ii) The sequence of local maximizers $\hat{\theta}$ satisfies

$$ \hat{\theta} - \theta^* = O_{\mathbb{P}_u} \left( v_n^{1/2} |W_n|^{-1} \right) = O_{\mathbb{P}_u} \left( |W_n|^{-1/2} \right). $$

Let us comment on the rate of convergence obtained in (ii). As mentioned in the introduction, standard continuous Gibbs point processes are non hyperuniform and results regarding the Takacs-Fiksel estimator for these point processes, see [6], show that the rate of convergence is $|W_n|^{-1/2}$. As a comparison, we obtain an upper-bound for the rate of convergence with similar rate. However, we also show that this rate can be reduced if one proves that $v_n = o(|W_n|)$. If $\Gamma$ is hyperuniform, it can be deduced (immediate if $|W_n|$ is an increasing ball) that $\mathrm{Var}[N(W_n)] = o(|W_n|)$. However, $N(W_n)$ is only the left-hand side of one particular DLR test function with $f=1$. It is an open and probably complex question, to derive conditions on a test function, say $f_l$, under which the hyperuniformity property of $\Gamma$ implies that $\mathrm{Var}_u[\widetilde{\mathrm{DLR}}_n(f_l, \Gamma; \theta^*)] = o(|W_n|)$. This question, of great importance from a statistical point of view, is left for future research.

**5. Simulation and numerical study.** In this section, we present simulation results and evaluate the performance of the estimator based on the Takacs-Fiksel method. The simulated point patterns were generated using a Metropolis-Hastings algorithm with a Gaussian move proposal and a Strauss interaction. The point patterns are generated on $[-30, 30]^2$ and then
clipped to the observation domain $[-\ell, \ell]^2$ (for $\ell = 8, 12, 16$). Replicated point patterns are generated from the same seed, for different sets of parameters $\theta_1 := \frac{1}{\sigma^2}$, where $\sigma$ is the standard deviation of the Gaussian move, $\theta_2$ is the strength of the Strauss interaction, and $R$ is the range of the Strauss interaction.

<figure><img src="image_3.png" /><figcaption>Fig 2: Simulated Gibbs perturbed lattices with Strauss interaction and Gaussian moves.</figcaption></figure>
In these simulations, we observe an interesting phenomenon regarding the effect of $\theta_2$ and $R$ on clustering behavior. For regular continuous Gibbs point processes with Strauss interaction, the larger the value of $\theta_2$, the more repulsive the interaction, and similarly, the larger $R$ is, the more repulsive the point process. For a fixed value of $\theta_1$, we observe that when $\theta_2$ is small (i.e., weak interaction), the influence of $R$ is negligible and the resulting point patterns are nearly indistinguishable. In this regime, the configurations display clusters of points separated by large empty regions, reminiscent of a Poisson point process. In contrast, when $\theta_2$ is large (i.e., strong interaction), the behavior changes markedly depending on the value of $R$. For $R = 0.5$, the configuration exhibits repulsive behavior: the points are more evenly spaced, with only a few pairs located close to one another. Conversely, for $R = 1.5$, a clear clustering phenomenon emerges, with groups of three to four interacting points occurring more frequently. Interestingly, these clusters themselves appear to be regularly spaced, as though they repel each other. According to Proposition 3.6, at least one Gibbs perturbed lattice with Strauss interaction and Gaussian moves is a hyperuniform point process. Figure 2 shows that we can obtain hyperuniform processes with different clustering behaviors simply by tuning the model parameters.

<table><caption>TABLE 1</caption><thead><tr><th rowspan="3">(θ<sub>1</sub>, θ<sub>2</sub>)</th><th colspan="9">W<sub>n</sub> =: [-ℓ, ℓ]<sup>2</sup></th></tr><tr><th colspan="3">ℓ = 8</th><th colspan="3">ℓ = 12</th><th colspan="3">ℓ = 16</th></tr><tr><th>mean</th><th>sd</th><th>RMSE</th><th>mean</th><th>sd</th><th>RMSE</th><th>mean</th><th>sd</th><th>RMSE</th></tr></thead><tbody><tr><td colspan="10"><b>R = 0.5</b></td></tr><tr><td>θ<sub>1</sub> = 4</td><td>3.88</td><td>0.19</td><td>0.23</td><td>4.01</td><td>0.16</td><td>0.16</td><td>4.10</td><td>0.18</td><td>0.21</td></tr><tr><td>θ<sub>2</sub> ≈ 0.11</td><td>0.11</td><td>0.19</td><td>0.19</td><td>0.03</td><td>0.09</td><td>0.12</td><td>0.12</td><td>0.11</td><td>0.11</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>1.10</td><td>0.07</td><td>0.12</td><td>0.99</td><td>0.07</td><td>0.07</td><td>1.00</td><td>0.05</td><td>0.05</td></tr><tr><td>θ<sub>2</sub> ≈ 0.11</td><td>0.15</td><td>0.09</td><td>0.10</td><td>0.07</td><td>0.11</td><td>0.12</td><td>0.11</td><td>0.06</td><td>0.06</td></tr><tr><td>θ<sub>1</sub> = 4</td><td>4.19</td><td>0.41</td><td>0.45</td><td>4.04</td><td>0.20</td><td>0.20</td><td>4.05</td><td>0.19</td><td>0.19</td></tr><tr><td>θ<sub>2</sub> ≈ 0.69</td><td>0.84</td><td>0.28</td><td>0.32</td><td>0.69</td><td>0.17</td><td>0.17</td><td>0.71</td><td>0.08</td><td>0.08</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>1.01</td><td>0.08</td><td>0.08</td><td>0.98</td><td>0.04</td><td>0.05</td><td>1.02</td><td>0.03</td><td>0.04</td></tr><tr><td>θ<sub>2</sub> ≈ 0.69</td><td>0.71</td><td>0.16</td><td>0.16</td><td>0.65</td><td>0.17</td><td>0.17</td><td>0.68</td><td>0.15</td><td>0.15</td></tr><tr><td>θ<sub>1</sub> = 4</td><td>3.98</td><td>0.37</td><td>0.37</td><td>3.93</td><td>0.15</td><td>0.17</td><td>4.01</td><td>0.15</td><td>0.15</td></tr><tr><td>θ<sub>2</sub> ≈ 2.3</td><td>2.16</td><td>0.30</td><td>0.33</td><td>2.33</td><td>0.20</td><td>0.21</td><td>2.50</td><td>0.25</td><td>0.32</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>0.99</td><td>0.07</td><td>0.07</td><td>1.01</td><td>0.05</td><td>0.05</td><td>0.98</td><td>0.03</td><td>0.04</td></tr><tr><td>θ<sub>2</sub> ≈ 2.3</td><td>2.30</td><td>0.52</td><td>0.52</td><td>2.26</td><td>0.23</td><td>0.23</td><td>2.24</td><td>0.18</td><td>0.19</td></tr><tr><td colspan="10"><b>R = 1.5</b></td></tr><tr><td>θ<sub>1</sub> = 4</td><td>4.03</td><td>0.50</td><td>0.50</td><td>4.16</td><td>0.30</td><td>0.34</td><td>4.06</td><td>0.27</td><td>0.27</td></tr><tr><td>θ<sub>2</sub> ≈ 0.11</td><td>0.11</td><td>0.06</td><td>0.07</td><td>0.11</td><td>0.07</td><td>0.07</td><td>0.13</td><td>0.05</td><td>0.06</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>1.06</td><td>0.10</td><td>0.12</td><td>1.05</td><td>0.07</td><td>0.09</td><td>0.99</td><td>0.03</td><td>0.03</td></tr><tr><td>θ<sub>2</sub> ≈ 0.11</td><td>0.11</td><td>0.09</td><td>0.09</td><td>0.10</td><td>0.06</td><td>0.06</td><td>0.11</td><td>0.04</td><td>0.04</td></tr><tr><td>θ<sub>1</sub> = 4</td><td>4.16</td><td>0.27</td><td>0.31</td><td>4.07</td><td>0.31</td><td>0.32</td><td>4.00</td><td>0.18</td><td>0.18</td></tr><tr><td>θ<sub>2</sub> ≈ 0.69</td><td>0.77</td><td>0.09</td><td>0.11</td><td>0.76</td><td>0.05</td><td>0.08</td><td>0.69</td><td>0.05</td><td>0.05</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>1.09</td><td>0.09</td><td>0.12</td><td>1.01</td><td>0.05</td><td>0.05</td><td>0.99</td><td>0.06</td><td>0.06</td></tr><tr><td>θ<sub>2</sub> ≈ 0.69</td><td>0.65</td><td>0.10</td><td>0.11</td><td>0.70</td><td>0.07</td><td>0.08</td><td>0.74</td><td>0.03</td><td>0.06</td></tr><tr><td>θ<sub>1</sub> = 4</td><td>92.97</td><td>49.08</td><td>101.61</td><td>6.27</td><td>0.47</td><td>2.32</td><td>5.33</td><td>0.61</td><td>1.46</td></tr><tr><td>θ<sub>2</sub> ≈ 2.3</td><td>101.57</td><td>55.77</td><td>113.86</td><td>3.21</td><td>0.46</td><td>1.01</td><td>3.43</td><td>0.37</td><td>1.19</td></tr><tr><td>θ<sub>1</sub> = 1</td><td>1.21</td><td>0.12</td><td>0.24</td><td>1.02</td><td>0.07</td><td>0.07</td><td>1.10</td><td>0.02</td><td>0.10</td></tr><tr><td>θ<sub>2</sub> ≈ 2.3</td><td>2.60</td><td>0.36</td><td>0.46</td><td>2.38</td><td>0.22</td><td>0.24</td><td>2.43</td><td>0.14</td><td>0.19</td></tr></tbody></table>

Mean, standard deviation and root-mean squared error (RMSE) of estimates based on 200 simulations from Gibbs perturbed lattices with Strauss interaction and Gaussian moves, and observed on [-ℓ, ℓ]<sup>2</sup> for ℓ = 8, 12, 16.

For each parameter configuration, 200 point patterns are simulated, and the parameters are estimated for each realization using the Takacs-Fiksel procedure with the edge corrections
described in Section 4. Due to edge correction effects, when $W_n = [-\ell, \ell]^2$ with $\ell = 8, 12, 16$, the estimator uses on average 100, 256, and 576 points, respectively. As shown in Table 1, the Takacs-Fiksel estimator performs excellently well. Both estimates of $\theta_1$ and $\theta_2$ exhibit very small biases, with standard deviations decreasing as $\ell$ increases. Consequently, the RMSEs also decrease with $\ell$ and remain low even for the smallest observation domain $[-8, 8]^2$. One notable exception occurs for the parameters $\theta_1 = 4$, $\theta_2 = -\log 0.1 \approx 2.3$, and $R = 1.5$, corresponding to the clustering regime described earlier with a more constrained move. In this regime, the estimator performs poorly: enlarging the observation window improves performance, but it still fails to reach the estimation accuracy observed in other parameter regimes. The system is highly rigid, with strong repulsion between particles and a strong attraction to the lattice sites. This competition between antagonistic forces may explain the failure of the estimation procedure. Further tests on larger observation windows are necessary, but they are computationally expensive due to the cost of generating point patterns and performing the optimisation. For weaker interactions (i.e., $\theta_2 = -\log 0.9 \approx 0.11$), the estimation procedure is consistently reliable, regardless of $\theta_1$ and $R$. In most parameter regimes, the RMSEs decay at a rate consistent with $|W_n|^{-1/2}$.

## APPENDIX A: AUXILIARY RESULTS FOR THE PROOF OF THEOREM 1

In this section, we prove auxiliary results needed for the proof of Theorem 1, which is that under assumptions (H1)-(H4) and (H5.1)-(Q1)[1] or (H5.2)-(Q2), $P_{\mathcal{L}}(H, \mathbb{Q}) \neq \emptyset$. Let $\Lambda_n = \mathcal{L} \cap [-n, n]^d$, we define the empirical field $\bar{P}_n$ as the probability measure

$$ \bar{P}_n := \frac{1}{|\Lambda_n|} \sum_{i \in \Lambda_n} \hat{P}_n(\cdot - i) $$

where $\hat{P}_n = \bigotimes_{k \in \mathcal{L}} \mathbb{P}_{\Lambda_n + 2kn/|k|}$. The empirical field is invariant under translation in $\mathcal{L}$. Our aim is to prove that the empirical field exhibits an accumulation point under the local topology. The local topology is defined as the coarsest topology such that for any cylinder event, also called local event, A the evaluation map $P \mapsto P(A)$ is continuous. The main tightness tool is the specific entropy. Let P be an invariant under translation probability measure on $\Omega$,

$$ (A.1) \qquad I_s(P) = \lim_{\Lambda \to \mathcal{L}} \frac{1}{|\Lambda|} I(P|_{\Lambda} | \mathbb{Q}^{\otimes \Lambda}) $$

where $I(\cdot|\cdot)$ is the Kullback-Leibler divergence and this limit exists due to a sub-additivity argument. The following proposition gives the properties of the specific entropy and the arguments that give the tightness of probability measures of lattice systems.

**PROPOSITION A.1 (Theorem 15.12 and Proposition 15.14 in [20]).**

For any probability measure P on $\Omega$ invariant under any translation on $\mathcal{L}$ we define the specific entropy of P as,

$$ (A.2) \qquad I_s(P) = \sup_{\Lambda \subset \mathcal{L}, |\Lambda| < \infty} \frac{1}{|\Lambda|} I(P|_{\Lambda} | \mathbb{Q}^{\otimes \Lambda}). $$

Therefore, the specific entropy does not depend on the sequence of finite sub-lattices $(\Lambda_n)_{n \in \mathbb{N}}$. Moreover, the specific entropy, $I_s$, is affine and the level sets $\{I_s(\cdot) \le C\}$ are compact and sequentially compact for the local topology.

Now we prove that the specific entropy for the empirical fields is bounded in order to exhibit a candidate for a Gibbs measure on the lattice.
**PROPOSITION A.2.** The empirical field $(\bar{\mathbb{P}}_n)_{n \in \mathbb{N}}$ has an accumulation point in the local topology that we denote by $\mathbb{P}_\infty$.

**PROOF.** Since the specific entropy is affine, we have that

$$I_s(\bar{\mathbb{P}}_n) = \frac{1}{|\Lambda_n|} \sum_{i \in \Lambda_n} I_s(\hat{\mathbb{P}}_n(\cdot - i)).$$

Furthermore, since the specific entropy does not depend on the sequence of sub-lattices, then for $i \in \Lambda_n$ we choose the sequence $(\Lambda_{(2k-1)n} + i)_{k \ge 1}$ and therefore

$$ 
\begin{align*}
I_s(\hat{\mathbb{P}}_n(\cdot - i)) &= \lim_{k \to \infty} \frac{1}{(2k-1)^d |\Lambda_n|} I(\hat{\mathbb{P}}_{n|\Lambda_{(2k-1)n}} | \mathbb{Q}^{\otimes \Lambda_{(2k-1)n}}) \\
&= \lim_{k \to \infty} \frac{1}{(2k-1)^d |\Lambda_n|} (2k-1)^d I(\mathbb{P}_n | \mathbb{Q}^{\otimes \Lambda_n}) \\
&= \frac{1}{|\Lambda_n|} I(\mathbb{P}_n | \mathbb{Q}^{\otimes \Lambda_n}).
\end{align*}
$$

Combining everything we have

$$I_s(\bar{\mathbb{P}}_n) = \frac{1}{|\Lambda_n|} I(\mathbb{P}_n | \mathbb{Q}^{\otimes \Lambda_n}).$$

Finally, according to (H1) and (H2), for $\mathbf{X} \sim \mathbb{P}_n$, we have

$$ 
\begin{align*}
I(\mathbb{P}_n | \mathbb{Q}^{\otimes \Lambda_n}) &= -\log Z_{\Lambda_n} - \mathbb{E}(\tilde{H}(\mathbf{X})) \\
&\le (A + B - \log \mathbb{Q}(B(0, \varepsilon)))|\Lambda_n|.
\end{align*}
$$

As a consequence, for all $n \in \mathbb{N}$

$$(A.3) \qquad I_s(\bar{\mathbb{P}}_n) \le A + B - \log \mathbb{Q}(B(0, \varepsilon)),$$

and thus, according to Proposition A.1, $(\bar{\mathbb{P}}_n)_{n \in \mathbb{N}}$ has an accumulation point. $\square$

A function $f$ is local if there is $\Delta \subset \mathcal{L}$ finite such that for any $\mathbf{x} \in \mathbb{Q}^\mathcal{L}$, $f(\mathbf{x}) = f(\mathbf{x}_\Delta)$. Moreover, $f$ is tame if $\mathbb{E}(e^{f(\tilde{\mathbf{X}}_\Delta)}) < \infty$ for $\tilde{\mathbf{X}}_\Delta \sim \mathbb{Q}^{\otimes \Delta}$.

**LEMMA A.3.** Let $f$ be a local and tame function on $\Delta \subset \mathcal{L}$ finite, for any stationary probability measure $P$ on $\Omega$, $\mathbf{X} \sim P$ and $\tilde{\mathbf{X}}_\Delta \sim \mathbb{Q}^{\otimes \Delta}$ we have

$$\mathbb{E}[f(\mathbf{X})] \le I_s(P)|\Delta| + \log \mathbb{E}[e^{f(\tilde{\mathbf{X}}_\Delta)}].$$

**PROOF.** By Jensen's inequality we have that

$$(A.4) \qquad \mathbb{E}_P[f(\mathbf{X})] = \mathbb{E}[f(\mathbf{X}_\Delta)] \le I(P_\Delta | \mathbb{Q}^{\otimes \Delta}) + \log \mathbb{E}[e^{f(\tilde{\mathbf{X}}_\Delta)}].$$

On the other hand, we know that

$$I_s(P) = \sup_{\Lambda \subset \mathcal{L}, |\Lambda| < \infty} \frac{1}{|\Lambda|} I(P_\Lambda | \mathbb{Q}^{\otimes \Lambda}).$$

Consequently, we have

$$(A.5) \qquad I(P_\Delta | \mathbb{Q}^{\otimes \Delta}) \le I_s(P)|\Delta|.$$

By combining (A.4) and (A.5), we obtain the desired inequality. $\square$
APPENDIX B: PROOF OF THEOREM 1

PROOF. We split the proofs into two parts. The first two show that the measure $P_\infty$ satisfies DLR equations under the two different sets of assumptions.

$\to \text{PART A. } P_\infty$ SATISFIES DLR EQUATIONS UNDER THE ASSUMPTIONS (H1)-(H4) AND (H5.1)-(Q1)[1]. Let $f$ be a bounded local function, $\Delta \subset \mathcal{L}$ finite and $\epsilon > 0$. We split the proof into 7 steps.

**Step 1.** Since $f$ is a local bounded function by convergence in the local topology, we have for $n$ large enough

$$ (B.1) \qquad \left| \int f d\mathbb{P}_\infty - \int f d\bar{\mathbb{P}}_n \right| \le \epsilon \|f\|_\infty. $$

**Step 2.** We define $B_k = [-k, k]^d \oplus \Delta$, $\Delta_k = \Lambda_k \cap \mathcal{L}$ and we consider the following events

$$ \begin{align*} I_{\Delta,l} &:= \{\exists i \in \Delta, i + x_i \in B_l^c\}, \\ T_{\Delta,k,l,R} &= \{\exists j \in \Delta_k^c, j + x_j \in B_{l+R}\}, \\ E_{k,l,R,\Delta} &= I_{\Delta,l} \cup T_{\Delta,k,l,R}. \end{align*} $$

We have that

$$ (B.2) \qquad \sup_{n \in \mathbb{N}} \left| \int f d\bar{\mathbb{P}}_n - \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n \right| \le \|f\|_\infty \epsilon. $$

Using the fact that $(X_i)_{i \in \mathcal{L}}$ are identically distributed under $\bar{\mathbb{P}}_n$ for any $n \ge 1$ and Markov inequality

$$ \begin{align*} \bar{\mathbb{P}}_n(E_{k,l,R,\Delta}) &\le \sum_{i \in \Delta} \bar{\mathbb{P}}_n(i + X_i \in B_l^c) + \sum_{j \in \Delta_k} \bar{\mathbb{P}}_n(j + X_j \in B_{l+R}) \\ &\le |\Delta| \bar{\mathbb{P}}_n(|X_0| \ge l) + |\Lambda_{l+R}| \bar{\mathbb{P}}_n(|X_0| \ge k - l - R) \\ &\le \left( \frac{|\Delta|}{l} + \frac{|B_{l+R}|}{k - l - R} \right) \mathbb{E}(|X_0|). \end{align*} $$

Furthermore, using Lemma A.3 and assumption (Q1)[1] and (A.3), there exists $c > 0$ such that for any $n \in \mathbb{N}$

$$ \mathbb{E}[c|X_0|] \le I_s(\bar{\mathbb{P}}_n) + \log \mathbb{E}[e^{|\tilde{X}|}] < \infty $$

for $\tilde{X} \sim \mathbb{Q}$. Consequently, there exists $C > 0$ such that

$$ (B.3) \qquad \bar{\mathbb{P}}_n(E_{k,l,R,\Delta}) \le C \left( \frac{|\Delta|}{l} + \frac{|B_{l+R}|}{k-l-R} \right). $$

As a consequence for $k$ and $l$ large enough we obtain (B.2).

**Step 3.** For $n$ large enough we have,

$$ (B.4) \qquad \left| \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n - \int \tilde{f}_{\Delta,l} \mathbb{1}_{T_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n \right| \le 2\epsilon \|f\|_\infty, $$

where

$$ \tilde{f}_{\Delta,l}(\mathbf{x}_{\Delta^c}) = \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \int f(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c}) \mathbb{1}_{I_{l,\Delta}^c}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} \mathbb{Q}^{\otimes \Delta}. $$
We consider $n$ large enough such that $\Delta \subset \Lambda_n$ and we define $\Lambda_n^* := \{i \in \Lambda_n, \Delta - i \subset \Lambda_n\}$. By definition, we have

$$ \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n = \frac{1}{|\Lambda_n|} \sum_{i \in \Lambda_n} \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} (\cdot + i) d\hat{\mathbb{P}}_n, $$

and as a consequence we have

$$ \left| \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n - \frac{1}{|\Lambda_n|} \sum_{i \in \Lambda_n^*} \int f \mathbb{1}_{E_{k,l,R,\Delta}^c} (\cdot + i) d\hat{\mathbb{P}}_n \right| \le \frac{|\Lambda_n \setminus \Lambda_n^*|}{|\Lambda_n|} \|f\|_{\infty}. $$

Let us consider

$$ \tilde{f}_{\Delta, E^c}(\mathbf{x}) = \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \int f \mathbb{1}_{E_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta^c} \cup \mathbf{x}'_{\Delta}) e^{-H_{\Delta}(\mathbf{x}_{\Delta^c} \cup \mathbf{x}'_{\Delta})} Q^{\otimes \Delta} $$

For $i \in \Lambda_n^*$, we have

$$ \int \tilde{f}_{\Delta, E^c}(\cdot + i) d\hat{\mathbb{P}}_n = \iint f \mathbb{1}_{E_{k,l,R,\Delta-i}^c} (\mathbf{x}_{(\Delta-i)^c} \cup \mathbf{x}'_{\Delta-i}) \frac{e^{-\tilde{H}_{\Delta-i}(\mathbf{x}_{(\Delta-i)^c} \cup \mathbf{x}'_{\Delta-i})}}{Z_{\Delta-i}(\mathbf{x}_{(\Delta-i)^c})} Q^{\otimes (\Delta-i)} d\hat{\mathbb{P}}_n. $$

Since $\Delta - i \subset \Lambda_n$ and that $\hat{\mathbb{P}}_n$ satisfies the DLR equations for $\Delta + i$, we have

$$ \int \tilde{f}_{\Delta, E^c}(\cdot + i) d\hat{\mathbb{P}}_n = \int f \mathbb{1}_{E_{k,l,R,\Delta}^c}(\cdot + i) d\hat{\mathbb{P}}_n $$

Similarly we have

$$ \left| \int \tilde{f}_{\Delta, E^c} d\bar{\mathbb{P}}_n - \frac{1}{|\Lambda_n|} \sum_{i \in \Lambda_n^*} \int \tilde{f}_{\Delta, E^c}(\cdot + i) d\hat{\mathbb{P}}_n \right| \le \frac{|\Lambda_n \setminus \Lambda_n^*|}{|\Lambda_n|} \|f\|_{\infty}. $$

Finally, we have that $\mathbb{1}_{E_{k,l,R,\Delta}^c} = \mathbb{1}_{I_{l,\Delta}^c} \mathbb{1}_{T_{k,l,R,\Delta}^c}$ and since $T_{k,l,R,\Delta}^c$ is an event that depends only on the configuration outside $\Delta$ we have

$$ \tilde{f}_{\Delta, E^c} = \tilde{f}_{\Delta, l} \mathbb{1}_{T_{k,l,R,\Delta}^c}. $$

Therefore, by combining everything above we obtain (B.4) for $n$ large enough.

**Step 4.** In this step we consider the following

$$ \hat{f}_{\Delta,l}(\mathbf{x}) = \frac{1}{Z_{\Delta,l}(\mathbf{x}_{\Delta^c})} \int f(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c}) \mathbb{1}_{I_{l,\Delta}^c}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} Q^{\otimes \Delta} $$

where

$$ Z_{\Delta,l}(\mathbf{x}_{\Delta^c}) = \int \mathbb{1}_{I_{l,\Delta}^c}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} Q^{\otimes \Delta}. $$

For $k, l$ and $n$ large enough we have,

$$ (B.5) \quad \left| \int \tilde{f}_{\Delta,l}(\mathbf{x}_{\Delta^c}) \mathbb{1}_{T_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta_k^c}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta,l}(\mathbf{x}_{\Delta^c}) \mathbb{1}_{T_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta_k^c}) d\hat{\mathbb{P}}_n \right| \le 3\epsilon \|f\|_{\infty}. $$

Since the difference between the normalisation constants is

$$ Z_{\Delta}(\mathbf{x}_{\Delta^c}) - Z_{\Delta,l}(\mathbf{x}_{\Delta^c}) = \int \mathbb{1}_{I_{l,\Delta}^c}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} Q^{\otimes \Delta}, $$
we have the following upper bound

$$ \left| \int \tilde{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n \right| \le \|f\|_{\infty} \times \\ \qquad \qquad \qquad \qquad \qquad \qquad \int \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \int \mathbf{1}_{I_{l,\Delta}}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} Q^{\otimes \Delta} d\bar{\mathbb{P}}_n $$

Furthermore, for $n$ large enough we can apply (B.4) for $f \equiv 1$, and obtain that

$$ \int \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \int \mathbf{1}_{I_{l,\Delta}}(\mathbf{x}'_{\Delta}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} Q^{\otimes \Delta} d\bar{\mathbb{P}}_n \le 2\epsilon + \bar{\mathbb{P}}_n(E_{k,l,R,\Delta}), $$

and by (B.3) we know that for $k, l$ large enough we have,

$$ \sup_{n \in \mathbb{N}} \bar{\mathbb{P}}_n(E_{k,l,R,\Delta}) \le \epsilon, $$

and therefore we have (B.5).

Step 5. For $k, l$ and $n$ large enough we have

$$ (B.6) \quad \left| \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta^c}) \mathbf{1}_{T_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta_k^c}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta^c}) \mathbf{1}_{T_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta_k^c}) d\bar{\mathbb{P}}_{\infty} \right| \le 3\epsilon \|f\|_{\infty}. $$

We can remark that under the event $T_{k,l,R,\Delta}^c$, $\hat{f}_{\Delta, l}$ depends only on the configurations in $\Delta_k \setminus \Delta$ since the points originating from $\Delta_k^c$ are too far away to interact with the points from $\Delta$. Therefore, we have

$$ \begin{aligned} & \left| \int \hat{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c} dP_{\infty} \right| \\ & \le \left| \int \hat{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n \right| \\ & \qquad + \left| \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) dP_{\infty} \right| \\ & \qquad + \left| \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) dP_{\infty} - \int \hat{f}_{\Delta, l} \mathbf{1}_{T_{k,l,R,\Delta}^c} dP_{\infty} \right| \\ & \le \|f\|_{\infty} (\bar{\mathbb{P}}_n(T_{k,l,R,\Delta}) + P_{\infty}(T_{k,l,R,\Delta})) \\ & \qquad + \left| \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) dP_{\infty} \right|. \end{aligned} $$

Using similar arguments as the ones used in Step 2, there exists a constant $C > 0$ such that

$$ \sup_{n \in \mathbb{N}} \bar{\mathbb{P}}_n(T_{k,l,R,\Delta}) \le \frac{C|B_{l+R}|}{k-l-R}. $$

As a consequence for $k$ large enough we have

$$ \bar{\mathbb{P}}_n(T_{k,l,R,\Delta}) + P_{\infty}(T_{k,l,R,\Delta}) \le 2\epsilon. $$

Once $k$ is fixed, by convergence in the local topology we have for $n$ large enough

$$ \left| \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n - \int \hat{f}_{\Delta, l}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_{\infty} \right| \le \epsilon \|f\|_{\infty}. $$
Step 6. For $l$ large enough we have

$$ (B.7) \qquad \left| \int \hat{f}_{\Delta,l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\mathbb{P}_{\infty} - \int \tilde{f}_{\Delta,l} d\mathbb{P}_{\infty} \right| \le 2\epsilon \|f\|_{\infty}. $$

Following similar computations as in Step 4 we show that

$$ \left| \int \hat{f}_{\Delta,l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\mathbb{P}_{\infty} - \int \tilde{f}_{\Delta,l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\mathbb{P}_{\infty} \right| \le \|f\|_{\infty} \int \int \frac{e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})}}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \mathbf{1}_{I_{l,\Delta}} \mathbb{Q}^{\otimes \Delta} d\mathbb{P}_{\infty}. $$

By dominated convergence theorem, we have for $l$ large enough

$$ \iint \frac{e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})}}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \mathbf{1}_{I_{l,\Delta}} \mathbb{Q}^{\otimes \Delta} d\mathbb{P}_{\infty} \le \epsilon. $$

Furthermore, for $k$ large enough we have

$$ \left| \int \tilde{f}_{\Delta,l} \mathbf{1}_{T_{k,l,R,\Delta}^c} d\mathbb{P}_{\infty} - \int \tilde{f}_{\Delta,l} d\mathbb{P}_{\infty} \right| \le \|f\|_{\infty} \mathbb{P}_{\infty}(T_{k,l,R,\Delta}) \le \epsilon \|f\|_{\infty}, $$

and this finishes the proof (B.7).

Step 7. Finally by dominated convergence theorem, for $l$ large enough we have

$$ (B.8) \qquad \left| \int \tilde{f}_{\Delta,l} d\mathbb{P}_{\infty} - \int \hat{f}_{\Delta} d\mathbb{P}_{\infty} \right| \le \epsilon \|f\|_{\infty}. $$

where

$$ \tilde{f}_{\Delta}(\mathbf{x}) = \frac{1}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \int f(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c}) e^{-\tilde{H}_{\Delta}(\mathbf{x}'_{\Delta} \cup \mathbf{x}_{\Delta^c})} \mathbb{Q}^{\otimes \Delta}. $$

Conclusion. By combining (B.1)-(B.8), we have

$$ \left| \int f d\mathbb{P}_{\infty} - \int \tilde{f}_{\Delta} d\mathbb{P}_{\infty} \right| \le 13\epsilon \|f\|_{\infty}. $$

This inequality, being true for all $\epsilon > 0$, shows that $P_{\infty}$ satisfies the DLR equations for bounded local functions. Then, these DLR equations can be extended to bounded measurable functions by a monotone class argument, which ends this first part.

→ PART B. $P_{\infty}$ SATISFIES DLR EQUATIONS UNDER THE ASSUMPTIONS (H1)-(H4) AND (H5.2)-(Q2). We consider the same notation and 7 steps as for the first part and detail only those that differ.

Step 2.bis. For $n$ large enough, we have

$$ (B.9) \qquad \left| \int f d\bar{\mathbb{P}}_n - \int \tilde{f}_{\Delta} d\bar{\mathbb{P}}_n \right| \le 2\epsilon \|f\|_{\infty}. $$

This step follows from the same arguments as in Step 2 where we replace $f \mathbf{1}_{E_{k,l,R,\Delta}^c}$ with $f$.

Step 3.bis. For $k > 0$ large enough, we have

$$ (B.10) \qquad \sup_{n \in \mathbb{N}} \left| \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta^c}) d\bar{\mathbb{P}}_n - \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n \right| \le 3\epsilon \|f\|_{\infty}. $$

By definition, the difference is equal to

$$ \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta^c}) - \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n = \iint f(\mathbf{x}) \frac{e^{-\tilde{H}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta} \cup \mathbf{x}'_{\Delta})}}{Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta})} \times \\ \left( \frac{e^{-\sum_{i \in \Delta} \sum_{j \in \Delta_k^c} \phi(|i+x'_i-j-x_j|)} Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) - Z_{\Delta}(\mathbf{x}'_{\Delta})}{Z_{\Delta}(\mathbf{x}'_{\Delta})} \right) d\mathbb{Q}^{\otimes \Delta} d\bar{\mathbb{P}}_n. $$
Since $\mathcal{Q}$ is bounded by (Q2) and by monotonicity of $\phi$ at infinity, we have that for $k$ large enough, any $x'_i, x_j \in \mathcal{Q}$

$$ \sum_{i \in \Delta} \sum_{j \in \Delta_k^c} |\phi(|i+x'_i-j-x_j|)| \le \sum_{i \in \Delta} \sum_{j \in \Delta_k^c} |\phi(|i-j|-D)|, $$

where $D = \text{diam}(\mathcal{Q})$. According to assumption (H5.2), we know that $\phi$ is summable and thus for any $\eta > 0$, there exists $k$ large enough such that

$$ \sum_{i \in \Delta} \sum_{j \in \Delta_k^c} |\phi(|i+x'_i-j-x_j|)| \le \eta. $$

Therefore for $k$ large enough we have

$$ \begin{aligned} & \left| \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta^c}) d\bar{\mathbb{P}}_n - \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n \right| \\ & \le \|f\|_{\infty} \int \int \frac{e^{-\tilde{H}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta} \cup \mathbf{x}'_{\Delta})}}{Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta})} \times \\ & \qquad \left( \left| \frac{(e^{\eta} - 1) Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta})}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \right| + \left| \frac{Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) - Z_{\Delta}(\mathbf{x}_{\Delta^c})}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \right| \right) d\mathbb{Q}^{\otimes \Delta} d\bar{\mathbb{P}}_n \\ & \le \|f\|_{\infty} \int \left( \left| \frac{(e^{\eta} - 1) Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta})}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \right| + \left| \frac{Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) - Z_{\Delta}(\mathbf{x}_{\Delta^c})}{Z_{\Delta}(\mathbf{x}_{\Delta^c})} \right| \right) d\mathbb{Q}^{\otimes \Delta} d\bar{\mathbb{P}}_n. \end{aligned} $$

Furthermore, for $k$ large enough we have,

$$ \begin{aligned} Z_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) &= \int e^{-\tilde{H}_{\Delta}(\mathbf{x}_{\Delta^c} \cup \mathbf{x}'_{\Delta})} + \sum_{i \in \Delta} \sum_{j \in \Delta_k^c} \phi(|i+x'_i-j-x_j|) d\mathbb{Q}^{\otimes \Delta} \\ &\le e^{\eta} Z_{\Delta}(\mathbf{x}_{\Delta^c}), \end{aligned} $$

and thus

$$ \left| \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta^c}) d\bar{\mathbb{P}}_n - \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n \right| \le \|f\|_{\infty} (e^{\eta} - 1)(1 + e^{\eta}). $$

Hence by choosing $\eta > 0$ small enough and $k$ large enough we obtain (B.10).

**Step 4.bis.** For any $k$ fixed, by convergence in the local topology we have for $n$ large enough

$$ (B.11) \qquad \left| \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_n - \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_{\infty} \right| \le \epsilon \|f\|_{\infty}. $$

**Step 5.bis.** By following the same argument as in Step 3.bis, we have for $k$ large enough

$$ (B.12) \qquad \left| \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta^c}) d\bar{\mathbb{P}}_{\infty} - \int \tilde{f}_{\Delta}(\mathbf{x}_{\Delta_k \setminus \Delta}) d\bar{\mathbb{P}}_{\infty} \right| \le 3\epsilon \|f\|_{\infty}. $$

**Conclusion.** By combining (B.1), (B.9), (B.10), (B.11) and (B.12), we have

$$ \left| \int f d\bar{\mathbb{P}}_{\infty} - \int \tilde{f}_{\Delta} d\bar{\mathbb{P}}_{\infty} \right| \le 10\epsilon \|f\|_{\infty} $$

and we conclude in the same way as for the first part.

Let $\Gamma$ be a stationary point of $\Gamma$ and let $U = u$. Parts A-B prove that $G_{\mathcal{L}}(H, \mathbb{Q}) \neq \emptyset$ and as a consequence we have that $P_{\mathcal{L}}(H, \mathbb{Q}) \neq \emptyset$. All that is left is to prove that the Gibbs perturbed lattice $\Gamma$ is a stationary point process. We start with a first observation, by definition
for $\mathbb{P} \in \mathcal{G}_{\mathcal{L}}(H, \mathbb{Q})$ and $\mathbf{X} \sim \mathbb{P}$ the law of $\mathbf{X}$ is invariant under translation and therefore for all $i \in \mathcal{L}$, $X_i \stackrel{d}{=} X_0$. We consider $W \subset \mathbb{R}^d$ bounded, since $(X_i)_{i \in \mathcal{L}}$ are identically distributed we have that

$$ 
\begin{aligned}
\mathbb{E}[\#(\Gamma \cap W)] &= \sum_{i \in \mathcal{L}} \mathrm{P}(i + X_i + U \in W) \\
&\leq \sum_{i \in \mathcal{L}} \mathrm{P}(X_i \in (W \oplus \mathcal{D}) - i) \\
&\leq \sum_{i \in \mathcal{L}} \mathrm{P}(X_0 \in (W \oplus \mathcal{D}) - i) \\
&\leq |W \oplus \mathcal{D}|.
\end{aligned}
$$

Since $\mathbb{E}(\#\Gamma \cap W)$ is finite then almost surely $\#\Gamma \cap W$ is finite and thus $\Gamma$ is almost surely locally finite. Now we consider $t \in \mathbb{R}^d$, then there is a unique $k \in \mathcal{L}$ and $V \in \mathcal{D}$ such that $t+U=k+V$ and we know that $V \sim U_{\mathcal{D}}$. Therefore, for any bounded measurable function $f$ on $\mathcal{C}$ we have

$$ 
\begin{aligned}
\mathbb{E}[f(\Gamma + t)] &= \mathbb{E}[f(\Gamma(\mathbf{X}, t+U))] \\
&= \mathbb{E}[f(\Gamma(\mathbf{X}, k+V))] \\
&= \mathbb{E}[f(\Gamma(\mathbf{X}, V))] = \mathbb{E}[f(\Gamma)].
\end{aligned}
$$

Therefore, for any $t \in \mathbb{R}^d$ we have $\Gamma = \Gamma + t$. Since $\Gamma$ is almost surely locally finite and its distribution is invariant under translations of $\mathbb{R}^d$, it follows that $\Gamma$ is a stationary point process. Finally, for any $u \in \mathcal{D}$ and $k \in \mathcal{L}$, we also have

$$ \mathbb{E}_u[f(\Gamma + k)] = \mathbb{E}_u[f(\Gamma(\mathbf{X}, k + u))] = \mathbb{E}_u[f(\Gamma(\mathbf{X}, u))] = \mathbb{E}_u[f(\Gamma)]. $$

Therefore, for any $u \in \mathcal{D}$ the law of $\Gamma | U = u$ is invariant under translation by any element of $\mathcal{L}$. $\square$

## APPENDIX C: PROOF OF PROPOSITION 3.6

**PROOF.** (i) Let $\mathbb{P} \in \mathcal{L}(H, \mathbb{Q})$. According to Lemma A.3 we have that

$$ (C.1) \qquad \mathbb{E}[|X_0|^d] \leq c (I_s(\mathbb{P}) + \log \mathbb{E}[e^{c|X|^d}]) $$

for $X \sim \mathbb{Q}$. It is not clear whether for any $\mathbb{P} \in \mathcal{PL}(H, \mathbb{Q})$, $I_s(\mathbb{P}) < \infty$. However we know by (A.3) that we have at least one measure $\mathbb{P}_\infty$ with a finite specific entropy. Hence, $\Gamma \sim \mathbb{P}_\infty$ is hyperuniform.

(ii) By DLR equations, for any $\mathbb{P} \in \mathcal{PL}(H, \mathbb{Q})$, we have

$$ \mathbb{E}[|X_0|^d] = \mathbb{E} \left[ \int |x|^d \frac{e^{-h(x+U, \Gamma_{0c})}}{Z_0(\Gamma_{0c})} \mathbb{Q}(dx) \right] $$

where $\Gamma_{0c} = \{i + X_i + U, i \in \mathcal{L} \setminus \{0\}\}$. Under the assumption (H5.3), there exists $M \geq 0$ such that $|h| \leq M$ and thus we have

$$ \mathbb{E}[|X_0|^d] \leq e^{2M} \mathbb{E}[|X|^d] < \infty $$

for $X \sim \mathbb{Q}$. $\square$

## APPENDIX D: PROOFS OF RESULTS RELATED TO DLR TYPE EQUATIONS

This Appendix focuses on the proofs of Propositions 3.7-3.9.
### D.1. Proofs of Proposition 3.7.

PROOF. Proposition 3.7(i) follows directly from summing (3.8) over all sites. Statement (ii) is then a consequence obtained by using the stationarity of $\Gamma | U = u$ and the invariance by translation of $f$. $\square$

### D.2. Proof of Proposition 3.8.

PROOF. This result is obtained by applying Proposition 3.7 twice. $\square$

### D.3. Proof of Proposition 3.9.

PROOF. Let $u \in \mathcal{D}$ and $i \in \mathcal{L}_u$. From (3.8) we have

$$ \mathbb{E}_u [\nabla \check{f}(i + X_i, \mathbf{\Gamma}_{ic})] = \mathbb{E}_u \left[ \frac{1}{Z_i(\mathbf{\Gamma}_{ic})} \int_Q \nabla \check{f}(i + x, \mathbf{\Gamma}_{ic}) \Lambda(i, x, \mathbf{\Gamma}_{ic}) dx \right]. $$

Using an integration by parts on the right-hand side term, we have

$$ \begin{aligned} \mathbb{E}_u [\nabla \check{f}(i + X_i, \mathbf{\Gamma}_{ic})] ={}& \mathbb{E}_u \left[ \frac{1}{Z_i(\mathbf{\Gamma}_{ic})} \oint_{\partial Q} \check{f}(i + x, \mathbf{\Gamma}_{ic}) e^{-h(i+x, \mathbf{\Gamma}_{ic}) - q(x)} dx \right] \\ & + \mathbb{E}_u \left[ \frac{1}{Z_i(\mathbf{\Gamma}_{ic})} \int_Q \check{f}(i + x, \mathbf{\Gamma}_{ic}) (\nabla h(i + x, \mathbf{\Gamma}_{ic}) + \nabla q(x)) e^{-h(i+x, \mathbf{\Gamma}_{ic}) - q(x)} dx \right]. \end{aligned} $$

Since $\check{f}$ belongs to $\mathcal{F}_R(H, \mathbb{Q}, \mathcal{L}_u)$, the first term of the previous equation equals 0. Then, by applying the DLR equations to the second term, we obtain

$$ \mathbb{E}_u [\nabla \check{f}(i + X_i, \mathbf{\Gamma}_{ic})] = \mathbb{E}_u [\check{f}(i + X_i, \mathbf{\Gamma}_{ic}) (\nabla h(i + X_i, \mathbf{\Gamma}_{ic}) + \nabla q(X_i))] $$

The result is proved by summing over all sites of $\mathcal{L}_u$ and by applying Fubini's theorem, justified by (3.14). $\square$

## APPENDIX E: PROOF OF THEOREM 2

The proof of this result is mainly based on the following ergodic theorem obtained by [29], which we recall here for the sake of self-consistency.

THEOREM 4 (Ergodic theorem [29]). Let $\mathcal{L}$ be a full rank lattice and $\mathbb{P}$ be an ergodic measure on $\mathcal{C}$. Let $F_W$ be a family of measurable functions, indexed by bounded sets $W \in \mathbb{R}^d$, which are additive, i.e. $F_{W \cup W'} = F_W + F_{W'} - F_{W \cap W'}$, shift invariant on $\mathcal{L}$, i.e. $F_W(\gamma) = F_{W+i}(\gamma + i)$ for any $i \in \mathcal{L}$, and integrable, i.e. $\mathbb{E}[|F_{[0,1]^d}|] < \infty$. Then, we have that for $\mathbb{P}$-a.e. $\gamma$,

$$ (E.1) \qquad \lim_{n \to \infty} |W_n|^{-1} F_{W_n}(\gamma) = \mathbb{E}[|F_{[0,1]^d}|] $$

where $W_n$ is a sequence of regular bounded domains of $\mathbb{R}^d$ such that $W_n \to \mathbb{R}^d$ as $n \to \infty$.

PROOF. (i) Let $u \in \mathcal{D}$. Due to the decomposition of stationary measures as a mixture of ergodic measures, see [30], and by Theorem 1, it suffices to prove Theorem 2 by assuming that $\mathbb{P}_u$ is ergodic. We introduce the following random variables.

$$ A_{1,W_n} = \sum_{i \in \mathcal{L}_u} 1(i + X_i \in W_n \ominus R) f(i, X_i, \mathbf{\Gamma}_{ic}; \theta) $$
$$
A_{2,W_n} = \sum_{i \in \mathcal{L}_U} \mathbf{1}(i + X_i \in W_n \ominus R, X_i \in \underline{M}_n, i \in (W_n \ominus (m_n + R))^c) f(i, X_i, \mathbf{\Gamma}_{i^c}; \theta)
$$

$$
A_{3,W_n} = \sum_{i \in \mathcal{L}_U} \mathbf{1}(i + X_i \in W_n \ominus R, X_i \in \bar{M}_n, i \in (W_n \ominus (m_n + R))^c) f(i, X_i, \mathbf{\Gamma}_{i^c}; \theta)
$$

and

$$
\tilde{B}_{1,W_n} = \sum_{i \in \mathcal{L}_U} \int_Q \mathbf{1}(i \in W_n \ominus (m_n + R) f(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx
$$

$$
\begin{align*}
\tilde{B}_{2,W_n} &= \sum_{i \in \mathcal{L}_U} \int_Q \mathbf{1}(i \in W_n \ominus (m_n + R), x \in \underline{M}_n, i + x \in (W_n \ominus R)^c) \times \\
&\qquad f(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx
\end{align*}
$$

$$
\begin{align*}
\tilde{B}_{3,W_n} = \sum_{i \in \mathcal{L}_U} \int_Q \mathbf{1}(i \in W_n \ominus (m_n + R), x \in \overline{M}_n, i + x \in (W_n \ominus R)^c) \times \\
&\qquad f(i, x, \mathbf{\Gamma}_{i^c}; \theta) \Lambda(i, x, \mathbf{\Gamma}_{i^c}; \theta) dx
\end{align*}
$$

where for the sequence $(m_n)_n$ defined in (M1), we let $\overline{M}_n = \{x \in \mathcal{Q} : |x| \ge m_n\}$ and $\underline{M}_n = \mathcal{Q} \setminus \overline{M}_n$. We can observe, by shortening a little bit the notation, that

$$
A_{W_n} := A_{W_n}(f, \Gamma; \theta) = A_{1,W_n} - A_{2,W_n} - A_{3,W_n}
$$

$$
\tilde{B}_{W_n} := \tilde{B}_{W_n}(f, \Gamma; \theta) = \tilde{B}_{1, W_n} - \tilde{B}_{2, W_n} - \tilde{B}_{3, W_n}
$$

$$
\mathrm{DLR}_n := \mathrm{DLR}_n(f, \Gamma; \theta) = A_{W_n} - B_{W_n}
$$

$$
\widetilde{DLR}_n := \widetilde{DLR}_n(f, \Gamma; \theta) = A_{W_n} - \widetilde{B}_{W_n}
$$

$$
DLR_n - \widetilde{DLR}_n = B_{W_n} - \widetilde{B}_{W_n}.
$$

We point out the following facts: for any $i \in \mathcal{L}_U$ and $x \in \mathcal{Q}$,

$$
(E.2) \quad \mathbf{1}(i + X_i \in W_n \ominus R, X_i \in \underline{M}_n, i \in (W_n \ominus (m_n + R))^c) = 0
$$

$$
(E.3) \quad 1(i \in W_n \ominus (m_n + R), x \in \underline{M}_n, i + x \in (W_n \ominus R)^c) = 0
$$

$$
(E.4) \quad 1(i \in W_n \ominus (m_n + R), x \in \overline{M}_n, i + x \in (W_n \ominus R)^c) \le 1(i \in W_n)
$$

and that as $n \to \infty$

$$
(E.5) \quad |W_n| \sim |W_n \ominus R| \sim |W_n \ominus (m_n + R)| \text{ and } |W_n \oplus m_n \setminus W_n| \sim 2dm_n|W_n|^{1-1/d}.
$$

In particular, we can observe that $A_{2,W_n} = \tilde{B}_{2,W_n} = 0$.

The proof of Theorem 2(i) reduces to showing that, as $n \to \infty$ and given $U = u$, al-
most surely $|W_n|^{-1}A_{1,W_n}$ and $|W_n|^{-1}\tilde{B}_{1,W_n}$ tend to (4.13)-(4.14) respectively, and that
$|W_n|^{-1}C_{W_n} \to 0$ for $C_{W_n} = A_{3,W_n}, \tilde{B}_{3,W_n}$ and $B_{W_n} - \tilde{B}_{W_n}$.

Proof that $|W_n|^{-1}A_{1,W_n}$ and $|W_n|^{-1}\tilde{B}_{1,W_n}$ tend to (4.13)-(4.14) respectively. We first fo-
cus on the term $A_{1,W_n}$. Let $\tilde{\Gamma}$ be the marked point process defined by $\tilde{\Gamma} = \{(i+X_i, X_i), i \in$
$\mathcal{L}_U\}$. By assumption on $\Gamma$ and $f$, one can rewrite $A_{1,W_n}$ as

$$
A_{1,W_n} = \sum_{z=(y,x) \in \tilde{\Gamma} \cap ((W_n \ominus R)) \times \mathcal{Q})} f_1(x) f_2(y, \tilde{\Gamma} \setminus z) = \sum_{z \in \tilde{\Gamma} \cap (W_n \ominus R) \times \mathcal{Q}} \tilde{f}(z, \tilde{\Gamma})
$$
where we observe that $\tilde{f}$ is shift invariant with respect to the first coordinate of $z$. By assumption, $A_1$ is furthermore additive and such that $\mathbb{E}_u[|A_{1,[0,1]^d}|] < \infty$. We are in a position to apply Theorem 4, which leads to $\mathbb{P}_u$-a.s. as $n \to \infty$, $|W_n|^{-1}A_{1,W_n} \to \mathbb{E}_u[A_{1,[0,1]^d}]$. Finally, using Proposition 3.7(ii), observe that

$$ 
\begin{align*}
\mathbb{E}_u[A_{1,[0,1]^d}] &= \sum_{i \in \mathcal{L}_u} \int_Q \mathbb{1}_{[0,1]^d}(i+x) \mathbb{E}_u[f(i,x, \mathbf{\Gamma}_{ic}; \theta) \Lambda(i,x, \mathbf{\Gamma}_{ic}; \theta^*)] dx \\
&= \sum_{i \in \mathcal{L}_u} \int_Q \mathbb{1}_{[0,1]^d}(i+x) \mathbb{E}_u[f(0,x, \mathbf{\Gamma}_{0c}; \theta) \Lambda(0,x, \mathbf{\Gamma}_{0c}; \theta^*)] dx \\
&= \int_Q \mathbb{E}_u[f(0,x, \mathbf{\Gamma}_{0c}; \theta) \Lambda(0,x, \mathbf{\Gamma}_{0c}; \theta^*)] dx.
\end{align*}
$$

For the term, $\tilde{B}_{1,W_n}$, assumptions of Theorem 2 allow us to apply Theorem 4 to $\tilde{B}_{1,W_n}$. Hence, as $n \to \infty$ and $\mathbb{P}_u$-a.s.

$$ |W_n|^{-1} \tilde{B}_{1,W_n} \to \mathbb{E}_u \left[ \int_Q f(0, x, \mathbf{\Gamma}_{0c}; \theta) \Lambda(0, x, \mathbf{\Gamma}_{0c}; \theta) dx \right]. $$

*Proof that* $|W_n|^{-1}A_{3,W_n} \xrightarrow{a.s.} 0$. Using the stationarity of $\mathbf{\Gamma}$ | $U = u$, (A)[f] and (E.4), there exists $c > 0$ such that

$$ 
\begin{align*}
|W_n|^{-1} \mathbb{E}_u |A_{3,W_n}| &\le c |W_n|^{-1} \#(\mathcal{L}_u \cap W_n) \int_{\overline{M}_n} \lambda_1(x; \theta^*) dx \\
&\le c \mathbb{P}(X \ge m_n).
\end{align*}
$$

for $X \sim \mathbb{Q}$. Using (M1) and Markov's inequality, it is clear that $\mathbb{P}(X \ge m_n) = O(e^{-m_n}) = O(|W_n|^{-\beta})$. Finally, Markov's inequality, (M1) and Borel-Cantelli's lemma conclude this step.

*Proof that* $|W_n|^{-1} \tilde{B}_{3,W_n} \xrightarrow{a.s.} 0$. This result is proved along similar lines as the previous result, except we now combine (B)[f], (M1) and (E.3)-(E.5).

*Proof that* $|W_n|^{-1}(B_{W_n} - \tilde{B}_{W_n}) \xrightarrow{a.s.} 0$. We first observe that

$$ (E.6) \qquad |B_{W_n} - \tilde{B}_{W_n}| \le \int_Q \int_Q Z_n(i, x, y, \mathbf{\Gamma}_{ic}; \theta) (1 - b_n(i, y)) dx dy $$

with

$$ Z_n(i, x, y, \mathbf{\Gamma}_{ic}; \theta) = b_n(i, x) |f(i, x, \mathbf{\Gamma}_{ic}; \theta)| \Lambda(i, x, \mathbf{\Gamma}_{ic}; \theta) \Lambda_n(i, y, \mathbf{\Gamma}_{ic}; \theta) $$

Now, we decompose the right-hand side of (E.6) as $\Delta_{\underline{M}_n} + \Delta_{\overline{M}_n}$, where

$$ \Delta_{\underline{M}_n} = \int_Q \int_Q Z_n(i, x, y, \mathbf{\Gamma}_{ic}; \theta) \times \mathbb{1}(i \in W_n \ominus (m_n + R), y \in \underline{M}_n, i + y \in (W_n \ominus R)^c) dx dy. $$

$$ \Delta_{\overline{M}_n} = \int_Q \int_Q Z_n(i, x, y, \mathbf{\Gamma}_{ic}; \theta) \times \mathbb{1}(i \in W_n \ominus (m_n + R), y \in \overline{M}_n, i + y \in (W_n \ominus R)^c) dx dy. $$
Using (E.3), $\Delta_{\underline{M}_n} = 0$. We leave the reader to check that using (B)[f] and the same arguments as before, we have that

$$|W_n|^{-1} \mathbb{E}_u |\Delta_{\overline{M}_n}| = O(|W_n|^{-\beta})$$

whereby the result is deduced using Markov's inequality and Borel-Cantelli's lemma applied to $|W_n|^{-1} \Delta_{\overline{M}_n}$.

(ii) Finally, to prove (4.17), we observe that the previous developments also prove that

$$\Delta_{\underline{M}_n} = O_{\mathbb{P}_u}(\mathbb{E}_u |\Delta_{\overline{M}_n}|) = O_{\mathbb{P}_u}(|W_n|^{1-\beta}).$$

Therefore, $B_{W_n} - \tilde{B}_{W_n} = O_{\mathbb{P}_u}(|W_n|^{1-\beta})$ which yields the result. $\square$

# APPENDIX F: PROOF OF THEOREM 3

## F.1. Auxiliary result.

LEMMA F.1. Assume (M1)-(M3). For $\ell = 0, 1, 2$, $l = 1, \dots, L$ and $\theta \in \Theta$, we have $\mathbb{P}_u$-a.s. as $n \to \infty$

$$(F.1) \qquad |W_n|^{-1} \mathrm{DLR}_n^{(\ell)}(f_l, \Gamma; \theta) \xrightarrow{a.s.} \mathcal{D}^{(\ell)}(f_l; \theta)$$

where

$$
\begin{align*}
\mathcal{D}^{(0)}(f_l; \theta) &= \mathcal{D}(f_l; \theta) = \mathbb{E}_u \left[ \int_Q f_l(0, x, \mathbf{\Gamma}_{0^c}; \theta) (\Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta^*) - \Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta)) \mathrm{d}x \right] \\
\mathcal{D}^{(1)}(f_l; \theta) &= \mathbb{E}_u \left[ \int_Q f_l^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta) (\Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta^*) - \Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta)) \mathrm{d}x \right] \\
&\quad - \mathbb{E}_u \left[ \int_Q f_l(0, x, \mathbf{\Gamma}_{0^c}; \theta) \Lambda^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta) \mathrm{d}x \right] \\
\mathcal{D}^{(2)}(f_l; \theta) &= \mathbb{E}_u \left[ \int_Q f_l^{(2)}(0, x, \mathbf{\Gamma}_{0^c}; \theta) (\Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta^*) - \Lambda(0, x, \mathbf{\Gamma}_{0^c}; \theta)) \mathrm{d}x \right] \\
&\quad - \mathbb{E}_u \left[ \int_Q f_l^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta) \Lambda^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta)^\top \mathrm{d}x \right] \\
&\quad - \mathbb{E}_u \left[ \int_Q \Lambda^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta) f_l^{(1)}(0, x, \mathbf{\Gamma}_{0^c}; \theta)^\top \mathrm{d}x \right] \\
&\quad - \mathbb{E}_u \left[ \int_Q f_l(0, x, \mathbf{\Gamma}_{0^c}; \theta) \Lambda^{(2)}(0, x, \mathbf{\Gamma}_{0^c}; \theta)^\top \mathrm{d}x \right].
\end{align*}
$$

In particular, we have $\mathcal{D}(f_l; \theta^*) = 0$ and $\mathcal{D}^{(1)}(f_l; \theta^*) = -\mathcal{I}(f_l; \theta^*)$, where $\mathcal{I}(f_l; \theta^*)$ is given by (4.11).

PROOF. The proof of this result is a direct application of Theorem 2, valid under the assumptions (M1)-(M3). $\square$

## F.2. Proof of Theorem 3.
PROOF. (i) We have, by Proposition 3.7, that

$$
\begin{align*}
\mathrm{Var}_u[\widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}; \theta^*)] &= \mathbb{E}_u[\widetilde{\mathrm{DLR}}_n(f_l; \theta^*)^2] \\
&= \mathbb{E}\left[\left\{\sum_{i \in \mathcal{L}_U} b_n(i, X_i) f_l(i, X_i, \mathbf{\Gamma}_{i^c}) \right.\right. \\
&\qquad \left. \left. - \int b_n(i, x) f_l(i, x, \mathbf{\Gamma}_{i^c}) \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx \right\}^2\right]
\end{align*}
$$

where to abuse and ease notation we let $f_l(i, x, \gamma) = f_l(i, x, \gamma; \theta^*)$, $\Lambda(i, x, \gamma) = \Lambda(i, x, \gamma; \theta^*)$
and $\int = \int_Q$. Using Propositions 3.7 and 3.8, we decompose the conditional variance

$$
\mathrm{Var}_u[\widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}; \theta^*)] = T_1 + T_2
$$

where

$$
T_1 = \mathbb{E}_u \left[ \sum_{i \in \mathcal{L}_U} \int b_n(i, x) f_l(i, x, \mathbf{\Gamma}_{i^c})^2 \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx \right. \\
\left. - \left\{ \int b_n(i, x) f_l(i, x, \mathbf{\Gamma}_{i^c}) \Lambda(i, x, \mathbf{\Gamma}_{i^c}) dx \right\}^2 \right]
$$

$$
T_2 = \mathbb{E}_u \left[ \sum_{\substack{i,j \in \mathcal{L}_U \\ \neq}} \iint b_n(i,x)b_n(j,y)f_l(i,x,\mathbf{\Gamma}_{i^c})f_l(j,y,\mathbf{\Gamma}_{j^c})\Lambda(i,x,\mathbf{\Gamma}_{i^c})\Lambda(j,y,\mathbf{\Gamma}_{j^c})\Delta x dy \right]
$$

and

$$
\Delta = \Delta(i, x, j, y, \Gamma_{i^c}, \Gamma_{j^c}, \Gamma_{(ij)^c}) = \frac{\Lambda(i, x, \Gamma_{(ij)^c})\Lambda(j, y, \Gamma_{(ij)^c} \cup \{i+x\})}{\Lambda(i, x, \Gamma_{i^c})\Lambda(j, y, \Gamma_{j^c})} - 1.
$$

Regarding the term $T_1$, using the stationarity of $\mathbf{\Gamma} \mid U = u$, Fubini's theorem and by Assumption (M3), we have the simple upper-bound

$$
(F.2) \quad T_1 \le \int \#(\mathcal{L}_u \cap (W_n \ominus (m_n + R))) \mathbb{E}_u [|f_l(0, x, \mathbf{\Gamma}_{0^c})| \Lambda(0, x, \mathbf{\Gamma}_{0^c})] dx = O(|W_n|).
$$

Let us consider the term $T_2$. Let $a, b \in \mathbb{R}^d$ and $\gamma \in \mathcal{C}$. By definition of $\lambda_2$, we have

$$
\lambda_2(a \cup b, \gamma) = \lambda_2(a, \gamma)\lambda_2(b, \gamma \cup a) = \lambda_2(b, \gamma)\lambda_2(a, \gamma \cup b)
$$

which can be rewritten as

$$
(F.3) \qquad \frac{\lambda_2(a, \gamma)}{\lambda_2(a, \gamma \cup b)} = \frac{\lambda_2(b, \gamma)}{\lambda_2(b, \gamma \cup a)}.
$$

Using this, we rewrite $\Delta + 1$ as

$$
\Delta + 1 = A(i, x, X_i, j, y, X_j, \mathbf{\Gamma}_{i^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(ij)^c}) \times \\
\frac{\iint \lambda(i, \tilde{x}, \mathbf{\Gamma}_{i^c}) \lambda(j, \tilde{y}, \mathbf{\Gamma}_{j^c}) d\tilde{x} d\tilde{y}}{\iint \lambda(i, \tilde{x}, \mathbf{\Gamma}_{i^c}) \lambda(j, \tilde{y}, \mathbf{\Gamma}_{j^c}) A(i, \tilde{x}, X_i, j, \tilde{y}, X_j, \mathbf{\Gamma}_{i^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(ij)^c}) d\tilde{x} d\tilde{y}}
$$

where for any $i, j \in \mathcal{L}_U$, $x, x', y, y' \in \mathcal{Q}$, we have, using the definition of $\Lambda$

$$
A(i, x, X_i, j, y, X_j, \mathbf{\Gamma}_{i^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(ij)^c}) = \frac{\lambda_2(i + x, \mathbf{\Gamma}_{(ij)^c}) \lambda_2(j + y, \mathbf{\Gamma}_{(ij)^c} \cup \{i + x\})}{\lambda_2(i + x, \mathbf{\Gamma}_{i^c}) \lambda_2(j + y, \mathbf{\Gamma}_{j^c})}
$$
Using (F.3), this term can be rewritten as

$$
\begin{equation}
\begin{aligned}
A(i, x, X_i, j, y, X_j, \mathbf{\Gamma}_{ic}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(ij)c}^c) = & \frac{\lambda_2(j + X_j, \mathbf{\Gamma}_{(ij)c}^c)}{\lambda_2(j + X_j, \mathbf{\Gamma}_{(ij)c}^c) \cup \{i+x\})} \\
& \times \frac{\lambda_2(i + X_i, \mathbf{\Gamma}_{(ij)c}^c \cup \{i+x\})}{\lambda_2(i + X_i, \mathbf{\Gamma}_{(ij)c}^c \cup \{i+x\} \cup \{j+y\})} \\
& \times \frac{\lambda_2(j + y, \mathbf{\Gamma}_{je}^c \cup \{i+x\})}{\lambda_2(j + y, \mathbf{\Gamma}_{je}^c)}.
\end{aligned}
\tag{F.4}
\end{equation}
$$

Hence, we view the term $T_2$ as a general functional of $i, X_i, j, X_j$ and $\Gamma$, to which we can apply Proposition 3.8 once more. We obtain that

$$
T_2 = \mathbb{E}_u \sum_{i,j \in \mathcal{L}_U} \int \int \int \int b_n(i,x)b_n(j,y) F_l(i,x,x',j,y,y', \mathbf{\Gamma}_{ic}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(ij)c}^c) \times \\
\Delta(i,x,x',j,y,y', \mathbf{\Gamma}_{ic}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(ij)c}^c) \lambda_1(x) \lambda_1(x') \lambda_1(y) \lambda_1(y') dx dx' dy dy'
$$

where

$$
\begin{align*}
F_l(i, x, x', j, y, y', \mathbf{\Gamma}_{ic}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(ij)c}^c) ={}& f_l(i, x, \mathbf{\Gamma}_{ic}^c) \lambda_2(i + x, \mathbf{\Gamma}_{ic}^c) \\
& \times f_l(j, y, \mathbf{\Gamma}_{je}^c) \lambda_2(j + y, \mathbf{\Gamma}_{je}^c) \\
& \times \lambda_2(i + x', \mathbf{\Gamma}_{(ij)c}^c) \lambda_2(j + y', \mathbf{\Gamma}_{(ij)c}^c \cup \{i + x'\}).
\end{align*}
$$

Similarly to the term $T_1$, we use Fubini's theorem and the stationarity of $\Gamma | U = u$ to derive

$$
\begin{align*}
T_2 ={}& \#(\mathcal{L}_u \cap (W_n \ominus (m_n + R))) \times \\
& \mathbb{E}_u \sum_{j \in \mathcal{L}_U \setminus \{0\}} \int \int \int \int 1(x \in W_n \ominus R) b_n(j, y) F_l(0, x, x', j, y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \times \\
& \qquad \Delta(0, x, x', j, y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \lambda_1(x) \lambda_1(x') \lambda_1(y) \lambda_1(y') dx dx' dy dy' \\
& = \#(\mathcal{L}_u \cap W_n \ominus (m_n + R)) \times (T'_1 + T'_2)
\end{align*}
$$

where

$$
T'_1 = \mathbb{E}_u \sum_{j \in \mathcal{L}_U} \int \int \int \int 1(x \in W_n \ominus R) b_n(j, y) \\
F_l(0, x, x', y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \times R_1(0, x, x', j, y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \times \\
\{A(0, x, x', j, y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) - 1\} \, dx \, dx' \, dy \, dy'
$$

$$
T'_2 = \mathbb{E}_u \sum_{j \in \mathcal{L}_U} \int \int \int \int \int 1(x \in W_n \ominus R) b_n(j, y) \\
F_l(0, x, x', y, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \times R_2(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) \times \\
\{A(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0c}^c, \mathbf{\Gamma}_{je}^c, \mathbf{\Gamma}_{(0j)c}^c) - 1\} \lambda_1(\tilde{x}) \lambda_1(\tilde{y}) dx dx' d\tilde{x} dy d\tilde{y}
$$

where

$$
R_1 = \frac{\iint \lambda(0, \tilde{x}, \mathbf{\Gamma}_{0c}) \lambda(j, \tilde{y}, \mathbf{\Gamma}_{je}) d\tilde{x} d\tilde{y}}{\iint \lambda(0, \tilde{x}, \mathbf{\Gamma}_{0c}) \lambda(j, \tilde{y}, \mathbf{\Gamma}_{je}) A(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0c}, \mathbf{\Gamma}_{je}, \mathbf{\Gamma}_{(0j)c}) d\tilde{x} d\tilde{y}}
$$
and

$$R_2 = \frac{\lambda_2(\tilde{x}, \mathbf{\Gamma}_{0^c})\lambda_2(j + \tilde{y}, \mathbf{\Gamma}_{j^c})}{\iint \lambda(0, \tilde{x}, \mathbf{\Gamma}_{0^c})\lambda(j, \tilde{y}, \mathbf{\Gamma}_{j^c})A(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c})d\tilde{x}d\tilde{y}}$$

We are now in a position to finally use the finite range property (H5.1). Indeed, for any $x, x', y, y' \in \mathcal{Q}$, any $j \in \mathcal{L}_U$, we can observe that $A(0, x, x', j, y, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) = 1$ if $j \notin \mathcal{B}_R(x, x', y, y')$ defined by

$$\mathcal{B}_R(x, x', y, y') = B(x - y', R) \cup B(x' - y, R) \cup B(x - y, R).$$

Now, let $\mathcal{E}_u$ be given by

$$\mathcal{E}_u = \sup_{l=1,...,L} \sup_{x,x',\tilde{x},y,y',\tilde{y} \in \mathcal{Q}} \sup_{j \in \mathcal{L}_u} \mathbb{E}_u \left[ \begin{array}{l} \left| F_l(0, x, x', y, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) \right| \\ \qquad R_1(0, x, x', j, y, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) \left| A(0, x, x', j, y, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) - 1 \right| + \\ \qquad (F.5) \quad R_2(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) \times \left| A(0, \tilde{x}, x', j, \tilde{y}, y', \mathbf{\Gamma}_{0^c}, \mathbf{\Gamma}_{j^c}, \mathbf{\Gamma}_{(0j)^c}) - 1 \right| \end{array} \right].$$

Now, we have on the one hand

$$\begin{aligned} T'_1 &\le \mathcal{E}_u \iint \int \sum_{j \in \mathcal{B}_R(x, x', y, y')} \lambda_1(x)\lambda_1(x')\lambda_1(y)\lambda_1(y')dx dx' dy dy' \\ &\le 3\mathcal{E}_u \iint \left\{ \#(\mathcal{L}_u \cap B(x-y, R)) \right\} \lambda_1(x)\lambda_1(y)dx dy \\ &\le 3\mathcal{E}_u \left\{ \sup_{x,y \in \mathcal{Q}} \#(\mathcal{L}_u \cap B(x-y, R)) \right\} \iint \lambda_1(x)\lambda_1(y)dx dy = O(1). \end{aligned}$$

On the other hand,

$$\begin{aligned} T'_2 &\le \mathcal{E}_u \iint \int \sum_{j \in \mathcal{B}_R(\tilde{x}, x', \tilde{y}, y') \cap (W_n \ominus (m_n + R))_{-y}} \lambda_1(x)\lambda_1(x')\lambda_1(\tilde{x})\lambda_1(y)\lambda_1(y')\lambda_1(\tilde{y})dx dx' d\tilde{x} dy dy' d\tilde{y} \\ &\le 3\mathcal{E}_u \left\{ \sup_{x,y \in \mathcal{Q}} \#(\mathcal{L}_u \cap B(x-y, R)) \right\} \iint \lambda_1(x)\lambda_1(y)dx dy = O(1) \end{aligned}$$

whereby we deduce, recalling (F.2), that $T_1 + T_2 = \#(\mathcal{L}_u \cap W_n \ominus (m_n + R)) \times (O(1) + T'_1 + T'_2) = O(|W_n|)$, which finally leads to the result.

(ii) We now focus on the consistency of $\hat{\theta}$. Let $k \in \mathbb{R}^p$, let $u_n = cv_n^{1/2}|W_n|^{-1}$ for some $c > 0$, we can observe that $u_n \to 0$ as $n \to \infty$. Let also $\Delta_n(k) = \text{TF}_n(\mathbf{f}, \mathbf{\Gamma}; \theta^*) - \text{TF}_n(\mathbf{f}, \mathbf{\Gamma}; \theta^* + u_n k)$. Since $\hat{\theta} = \text{argmax}_{\theta}(-\text{TF}_n(\mathbf{f}, \mathbf{\Gamma}; \theta))$, the result is proved if for any $\varepsilon > 0$, there exists $K = K(U)$ such that

$$P\left(\sup_{|k|=K} \Delta_n(k) > 0 \mid U = u\right) \le \varepsilon.$$

We shorten and introduce a few notation.

$$e_n(\theta) = \text{TF}_n^{(1)}(\mathbf{f}, \gamma; \theta) = 2 \sum_{l=1}^{L} \text{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}, \theta) \text{DLR}_n(f_l, \mathbf{\Gamma}, \theta)$$
$$ \tilde{e}_n(\theta) = 2 \sum_{l=1}^{L} \mathrm{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}, \theta) \widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}, \theta) $$

$$ \tilde{e}_n(\theta) = 2|W_n| \sum_{l=1}^{L} \mathcal{I}_l(\theta^*) \widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}, \theta) $$

Under assumption (M2), there exists $t \in (0, 1)$ such that, using a Taylor expansion, we can decompose $\Delta_n(k)$ as

$$ \begin{aligned} \Delta_n(k) &= -u_n k^\top \mathrm{TF}_n^{(1)}(\mathbf{f}, \mathbf{\Gamma}, \theta^* + t u_n k) \\ &= -u_n k^\top e_n(\theta^* + t u_n k) \\ &= T_1 + T_2 + T_3 + T_4 \end{aligned} $$

where, using the short notation $\theta_n^* = \theta^* + t u_n k$,

$$ T_1 = u_n k^\top (\tilde{e}_n(\theta_n^*) - e_n(\theta_n^*)) $$

$$ T_2 = u_n k^\top (\tilde{e}_n(\theta^*) - \tilde{e}_n(\theta_n^*)) $$

$$ T_3 = u_n k^\top (\tilde{e}_n(\theta^*) - \tilde{e}_n(\theta_n^*)) $$

$$ T_4 = -u_n k^\top \tilde{e}_n(\theta^*) $$

We now treat these four terms separately. In the following, we abuse notation by denoting $o_{a.s.}$ the almost sure behaviour with respect to $\mathbb{P}_u$.

Term $T_1$. We have that

$$ \tilde{e}_n(\theta_n^*) - e_n(\theta_n^*) = 2 \sum_l \mathrm{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta_n^*) \left( \widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}; \theta_n^*) - \mathrm{DLR}_n(f_l, \mathbf{\Gamma}; \theta_n^*) \right) $$

Let $\theta \in \Theta$. On the other hand, we can apply Lemma F.1 and obtain

$$ \mathrm{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta) = O_{\mathrm{a.s.}}(|W_n|). $$

We may also apply, under assumptions (M1)-(M3), Theorem 2(ii) to $f = f_l$ to derive

$$ \widetilde{\mathrm{DLR}}_n(f_l, \mathbf{\Gamma}; \theta_n^*) - \mathrm{DLR}_n(f_l, \mathbf{\Gamma}; \theta_n^*) = O_{\mathbb{P}_u}(|W_n|^{1-\beta}). $$

Therefore, for any $\varepsilon > 0$, we have, as $n \to \infty$

$$ |W_n|^{-2+\beta-\varepsilon} (\tilde{e}_n(\theta_n^*) - e_n(\theta_n^*)) \xrightarrow{\mathbb{P}_u} 0 $$

whereby we deduce by the continuous mapping theorem that

$$ \begin{aligned} u_n (\tilde{e}_n(\theta_n^*) - e_n(\theta_n^*)) &= O_{\mathbb{P}_u} (u_n |W_n|^{2-\beta+\varepsilon}) \\ &= O_{\mathbb{P}_u} \left( v_n \frac{|W_n|^{1-\beta+\varepsilon}}{\sqrt{v_n}} \right) \\ &= o_{\mathbb{P}_u} (v_n) \end{aligned} $$

for some small $\varepsilon > 0$ and as soon as $\beta > 1/2$. Hence, $T_1 \le |k|T'_1$ with $T'_1 = o_{\mathbb{P}_u}(v_n)$.
Term $T_2$. Under assumption (M2), there exists $s \in (0, 1)$ such that

$$T_2 = -tu_n^2 k^\top \tilde{e}_n^{(1)}(\tilde{\theta}_n^\star)k$$

where $\tilde{\theta}_n^\star = \theta^\star + \text{stku}_n \to \theta^\star$ as $n \to \infty$, and where for any $\theta \in \Theta$

$$\tilde{e}_n^{(1)}(\theta) = 2 \sum_{l=1}^{L} \left\{ \begin{aligned} & \text{DLR}_n^{(2)}(f_l, \mathbf{\Gamma}; \theta) \widetilde{\text{DLR}}_n(f_l, \mathbf{\Gamma}; \theta) + \\ & \text{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta) \widetilde{\text{DLR}}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta)^\top \end{aligned} \right\}.$$

Lemma F.1 and Theorem 2 show that for any $\theta \in \Theta$, $\mathbb{P}_u$-a.s. as $n \to \infty$

$$|W_n|^{-2} \tilde{e}_n^{(1)}(\theta) \xrightarrow{a.s.} 2 \sum_{l=1}^{L} \left\{ \mathcal{D}^{(2)}(f_l; \theta) \mathcal{D}(f_l; \theta) + \mathcal{D}^{(1)}(f_l; \theta) \mathcal{D}^{(1)}(f_l; \theta)^\top \right\}.$$

By the continuous mapping theorem, we deduce that $\mathbb{P}_u$-a.s. as $n \to \infty$

$$|W_n|^{-2} \tilde{e}_n^{(1)}(\check{\theta}_n^\star) \xrightarrow{a.s.} 2 \sum_{l=1}^{L} \left\{ \begin{aligned} & \mathcal{D}^{(2)}(f_l; \theta^\star) \mathcal{D}(f_l; \theta^\star) + \mathcal{D}^{(1)}(f_l; \theta^\star) \mathcal{D}^{(1)}(f_l; \theta^\star)^\top \end{aligned} \right\} \\ =: 2 \sum_{l=1}^{L} \mathcal{I}_l(\theta^\star) \mathcal{I}_l(\theta^\star)^\top = 2\mathcal{I}(\theta^\star).$$

We finally obtain that

$$\begin{aligned} T_2 &= -tu_n^2 |W_n|^2 |k|^2 \frac{k^\top(2\mathcal{I}(\theta^\star))k}{|k|^2} + \\ &tu_n^2 |W_n|^2 k^\top \{2\mathcal{I}(\theta^\star) - |W_n|^{-2}\tilde{e}_n^{(1)}(\check{\theta}_n^\star)\}k \\ &\leq -\kappa |k|^2 u_n^2 |W_n|^2 + |k| T_2' \end{aligned}$$

where $\kappa = 2t\lambda_{\min}(\mathcal{I}(\theta^\star)) > 0$ and

$$\begin{aligned} T_2' &= \sup_{k, \|k\|=K} |k| u_n^2 |W_n|^2 \left| 2\mathcal{I}(\theta^\star) - |W_n|^{-2} \tilde{e}_n^{(1)}(\check{\theta}_n^\star) \right| \\ &= o_{\text{a.s.}}(u_n^2 |W_n|^2) = o_{\text{a.s.}}(v_n). \end{aligned}$$

Term $T_3$. We can observe that

$$\tilde{e}_n(\theta^\star) - \check{e}_n(\theta^\star) = 2|W_n| \sum_{l=1}^{L} \left\{ |W_n|^{-1} \text{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta^\star) - \mathcal{I}_l(\theta^\star) \right\} \widetilde{\text{DLR}}_n(f_l, \mathbf{\Gamma}; \theta^\star).$$

Under assumptions (M1), (M2) and (M3), one can apply Lemma F.1 and Theorem 3(i) to show that

$$|W_n|^{-1} \text{DLR}_n^{(1)}(f_l, \mathbf{\Gamma}; \theta^\star) - \mathcal{I}_l(\theta^\star) = o_{\text{a.s.}}(1) \quad \text{and} \quad \widetilde{\text{DLR}}_n(f_l, \mathbf{\Gamma}, \theta^\star) = O_{\mathbb{P}_u}(v_n^{1/2}).$$

Since $u_n|W_n| = O(v_n^{1/2})$, we deduce that $T_3 \le |k|T_3'$ with $T_3' = o_{\mathbb{P}_u}(v_n)$.

Term $T_4$. By definition and Proposition 3.7, the random variable $\tilde{e}_n(\theta^\star)$ is centered. Therefore, we have in particular from Theorem 3(i) and conditional Markov's inequality, that $\tilde{e}_n(\theta^\star) = O_{\mathbb{P}_u}(\sqrt{\text{Var}_u[\tilde{e}_n(\theta^\star)]}) = O_{\mathbb{P}_u}(v_n^{1/2})$
**Conclusion.** We have that

$$ \Delta_n(k) \le u_n |W_n| |k| |\tilde{e}_n(\theta^*)| - \kappa |k|^2 u_n^2 |W_n|^2 + |k|T $$

with $T = T'_1 + T'_2 + T'_3$. Let $u \in D$, we then have

$$ P\left( \sup_{|k|=K} \Delta_n(k) > 0 \mid U = u \right) \le P\left( |\tilde{e}_n(\theta^*)| + u_n^{-1} |W_n|^{-1} T \ge u_n |W_n| \kappa K \mid U = u \right). $$

Now, the choice of the sequence $u_n$ leads to $u_n|W_n| = cv_n^{1/2}$ and to

$$ |\tilde{e}_n(\theta^*)| + u_n^{-1} |W_n|^{-1} T = O_{\mathbb{P}_u}(v_n^{1/2}) + v_n^{-1/2} o_{\mathbb{P}_u}(v_n) = O_{\mathbb{P}_u}(v_n^{1/2}) $$

whereby we deduce the result. □

**Acknowledgments.** The research of C. Renaud Chan and J.-F. Coeurjolly is supported by Persyval-lab (ANR-11-61 LABX-0025-01). The authors would like to thank David Dereudre for insightful discussions and Rémy Drouilhet for assistance with the implementation of the simulation algorithm.

## REFERENCES

[1] ALMEIDA, M. P. and GIDAS, B. (1993). A variational method for estimating the parameters of MRF from complete or incomplete data. *The Annals of Applied Probability* 103–136.

[2] BADDELEY, A. and DEREUDRE, D. (2013). Variational estimators for the parameters of Gibbs point process models. *Bernoulli* 19 905–930.

[3] BADDELEY, A. and TURNER, R. (2000). Practical Maximum Pseudolikelihood for Spatial Point Patterns: (with Discussion). *Australian & New Zealand Journal of Statistics* 42 283–322.

[4] BESAG, J. (1974). Spatial interaction and the statistical analysis of lattice systems. *Journal of the Royal Statistical Society: Series B (Methodological)* 36 192–225.

[5] BUISINE, R., VANBESIEN, O., DEREUDRE, D., BURGNIES, L. and LHEURETTE, É. (2025). Absorption properties of hyperuniform and Poisson Point Process distributed metasurfaces at infrared wavelengths. *Journal of Applied Physics* 138.

[6] COEURJOLLY, J.-F., DEREUDRE, D., DROUILHET, R. and LAVANCIER, F. (2012). Takacs-Fiksel method for stationary marked Gibbs point processes. *Scandinavian Journal of Statistics* 39 416–443.

[7] COEURJOLLY, J.-F. and DROUILHET, R. (2010). Asymptotic properties of the maximum pseudo-likelihood estimator for stationary Gibbs point processes including the Lennard-Jones model. *Electronic Journal of Statistics* 4 677–706.

[8] COEURJOLLY, J.-F. and LAVANCIER, F. (2017). Parametric estimation of pairwise Gibbs point processes with infinite range interaction. *Bernoulli* 23 1299–1334.

[9] DEREUDRE, D. (2019). Introduction to the theory of Gibbs point processes. In *Stochastic Geometry: Modern Research Frontiers* 181–229. Springer.

[10] DEREUDRE, D., DROUILHET, R. and GEORGII, H.-O. (2012). Existence of Gibbsian point processes with geometry-dependent interactions. *Probability Theory and Related Fields* 153 643–670.

[11] DEREUDRE, D. and FLIMMEL, D. (2024). Non-hyperuniformity of Gibbs point processes with short-range interactions. *Journal of Applied Probability* 61 1380–1406.

[12] DEREUDRE, D., FLIMMEL, D., HUESMANN, M. and LEBLÉ, T. (2024). (Non)-hyperuniformity of perturbed lattices. *arXiv:2405.19881*.

[13] DEREUDRE, D. and LAVANCIER, F. (2011). Practical simulation and estimation for Gibbs Delaunay–Voronoi tessellations with geometric hardcore interaction. *Computational Statistics & Data Analysis* 55 498–519.

[14] DEREUDRE, D., LAVANCIER, F. and HELISOVÁ, K. S. (2014). Estimation of the Intensity Parameter of the Germ-Grain Quermass-Interaction Model when the Number of Germs is not Observed. *Scandinavian Journal of Statistics* 41 809–829.

[15] FIKSEL, T. (1984). Estimation of parameterized pair potentials of marked and non-marked Gibbsian point processes. *Elektronische Informationsverarbeitung und Kybernetik* 20 270–278.

[16] FIKSEL, T. (1988). Estimation of interaction potentials of Gibbsian point processes. *Statistics* 19 77–86.

[17] FRIEDLI, S. and VELENIK, Y. (2017). *Statistical mechanics of lattice systems: a concrete mathematical introduction*. Cambridge University Press.
[18] GAETAN, C. and GUYON, X. (2010). *Spatial statistics and modeling* **90**. Springer.

[19] GEORGII, H.-O. (1979). *Canonical Gibbs Measures, 760 of Lecture Notes in Mathematics*. Springer, Berlin, Heidelberg, New York.

[20] GEORGII, H.-O. (2011). *Gibbs Measures and Phase Transitions*. De Gruyter, Berlin, New York.

[21] HOUGH, J. B., KRISHNAPUR, M., PERES, Y. et al. (2009). *Zeros of Gaussian analytic functions and determinantal point processes* **51**. American Mathematical Society.

[22] ILLIAN, J., PENTTINEN, A., STOYAN, H. and STOYAN, D. (2008). *Statistical analysis and modelling of spatial point patterns*. John Wiley & Sons.

[23] JENSEN, J. L. and MØLLER, J. (1991). Pseudolikelihood for exponential family models of spatial point processes. *The Annals of Applied Probability* **1** 445–461.

[24] KENDALL, W. S., VAN LIESHOUT, M. N. M. and BADDELEY, A. J. (1999). Quermass-interaction processes: conditions for stability. *Advances in Applied Probability* **31** 315–342.

[25] KIDERLEN, M. and HÖRIG, M. (2013). Matérn's hard core models of types I and II with arbitrary compact grains CSGB Research Reports No. 05, Centre for Stochastic Geometry and Advanced Bioimaging, Aarhus University.

[26] LACHIÈZE-REY, R. (2025). Hyperuniform random measures, transport and rigidity. *arxiv:2510.18392*.

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

[28] MØLLER, J. and WAAGEPETERSEN, R. P. (2003). *Statistical inference and simulation for spatial point processes*. CRC press.

[29] NGUYEN, X. X. and ZESSIN, H. (1977). Integral and differential characterizations of the Gibbs process. *Advances in Applied Probability* **9** 446–447.

[30] PRESTON, C. (1976). *Random fields* **534**. Springer.

[31] RUELLE, D. (1970). Superstable interactions in classical statistical mechanics. *Communications in Mathematical Physics* **18** 127–159.

[32] STRAUSS, D. J. (1975). A model for clustering. *Biometrika* **62** 467–475.

[33] TAKACS, R. (1986). Estimator for the pair-potential of a Gibbsian point process. *Statistics: A Journal of Theoretical and Applied Statistics* **17** 429–433.

[34] TORQUATO, S. (2016). Hyperuniformity and its generalizations. *Physical Review E* **94** 022122.

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

[36] TORQUATO, S. and STILLINGER, F. H. (2003). Local density fluctuations, hyperuniformity, and order metrics. *Physical Review E* **68** 041113.

[37] VALKÓ, B. and VIRÁG, B. (2009). Continuum limits of random matrices and the Brownian carouse. *Inventiones mathematicae* **177** 463–508.