D Dereudre
email: david.dereudre@univ-valenciennes.fr

F Lavancier
email: frederic.lavancier@univ-nantes.fr

Practical simulation and estimation for Gibbs Delaunay-Voronoi tessellations with geometric hardcore interaction

Keywords: Gibbs point process, random tessellations, stochastic geometry, pseudo-likelihood estimator, spatial statistics

HAL is

Introduction

In the domains of physics and biology, some large-scale random geometric structures can be mathematically modeled using Poisson-Voronoi or Poisson-Delaunay tessellations. In cosmology for instance, since [START_REF] Matsuda | Topology of supercluster-void structure[END_REF], modeling the large-scale galaxy distribution generally relies on Voronoi tessellations (see [START_REF] Icke | Fragmenting the universe[END_REF] and [START_REF] Icke | Fragmenting the universe. II -Voronoi vertices as Abell clusters[END_REF]). In biology, Voronoi tessellations are often used to model the cellular configuration of a tissue (since the seminal work of [START_REF] Honda | Description of cellular patterns by Dirichlet domains: the two-dimensional case[END_REF]). This tool is also relevant to model the geometrical structure of proteins (cf. [START_REF] Poupon | Voronoi and Voronoi-related tessellations in studies of protein structure and interaction[END_REF] for a state of the art) or microstructures like foams. Mathematical properties of Poisson-Voronoi and Poisson-Delaunay tessellations have been widely studied (see [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF] for instance).

Unfortunately, these models have the disadvantage of yielding strong independence properties due to the Poissonian nature of the underlying point process. In different biological studies, the necessity to introduce an interaction between the cells of the tessellation to achieve greater realism has indeed been emphasized. In [START_REF] Farhadifar | The influence of cell mechanics, Cell-Cell interactions, and Proliferation on Epithelial Packing[END_REF] for instance, the interaction between neighbouring epithelial cells is dealt with using a Hamiltonian energy. This Hamiltonian is a function of the area of each cell of the Voronoi tessellation, but it involves also a pair-interaction that depends on the length of the common edge of two cells. The same kind of interaction (but between two types of cells) is also considered in [START_REF] Emily | A statistical approach to estimating the strength of cell-cell interactions under the differential adhesion hypothesis[END_REF]. Moreover, some geometric hardcore interactions are sometimes demanded. As an example, Lautensack and Sych ( [START_REF] Lautensack | 3D image analysis of open foams using random tessellations[END_REF]) modeled foams by a tessellation built from a Matern model with hardcore interaction. The resulting tessellation is then constrained to reach a desired regularity. The study of the regularity of the tessellation is also at the heart of the article of Eglen and Willshaw ( [START_REF] Eglen | Influence of cell fate mechanisms upon retinal mosaic formation: a modelling study[END_REF]): their work shows the relevance of forcing the geometry of cells in order to model retinal neurons.

It is thus natural to consider Gibbsian modifications of the Poisson-Voronoi or Poisson-Delaunay tessellation, involving a smooth interaction but also a hardcore interaction (in a general sense, see Definition 1), in order to produce more realistic models of interacting random structures.

A first mathematical model has been proposed by Ord (see the discussion in [START_REF] Ripley | Modelling spatial patterns (with discussion)[END_REF]). In this model, the interaction relies on each cell of the Voronoi tessellation. In particular, a classical hardcore interaction forces the cells not to be too small. This model can be viewed as a nearest neighbour Gibbs point process and was studied in [START_REF] Baddeley | Nearest-neighbor Markov point processes and random sets[END_REF]. Its existence on the infinite support R d is implied by the results in [START_REF] Bertin | Existence of "nearest-neighbour" spatial Gibbs models[END_REF] and [START_REF] Bertin | Existence of Delaunay Pairwise Gibbs Point Processes with Superstable Component[END_REF]. A Birth-Death simulation algorithm for simulating such nearest neighbour Gibbs processes is presented in [START_REF] Bertin | Spatial Delaunay Gibbs point processes[END_REF]. However, tessellations involving geometric hardcore interactions do not generally belong to a classical theoretical framework as the previous one. They are in general not hereditary in the sense that, when removing a point from an allowed tessellation, the resulting tessellation may become forbidden (see Section 4, or [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF], about this property). Consider for instance a generalization of the Ord process where the cells are forced not to be too large: this natural model is not hereditary. The existence of a Gibbs Delaunay-Voronoi tessellation on R d associated to a large class of possible non-hereditary interactions has been proved recently in [START_REF] Dereudre | Gibbs Delaunay tessellations with geometric hardcore conditions[END_REF] and [START_REF] Dereudre | Existence of Gibbsian point processes with geometry-dependent interactions[END_REF]. For these processes, no simulation algorithm has been presented so far.

From a statistical point of view, the issue is the estimation of the interaction. Assuming a parametric form, this may be achieved through the maximum likelihood or the pseudo-likelihood procedure. The maximum likelihood estimator suffers from a lack of theoretical justifications, except for restrictive examples of interacting point processes (see [START_REF] Møller | Statistical Inference and simulation for Spatial point Processes[END_REF] for a review), which do not concern tessellation models. On the other hand, some theoretical results are available for the pseudo-likelihood estimator. Consistency and asymptotic normality are proved in [START_REF] Billiot | Maximum pseudolikelihood estimator for exponential family models of marked Gibbs point processes[END_REF] in a general framework including some Gibbs tessellation models, but without any hardcore interaction. A generalization to interactions involving a possible non-hereditary hardcore part is considered in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]. In this article, the consistency of the estimation of both the hardcore part and the smooth part of the interaction is proved, in a setting concerning a large class of tessellation models.

In the present article, we rewrite these theoretical results, sometimes established in an abstract setting, to the framework of Gibbs Delaunay-Voronoi tessellations. Moreover, some theoretical complements are given. In particular, a Birth-Death-Move algorithm is presented to simulate Gibbs tessellations with non-hereditary interactions, and a convergence result is proved. We also extend the recent concept of residuals introduced in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF] to the non-hereditary setting. Nevertheless, the aim of this article is mainly to clarify how to handle (non-hereditary) Gibbs Delaunay-Voronoi tessellations in practice: what kinds of models are available? How should we simulate these tessellations? How should we fit them to a data set and validate the fitted model?

In the first part, the formal definition of Gibbs Delaunay-Voronoi tessellations is given. We restrict ourselves to tessellations on the plane R 2 for simplicity. Three example models are then considered: a non stationary crystallized triangulation model, a stationary interacting Delaunay model and a Voronoi tessellation model. We think that these models could be relevant for the biological applications cited before. Moreover, they are used throughout the article to illustrate the proposed methods. In Section 3, we explain how to simulate (non-hereditary) Gibbs Delaunay-Voronoi tessellations thanks to a Birth-Death-Move Metropolis-Hastings algorithm. Some simulations of the three above examples are presented. In Section 4, we consider the estimation issue. As explained there, the pseudo-likelihood approach is preferred to the maximum likelihood procedure for practical reasons. As a matter of fact, the maximum likelihood estimator is prohibitively time-consuming in our setting. However, if possible, maximum likelihood could be used in a second step to refine the pseudo-likelihood estimation. A procedure is presented to estimate both the hardcore parameters and the interaction, as considered in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]. Finally, the concept of residuals as recently introduced in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF] is generalized, which gives a method to validate the fitted model. In the appendix, we present some theoretical justifications. They concern the existence of Delaunay-Voronoi tessellations, the convergence of the simulation algorithm, and the consistency of the estimation procedure.

2. The Gibbs Delaunay-Voronoi tessellations model.

The Poisson Delaunay-Voronoi tessellations.

In paragraph 2.1.1, we recall the basic definition of point configurations. In 2.1.2, some regularity assumptions are given to ensure that the Delaunay-Voronoi tessellations are well-defined. Randomness is introduced in paragraph 2.1.3, via Poisson point processes, to define the wellknown Poisson Delaunay-Voronoi tessellations which are models of random tessellations without interaction between the cells. The interaction is introduced in Section 2.2.

Point configurations.

Let us denote by R 2 the 2-dimensional Euclidean real space. B(R 2 ) is the set of bounded Borel sets in R 2 . The state space M(R 2 ) is the set of regular locally finite point configurations γ in R 2 defined by

M(R 2 ) =    γ ⊂ R 2 such that -a) for all Λ in B(R 2 ), Card(γ ∩ Λ) < +∞ -b) four points of γ are not on a same circle -c) for every half plane H in R 2 , Card (γ ∩ H) > 0    , (1)
where Card(γ ∩ ∆) denotes the number of points from γ in the set ∆. Let γ be in M(R 2 ) and Λ a Borel set in R 2 , we denote by γ Λ the restriction of γ on Λ which is just the set γ ∩ Λ. For a point x in R 2 , we denote by γ + x the configuration γ ∪ {x} and if x belongs to γ, γx denotes the set γ\{x}.

Delaunay-Voronoi tessellations.

Let us recall the definition of Delaunay-Voronoi tessellations, which are given for example in [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF] Concerning the Voronoi tessellation coming from γ, for every x in γ, we define the Voronoi cell C(x, γ) by

C(x, γ) = z ∈ R 2 such that ∀y ∈ γ\{x} |z -x| ≤ |z -y| .
From points a) and c) in (1), we remark that C(x, γ) is a bounded closed convex set in R 2 . The Voronoi tessellation Vor(γ) is defined by the set of all C(x, γ) for x in γ. Vor(γ) is also a partition of the plane R 2 .

There are some relations between these two tessellations. Indeed, T = {x, y, z} in γ is a Delaunay triangle if and only if C(x, γ) ∩ C(y, γ) ∩ C(z, γ) = ∅.

Poisson Delaunay-Voronoi tessellations.

In this paragraph we define the Poisson Delaunay-Voronoi tessellations as in [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF] or [START_REF] Stoyan | Stochastic Geometry ans its applications[END_REF]. Let us recall that the space of point configurations M(R 2 ) is endowed with the σ-algebra σ(M(R 2 )) generated by the sets

