Anisse Ismaili
email: anisse.ismaili@lip6.fr

Evripidis Bampis
email: evripidis.bampis@lip6.fr

Nicolas Maudet
email: nicolas.maudet@lip6.fr

Patrice Perny
email: patrice.perny@lip6.fr

A Study on the Stability and Efficiency of Graphical Games with Unbounded Treewidth

Keywords: I.2.11 [Distributed Artificial Intelligence]: Multiagent systems Algorithms, Theory, Experimentation Graphical Games, Approximation

Graphical games (GG) provide compact representations of multiplayer games involving large populations of agents when influences among them have some locality property. The notion of pure Nash equilibrium (PNE), not requiring randomized strategies, is a fundamental stability concept. However, recent results show that for many natural topologies, a PNE is very unlikely to exist when the number of agents is large, which challenges the relevance of PNE in large GG.

In this paper, we investigate how far we can get from the notion of individual stability captured by the concept of PNE, by only requiring agents to be almost in best-response ( -Nash), or by requiring almost all agents to be in bestresponse. We study these approximated notions of PNE for different topologies, including graphs with unbounded treewidth, like grids. This makes the problem computationnally very challenging and requires the comparison and use of several algorithmic solutions. Our results reveal surprisingly good asymptotic properties, tempering the claim that individual stability is not a relevant notion for large GG. Finally, as approximated PNE provide various tradeoffs between stability and social utility maximization, we propose an approach to construct a minimal-size e-covering of all feasible Pareto-dominant tradeoffs.

INTRODUCTION

Games involving a large number of players require a prohibitive amount of space in order to be represented in normal form. However, in many cases, the payoff functions of agents do not depend on all other agents. Graphical games (GG) provide compact representations of multiplayer games involving large populations of agents when influences among them exhibit such locality property [START_REF] Kearns | Graphical models for game theory[END_REF][START_REF] Kearns | Graphical games[END_REF]. In a GG, each vertex is an agent, and an agent's payoff is a function of his action and the actions of the agents in the innerneighborhood.

This paper investigates the stability and efficiency of graphical games, with agents playing pure strategies (i.e. not randomized). In this context, the notion of Pure Nash Equilibria (PNE) is arguably the most fundamental stability notion: no agent has an individual incentive to deviate. Stability usually comes at the price of efficiency: the ratio of the best (resp. worst) NE over the social welfare optimum is known as the price of stability [START_REF] Anshelevich | The price of stability for network design with fair cost allocation[END_REF] (resp. price of anarchy).

Unfortunately, there is no guarantee for such a PNE to always exist (in contrast to the case of mixed Nash equilibria). The existence problem of PNE in general GG has been known to be hard for a long time [START_REF] Gottlob | Pure Nash equilibria: Hard and easy games[END_REF].

By exploiting non-serial dynamic programming [START_REF] Bertele | Nonserial dynamic programming[END_REF] as used in graphical models for probabilistic inference, [START_REF] Daskalakis | Computing pure Nash equilibria in graphical games via markov random fields[END_REF] provides an algorithm which handles efficiently graph structures with bounded (or logarithmic, when the number of actions is bounded) treewidth. One can also cite [START_REF] Soni | Constraint satisfaction algorithms for graphical games[END_REF] which assumes a bounded treewidth too. Recently, [START_REF] Jiang | Pure Nash equilibria: complete characterization of hard and easy graphical games[END_REF] even showed that this criteria suffices to fully characterize hard vs. easy instances in GG. Essentially, GG are easy if and only if they have a bounded treewidth (after their sinks have been iteratively removed). The case of unbounded treewidth is nevertheless of great interest, since many natural graphs fall in this category (e.g. grids).

Another important issue when there are several possible PNE is to search for the best one in terms of some welfare function [START_REF] Chapman | A distributed algorithm for optimising over pure strategy Nash equilibria[END_REF]. However, recent results show that when the number of agents becomes large, the likelihood to find at least one PNE converges to 0 for many graph topologies. Indeed, in [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF] different natural topologies are investigated, both theoretically and experimentally. They show that the structure heavily affects the probability of existence of a PNE: for instance, while the probability converges to 0 in the case of trees with an unbounded diameter, it remains as high as 1 -1 e for the standard (random) normal form, or bipartite graphs. They suggest that the probability of existence increases as the length of paths diminishes. From the results of [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF], there is no evidence that the probability of PNE existence in general graph topologies tends to a strictly positive value when the diameter diverges more quickly than Ω(log(n)). The same problem is addressed in [START_REF] Daskalakis | Connectivity and equilibrium in random games[END_REF] for graphs drawn from the Erdös-Renyi G(n, p) model (where n is the number of vertices and p ∈ [0, 1] is the probability for an unoriented edge to exist). They essentially show that when c n 2 < p < ln(n) 2n (medium connectivity) the probability for a PNE to exist converges to 0, and when p > 2 ln(n) n (high connectivity) the number of PNEs converges to a poisson(1) distribution (the probability of existence is again 1 -1 e ). These results may challenge the relevance of the notion in GG involving a large number of agents (at least for the aforementioned topologies). However, they provide a partial picture of the problem. In particular, they mostly pay attention to the 2-action case and remain constrained by the exact definition of PNE. In particular, it does not address the question of whether it is possible to find a profile of strategies which is "almost" a PNE. This is the first problem that we study in this paper. This requires to define what an approximate PNE is. As a PNE requires all agents to strictly play a best-response, we shall consider two possible approximations in this paper. The first one is to modify the condition under which agents may deviate. Specifically, we consider the so-called multiplicative -Nash equilibrium. The rationale of this notion is that the relative incentive to deviate for agents may be small in practice. The second one that we call k-Nash consists in relaxing the constraint that all agents are in best-response. The rationale of this approximation is that it may only require, for the system designer, to "convince" or "impose" the strategy of a small number of agents. This notion bears some similarity with the notion of Stackelberg threshold [START_REF] Sharma | Stackelberg thresholds in network routing games or the value of altruism[END_REF], with the difference that this measures the number of agents whose strategy must be fixed in order to guarantee a social optimum at equilibrium, while in our case, we seek to guarantee that an equilibrium exists. We also note that the study of k-Nash approximations is very fundamental, since it does not rely on the payoff random distributions, but just on best response tables.

Unfortunately, as previously mentioned, instances involving unbounded treewidth are computationally very challenging to handle. We deploy a range of algorithmic solutions (junction tree algorithms [START_REF] Daskalakis | Computing pure Nash equilibria in graphical games via markov random fields[END_REF], a SAT formulation due to [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF], and Mixed Integer Linear Programming), and compare their respective efficiency on the different problems we face. This allows us to perform different experiments involving a significant number of agents (typically at least 100 agents) and action-sets of size greater than 2. Our results shed a new light on the PNE problem in large GG: in a nutshell, we show that in various GG whose PNE probability of existence converges to 0, we actually get very close to PNE (in particular, a very small constant fraction of agents may be unsatisfied), and that the action-set size may greatly affect the results (in particular, it contradicts the assumption that the divergence of the graph diameter may be a sufficient condition to conclude on the asymptotic inexistence of PNE).

We next turn our attention to the tradeoff between stability and utilitarian optimality. Equipped with a suitable approximation of PNE, it is actually possible to generalize the notion of price of stability by seeing it as a bi-criteria problem. The set of non-dominated tradeoffs would provide precious indications for the system supervisor ("what gain in social welfare can be expected by relaxing this much of the stability requirement?"). Specifically, we question how the price of stability evolves as a function of the best response -relaxation, by computing the values of an (1 + e)approximated bi-criteria Pareto-set ( , v) of minimal size, with a minimal number of calls to an utilitarian optimization problem and an -stability optimization problem, thanks to a greedy algorithm defined in [START_REF] Diakonikolas | Small approximate Pareto sets for biobjective shortest paths and other problems[END_REF].

The remainder of this paper is as follows. In Section 2 we introduce the necessary background on graphical games. In Section 3, we provide the different algorithmic approaches we use, and discuss how they perform for the different problems we have to deal with. Approximated notions of PNE are introduced and investigated in Section 4. Finally, Section 5 details the methodology to obtain an approximated Pareto-set for the stability vs. efficiency tradeoff.

PRELIMINARIES

We remind common notations of multiagent games. The set of agents is N = {1, . . . , i, . . . , n}. Each agent i chooses his action ai in his particular action-set Ai. An action-profile is any element: (a1, . . . , ai, . . . , an) ∈ n i=1 Ai denoted by a ∈ A. An agent i's adversary-action-profile is denoted a-i ∈ j =i Aj and his payoff function is vi :

Ai × A-i → R. Definition 1 A multiagent game (or n-agent game when |N | = n) is a t-uple Γ = N, {Ai} i∈N , {vi} i∈N .
Nash equilibria formulate the "individual stability" of an action-profile (ai, a-i) ∈ A: each agent i, given the fixed adversary-action-profile a-i, has no incentive to change unilaterally his action ai, to get an individually prefered actionprofile.

Definition 2 Given a multiagent game Γ and a player i, a best response function BRi : A-i → P (Ai) maps adversaryaction-profiles to the set of actions which "satisfy" i.

Definition 3 An action-profile a ∈ A is a Nash equilibrium if and only if for every agent i, ai ∈ BRi (a-i).

If we consider BRi (a-i) = argmax a i ∈A i {vi (ai, a-i)}, then we are working on Pure Nash Equilibria (PNE). If we con-

sider BR i (a-i) = {ai ∈ Ai|∀a i ∈ Ai : v i( a i ,a -i) v i( a i ,a -i) ≤ (1 + )}
(assuming positive utilities) , then we are working on multiplicative -approximated Nash Equilibria ( -PNE).

In the most general normal forms, each payoff function requires O i∈N |Ai| space, which is not compact if the number of agents is unbounded. A graphical game (GG) catches the locality of payoff functions. Its support graph is an oriented graph G = (N, E), where N is the set of agents, E are the edges which define dependencies between payoffs, and agent i's inner-neighborhood (including himself) is denoted ν(i). Player i's payoff function depends only on the actions of players in ν(i). A local-action-profile is any element: (ai, . . .) ∈ j∈ν(i) Aj denoted by: a (i) ∈ A (i) . In a graphical game, agent i's payoff is a function vi : A (i) → R, which by the way requires O j∈ν(i) |Aj| space. Hence GG are a compact multiagent games representation when ∆ = maxi∈N {|ν(i)|} is bounded and small.

Definition 4 A graphical game (GG) with support graph G = (N, E) is a t-uple Γ = G, {Ai} i∈N , {vi} i∈N .
What we are mostly going to consider are unoriented graphical games, where (i, j) ∈ E ⇔ (j, i) ∈ E (payoff dependencies between agents are mutual). We might also want to find an action-profile which maximizes the average payoff, like in the VCG mechanism and Clarke's prices.

Definition 5 Given a GG Γ = G, {Ai} i∈N , {vi} i∈N the utilitarian optimization problem (denoted as max -GG) is to find a * ∈ argmax a∈A { i∈N vi(a (i) )}.
Given a fixed support graph G = (N, E), we study randompayoff graphical games with respect to the following distribution: each payoff values vi(a (i) ) of the functions vi : A (i) → R are drawn uniformly and independently from the interval [0, 1]. For the PNE and the k-Nash problems, it is equivalent to suppose that each adversary-local-action-profile has a singleton best-response. But for the -Nash problem, it is an important (but common) hypothesis that can influence the results.

