Guy Bouchitt

Rajesh Mahadevan

MEAN FIELD THEORY FOR A GENERAL CLASS OF SHORT-RANGE INTERACTION FUNCTIONALS

Keywords: empirical measures, non-local functionals, Γ-convergence, mean-field energy, sub-additivity Mathematics Subject Classification: 49J45, 49K21, 49N15, 60B10, 70-10, 82B21

 obtained in the homogeneous case (r) = r -s where s > d.

Introduction

We consider a repulsive interaction function on (R d ) N of the kind (1.1)

c ε N (x 1 . . . , x N ) = i =j |x i -x j | ε .
where:

• N is the number of particles in R d ;

• ε > 0 scales the interaction distance between particles.

• the two-particle cost : [0, +∞] → [0 + ∞] satisfies:

(H1) is l.s.c. and (0) > 0 ( (0) = +∞ is allowed) (H2) ∃r 0 ≥ 0 such that is finite and non increasing on [r 0 , +∞) and lim r→∞ (r) = 0.

Date: October 24, 2023. 1991 Mathematics Subject Classification. 00X99. Part of this work was prepared at FCFM, Univ. Concepci'on. The first author warmly thanks this institution for its financial support and hospitality during the fall of 2022.

The second author acknowledges the support of Dirección de Postgrado, UdeC for the funding received through the project UCO 1866 .

In the whole paper, we denote by Ω the closure of a smooth domain Ω ⊂ R d where the N particles are located. We speak of a confined system when Ω is compact (container). Given a continuous exterior potential U : Ω → R, we consider the finite dimensional problem:

(1.2) E ε N (Ω, U ) := inf h N c ε N (x 1 . . . , x N ) + 1 N N i=1 U (x i ) : x i ∈ Ω ,
where h N is a suitable chosen normalization factor. Since the seminal work of Choquet in 1958 [START_REF] Choquet | Diamètre transfini et comparaison de diverses capacités[END_REF] and the growing interest of the quantum and statistical mechanics community, a lot of work has been devoted to the limit behavior of E ε N (Ω, U ) as N → ∞ (ε fixed) as well as the characterization of the weak cluster points of the empirical measures associated with N -point configurations of minimal energy. The cornerstone of the mean field theory consists in identifying a limit energy functional on measures whose minimizers are precisely these cluster points.

1.1. State of the art. The scaling factor h N in (1.2) must be selected so that the limit of the infimum belongs to (0, +∞). In turn this issue relies heavily on the integrability properties of the function g(x) = (|x|). Let us report on two cases of major interest:

1.1.1. Long range interaction case. Here Ω = R d and we take ε = 1. Moreover in addition to (H1)(H2), we assume that g ∈ L 1 loc (R d ) i.e.:

(1.3) 1 0 r d-1 (r) dr < +∞.

In that case, a relevant choice is h N = 1 N 2 meaning roughly that the interaction energy is averaged over all pairs of distinct points in {x 1 , x 2 , . . . x N }. The identification of the mean-field energy is well known in the case of Riesz potentials (r) = 1 r s for 0 < s < d, in the Logarithmic case (r) = -log(r) for d = 2 and more generally for of positive type i.e. such that the Fourier transform of g(x) = (|x|) is positive in R d (see for instance the monograph by S. Serfaty [START_REF] Serfaty | Coulomb gases and Ginzburg-Landau vortices[END_REF]). It is given by a non-local functional, the so called Direct energy:

(1.4) D (ρ) := (|x -y|) ρ ⊗ ρ(dxdy).

Accordingly the limit problem associed with (1.2) reads:

E ∞ (Ω, U ) = inf D (ρ) + U dρ : ρ ∈ P(R d ) ,
where the infimum is reached at a unique configuration provided U growths suitably at infinity. At this stage, a few comments are in order:

-in the case of a confining external potential U , there are several impressive works devoted to the next order asymptotics [START_REF] Serfaty | Coulomb gases and Ginzburg-Landau vortices[END_REF][START_REF] Petrache | Next order asymptotics and renormalized energy for Riesz interactions[END_REF][START_REF] Serfaty | Systems of points with Coulomb interactions[END_REF][START_REF] Cotar | Next-order asymptotic expansion for nmarginal optimal transport with coulomb and riesz costs[END_REF] in the case of Riesz potentials (r) = r -s for d ≥ 3 and d -2 ≤ s < d revealing an asymptotic behavior as N → ∞ of the form:

(1.5) lim

N →∞ N 1-s d (E N (Ω, U ) -E ∞ (Ω, U )) = C(s, d) (ρ U ) 1+ s d ,
where ρ U is the unique minimizer realizing E ∞ (Ω, U ).

-if U remains bounded at infinity (for instance a Coulomb potential vanishing at infinity), the existence of a minimizer ρ U may fail due to a loss of mass at infinity along minimizing sequences. A relaxation procedure leads to consider minimizers in the class of sub-probabilities ρ ∈ P -(R d ) and involves the weak* lower semicontinuous convexification of the Direct energy D. If is of positive type, this relaxed energy coincides with the natural 2-homogeneous extension of D to P -(R d ) while almost nothing is known if is merely locally integrable. For further details and examples of relaxed minimizers, we refer to the recent paper [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF].

1.1.2. Short range interaction case. Following an idea developed for the hard spheres model [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF], we look now at ε as a small parameter tending to zero with a prescribed speed as N → ∞. If one thinks to a container Ω of unit volume and ε to be the average distance of a particle to the others, it is natural to consider an asymptotic analysis where the product ε d N remains constant or converges to a given intensity factor κ ∈ (0, +∞). In a crowd model, this factor κ is related to a congestion ratio (see the hard speres model in Section 3.4 and Remark 3.13). This of course means that we need to assume that ε ∼ N -1 d . Accordingly, in order to obtain a precise scaling for h N ensuring a non-trivial behavior of the infimum (1.2), it is crucial to make an additional integrability on at infinity namely

(H3) +∞ r 0 (r)r d-1 dr < +∞.
It turns out that, under (H3), the right scaling factor in (1.2) is h N = 1 N in contrast with the long range case. This covers the case of hyper singular Riesz potentials (r) = r -s with s > d. For such potentials the parameter ε can be dropped thanks to the homogeneity and the normalized interaction energy becomes 1

N 1+ s d c 1 N (x 1 . . . , x N ).
Under the latter scaling, it was proved recently [START_REF] Hardin | Asymptotic properties of short-range interaction functionals[END_REF][START_REF] Douglas | Large deviation principles for hypersingular Riesz gases[END_REF] that the mean field energy is a local funtional defined on absolutely continuous measures

ρ = u L d Ω by F (ρ) = C(s, d) Ω u 1+ s d dL d , being C(s, d
) a universal constant. However, extending this result to more general costs seems to be difficult in the framework developed in [START_REF] Hardin | Asymptotic properties of short-range interaction functionals[END_REF], except possibly if is very close to a power potential. 1.2. Our contribution. This paper proposes a significant simplification of asymptotic analysis in the short range case. The approach is based on two components: first, treating the interaction distance ε in (1.1) as an infinitesimal parameter, and second, using an -counterpart of the traditional empirical measure frequently utilized in mean-field theory.

Thus, for every cost satisfying (H3), we can determine the mean field energy in terms of a local integral functional of the type Ω f (u) dx, where u = dρ dx denotes the local particle density and f is a convex integrand that exhibits super-linear growth at infinity. This expands upon previous findings [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF] that were limited to the hard-spheres model (where (r) = +∞ when r < 1 and (r) = 0 otherwise). Similarly, this permits to handle the case of hyper-singular Riesz potentials (r) = r -s for s > d analyzed in [START_REF] Hardin | Asymptotic properties of short-range interaction functionals[END_REF][START_REF] Douglas | Large deviation principles for hypersingular Riesz gases[END_REF].

It is noteworthy that fulfilling the integrability condition (H3) is crucial and cannot be sidestepped. When assuming that Ω is bounded, a cost that satisfies +∞ r 0 r d-1 (r) dr = +∞ would give rise to an infinite limit in (1.2) if the scaling is by

h N = 1 N and ε d N ∼ 1 (see Remark 3.14).
1.3. Setting of the asymptotic problem and notations. From now on, Ω will be a bounded domain of R d with Lipschitz boundary (∂Ω needs to be L d -negligible) and we consider a cost function which satisfies the standing assumptions (H1), (H2) and (H3).

For purposes of presentation, we utilize the infinitesimal length ε as the main parameter while the number of particles N = N ε approaches infinity as ε tends towards zero, following the scale N ε ∼ κ ε -d where κ is a positive constant. Later, we will establish that assigning a value to κ is not required because a uniform bound on the N ε -point interaction energy will automatically result in lim sup ε→0 N ε ε d < +∞. We can now incorporate measures in Ω into a variational framework for addressing the mean field problem. For each finite subset S ⊂ Ω, we define its ε-scaled empirical measure as follows:

(1.6) ρ ε S := ε d x∈S δ x .
This measure belongs to the set of non-negative Borel measures on Ω, denoted by M + (Ω). Here, ρ represents the total mass, which may be infinite, of any element ρ ∈ M + (Ω). Through this, we can observe that ρ ε S = ε d (S), which could deviate from the classical empirical measure of S ε with a total mass equal to 1. A key avantage of this approach is that it allows to avoid the non-local constraint that all competitors must belong to the subclass P(Ω) of probability measures.

Next, we define the ε-scaled interaction energy of a discrete set S ⊂ R d as follows:

(1.7)

ξ ,ε (S) = (x,y)∈S 2 \∆ |x -y| ε where ∆ := {(x, x) : x ∈ R d }.
We will refer to the interaction energy corresponding to ε = 1 as the "ground interaction energy", denoted by ξ . When the cost function is fixed, we will use ξ ε instead of ξ ,ε . With this in mind, we can define a scaled energy functional F ε : M + (Ω) → [0, +∞] for every ε > 0 in the following way:

F ε (ρ) = ε d ξ ε (S) if ∃S ⊂ Ω such that ρ = ρ ε S +∞ otherwise. (1.8)
The discrete problem (1.2) for h N = 1 N and N ε d ∼ κ can then be expressed through the relation:

κ E ε N (Ω, U ) ∼ inf F ε (ρ) + Ω U dρ : ρ ∈ M + (Ω) as ε → 0.
Accordingly the mean-field energy will be represented by a functional F : M + (Ω) → [0, +∞] characterized by the property that, for every U ∈ C(Ω), one has the convergence of infima

inf F ε (ρ) + Ω U dρ → inf F (ρ) + Ω U dρ
accompanied by the tight convergence of minimizers. This falls squarely within the Γ-convergence theory ([2, 11, 7]) on which we rely to support our results.