{γ ∈ M(R 2 ), N Λ (γ) = n}, n ∈ N, Λ ∈ B(R 2 )
, where N Λ (γ) denotes the number of points of γ in Λ. The most prominent probability measures on M(R 2 ) are the Poisson processes. Let us denote them by π ν , where ν is a locally finite measure on R 2 and stands for the intensity measure (see [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF] page 83). When ν is equal to zλ (z > 0, λ the Lebesgue measure) we simply write π z which represents the classical stationary Poisson Point Process with intensity z. Let us remark that π ν is not necessary stationary but obviously π z is.

For every Λ in B(R 2 ), π ν Λ (respectively π z Λ ) denotes the Poisson process π ν (respectively π z ) restricted on Λ.

The law of Del(γ) (respectively Vor(γ)) under the process π ν is called the Poisson Delaunay (respectively Voronoi) tessellation with intensity ν. These Poisson Delaunay-Voronoi tessellations are well-studied in [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF], Section 4.

Random Delaunay-Voronoi tessellations with interaction.

This section is devoted to the presentation of interacting random Delaunay-Voronoi tessellations. The interaction is introduced, via Gibbs modifications of the Poisson Delaunay-Voronoi tessellations, by specifying the conditional densities. This is the classical strategy used in physics and biology (see for example [START_REF] Ruelle | Superstable interactions in classical statistical mechanics[END_REF]).

For every Λ in B(R 2 ), we consider the conditional density f Λ with respect to the Poisson process

π ν Λ defined by f Λ (γ Λ , γ Λ c ) = 1 Z Λ (γ Λ c ) e -EΛ(γΛ,γ Λ c ) , (2)
where γ Λ is a point configuration inside Λ, γ Λ c is a point configuration outside Λ and E Λ (γ Λ , γ Λ c ) is the energy of γ Λ given the outside configuration γ Λ c . E Λ is a functional from Vor(γ) or Del(γ) to R ∪ {+∞} which we will precise later. Z Λ (γ Λ c ) := e -EΛ(γΛ,γ Λ c ) π ν Λ (dγ Λ ) is the normalization constant in order to have a probability density under π ν Λ . Let us remark that the conditional densities favor the point configurations with low energy and conversely penalize the point configurations with high energy. If the energy E Λ (γ Λ , γ Λ c ) is equal to infinity then, with probability one, the configuration is even forbidden and does not appear. One classical example of such a situation is when the points of γ are prohibited from being closer than a distance R apart, id est E Λ (γ Λ , γ Λ c ) = +∞ if there exist x in γ Λ and y in γ such that |x -y| ≤ R. This constraint is usually designated as a hardcore interaction. In this paper, we generally call hardcore interaction any situation where E Λ (γ Λ , γ Λ c ) = +∞. Definition 1. An energy (or an interaction) is said to contain a hardcore part if E Λ (γ Λ , γ Λ c ) = +∞ for some γ and some Λ.

We denote by M ∞ (R 2 ) the set of allowed configurations which is defined by

M ∞ (R 2 ) = γ ∈ M(R 2 ) such that for all Λ in B(R 2 ), E Λ (γ Λ , γ Λ c ) < +∞ . (3)
Now let us define the model of random Delaunay-Voronoi tessellations with interaction.

Definition 2. A probability measure P on M ∞ (R 2 ) is a Gibbs Delaunay-Voronoi tessellation for the energies (E Λ ) Λ∈B(R 2 ) and the intensity measure ν if for every Λ in B(R 2 ) and for P -almost every outside configuration γ Λ c , the law of P given γ Λ c is absolutely continuous with respect to π ν with the density f Λ (., γ Λ c ).

This definition of Gibbs measures is the classical one that can be found for example in [START_REF] Preston | Random fields[END_REF]. Let us point out several problems about the existence of these Gibbs Delaunay-Voronoi tessellations. First of all there are some conditions on the energies (E Λ ) Λ∈B(R 2 ) to ensure that the conditional densities (f Λ ) Λ∈B(R 2 ) are well-defined and compatible. Moreover, even if it is the case, it is not obvious that Gibbs Delaunay-Voronoi tessellations exist and are unique (the non unicity of Gibbs processes is called phase transition in statistical mechanics). In Section 5.1, we give some conditions to ensure the existence of stationary Gibbs Delaunay-Voronoi tessellations for energy functions given by ( 4) and ( 5) below. As far as we know, uniqueness or phase transition for such Gibbs measures have not been proved. For similar Gibbs models with multi-type particles, a phase transition result is given in [START_REF] Bertin | Phase transition in the nearest-neighbor continuum Potts models[END_REF].

Since we are interested by models of random Delaunay-Voronoi tessellations, the energy functions have to depend on the local geometry of the Delaunay or Voronoi tessellations. We need some notations. Let us first define the neighbour relations ∼ Del , ∼ Vor between the cells in Del(γ) or Vor(γ).

For all T, T ′ ∈ Del(γ), T ∼ Del T ′ if Card(T ∩ T ′ ) ≥ 2. For all C, C ′ ∈ Vor(γ), C ∼ Vor C ′ if C ∩ C ′ = ∅.
In

Λ in B(R 2 ). A triangle T ∈ Del(γ) (respectively a cell C ∈ Vor(γ)) is outside Λ if for every configuration γ ′ Λ in Λ, T (respectively C) is in Del(γ Λ c ∪ γ ′ Λ ) (respectively Vor(γ Λ c ∪ γ ′ Λ )
). In other words, T (or C) is outside Λ if T (or C) remains in Del(γ) (or Vor(γ)) for any modification of the configuration γ inside Λ. T (or C) is inside Λ if it is not outside Λ. We denote by Del Λ (γ) (respectively Vor Λ (γ)) the cells T in Del(γ) (respectively C in Vor(γ)) which are inside Λ. a) A general form for the energy of a Delaunay tessellation We define the energy E Λ of the Delaunay tessellation by

E Λ (γ Λ , γ Λ c ) = T ∈DelΛ(γ) V 1 (T ) + {T,T ′ }⊂Del(γ) T ∼ Del T ′ T or T ′ in DelΛ(γ) V 2 (T, T ′ ), (4)
where V 1 is a function from the space of triangles T to R ∪ {+∞} and V 2 is a symmetric function from T 2 to R ∪ {+∞}. In Section 2.3, we give precise examples of functions V 1 and V 2 .

b) A general form for the energy of a Voronoi tessellation Similarly, we define the energy E Λ of the Voronoi tessellation by

E Λ (γ Λ , γ Λ c ) = C∈VorΛ(γ) V 1 (C) + {C,C ′ }⊂Vor(γ) C∼VorC ′ C or C ′ in VorΛ(γ) V 2 (C, C ′ ), (5)
where V 1 is a function from the space of bounded convex sets C to R ∪ {+∞} and V 2 is a symmetric function from C 2 to R∪{+∞}. In Section 2.3, a precise example of functions V 1 and V 2 is provided.

Three explicit reference models

We present in this section three explicit examples of Gibbs Delaunay-Voronoi models, that will be our reference models until the end of the paper. All the functions V 1 and V 2 given in this section satisfy the assumptions in Section 5.1 and so the associated Gibbs Delaunay-Voronoi tessellations exist.

Model 1: a non stationary crystallized triangulations model.