Throughout this paper, B (i) = {a (i) ∈ A (i) |ai ∈ BRi (a-i)} will denote the best-responses local-action-profiles, and

B (i) = A (i) \ B (i) the non-best-responses local-action-profiles. For -best-responses, < ⇒ B (i) ⊃ B (i) .

COMPUTATION IN GRAPHICAL GAMES

In this section, we present various ways to compute a PNE in a general GG. Recall that this is an NP-Hard problem. It is known that this problem's structure is close to a logical problem [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF]. In what follows, we first remind how junction trees can be used when (hyper)treewidth is bounded. Then we detail the SAT formulation that [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF] briefly proposes and massively uses. We also propose a MIP formulation of the same problem, and we show how this MIP can be adapted to the case of an unbounded maximal degree ∆, with a cutting plane algorithm polynomial in ∆, for non-best-responses.

Junction Tree Algorithm. Computing a PNE in a GG Γ = G, {Ai} i∈N , {vi} i∈N can be reduced in various ways to the Junction Tree algorithm. For instance, [START_REF] Daskalakis | Computing pure Nash equilibria in graphical games via markov random fields[END_REF] builds a maximum likelihood problem in a Markov random field (ML-MRF) over A1 × . . . × An which is clique-factorized and whose cliques roughly correspond to the agents neighborhood: Ci = ν(i) and P (a) = i∈N φC i (a (i) ). The potential functions are

φC i (a (i) ) = 1 if a (i) ∈ B (i) otherwise
in such a way that the maximum likelihood are 1 if and only if there is a PNE. With the same reduction, one can prove the NP-Hardness of max -GG. It is sufficient to change the multiplicative algebra into an additive algebra.

One can remark that if we choose = e -1 = exp(-1), this ML-MRF not only decides on the existence of a PNE, but also maximizes the number of agents in a best-response local-action-profile. Indeed, if π * is the maximal log-likelihood of this ML-MRF, then -π * is precisely the minimal number of unsatisfied agents. Unfortunately, these algorithms assume a bounded degree (or at most logarithmic with a fixed number of actions) treewidth, for tractability. They behave very badly in the general case, including the St Paul Town-Game.

SAT. Let us now introduce a decision variable xi,a i ∈ {0, 1} for each agent i and each one of his action ai. These variables can be both interpreted as a boolean or a binary variable, whose meaning is "agent i chooses action ai". From this, given a GG Γ = G, {Ai} i∈N , {vi} i∈N , the SAT formulation proposed in [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF] is obtained by two sets of clauses. The first one ensures that each agent chooses one and only one action: ∀i,

a i ∈A i xi,a i ∀i, ∀{ai, a i } ⊂ Ai, ¬xi,a i ∨ ¬x i,a i
The second set of clauses forbids the non-best-responses local-action-profiles :

∀i, ∀a (i) ∈ B (i) , j∈ν(i) ¬x j,a (i) j
One can remark that, this formulation can also be used to answer for a fixed the decision problem "Is there an -Nash Equilibrium?", by simply prohibiting B (i) .

MIP. For the first MIP formulation that we propose, we use the same variables and the same constraints, and we introduce the relaxation variables ρi ∈ [0, 1] the sum of which must be minimized so that a PNE exists if and only if the optimal value is 0. First of all we need to add constraints enforcing that i chooses exactly one action:

∀i ∈ N, ∀ai ∈ Ai : xi,a i ∈ {0, 1} (1)
∀i ∈ N :

a i ∈A i xi,a i = 1 (2)
Then we impose that for each non-best-response local-actionprofile, at least one decision of the inner-neighborhood ν(i) is 0. These contraints are relaxed by using continuous variables ρi.

∀i ∈ N : ρi ∈ [0, 1] (3) ∀i ∈ N, ∀a (i) ∈ B (i) : j∈ν(i) x j,a (i) j ≤ |ν(i)| -1 + ρi (4)
We can then write the MIP addressing the PNE problem in GG: minimize

n i=1 ρi subject to (1) ∧ (2) ∧ (3) ∧ (4)
By introducing a cutoff if maxi ρi > 0, we avoid many exploration nodes that are useless for the Pure NE problem. The number of non-best-response local-action-profile constraints is an exponential of ∆. Hence this is a compact formulation if ∆ is bounded.

Techniques. We tested three computation techniques: maximum likelihood in a Markov random field with the junction tree algorithm implemented in libDAI [START_REF] Mooij | libDAI: A free and open source C++ library for discrete approximate inference in graphical models[END_REF], the SAT formulation with the miniSAT solver [START_REF] Een | An extensible SAT-solver[END_REF], and the MIP formulation in GuRoBi [START_REF]Gurobi Optimization[END_REF]. We shall summarize how these approaches fit the different problems considered in the following section.

Figure 2 compares the computation times in unbounded (hyper)treewidth GGs, using our MIP formulation (ILP) and junction-trees (ML MRF). Both techniques are surpassed by miniSAT (which required an almost 0 time, here).

APPROXIMATING PURE NE

In this section we investigate two different kinds of approximation of pure equilibria in GG. We conduct experiments with these approximations on support graphs whose random payoff GG are known to be unfortunately not likely to have a PNE when n is high: for paths, binary trees [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF], and medium connected Erdös-Renyi [START_REF] Daskalakis | Connectivity and equilibrium in random games[END_REF], where we know that P (∃a : PNE ) → 0, when n → ∞, and for square grids, for which we experimentally measured (see Figure 4) that the PNE probability tends to 0 for small action-sets.

-Nash.The first notion we study is the multiplicative approximated NE, which minimizes (a) = maxi∈N { i(a (i) )}

(where i(a (i) ) is defined by (1 + i(a (i) ))vi(a (i) ) = βi(a (i) ) and βi(a (i) ) is the best individual value that i can achieve with an individual deviation). This general idea of this approximation is standard in the literature (although we note that it is usually used in the additive sense), and reflects the intuition that "agents are indifferent to sufficiently small gains" [START_REF] Shoham | Multiagent systems: Algorithmic, game-theoretic, and logical foundations[END_REF]. We emphasize that this notion is often used as a discretization to facilitate the computation of mixed equilibria. In this case, a common criticism is to say that it is not guaranteed to be close to an exact equilibrium (which does exist for sure). This argument does not hold here, as our objective is to minimize : when a PNE does exist, we indeed get = 0, otherwise we find an action profile which is the less sensible to individual deviation: the smallest such that there is an action profile where agents are not likely to individually deviate when the gain is less than (1 + ) times what they already have.

k-Nash. The second kind of approximation is to find an action-profile, which minimizes the number k of agents not playing a pure best-response (k-Nash), and hence to find the minimal fraction of unsatisfied agents that it is necessary to "convince" or "constrain" to obtain stability. This notion has a great advantage over the -Nash approximation: it does not rely on the payoff distributions, but just on singleton best response tables. We note that, in principle, the fraction of unsatisfied agents may converge to 0 even though the probability of PNE does so.

Methodology. Remark that both the maximum log likelihood and the MIP presented in Section 3 can be directly adapted to k-Nash: maximizing the number of agents playing a best response. If the support graph is of polynomial size and has a bounded treewidth, then one should use a junction tree algorithm. Otherwise, when the treewidth is unbounded, or when ∆ is unbounded, one should use the MIP formulation when addressing k-Nash.