The paper is organized as follows: in Section 2, we establish a lower bound which allows to obtain the strong equi-coercivity of the sequence (F ε ); in addition we show that any weak* cluster point of a sequence (ρ ε ) with uniformly bounded energy is absolutely continuous with respect to the Lebesgue measure; in Section 3, we state the Γ-convergence of F ε as ε → 0 to a convex functional of the form F (ρ) = Ω f ( dρ dx ) dx, where the effective integrand f growths at least quadratically at infinity. It is given by the thermodynamical limit of a subadditive set funtion (Krengel's theorem). Some examples and applications are given. The Section 4 is devoted to the proof of the main theorem.

Notations:

-B(x, r) is the open ball of the Euclidean space R d centered at x and of radius r ; if x = 0, we simply denote B r ; -

Q k denotes the hypercube [-k/2, k/2) d , Q(x 0 , r) := x 0 + r Q 1 ; -∆ = {(x, x) : x ∈ R d }
µ h = 1, µ h → 1, µ = µ .
-To any non-empty set A, we associate the functions:

1 1 A (x) = 1 if x ∈ A 0 otherwise , χ A (x) = 0 if x ∈ A +∞ otherwise .

Energy estimates and compactness.

We begin with some elementary properties of the set function ξ (S) (defined by (1.7) for ε = 1). Lemma 2.1. Let S 1 , S 2 be finite disjoint subsets of R d . Then we have:

(i) (super-additivity) ξ (S 1 ∪ S 2 ) ≥ ξ (S 1 ) + ξ (S 2 ). (i) (sub-additivity at large distance) ξ (S 1 ∪ S 2 ) ≤ ξ (S 1 ) + ξ (S 2 ) + 2 + (η) (S 1 ) (S 2 ),
where η := dist(S 1 , S 2 ).

Proof. Since S 1 and S 2 are non-intersecting, we can split (S 1 ∪ S 2 ) 2 in four disjoint pieces as follows:

(S 1 ∪ S 2 ) 2 = (S 1 × S 1 ) ∪ (S 2 × S 2 ) ∪ (S 1 × S 2 ) ∪ (S 2 × S 1 ).
The inequality (i) is then straightforward whereas, for the (ii), we simply majorize by + (η) the contribution (|x -y|) of each pair (x, y) in (S 1 × S 2 ) ∪ (S 2 × S 1 ) where |x -y| ≥ η holds.

Recalling the notations (1.6) and (1.8) given in the introduction, we may rewrite a ε-rescaled version of Lemma 2.1 as follows: for every ε > 0 and all pairs of non-intersecting subsets (S ε , S ε ), one has (2.1)

F ε (ρ ε ) + F ε (ρ ε ) ≤ F ε (ρ ε + ρ ε ) ≤ F ε (ρ ε ) + F ε (ρ ε ) + 2 + (η ε ) ε d ρ ε ρ ε ,
where

ρ ε := ρ ε S ε , ρ ε := ρ ε S ε and η ε = ε -1 dist(S ε , S ε ) .

2.1.

A fundamental lower-bound. Since we are only assuming that (r) is decreasing for suitably large r, we need to define :

+ (r) := sup { (s) : s ≥ r} , - (r)
:= inf x,y∈[0,r] d (|x -y|) (2.2)
which are monotone non-increasing and satisfy:

+ (r) = (r) ∀r ≥ r 0 , -(r) = inf (s) : s ≤ r √ d ≤ + (r √ d).
Next, in the same line as in the survey [START_REF] Mathieu Lewin | Coulomb and Riesz gases: the known and the unknown[END_REF], we derive a very simple but fundamental lower bound for ξ (S) when S is a N -point system contained in a Borel subset B of finite volume in R d . For such a B and any δ > 0, we denote by m δ (B) the minimal number of disjoint δhypercubes

Q j = x j + [-δ/2, δ/2[ d such that B ⊂ ∪ m δ (B) j=1 Q j . Then it's easy to check that m δ (B) ∼ δ -d L d (B) as δ → 0.
For such subsets B, we will often use the following equivalent version obtained by keeping δ fixed while using large dilations:

(2.3) lim ε→0 ε d m δ B ε = δ -d L d (B). Lemma 2.2. Let B ⊂ R d be a Borel subset such that L d (B) < +∞.
Then for any N -point system S ⊂ B and any δ > 0, we have

(2.4) ξ (S) ≥ N -(δ) (ζ -1) + where ζ = N m δ (B) . Proof. Let {Q j , 1 ≤ j ≤ m δ (B)} be a covering of B by disjoint δ- hypercubes Q j and denote n j = (S ∩ Q j ).
We have N = n j while, by the super additivity of ξ (see Lemma 2.1) and the definition of -(δ), we have:

ξ (S) ≥ m δ (B) j=1 ξ (S ∩ Q j ) ≥ -(δ) m δ (B) j=1 n j (n j -1). If N > m δ (B) (i.e. ζ > 1)
, the desired lower bound (2.4) follows by noticing that the infimum inf

   m δ (B) j=1 t j (t j -1) : m δ (B) j=1 t j = N, t j ∈ R    is reached for t j = N m δ (B) = ζ , ∀j. If ζ ≤ 1 the inequality (2.4) is trivial since ξ (S) ≥ 0.
Remark 2.3. By the lower semicontinuity assumption (H1), since -(0) = (0) ∈ (0, +∞], we may always find a δ > 0 such that -(δ) > 0. Note that the inequality (2.4) is still valid if -(δ) = +∞ provided we agree that 0×+∞ = 0 (this situation occurs when = +∞ on an interval [0, δ 0 )).

Strong coercivity and compactness.

We recall the definition of the functional F ε given in (1.8) whose domain consists of ε-empirical measures of discrete subsets S ε ⊂ Ω (see the definition (1.6)).

Lemma 2.4 (strong coercivity). Assuming that satisfies (H1) (H2), let δ > 0 be such that -(δ) ∈ (0, +∞] and set α :=

δ d -(δ) L d (Ω)
.

Then we have:

(2.5) lim inf ε→0 F ε (ρ ε ) ρ ε 2 ≥ α , whenever (ρ ε ) is a sequence such that ρ ε → +∞.
Proof. Without loss of generality, we may assume that

F ε (ρ ε ) < +∞. Thus ρ ε = ρ ε Sε for a suitable N ε -point system S ε ⊂ Ω while: F ε (ρ) = ε d ξ ε (S ε ) = ε d ξ S ε ε , ρ ε = N ε ε d .
By applying the lower bound (2.4) to the subset ε -1 S ⊂ ε -1 Ω, we get:

(2.6) F ε (ρ ε ) ≥ -(δ) ρ ε ρ ε β ε -1 + where β ε = ε d m δ Ω ε .
From (2.3), we know that β ε → δ -d L d (Ω). Hence (2.5) follows by dividing (2.6) by ρ ε 2 and passing to the limit as ε → 0.

Proposition 2.5 (compactness). Assuming that satisfies (H1) and (H2), let U : Ω → R be a bounded Borel function and

(ρ ε ) a sequence in M + (Ω) such that (2.7) sup ε F ε (ρ ε ) + U dρ ε < +∞.
Then :

(i) there exists a constant C such that

ρ ε + F ε (ρ ε ) ≤ C < +∞ for every ε > 0. (ii) any weak* cluster point of (ρ ε ) is of the form ρ = u L d Ω with u ∈ L 1 (Ω).
Remark 2.6. The proposition above implies that sequences of point configurations S ε with equi-bounded energies admit a finite limiting intensity factor κ := lim sup ε ε d (S ε ) while, by the assertion (ii), S ε is not allowed to concentrate anywhere as ε → 0.

Proof. Suppose that ρ ε has no upper bound. Then the uniform energy upper bound (2.7) implies that

lim inf ε→0 F ε (ρ ε ) ρ ε ≤ sup Ω |U |,
while by (2.5) the left hand member of the previous inequality is infinite. So there is a contradiction and we can conclude that (ρ ε ) is bounded. Then it follows from (2.7) that F ε (ρ ε ) is bounded as well, whence the assertion (i) .

Let us now prove the assertion (ii); we know that the sequence (ρ ε ) is bounded and therefore admits weak* cluster points. Given such a cluster point, we can assume, without loss of generality, that ρ ε * ρ in M + (Ω). Let us introduce for every t > 0 the set

E t := x ∈ Ω : lim inf r→0 ρ(B(x, r)) ω d r d > t .
Thanks to the upper semi-continuity of the map x → ρ(B(x, r), we infer that E t is a Borel subset of Ω. We are going to prove that

(2.8) lim t→+∞ ρ(E t ) = 0.
To that aim, we consider the family of closed balls in R d defined by

F t := B(x, r) : x ∈ E t , r < r x , ρ(∂B(x, r)) = 0 ,
where r x > 0 is chosen so that ρ(B(x, r)) > tω d r d for every r < r x . Since F t determines a fine covering of the bounded Borel set E t , we may invoke the Vitali-Besicovitch covering theorem (see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]Thm 2.19]) which provides the existence a countable subfamily (B n ) such that

(2.9) ρ(B n ) > t L d (B n ) ∀n , ρ(E t \ ∪ n B n ) = 0.
Next we associate with the weak* convergent sequence (ρ ε ), two set functions defined on Borel subsets A ⊂ Ω:

η ε (A) := F ε (ρ ε A) , η(A) := lim inf ε→0 η ε (A).
We can readily check that η ε and η are monotone with respect to the inclusion while η ε (Ω) = F ε (ρ ε ) implies that η(Ω) ≤ β < +∞. Moreover, by the first inequality in (2.1), η ε is super-additive on disjoint Borel subsets. Obviously this holds true also for the set function η. By applying this property to the sequence of disjoint balls B n , we get the upperbound:

(2.10)

n η(B n ) ≤ C .
On the other hand, thanks to the coercivity inequality (2.6) that we apply with Ω = B n , we obtain:

η ε (B n ) = F ε (ρ ε B n ) ≥ -(δ) ρ ε (B n ) ρ ε (B n ) β ε (B n ) -1 + , (2.11) where β ε (B n ) := ε d m δ B n ε . Since ρ(∂B n ) = 0 by construction and
thanks to (2.3) and (2.9), we infer that

lim ε→0 ρ ε (B n ) β ε (B n ) = ρ(B n ) δ -d L d (B n ) ≥ t δ d .
Therefore, passing to the limit ε → 0 in (2.11) , we deduce that:

(2.12)

η(B n ) ≥ -(δ) ρ(B n ) (t δ d -1) + .
All in all, after collecting the second equality of (2.9), (2.10) and (2.12), we are led to:

ρ(E t ) ≤ n ρ(B n ) ≤ C -(δ)(tδ d -1) + .
Our claim (2.8) follows by sending t → +∞. The absolute continuity property ρ L d stated in the assertion (ii) is a consequence of the Besicovitch differentiation theorem [1, theorem 2.22]), which states that the singular part ρ s in the Lebesgue-Nikodym decomposition of ρ with respect to the Lebesgue measure coincides with ρ E ∞ being E ∞ = ∩ t>0 E t . In our case ρ s = ρ(E ∞ ) = 0 due to (2.8).

2.3.

Upper-bound of energies. In the same way as in [START_REF] Mathieu Lewin | Coulomb and Riesz gases: the known and the unknown[END_REF], we will be using an upper bound of ξ (S) when S is arranged on a d-dimensional periodic Bravais lattice G1 . To such a lattice we associate the -Epstein zeta function defined for every r > 0 by:

(2.13) Λ ,G (r) := x∈G\{0} (r|x|).
In the case where the cartesian lattice G = Z d is used, we will write simply Λ (r). The finiteness of this function for large r, under the condition (H3), turns out to be be crucial for deriving an uniform upper bound for the scaled energy F ε given in (1.8).

Lemma 2.7. Under (H1) -(H3), there exists

C G > 0 such that (2.14) Λ ,G (r) ≤ C G r d (r 0 ) r d 0 + d +∞ r 0 t d-1 (t) dt ∀r ≥ r 0 max{1, a -1 G } ,
where a G := min{|y| : y ∈ G \ {0}}.

Proof. To simplify, we chose the lattice G so that a G = 1. Up to substituting with + which satisfies + ≥ and + = on [r 0 , +∞), we may also assume that is non-increasing on R + . Accordingly, for any s ∈ [0, (0 + )), the set of values { > s} forms a non-empty interval [0, -1 (s)). The pseudo inverse -1 (s) is the supremum of all t ≥ 0 such that (t) > s, and it is a monotone non-increasing function on [0, +∞). Therefore, we have the following equivalence:

-1 (s) > t ⇐⇒ (t) > s.
By applying the layer cake formula to the counting measure on G, we get

Λ ,G (r) = ∞ 0 N r (s) ds,
where the integer function N r (s) := ({x ∈ G \ {0} : (r|x|) > s}) satifies N r (s) = 0 if r ≥ r 0 and s ≥ (r 0 ) (we assumed that a G = 1). On the other hand, for any periodic Bravais lattice G ⊂ R d , there exists a constant C G > 0 such that

(B r ∩ G) ≤ C G r d , ∀r > 0.
This implies the inequality:

N r (s) = ({x ∈ G \ {0} : r|x| < -1 (s)}) ≤ C G r d ( -1 (s)) d .
Therefore, for every r ≥ r 0 , we are led to:

Λ ,G (r) ≤ C G r d (r 0 ) 0 ( -1 (s)) d ds.
Then, after noticing that {s ∈ [0, (r 0 )] : -1 (s) > t} = [0, (r 0 ) ∧ (t)] holds for any t ≥ 0, we obtain the desired inequality by applying once again the layer cake formula:

(r 0 ) 0 ( -1 (s)) d ds = d ∞ 0 ( (r 0 ) ∧ (t)) t d-1 dt = (r 0 )r d 0 + d +∞ r 0 t d-1 (t)dt.
Remark 2.8. In view of Lemma 2.7, the -Epstein zeta function Λ ,G (r) vanishes at infinity. However the behavior in O(r -d ) as r → ∞ suggested by (2.14) is not optimal as we can see in the case of a Riesz potential (r) = r -s with s > d, where Λ ,G (r) = C r -s . On the other hand, it is noteworthy that Λ ,G is not continuous in general. A very simple example to see this is given by the step function = 1 2 1 1 [0,1) which satisfies (H1) -(H3). For d = 1 and the lattice G = Z, we find that Λ (r) = {n ∈ N : 0 < nr < 1} = [r -1 ] , where [•] denotes the integer part.

Next, by applying Lemma 2.7 in the case of the Cartesian lattice G = Z d , we derive a fundamental upper bound for the short-range interaction energy. Lemma 2.9. Let r > 0 and S be a finite subset of the lattice r Z d . Then

(2.15) ξ (S) ≤ (S) Λ (r)
As a consequence, for every a > 0, there exists

S ε ⊂ Ω such that ρ ε = ρ ε Sε satisfies (2.16) ρ ε * a L d Ω and lim sup ε F ε (ρ ε ) ≤ a Λ (a -1 d ) |Ω|
where the right hand side upper bound is finite whenever 0 ≤ a ≤ r -d 0 . Proof. Let S = {r x i : 1 ≤ i ≤ N } where N = (S) and x i ∈ Z d . Noticing that, for every i, the set {x i -x j : j = i} consists of N -1 distinct elements of Z d \ {0}, we infer that j =i

(r|x i -x j |) ≤ z∈Z d \{0} (r|z|) = Λ (r),
hence the desired inequality (2.15) by summing with respect to i.

Taking now ρ ε = ρ ε Sε where S ε = Ω ∩ (r ε Z) d and r ε = ε a -1 d
, we obtain a sequence such that ρ ε * a L d Ω as ε → 0. Indeed, by the periodicity of the Euclidean lattice, ρ ε converges to a uniform density on Ω while its total mass

ρ ε = ε d (S ε ) ∼ ε d (r -d ε |Ω|) converges to a |Ω| as ε → 0. Eventually, by applying (2.15), we get F ε (ρ ε ) = ε d ξ ( Sε ε ) ≤ ε d N ε Λ (a -1/d ), whence: lim sup ε F ε (ρ ε ) ≤ a Λ (a -1 d ) |Ω|.
The finiteness of Λ (a - (thereby fixing the volume of the so called fundamental domain of G).

The existence and the determination of an optimal lattice G for the following minimization problem:

inf{Λ ,G (r) : G = F Z d , |det(F )| = 1}
touches on a very hard and famous problem related to crystallisation conjectures (see [START_REF] Mathieu Lewin | Coulomb and Riesz gases: the known and the unknown[END_REF][START_REF] Blanc | The crystallization conjecture: a review[END_REF]). Note that, for a general cost , the answer to this problem will depend of the value r (thus of the local density a = dρ dx of the limiting measure ρ).

3. The Γ-convergence result.

3.1.

A quick overlook. The notion of Γ-convergence is popular in the community of calculus of variations and very much used in the analysis of sharp-interface models, dimension reduction for problems in mechanics, optimal design and homogenization. As pointed out in the introduction, this tool is also perfectly suited to justify a meanfield approach for large particle systems subject to a minimum energy criterion (see [START_REF] Serfaty | Systems of points with Coulomb interactions[END_REF]). For the convenience of the reader, let us give here some basic definitions and main properties. For futher details, we refer to the monographs [START_REF] Attouch | Variational convergence for functions and operators[END_REF][START_REF] Dal | An introduction to Γ-convergence[END_REF][START_REF] Braides | Gamma-convergence for Beginners[END_REF].

Let (E, τ ) be a metrizable topological space and consider a sequence of functionals F n : E → (-∞, +∞]. Then the lower Γ-limit F -and the upper Γ-limit F + of F n are defined by:

F -(u) := inf un→u lim inf n→∞ F n (u n ) , F + (u) := inf un→u lim sup n→∞ F n (u n ).
Both are τ -lower semicontinuous (see [START_REF] Attouch | Variational convergence for functions and operators[END_REF]), whereas in general it holds F -≤ F + . In practice it is useful to check that these functionals are proper i.e. that they range into R ∪ {+∞} being not identically +∞.

If F n admits a lower bound independent of n, this amounts to checking that the existence of u 0 ∈ E such that (3.1)

F + (u 0 ) = inf un→u 0 lim sup n→∞ F n (u n ) < +∞.
We say that

F n Γ-converges to F (denoted F n Γ - → F ) if F = F -= F +
or equivalently if the two following conditions are fulfilled: a) (lowerbound) For any sequence u n converging to u, we have the inequality

lim inf n→∞ F n (u n ) ≥ F (u);
b) (recovering sequence) For every u ∈ X such that F (u) < +∞, there exists (u n ) such that