In this model, we propose to define a non stationary random triangulation in which the angles of triangles are forced to be larger than a fixed real α in [0, π 3 [. If α is chosen close to π 3 , the model produces rigid random triangulations. Moreover, it is possible to have a non stationary density of points. We assume that the intensity measure ν is absolutely continuous with respect to the Lebesgue measure λ and the energy E Λ is defined by (4) with

V 1 (T ) = +∞ if α(T ) ≤ α, 0 otherwise, and V 2 = 0,
where α(T ) is the minimal angle inside T . In fact, the energy E Λ is equal to plus infinity if there exists a triangle inside Λ which is too flat.

Model 2: a stationary interacting Delaunay model.

In this model, we propose to study an example of stationary Delaunay triangulation with interaction. It is a simple model in order to present different practical and theoretical aspects in this work (modeling, simulation, estimation). This model has not the ambition to be directly usable in physics or biology.

First of all, we fix the intensity measure ν to be equal to zλ. Via a geometric hardcore interaction, we force the edges not to be too small, the triangles not to be too large and via a smooth interaction the large perimeters of triangles are favored or penalized (depending on the sign of θ). More precisely, let 0 < ε < α and θ be in R, the energy E Λ is defined by (4) with

V 1 (T ) =    +∞ if l(T ) ≤ ε, +∞ if R(T ) ≥ α, θPer(T ) otherwise,
and

V 2 = 0,
where l(T ) is the minimal length of sides of T , R(T ) is the radius of the circumscribed ball of T and Per(T ) is the perimeter of T .

Model 3: a Voronoi tessellation model.

In this third example, the interaction is defined as far as possible to fit with the biological applications evoked in the introduction (geometric regularities of cells, interaction between cells), although other interactions could be chosen. We suppose that the model is stationary so we fix the intensity ν to be equal to zλ. The geometry of cells is controlled by a hardcore interaction V 1 which forces the cells not to be too small, too large or too flat. Moreover, we add a smooth interaction involving a competition between the volumes of neighbours cells.

Let 0 < ε < α, B > 1/(2 √
3) and θ be in R. The energy E Λ is defined by ( 5) with

V 1 : C → V 1 (C) =        +∞ if h min (C) ≤ ε, +∞ if h max (C) ≥ α, +∞ if h 2 max (C) ≥ BVol(C), 0
otherwise, and

V 2 : (C, C ′ ) → V 2 (C, C ′ ) = θ max(Vol(C), Vol(C ′ )) min(Vol(C), Vol(C ′ )) -1 1 2
, where h min (C) is the minimal distance between the center x of the cell C and the edges of the boundary of C. Similarly, h max (C) is the maximal distance between x and the edges of C (see Figure 1). The choice of the power 1 2 is arbitrary and may be changed. Nevertheless, it seems to lead to more realistic simulations. The parameter B controls the form of the cell: the smaller B, the more regular the cell. For instance, for a hexagonal cell, B = 1/(2 √ 3) ≈ 0.29. Let us remark that if θ is positive the interaction V 2 forces the neighbour cells to have the same volume. Conversely, if θ is negative it forces the neighbour cells to have different volumes. The sign of θ is crucial in this model.

Simulations

Gibbs tessellations on a finite window

In this section, we deal with the simulation of our models. First of all, let us remark that Gibbs Delaunay-Voronoi tessellations are processes on R 2 and so one has to restrict or approximate them inside a fixed window Λ. The most natural choice would be to simulate the restriction on Λ but it is almost impossible to do it. Therefore the common method is to consider the finite volume Gibbs approximations on Λ, which is the probability measure absolutely continuous with respect to the Poisson process π ν Λ with the density f Λ given in [START_REF] Baddeley | Nearest-neighbor Markov point processes and random sets[END_REF]. The outside point configuration γ Λ c is then fixed and chosen arbitrarily. Results in statistical mechanics show, in general, that the thermodynamic limits of these finite volume Gibbs measures, when the volume Λ goes to R 2 , are Gibbs measures (see [START_REF] Preston | Random fields[END_REF]). Therefore, finite volume Gibbs measures are good approximations of our models.

There are essentially three possibilities to fix the outside configurations. The first possibility is to consider the empty outside configuration but it is not usable in our context because it produces non bounded Delaunay-Voronoi cells and so non computable energies. The second possibility is to fix an outside point configuration which has to be specified explicitly. Again, this is not practical because the strong hardcore interaction, which appears in our three reference models, makes it difficult to find such a configuration. The third possibility is the periodic outside configuration which is built by periodization in the full plane R 2 of the random configuration inside Λ. We choose this approach because it seems relevant to deal with the hardcore problems coming from the boundary effects. Now, let us give the precise construction of the periodic finite volume Gibbs measure. The simulation window Λ is chosen as the square [0, 1] 2 . Any other window in R 2 can be considered in the same way. A simple rescaling procedure enables us to reduce to this case. For every point configuration γ in [0, 1] 2 , we denote by γ the periodic configuration on R 2 defined by γ =

i∈Z 2 τ i (γ), (6)
where τ i is the translation in R 2 with respect to the vector i.

In the sums over the cells in (4) or (5), we must ensure that each collection of periodic cells has a unique contribution in the computation of the periodic energies. A solution consists in selecting only the cells whom barycenters are in [0, 1] 2 . Therefore, for any Voronoi cells C (respectively Delaunay triangle T ), we denote by < C > (respectively < T >) the barycenter of the cell C (respectively triangle T ). Similarly, for any couple of Voronoi cells (C, C ′ ) (respectively triangles (T, T ′ )) we denote by < C, C ′ > (respectively < T, T ′ >) the barycenter of the set C ∪ C ′ (respectively T ∪ T ′ ) .

The periodic energy Ē(γ) associated to the energy of a Delaunay tessellation (4) is defined by

Ē(γ) = T ∈Del(γ) <T >∈[0,1] 2 V 1 (T ) + {T,T ′ }⊂Del(γ) C∼ Del C ′ <T,T ′ >∈[0,1] 2 V 2 (T, T ′ ). (7)
Similarly the periodic energy Ē(γ) associated to the energy of a Voronoi tessellation ( 5) is defined by

Ē(γ) = C∈Vor(γ) <C>∈[0,1] 2 V 1 (C) + {C,C ′ }⊂Vor(γ) C∼VorC ′ <C,C ′ >∈[0,1] 2 V 2 (C, C ′ ). ( 8
)
Now let us define the Gibbs process on [0, 1] 2 with periodic outside configuration.

Definition 3. The periodic Gibbs Delaunay-Voronoi tessellation P is the point process on [0, 1] 2 which is absolutely continuous with respect to the Poisson point process on [0, 1] 2 (denoted by π ν 0 ), where the density f is defined, for every γ in [0, 1] 2 , by

f (γ) = 1 Z e -Ē(γ) and Z = e -Ē(γ ′ ) π ν 0 (dγ ′ ).
Remark 1. In the case where the intensity measure and the energy functions are stationary, we know there exists at least one stationary Gibbs measure (see Theorems 1 and 2 in the appendix). However, non stationary Gibbs measures may exist too. This phenomenon is called the breakdown of symmetry (see [START_REF] Preston | Random fields[END_REF] 4.1, for instance). It is not proved theoretically that symmetry breakdown can occur for the models studied in this paper. Sometimes, this phenomenon can be observed via simulations. Nevertheless let us remark that this will not be possible with the periodic Gibbs models simulated here. Indeed, if the intensity measure and the energy functions are stationary, as in models 2 and 3, then the associated periodic Gibbs models are stationary too and the symmetry breakdown is not observable. Free or configurational boundary conditions should be investigated.

Algorithm of simulation

As explained above, to deal with the boundary effects, we choose to simulate periodic tessellations. For sake of brevity, we confuse γ and γ in this section, similarly we will use f instead of f . The window of simulation is Λ = [0, 1] 2 and we omit the indexation by Λ in the sequel. Moreover, the notation M ∞ (R 2 ) is abusively extended in this section to the periodic tessellations with finite energies.

There exist different algorithms to simulate finite volume Gibbs processes. Some perfect simulation algorithms have been developed (see [START_REF] Geyer | Constrained Monte Carlo maximum likelihood for dependent data[END_REF], [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF]) but they seem not to be really implementable in our context due to the strong rigidity of our hardcore models. So we make the choice to simulate them via the classical Birth-Death-Move Metropolis-Hastings algorithm, which we recall below (see also [START_REF] Møller | Statistical Inference and simulation for Spatial point Processes[END_REF]).

For x ∈ [0, 1] 2 , N (x, σ 2 ) denotes the Gaussian distribution on R 2 centered at x with covariance matrix diag(σ 2 , σ 2 ), where σ > 0. This law is the proposal density for moving a point. Note that if the moved point falls outside [0, 1] 2 (in step 5. below), it is replaced inside [0, 1] 2 by the periodic property. We assume that the intensity measure ν is absolutely continuous with respect to the Lebesgue measure, and we denote by g its density. In particular, in the stationary case, i.e. when ν = zλ with z > 0, then g is identically equal to z.

1. For γ 0 ∈ M ∞ (R 2 ), let n = card(γ 0 ).
γ 1 = γ 0 + x if b < f (γ0+x)g(x) (n+1)f (γ0) , γ 0 otherwise.
4. If a > 2/3 then generate x uniformly on γ 0 and set

γ 1 = γ 0 -x if b < nf (γ0-x) f (γ0)g(x) , γ 0 otherwise.
5. If 1/3 < a < 2/3 then generate x uniformly on γ 0 , generate y ∼ N (x, σ 2 ) and set

γ 1 = γ 0 -x + y if b < f (γ0-x+y) f (γ0)
, γ 0 otherwise.

6. Iterate from 1. where γ 0 ← γ 1 .

This algorithm can be refined: The probability of move, birth or death proposals may differ and may depend on x, similarly the law for choosing a point before killing it, adding it or moving it may be chosen properly (e.g. according to the intensity law ν). The idea of this procedure is that, starting from an allowed configuration γ 0 , the iterations converge to the realization of an invariant measure which is the Gibbs process we want to simulate. For classical Gibbs point processes, this convergence is proved for example in Section 7.3 of [START_REF] Møller | Statistical Inference and simulation for Spatial point Processes[END_REF]. In our setting, the convergence is not obvious. It mainly relies on the following connectivity property (see also Definition 5 in the appendix): from any allowed configuration γ, it is possible to reach another allowed configuration γ ′ thanks to an iterative birth-death-move procedure as above. Since a hardcore Gibbs tessellation may be very rigid, this property does not always hold (consider for instance the Delaunay tessellations where all the triangles are imposed to be almost equilateral). Yet, in most situations, the connectivity exists. Let us note that the moving step is crucial here, because it allows the connectivity of rigid tessellations that a simple birth-death procedure would not.

We show in the appendix that, under some assumptions, the algorithm converges. These assumptions are fulfilled for Model 2 when 2ε < α < 1 2 . For the Voronoi tessellation presented in Model 3, we claim that the convergence holds for a large reasonable set of parameters (ǫ, α, θ). However, in this case theoretical justifications become tedious and are not achieved in this paper.

Practical implementation

In the above algorithm, the choice of the initial configuration is crucial. We must start from an allowed configuration γ 0 in [0, 1] 2 , i.e. γ0 ∈ M ∞ (R 2 ). For this reason, we cannot start from the empty configuration. In our simulations, we chose to start from the point configuration whom Delaunay tessellation is a regular lattice of triangles (see its plot on top left of Figure 2). This starting configuration is allowed by all our hardcore models, provided the distance between points is properly chosen.

The computation of the ratio in steps 3.-5. of the algorithm is time-consuming since it supposes the computation of two tessellations plus a calculus on their cells to obtain their respective density. But it is possible to simplify this computation by focusing on a smaller window. Indeed, consider for instance step 3, i.e. the birth case (the same approach remains true for the death or the move step). When one adds a point x in a configuration, the new tessellation differs from the previous one only in a neighbourhood of x. Thus, the ratio of the densities reduces to a difference of energies in this neighbourhood. The size of this neighbourhood is determined by the size of cells around x. For instance in Model 2, the diameter of cells is forced to be smaller than 2α, so it suffices to focus on a ball with radius 2α around x to compute the difference of energies.

When there is a hardcore interaction, the convergence of the algorithm may be slow. Indeed, when one adds a point, the new tessellation can be forbidden. Hence, in presence of a strong hardcore interaction, most of the new tessellations proposed by the birth or death step in the algorithm will be refused. For this reason, we check the progress of the algorithm by monitoring control. Every one thousand iterations, we count the total number of points in the configuration as well as the number of accepted birth steps (among one thousand steps), the number of accepted death steps and the number of accepted move steps. The plot of theses numbers all along the iterative algorithm helps us to check the stabilization of the iterative process (though this is not a proof of the convergence to the invariant measure).

Some examples

We present some simulations of the models introduced in Section 2.3. We use the Birth-Death-Move algorithm presented before with σ = 0.015.

Figure 2 shows a simulated tessellation from Model 1 where α = π/6. Let us denote by g the density of the non-stationary intensity measure ν. It is 1-periodic with respect to each component, i.e. for all (x, y) ∈ R 2 , g(x + 1, y) = g(x, y) and g(x, y + 1) = g(x, y), and for (x, y) ∈ [0, 1] 2 , we assume g(x, y) = 100[(x -0.5) 2 + (y -0, 5) 2 ] -0.75 .

The triangles of such a tessellation are forced not to be too flat and to be more dense around the point (0.5, 0.5). The monitoring control seems to justify the convergence of the algorithm after 5.10 4 iterations. In Figure 3, two tessellations from Model 2 have been simulated with α = 0.08, z = 1000, and θ = ±5. We did not introduce in these simulations the hardcore parameter ǫ. When θ > 0 (on bottom), the tessellation is more likely to exhibit a small number of vertices since the total sum of perimeters is then low, which minimizes the energy. For θ < 0 (on top), this is the contrary: the energy is minimal when the total sum of perimeters is high, inducing a lot of vertices in the tessellation. In these two cases, the size of the triangles is controlled by the hardcore parameter α, which might be unnecessary when θ < 0 but is certainly a big constraint when θ > 0.

Model 3, involving Voronoi tessellation, has been simulated for α = 0.05, B = 0.625, z = 100, and θ = -0.8, -0.5, 0.5, 0.8. The hardcore parameters α and B force the cell not to be too large or too flat. We did not impose a minimal length of sides through the hardcore parameter ǫ. The parameter θ quantifies the dependence between the size of neighbour cells: when θ < 0 they are most likely to have different sizes, whereas for θ > 0 they tend to exhibit the same volume. These two opposite behaviors are clearly observable in the two extreme cases θ = -0.8 and θ = 0.8 in Figure 4. When |θ| = 0.5, this difference is more difficult to distinguish. A challenging task will be to properly estimate the parameters α, B, z, and θ in both the apparently closed situation θ = -0.5 and θ = 0.5. This problem is addressed in the next section.

In Figure 5, one can check the convergence of the algorithm for the four above simulations. Note that in the very rigid case θ = 0.8, the iteration process needs a lot of time (tens of thousands iterations) before starting a birth-death step. This shows the importance of the moving step.

More simulations of Model 3 are presented in Figure 6 with their monitoring control in Figure 7. They show the impact of B on the geometry of the cells. It plays a bigger role when θ < 0, since in this case we note the cells may be very flat without this hardcore parameter (B = +∞). Let us remark that some clusters of small cells appear when θ < 0. This is more visible when B is not too constraining.

Estimation

In this section, we focus on stationary Gibbs tessellations, in particular we suppose that the underlying Poisson process is π z , the stationary Poisson point Process with intensity zλ where λ is the Lebesgue measure. This is the case of Model 2 and 3 presented before. We apply our estimation procedure to the parameters z, ǫ, α, B and θ involved in these models. An interesting generalization would be to include the estimation of a non-stationary intensity (as in Model 1) to the estimation of the parameters of the interaction. We do not deal with this task in this paper.

Let us specify some notations. We assume that we observe a tessellation coming from a point configuration γ on a window Λ n = [-n, n] 2 . This tessellation is defined through an energy function E Λ (γ Λ , γ Λ c ) as in (4) or [START_REF] Bertin | Spatial Delaunay Gibbs point processes[END_REF]. In the following, we denote it by E β,θ Λ (γ Λ , γ Λ c ), since we assume a parametric form. We actually suppose that it depends on parameters (β, θ): β is the hardcore parameter (parameterizing the finiteness of the energy function), while θ is the smooth interaction parameter. For instance, in Model 2: β = (ǫ, α), θ = θ; in Model 3: β = (ǫ, α, B), θ = θ. This distinction is presented clearly in AppendixA.3, where some theoretical results for the asymptotic consistency of our estimators are given.

We consider a classical two-step estimation procedure. We first estimate β, then we estimate θ and z by pseudo-likelihood where β is replaced by its estimator.

The choice of the pseudo-likelihood approach instead of the classical maximum likelihood estimator is mainly imposed by practical reasons. Indeed, maximum likelihood requires the estimation, by simulations, of an unknown normalizing constant. This approach demands to simulate several tessellations according to the model, which is extremely time-consuming in the situation when a hardcore interaction is involved (see previous section). Moreover, the pseudo-likelihood procedure has the advantage of being asymptotically consistent for a large class of models (see [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]), which has not been proved for the maximum likelihood estimator in such a general setting. However, when the hardcore interaction is not too strong, the maximum likelihood estimation may constitute a second step to refine a pseudo-likelihood approach.

There is a major difficulty to overcome in order to implement the pseudo-likelihood estimation in our case: the hardcore interactions are not necessarily hereditary. An interaction is hereditary if, for every forbidden point pattern γ, then, for every point x, the configuration γ + x remains forbidden. This is equivalent to: for every allowed point configuration γ, then for every point x ∈ γ, the configuration γx remains allowed. In other words, an interaction is hereditary if one can remove any point from γ. This property concerns only the hardcore interaction. So every interaction involving no hardcore part is necessarily hereditary. The models presented in Section 2.3 are not hereditary. Indeed, if one removes a point from an allowed tessellation, the new tessellation may contain cells that are too large (for instance). As a consequence, we must modify the classical pseudo-likelihood contrast to take into account the so-called removable points as introduced in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF] (see Definition 4).

The two-step procedure

The first step consists in estimating the hardcore parameter β. Let us first assume that β is a one-dimensional parameter. We suppose the following inclusion

if β < β ′ then ∀Λ, E β,θ Λ (γ Λ , γ Λ c ) < +∞ ⇒ E β ′ ,θ Λ (γ Λ , γ Λ c ) < +∞. ( 9
)
In this case, a consistent estimator of β is

β = inf{β > 0, E β,θ Λn (γ Λn , γ Λ c n ) < +∞}. ( 10
)
If instead of ( 9), the converse implication holds, then it suffices to replace the infimum by a supremum in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF].