Now for the computation of -Nash, since for a fixed , the SAT formulation is the fastest one by a huge margin (regardless of the treewidth or the support graph's maximal degree ∆), it turns out to be faster than MIP to make a dichotomy over with the SAT formulation at each iteration.

We deduced the proportions of instances having a PNE from the fastest previous approximations. For the convenience of the reader, we summarize the characteristics of the graphs and the algorithmic approaches employed (for the different problems) in Table 1.

Results on paths and binary trees. We know from [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF] that the probabilities of PNE to exist in paths and binary trees tend to 0 (they actually reach 0 with about 100 agents). Our experimental results are presented in Figure 3. We observe that these results are insensitive to the size of the action-set. On these structures, more actions for agents do not help to reach a PNE. Now regarding the approximations of PNE, one may intuitively expect the situation to worsen as the probability of PNE decreases. This is not the case: in both topologies, the k-Nash equilibria we found reveal that the average minimal fractions of unsatisfied agents tend to small constants, between 3% and 6%. It means that the average number of unsatisfied agents is strictly proportional to the total number of agents. Moreover, this convergence is very fast: even when the probability of PNE is as high as 50% (path, 2 actions), the average fraction of unsatisfied agents is already at 6%. In this sense, we are not far from individual rationality equilibrium (even though the number of Unlike for k-Nash, the experiments for -Nash show that in paths and trees, the results are affected by the size of the action-set. More precisely, the smaller the action-sets are, the harder it is to guarantee a small incentive for all agents in the same action-profile. To conclude, although more actions available for agents do not increase the probability of PNE, they do increase significantly the quality of the -Nash approximation (compare 200 agents on a path: is less than 10% for 4 actions vs. 50% for 2 actions). Results on grids. These results are in Figure 4. The probability of PNE existence is greatly affected by the actionset size in this topology (recall this was not the case for paths and binary trees). For 2 or 3 actions for each agent, these probabilities tend to 0 (at different rate) as the grid sizes increase. With action-sets of size 4, the influence of the action-set size can be observed even more dramatically: the probabilities of pure NE to exist do not tend to 0 any longer (they seem to tend to the 1 -1 e value). This shows that the diameter of the graph cannot be a sufficient parameter to conclude on the fact that the probability of PNE existence converges asymptotically to 0 or not. To the best of our knowledge, this is the first time that this is observed.

Turning our attention to k-Nash, we see that we are not so far from individual stability since the minimal fraction of unsatisfied agents are below 2%. About -Nash, tends to 0 with 3 or 4 actions, but does not with 2. We also note that this topology allows to reach very close to stable -Nash approximations: for instance, remains negligible for 250 agents with a 3-action set, when the probability of PNE has already decreased to 20%. ) Erdös-Renyi graphs, we see that the convergence to 0 of PNE existence is not affected by the action-set size (not shown here). We can also observe that the average minimal fraction of unsatisfied agents decreases below 2%. As for the -Nash, the values obtained show that it is more difficult than with grids to get a close -approximation of PNE. For completeness, we mention that we also performed some experiments on small worlds (Watts-Stogratz) graphs, which are free-scale networks capturing some properties of biological, social, and computer networks. To sum up, we observed that the actionset size has a strong influence on both the PNE probability and the minimal values. By focusing our attention to p ∈ [0.03, 0.375], we observe that we obtain very small values (from 3% to less than 1%). Note that when p = 0.375 the PNE probability has reached the 1 -1 e value [START_REF] Dilkina | The impact of network topology on pure Nash equilibria in graphical games[END_REF]. Conclusions. Some general conclusions can be drawn from these experiments. First of all, even though the PNE probability tends to 0, in general, the approximated equilibria are very good. This is in line with some observations of [START_REF] Vickrey | Multi-agent algorithms for solving graphical games[END_REF] who mentioned that "it appears that, with random payoffs, fairly good equilibria often exist in pure strategies" after some experiments on grids (with 3-actions games, and under the additive -Nash interpretation). However, we also see that the picture is more complex. While k-Nash approximations always require only a small (decreasing or constant) fraction of agents to be convinced, it may sometimes be difficult to get small -approximations when the action-set is 2. But augmenting the size of the action-set always improves this approximation. As approximated PNE may provide various solutions, this gives us the flexibility to optimize the efficiency (in the utilitarian sense).

STABILITY VS SOCIAL UTILITY

Given a GG, we investigate the bi-criteria problem that consists in both maximizing "individual stability" and "overall utility". Formally, concerning stability, we will focus on the maximization of a strictly decreasing function of , ϕ(a) = 1/(1 + (a)), which corresponds to the guaranteed fraction of individual rationality ϕi(a (i) ) = vi(a (i) )/βi(a (i) ) where βi(a (i) ) is the best individual value that i can achieve with an individual deviation. Note that maximizing this function means that we minimize (a), which represents the worst relative opportunity loss among agents for actionprofile a. Concerning utility, we maximize the standard utilitarian criterion defined by v(a) = i∈N vi(a (i) ). Hence, to any a in A is assigned a performance vector (ϕ(a), v(a)) in the bi-criteria valuation space. The set of feasible performance vectors is then defined by F (A) = {(ϕ(a), v(a)), a ∈ A} and the components are assumed to be bounded.

Since the two objectives considered are possibly conflicting, we are interested in finding action-profiles achieving different tradeoffs. To this end, determining a representative subset of feasible Pareto-optimal tradeoffs within F (A) would provide useful information for the system supervisor. We recall that a vector (ϕ, v) ∈ F (A) is Pareto-optimal if we cannot improve one dimension without downgrading the other. Note that the set P of Pareto-optimal vectors in F (A) possibly includes a huge number of elements. Not only would its exact determination induce a prohibitive computational time, it would also provide the supervisor with an unnecessary large sample of possibilities. We propose instead to determine an approximation of P using the notions of e-dominance and e-covering. In general, multiple e-covering of the Pareto set P exist. The set P itself in an obvious example of e-covering of P. More interestingly, the existence of e-covering of polynomial size is established for general multiobjective combinatorial problems with a bounded number of criteria [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF]. Here, for any given e = (eϕ, ev) ∈ R 2 + we would like to be able to determine an e-covering of P as small as possible, ideally an e-covering of P of minimal cardinality. Since we have only two objectives we can resort to a greedy approach proposed in [START_REF] Diakonikolas | Small approximate Pareto sets for biobjective shortest paths and other problems[END_REF] for computing an e-covering of minimal cardinality in general bi-criteria optimization problems. Let us explain how this approach can be further specified to explore the feasible tradeoffs attached to the action-profiles in a GG.

The construction proposed in [START_REF] Diakonikolas | Small approximate Pareto sets for biobjective shortest paths and other problems[END_REF] relies on solving a sequence of optimization problems alternating two complementary subproblems:

Restrict-v(α). For any given value α, we want to maximize v subject to the constraint ϕ ≥ α, or answer no when no such solution exists.

Restrict-ϕ(γ). For any given value γ, we want to maximize ϕ subject to the constraint v ≥ γ, or answer no when no such solution exists.

We present now how we can solve these two subproblems in our context. Since these two problems are NP-hard, and since they include too much numerical information to resort to SAT reformulation, we focus on Mixed Integer Programming formulations (MIP).

Solving Restrict-v(α). With no loss of generality, we assume that payoffs are positive (vi(a) ≥ 0). If it is not the case, a simple translation of payoffs is sufficient. We remind players action variables xi,a i ∈ {0, 1} whose value is 1 iff player 1 chooses action ai. First of all we need to add constraints enforcing that i chooses exactly one action :

∀i ∈ N, ∀ai ∈ Ai : xi,a i ∈ {0, 1} (5) ∀i ∈ N : a i ∈A i xi,a i = 1 (6)
Then we define the variables of local action profiles z i,a (i) , in such a way that we have z i,a (i) = j∈ν(i) x j,a (i)

j : ∀i ∈ N, ∀a (i) ∈ A (i) : 0 ≤ z i,a (i) ≤ 1 (7) ∀i ∈ N, ∀a (i) ∈ A (i) , ∀j ∈ ν(i) : z i,a (i) ≤ x j,a (i) j (8)
And make players payoffs vi consistent with the activated local action profile:

∀i ∈ N : vi = a (i) ∈A (i) vi(a (i) )z i,a (i) (9)
To complete the formulation of Restrict-v(α), we must also require that ϕ ≥ α:

∀i ∈ N, ∀a (i) ∈ B (i) α , z i,a (i) = 0 (10)
To avoid very large integrality gaps, the following reinforcement is useful:

∀i ∈ N, a (i) ∈A (i) z i,a (i) = 1 (11)
Moreover constraints [START_REF]Gurobi Optimization[END_REF] become necessary when the problem is empty due to constraints [START_REF] Gottlob | Pure Nash equilibria: Hard and easy games[END_REF]. We can then write our MIP formulation solving Restrict-v(α): maximize n i=1 vi subject to (5) ∧ (6) ∧ (7) ∧ (8) ∧ (9) ∧ (10) ∧ [START_REF]Gurobi Optimization[END_REF] Since all vi(a (i) ) are positive, and we maximize, the local action profiles variables are consistent.

Solving Restrict-ϕ(γ). In this problem, we want to maximize ϕ(a) = 1/(1 + (a)) with (a) = maxi∈N { i(a)} and i(a) = a (i) ∈A (i) i(a (i) )z i,a (i) . We recall that our actual objective is to minimize i, i ∈ N under the constraint that the total payoff is at least γ. However, it is not possible to proceed exactly as above to formulate the restrict-ϕ problem using variables z i,a (i) because now the linearizations as well as the consistencies of z i,a (i) = j∈ν(i) x j,a (i) j are lost. To recover the necessary monotonicity of the objective function with respect to z i,a (i) variables we maximize a strictly decreasing function of . We therefore define ϕi(a

(i) ) = 1 1+ i (a (i) ) = v i (a (i) ) β i (a (i)
) . Hence we have :

∀i ∈ N : ϕi = a (i) ∈A (i) v i (a (i) ) β i (a (i) ) z i,a (i) (12)
The total payoff must be at least γ, therefore :

i∈N a (i) ∈A (i) vi(a (i) )z i,a (i) ≥ γ (13)
Finally, our MIP adressing Restrict-ϕ(γ) is :

maximize mini∈N {ϕi} subject to (5) ∧ (6) ∧ (7) ∧ (8) ∧ (11) ∧ (12) ∧ (13)
If we assume v i (a (i) ) β i (a (i) ) ≥ 0 (without loss of generality), then the local-action-profile variables are consistent, and the MIP is valid. We therefore obtain MIP formulations to solve the two basic sub-problems Restrict-v(α) and Restrict-ϕ(γ). Now, the greedy construction of an e-covering starts with the initial call q0 =Restrict-v(0). Then we compute the following alternated sequences for n ≥ 1:

pn = Restrict-ϕ(qn-1/(1 + ev)) (14)
qn+1 = Restrict-v((1 + eϕ)ϕ(pn)) (15)
We let n increase until the feasible domain of Restrict becomes empty. The resulting set {p1, . . . , pK } provides an e-covering set of minimal cardinality. This procedure makes only 2K calls to Restrict where K is the size of the output. Further details on this greedy approach and its optimality for general bi-objective problems are given in [START_REF] Diakonikolas | Small approximate Pareto sets for biobjective shortest paths and other problems[END_REF].

We have performed numerical tests to check the practical feasibility of this approach. The results can be seen on Figure 6 that provides e-coverings for instances of size 7×7 and e ∈ {(0.05, 0.25), (0.015, 0.075), (0.002, 0.01)}.

As expected, we observe that points belonging to the ecovering set spread throughout the approximate Pareto front with a density of points that increases as the norm of e = (eϕ, ev) diminishes. Note that the supervisor can control independently the density of points on the v-axis and on the ϕ-axis by choosing ev = eϕ. Remark also that some of these points do not belong to the boundary of the convex-hull of the e-covering set. Such tradeoffs are known as unsupported solutions and cannot be obtained by optimization of a linear combination of the two criteria. The numerical tests confirm the adequacy of the e-covering concept to provide representative approximations of the Pareto ensuring some diversity (inherent to the covering property) while avoiding any redundancy (due to minimality). Such coverings provide the supervisor with a collection of typical action-profiles achieving various efficient tradeoffs between stability and social utility maximization. The final choice within the e-covering can then easily be made depending on the aspirations in terms of stability and social welfare.

CONCLUSION

In this paper we investigated stability and efficiency issues in GGs, without confining ourselves to the tractable case of bounded treewidth. This is made possible through the use of a variety of algorithmic solutions. Previous works showed the influence of the network topology on the asymptotic probability of PNE, exhibiting many cases where this probability tends to 0 with n. What our experimental results first exhibit is that there is an interplay between the topology and the action-set size: for some topologies only, the action-set makes a crucial difference in terms of convergence to 0 for the PNE existence. Thus, we observed a (1 -1 e ) probability of PNE existence in square-grids with 4 actions. To the best of our knowledge, this is so far the topology with the fastest diameter's divergence ( √ n) observed not to tend to 0. Our second main finding is that, despite the low probability of PNE existence, approximations of equilibria are in most cases close to individual stability. The convergences of k-Nash to particular values calls for a deeper understanding of how graph topologies influence the minimal number of unsatisfied agents. Eg., it is an open question whether there are any topology where the k-Nash over random-payoff GG may increase. While k-Nash does not rely on the payoff random uniform [0, 1] distribution, -Nash does, so it would also be of interest to study non-uniform distributions.

Finally, we have seen that these approximations of PNE provide many solutions. This allows the computation of tradeoffs between almost individually stable solutions, and efficient (socially optimal) one. Knowing how the price of stability would evolve as a function of the would certainly be useful for the system designer. Our last section shows how this can be performed in practice without returning a set of undominated trade-off solutions of prohibitive size. We want to emphasize that the MIP formulations we used for the tradeoffs determinations can be generalized to other criteria. Eg., PNE can be -approximated in an additive way, or the overall utility can be replaced by a fair criterion.

Figure 1 :

 1 Figure 1: St Paul Town-Game

Figure 2 :

 2 Figure 2: Computation average times for the PNE problem in square grids made of 2 actions agents.

Table 1 :

 1 Dichotomy of SAT [1500] medium connect. Erdös-Renyi Ω(βn), with β > 0 [9] log n MIP [500] Dichotomy of SAT [500] Summary of graph parameters and algorithmic solution for the different problems. The number in brackets indicates the number of payoffs GG per graph and action-size

Figure 3 :

 3 Figure 3: Indicators of PNE proximity in randompayoff GG : path and binary trees

Figure 4 :

 4 Figure 4: Indicators of PNE proximity in randompayoff GG : square grid graph

Figure 5 :

 5 Figure 5: Indicators of PNE proximity in randompayoff GG: medium connectivity Erdös Renyi graph

  e-dominance. For any e = (eϕ, ev) ∈ R 2 + the e-dominance relation is defined on performance vectors x = (xϕ, xv), y = (yϕ, yv) ∈ F (A) by x e y ⇔ [(xϕ(1 + eϕ) ≥ yϕ) and (xv(1 + ev) ≥ yv)]. e-covering. For any e = (eϕ, ev) ∈ R 2 + and any set S ⊆ F (A), the set T ⊆ S is said to be an e-covering of S if ∀s ∈ S, ∃t ∈ T, t e s.

Figure 6 :

 6 Figure 6: e-Covering of Pareto-optimal tradeoffs


Riadh Dhaou
email: dhaou@enseeiht.fr

Benoît Escrig
email: escrig@enseeiht.fr

Béatrice Paillassa
email: paillassa@enseeiht.fr

Caroline Bès
email: caroline.bes@cnes.fr

Extending satellite service availability through energy efficient cooperation

Keywords: satellite communications, hybrid systems, emergency communications, cooperative protocols, relay selection

In this paper, we address the design of a cooperative protocol for a hybrid satellite/terrestrial emergency system. We want to perform energy savings compared to the case where all the terrestrial relay nodes are forwarding satellite messages to ground receivers. This is done via the selection of relevant relay nodes. The parameterization of the protocol phases has been done through simulations and takes into account the duration of the selection process, the number of selected nodes, and the signaling overhead. The selection process based on a node identifier (ID) appears to provide greater energy savings compared to the selection process based on the signal to interference and noise ratio (SINR). The solutions have been implemented in the real case scenario of forest fire that has been thoroughly documented by the US administration. According to the scenario parameters, 100% of the masked nodes are reached after cooperation.

I. INTRODUCTION

Satellite Communications (SCs) are known to provide services over large coverage areas without the need of a ground infrastructure [START_REF] Chini | A survey on mobile satellite systems[END_REF]. This benefit is paramount in case of a natural disaster (tsunami, earthquake) when rescue teams must be deployed over devastated areas [START_REF] Del Re | Salicesatellite-assisted localization and communication systems for emergency services[END_REF] [3] [START_REF] Iapichino | Advanced hybrid satellite and terrestrial system architecture for emergency mobile communications[END_REF]. However, the availability of a satellite data service is limited by the masking effect when there are obstacles between the satellite and the ground receivers, i.e., when the received signal to interference and noise ratio (SINR) is too low to achieve reliable decoding of the satellite message.

Diversity techniques such as time, frequency, and satellite diversity techniques [START_REF] Cocco | Cooperative satellite to land mobile gap-filler-less interactive system architecture[END_REF] can tackle this problem by providing ground receivers with different versions of a satellite signal. 1 In the context of hybrid satellite/terrestrial systems, cooperative communications can also enhance the service availability by allowing relay nodes to forward satellite messages to masked nodes [START_REF] Paillassa | Improving satellite services with cooperative communications[END_REF] [7] [START_REF] Kim | Evaluation of cooperative techniques for hybrid/integrated satellite systems[END_REF].

Cooperation schemes have already been implemented in broadcasting systems such as Digital Video Broadcasting Satellite Handheld (DVB-SH) by means of fixed relays, called gap-fillers. But these relays cannot handle dynamic situations, i.e., situations when ground receivers are mobile nodes [START_REF] Etsi | Digital video broadcasting[END_REF]. On 1 Increasing the number of versions reduces the probability that all the versions exhibit a low SINR at the same time.

the contrary, it has been shown that non masked (NoM) mobile nodes can serve as relay nodes in a dynamic environment. 2For instance, in [START_REF] Paillassa | Improving satellite services with cooperative communications[END_REF], satellite downlink (DL) communications are divided into two time-slots. During the first time-slot, the satellite is transmitting its message toward ground receivers. During the second time slot, the NoM nodes broadcast the satellite message. 3 Hence, the masked nodes are given a second opportunity to get the satellite message. No extra signaling is needed here to set the cooperation. No extra infrastructure, such as gap-fillers, is required. However, simplicity comes at the price of extra energy consumption. A masked node may receive two or more versions of the satellite signal even though only one version is needed. Moreover, an NoM node may broadcast the satellite message even if it is not in the range of a masked node.

The main contribution of the paper consists in designing a cooperative transmission scheme for hybrid satellite/terrestrial systems that both (i) increases the service availability, by tackling the masked node problem, and (ii) allows for some energy savings compared to the case when all NoM nodes serve as relay nodes, by selecting the best relay for each masked node. A second contribution consists in evaluating the proposed solutions on a real case. The context corresponds to a forest fire that has been investigated in [10].

The rest of the paper is organized as follows. The second section presents the wildfire scenario. The third section describes the protocol and the relay selection procedures. The fourth section addresses the simulation results and we conclude in section five.

II. WILDFIRE SCENARIO

A hybrid terrestrial/satellite system is implemented in the context of a wildfire. SCs are established between a remote operations center and emergency workers (firefighters, foresters) that are deployed on the fireground. The idea here is to propose a cooperation protocol between ground terminals so that masked nodes could receive the information from the operations center via NoM nodes. The terrain is so steep that no fire trucks with satellite transmission equipment can access the area.

The scenario is based on a real case that has been thoroughly documented by the U.S. Department of Agriculture [10]: a forest fire in a mountain area. In particular, the investigation showed that better communication between the fire control center and firefighters would have lead to a reduced number of deaths during the incident.

The DL traffic is constituted by images/maps transmissions (CBR 500 KB every 10 minutes). For sake of simplicity, the uplink (UL) traffic, and hence voice calls, will not be considered in this paper.

Twenty-one half-duplex nodes are deployed on the fireground (a first group of seven nodes and an additional group of fourteen nodes). Nodes are able to transmit on both a satellite link and a terrestrial WIFI link. Cooperation takes place in adhoc mode, i.e., there is no fixed relays. Emergency workers are deployed on a 400 hectares area (1000 meters x 4000 meters) and the deployment follows the fire line [10] according to three concentric and elliptic areas (Z 1 , Z 2 , and Z 3 ) depicted in Fig. 1. In the context of a forest environment, propagation models comprise two components: a plane earth path loss component and an additional component that varies according to studies, contexts, forest types, and frequency ranges. We have selected the Fitted ITU-R model developed in [START_REF] Al-Nuaimi | Measurements and prediction model optimization for signal attenuation in vegetation media at centimetre wave frequencies[END_REF] [START_REF] Meng | Empirical near ground path loss modeling in a forest at vhf and uhf bands[END_REF]. The multipath fading margin is neglected here but multipath effects are taken into account by using OFDM (Orthogonal Frequency Division Multiplex) waveforms.

III. PROTOCOL

The proposed protocol is called OHR for One Hop Relay. The time axis for the transmission of one frame is divided into two time-slots. The first time-slot (TS1) is for the satellite transmission. Data and control frames are transmitted during the second time-slot (TS2) in order to detect and select relay nodes, and to forward the information transmitted by the satellite (see Fig. 2).

During the detection step, all NoM nodes that have successfully decoded the satellite message, send a detection (DET) frame. The same frame is used for all nodes. The NoM nodes are considered as potential relays or relay candidates. Nodes that have not received the message from the satellite, but have successfully decoded the DET frame, detect themselves as being masked nodes. Considering OFDM signals on the terrestrial segment, DET frames from different NoM nodes are arriving at masked nodes within the OFDM guard interval, so they can be combined constructively. Basically, the During the selection step, each masked node aims at selecting its best relay among the NoM nodes in its range. Limiting the number of relays reduces considerably the energy consumption caused by data forwarding. When a masked node successfully decodes the DET frame, it sends a CFC frame (Call For Cooperation) to initiate the selection process. The CFC frames are transmitted after a random timout in order to minimize the risk of a collision. In case of a collision, the masked nodes wait for a random duration before trying to access the channel again. At that point, NoM nodes know that they have to help masked nodes. We now proceed to the competition step.

Potential relays are evaluated according to a given metric. In this paper, we compare two different metrics. The first metric is based on a node identifier (ID). The best relay corresponds to the node with the lowest address. The second metric is based on an SINR measurement. The best relay corresponds to the node that exhibits the highest SINR at the masked node. The metrics are collected by masked nodes through RTC messages (Ready To Cooperate) sent by the potential relays. Each relay candidate triggers a timer upon receiving the CFC frame. When the timer expires, an RTC frame is sent. So, when a node sends an RTC frame, it signals its presence to masked nodes and indicates its willingness to cooperate. Hence, masked nodes get a knowledge of their neighboring nodes and can select their best relay.

Note that the RTC frame of the NoM node N i contains its ID when the ID metric is used. When the SINR metric is considered, N i send the SINR values corresponding to each received CFC frame, i.e, the SINR values between itself and all the masked nodes in its range that sent a CFC frame. We can go a step further. An RTC frame of an NoM node N i could contain the list of all the masked nodes in the range of N i and their corresponding metrics. When the RTC frame from N i is overheard by another potential relay N j , this node could waive its participation to the selection process as soon as it discovers that N i has better metrics with respect to the masked nodes N j is about to help. This optimization has been considered and evaluated. It results that the duration of the cooperation set-up and the signaling overhead increased dramatically, thus precluding the implementation of such an approach. Due to the lack of space, the corresponding results will not be reported here.

We now address the way timers are set. Timer values depend on the metric. When the selection is based on the lowest ID (resp. highest SINR), timer values are proportional (resp. inversely proportional) to the metric [START_REF] Escrig | On-Demand Cooperation MAC Protocols with Optimal Diversity-Multiplexing Tradeoff[END_REF] [START_REF] Bletsas | A Simple Cooperative Diversity Method Based on Network Path Selection[END_REF]. So the best relay always signals its presence first. This contributes to minimize the duration of the distributed selection process. When IDs are used, timers can be configured so that there are no collision between different RTC frames. For instance, a timer τ i associated to node N i can be computed using the formula τ i = N k=1 a ik × τ k where the a ik denote the hexadecimal digits of the MAC address of N i and τ k denotes the time increment associated to the k th digit. When the SINRs are used, there is a risk of a collision since the metrics belong to the set of real numbers. In the case of a collision, NoM nodes access the channel later by triggering a timer with a random duration. Once the masked nodes have collected all the CFC frames, they select the best relay. The result of the selection step is transmitted in a CTS (Clear To Send) control frame that contains the identifier of the best relay. Here also, when a collision occurs, the masked nodes wait for a random duration before trying to access the channel again. In this paper, the values of the timers and the duration of the selection step are adjusted by simulation experiments based on the studied scenarios.

During the relaying step, NoM nodes are forwarding the satellite message. Since the signals are to be combined at the masked terminals, the forwarding is performed using the same modulation and coding scheme. There is no adaptive modulation and coding (AMC) mechanism. Hence, the energy savings provided by the SINR-based method are limited since it is not possible to transmit the message using a higher order modulation. This constraint limits the interest of selecting the best relay based on the SINR criterion. We assume that the relay terminals are synchronized, so they all transmit quasi-simultaneously and the forwarded signals can be added constructively at the receiver using a maximum ratio combiner. Terrestrial signals are not combined with the satellite one, since masked nodes did not receive a usable signal. When the node speed is low, the detection and selection phases can be omitted for subsequent frames (see Fig. 2).

Note also that the cooperative transmission scheme can be combined with a repetition scheme from the satellite. We consider that a single repetition scheme from the satellite is not enough to address any shadowing profile, any mobility model. The service availability is increased here through the help of surrounding nodes.

IV. PROTOCOL REFINEMENTS THROUGH SIMULATIONS

The choice for a simulation-based approach comes from the fact that there is no fully satisfactory mathematical solutions to our problem. For instance, minimizing the duration of the timer-based selection process leads to timers with smaller values. This increases the probability of a collision between relay candidates when random timers are used. That, in turns, increases the duration of the selection process [START_REF] Bletsas | A Simple Cooperative Diversity Method Based on Network Path Selection[END_REF]. So, instead of considering all possible options, we present here simulations results to attest the feasibility of the proposed solution.

The implementation of the proposed protocol is performed under the ns-2 network simulator. The purpose of the simulations is to refine the design of the protocols. Simulations aim to evaluate the duration of the protocol phases (detection, selection, and cooperation), the performance of the OHR protocol in terms of number of potential relays, number of selected relays according to the selection criterion: lowest ID or maximum SINR. The validation step of the protocol requires the development of specific modules in ns-2: a new propagation model (Fitted ITU-R) adapted to the forest environment, a technique for the coherent combining of multiple signals, two selection methods based on the lowest ID or the maximum SINR.

In the simulated scenarios, the terminals are placed on an ellipse. A given percentage of nodes is masked. The choice of masked nodes is different from one scenario to another: either the masked nodes are randomly selected, or an entire geographical area is shadowed. To synchronize the nodes and separate the protocol steps, each protocol phase has a specific duration. The values of these durations have been adjusted so the risk of collision is minimized during the cooperation setup and have been determined by ns-2 simulations. The results are depicted in Table I. These results have been obtained for a scenario in which the duration of the DL time-slot is 0.1 s, so the total duration of the three-phase protocol is also 0.1 s.

Several parameters have been evaluated in the simulations, including the duration of the cooperation set-up, i.e., the average time for a masked node to establish a cooperation, and the percentage of lost packets. These parameters have been estimated based on the density of the nodes, i.e., the total number of nodes in a specific geographical area, and according to the transmitted power. Each parameter has been evaluated taking into account the two relay selection methods, and the percentage of masked nodes, 25%, 50% and 75% of the total number of nodes, with a random spatial distribution.

A. Cooperation set-up time

First, we vary the node density and we plot the time required to establish cooperation for different percentages of masked nodes. The results are depicted in Fig. 3. Since the intervention zone has a fixed surface, the parameter "Density" in the figures depicts the number of nodes in the intervention zone. In Fig. 3, the set-up duration converges rapidly to 1.25 s whatever is the selection method (ID or SINR). So we choose to set the maximal duration for the selection step to 1.5 s. Moreover, the transmission power is set so that the node graph is a connected graph, i.e., there is at least one path between any two nodes. This minimizes the probability of a set-up fail.

B. Choice of the selection method

First, the overhead induced by the selection step is investigated (see Fig. 4 and Fig. 5). Here also, the parameter "Density" depicts the number of nodes in the intervention zone. Since there is no great difference between minimum/maximum values and the mean values, only mean value curves are plotted in the following figures. We note that the overhead decreases with the percentage of masked terminals. The total number of transmitted bits increases with the density of nodes, irrespective of the selection mode, ID or SINR. Second, the selection methods have been evaluated according to the number of selected relays. Fig. 6 shows the number We notice a clear difference between the lowest ID and the maximum SINR methods. When the number of nodes is increasing, it is always possible to find a best relay according to the SINR criterion. So the number of selected relays is increasing when the SINR criterion method is used. On the contrary, when the ID solution is used and the number of nodes is increasing, the probability to find a relay with the lowest ID approach increases with the density, leading the a decreased number of selected relays, and hence to a lower energy consumption.

Several other outcomes were observed, which are rather straightforward:

• When the number of masked nodes increases, the time required for all nodes to be able to cooperate increases.

• The SINR increases as the density of nodes increases.

• The percentage of potential relay decreases with the density of terminals.

• Conversely, the percentage of selected relay tends to decrease when the density of terminals increases.

More importantly, according to the forest fire scenario, all the masked nodes receive the satellite message after the cooperation step.

At first glance, we can say that the selection mechanism based on the SINR provides a better SINR on each link but the energy consumption at the network level is higher because there is potentially one single best relay for each masked node. On the contrary, the SINR on each link is not minimized when best relays are selected according to the lowest ID. However, the energy consumption is lower since one single NoM node can serve several masked nodes. We can go a step further using the following example.

Note first that, in all cases, all the relays are using the same waveform, i.e., the same modulation and coding scheme (MCS). The selected MCS is the one that provides the best robustness to channel impairments. We analyze the scenario presented in Fig. 7 where M 1 and M 2 are two masked nodes and R 1 and R 2 are two potential relays. We assume that M 1 and M 2 are in the range of both R 1 and R 2 . When the IDbased solution is considered, one relay can serve both M 1 and M 2 . Let R 2 be the selected relay. When R 1 is closer to M 1 than R 2 , one can claim that, in the case of the SINR-based selection, the transmission between R 1 and M 1 could use a little less energy than in the case of the ID-based solution. However, in the case of the SINR-based solution, we also must take into account that two transmissions will take place instead of one. So the energy saving that is earned on the R 1 -M 1 link is counterbalanced by the fact that two transmissions must take place, the one from R 1 to M 1 and the one from R 2 to M 2 . That leads us to select the ID-based solution. V. CONCLUSION SCs provide a relevant contribution to emergency systems since they allow the provision of data services over wide coverage areas. In this paper, we propose a cooperative hybrid satellite/terrestrial system in the context of a forest fire. The satellite segment is used in a first time-slot to broadcast emergency data. During a second-time slot, relay nodes are selected among the NoM nodes in the neighborhood of masked nodes. Then the selected nodes forward the satellite message they have received. The issue is to perform some energy savings compared to the case where all the NoM nodes forward the satellite messages to their neighboring nodes. For that purpose, two selection mechanisms have been proposed: one is based on the lowest ID and the other one is based on the maximum SINR.

From the simulation results, it appears that the ID-based selection process provides greater energy savings than the one based on the SINR. Even if the SINR is not maximized for all wireless links, this is largely counterbalanced by the fact that a relay can serve several masked nodes. The second contribution of the paper consists of applying the proposed protocol to the real case scenario of a forest fire. In the specific case that has been studied in this paper, the satellite messages could be delivered to 100% of the masked nodes after a cooperation step.

Further studies could address the case of multi-hop cooperation. Masked nodes, once they have received the satellite message from NoM nodes, could in their turn forward the message. Furthermore, fixed nodes have been considered in this paper. The case of mobile nodes should also be addressed to tackle more compelling issues and a special attention should be paid to mobility models. Lastly, synchronization at the signal level is also a major issue. This is left for future work.

Fig. 1 .

 1 Fig. 1. Fire spread diagram according to [10].

Fig. 2 .

 2 Fig. 2. Frame exchange sequence in the OHR protocol.

Fig. 3 .

 3 Fig. 3. Cooperation set-up duration for the OHR protocol. Mean values are represented by bold lines, minimum and maximum values by dotted lines. Labels ID and SINR denote the selection methods. The percentages represent the percentages of masked nodes.

Fig. 4 .

 4 Fig. 4. Mean values for the total number of emitted bits during the selection phase for the OHR protocol. Labels ID and SINR denote the selection methods. The percentages represent the percentages of masked nodes.

Fig. 5 .

 5 Fig. 5. Mean values for the total number of received bits during the selection phase for the OHR protocol. Labels ID and SINR denote the selection methods. The percentages represent the percentages of masked nodes.

Fig. 6 .

 6 Fig.6. Mean values for the number of selected relays as a percentage of detected relays for the OHR protocol. Labels ID and SINR denote the selection methods. The percentages represent the percentages of masked nodes.

Fig. 7 .

 7 Fig. 7. Relay transmissions in the case of the ID-based selection (dotted line) and the SINR-based selection (solid line).

NoM nodes are the nodes that have successfully decoded the satellite message.

We assume that not all the ground receivers are masked. This assumption makes sense since the sizes of obstacles are small compared to the sizes of the coverage areas.


LOSSY COMPRESSION OF UNORDERED TREES

ROMAIN AZAÏS, JEAN-BAPTISTE DURAND AND CHRISTOPHE GODIN OBJECTIVES DAG compression is a classical technique that consists in building a Directed Acyclic Graph (DAG) that represents an unordered tree without the redundancy of its subtrees. Only trees with a high level of redundancy are efficiently compressed by this method. We introduce a lossy compression method that consists in computing the DAG of a self-nested structure that both presents the highest level of redundancy and approximates the initial data to compress.

EDIT DISTANCE

This study requires to introduce a distance on the space of unordered trees in order to quantify the quality of a given self-nested estimate of a tree structure. We consider an editing distance defined from the following tree edit operations: insertion and deletion of a leaf. An edit script is a sequence of edit operations. The result of applying an edit script s to a tree τ is the tree τ s obtained by applying the component edit operations to τ in the order they appear in the script.

The cost of an edit script is only the number of edit operations. Finally, given two trees τ 1 and τ 2 , the edit distance δ(τ 1 , τ 2 ) is the length of the minimum edit script that transforms τ 1 into a tree that is isomorphic to τ 2 .

EXAMPLE

The tree data to compress and its DAG version (top) and the solutions provided by the lossy compression algorithms NEST (middle) and LP-RFC (bottom). The NEST solution has too many nodes for looking like the initial tree (107 nodes added to obtain a self-nested tree). The LP-RFC solution is visually close to the initial tree, which is confirmed by the error rate of 35%. Both algorithms achieve a compression rate greater than 80%.

[1] Azaïs, Durand, and Godin. Lossy compression of trees via linear DAGs. Preprint, 2016.

[2] Godin and Ferraro. Quantifying the degree of selfnestedness of trees. Application to the structural analysis of plants. IEEE TCBB, 2010.

WORST-CASE APPROXIMATION ERROR

Among trees of height h and outdegree less than m, the tree that is the farthest to a self-nested tree may be identified. The editing distance to its best self-nested approximation is of order 0.25 m h . The diameter of the space of unordered trees being of order m h , it follows that the largest area without any self-nested tree is of relative radius 0.25. In addition, the error rate for this tree (and thus the worst error rate that must be expected from any lossy compression algorithm) is of order 0.33.

SELF-NESTED TREES

Self-nested trees present remarkable compression properties because of the systematic repetition of subtrees. They are defined as trees such that all the subtrees of a given height are isomorphic.

The DAG related to a self-nested tree τ is linear (there exists a path going through all its vertices) and has thus height(τ ) + 1 nodes. In other words, self-nested trees achieve maximum compression rates. In addition, self-nested trees are very rare, while still being quite close to any unordered tree.

ALGORITHM

The preexisting algorithm (NEST) only adds nodes to a tree τ to transform it into a selfnested structure. Instead of adding nodes, we propose to replace some internal structures of τ by their centroid (RFC). In particular, this allows us to delete some nodes and thus gives flexibility. These self-nested estimates may be computed in polynomial time. Nevertheless, this procedure can only modify subtrees and not delete them. That is why we introduce a new algorithm that exploits local pruning of τ (LP-RFC).

Figure 1 :

 1 Figure 1: Error and compression rates on the example.

Figure 3 :

 3 Figure 3: Compression rates on the simulated dataset.

Table 1 :

 1 Frequency of self-nested trees with given maximal height and ramification number with respect to unordered trees under the same constraint.

		deg ≤ 2	deg ≤ 3	deg ≤ 4
	h ≤ 2	0.88	6.18.10 -1	3.52.10 -1
	h ≤ 3	0.49	3.38.10 -2	7.43.10 -5
	h ≤ 4	0.07	2.90.10 -8	4.16.10 -23
	h ≤ 5 3.36.10 -4 3.56.10 -28 1.66.10 -100

CONTACT INFORMATION

Romain Azaïs

Inria team BIGS -Institut Élie Cartan de Lorraine Université de Lorraine in Nancy, France romain.azais@inria.fr

SIMULATION STUDY

Our simulations are performed on 500 small binary trees generated from a stochastic model. The three algorithms are equivalent in terms of compression rates. The key parameter is thus the error rate that is much better for RFC and LP-RFC algorithms than for NEST procedure (substantial gain of 20%). Local pruning is useful in RFC 24.2% of the time and makes the error decrease of 1.8% (7.3% when useful).


F Ladieu

C Brun

D L'hôte

Nonlinear dielectric susceptibilities in supercooled liquids : a toy model

The dielectric response of supercooled liquids is phenomenologically modeled by a set of Asymmetric Double Wells (ADW), where each ADW contains a dynamical heterogeneity of Ncorr molecules. We find that the linear macroscopic susceptibility χ1 does not depend on Ncorr contrary to all higher order susceptibilities χ 2k+1 . We show that χ 2k+1 is proportional to the k th moment of Ncorr, which could pave the way for new experiments on glass transition. In particular, as predicted by Bouchaud and Biroli on general grounds [Phys. Rev. B, 72, 064204 (2005)], we find that χ3 is proportional to the average value of Ncorr. We fully calculate χ3 and, with plausible values of few parameters our model accounts for the salient features of the experimental behavior of χ3 of supercooled glycerol.

Upon fast enough cooling, most liquids do not cristallize but enter into a supercooled liquid state [1][START_REF] Donth | The Glass Transition[END_REF][3][4][5], where the viscosity η dramatically increases with lowering the temperature T . Below the glass transition temperature T g , η is so high that the system is in practice a solid -the glass-, yet, no structural difference between the glass and the liquid state has ever been detected [1]. Over the past fifteen years, a major breakthrough was the discovery of Dynamical Heterogeneities (D.H.) in supercooled liquids [6][7][8][9][10][11] ; i.e., relaxation happens through collective events gathering N corr molecules, and some groups are relaxing much faster than others. As it is expected that an increase of N corr when lowering T could increase dramatically η, an significant effort was made for measuring the T -dependence of N corr [8,12,13].

It has been argued that the most direct way to draw accurately the T -dependence of N corr from experimental data is based on the a.c. nonlinear susceptibility χ 3 [14][15][16], where χ 3 is the third order response of the fluid to a field with an angular frequency ω. This field can be of any nature, e.g. electric as in [15,16]. More precisely two nonlinear susceptibilities are related to χ 3 : χ 3 , which correspond to the third order nonlinear response at the third harmonics (i.e., at "3ω") and at the first harmonics (i.e., at "1ω") respectively. Bouchaud and Biroli (BB) have shown [14,16] that χ

Here k B is the Boltzman constant, and τ α (T ) is the typical relaxation time at temperature T corresponding to the relaxation frequency f α = 1/(2πτ α ) where the imaginary part of the linear response is maximum.

∆χ 1 = χ 1 (ω = 0) -χ 1 (ω → ∞)
is the part of the static linear susceptibility corresponding to the slow relaxation process of interest, a 3 is the volume occupied by one molecule, and H and K are two complex scaling functions that approach zero for both small and large ωτ α . Note that the humped shapes of |H (ωτ α ) | and |K (ωτ α ) | are distinctive features of the glassy correlations. BB's prediction relies on very general grounds, such as a generalised fluctuation dissipation relation, and was inspired by spin glass physics [17], where a true second order phase transition happens at T c , accompanied by a critical divergence of χ 3 (while the linear suceptibility χ 1 does not diverge). A consequence of this generality is that the detailed expressions of the scaling functions H and K remain unkown. Here we present a phenomenological "toy" model where, for the first time, BB's predictions are recovered with an explicit expression for the functions H and K. By using plausible values of free parameters, the most salient experimental features of [15,16] can be accounted for. Moreover we obtain new predictions on higher order nonlinear susceptibilities χ 2k+1≥5 . This could motivate new experiments deepening our understanding of the glass transition.

Model : We assume that all D.H.'s are independent from each other and that a given D.H. is a group of N corr molecules evolving in an Asymmetric Double Well potential (ADW), depicted in Fig. 1. Each ADW is characterised by the height of its barrier V and by an asymmetry energy ∆. We neglect internal field effects. On Fig. 1,z represents the axis of the external electric field E(t) = E cos(ωt), and θ 1 is the angle between the field and the well which has the deepest energy at E = 0. For simplicity we assume that θ 2 = θ 1 + π. With respect to earlier versions [19,[START_REF] Böhmer | Chapter 14 in Broadband Dielectric Spectroscopy[END_REF], a key refinement is the assumption that the magnitude of the net dipolar moment µ, in either of the two wells, is given by µ = µ molec √ N corr where µ molec is the molecular moment : This estimator of µ is assumed here because there should not exist any geometrical ordering among the molecules contributing to a given D.H. [6]. With v DH = N corr a 3 the volume of a D.H., the simplest approach, for θ 1 = 0 and ∆ = 0, yields a static polarisation given by (µ/v DH ) tanh(µE/k B T ).

Expanding in E gives ∆χ 1 ∝ µ 2 /v DH , which is independent of N corr since the N corr dependence of µ 2 cancels that of v DH . For all higher orders such a cancellation does not happen, e.g. χ 3 ∝ µ 4 /v DH ∝ N corr . This is the main reason why we find below that χ 1 is blind to N corr contrarily to all higher order susceptibilities.

Let us now consider a set of N identical ADW's. With Π j,k the transition rate from the well k to the well j, we obtain the number n 1 -resp. where e(t)

≡ F cos(ωt) , F = µ molec √ N corr cos(θ 1 ) k B T E and M = µ molec cos θ 1 √ N corr a 3 , δ = tanh( ∆ 2k B T
).

(

) 2
Setting E = 0 in Eq. 2 yields P = P 0 = Mδ. As M ∼ cos θ 1 , we obtain < P 0 >= 0 where the brackets denote the average over the isotropically distributed values of θ 1 . In the limit of small fields (i.e., e → 0), expanding Eq. 2 to the first order in e yields :

< P 1 (t) >= < MF > (1 -δ 2 ) 1 + (ωτ ) 2 cos(ωt -arctan ωτ ), (3)
i.e., a Debye response, as expected in any double well model [START_REF] Frölich | Theory of Dielectrics[END_REF]. As the linear dielectric spectra of supercooled liquids are asymmetric in frequency, we assume, as in other phenomenological models [11,[START_REF] Brun | [END_REF], that the values of τ are distributed according to G(τ ), [23]. G(τ ), given in [18] and in Fig. 1, is chosen to recover accurately the experimentally well known linear susceptibility χ 1 (ω, T ) by weighting Eq. 3 with G(τ ) and summing over all values of τ . More precisely, G(τ ) determines the shape of χ 1 (ω, T ), but not its overall magnitude ∆χ 1 . We use the experimentally well known value of ∆χ 1 as an additional constraint in our model : from Eq. 3, we obtain µ 2 molec = 3k B T ǫ 0 a 3 ∆χ 1 /(1 -δ 2 ) ; i.e., µ molec is no longer a free parameter.

Computing χ

3 and χ

(1)

3 : For a given value of N corr and of δ, we consider all ADW's having the same τ and θ 1 . By using Eq. 2, we compute the polarisation up to the third order in field. We first average the result over θ 1 , then sum over τ with weight G(τ ), and finally average over the values of N corr existing among various D.H.'s. The latter average is denoted by [ ] av . This yields quantities -labelled below by an index "m" standing for "model"-which are comparable to experiments [15,16]. Note that this method, by using the values of a 3 , ∆χ 1 , G(τ ) drawn from standard experiments, eliminates θ 1 , ν 0 , V and µ molec . Thus, when comparing our model to the experimental values of χ 3 at a given T , the two remaining free parameters are [N corr ] av and δ. For simplicity we take a single value for δ, and postpone the possible averaging over δ to Ref. [18].

In practice, we solve Eq. 2 by assuming e ≪ 1, and develop Eq. 2, as well as P , in series of e, up to the third order. As the polarisation of a given set of ADW's sharing the same θ 1 , τ, N corr is not symmetric with respect to field reversal E → -E, we set P (t) = q=3 q=0 P q (t) where P q ∝ E q . Since e ≪ 1, one has |P q | ≫ |P q ′ >q | ; i.e., all P q ′ >q can be neglected when looking for P q . Thus, P q is obtained by keeping only the terms ∝ e q in Eq. 2. This was illustrated above to get first P 0 and then P 1 (t) -see Eq. 3-. Repeating the procedure to the order e 2 yields P 2 , [18]. Finally going to the order e 3 gives :

τ d(P 3 ) dt + P 3 = M(1 -δ 2 ) 6 e 3 - 1 2 P 1 e 2 -δP 2 e. (4)
As P 1 and P 2 are known, the analytical expression of P 3 (t) is readily obtained from Eq. 4. After averaging over θ 1 , τ, N corr , one obtains P 3,m that must be identified with the third order term P 3 of the experimental polarisation. As e 3 ∝ E 3 (3/4 cos(ωt) + 1/4 cos(3ωt)), we recall that P 3 naturally defines the first and third harmonics cubic susceptibilities (with phases -δ (k) 3 , k = 1, 3) as [24] :

P 3 (t) ǫ 0 = 3E 3 4 |χ (1) 3 | cos(ωt -δ (1)
3 ) +

E 3 4 |χ (3) 3 | cos(3ωt -δ (3)
3 ).

(

) 5
Results of the ADW model : A dimensional analysis shows that our model yields P q,m ∝ [< MF q >] av , which has two important consequences. First it yields P q,m ∝< (cos θ 1 ) q+1 >. This implies that the even terms P 2k,m ≡ 0, which ensures that the macroscopic polarisation reverses exactly upon the E(t) → -E(t) reversal, as required by macroscopic symmetry considerations. Second, all odd terms P 2k+1,m are non zero, yielding for the susceptibilities : χ 2k+1,m ∝ N k corr av . This shows that the linear susceptibility χ 1,m is blind to the value of N corr , contrary to higher order susceptibilities which are directly proportionnal to the k th moment of N corr . This first important result is reminiscent of the spin-glass transition [17] which has inspired BB's prediction.

The above mentionned analysis yields χ

3,m that we convert into its dimensionless form X

(3) 3,m = χ (3) 3,m k B T /[ǫ 0 a 3 (∆χ 1 ) 2 ]. Writing X (3) 3,m = |X (3) 3,m | exp[-iδ (3) 3,
m ] where i 2 = -1, we get finally [18], with x = ωτ :

3,tot is the weighted sum (see text) of X

3,m (see Eq. 6) and of X

3,trivial corresponding to the cubic susceptibility of independant molecules undergoing rotational brownian motion [16,25]. For f /fα ≥ 1 one has X is displayed here, and that [Ncorr]av = 15. δ and fex have the same values as in Fig. 3.

X (3) 3,m = 9[N corr ] av 5(1 -δ 2 ) ∞ 0 G(τ ) D (3) 3 (x)e i Ψ (3) 3 (x)-arctan(3x) 1 + (3x) 2 dτ D (3) 3 (x)e iΨ (3) 3 (x) ≡ 1 6 - e -i arctan x √ 1 + x 2 1 2 -δ 2 e -i arctan(2x) √ 1 + 4x 2 . ( 6
)
Note that G nearly obeys Time-Temperature Superposition (TTS) ; i.e., it is nearly [START_REF] Brun | [END_REF] independent on T when plotted as a function of λ = τ /τ α . As ωτ = λωτ α , Eq. 6 shows that X

(3) 3,m equals [N corr ] av times a function which does not depend on T -we take δ as a constant in T -, when plotted as a function of ωτ α . Thus Eq. 6 gives the first phenomenological expression of the function H(ωτ α ) of Eq. 1 -we recall that according to BB's prediction X

(3) 3 is [N corr ] av H-. Eq. 6 thus shows explicitly that the T dependence of X

(3) 3 is directly that of [N corr ] av , up to small effects coming from small violations of TTS in G(τ ). This is the second important result of our model. Fig. 2 shows the frequency behavior of |X

3,m |. For most values of δ, the spectrum has a low pass character. In the vicinity of δ ⋆ = 1/ √ 3 the spectrum has a humped shape. To understand this, let us note P stat the solution of Eq. 2 at ω = 0. One gets P stat = M tanh[e + ∆/(2k B T )]. Expansion to order e 3 yields X stat 3,m (δ ⋆ ) = 0. Around δ ⋆ , X stat 3,m moves from a negative "Ising-like" value (low δ's), to a positive value for very asymmetric ADW's (high δ's). When ω = 0, the effective relaxation time τ /(cosh e + δ sinh e) comes into play, which contributes also to X

(3) 3,m . This is why close to δ ⋆ , |X

3,m | has a humped shape in frequency. A deeper, i.e., much less model dependent, reason for this humped shape is given below.

To compare our model to the nonlinear susceptibilities of glycerol reported in [15,16], we first focus on the case f ≥ f α . Fig. 3 shows that choosing [N corr ] av = 5 and δ = 0.60 yields a very good agreement between our model and the values of X

(3) 3 (f ≥ f α ) measured at T = 204.7K ≃ T g + 16K. We emphasize that the agreement is good for both the modulus and the phase of X

(3) 3,m . Fig. 4 shows the same kind of comparison for X (1) 3 , for which an expression similar to Eq. 6 is given in [18]. On Fig. 4 the best agreement between our model and the data reported in [16] is obtained with [N corr ] av = 15 and δ = 0.60 : With respect to the data, our model underestimates the phase by ≃ 20 • and yields a maximum for the modulus at f ⋆ ≃ 1.6f α not far from the experimental value of 2.5f α . The fact that the optimal [N corr ] av is not the same in Fig. 3 and in Fig. 4 may come from interferences between the nonlinear responses of the D.H.'s with different τ , see Eq. 6 and Ref. [18]. These interferences have different effects on X 3,m , see [18], and this is not fully captured by our toy model, due to its simplicity. We emphasize, on the other hand, that [N corr ] av = 5 -15 is the right order of magnitude when comparing to the values given by 4D-NMR experiments [8] or by Ref. [13]. Moreover our model accounts for the fact that |X 3 | in glycerol. Finally, δ = 0.60 amounts to ∆ ≃ 1.4k B T ≈ 1.4k B T g , i.e. it does not introduce a new energy scale. Now, let us move to the case f < f α . Here we must take into account the finite lifetime τ ex of D.H.'s ; i.e., the fact that the liquid flows at large times [10]. The effective value of [N corr ] av decreases with frequency when f τ ex ≤ 1, since a given molecule is involved in various DH's at large times ; i.e., it becomes independent of other molecules in the long run : as X (k=1,3) 3,m ∝ [N corr ] av , this will give a humped shape to the nonlinear susceptibility even for the values of δ where |X (k) 3,m | has a low pass character. To take this idea into account, we simply use the -well knownnonlinear response X (k) 3,trivial of independent molecules (see [15,16,25]) and assume that it dominates the measured X (k) 3 when f τ ex ≪ 1. In practice, we write heuristically the total cubic susceptibility X (k) 3,tot (with once again k = 1, 3) as :

X (k) 3,tot = pX (k) 3,m + (1 -p)X (k) 3,trivial with p = exp (-f ex /f ), see [18]. For f > f α , X (k) 3,tot is of course very close to X (k) 3,m , since p ≃ 1. For f ≤ f α , Figs. 3-4 show that, with f ex = 0.14f α , X (k)
3,tot has the same global qualitative trends as the measured X (k) 3 in glycerol. We note that f ex /f α = 0.14 amounts to Q = τ ex /τ α ≃ 7, which is compatible with the values Q ≃ 3 -10 reported before [10] albeit still debated [18]. We think that the oscillation of |X (k) 3,tot | around 0.1f α is unphysical and comes from the very naive way of including τ ex in our analysis.

To conclude, we have developped a very simple toy model for the nonlinear susceptibilities in supercooled liquids. We find that χ 2k+1,m ∝ [N k corr ] av ; i.e., that χ 1,m is blind to the value of N corr contrary to all higher order susceptibilities. This yields the first phenomenological expression of the scaling functions involved in BB's predictions. With reasonnable values of parameters, the main trends of nonlinear experimental data are recovered. Our model explains very simply why the nonlinear responses yield brand new information on the glassy dynamics. This simplicity may trigger more experiments deepening our understanding of the glass transition. where e(t) ≡ F cos(ωt) ,

F = µ molec √ N corr cos(θ 1 ) k B T E and M = µ molec cos θ 1 √ N corr a 3 , δ = tanh( ∆ 2k B T ). (8)
As explained in the article, the two sources of nonlinearity in Eq. 8 are : (i) the nonlinear character of the equilibrium value P stat = M(δ cosh e + sinh e)/(δ sinh e + cosh e) = M tanh[e + ∆/(2k B T )] ; and (ii) the nonlinear character of the instantaneous relaxation time τ ef f = τ /(δ sinh e + cosh e).

We expand P (t) in series of powers of the field E up to third order P (t) = P 0 + P 1 (t) + P 2 (t) + P 3 (t) where P q ∝ E q . As E(t) 2 = E 2 (1 + cos(ωt))/2 and E(t) 3 = E 3 (3 cos(ωt) + cos(3ωt))/4, P 2 (t) and P 3 are the sum of two terms :

P 2 (t) = P (0) 2 + P (2) 2 (t) P 3 (t) = P (1) 3 (t) + P (3) 3 (t), (9)
where the superscript in parentheses indicates the index of the relevant harmonics. For example, P 3 (t) is given by a term oscillating at the fundamental frequency, and by a term oscillating at three times the fundamental frequency.

As the condition e ≪ is well obeyed experimentally, one can neglect all P q ′ >q terms when computing P q . Therefore P q is obtained by keeping only the terms ∝ e q in Eq. 8 above.

To the order e 0 it is found that :

P 0 = Mδ. (10)
Now, going to the order e 1 , one has (by using the result for P 0 in Eq. 10) :

τ dP 1 dt + P 1 = M 1 -δ 2 × F cos(ωt), (11)
which yields :

P 1 (t) = M(1 -δ 2 ) 1 + (ωτ ) 2 F cos (ωt -arctan(ωτ )) . (12)
We now go to the order e 2 and get :

τ d(P 2 ) dt + P 2 = -δF P 1 (t) cos(ωt). (13)
As P 1 (t) oscillates at frequency ω, the right hand side of Eq. 13 contains one constant term and another term oscillating at 2ω. Therefore, one finds :

P (0) 2 = M(δ -δ 3 )F 2 2 1 + (ωτ ) 2 cos [π + arctan(ωτ )] P (2) 2 (t) = M(δ -δ 3 )F 2 2 1 + (ωτ ) 2 1 + (2ωτ ) 2 cos [2ωt + π -arctan(ωτ ) -arctan(2ωτ )] . (14)
Finally, we reach the order e 3 and get :

B. Averaging over θ1,τ and Ncorr.

As explained in the main article, three kinds of averages must be done in our ADWP model : (i) First, we have to average over the angle θ 1 , the values of which are assumed to be isotropically distributed. Denoting this average by < >, one finds < (cosθ 1 ) 2k+1 >= 0 and < (cosθ 1 ) 2k >= 1/(2k + 1), for any integer k. As we have found above that P q ∝ MF q , one obtains < P q >∝< MF q >∝< (cosθ 1 ) q+1 >. Therefore, all the even integer harmonics vanish, contrarily to all odd harmonics which are found to be

< P 2k+1 >∝ N k corr . (22)
(ii) Second we have to average over various relaxation times τ , with weight G(τ )dτ . The distribution function G is chosen so as to recover accurately the experimental linear response, χ 1 . Therefore G must simultaneously solve the two following equations for the real part, χ ′ 1 , and the imaginary part, χ ′′ 1 :

χ ′ 1 (ω) -χ ′ 1 (∞) ∆χ 1 = ∞ -∞ G(ln τ ) × 1 1 + (ωτ ) 2 × d ln τ χ ′′ 1 (ω) ∆χ 1 = ∞ -∞ G(ln τ ) × ωτ 1 + (ωτ ) 2 × d ln τ, (23)
where we have used the fact that G(τ )dτ = G(ln τ )d ln τ . In practice one uses [START_REF] Brun | [END_REF]23] :

G(ln(τ )) = N GGE e -( β
Here Γ(x) is the Euler gamma function and α, β, σ, γ, τ 0 are T dependent parameters. For glycerol, a good set of parameters is given by :

with T expressed in Kelvins. Note that τ 0 is nearly proportionnal to the typical relaxation time τ α defined by τ α = 1/(2πf α ) where f α is the frequency of the peak of χ ′′ 1 . Additionally one finds from Eq. 23 :

∆χ 1 = (1 -δ 2 )(µ molec ) 2 < [cos(θ 1 )] 2 > k B T ǫ 0 a 3 , (26)
and with < (cos θ 1 ) 2 >= 1/3, we obtain 3,m as those obtained in the main article with a single value of δ. The line is a guide to the eyes.

1 ) 3

 13 related to the average value of N corr over the various D.H.'s existing at a given T -noted [N corr (T )] av -by : χ (3) 3 (ω, T ) ≈ ǫ 0 (∆χ 1 ) 2 a 3 k B T [N corr (T )] av H (ωτ α ) χ ((ω, T ) ≈ ǫ 0 (∆χ 1 ) 2 a 3 k B T [N corr (T )] av K (ωτ α ) .

Figure 1 :

 1 Figure 1: Inset : ADW model where each D.H. of Ncorr molecules evolves in an asymmetric double well making an angle θ1 with respect to the applied field E. Main graph : distribution of relaxation times for Glycerol at T = 204.7K, [18, 23].

  n 2 -of ADW's in state 1 -resp. state 2-by solving the two master equations : ∂n 1 /∂t = -Π 21 n 1 + Π 12 n 2 and ∂n 2 /∂t = -Π 12 n 2 + Π 21 n 1 . Assuming thermally activated barrier hoppings, one gets [19, 20] : Π 12 = W exp[(∆/2 + µE cos θ 1 )/k B T ], Π 21 = W exp[-(∆/2 + µE cos θ 1 )/k B T ] where W = ν 0 exp[-V /k B T ]. Here ν 0 = 1/τ m where τ m is the microscopic characteristic time of the thermal fluctuations within each well. The polarisation P of the set of N identical ADW's is given by P = µ cos(θ 1 )(n 1 -n 2 )/(N v DH ). The two master equations yield the dynamical equation for P , which involves the relaxation time τ = 2W cosh(∆/2k B T ) of the identical ADW's : τ dP dt + P (δ sinh e + cosh e) = M (δ cosh e + sinh e)

Figure 2 :

 2 Figure 2: (Color Online) Effect of the dimensionless asymmetry energy δ on |X (3) 3,m | : the spectra have a low pass character excepted close to δ ⋆ = 1/ √ 3. Here [Ncorr]av = 5.

Figure 3 :

 3 Figure 3: (Color Online) For [Ncorr]av = 5, δ = 0.60 comparison of the ADW model with the experiments of Ref. [16] at T = 204.7K. X

  and the experiments are very well accounted for by the model. For f /fα ≤ 1, only the global trends of the data are restored by the model with fex/fα = 0.14. Inset : Phases corresponding to the main graph, same symbols.

Figure 4 :

 4 Figure 4: (Color Online) Same symbols as in Fig. 3, excepted that X (1) 3

( 3 ) 3 |

 33 is peaked at a frequency ten times smaller than |X(1)

  δ sinh e + cosh e) = M (δ cosh e + sinh e)

  α = 10 β = -5.5996 × 10 -1 + 4.0900 × 10 -3 T + 1.50795 × 10 -5 T 2 σ = 1.57 × 10 -1 exp 407.525 T -141 γ = -7.826920 + 1.015 × 10 -1 T -4.32345 × 10 -4 × T 2 + 6.34415 × 10 -7 × T 3 τ 0 = 1.1511 × 10 -15 × exp 19.08905 × 127.38588 T -127.38588 ,

(µ molec ) 2 Fig. S 5 :

 25 Fig. S 5: Values of δ1, δ2 yielding, with a gaussian weight distribution, the same values of X (k)

We thank R. Tourbot for his outstanding technical help, S. Nakamae for carefully reading the paper. We thank G. Diezemann for encouraging discussions in dec. 2011, and C.Alba-Simionesco, G. Biroli, J.-P. Bouchaud, J.-P. Carton, P.M. Déjardin for long lasting help.

Note added : See also on the same subject the paper of Gregor Diezemann to appear on Condmatt TODAY plus or minus a few days !...

Electronic Physics Auxiliary Publication Service : Supplementary Information for :

Nonlinear dielectric susceptibilities in supercooled liquids : a toy model.

F. Ladieu, C. Brun, D. L'Hôte

We detail hereafter the calculations for X 3,m summarized in our main article. We then give a short justification of our assumption p = exp[-f ex /f ] made in the end of the main article. Finally, we give more informations about what happens when averaging over the dimensionless asymmetry parameter δ.

I. THE NONLINEAR SUSCEPTIBILITIES IN THE ASYMMETRIC DOUBLE WELL POTENTIAL MODEL.

A. Calculations for one set of identical Asymmetric Double Wells.

In this section we consider a set of N identical Asymmetric Double Wells (ADW) ; i.e., a set of ADW's sharing the same values for all microscopic parameters of the model. Denoting n 1 (respectively n 2 ) the number of ADW's in state θ 1 (respectively θ 2 = θ 1 + π), the polarisation P of the considered set of ADW's is given by :

where it was assumed that the net dipolar moment in either of the two states of a given ADW is given by µ = µ molec √ N corr , with µ molec the molecular dipole moment (see the main article). Combining the two master equations for n 1 (t) and n 2 (t), with n 2 (t) = N -n 1 (t), one gets the equation for the dynamics of P :

We separate the terms oscillating at ω from those oscillating at 3ω. Denoting by

2 |-the amplitude of P 1 (t) -respectively P (2) 2 (t)-, one obtains :

3 ) dt + P

(1) 3

as well as

3 ) dt + P