u n → u and F n (u n ) → F (u).
This is the case in particular if F n = F does not depend of n; then F n Γ -→ cl(F ) where cl(F ) denotes the τ -lower semicontinuous envelope of F . Among all properties of Γ-convergence, we give some which will be used in this paper. Proposition 3.1. Let F n : E → (-∞, +∞] and assume (3.1). Then:

(i) F n Γ - → F ⇐⇒ cl(F n ) Γ - → F ; (ii) ( Kuratowski compactness Theorem) If (E, τ ) is a second
countable topological space (for instance a separable metric space), then any sequence (F n ) admits a Γ-convergent subsequence; (iii) ( convergence of infima) Suppose that F n Γ -→ F and that the following equi-coercivity property holds:

sup n F n (u n ) < +∞ =⇒ {u n }is τ -relatively compact.
Then lim n→∞ inf X F n = min X F and the minimum set for F coincides with the cluster points of all sequences (u n ) such that

F n (u n ) -inf F n → 0 ; (iv) ( stability) F n Γ - → F =⇒ F n +G Γ - → F +G for every continuous perturbation function G : E → R.
Remark 3.2. The continuity requirement for G in the assertion (iv) is often too restrictive. Actually the same conclusion holds under the following milder condition:

(3.2) inf K G > -∞ for any compact K ⊂ E G(u n ) → G(u) whenever u n → u and F (u) < +∞
For the convenience of the reader, a brief proof of the sufficiency of this condition is given below.

Proof. To check condition a), we consider a sequence (u n ) such that u n → u. Without loss of generality, assume that

F n (u n ) + G(u n ) ≤ C for a suitable constant C. Since G(u n ) is lower bounded (take K to be {u n , n ∈ N} ∪ {u}), we infer that F n (u n ) ≤ C for another constant C . By the Γ-convergence F n Γ - → F , it follows that F (u) ≤ lim inf n F n (u n ) < +∞. Therefore G(u n ) → G(u) and lim inf n (F n + G)(u n ) ≥ (F + G)(u).
For checking condition b), we may restrict to elements u ∈ E such that (F +G)(u) < +∞. Then F (u) < +∞ and any recovering sequence

u n → u such that F n (u n ) → F (u) will satisfy G(u n ) → G(u).

3.2.

The main result. In our context, the Γ-convergence issue applies to the sequence (F ε ) defined in (1.8) and to the ambiant topological space M + (Ω) embedded with the weak* topology (tight convergence). Thanks to the equi-coercivity property established in Proposition 2.5, there is no loss of generality in working in a fixed closed ball of M + (Ω) which is metrizable and compact. Therefore all the properties mentioned in the former subsection are applicable (after substituting the index n → ∞ with the continuous parameter ε → 0).

F (ρ) := Ω f (u) dx if ρ = u L d Ω +∞ otherwise
and f : R + → [0, +∞] is convex, l.s.c. and satisfies

(3.3) f (0) = f (0 + ) = 0 , lim inf t→+∞ f (t) t 2 > 0.
The proof of Theorem 3.3 is postponed to Section 5 while the property (3.3) is established in the next Sub-section. Let us now consider the Fenchel conjugate of f (implicitly extended by +∞ on (-∞, 0)) given by (3.4) f * (λ) := sup {λ t -f (t) : t ∈ R + } From (3.3), one can check that the supremum in (3.4) is actually a maximum (which is attained at t = 0 if λ ≤ 0). Therefore, f * is convex, continuous and vanishes on (-∞, 0]. As a consequence it admits left and right derivatives (f * ) (λ -) ≤ (f * ) (λ + ) for any λ so that the subdifferential ∂f * (λ) = [(f * ) (λ -), (f ) * ) (λ + )] is non empty.

Corollary 3.4. Let U ∈ C(Ω) be an external potential. Then, for every ε > 0, there exists a finite set S ε ⊂ Ω minimizing

I (ε) (U ) := inf S⊂Ω (x,y)∈S 2 \∆ |x -y| ε + x∈S U (x) .
Moreover sup ε ε d (S ε ) < +∞ and lim ε→0 ε d I (ε) (U ) = I (U ), where

(3.5) I (U ) := min u∈L 1 (Ω) Ω (f (u) + u U ) dx = - Ω f * (-U ) dx.
Furthermore any weak* cluster point of ρ ε Sε belongs to the minimum set of (3.5) given by

(3.6) S := u L d Ω : u(x) ∈ ∂f * (-U (x)) a.e. x ∈ Ω .
Remark 3.5 (Non-uniqueness). In general f is not strictly convex and ∂f * can be multi-valued (see for instance the hard spheres case depicted in Subsection 3.4 or the example of the step function given in (3.19) where f is piecewise affine). Note that, since f * vanishes on R -, it holds I (U ) = 0 for every non-negative potential U .

Remark 3.6. A natural variant of I (ε) (U ) consists in prescribing the total number N ε of particles to satisfy N ε ∼ κ ε -d for some given intensity factor κ ∈ (0, +∞). Accordingly, Corollary 3.4 can be restated by adding a total mass constraint in the limit problem, that is with I (U ) in (3.5) replaced by

I ,κ (U ) := inf Ω (f (u) + u U ) dx : u ∈ L 1 (Ω; R + ) , Ω u dx = κ .
The associated minimum S ,κ can be determined by selecting a suitable Lagrange multiplier depending implicitly on κ; as a consequence, an explicit form for S ,κ of the kind (3.6) is not available. Remark 3.7 (Clustering). For ε > 0 fixed, the lower semicontinuity property of F ε requires that (0+) = +∞. Otherwise, if (0 + ) < +∞, a sequence of subsets S n ⊂ Ω such that sup n F ε (ρ ε Sn ) < +∞ (thus retaining a finite number of points) can collapse into several clusters while retaining finite energy. In this case, the relaxed functional cl(F ε ) can be obtained directly by extending its domain to ε-empirical measures associated with multisets (instead of sets) and by extending the definition (1.6) accordingly taking into account the multiplicity of each cluster of particles. Note however that considering cl(F ε ) instead of F ε will not change the mean-field energy F given in Theorem 3.3 in virtue of the assertion i) of Proposition 3.1.

3.3.

Characterization and properties of f . The convex integrand f will be characterized indirectly through its Fenchel conjugate. For every λ ∈ R and any Borel subset B ⊂ R d , we define:

(3.7) Γ (λ, B) := sup {λ (S) -ξ (S) : S finite ⊂ B} .
The key properties of this bivariate function are summarized in the two following lemmas. Let ϕ : R → R be the convex, continuous function defined by:

(3.8) ϕ(t) =      (1+t) 2 4 if t ≥ 1 t if 0 ≤ t < 1 0 if t < 0 .
A straightforward computation shows that its Fenchel conjugate is given by:

ϕ * (ζ) = ζ(ζ -1) + if ζ ≥ 0 , ϕ * (ζ) = +∞ if ζ < 0.
Lemma 3.8. Assume that L d (B) < +∞. Then the map λ → Γ (λ, B) is convex, continuous, vanishes on R -and satisfies:

(3.9) 0 ≤ Γ (λ, B) ≤ m δ (B) -(δ) ϕ λ -(δ)
, for every λ ∈ R, being δ > 0 such that -(δ) ∈ (0, +∞] and ϕ being defined by (3.8).

Notice that the right hand side of (3.9) is convex as a function λ. On the other hand, due to the linear behavior of ϕ on [0, 1], we can easily infer that Γ (λ, B) ≤ λ m δ (B) whenever -(δ) = +∞.

Proof. As a supremum of affine functions, Γ (λ, B) is convex and l.s.c. with respect to λ; it is non-negative (follows by taking S the empty set) and vanishes if λ ≤ 0. The continuity property follows classically from (3.9) which provides the finiteness of Γ (λ, B) since m δ (B) < +∞ once L d (B) < +∞. It remains to prove the upper bound in (3.9) which clearly follows from (2.4). Indeed, in term of ζ = (S) m δ (B) , we have for every finite S ⊂ B:

λ (S) -ξ (S) ≤ m δ (B) (λ ζ --(δ) ζ (ζ -1) + ) .
Taking the supremum of the right hand member with respect to ζ ≥ 0 and noticing that ζ (ζ -1) + = ϕ * (ζ), we derive (3.9)after some easy manipulations. Lemma 3.9. For every λ ≥ 0, the set function B → Γ (λ, B) is subadditive on disjoints Borel subsets and translation invariant.

Proof. Let B 1 , B 2 be such that B 1 ∩ B 2 = ∅ and let S be a finite subset of B 1 ∪ B 2 . Then, setting S i = S ∩ B i , we get a partition S = S 1 ∪ S 2 and by appling the super-additivity part of (2.1), we infer that:

λ (S) -ξ (S) ≤ 2 i=1 (λ (S i ) -ξ (S i )) ≤ Γ (λ, B 1 ) + Γ (λ, B 2 ),
hence the desired sub-additivity property by taking the supremum with respect to S ⊂ B 1 ∪ B 2 . The invariance by translation is trivial.

In virtue of Lemma 3.9, we may now apply to Γ(λ, •) a classical result by Krengel [START_REF] Krengel | Ergodic theorems[END_REF], which ensures the existence, for every λ, of a limit for the ratio Γ(λ,Q k ) k d as k → +∞ (thermodynamical limit). Let us define the function

(3.10) g (λ) := inf k∈N Γ(λ, Q k ) k d .
Then we have:

(3.11) g (λ) = lim k→+∞ Γ(λ, Q k ) k d = lim ε→0 ε d Γ(λ, Q 1/ε ) .
For a proof of (3.11), we refer for instance to [START_REF] Licht | Global-local subadditive ergodic theorems and application to homogenization in elasticity[END_REF]. Since the effective profile f in Theorem 3.3 will be identified through the relation f * = g (see the last step of the proof in Section 5), we now establish some useful bounds for g . Recalling the definition of Λ in (2.13), we introduce the function defined on R by:

H (t) := t Λ (t -1 d ) if t > 0 , +∞ if t ≤ 0 .
From (H 1 )(H 2 ), it is straightforward that this function H is l.s.c. and finite on (0, r -d 0 ] while H (0 + ) = H (0 + ) = 0. Note that H is not convex in general, even it could be discontinuous, as happens for which is the step function ) . Indeed, in this case and for d = 1, we get H (t) = t [t] (see Remark 2.8). However H is convex continuous on R + in many classical cases including that of hyper singular Riesz potentials (r) = r -s where H (t) = C t 1+s/d . Proposition 3.10. The fonction g (λ) is convex, continuous, nonnegative, vanishes on R -and, for λ > 0, satisfies the inequalities:

(r) = 1 2 1 1 [0,1
(3.12) H * (λ) ≤ g (λ) ≤ δ -d -(δ) ϕ λ - (δ)
, being ϕ given by (3.8) and δ > 0 choosen such that -(δ) > 0.2

Proof. In view of (3.11) and of the convexity of Γ(•, Q k ), the function g λ is convex, non-negative and vanishes for λ ≤ 0 as a pointwise limit as k → ∞ of the sequence of functions k -d Γ(•, Q k ).

To prove the right hand inequality in (3.12), it is enough to apply (3.9) with B = Q k and, after dividing by k d , pass to the limit taking into account that, in virtue of (2.3)

, m δ (Q k ) ∼ ( k δ ) d as k → ∞.
Since the majorant is a convex, continuous function of λ, we infer that g is continuous on R as well.

Eventually, let us apply, for every t > 0 anf k ∈ N * , the fundamental upper bound (2.15) to the finite subset S t,k

:= t -1 d Z d ∩ Q k .
Then, recalling (3.7), we have:

Γ (λ, Q k ) ≥ λ (S t,k ) -ξ (S t,k ) ≥ (S t,k ) λ -Λ (t -1 d ) .
Since (S t,k ) ∼ t k d as k → ∞ and in virtue of (3.11), we are led to:

g (λ) = lim k→∞ Γ (λ, Q k ) k d ≥ t(λ -Λ (t -1 d )) = λt -H (t).
The left hand side inequality in (3.12) follows by taking the supremum with respect to t.

Corollary 3.11. The integrand f := g * is convex, l.s.c. and satisfies f (0) = 0 while f (t) = +∞ for t < 0. Moreover, for every t > 0 and δ > 0, we have the inequalities:

(3.13) -(δ) t (t δ d -1) + ≤ f (t) ≤ H * * (t). Accordingly, f is finite on [0, r -d 0 ]
, monotone non-decreasing on [0, +∞) and satisfies:

(3.14) f (0 + ) = 0 , lim inf t→+∞ f (t) t 2 ≥ sup δ>0 -(δ)δ d > 0
Proof. Passing carefully to Fenchel conjugates in inequalities (3.12), we are led to (3.13) from which the other statements follow directly.

In particular, as f (0 + ) = 0, we infer that the convex function f is monotone non-decreasing on R + .

Remark 3.12 (growth conditions). The inequality in (3.14) confirms that f grows at least quadratically at infinity, as announced in the introduction (see Theorem 3.3). More specifically, we can highlight two subcases for a cost satisfying1 (H1)

-(H 3 ). a) k := sup δ>0 -(δ)δ d = +∞. Then lim inf t→+∞ f (t)
t 2 = +∞ and f has a super quadratic growth. Note that this conclusion is consistent with the case (r) = r -s for s > d (see the next subsection).

b)

∞ 0 + (t)t d-1 dt < +∞. In this case k < +∞ and, thanks to (3.13) and to the estimate given in Lemma 2.7 (that we can apply to + with r 0 = 0), we obtain the lower and upper bounds:

0 < k ≤ lim inf t→+∞ f (t) t 2 ≤ lim sup t→+∞ f (t) t 2 ≤ C ∞ 0 + (t)t d-1 dt.
It follows that, under the integrability condition R d + (|x|) dx < +∞, f enjoys a quadratic growth from above and from below.

(r) = +∞ if r < 1 0 if r < 1
The computation of g through (3.10) and (3.11) leads to a linear function on R + namely g (λ) = γ d λ, where γ d denotes the densest spheres packing volume fraction in R d . This famous universal constant can be defined as

(3.15) γ d := inf k∈N * S(Q k ) k d = lim k→∞ S(Q k ) k d ,
where, for any Borel set A ⊂ R d , S(A) denotes the maximal number of points in A with mutual distance larger or equal to 1. The meanfield energy density f given by Theorem 3.3 is therefore the indicator function of the interval [0, γ d ]

f (t) = 0 if t ≤ γ d , f (t) = +∞ otherwise .
Furthermore, for every continuous external potential U ∈ C(Ω), we recover from Corollary 3.4 the convergence:

min

S∈Fε(Ω) ε d z∈S U (z) → γ d Ω U (z) dz ,
where F ε (Ω) is the family of finite subsets S ⊂ Ω satisfying |x -y| ≥ ε for all (x, y) ∈ S 2 \ ∆.

Remark 3.13. A variant of the previous result was obtained recently in [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF] in the case where the total number of particles N ε is prescribed to satisfy N ε ε d → κ as ε → 0 where κ > 0 is a given real parameter.

With our notations this condition amounts to restrict the Γ-limit F to measures ρ such that Ω udx = κ. Since the domain of F consists of density measures ρ = u L d Ω such that u ≤ γ d a.e., the latter integral condition requires that κ ≤ γ d |Ω| hence a congestion ratio θ := κ γ d |Ω| not larger than 1. In this case and if, following the classical empirical measure representation, u is normalized to be a probability density by setting ũ := u Ω udx = u κ , we recover a mean-field energy vanishing for ũ ≤ γ d κ and infinite otherwise, exactly as stated in [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF]Thm 6.1]. Note that the duality technique used there could only handle cost funtions taking values in the discrete set {0, +∞}.

3.4.2. The case of Riesz potentials. Short range potential of Riesz type corresponds to fixing s > d and taking

l(r) = r -s on R d + .
In this case, it's easy to establish from the homogeneity of the cost that the bivariate function Γ defined in (3.7) satisfies, for every t ≥ 0, the following scaling law :

(3.16) Γ (tλ, B) = tΓ (λ, t 1/s B) .
It follows from (3.16) that:

Γ (tλ, Q k ) k d = tΓ (λ, t 1/s Q k ) k d = t 1+d/s Γ (λ, Q t 1/s k ) (t 1/s k) d .
Sending k → ∞ and applying (3.11) two times, we get:

g (tλ) = lim k→∞ Γ (tλ, Q k ) k d = lim k→∞ t 1+d/s Γ (λ, Q t 1/s k ) (t 1/s k) d = t 1+d/s g (λ).
In virtue of the equality f = g * , we deduce that:

(3.17)

f (t) = C(s, d) t 1+s/d
where C(s, d) = f (1) is a universal constant. We thus recover the Γ-convergence result proved in [START_REF] Hardin | Asymptotic properties of short-range interaction functionals[END_REF][START_REF] Douglas | Large deviation principles for hypersingular Riesz gases[END_REF].

Remark 3.14. If we chose (r) = r -s where s < d, then condition (H3) is violated, and the scaling defined in equation (1.8) that we used to define F ε will result in an infinite Γ-limit. This means that F (ρ) will be equal to +∞ whenever ρ is not equal to zero, and F (0) will be equal to zero. We can observe this when we consider a system of N ε particles in S ε ⊂ Ω such that ρ ε := ρ ε Sε converges weakly to ρ, and sup ε F ε (ρ ε ) < +∞. Assuming ρ = 0, then we have N ε ∼ ρ ε -d as ε → 0. Moreover, due to the power law property of , we can write:

(3.18) F ε (ρ ε ) = ε s+d ξ (S ε ) ∼ ρ 2 ε d-s ξ (S ε ) N 2 ε .
As satisfies (1.3), the convergence result of the long range case holds (with h N = N -2 , see [START_REF] Bindini | Relaxed many-body optimal transport and related asymptotics[END_REF], [START_REF] Serfaty | Coulomb gases and Ginzburg-Landau vortices[END_REF]). Therefore, based on (1.4), and given that the standard empirical measure linked to S ε converges to ρ =

(3.19) (r) = M 2 if r < 1 0 if r < 1 (M positive parameter) .
Applying the lower bound (3.13) with δ = 1, we deduce that f ≥ h where h(t) := M 2 t(t -1) + .

In turn this lower bound is optimal for integer values of t since, as proved below, f coincides on R + with the piecewise affine interpolation of h given by:

(3.20) f (t) = h(k) + (t -k)(h(k + 1) -h(k)) ∀t ∈ [k, k + 1] , k ∈ N.
Proof. Owing to (3.11), the Fenchel conjugate of f is given by

g (λ) = lim K→+∞ 1 K sup S⊂[0,K] λ (S) - M 2 ({(x, y) ∈ S 2 \ ∆ 1 }) ,
where ∆ 1 := {(x, y) ∈ R d × R d : |x -y| ≥ 1}. Let S ⊂ [0, K] be an optimal set which we split in K disjoint pieces namely S = K i=1 S i where S i = S ∩ [i, i + 1). Let us denote n i the number of points in S i . By pushing them to the center of the interval [i, i + 1], we see that the number of pairs in S 2 \ ∆ 1 decreases to K i=1 n i (n i -1). It follows that:

g (λ) = lim K→∞ sup n i ∈N 1 K λ K i=1 n i - M 2 K i=1 n i (n i -1) = sup n∈N λ n - M 2 n(n -1) = (h + χ N ) * (λ),
where χ N denotes the indicator function of the integers. Therefore f = (g ) * is nothing else but the convexification of h + χ N given by the interpolation formula (3.20).

As demonstrated above, optimal point configurations for a constant external potential are obtained by periodically grouping a suitable number of points. Therefore, optimal sets are essentially multisets (see Remark 3.7). It is probable that a similar phenomenon occurs in higher dimensions.

The situation will vary if we consider a non-monotonic step function, such as the following:

(r) = 1 for r ∈ [0, 1], (r) = 4 for r ∈ (1, 2), (r) = 0 for r ≥ 2.

In this case, we anticipate that optimal configurations may be periodic, but associated with a non-uniform Voronoi tessellation, consisting of patterns of different sizes, as observed in the context of optimal location problem (see [START_REF] Bouchitté | Asymptotic analysis of a class of optimal location problems[END_REF]Sec 3.4]).

Proof of the main Theorem