In the case of a multi-dimensional hardcore parameter β, we estimate each of its components as above.

For instance, for Models 2 and 3 presented in Section 2.3, Property ( 9) is satisfied by the hardcore parameters α and B, while the converse holds for ǫ. As a consequence, following [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF], natural estimators for these examples are (the notations are the same as in 2.3):

• For Model 2:

ǫ = min{l(T ), T ∈ Del Λn (γ)}, α = max{R(T ), T ∈ Del Λn (γ)}.
• For Model 3:

ǫ = min{h min (C), C ∈ Vor Λn (γ)}, α = max{h max (C), C ∈ Vor Λn (γ)}, B = max{h 2 max (C)/Vol(C), C ∈ Vor Λn (γ)}.
The second step consists in estimating the smooth interaction parameter θ and the intensity parameter z. We use the pseudo-likelihood procedure for the reasons explained before. To deal with the non-hereditary problem, we must introduce the concept of removable points. Definition 4. Let γ be in M(R 2 ) and x be a point of γ, then x is removable from γ if there exists

Λ ∈ B(R 2 ) such that x ∈ Λ and E β,θ Λ (γ Λ -x, γ Λ c ) < +∞. ( 11
)
The following proposition, proved in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF], gives a more intuitive approach and justifies the name of removable points.

Proposition 1. Let γ be in M ∞ (R 2 ) and x be a point of γ, then x is removable from γ if and only if γ -x is in M ∞ (R 2 ).
From the definition, it is clear that the property of being a removable point from γ depends only on the hardcore parameter β and not on θ or z. Thus, we denote by R β (γ) the set of removable points in γ.

The more rigid the tessellation, the less removable points there are. In particular, if there is no hardcore part in the interaction function, every point of γ is removable. In Figure 8, the removable points of previous simulations are encircled.

We are now in position to introduce the pseudo-likelihood contrast function, adapted to the non-hereditary case:

P LL Λn (γ, z, β, θ) = Λn z exp -h β,θ (x, γ) dx + x∈R β (γ)∩Λn h β,θ (x, γ -x) -ln(z) ,
where h β,θ (x, γx) is the local energy of x in γ, defined for every x ∈ R β (γ) by

h β,θ (x, γ -x) = E β,θ Λ (γ Λ , γ Λ c ) -E β,θ Λ (γ Λ -x, γ Λ c ), ( 12
)
where Λ is a set containing x as in Definition 4. Let us point out that h β,θ (x, γ) is just equal to h β,θ (x, (γ + x)x) and is always well-defined for γ in M ∞ (R 2 ). The parameters θ and z are estimated by minimizing P LL Λn , where the hardcore parameter β is replaced by its estimator β obtained in the first step:

(ẑ, θ) = argmin z,θ P LL Λn (γ, z, β, θ). ( 13
)
The consistency of this estimation procedure is considered in AppendixA.3.

Practical implementation

The optimization of P LL Λn requires the calculus of the local energy h β,θ (x, γ) for any x ∈ Λ n . This is the same calculus as the one needed in step 3 of the algorithm presented in Section 3.2 and, as explained there, it can be achieved by focusing on a window around x. Moreover, this computation requires the knowledge of γ Λ c , the configuration outside this window. To prevent boundary problems, it is actually necessary to compute P LL on a sub-window of the initial observation window Λ n . We denote abusively Λ n this sub-window in the following.

The derivative of P LL Λn with respect to z yields the following estimator for z:

ẑ = Λn exp -h β,θ (x, γ) dx card(R β (γ) ∩ Λ n ) , ( 14
)
where card(R β (γ) ∩ Λ n ) is the number of removable points from the observed point pattern γ in Λ n .