(3) 3

By using Eqs. 12-14 and the two previous equations, one finds :

3 ) dt + P

3 (ωτ ) cos ωt + Ψ

where D

3 (ωτ ) cos(ωt + Ψ

3 (ωτ )) ≡

as well as

3 ) dt + P

3 (ωτ ) cos 3ωt + Ψ

with D

3 (ωτ ) cos(3ωt + Ψ

(3)

Note that the above definitions of D

3 and Ψ

3 are consistent with those in the main article. The solution of Eq. 18 is given by

The solution of Eq. 19 is given by

At this point the two free parameters of our ADWP model are N corr and δ. In particular Eq. 27 sets the value of µ molec , since ∆χ 1 and a 3 are experimentally well known.

(ii) The third and the last average to be taken is over the values of N corr . Indeed the proportionality expressed by Eq. 22 remains true when averaging over the τ 's. Therefore, it is very easy in our model to take into account the fact that there exists a distribution of the values of N corr among various dynamical heterogeneities of a real supercooled liquid. As all above equations have been derived for given free parameters N corr and δ, we superpose the ensemble of models with the same δ but with different values of N corr . Denoting the average over N corr by [ ] av , we obtain from Eq. 22 :

where χ 2k+1,m is the macroscopic nonlinear susceptibility of the order 2k+1, and where the index m stands for "model", so as to avoid any confusion between the nonlinear susceptibilities produced by the model and those corresponding to what is experimentally measured (denoted χ 2k+1 ). Note that χ 2k+1,m represents generically the set of components of the macroscopic polarisation which is proportionnal to E 2k+1 and oscillates at one of the odd harmonics between 1ω and (2k + 1)ω. For example, χ 3,m corresponds to two terms : one is proportionnal to χ

3,m -note the presence of exponent (1)-and oscillates at 1ω, and the other one is proportionnal to χ The macroscopic polarisation P is given by [24] :

where the function χ 1 (t) corresponds to the experimental macroscopic linear response while χ 3 (t 1 , t 2 , t 3 ) is the experimental macroscopic nonlinear response. It is shown in ref. [24], that for a field E(t) = E cos(ωt) one gets :

3 ) + 1/4E 3 χ