First, we check the properness property (2.16) to make sure that the upper Γ-limit of F ε is not trivial. To do this, it is sufficient to apply Lemma 2.9 by choosing u 0 = a 1

1 Ω for a ∈ [0, r -d 0 ]. Next, by virtue of the equi-coercivity property of F ε proved in Proposition 2.5 and of the Kuratowski compactness theorem (see Proposition 3.1 and the introductory comment of Subsection 3.2), we can find a sequence ε k → 0 and a weak* lower semicontinuous functional

F : M + (Ω) → [0, +∞] such that F ε k Γ - → F as k → ∞.
Note that the limit F may a priori depend on the chosen sequence ε k → 0. Accordingly, we will complete the proof of Theorem3.3 in two steps which are outlined below:

Step 1: we show that F is a local functional of the form (4.1)

F (ρ) := Ω j(x, u) dx if ρ = u L d Ω +∞ otherwise
where j : Ω × R → [0, +∞] is a suitable convex normal integrand such that j(•, 0) = 0 a.e. in Ω. As a consequence F is convex, weak* l.s.c. and coincides with its Fenchel biconjugate, i.e. :

F (ρ) = F * * (ρ) = sup v∈C(Ω) v dρ -F * (v) .
Step 2: we identify the Fenchel conjugate F * in terms of the convex function g defined in(3.10), namely:

F * (v) = Ω g (v) dx , for every v ∈ C(Ω).
It follows that the limit functional F does not depend on the sequence (ε k ). Also, since (3.14) the convex function f = g * has a superlinear growth at infinity, by applying a classical result on convex functionals on measures (see for instance [START_REF] Bouchitté | Integral representation of convex functionals on a space of measures[END_REF]), we obtain the equalities F = F * * = F where:

F (ρ) := Ω f (u) dx if ρ = u L d Ω +∞ otherwise .
This will conclude the proof of the Γ-convergence of the whole sequence (F ε ) as stated in Theorem 3.3.

Proof of Step 1. Let ρ ∈ M + (Ω) be such that F (ρ) < +∞. Then, there exists a recovering sequence ρ k * ρ such that lim sup k→∞ F ε k (ρ k ) = F (ρ) < +∞. By the assertion ii) of Proposition 2.5, we infer that ρ is an absolutely continuous measure. Accordingly, there exists a functional J : L 1 (Ω) → [0, +∞] such that

F (ρ) := J(u) if ρ = u L d Ω +∞ otherwise .
The following result will be crucial for deriving the integral representation and the convexity of J. Its delicate proof is postponed to the end of this section.

Lemma 4.1. The functional J : L 1 (Ω) → [0, +∞] defined above satisfies the following: (i) J is weakly lower semicontinuous and satisfies J(0) = 0;

(ii) The domain of J is a subset of L 1 (Ω; R + ) and J(u 1 1 A ) ≤ J(u) holds for every u ∈ L 1 (Ω, R + ) and every Borel subset A ⊂ Ω; (iii) It holds J(u + v) = J(u) + J(v) whenever u v = 0 is satisfied almost everywhere in Ω.

In view of the assertions (i) and (iii) of Lemma 4.1 and since the Lebesgue measure on Ω is atomless, we may apply a classical integral representation (see for instance Hiai and Umegaki [START_REF] Hiai | Representation of additive functionals on vector-valued normed Köthe spaces[END_REF][START_REF] Hiai | Integrals, conditional expectations, and martingales of multivalued functions[END_REF] or the monograph [START_REF] Buttazzo | Semicontinuity, relaxation and integral representation in the calculus of variations[END_REF]) according to which there exists a suitable convex normal integrand j such that (4.1) holds. Moreover, as J ≥ 0 and J(0) = 0, we have j(•, 0) = 0 a.e. while, due to the assertion (ii), the integrand j satisfies j(x, t) = +∞ if t < 0.

Proof of Step 2. From Step 1 and by a classical result on integral functionals (see for instance [START_REF] Bouchitté | Integral representation of convex functionals on a space of measures[END_REF]), the Fenchel conjugate of F is given for every v ∈ C(Ω) by:

F * (v) = sup u∈L 1 (Ω) Ω v u dx - Ω j(x, u(x)) dx = Ω j * (x, v(x)) dx.
Obviously we may extend this equality to all functions v ∈ L ∞ (Ω). Noticing that j * (x, 0) = -inf j(x, •) = -j(x, 0) = 0, we observe that, for every λ ∈ R and for every hypercube Q(x 0 , a) ⊂ Ω, we have

F * (λ 1 1 Q(x 0 ,a) ) = Q(x 0 ,a) j * (x, λ) dx.
Next we claim that, for any such an hypercube Q(x 0 , a) ⊂ Ω, the following holds:

(4.2) F * (λ 1 1 Q(x 0 ,a) ) = a d g (λ)
Suppose that this claim is true. Then, by considering Lebesgue points of j * (•, λ) for λ in a dense countable subset D of R, we can find a Lebesgue negligible subset N ⊂ Ω such that j * (x, λ) = g (λ) for all (x, λ) ∈ (Ω \ N ) × D. Thanks to the continuity of g proved in Proposition 3.10 and to the convexity of j * (x, •), the latter equality can be then extended to all (x, λ) ∈ (Ω\N )×R, so that we have F * (v) = Ω g (v) dx for every v ∈ C(Ω). Hence the conclusion of Step 2 is reached and the proof of Theorem 3.3 is complete provided we can confirm (4.2).

We now focus on the proof of the equality (4.2). For λ ≤ 0, this equality is trivial since g (λ) = j * (x, λ) = 0. Next we observe that, for every λ ≥ 0 and Q(x 0 , a) ⊂ Ω, we have:

F * ε k (λ 1 1 Q(x 0 ,a) ) := sup λ ρ(Q(x 0 , a) -F ε k (ρ) : ρ ∈ M + (Ω) = ε d k sup S⊂Ω {λ (S ∩ Q(x 0 , a)) -ξ ε k (S)} = ε d k sup S⊂Q(x 0 ,a) {λ (S) -ξ ε k (S)} = ε d k sup S ⊂Q(x 0 , a ε k ) {λ (S ) -ξ (S /ε k )} = ε d k Γ (λ, Q(x 0 , a ε k ))
where:

-to pass from the second to the third line, we substitute any competitor S ⊂ Ω with S ∩ Q(x 0 , a) which has larger energy;

-to pass from the third line to the two last lines, we set S = S /ε k for going back from the ε k -scaled energy (1.7) to the ground interaction energy ξ and ultimately recover the set function Γ defined in (3.7). Therefore, thanks to (3.11), we can pass to the limit k → +∞ (the position of x 0 is irrelevant) and obtain the equality

lim k→+∞ F * ε k (λ 1 1 Q(x 0 ,a) ) = a d g (λ).
So, proving (4.2) reduces to checking the equality

lim k→+∞ F * ε k (λ 1 1 Q(x 0 ,a) ) = F * (λ 1 1 Q(x 0 ,a) ),
that we rewrite in the equivalent form: The left hand side infimum in (4.3) being non-positive (easily seen by taking ρ = 0 as a competitor), we may apply Proposition 2.5 with the choice U = -λ 1 1 Q(x 0 ,a) . Therefore any minimizing sequence (ρ k ) for the left hand side of (4.3) is bounded in M + (Ω) hence weakly* relatively compact. By the assertion iii) of Proposition 3.1, we will be able to conclude the convergence of infima in (4.3) if we can show that (4.4)

F ε k + G Γ - → F + G being G(ρ) := -λ ρ(Q(x 0 , a).
In virtue of Theorem 3.3, we already know that F ε k Γ -→ F . Then it is enough to invoke the stability property of the assertion iv) of Proposition 3.1. However the functional G given above is not weak* continuous on M(Ω) and therefore, we need to verify the less stringent requirements set out in (3.2). The first one is satisfied since |G(ρ)| ≤ λ ρ . For the second one, we observe that, if F (ρ) < +∞, then ρ is of the form ρ = u L d Ω (see Proposition 2.5), hence ρ(∂Q(x 0 , a)) = 0 and every sequence ρ n * ρ satisfies ρ n (Q(x 0 , a)) → ρ(Q(x 0 , a)) = Q(x 0 ,a) u dx.

This confirms the validity of (4.4), hence that of (4.3). As a result the equality (4.2) is proved and, as announced, this achieves Step 2 and the proof of Theorem 3. To show that J(0) = 0, we consider as S k a singleton {x 0 } so that the associated measure ρ k = (ε k ) d δ x 0 satisfies ρ k * 0 while F ε k (ρ k ) = 0.

(ii) By the definition of the Γ-limit, F (ρ) < +∞ implies that ρ is a weak* limit of a sequence (ρ k ) in M + (Ω), hence of the form ρ = u L d Ω with u ≥ 0. Let now u ∈ L 1 (Ω; R + ) and A a Borel subset of Ω. We show first that J(u 1 1 A ) ≤ J(u) if L d (∂A) = 0. We may assume that J(u) < +∞. Hence, there exists a family of subsets S k ⊂ Ω such that ρ k = ρ indeed the convergence ρ k * ρ is tight while ρ(∂A) = 0. Therefore, F ε k (ρ k ) ≤ F ε k (ρ k ). By passing to the limit, as k → +∞, we deduce that

J(u 1 1 A ) ≤ F (ρ 1 1 A ) ≤ lim inf k→∞ F ε k (ρ k ) ≤ lim sup k→∞ F ε k (ρ k ) = J(u).
To extend the inequality to any Borel subset A of Ω, it is enough to consider an approximating sequence (A n ) such that L d (∂A n ) = 0 , L d (A n ∆A) → 0, and then pass to the limit in the inequality J(u) ≥ J(u1 1 An ) while letting n → +∞. Indeed, the conclusion will then follow from the lower semicontinuity of J with respect to the norm convergence in L 1 (Ω). Now, to construct such a sequence (A n ), we consider a compact subset K n ⊂ A and an open subset ω n ⊃ A such that L d (ω n \ K n ) ≤ 1 n . For every n, we can choose a suitable r n > 0 such that the enlarged open set A n = K n + B(0, r n ) satisfies L d (∂A n ) = 03 while A n ⊂ ω n . Then clearly L d (A n ∆A) ≤ L d (ω n \ K n ) → 0.

Let us now prove now the assertion (iii). In a first step, we assume that spt(u) ∩ spt(v) = ∅ so that there exists open subsets A ⊃ spt(u) and B ⊃ spt(v) such that dist(A, B) := η > 0.

We begin by proving the inequality J(u + v) ≥ J(u) + J(v). Without loss of generality, we may assume that J(u + v) < +∞ (hence u and v are non-negative). Then there exists a sequence of sets S k ⊂ Ω such that

ρ k = µ ε k S k * (u + v) L d Ω , F ε k (ρ k ) → J(u + v).
We write S k = S k ∪ S k where S k = S k ∩ A and S k = S k ∩ B are disjoint. Then, we have

ρ k = ρ k + ρ k where ρ k = µ ε k S k and ρ k = µ ε k S k . Clearly ρ k * u L d Ω while ρ k * v L d Ω.
Therefore, by applying the Γ -lim inf inequality to ρ k and ρ k while taking into account the super-additivity property of F ε k (see (2.1)), we deduce that:

J(u + v) = lim k→∞ F ε k (ρ k ) ≥ lim inf k→∞ F ε k (ρ k ) + lim inf k→∞ F ε k (ρ k ) ≥ F (u L d Ω) + F (v L d Ω) = J(u) + J(v).
To show the converse inequality J(u + v) ≤ J(u) + J(v), we assume without any loss of generality that J(u) < +∞ and J(v) < +∞. Then

Theorem 3 . 3 .

 33 Let satisfy the standing assumptions -(H3) and let F ε : M + (Ω) → [0, +∞] be given by (1.8). Then F ε Γ -→ F (for the weak* topology) where

3. 4 .. 3 . 4 . 1 .

 4341 ExamplesThe hard spheres model. The hard spheres potential is given by

2 ε≥ 3 . 4 . 3 .

 2343 D (ρ) > 0.This contradicts(3.18) since sup ε F ε (ρ ε ) < ∞. Therefore ρ = 0. The case of finite costs. Many examples of finite costs can be considered as, for instance, being a step function with compact support. Owing to Corollary 3.11, the effective convex integrand f has a quadratic growth on R + . The simplest one is the penalized version of the hard spheres potential defined by:

( 4

 4 .3) inf ρ∈M + (Ω) {F ε k (ρ) -λ ρ(Q(x 0 , a))} → inf ρ {F (ρ) -λ ρ(Q(x 0 , a))} .