In the simple case where θ is a one-dimensional parameter and h β,θ (x, γ) is a sufficiently regular function, the minimization of P LL Λn in θ can be reduced to the determination of the root of an equation. Indeed, from the derivative of P LL Λn with respect to θ, we obtain in this case that θ is the solution of

z Λ ′ n ∂h β,θ ∂θ (x, γ) exp -h β,θ (x, γ) dx = x∈R β (γ)∩Λn ∂h β,θ ∂θ (x, γ -x), (15)
where

Λ ′ n = {x ∈ Λ n , γ + x ∈ M ∞ (R 2 )}.
Moreover, when h β,θ (x, γ) depends linearly on θ (as in Models 2 and 3), for all x such that γ

+ x ∈ M ∞ (R 2 ), ∂h β,θ ∂θ (x, γ) = θh β,1 (x, γ),
which simplifies equation ( 15) above.

From a practical point of view, we first estimate θ thanks to [START_REF] Honda | Description of cellular patterns by Dirichlet domains: the two-dimensional case[END_REF], where z is replaced by ẑ given by [START_REF] Geyer | Constrained Monte Carlo maximum likelihood for dependent data[END_REF]. Then we deduce ẑ by plugging θ into [START_REF] Geyer | Constrained Monte Carlo maximum likelihood for dependent data[END_REF]. In both these estimations, the involved integrals are approximated by Monte Carlo (this is the most time-consuming step of the estimation procedure).

Some examples 4.3.1. For Model 2

We implement the estimation procedure on simulations of Model 2. We do not introduce the hardcore parameter ǫ here. The estimation of α, θ and z has been done from 200 replications of Model 2 when α = 0.08, z = 1000 and θ = ±5, simulated as in Section 3.2. The results are shown in Figure 9 and 10. We have distinguished two cases: first estimating θ by supposing z = 1000 known, then estimating both θ and z. In this last case, one can note in the bottom right plot of these figures the closed relation between ẑ and θ. Although the models that we consider are well identifiable, it is not surprising to observe this closed relation: it is implied by the Euler's formula, which connects linearly the number of cells and the number of vertices in a tessellation (see 3.2.11 in [START_REF] Møller | Lectures On Random Voronoi Tessellations[END_REF]). Therefore, if z is under-estimated, θ will tend to be under-estimated as well, in order to respect this linear relation.

When θ = -5 (Figure 9), the simulated tessellations rely on about 1500 points and all of them are removable. The estimation of α when θ = -5 actually shows that this hardcore parameter is useless in this case: the cells of the tessellation naturally satisfy the hardcore condition. It is interesting to note that this misspecification does not affect the estimation of the smooth interaction parameter θ. The average of θ is about -5, while its standard deviation is 0.4 when z is known and 1.6 when z is estimated. The average of ẑ is 1002 and its standard deviation 145.

Example of tessellation α θ when z is known θ when z is estimated ẑ Scatterplot of ( θ, ẑ) When θ = 5 (Figure 10), the hardcore plays an important role in the model. It is well estimated with a standard deviation of 3.10 -2 . The standard deviation of θ is 0.3 when z is known and 1.9 when z is estimated. The average of ẑ is 1049 and its standard deviation 313. These estimations seem less accurate than when θ = -5. This certainly comes from the fact that, when θ = 5, our simulated tessellations on [0, 1] × [0, 1] rely only on 500 points. Most of these points are removable (more than 90%), as showed in the left example of Figure 8.

For Model 3

Two hundred replications of Model 3 where α = 0.05, B = 0.625, z = 100 and θ = ±0.5 have been simulated according to the algorithm presented in Section 3.2 (see Figure 4 for an example). As above, the hardcore parameter ǫ was not introduced here. The results of the estimations are shown in Figure 11 when θ = -0.5 and in Figure 12 when θ = 0.5. Two situations are considered, assuming z = 100 is known or not. The particularity of these simulated Voronoi tessellations is their rigidity. The hardcore interactions are strong, forcing the cells not to be too large (through α) neither too flat (through B). This is confirmed by the accuracy of their estimation in both cases θ = ±0.5 (see the histograms in Figures 11 and12). But, as a consequence, there are only a few removable points, making the estimation of the smooth interaction parameters more difficult. Yet, it appears from these simulations that, in spite of the apparent similarity of the tessellations when θ = -0.5 and θ = 0.5 (see Figure 4) and in spite of the few number of removable points, the estimation procedure is mostly available to properly distinguish them.

When θ = -0.5 (Figure 11), there are in average 45 removable points on 265 points. The estimation of θ remains correct: the average and the standard deviation of θ are respectively -0.52 and 6.4 10 -2 when z is known, and -0.56 and 14.5 10 -2 when z is estimated. The average of ẑ is 94 while its standard deviation is 45.

When θ = 0.5 (Figure 12), there are only 3.5 removable points in average on about 215 points. In this latest extreme case, some estimations of θ and z were even impossible since there were no removable points at all (in 5 percent of the simulations). This shows the limit of the estimation procedure in presence of a very rigid tessellation. The average of θ is 0.55 and its standard deviation is 22.8 10 -2 when z is not estimated. This is surprisingly reasonable in view of the few numbers of removable points. When both z and θ are estimated, the results become bad: the average of θ is 0.53 with a standard deviation of 48 10 -2 and the average of ẑ is 189 with a standard deviation of 345. Their joint distribution is plotted on bottom left of Figure 12. A zoom in is plotted on bottom middle, where more than 90% of the points are remaining. The last plot on bottom right shows the repartition of θ according to the number of removable points in the tessellation, when z = 100 is assumed to be known. There is a clear bias when the number of removable points is low. Since z = 100 is fixed, this low number of removable points is associated with a strong rigidity, so θ is most likely to be high. Moreover, the standard deviation of θ seems to decrease with the number of removable points. This is confirmed by Table 1 which contains, for a fixed number of removable points card(R β (γ)), the number of tessellations from our simulations having this number of removable points (named replications) and the standard deviation of θn calculated from these tessellations (denoted sd( θ)).

card(R β (γ)) 1 2 3

Analysis of residuals

When fitting a model to a data set, the analysis of the residuals is a standard way to check the quality of the model. The concept of residuals for spatial point processes is not simple. A general definition is proposed in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF], where the authors also present several diagnostic tools based on residuals. The definition relies on the Campbell equilibrium equation due to Nguyen and Zessin (see [START_REF] Nguyen | Integral and differential characterizations of the Gibbs process[END_REF]), where the Papangelou conditional intensity is involved. In our context, the Papangelou conditional intensity does not always exist, because the hardcore interactions are not necessarily hereditary (see Remark 2 in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]). Yet, a Campbell equilibrium equation still holds, provided we restrict the support to the set of removable points. 1).

Proposition 2. Let P be a stationary Gibbs Delaunay-Voronoi tessellation as defined in Definition 2. For every bounded non negative measurable function ψ from R 2 × M(R 2 ) to R, we have

E P   x∈R β (γ) ψ(x, γ -x)   = E P R 2 ψ(x, γ)e -h β,θ (x,γ) ν(dx) ,
where h β,θ is defined in ( 12) and E P denotes the expectation under P .

This proposition is proved in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]. From this equation, following [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF], we can define the innovation process, for any bounded set ∆ in R 2 and for every function ψ as above:

I ∆, ψ, h β,θ , ν = x∈R β (γ)∩∆ ψ(x, γ -x) - ∆ ψ(x, γ)e -h β,θ (x,γ) ν(dx).
The residuals are then defined as an estimation of the innovations:

R ∆, ψ, h β, θ, ν = x∈R β (γ)∩∆ ψ(x, γ -x) - ∆ ψ(x, γ)e -h β, θ (x,γ) ν(dx),
where ν is an estimation of the intensity measure ν, which is simply ẑλ in the stationary case.

This generalization of the residuals to the setting of possible non-hereditary interactions allows to perform several diagnostic plots. We refer to [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF] for a presentation of different relevant choices for ψ, and for some diagnostic tools. A smoothed version of the residuals is also proposed, leading to more appealing graphics. The main purpose of the residuals analysis is to check whether the fitted model is misspecified.

As an illustration, let us assess the effect of a misspecified model to the tessellation simulated in top left of Figure 6. It actually corresponds to a sample from Model 3 where α = 0.05, B = +∞, z = 100 and θ = -0.5. But we will improperly fit a stationary Poisson process to this sample, then we will fit Model 2 (where the interaction relies on the Delaunay triangulation). Finally the correct Model 3 will be fitted for a sake of comparison. Figure 13 represents the sample according to these three points of view.

Voronoi Points Delaunay

Figure 13: Voronoi tessellation (left) and Delaunay tessellation (right) from the same point configuration (middle), coming from a simulation of Model 3 (top right of Figure 6).

We consider the simple case when ψ = 1. This corresponds to the so-called raw residuals, which have the following form in the stationary case:

R ∆, 1, h β, θ, ẑ = card R β (γ) ∩ ∆ -ẑ ∆ e -h β, θ (x,γ) dx.
To check the fitted model, we use the QQ-plot diagnostic presented in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF]. It consists in comparing the empirical quantiles of the fitted residuals to the empirical quantiles of bootstrapped residuals.

If we fit a stationary Poisson process to the sample, we obtain an estimated intensity ẑ = 833. The raw residuals, computed on squares ∆ with side 0.01, are shown in top left of Figure 14. The same kind of residuals have then been computed on 100 simulated Poisson process with intensity ẑ. A QQ-plot of these residuals with a 95%-confidence interval is shown on top right of Figure 14, where the residuals of the original sample have been added (crosses). An example of raw residuals from a simulated Poisson process is represented in the middle of this plot. It appears that the residuals of our sample do not behave as those from the simulated Poisson point processes. The stationary Poisson model is then misspecified.

Similarly, if we fit Model 2 to the same sample, we obtain α = 0.055 and θ = 4.49 when z = 1000. The raw residuals for this model are shown in bottom left of Figure 14. We have bootstrapped residuals from 100 simulated samples from Model 2 with the same parameters. One example of such residuals is shown on bottom middle. The QQ-plot, in bottom right, shows that the original sample does not seem to follow Model 2.