3 ) + ...

We now must identify the result of our model with above relations giving the experimental macroscopic polarisation. We start from Eqs. 20-21 and average over θ 1 which yields :

We then average over the τ 's, as in Eq. 23, and then over N corr . With Eqs. 21-30, we obtain :

where the last equality was obtained by replacing χ

As in the main article one defines the dimensionless nonlinear suscceptibility as X

3,m = χ

3,m ] one obtains :

3,m ≡ phase of X

3,m = arctan S

(3) SIN

S

(3) COS with S

(3)

A similar calculation for χ

3,m yields :

Note that in the first equality of Eq. 35 there is a factor 3 instead of 9 found in Eq. 34. This comes from Eq. 30 where there is a factor 3/4 for the cubic term oscillating at ω while it is only 1/4 for the cubic term oscillating at 3ω.

II. MORE ON THE FREQUENCY DEPENDENCE OF THE WEIGHT p = exp[-fex/f ].

All above calculations have been made as if the lifetime τ ex of the considered Asymmetric Double Wells is infinite. As a supercooled liquid is ergodic above T g , the heterogeneity of the dynamics implies that τ ex must be finite. This comes from the fact that a region of space relaxing faster than the average must become a region relaxing slower than the average, to restore ergodicity. We shall assume, for simplicity, that any ADW is reshuffled with the same characteristic time τ ex , whatever the value of τ it had just before.