3 .

 3 Proof of Lemma 4.1.(i) Let u k → u in L 1 (Ω). Then ρ k = u k L d Ω and ρ = u L d Ω are such that ρ k * ρ in M + (Ω). Since F = Γ-lim F ε k is wzak* lower semicontinuous, we infer that lim inf k→∞ J(u k ) = lim inf k→∞ F (ρ k ) ≥ F (ρ) = J(u).

, we have ρ k * u 1 1 A

 1 ε k S k satisfies ρ k * ρ := u L d Ω and F ε k (ρ k ) → J(u). If we let S k = S k ∩ A and ρ k = ρ ε k S k L d Ω, since

  stands for the diagonal of R d ;

	-(S) denotes the counting measure of a subset S ⊂ R d (+∞ if
	S is infinite);
	-L d is the Lebesgue measure in R d ; given any Borel set B, |B|
	is a short notation for L d (B) ; ω d is such that |B(x, r)| = ω d r d ;
	-C(Ω) denotes the Banach space of continuous functions on the
	compact subset Ω equipped with the uniform norm;
	-M(Ω) stands for the space of signed Radon measures on Ω
	equipped with the total variation norm;
	-P -(Ω) (resp. P(Ω)) is the subset of Borel measures µ ∈ M + (Ω)
	such that µ := µ(Ω) ≤ 1 (resp. µ = 1).
	-The topological support of µ ∈ M + (Ω) is denoted supp(µ)
	while µ A represents its restriction to a Borel subset A ⊂ Ω;
	-The bracket •, • will denote the duality between C(Ω) and
	M(Ω):
		v, µ = vdµ,
	This duality induces the weak* topology on M(Ω) which can
	be identified with the dual of C(Ω); as Ω is compact, the weak*
	convergence µ h	* µ in M + (Ω) implies the tight convergence
	since

  Remark 2.10. The upper bound (2.15) and (2.16) obtained by choosing G = Z d as the reference lattice are in general not optimal since the Epstein-zeta function of another lattice could provide better ones.

	1 d ) for a ∈ [0, r -d 0 ] follows from Lemma 2.7.

Obviously the inequality (2.15) holds true after replacing Λ with the Epstein function Λ ,G of any Bravais lattice G = F Z d , while for the valisity of (2.15), we need to add the normalization condition |det(F )| = 1

i.e. of the form G = F Z d for some invertible matrix F ∈ R d×d

The second inequality becomes g (λ) ≤ δ -d λ in the case where l -(δ) = +∞.

Here we use the fact that the functionα n (r) = L d ({x ∈ Ω : dist(x, K n ) > r})is bounded monotone non increasing so that it is continuous except possibly on a finite or countable subset of R + .

we consider recovering sequences

. Up to dropping the elements of S k which ar not in A and the elements of S k which are not in B, we may assume that spt(ρ k ) ⊂ A and spt(ρ k ) ⊂ B. Indeed, removing these points will not affect the weak* convergence to u L d Ω and v L d Ω respectively while the total energy ξ ε k will not increase. Therefore, by exploiting the right hand inequality in (2.1) and since ρ k = ρ k + ρ k converge weakly* to ρ = (u + v) L d Ω, we are led to the following set of inequalities:

where, in the second line, we used the tight convergence of ρ k , ρ k . The conclusion follows by noticing that + coincides with for large values where it is non-increasing. Thus the integrability condition (H3) implies that r d (r) → 0 as r → +∞. The desired sub-additivity inequality follows.

In a second step, we remove the strict separation condition on the supports of u and v by simply assuming that the upper-level sets A := {u > 0} and B := {v > 0} satisfy L d (A ∩ B) = 0 (which, for u, v nonnegative, is equivalent to say that u v = 0 a.e.). To that aim, possibly after substituting A, B with non-intersecting Borel representatives, we consider increasing sequences of compact subsets K n ⊂ A, K n ⊂ B such that:

In virtue of assertion (ii), we infer that, for every n:

In virtue of the lower semicontinuity of J, since u n → u and v n → v in L 1 (Ω), we deduce from above the inequalities J(u + v) ≤ J(u) + J(v) and J(u) + J(v) ≤ J(u + v), hence the desired additivity property.

Imath, Université de Toulon, BP 20132, 83957 La Garde Cedex-FRANCE Email address: bouchitte@univ-tln.fr Depto. de Matemática, FCFM, Univ. de Concepción, Concepción, Chile Email address, Mahadevan: rmahadevan@udec.cl


Nicolas Basset
email: bassetni@univ-grenoble-alpes.fr

Alexandre Donzé

Wordgen : a Timed word Generation Tool

Sampling timed words out of a timed language described as a timed automaton may seem a simple task: start from the initial state, choose a transition and a delay and repeat until an accepting state is reached. Unfortunately, simple approach based on local, on-thefly rules produces timed words from distributions that are biased in some unpredictable ways. For this reason, approaches have been developed to guarantee that the sampling follows a more desirable distribution defined over the timed language and not over the automaton. One such distribution is the maximal entropy distribution, whose implementation requires several non-trivial computational steps. In this paper, we present Wordgen which combines those different necessary steps into a lightweight standalone tool. The resulting timed words can be mapped to signals used for modelbased testing and falsification of cyber-physical systems thanks to a simple interface with the Breach tool.

INTRODUCTION

Timed automata (TA) is a well-known formalism to describe realtime systems. They are used in many verification tools such as Uppaal [START_REF] David | Uppaal SMC tutorial[END_REF][START_REF] Kim G Larsen | UPPAAL in a nutshell[END_REF] or Prism [START_REF] Kwiatkowska | PRISM: Probabilistic symbolic model checker[END_REF] since the algorithms behind their analysis are well understood. Verification can be conducted either exactly, using regions or zone graph abstractions, or approximately using statistical model checking (SMC) when the size of the system is too large. SMC can also deal with more complex systems, like hybrid systems and/or probabilistic systems.

Statistical model checking relies on the sampling of large numbers of trajectories of the system under test, which can be difficult. In [START_REF] Bohlender | A review of statistical model checking pitfalls on real-time stochastic models[END_REF], many examples are shown where the sampling of TA using current tools (e.g., UPPAAL SMC [START_REF] David | Uppaal SMC tutorial[END_REF] or Modes [START_REF] Carlos E Budde | A statistical model checker for nondeterminism and rare events[END_REF]) is not satisfactory for various reasons. The main observation of the authors is that these tools sample TA differently leading to different and often contradicting results. Indeed when sampling a TA, discrete and continuous choices must be made in each state. The order of these choices lead to different algorithms with different results. Moreover the distribution used for both the discrete and the continuous choice should not be arbitrary (see section 2.2 or [START_REF] Bohlender | A review of statistical model checking pitfalls on real-time stochastic models[END_REF] for more details).

To overcome these shortfalls, probability distribution over delay and transitions must be more formally specified. Distributions that maximize the entropy of the sampling process are a good choice because they are known to be the least biased [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF]. They are uniform distributions (UD) over timed languages which have the following characteristics:

(1) UD being defined over the timed language and not over the automaton, two automata recognizing the same timed language will have the same UD; (2) UD guarantee that given two timed words of the same length they are chosen with the same probability; (3) UD can be computed using time-language volume [START_REF] Asarin | Entropy of regular timed languages[END_REF] for which tractable algorithms have been described in [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF]; (4) However, UD computation relies on an analysis of the automaton which can be costly, in particular it can be significantly more costly than more "naïve" statistical model checking methods.

To our knowledge the only implementation of the computation of the uniform distribution over timed language is the tool-chain described in [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF], which is not satisfactory for several reasons which we detail in Section 3.4. The computation of UD requires the computation of the zone graph and thus is more involved than most classical computations on time automata (reachablility, emptiness, TCTL model checking). As a consequence there is no point to use UD sampling to solve these problems. UD sampling can, however, be used to construct probabilistic semi-algorithms to undecidable problems like the inclusion of time language [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF].

Another application of TA sampling is the validation of Cyber Physical System (CPS). In this setting we are interested in hybrid systems containing both a discrete part (e.g., modeled by finite state automata) and a continuous part (e.g., modeled by a system of ordinary differential equations) which interact with their environment through real-valued time signals. The formal verification of such systems is known to be difficult. A commonly used alternative is, similarly to SMC, to rely on generating time signals and check that the system behaves according to its specification over these signals.

It is then crucial to sample signals that cover the space of admissible input signals well. This space is not always easily described, and by extension sampled, using standard numerical constraints or even advanced specification languages such as signal temporal logics (STL). In [START_REF] Barbot | Generation of Signals Under Temporal Constraints for CPS Testing[END_REF][START_REF] Barbot | Falsification of Cyber-Physical Systems with Constrained Signal Spaces[END_REF], the authors show that TA can provide concise descriptions of the spaces of input signals for complex CPS models, and validate them using UD sampling. The cost of computing and sampling UD for a TA in these cases is negligible compared to other methods to generate admissible input signals (e.g., using STL monitoring-based rejection methods) to run and simulate the CPSs, making the method a good fit.

In this paper we presents the tool Wordgen which takes as input a TA, compute the UD over the timed language of the TA and outputs timed words sampled from this language using UD. We also present a simple interface for timed word based signal generation with the tool Breach [START_REF] Donzé | Breach, A Toolbox for Verification and Parameter Synthesis of Hybrid Systems[END_REF], which is a well-known tool for CPS validation. The paper is structured a follows: Section 2 describes the necessary definition, Section 3 describes the computational steps of Wordgen, Section 4 presents some implementation details, Section 5 presents a case study, and Section 6 concludes.

DEFINITIONS 2.1 Timed Automata