Finally, Model 3 is fitted. The estimation gives α = 0.049, B = 97.8 and θ = -0.56 when z = 100. The same plots as before are represented in Figure 15. According to the QQ-plot, one should not reject the fitted model. Let us remark that the estimation of θ is rather bad for our sample: the error is -0.06 although other simulations shows that the standard deviation of the errors is about 0.02. This is the reason why the distribution of the residuals is on the edge of the confidence interval in the QQ-plot. The two residuals images represented on the left show that some residuals may be very negative on some squares (the black ones). This is confirmed by the dispersion of the lowest quantiles in the QQ-plot. Thus, the distribution of the residuals can not be Gaussian in this example. This differs from the asymptotic gaussianity of most residuals conjectured in [START_REF] Baddeley | Residual analysis for spatial point processes[END_REF].

AppendixA.

AppendixA.1. Existence of Gibbs Delaunay-Voronoi tessellations.

The existence results presented here are published in a more general setting in [START_REF] Dereudre | Gibbs Delaunay tessellations with geometric hardcore conditions[END_REF] and [START_REF] Dereudre | Existence of Gibbsian point processes with geometry-dependent interactions[END_REF]. They are slightly modified and simplified so that they suit better the setting of Gibbs Delaunay-Voronoi tessellations. We suppose that the energy functions have the forms (4) or [START_REF] Bertin | Spatial Delaunay Gibbs point processes[END_REF]. The three following assumptions H1-H3 are sufficient to define the conditional densities f Λ in [START_REF] Baddeley | Nearest-neighbor Markov point processes and random sets[END_REF].

H1 M ∞ (R 2 ) = ∅. For every γ in M ∞ (R 2 ) and every Λ in B(R 2 ), γ in M(R 2 ) is called a (r, Λ)-modification of γ (with r > 0) if there exist distinct y 1 , y 2 , . . . , y n in Λ satisfying |y i -x i | < r for every 1 ≤ i ≤ n (with γ Λ = {x 1 , x 2 , . . . , x n }) such that γ = {y 1 , y 2 , . . . , y n } ∪ γ Λ c . H2 M ∞ (R 2 ) is a locally open set in M(R 2 ) which means that for every γ in M ∞ (R 2 ), every Λ in B(R 2 ) there exists r Λ (γ) > 0 such that any (r Λ (γ), Λ)-modification of γ is in M ∞ (R 2 ).
H3 The interactions V 1 and V 2 are stable which means that there exists a constant

K > 0 such that V 1 ≥ -K and V 2 ≥ -K.
Now let us give a collection of assumptions used in the proof of the existence of Gibbs Delaunay-Voronoi tessellations. For every R > 0, we denote by γ R an infinite configuration in R-equilateral position. That means that any triangle in Del(γ R ) is equilateral with length of sides equal to R (see the initial configuration presented in Figure 2 for an example). H5 There exist R > 0, 0 < r < R/2 and A ≥ 0 such that for every Λ ∈ B(R 2 ), every (r, Λ)modification γ of γ R , and every T, T

H4 There exist K 1 > 0 and K ′ 1 > 0 such that, for all γ ∈ M ∞ (R 2 ) and all Λ ∈ B(R 2 ), Card(γ Λ ) ≥ K 1 Vol(Λ) -K ′ 1 .

Residuals when fitting a Poisson

′ ∈ Del(γ) with T ∼ Del T ′ , V 1 (T ) ≤ A and V 2 (T, T ′ ) ≤ A; (A.1) Respectively, for every C, C ′ ∈ Vor(γ) with C ∼ Vor C ′ , V 1 (C) ≤ A and V 2 (C, C ′ ) ≤ A. (A.2)
Now we are able to give a first existence theorem Theorem 1. There exists a stationary Gibbs Delaunay-Voronoi tessellation for any intensity ν = zλ (z > 0) and any energy functions (E Λ ) Λ∈B(R 2 ) satisfying assumptions H1, H2, H3, H4 and H5.

Assumptions H4 and H5 can be substituted by the following one.

H6 There exists A > 0 such that for every r > 0 we can find R > 2r such that for every Λ ∈ B(R 2 ) and every (r, Λ)-modification γ of γ R the property (A.1) or (A.2) holds.

We have the second following existence theorem.

Theorem 2. There exists a stationary Gibbs Delaunay-Voronoi tessellation for any intensity ν = zλ (z > 0) and any energy functions (E Λ ) Λ∈B(R 2 ) satisfying the Assumptions H1, H2, H3 and H6.

The proofs of Theorems 1 and 2 can be found in [START_REF] Dereudre | Existence of Gibbsian point processes with geometry-dependent interactions[END_REF]. They rely on entropy tools which are only available in the setting of stationary processes (i.e. ν = zλ). Concerning our three example models, the following corollary holds (the existence of Model 2 is also proved in [START_REF] Dereudre | Gibbs Delaunay tessellations with geometric hardcore conditions[END_REF]).

Corollary 1. In the stationary case, i.e. when the intensity measure ν is equal to zλ, Gibbs Delaunay-Voronoi tessellations for models 1, 2 and 3 exist.

Proof. First of all, assumptions H1, H2, H3 are obviously satisfied for the three models. Concerning Model 1, the existence is given by Theorem 2. Assumption H6 is proved by taking A = 0 and R large enough with respect to r such that any (r, Λ)-modification γ of γ R have Delaunay triangles with angles larger than α. Concerning Models 2 and 3, the existence is given by Theorem 1. Assumption H4 comes from the hardcore interaction which forces the cells to be not too large. The uniform bound in H5 is obvious if R and r are chosen such that any (r, Λ)-modification γ of γ R is in M ∞ (R 2 ).

AppendixA.2. Convergence of the algorithm

The Birth-Death-Move algorithm used in this paper is presented in [START_REF] Møller | Statistical Inference and simulation for Spatial point Processes[END_REF] page 115 where the convergence is proved in Proposition 7.7 if the associated Markov Chain is irreducible and aperiodic (see also [START_REF] Meyn | Markov chains and stochastic stability[END_REF]). In our setting, there is no problem with the aperiodicity since the probability that nothing happens during one step of the algorithm is positive. In general to prove the irreducibility, it is sufficient to point out that every configuration γ is linked by a finite number of algorithm steps to the empty configuration. In our case, it is not possible because there is a strong hardcore interaction and so the connection with the empty configuration is in general false. So we need the connectivity assumption H7 based on the following definition. Definition 5. γ and γ ′ in M ∞ (R 2 ) are connected if there exist n ≥ 0 and a sequence of configurations γ 0 , γ 1 , . . . , γ n-1 , γ n in M ∞ (R 2 ) such that γ 0 = γ, γ n = γ ′ and, for each 0 ≤ i ≤ n -1, γ i and γ i+1 differ only by one step of the algorithm (a birth, a death or a move).

H7 For every γ and γ ′ in [0, 1] 2 such that γ and γ′ , defined in [START_REF] Bertin | Phase transition in the nearest-neighbor continuum Potts models[END_REF], are in M ∞ (R 2 ), then γ and γ ′ are connected.

The deterministic connectivity assumption H7 and the flexibility assumption H2 on the space M ∞ (R 2 ) ensure that for all configurations γ, γ ′ in M ∞ (R 2 ) and every r > 0, the algorithm may generate from γ, with a positive probability and a finite number of steps, a (r, [0, 1] 2 )-modification of γ ′ (see the definition after H1). Irreducibility of the Markov chain follows and we have the following proposition.

Proposition 3. Under the assumptions H2 and H7, the Birth-Death-Move algorithm presented in Section 3.2 converges to P (see Definition 3) in total variation norm for P -a.s. every initial condition.

It seems difficult in our context to obtain rates of convergence, because the energy functions are not locally stable (the local stability means that |E Λ (γ ∪ {x}) -E Λ (γ)| is uniformly bounded with respect to Λ, γ and x). Moreover the space M ∞ (R 2 ) may be very complicated since there is no upper bound in general for the number of steps n in assumption H7.

Let us remark that assumption H7 is not easy to check. If H7 is not satisfied then the algorithm converges to the restriction of P on the connected component of the initial configuration γ 0 for the connection relation defined below (see definition 5). In this case, the limiting distribution may depend on the initial configuration. For Model 1, we can show that H7 is satisfied if α is small enough. We don't give the proof here but the scheme is essentially the same than in the following Proposition 4 which deals with Model 2. For Model 3, it is more complicated, we have not proved it but it seems satisfied if α and B are large enough. Proposition 4. In Model 2, if 2ε < α < 1 2 then the algorithm presented in Sections 3.2 converges to P .

Proof. According to Proposition 3, it suffices to show H2 and H7. Since H2 is obviously satisfied for Model 2, it remains to show H7.

Let γ and γ ′ be in [0, 1] 2 such that γ and γ′ are in M ∞ (R 2 ). To simplify the notations, we say that γ is in M ∞ (R 2 ) if γ is in M ∞ (R 2 ). We start the sequence by putting γ 0 = γ and we construct the sequence γ i by an algorithmic procedure.

In a first step (called saturation) we add points until there does not exist any ball with radius ε without points. More precisely, we test if there exists x in [0, 1] 2 such that γ0 ∩ B(x, ε) = ∅. If it is not the case, the saturation is finished. If it is the case we add the point x by a birth-step action and we put γ 1 = γ 0 + x. Then, we test again if there exists x in [0, 1] 2 such that γ1 ∩ B(x, ε) = ∅. If it is not the case the saturation is finished otherwise we put γ 2 = γ 1 + x. We go on like this until the saturation procedure stops which is always the case since [0, 1] 2 may contain only a finite number of points with a distance between them bigger than ε. Let us remark that this construction produces configurations in M ∞ (R 2 ). We denote by γ m1 the saturated configuration of γ.