After reshuffling, the glassy correlations are different from those established before. Thus, if one performs an average over time longer than τ ex , a given molecule is no longer correlated to any other molecule. This is why, one expects any molecule to become effectively independent of all other molecules in the limit of large times t ≫ τ ex . Therefore one expects, at large times, the measured nonlinear dimensionless susceptibilities X (k) 3 to be dominated by the corresponding susceptibilities X (k) 3,trivial of independent molecules undergoing Brownian rotational motion. Note that X (k) 3,trivial has been fully calculated in Ref. [25]. Very few things are quantitatively established concerning the reshuffling phenomenon ; even the value of Q = τ ex /τ α remains a subject of discussions [START_REF] Richert | Experimental approaches to Heterogeneous Dynamics[END_REF]. Therefore a detailed description of its impact on nonlinear susceptibilities is not available at present. This is why, we heuristically add the nonlinear susceptiblities X (k) 3,m given by our ADWP model (multiplied by the weight p) to X (k) 3,trivial (multiplied by the complementary weight (1 -p)). As the limit of large times t ≫ τ ex corresponds to low frequencies f τ ex ≤ 1, we physically expect that p vanishes in this limit. The simplest way to express this idea quantitatively is to state that the weight (1-p) of the trivial response is given by the probability of a reshuffling event happening during one E oscillation period of 2π/ω. It is reasonnable to assume that the probability of the reshuffling events are given by a Poissonian distribution (1/τ ex ) exp (-t/τ ex ), and therefore :