A timed automaton (TA) A is defined as a tuple (Σ, 𝑋, 𝑄, 𝑖 0 , F , Δ) where Σ is a finite set of events; 𝑋 is a finite set of clocks; 𝑄 is a finite set of locations; 𝑖 0 is the initial location; F ⊆ 𝑄 is a set of final locations and Δ is a finite set of transitions. The finite set of clock 𝑋 is a finite set of non-negative real-valued variables . A guard is a finite conjunction of clock constraints of the form 𝑥 ⊲⊳ 𝑐 with 𝑥 ∈ 𝑋 , 𝑐 ∈ N + and ⊲⊳∈ {<, ≤, =, ≥, >} . For a clock vector ì 𝑥 ∈ [0, ∞[ 𝑋 and a non-negative real 𝑡, we denote by ì 𝑥 + 𝑡 the vector ì 𝑥 + (𝑡, . . . , 𝑡). A transition 𝛿 ∈ Δ has an origin 𝛿 -∈ 𝑄, a destination 𝛿 + ∈ 𝑄, a label 𝑎 𝛿 ∈ Σ, a guard 𝔤 𝛿 and a reset function 𝔯 𝛿 determined by a subset of clocks 𝐵 ⊆ 𝑋 . Fig. 1(left) depicts an automaton with two clocks 𝑥 and 𝑦, a single location 𝑙 1 and two transitions. The transition labelled by 𝑎 is guarded by 𝑥 < 10 ∧ 𝑦 < 3 and reset 𝑦 while transition labelled by 𝑏 is guarded by 8 ≤ 𝑥 ≤ 10 ∧ 1 ≤ 𝑦 ≤ 2.

Additional clock constraints can be added to each location, these are known as invariants. In our setting, an automaton without invariant is obtained by adding the constraint of the invariant to all outgoing guard of a location.

A configuration 𝑠 = (𝑞, ì 𝑥) ∈ 𝑄 × [0, ∞) 𝑋 is a pair of location and a clock vector. The initial configuration of A is (𝑖 0 , ì 0). A timed transition is a pair (𝑡, 𝛿) of a time delay 𝑡 ∈ [0, ∞) followed by a discrete transition 𝛿 ∈ Δ. The delay 𝑡 represents the time before firing the transition 𝛿.

A run is an alternating sequence (𝑞 0 , ì 𝑥 0 )

𝑡 1 ,𝛿 1 ----→ (𝑞 1 , ì 𝑥 1 ) . . . 𝑡 𝑛 ,𝛿 𝑛 ----→ (𝑞 𝑛 , ì
𝑥 𝑛 ) of configurations where 𝑞 𝑖 is the successor of 𝑞 𝑖 -1 by 𝛿 𝑖 , and the vector ì 𝑥 𝑖 -1 + 𝑡 satisfies the guard 𝔤 𝛿 and ì 𝑥 𝑖 = 𝔯 𝛿 ( ì 𝑥 𝑖 -1 + 𝑡). This run is labelled by the timed word (𝑡 1 , 𝑎 1 ) • • • (𝑡 𝑛 , 𝑎 𝑛 ) where for every 𝑖 ≤ 𝑛, 𝑎 𝑖 is the label of 𝛿 𝑖 . The set of timed words labelling all the runs leading from the initial configuration (𝑖 0 , ì 0) to a final configuration (𝑞 𝑛 ∈ F ) is called the timed language of A.

Isotropic Sampling

Isotropic sampling is the sampling of timed words implemented in tools like UPPAAL-SMC or Modes with some variation. Given an automaton (Σ, 𝑋, 𝑄, 𝑖 0 , F , Δ) and a configuration of the automaton (𝑞, ì 𝑥) the next delay and next event is computed by choosing:

• The next transitions 𝛿 ∈ Δ such that 𝛿 -= 𝑞 and ∃𝑡 ≥ 0 s.t. ì 𝑥 + 𝑡 |= 𝔤 𝛿 , i.e., transition for which there exists a time where the guard is satisfied. All such transition are given the same probability to be chosen (discrete uniform law); • The next delay 𝑡 such that ∃𝛿 ∈ Δ s.t. 𝛿 -= 𝑞 ∧ ì 𝑥 +𝑡 |= 𝔤 𝛿 , i.e., delay for which there exists at least one acceptable transition at this time. Time 𝑡 is chosen uniformly (continuous uniform law) among acceptable times if this set is bounded and using an exponential distribution otherwise.

The order of these two choices will substantially change the sampling algorithm and its performance. Thus isotropic sampling is not well defined and as a consequence, different tools implement different strategies. As an example, if transitions are chosen first, an analysis of reachable guards must be conducted to select a satisfiable one, but then the time has to satisfy a single guard which is usually a convex set and thus is easy to sample. The other way around requires a similar trade-off. In [START_REF] Bohlender | A review of statistical model checking pitfalls on real-time stochastic models[END_REF] the authors show that these choices (and other choices of implementation) lead to tools giving contradictory results on some examples.

Another problem with isotropic sampling is that it is local to each step. More precisely along a timed word, a choice of time at step 𝑖 can make a guard at step 𝑗 > 𝑖 unsatisfiable. For example, in Fig. 1 the bottom plot shows isotropic sampling of the automaton on the left. Most transitions are sampled in the bottom right corner because the local rule cannot stir the simulation toward the green transition to avoid becoming stuck in a time lock.

COMPUTATIONAL STEPS

Given a timed automaton, Wordgen computes the UD over its time language and samples it. In this section, we explain the different operations performed by Wordgen to achieve this. Fig. 2 depicts an overview of the tool architecture and Fig. 3 shows a running example illustrating each computation step. The computation starts by parsing the input file which can be either in the Prism PTA or Uppaal file format. The result is a data structure containing a TA A. Fig. 3(a) shows an example of TA.

Forward Reachability and Splitting

First the automaton is transformed into a standard zone graph [START_REF] Bengtsson | Timed automata: Semantics, algorithms and tools[END_REF] by performing a forward reachability analysis as the traversal of a graph. Zones are sets of clock constraints implemented as Difference Bound Matrices (DBM) for which intersection time elapse (TE) and reset (𝑅𝑒𝑠𝑒𝑡) are efficient. The vertices of the graph are pairs of automaton location and zone, edges are automaton transitions with guard specified as zone. Starting from the initial location and the zone containing the clock vector ì 0, a graph is constructed by computing for each vertex of the graph (𝑙, 𝑧) and each transition 𝛿 s.t. 𝛿 -= 𝑙 the new vertex (𝛿 + , 𝑅𝑒𝑠𝑒𝑡 (𝑇 𝐸 (𝑧) ∩ 𝔤 𝛿 , 𝔯 𝛿 ). Fig. 3(b) shows the reachability graph of automaton Fig. 3(a) as computed by Wordgen.

Then, the zone graph is split so that for any vertex of the zone graph, for each available transition 𝛿, there exists two functions lb 𝛿 and ub 𝛿 of the clock vector ì 𝑥 such that, for all

𝑡 ∈ R + , ì 𝑥 + 𝑡 |= 𝔤 𝛿 is equivalent to 𝑡 ∈ [lb( ì 𝑥), ub( ì 𝑥)]
with lb and ub of the form ì 𝑥 ↦ → 𝑐 -𝑥 where 𝑐 ∈ N and 𝑥 is a clock or zero. This is always possible by splitting a zone in smaller one. Detail on the computation and proof of correctness can be found in [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF]. In Fig. 3(b) only one transition does not satisfy the splitting constraint (in black), Fig. 3(c) shows the result of the splitting algorithm, the vertex (1) have been split. In the worst case zone are spitted until only region remains which leads to an exponential blow up in the number of clocks.

𝑙 1 𝑎, {𝑦}, 𝑥 < 10 ∧𝑦 < 3 𝑏, {𝑥 }, 𝑥 ∈ [8; 10] ∧𝑦 ∈ [1; 2]

Weight and Distribution Computation

Each vertex (respectively, each edge) of the zone graph is equipped with a weight function 𝑣 𝑚 (𝑙, ì 𝑥) (respectively, 𝑣 𝑚 (𝛿, ì 𝑥)) which maps a clock vector to the volume of the timed language reachable from this clock vector in 𝑚 steps. The weight is computed with the following recursive equations following the theory of volume of timed language of [START_REF] Asarin | Entropy of regular timed languages[END_REF]:

𝑣 0 (𝑙, ì 𝑥) = 1 𝑙 ∈ F ; 𝑣 𝑚 (𝛿, ì 𝑥) = ∫ 𝑡 =ub 𝛿 𝑡 =lb 𝛿 𝑣 𝑚-1 𝛿 + , 𝔯 𝛿 ( ì 𝑥 + 𝑡) d𝑡; 𝑣 𝑚 (𝑙, ì 𝑥) = ∑ 𝛿 -=𝑙 𝑣 𝑚 (𝛿, ì 𝑥). (1)
Since lb 𝛿 and ub 𝛿 are linear functions of clock vectors, 𝑣 𝑚 (_) are polynomials with the clocks set as set of variable. The degree of the polynomials are 𝑚 and can be computed in polynomial time with respect to the graph size but exponential in 𝑚 and the number of clocks. Figure 3(d) shows the weights computed by Wordgen for the transition from 1 to 2 labeled by 𝑏.

The zone graph is finally augmented with a probability distribution to become a stochastic process. Discrete choices between transitions are computed using the weights of the transitions. Let us consider a vertex (𝑙, ì 𝑥) where 𝑘 transitions are available: 𝛿 1 , 𝛿 2 , . . . 𝛿 𝑘 . The probability to sample 𝛿 𝑖 is

𝑣 𝑚 (𝛿 𝑖 , ì 𝑥 ) 𝑘 𝑗 =1 𝑣 𝑚 (𝛿 𝑗 , ì 𝑥 )
. Continuous time distributions are computed as the derivative of the weight normalized by the total weights, that is from vertex (𝑙, ì 𝑥) for transition 𝛿 the probability density function (PDF) of the time 𝑇 is

𝑇 ↦ → ∫ 𝑡 =𝑚𝑖𝑛 (𝑇 ,ub 𝛿 ) 𝑡 =lb 𝛿 𝑣 𝑚-1 𝛿 + , 𝔯 𝛿 ( ì 𝑥 + 𝑡) d𝑡 1 𝑣 𝑚 (𝛿, ì
𝑥 ) . Note that the variable 𝑡 appears only in the numerator thus the PDF is a polynomial in 𝑡. The Cumulative Density function (CDF) is the integral of the PDF from 0 to 𝑡. It is denoted for transition 𝛿 and clock vector ì 𝑥 as (𝐶𝐷𝐹 𝑚 (𝛿, ì 𝑥)(𝑡)). It is also a polynomial in 𝑡. In Fig. 3(d) the volume and CDF of the transitions between 1 and 2 are reported.

Sampling

The last step consists of sampling the stochastic process obtained from the previous operation. To generate a timed word of length 𝑛, one starts with a sequence (𝑢 𝑖 ) 2𝑛

𝑖=1 ∈ [0, 1] of real values, which corresponds to a point in the unit box of dimension 2𝑛. These values can either be computed by Wordgen using pseudo-random number generators or low discrepancy sequences, they can also be computed by an external optimization algorithm like in [START_REF] Barbot | Falsification of Cyber-Physical Systems with Constrained Signal Spaces[END_REF] where Breach [START_REF] Donzé | Breach, A Toolbox for Verification and Parameter Synthesis of Hybrid Systems[END_REF] was used to drive the generation. The parameter 𝑚 of the distribution should be equal to 𝑛 -𝑖 at step 𝑖 for exact uniform sampling, however, as the computation of distribution for large value of 𝑚 became intractable, receding horizon can be used, a maximal value 𝑚 0 is computed then the distribution with 𝑚 = min(𝑚 0 , 𝑛 -𝑖) is used at each step. The correction of this approach and rate of convergence are discussed in [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF].

𝑙 1 [𝑎], (𝑥 < 2) ∧ (𝑦 < 4), {𝑦 } [𝑏 ], (𝑥 < 3) ∧ (𝑦 < 2), {𝑥 } (2) y in [0,2); x=0 b, [0 ; 2-y] (1) y=0; x in [0,2) a, [0 ; 2-x] b, y in [0,2); x < 3 a, [0 ; 2-x] (0) y=0; x=0 b, [0 ; 2-x] a, [0 ; 2-x] (2) y in [0,2); x=0 b, [0 ; 2-y] (3) y=0; x in [0,1) a, [0 ; 1-x] (1) y=0; x in (1,2) a, [1-x ; 2-x] b, [0 ; 2-y] a, [0 ; 1-x] a, [1-x ; 2-x] b, [0 ; 3-x] a, [0 ; 2-x] (0) y=0; x=0 b, [0 ; 2-x] a, [0 ; 1-x] a, [1-x ; 2-x] (a) (b) (c) (𝑑) 𝑚 1 2 3 4 𝑣 𝑚 (1 → 2, (𝑥, 𝑦)) 3 -𝑥 15-2𝑥 -𝑥 2 2 126-24𝑥 -3𝑥 2 -𝑥 3 6 1389-256𝑥 -48𝑥 2 -4𝑥 3 -𝑥 4 24 𝐶𝐷𝐹 𝑚 (1 → 2, (𝑥, 𝑦))(𝑡) 𝑡 3-𝑥 2𝑡 -𝑡 2 15-2𝑥 -𝑥 2 69𝑡 -12𝑡 2 +𝑡 3 126-24𝑥 -3𝑥 2 -𝑥 3 760𝑡 -138𝑡 2 +16𝑡 3 -𝑡 4 1389-256𝑥 -48𝑥 2 -4𝑥 3 -𝑥 4 (e)
Starting from the initial configuration of the automaton, the timed word is computed iteratively. At step 𝑖, in configuration (𝑙, ì 𝑥) the next transition 𝛿 is chosen by sampling the discrete probability distribution 𝑣 𝑚 (𝛿, ì 𝑥) using real 𝑢 2𝑖 . Then delay 𝑡 is distributed according to the PDF of 𝛿. This is done using the so-called inverse transform sampling method. This consists of computing the root of the CDF minus 𝑢 2𝑖+1 . Such root exists as the CDF is a strictly increasing polynomial taking values in [0, 1]. This root can be effectively computed using Newton's method. When a timed word of desired length is produced, it is outputted and the next timed word is computed. Fig. 3(e) shown an example of a timed word. We could have sampled first the time and then the transitions; the sampling distribution would be the same (see [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF]) but requires to sample from a piece-wise polynomial distribution.

Motivation for a New Implementation

In [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF][START_REF] Barbot | Generation of Signals Under Temporal Constraints for CPS Testing[END_REF] a tool-chain was implemented to conduct the same computation as Wordgen. The different parts of the computation were performed by specialized tools: the TA analysis was performed by a modified version of the model checker Prism, the probability distribution was computed by a script written in SageMath, a computer algebra system well suited for computation over polynomials, and the sampling was performed by the statistical model checker Cosmos. This work is a new standalone implementation in the OCaml programming language of all these steps. The main focus was to make this tool easier to use and to extend than the previous implementation. We describe some of the main advantages next.

Performance. While for their original purpose, each tool used in the legacy tool-chain are known to be very efficient, the specific operations they were used for in this work are not the ones they were designed and optimized for. Hence, a specialized re-implementation of these operations does not harm performance.

• Prism forward reachability algorithm is faster than our implementation over large automata. However, automaton size is not the bottle neck here: on all examples we tried, this time was small for both tools as compared to the other computational steps. • SageMath has an efficient library to handle polynomials formally and can handle various operations but only a few of these operations are relevant to our setting and those can be implemented with more specialized algorithms and data structures. In Table 1 we show on an example that Wordgen is faster, even when using exact rational numbers. 1: Performance of Wordgen compared to previous tools for distribution computation and sampling. All times are in seconds. All computations have been done on a single core of an Intel i7-6700 processor. Sampling generates 50000 trajectories of length 100, which is the number of trajectories used to plot Fig. 1 sampling algorithms over different kinds of automata: symbolic automata and extensions over the work of [START_REF] Barbot | Uniform Sampling for Timed Automata with Application to Language Inclusion Measurement[END_REF][START_REF] Barbot | Generation of Signals Under Temporal Constraints for CPS Testing[END_REF]. The modularity of Wordgen makes it easy to prototype new ideas.

DATA STRUCTURE AND IMPLEMENTATION

Wordgen is written in OCaml with some borrowing from other projects. It is freely available under the GPLv3 licence here [START_REF] Barbot | Wordgen git repository[END_REF]. The Uppaal file parser, time-automaton data structure and Differential Bound Matrix (DBM) library are taken from SymRob, a tool for the symbolic robustness analysis of time automata [START_REF] Sankur | Symbolic quantitative robustness analysis of timed automata[END_REF]. The automata data structure compute the synchronized automaton from a network of automata in the Uppaal formalism. The Prism file parser is taken from Cosmos [START_REF] Ballarini | HASL: A new approach for performance evaluation and model checking from concepts to experimentation[END_REF].

Rational numbers representation

The polynomials ring is built upon a rational field, two implementations are provided for rationals, either as floating point numbers or as exact rationals using arbitrary precision integers. Wordgen can use an exact rational number library based on GMP 1 . Floating point numbers are faster but prone to numerical errors.

Web application

The tool Wordgen features a graphical user interface available at [START_REF] Barbot | Wordgen Web Application[END_REF]. This webpage is a slightly modified version of the full tool compiled to JavaScript. It fully run in the browser without any installation. It featured two editors for the timed automaton either as a graph or as a Prism file and provide visualization of the results and internal representation using graphviz and gnuplot. Some functionalities are missing notably the Uppaal file parser and the speed of the tool is reduced compared to the native desktop version. It is very user-friendly compared to the native command-line version.

Usage of the command-line application

To run Wordgen on a unix commandline, one is to specify the automaton, the number of iteration of the distribution computation (the 𝑚 in equation 1) and the number of trajectories. Assuming the file twoears.prism contains the automaton describes in figure 3(a) we can run Wordgen as follows where --poly specify the parameter 𝑚 and --traj specify the number of trajectories: Wordgen features many options and parameters, the web application is a great way to explore them as they are displayed as a HTML form.

Signal Generation with Breach

To convert a timed word into a 𝑛-dimensional signal, i.e., a function from time in R + to R 𝑛 , the idea is to map the sequence of pairs of durations and labels to a sequence of control points which are vectors in the time domain R + ×R 𝑛 that Breach2 can use as a base for interpolation using different interpolation schemes, e.g., constant, linear, splines, etc. In practice, we need to define parameters for control points for each transition in a timed word of a given length, which is a tedious process. To make it easier, we implemented a wrapper class that only requires the user to map each event in the alphabet Σ to values or ranges of values for signals. To sample a signal, Wordgen sample a timed word, then Breach uses its builtin sampling methods to generate a value in the domain defined by each label of each transition.

As an example, consider the automata used in [START_REF] Barbot | Falsification of Cyber-Physical Systems with Constrained Signal Spaces[END_REF] to define inputs to an automatic transmission system for a car. We first indicate the name of the file defining the automaton, the timed word length and the names of the signals to generate, namely for throttle and braking. Then for each symbol in the alphabet of the automaton, we define either a single value or a range of values that the signals will be in when a label in a transition corresponds to this symbol. The resulting figure is shown on Fig. 4.

APPLICATION TO CPS FALSIFICATION

In [START_REF] Barbot | Falsification of Cyber-Physical Systems with Constrained Signal Spaces[END_REF] two case studies are proposed where Wordgen is used to falsify CPS by generating timed words interpreted as signals. We mentioned the driving automata in the previous section and give more details about the second one related to ΔΣ modulators. ΔΣ modulators are important components of analog-to-digital converters. Practical quantizers have limited input and output ranges, which may lead them to saturation, and we want to check whether the output ever saturates. Specifically, we consider the falsification of the absence of saturation of some quantizer signal Out under a certain class of nearly oscillatory inputs In. Formally In and Out must satisfy for some t 𝑠 ≥ 0 and ∀𝑡 ≥ 0,

|Out(𝑡)| < 2 (2) ∃𝑇 ∈ [8t 𝑠 , 12t 𝑠 ] such that In(𝑡 + 𝑇 ) = In(𝑇 ) (3)
Standard falsification approaches translates requirements in Signal Temporal Logic (STL). Encoding (2) as an STL formula is trivial: 𝜑 ¬sat = □|Out| < 𝑝 𝑠𝑎𝑡 . However, enforcing that In satisfies (3) is not so simple. For instance, unbounded periodic properties are known to be beyond STL expressivity [START_REF] Maler | Monitoring Temporal Properties of Continuous Signals[END_REF], and this is before considering that periods may be uncertain.

To solve this problem using Wordgen, an automata was designed as a model of a pseudo-periodic signal, shown on Fig. 5(a). The important part is the cycle of length 4, each transition of this TA has its own clock (𝑥 1 , 𝑥 2 , 𝑥 3 , 𝑥 4 ) which is reset when the transition is taken. where 𝑘 is the remainder of 𝑘/4. Since the discrete sequence 𝑢 𝑖 is periodic, the resulting signal satisfies (3) iff ∀𝑖, 8t 𝑠 ≤ 𝜏 𝑖 +𝜏 𝑖+1 +𝜏 𝑖+2 +𝜏 𝑖+3 ≤ 12t 𝑠 . Note that this constraint is satisfied by the delays of the timed words of our TA of Fig. 5(a). Hence by using Wordgen to generate timed words and mapping labels a, b, c, d to values 𝑢 0 , 𝑢 1 , 𝑢 2 , 𝑢 3 we obtain the desired signals.

For the saturation threshold 𝑝 𝑠𝑎𝑡 = 2 used in the model [START_REF] Brigati | Modeling Sigma-Delta modulator non-idealities in Simulink[END_REF], the property 𝜑 ¬sat was easily falsified in our optimization setting for arbitrary long timed word. Using other falsification methods required to use rejection sampling which was only tractable for signals of length up to three periods.

CONCLUSION AND FUTURE WORK

We have presented Wordgen, a uniform sampler for timed language. To our knowledge it is the first standalone tool able to perform this sampling. Compared to previous tool, Wordgen is faster, easier to use, feature a graphical user interface and have more features. It has been successfully used for signal generation in the context of falsification of cyber-physical systems where other approaches fail.

Wordgen was designed as a fast prototyping tool, in the future we plan to incorporate several new features. By using Boltzmann sampling, we will be able to sample words of different lengths. Using recent work on time automaton volume computation, we will be able to control the average duration of generated timed word and to deal with unbounded clocks.

Figure 1 :

 1 Figure 1: (Left) A TA with 2 clocks 𝑥 and 𝑦, 1 location and 2 transitions. The transition labeled by 𝑎 is guarded by upper-bounds on clocks, while the guard on the transition labeled by 𝑏 is a small box 𝑥 ∈ [8; 10] ∧ 𝑦 ∈ [1; 2]. The only way to reset the clock 𝑥 is to hit the box. (Top Right) discounted uniform sampling with 10 iterations (m=10), (Bottom Right) isotropic sampling, each sampling contains 50000 trajectories of length 100 where a pixel is drawn at each transitions at positions (𝑥, 𝑦) before applying resets.

Figure 2 :

 2 Figure 2: Architecture of Wordgen. Round part represents computation steps while square one represents data. Arrows represent flow of data. Lines are dependencies.

Figure 3 :

 3 Figure 3: Illustration of the different steps of computation. (a) a TA; (b) the forward reachability zone graph, red edges are already split; (c) the split zone graph; (d) weights and distribution computed on edge from vertices 1 to 2 with 0 to 4 iterations; (e) a timed word sampled by Wordgen.

>Figure 4 :

 4 Figure4: Three trajectories generated from the driving automata of[START_REF] Barbot | Falsification of Cyber-Physical Systems with Constrained Signal Spaces[END_REF].

  TA_filename = ' driving_TA . prism '; num_evts = 10; signals = { ' throttle ' , ' brake ' };

Figure 5 :A

 5 Figure 5: (a) A timed automaton for pseudo-periodic signal. To avoid overloading the figure, a global clock 𝑥 (reset to 0 at each transition) and the global invariant 𝑥 1 , 𝑥 2 , 𝑥 3 , 𝑥 4 < 12 and 𝑥 < 4 hold for each guard and are not depicted. (b) Five time words of length 40 sampled from (a) with the discounted uniform sampling with 𝑚 = 5, each action a,b,c,d are given a value shown in ordinates. % label init : Throttle range [0 , 100] but brake is 0 labels_ranges . init . brake = 0; labels_ranges . init . throttle = [0 100]; % label a: no acceleration and no braking labels_ranges .a. throttle =0; labels_ranges .a. brake =0; % label b: no acceleration , braking in range [0 , 350] labels_ranges .a. brake = [0 350]; labels_ranges .b. throttle =0; % label c: no braking labels_ranges .c. brake =0; (...) A BreachTASignalGen can then be created and used to generate signals. It is derived from Breach native class BreachSignalGen. B_ta = BreachTASignalGen ( signals , TA_filename , labels_ranges , num_evts ); B_ta . SampleDomain (5) % sample 5 timed words calling wordgen B_ta . Sim (0:0.01:50) ; % interpolate control points for t in [0 , 50] B_ta . PlotSignals () ;

  Fig. 5(b) depicts time words sampled from this automaton by Wordgen. The idea is to extend this automaton into a signal generator interpolating the signal values between points of a periodic discrete sequence of the form:𝑢 0 𝜏 0 𝑢 1 𝜏 1 𝑢 2 𝜏 2 𝑢 3 𝜏 3 𝑢 0 𝜏 4 𝑢 1 𝜏 5 𝑢 2 𝜏 6 𝑢 3 𝜏 7 𝑢 0 𝜏 𝑁 𝑢 N . . .. The value In(𝑡) is obtained by finding 𝑘 such that 𝑘 0 𝜏 𝑖 ≤ 𝑡 < 𝑘+1 0 𝜏 𝑖 and interpolating between 𝑢 𝑘 and 𝑢 𝑘+1

•

  Cosmos is designed to simulate large concurrent systems specified as Petri Nets; in our setting we are only simulating single automata and thus many optimizations of Cosmos are pure overhead. In Table1we can see that Wordgen and Cosmos sampling times are of the same order of magnitude, Wordgen being faster for small automata.

	𝑚 Sage Cosmos		Wordgen
			Sampling	Float Rationnal Sampling
	1	2.4	163	0.02	0.03	51
	5	4.8	179	0.1	0.1	122
	10 16.2	265	0.5	0.12	208
	15	48	383	2.2	3.1	374
	20 114	542	6.7	9.7	475
	Table
						Moreover, these tools have a large constant startup time (Loading
						the JVM for Prism, loading Python libraries for Sage, code genera-
						tion and compilation for Cosmos) making running the tool-chain a
						long and tedious process when designing a given model.
						Usability. Making the tools of the legacy tool-chain work together
						is quite cumbersome. Indeed, each of them is written in a different
						programming languages (Java, Python/Sage , C++) which requires
						to serialize data to files between each operation. The coordination
						between the three tools requires several scripts and adding new

features means implementing modifications in each individual tools as well as in the file format of the exchange file. New features. In Wordgen the data-structure to represent distributions can be replaced easily. On-going work uses this to implement

https://github.com/decyphir/breach

This work was financed by the ANR MAVeriQ (ANR-20-CE25-0012), the join ANR-JST project CyPhAI and the Auvergne-Rhône-Alpes region project DETAI..