In a second step, we add the points of γ ′ to γ m1 by the following way. Let x be a point of γ ′ . By definition of the saturation, the configuration γ m1 + x is not in M ∞ (R 2 ) since there exists at least one point y in γ m1 such that |x -ȳ| ≤ ε (ȳ is the periodic version of y such that ȳ ∈ B(x, ε)). If this point y is unique we use a move-step action to move y to x. So we put γ m1+1 = (γy) + x. If these points are non unique, they are removed (except one) by death-step actions and the last one is moved as above. We denote by γ m1+1 , . . . , γ m2 this sequence and we remark that these configurations are in M ∞ (R 2 ) since γ m1 is saturated and 2ε < α < 1 2 . Now we saturate again the configuration γ m2 as above and we add another point of γ ′ to γ m2 . We go on until we have added all the points of γ ′ and we denote by γ m3 the final configuration.

It remains to remove the points of γ m3 which are not in γ ′ . It is sufficient to apply several death-step actions since the obtained configurations are in M ∞ (R 2 ).

AppendixA.3. Consistency of the estimation procedure

Let us suppose that the energy function, defined in (4) and [START_REF] Bertin | Spatial Delaunay Gibbs point processes[END_REF], is parameterized by β and θ, and is denoted by E β,θ Λ . We first need to distinguish properly the hardcore parameter β from the other parameter θ. This is the purpose of the following assumption.

S1: For all γ ∈ M(R 2 ), for all β and for all θ and θ ′ ,

∀Λ ∈ B(R 2 ), E β,θ Λ (γ Λ , γ Λ c ) < ∞ ⇐⇒ E β,θ ′ Λ (γ Λ , γ Λ c ) < ∞.
Under S1, the support of the energy is parameterized by β only, and not by θ, which confirms that β is the hardcore parameter. This assumption is satisfied by Models 2 and 3 with β = (ǫ, α) and β = (ǫ, α, B) respectively.

The strong consistency (under any stationary Gibbs measure) of β defined in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF] and (ẑ, θ) defined in [START_REF] Farhadifar | The influence of cell mechanics, Cell-Cell interactions, and Proliferation on Epithelial Packing[END_REF] are established in Theorem 2 in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF], under some regularity assumptions. These assumptions have been checked for Model 2 in Proposition 5 in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]. Concerning Model 3, the assumptions could be checked in the same way excepted for assumption S3 involved in [START_REF] Dereudre | Campbell equilibrium equation and pseudo-likelihood estimation for non-hereditary Gibbs point processes[END_REF]. We have not succeeded to prove it but it seems true at least for B large enough.

  page 15. For a point configuration γ in M(R 2 ), a set of three points T = {x, y, z} belonging to γ is a Delaunay triangle in γ if the open circumscribed ball B(T ) of T does not contain any point of γ. The Delaunay tessellation Del(γ) is defined by the set of all Delaunay triangles T in γ. By points b) and c) in (1), Del(γ) is a partition of the space R 2 .

Figure 1 :

 1 Figure 1: Example of Voronoi cell C with center x and distances hmin, hmax.

2 .

 2 Draw independently a and b uniformly on [0, 1]. 3. If a < 1/3 then generate x uniformly on [0, 1] 2 and set

Figure 2 :

 2 Figure 2: Simulation of Model 1 with α = π/6 and g(x, y) = 100[(x-0.5) 2 +(y-0.5) 2 ] -0.75 (5.10 5 iterations). Top left: initial point configuration; Top right: monitoring control (from top to bottom: number of moved points, number of birth points, number of killed points and total number of points, pointed out every 1000 iterations); Bottom: final simulated tessellation.

θ = - 5 θ = 5 Figure 3 : 8 Figure 4 :

 55384 Figure 3: Simulation of Model 2 with α = 0.08, z = 1000 and θ = -5 (top), θ = 5 (bottom) after 2.10 5 iterations. Top left: final simulated tessellation when θ = -5; Top right: monitoring control when θ = -5 (with the same plots as for Figure 2); Bottom left: final simulated tessellation when θ = 5; Bottom right: monitoring control when θ = 5 (with the same plots as before).

Figure 5 : 5 B = 1 , 5 B 5 Figure 6 :

 551556 Figure 5: Monitoring control for the simulations of Model 3 presented in Figure 4, in the same order (from top left to bottom right: θ = -0.8, -0.5, 0.5, 0.8). They consist in the same plots as for Figure 2.

Figure 7 :

 7 Figure 7: Monitoring control for the simulations of Model 3 presented in Figure 6, in the same order (B = +∞ (left), B = 1 (right), θ = -0.5 (top), θ = 0.5 (bottom)). They consist in the same plots as for Figure 2.

Figure 8 :

 8 Figure 8: Removable points (encircled) from: The Delaunay tessellation simulated in Figure 3 where θ = 5 (left); The Voronoi tessellation simulated in Figure 4 where θ = -0.5 (right).

Figure 9 :

 9 Figure 9: Estimation of Model 2 when α = 0.08, θ = -5, z = 1000, from 200 replications.

Figure 10 :

 10 Figure 10: Estimation of Model 2 when α = 0.08, θ = 5, z = 1000, from 200 replications.

Figure 11 :

 11 Figure 11: Estimation of Model 3 when α = 0.05, B = 0.625, θ = -0.5, z = 100, from 200 replications.

Figure 12 :

 12 Figure12: Estimation of Model 3 when α = 0.05, B = 0.625, θ = 0.5, z = 100, from 200 replications. On top left: A typical tessellation which is estimated. Bottom right: Repartition of θ according to the number of removable points observed in the tessellation (see also Table1).

Figure 14 :Figure 15 :

 1415 Figure 14: Analysis of residuals for misspecified models.

  fact, T ∼ Del T ′ if T and T ′ have a common edge and C ∼ Vor C ′ if C and C ′ have a common edge at their boundary. Now let us define the cells in Del(γ) or Vor(γ) which are inside or outside a given bounded set

Table 1 :

 1 Standard deviation of θ according to the number of removable points, from replications of Model 3 with θ = 0.5.

	4	5	6	7	> 7

Acknowledgements

Simulations have been done on the computer servers supported by the research program MOSTAPAD-CPER 2007-2013, Project 7 "STIC et Calcul", of the Fédération de Mathématiques des Pays de Loire.


Frédéric Mouton

Non-tangential, radial and stochastic asymptotic properties of harmonic functions on trees * †

Keywords: harmonic functions -trees -Fatou theorem -random walks. † Math. Classif. : 05C05, 31C20, 31C35, 60J15, 60J50

For a harmonic function on a tree with random walk whose transition probabilities are bounded between two constants in (0, 1/2), it is known that the radial and stochastic properties of convergence, boundedness and finiteness of energy are all a.s. equivalent. We prove here that the analogous non-tangential properties are a.e. equivalent to the above ones.

We are interested in the comparison between some non-tangential asymptotic properties of harmonic functions on a tree and the corresponding radial properties, using analogous stochastic ones. We proved in a previous work [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF], under a reasonable uniformity hypothesis, the almost sure equivalence between different radial and stochastic properties: convergence, boundedness and finiteness of the energy. The probabilistic-geometric methods, adaptated from those we used in the setting of manifolds of negative curvature [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques en courbure négative[END_REF], were flexible and presumed to extend to the non-tangential case for trees.

A recent article [START_REF] Atanasi | The lusin area function and local admissible convergence of harmonic functions on homogeneous trees[END_REF] shows by combinatorial methods the equivalence of the three non-tangential corresponding properties in the particular case of homogeneous trees. It seems to be time to show explicitely that our methods give in a swift way the non-tangential results for general trees satisfying the uniformity hypothesis above.

We use our previous results to compare the non-tangential notions with the radial and stochastic ones: we prove on one hand that the stochastic convergence implies the non-tangential convergence in the section 3 and on the other hand that the non-tangential boundedness implies almost surely the finiteness of the non-tangential energy in the section 4. The notations are fixed in the section 1 and our main result is stated in the section 2.

Setting

Let us briefly fix the notations (for details see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]). We consider a tree (S, A) i.e. a non-oriented, locally finite, connected and simply connected graph with vertices in S and edges in A. We will use the usual notions of path, distance and geodesic path and note x ∼ y iff (x, y) ∈ A.

We also consider a transient random walk (X n ) n on S such that the transition probability p(x, y) > 0 iff x ∼ y. Denote by P x the distribution of the walk starting from x and by p n (x, y) the probability P x [X n = y] of reaching y from x in n steps.

The Green function G(x, y) = ∞ n=0 p n (x, y) is finite by transience. Denote by H(x, y) the probability of reaching y starting from x. If z is on the geodesic path [x, y], the simple connectivity implies H(x, y) = H(x, z)H(z, y) and G(x, y) = H(x, z)G(z, y).

(

) 1
If U ⊂ S, the Green function of U , defined on U × U , is the expectation of the number of times the walk starting from x hits y before exiting U . The Laplacian of a function

f on S is ∆f (x) = E x [f (X 1 )] -f (x). The function f is harmonic if ∆f = 0.
Let u be a fixed harmonic function. The stochastic energy of u is J * (u) = ∞ k=0 ∆u 2 (X k ) (non-negative terms). The events L * * , N * * and J * * are defined respectively by the convergence of (u(X n )) n , its boundedness and the finiteness of the stochastic energy. The Martingale theorem implies J * * ∼ ⊂ L * * (P x -almost sure inclusion) [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]. It is known since P. Cartier [START_REF] Cartier | Fonctions harmoniques sur un arbre[END_REF] that geometric and Martin compactifications agree and the random walk converges almost surely to a point of the boundary ∂S. The exit law starting from x is the harmonic measure µ x and µ = (µ x ) x is a familly of equivalent measures. Conditioning by Doob's method of h-processes gives probabilities P θ

x (ending at θ). Asymptotic events verify 0-1 law and we define sets

L * = θ ∈ ∂S|P θ x (L * * ) = 1 , N * = θ ∈ ∂S|P θ x (N * * ) = 1 , J * = θ ∈ ∂S|P θ x (J * * ) = 1
, which determine stochastic notions of convergence, boundedness and finiteness of the energy at θ ∈ ∂S. For θ ∈ L * , lim u(X n ) is P θ x -a.s. constant (independent from x) and called the stochastic limit at θ.