This is the weighting function that has been used in Figs. 34of the main article. Of course it plays a role only for the range f ≤ f α as one has p(f ≥ f α ) ≃ 1 since f ex /f α ≪ 1.

III. AVERAGING OVER δ.

For simplicity we have presented in the main article the results of our ADWP model obtained for a single value of the dimensionless asymmetry δ. One can generalise the results by averaging over δ, at the cost of additional parameters. To investigate this question, we computed the values of X (k) 3,m for 100 values of δ linearly distributed in the [0; 0.99] interval. We then averaged the complex values of X (k) 3,m by a weight w(δ). For simplicity we have used either a flat distribution which is non zero only between δ min ≥ 0 and δ max < 1 ; or a "gaussian" distribution where

Here C is the proper normalisation constant taking into account that δ is defined only on the [0; 1] interval. Note that δ 1 is close to, but not exactly equal to, the average of δ ; and similarly δ 2 is not exactly its standard deviation due to the fact that δ is restricted to the [0; 1] interval.

Two interesting features are worth noting in this averaging procedure over δ :

• First, the values of X • Second, the shape chosen for w(δ) can strongly change the resulting X

3,m values. To investigate this point, we have fixed the two first moments of δ, and chosen accordingly the parameters δ 1 , δ 2 , δ min and δ max . It is found that X (k) 3,m can be strongly different for a gaussian weight and for a flat weight distributions. This clearly shows the strong importance of the interference effects, evoked in the main article, between the nonlinear susceptibilities of the dynamical heterogeneities corresponding to different values of τ . These interference effects are strong enough to yield, e.g., a change in the log-log slope of X 3,m (f α ) by a factor significantly different from 1 (i.e., larger than 2, or smaller than 1/2). We emphasize that the changes of X 3,m (f ). This is the reason why it is not surprising that fitting the measured values of X 3 requires different values of [N corr ] av , as in the main article. Indeed, it is very likely that the extreme simplicity of our model cannot fully capture these complicated interference effects. However, relaxing only this constraint that the values of [N corr ] av should be the same when fitting X