Fix a base point o. For θ ∈ ∂S, γ θ is the geodesic ray from o to θ and for c ∈ N, Γ θ c = {y ∈ S|d(y, γ θ ) ≤ c} is a non-tangential tube. Let u be a harmonic function. For c ∈ N, its c-non-tangential energy at θ is J θ c (u) = y∈Γ θ c ∆u 2 (y) and its radial energy at θ is

J θ (u) = J θ 0 (u) = ∞ k=0 ∆u 2 (γ θ (k)).
There is radial convergence, boundedness or finiteness of the energy depending wether (u(γ θ (n))) n converges, is bounded or has finite radial energy. There is nontangential convergence of u at θ if for all c ∈ N, u(y) has a limit when y goes to θ staying in Γ θ c . There is non-tangential boundedness (resp. finiteness of the energy) if for all c ∈ N, u is bounded on Γ θ c (resp. J θ c (u) < +∞).

Main result

We now suppose (H): ∃ε > 0, ∃η > 0, ∀x ∼ y, ε ≤ p(x, y) ≤ 1 2η, a discrete analogue of the pinched curvature for manifolds. It also forces at least three neighbors for each vertex, and ensures transience. We proved in [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]: Theorem 2.1 For a harmonic function u on a tree with random walk satisfying (H), the notions of radial convergence, radial boundedness, radial finiteness of the energy, stochastic convergence, stochastic boundedness, stochastic finiteness of the energy, are µ-almost equivalent.

We prove here the following theorem: Theorem 2.2 Under the same hypotheses, the notions of non-tangential convergence, non-tangential boundedness and non-tangential finiteness of the energy are µ-almost equivalent to the notions above.

Considering the trivial implications, it is sufficient to prove that stochastic convergence implies non-tangential convergence and non-tangential boundedness implies almost surely non-tangential finiteness of the energy.

Stochastic implies NT convergence

The first implication needs the following lemma due to A. Ancona in a general setting [START_REF] Ancona | Théorie du potentiel sur les graphes et les variétés[END_REF], but easily proved here by simple connectivity: Lemma 3.1 If (x n ) n is a sequence converging non-tangentially to θ ∈ ∂S, the walk hits P θ o -a.s. infinitely many x n .

Let us see how this lemma helps. Assume that the harmonic function u has a stochastic limit l ∈ R at θ but does not converge non-tangentially towards l at θ. There exists δ > 0 and a sequence (x n ) n converging non-tangentially to θ such that |u(x n ) -l| ≥ δ for all n. As the random walk (X k ) k hits P θ o -a.s. infinitely many x n by the lemma, one can extract a subsequence (X kj ) j such that |u(X kj ) -l| ≥ δ for all j. Hence, P θ o -almost surely, the function u does not converge towards l along (X k ) k which leads to a contradiction.

Le us now prove the lemma. Recall that the principle of the method of Doob's h-processes is to consider a new Markov chain defined by p θ (x, y) = K θ (y) K θ (x) p(x, y) where the Martin kernel K θ (x) is defined as lim y→θ G(x,y) G(o,y) (see for example [START_REF] Dynkin | Boundary theory of markov processes (the discrete case)[END_REF]). This formula leads to analogous fomulae for the p θ n and the associated functions H θ and G θ . Consider for a fixed n the projection y n of x n on the geodesic ray γ θ (see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]). As the random walk starting from o and conditioned to end at θ hits almost surely y n due to the tree structure, the strong Markov property gives

H θ (o, x n ) = H θ (y n , x n ) = K θ (xn) K θ (yn) H(y n , x n ). By definition of the Martin kernel, K θ (xn) K θ (yn) = lim y→θ G(xn,y)
G(yn,y) and G(x n , y) = H(x n , y n )G(y n , y) as soon as

y n ∈ [x n , y], so H θ (o, x n ) = H(x n , y n )H(y n , x n ).
The distance between x n and y n is bounded as (x n ) n converges non-tangentially to θ, hence the last product is bounded from below by a constant C > 0 using (H). By Fatou's lemma, the probability conditioned to end at θ of hitting infinitely many x n is not smaller than C and the asymptotic 0-1 law ensures that it equals 1, which completes the lemma's proof.

NT boundedness implies finite NT energy

Denoting N c = {θ ∈ ∂S| sup Γ θ c |u| < +∞} and J c = {θ ∈ ∂S|J θ c (u)
< +∞}, we will show that for all c ∈ N, N c+1 ∼ ⊂ J c , which will give the wanted result by monotonous intersection. Let us write N c+1 = N ∈N N N c+1 , where

N N c+1 = θ ∈ ∂S sup Γ θ c+1 |u| ≤ N .
By countability it is sufficient to prove that for all N ,

N N c+1 ∼ ⊂ J c . Let us fix N ∈ N. Denote Γ = θ∈N N c+1
Γ θ c and τ the exit time from Γ. As

M n = u 2 (X n ) - n-1 k=0 ∆u 2 (X k )
is a martingale (see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]), Doob's stopping time theorem for the bounded exit time

τ ∧ n gives E o [M τ ∧n ] = E o [M 0 ] = u 2 (o) ≥ 0, hence E o τ ∧n-1 k=0 ∆u 2 (X k ) ≤ E o u 2 (X τ ∧n ) .
As X τ ∧n is at distance at most 1 from Γ, it lies in a tube Γ θ c+1 where θ ∈ N N c+1 and |u(X τ ∧n )| ≤ N . When n goes to ∞, monotonous convergence (∆u 2 ≥ 0) and the desintegration formula (see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]) give then, for µ-almost all θ ∈ ∂S,

E θ o τ -1 k=0 ∆u 2 (X k ) < +∞.
Let us use a conditioned version of formula 2 from [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF], which will be proved later :

Lemma 4.1 For a function ϕ ≥ 0 on Γ and τ the exit time of Γ,

E θ o τ -1 k=0 ϕ(X k ) = y∈Γ ϕ(y)G Γ (o, y)K θ (y).
This lemma implies that for µ-almost all θ ∈ ∂S, y∈Γ ∆u 2 (y)G Γ (o, y)K θ (y) is finite. In order to get an energy, we will show that G Γ (o, y)K θ (y) is bounded from below using the two following lemmas. The first one is due to A. Ancona [START_REF] Ancona | Théorie du potentiel sur les graphes et les variétés[END_REF] but has a very simple proof in the present context of trees. The second one enables comparison between G Γ and G. ⊂ J c . The proof of that fact is the same as in the analogous radial proof [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF] which completes the theorem's proof.

Let us now prove the lemmas. Concerning lemma 4.1, using Fubini,

E θ o τ -1 k=0 ϕ(X k ) = ∞ k=0 E θ o ϕ(X k )1 (k<τ ) .
The random variable ϕ(X k )1 (k<τ ) being measurable with respect to the σalgebra generated by (X i ) i≤k (see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]) and using formula 2 from [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF], the expectation above equals

∞ k=0 E o ϕ(X k )1 (k<τ ) K θ (X k ) = E o ∞ k=0 ϕ(X k )1 (k<τ ) K θ (X k ) = y∈Γ ϕ(y)G Γ (o, y)K θ (y),
which finishes the proof of lemma 4.1.

Let us prove lemma 4.2. Denote π(y) the projection of y on γ θ (see [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]) and remark that for z ∈ (π(y), θ), G(o, z) = H(o, π(y))G(π(y), z) and G(y, z) = H(y, π(y))G(π(y), z) by formula 1. Hence G(y,z) G(o,z) = H(y,π(y)) H(o,π(y)) does not depend anymore on z and its limit when z goes to θ is then K θ (y) = H(y,π(y)) H(o,π(y)) . By formula 1, G(o, y)K θ (y) = H(y, π(y)) G(o, y) H(o, π(y))

= H(y, π(y))H(π(y), y)G(y, y).

Lemma 4 . 2 Lemma 4 . 3

 4243 ∀c ∈ N, ∃α > 0, ∀θ ∈ ∂S, ∀y ∈ Γ θ c , G(o, y)K θ (y) ≥ α. For U ⊂ S containing Γ θ c and τ the exit time of U , lim y∈Γ θ c ,y→θ G U (o, y) G(o, y) = P θ o [τ = +∞]. By lemma 4.2, for µ-almost all θ ∈ N N c+1 , y∈Γ θ c ∆u 2 (y) G Γ (o, y) G(o, y) < +∞. If we show that for µ-almost all θ ∈ N N c+1 , P θ o [τ = +∞] > 0, lemma 4.3 gives N N c+1 ∼

But G(y, y) ≥ p 2 (y, y) ≥ 3ε 2 and H(y, π(y))H(π(y), y) ≥ ε 2c by (H) and d(y, π(y)) ≤ c, which finishes the proof of lemma 4.2.

Let us prove lemma 4.3 :

and by definition of Martin's kernel, if we could switch the limit and expectation, by a conditioning formula [START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF],

We now justify that inversion by Lebesgue's theorem. The idea is to bound, when τ is finite, G(Xτ ,y) G(o,y) by a multiple of K θ (X τ ). We compare for that purpose G(X τ , y) with K θ (X τ ). Denote again by π the projection function on γ θ . We

, by formula 1 and the remark that this formula also implies by definition of K θ and by taking the limit that

. We also have, by definition and formula 1,

We know that G is bounded (see [START_REF] Picardello | Finite truncations of random walks on trees[END_REF][START_REF] Mouton | Comportement asymptotique des fonctions harmoniques sur les arbres[END_REF]) and H is a probability, so it just remains to compare H(o, π(y)) with G(o, y). But H(o,π(y)) G(o,y) = (G(π(y), y)) -1 and 1 G is bounded by 1 3ε 2 . Merging the two cases gives a constant β such that G(Xτ ,y) K θ (Xτ ) ≤ βG(o, y), which enables to use Lebesgue's theorem and completes the proof of lemma 4.3.