Mahmoud El

Chamie Chadi Barakat

Giovanni Neglia

Geographically Fair In-Network Caching for Mobile Data Offloading

Data offloading from the cellular network to lowcost WiFi has been the subject of several research works in the last years. In-network caching has also been studied as an efficient means to further reduce cellular network traffic. In this paper we consider a scenario where mobile users can download popular contents (e.g., maps of a city, shopping information, social media, etc.) from WiFi-enabled caches deployed in an urban area. We study the optimal distribution of contents among the caches (i.e., what contents to put in each cache) to minimize users' access cost in the whole network. We argue that this optimal distribution does not necessarily provide geographic fairness, i.e., users at different locations can experience highly variable performance. In order to mitigate this problem, we propose two different cache coordination algorithms based on gossiping. These algorithms achieve geographic fairness while preserving the minimum access cost for end users.

I. INTRODUCTION

The reduction of the load on the 3G/4G infrastructure by offloading mobile traffic to low-cost WiFi has been the subject of several research works in the last years [START_REF] Dimatteo | Cellular Traffic Offloading through WiFi Networks[END_REF], [START_REF] Deshpande | Performance Comparison of 3G and Metro-scale WiFi for Vehicular Network Access[END_REF]. The main idea is to redirect content requests, when it is possible, from the mobile infrastructure to WiFi hotspots. If requests are delay tolerant (videos, images, documents, emails, etc.), it is indeed possible to delay them until a WiFi opportunity arises. In [START_REF] Mehmeti | Is it worth to be patient? analysis and optimization of delayed mobile data offloading[END_REF], the authors study the performance of delaying mobile data offloading with user-defined deadlines in the framework of queuing models. In [START_REF] Balasubramanian | Augmenting Mobile 3G Using WiFi[END_REF] the authors propose schemes for the prediction of future WiFi opportunities and advocate the gain in terms of 3G traffic reduction, while in [START_REF] Whitbeck | Fast Track Article: Push-and-track: Saving Infrastructure Bandwidth Through Opportunistic Forwarding[END_REF] requests are redirected to neighboring mobile devices where the operators push the content.

In addition to traffic redirection via low-cost WiFi, innetwork caching has also been studied to further reduce cellular network traffic [START_REF] Neves Dos Santos | CEDO: Content-centric Dissemination Algorithm for Delay-tolerant Networks[END_REF], [START_REF] Bhatia | UNAP: User-Centric Network-Aware Push for Mobile Content Delivery[END_REF]. Indeed, because of the skewed content popularity, even limited-size caches at the network edges can reduce the load on the core. The authors in [START_REF] Erman | To Cache or Not to Cache: The 3G Case[END_REF] illustrate this benefit by the help of a measurement-based study of 3G HTTP traffic. The idea of in-network caching is the key principle in Information-Centric Networks (ICN), like the Content-Centric Network (CCN) paradigm [START_REF] Jacobson | Networking Named Content[END_REF].

Our present work shares the same objective of mobile traffic offloading with this set of related work taking advantage of innetwork caching. In particular we consider a scenario where WiFi-enabled caches are deployed at different locations in a

The authors are with INRIA Sophia Antipolis, 2004 Route des Lucioles, B.P. 93, 06902 Sophia Antipolis Cedex, France. Emails: shamieh.m@gmail.com, {chadi.barakat, giovanni.neglia}@inria.fr city to provide mobile users (e.g., tourists) access to a set of relatively static popular contents (e.g., maps of the city, shopping information, social media, etc.). Caches are assumed to be only loosely connected to each other and to the Internet, e.g., through some low-rate or intermittent communication infrastructure. For example they could take advantage of a Delay Tolerant Network (DTN) backbone leveraging the public transportation system [START_REF] Gaito | On the impact of a road-side infrastructure for a dtn deployed on a public transportation system[END_REF]. Fig. 1 shows our user-caches scenario.

Users can retrieve the contents at any time through their 3G data connection or wait to enter in the transmission range of a cache storing the content. The user is assumed to incur a higher cost when s/he accesses the content through the cellular network than when s/he exploits the WiFi connection, e.g., because of direct data charges, data roaming, reduction of the monthly residual amount of data available, larger battery depletion, or lower download rates. At the same time, the user is assumed to be impatient, i.e., s/he may not be willing to delay content fruition more than a certain time and can then resort to use her/his cellular data connection.

Two questions arise in such scenario. First, how should the set of contents be distributed across the caches in order to minimize users' cost to access the contents? Second, can the system benefit from some form of cooperation among the caches despite their loose connectivity? In this paper we address both these issues. We start formulating an optimization problem for contents' placement in caches (i.e., which files to put in which cache) in the presence of impatient users, and then propose an efficient algorithm to find the placements that minimize user's content retrieval cost and maximize the amount of offloaded traffic, thus relieving the cellular network. We highlight the fact that the cost-optimal distribution of con-tents does not guarantee in general geographical fairness, i.e., users experience high variability in performance at different locations in the network. We propose then a new principle of "content shuffling by gossiping" that takes advantage of the limited communication infrastructure connecting the caches. We believe the idea of local content shuffling can have applications well beyond the scenario considered in this paper.

The rest of the paper is organized as follows. In Section II we describe how we modeled the scenario. Section III studies the optimal in-network caching distribution to minimize content retrieval costs for users. We then introduce the concept of geographic fairness in Section IV and propose two caches' coordination algorithms to improve it. Section V presents some simulation results, and Section VI concludes the paper.

II. MODEL A. Scenario

Consider a network of WiFi-enabled caches that are deployed in a region where a set of popular contents (maps of the city, shops information, etc.) is made available to mobile users (e.g., tourists). The users can use their WiFi supported devices (e.g., tablets, smartphones) while moving in that area to connect to the closest cache. When the user requests a content, s/he has two ways to get it: either to download the content using her/his 3G data connection with a cost C 3G , or to wait to meet a cache storing the content and download it through WiFi with a cost C W iF i . We assume that the cost per content using the WiFi is smaller than that of 3G, i.e.,

C (i) W iF i < C (i) 3G
for every content i.

A single cache cannot store the full content catalog, hence, while moving around the city, the user connects to different caches that may or may not have the requested content. Note that if a cache does not contain the content, it will not retrieve the content for the user, but s/he can either search for the content from another cache or download it using her/his 3G connection. This is coherent with our assumption that caches are poorly connected among themselves and to the Internet. We consider impatient users, i.e., each user has a maximum time (we call it the patience time) s/he is willing to wait to download a content from a cache. If a cache with the content has not been found by such time, the user will use the 3G to get it and will then incur a higher cost.

The choice of which content to store at each cache is fundamental to satisfy a large number of users' queries and reduce their costs. Intuitively, if the patience time is the same for all the contents, highly requested contents should be stored in more caches than the less requested ones. What happens if patience times are heterogeneous is not clear (even qualitatively), and in any case a model is required to determine the exact number of replicas.

B. Contents

We suppose that the catalog size for users' requests is K, i.e., the total number of different contents or files requested by the users is K (along this paper, we use the two terms contents and files interchangeably). A request for one of the catalog contents is a random variable R ∈ {1, . . . , K} with the following discrete probability distribution:

q i = Prob(R = i) i = 1, . . . , K. (1)
Without loss of generality, we suppose that the contents are ordered in decreasing order of popularity, i.e., the first content has the highest popularity and the K-th content has the lowest one, i.e., q 1 ≥ q 2 ≥ • • • ≥ q K . Contents are assumed to have the same storage size.

C. Caches

The service provider has a fixed number N of identical caches. Each cache can store up to B contents and provide them to nearby users upon request. The caches are distributed in different geographic places in the city.

Caches have the possibility to establish a connection with each other from time to time using a low cost network (a dedicated wired network, an ad-hoc network using buses or any other mobile devices as data ferries [START_REF] Zhao | A message ferrying approach for data delivery in sparse mobile ad hoc networks[END_REF], etc.). Once two caches are connected, we assume that they can exchange contents. Therefore, the caches form a network that can be represented by an undirected graph G = (V, E), where the caches are its vertices (V = 1, . . . , N ) and each edge (i, j) ∈ E identifies that there is a recurrent communication opportunity between cache i and j (a formal communication model is introduced in Section IV).

In the first part of this work, we will study the optimal content distribution among caches, in the sense that this distribution minimizes the users' cost to retrieve the content. To this regard, the cache network topology does not play any role, because the optimization problem can be solved in a centralized way and the caches do not need to exchange any content. However, when we study the problem of geographic fairness later in the paper, we will propose algorithms that rely on contents' shuffling among the caches. Here, the inter-caches communication infrastructure is important in determining the performance of the algorithms.

D. Users

We consider a homogeneous set of users, so in this section we will refer to one of them. As mentioned earlier, the user generates requests for contents to download according to the probability distribution introduced in Section II-B. The user keeps moving inside the city and can get in contact with a given WiFi cache multiple times. We say that a user is associated to a cache if her/his WiFi device is actually associated to the cache's Access Point and can download contents. During an association, the user will complete the download of all the contents s/he is interested in and which are stored at the cache. This holds if we assume that the user stops moving until the download is complete, or that the download time is smaller than the user's association time to the cache.

Let us consider a request generated by the user at time t 0 . We define the residual time between the user and cache j as the first time after t 0 when the user associates to cache j and we denote it by X j . Through the whole paper, we consider valid the following assumption: Assumption 1. The residual times X 1 , . . . , X N between the user and the caches in the network are independent and identically distributed (i.i.d.) random variables.

To highlight the generality of this mobility model for the user, we give the following mobility model example that satisfies our assumption.

Example: Assume that for each cache j the point process of the instants when the user associates to cache j is a renewal process, where the random time Z j between two consecutive associations has Cumulative Distribution Function (CDF) F Zj (t). If the N renewal processes are i.i.d. and stationary, then Assumption 1 is satisfied. The residual time X j is the residual time for t = 0 of the corresponding renewal process [12, p. 116] and its CDF can then be calculated as follows:

F Xj (x) = x 0 FZj (τ )dτ /E[Z j ], where FZj (τ ) = 1 -F Zj (τ )
is the complementary CDF.

While our results below hold under Assumption 1 for generic distributions for the variables X 1 , . . . , X N , we often consider the exponential distribution as an example in order to derive closed form expressions and provide further insights about the problem.

The user is impatient, i.e., he is willing to wait a time at most T i before being able to download content i. If the user finds content i at one cache by time T i , s/he downloads it from the cache using the Wifi connection at a cost C (i) W iF i . If the timer expires and the user has not found the content, s/he downloads it using 3G at a cost C (i) 3G . We refer to T i as the user's patience time.

Let R be the (random) content requested at time t 0 , then the cost C R incurred to download it is:

C R = C (R) 3G 1 file downloaded by 3G + C (R)
W iF i 1 file downloaded by WiFi , where 1 is the indicator function (equal to 1 if the condition is satisfied and 0 otherwise). The randomness in C R comes from two factors: 1) the identity of the requested content and 2) the instants when the user associates to the different caches after t = t 0 .

III. TUNING CONTENT DISTRIBUTION FOR OPTIMAL

NETWORK PERFORMANCE Let C be the expected cost for the request C R , i.e., C = E[C R ]. C depends obviously on the distribution of contents among the caches. On one hand, placing more replicas of popular contents in the network of caches reduces the costs because it allows more frequent downloads from WiFi, which is less costly than 3G. On the other hand, placing popular contents everywhere reduces the remaining available space for other, less popular, contents to be stored in the network without bringing necessarily the same level of cost reduction. Therefore, due to caching space constraints of the caching network, the number of copies should be selected wisely to reduce the costs. Let us start with the following proposition that gives the expression of C. Proposition 1. The expected cost per content request for any user is

C = α + K i=1 γ i q i (1 -F Y (i) (T i )) , (2)
where α and γ i are positive constants, Y (i) is the minimum of the residual times between the user and any of the caches having content i, and

F Y (i) (T i ) = Prob(Y (i) ≤ T i ) is its cumulative distribution function.
Proof: Suppose that there are n i copies of a content i in the network (i.e., n i is the number of caches storing content i). Then the minimum of the residual times is

Y (i) = min{X 1 , X 2 , . . . , X ni },
where X 1 , . . . , X ni are i.i.d. random variables having the same distribution (due to Assumption 1). The probability that the user connects to a cache having content i within time delay T i (and downloads it using WiFi) is:

Prob[content downloaded by WiFi | it is content i] = Prob[Y (i) ≤ T i ] = F Y (i) (T i ), (3)
and the probability it is downloaded using 3G is:

Prob[content downloaded by 3G | it is content i] = Prob[Y (i) > T i ] = 1 -F Y (i) (T i ). (4)
Let A 1 be the event that the content is downloaded by WiFi and A 2 the event that the content is downloaded by 3G. Then,

E[C R ] = K i=1 C (i) W iF i Prob[A 1 | R = i] × Prob[R = i] + K i=1 C (i) 3G Prob[A 2 | R = i] × Prob[R = i] = K i=1 C (i) W iF i (F Y (i) (T i ))q i + K i=1 C (i) 3G (1 -F Y (i) (T i ))q i = K i=1 q i C (i) W iF i + K i=1 q i (C (i) 3G -C (i) W iF i )(1 -F Y (i) (T i )), then the proposition follows taking α = K i=1 q i C (i) W iF i and γ i = C (i) 3G -C (i) W iF i .
For the special case when the residual time between the user and a cache is exponentially distributed with rate λ, then Y (i) = min{X 1 , X 2 , . . . , X ni } is also exponentially distributed with rate n i λ. Therefore the cost C can be written in closed form as follows:

C = α + K i=1 γ i q i exp(-λn i T i ). (5)
Proposition 1 shows that the average cost depends on the number of replicas each content has in the network. Notice that the average cost lies between

C W iF i = K i=1 q i C (i) W iF i and C 3G = K i=1 q i C (i) 3G , i.e., C W iF i < C ≤ C 3G , where C = C 3G if
the WiFi is not used or the caches do not store any content (i.e., n i = 0 for all i), and C can get arbitrary close to C W iF i if there are enough copies of every content (i.e., lim n1,...,n K →∞ C = C W iF i ).

A. Optimizing Content Distribution

Due to the storage size constraints at caches, the service provider can optimize the performance of the network by selecting the number of replicas of each content to be placed in the caches so that the cost is minimized. The optimal number of replicas for each content can be derived solving the following optimization problem:

minimize n1,n2,...,n K f (n 1 , n 2 , . . . , n K ) = C given in (2) subject to K i=1 n i ≤ B.N, n i ∈ {0, 1, . . . , N }, for i = 1, . . . , K. (6)
The first constraint states that the total number of replicas in the network should not exceed the available space in the caches. The second set of constraints limits the number of replicas of a content to the number of available caches and restricts this number to a non-negative integer.

The result of this optimization is the number of replicas of each content to be placed in the caches. Note that any placement of that number of contents under the condition that the same content is not placed more than once in the same cache (otherwise it just consumes space without adding any value) would provide the optimal cost. This placement is clearly not unique in general. We will get back to this issue later in the paper when talking about geographically fair placements in Section IV. Before providing an efficient algorithm to solve (6), we prove the following important property of its objective function: Proposition 2. The objective function f (n 1 , . . . , n K ) is a separable convex function in the variables of its argument.

Proof: The function f is clearly separable in its variables because

f (n 1 , . . . , n K ) = α + K i=1 γ i q i (1 -F Y (i) (T i )), and F Y (i) (T i ) is a function of n i only. Then it is sufficient to prove that each of the functions 1 -F Y (i) (T i ) is convex in n i . Denote Y (n i ) = min{X 1 , . . . , X ni } where X 1 , . . . , X ni are all i.i.d. (due to Assumption 1). Let us define FY (ni) (T ) = 1 -F Y (ni) (T ) = Prob(Y (n i ) > T ).
Then we have for every

n i > 0, FY (ni+1) (T ) -FY (ni+2) (T ) = Pr(min{Y (n i ), X} > T ) -Pr(min{Y (n i + 1), X} > T ) = Prob(Y (n i ) > T ) × Prob(X > T ) -Prob(Y (n i + 1) > T ) × Prob(X > T ) ≤ Prob(Y (n i ) > T ) -Prob(Y (n i + 1) > T ) = FY (ni) (T ) -FY (ni+1) (T ),
thus f is convex in the variables of its arguments.

B. A Greedy Optimal Algorithm

The optimization problem [START_REF] Neves Dos Santos | CEDO: Content-centric Dissemination Algorithm for Delay-tolerant Networks[END_REF] gives the number of copies of every content in the network guaranteeing minimal cost. Let n * 1 , n * 2 , . . . , n * K be this optimal allocation and let C * = f (n * 1 , n * 2 , . . . , n * K ) be the corresponding minimum cost value. Let us define the following function for a content i:

D i (n i ) = γ i q i F Y (ni+1) (T i )) -F Y (ni) (T i )
where Y (n i ) = min{X 1 , . . . , X ni } and, by convention, F min{X1,...,Xn i } (T i ) = 0 if n i = 0. In the particular case of an exponential mobility model with rate λ, then D i (n i ) has the following closed form:

D i (n i ) = γ i q i (exp(-n i λT i ) -exp(-(n i + 1)λT i )) .
The optimal number of replicas of each content can be given by the following greedy algorithm:

Algorithm 1 Calculate n * 1 , n * 2 , . . . , n * K of optimization (6) 1: n 1 = n 2 = • • • = n K = 0 2: S = {1, 2, . . . , K} 3: V i = D i (n i ) ∀i ∈ S 4: while i n i < B.N do 5: j = Argmax i∈S V i 6: n j ← n j + 1 7: V j ← D j (n j ) 8: if n j = N then 9:
Remove j from S 10: end if 11: end while 12: return n 1 , n 2 , . . . , n K Proposition 3. Algorithm 1 returns the optimal solution n * 1 , n * 2 , . . . , n * K of optimization (6). Proof: At every iteration, the algorithm is adding one content having the maximum difference V i (steps 5 and 6 in the algorithm) and removing the content if n i > N (step 9). Thus Algorithm 1 is eventually equivalent to calculating the elements D i (m), for i = 1, . . . , K and m = 0, . . . N -1, and then considering the contents corresponding to the largest BN values among those elements. The latter can be shown to provide the optimal solution by adapting the proof for a similar algorithm in the context of the discrete resource allocation problem with a separable convex function given in [START_REF] Ibaraki | Resource Allocation Problems: Algorithmic Approaches[END_REF]Chapter 4]. Also in our case, what is fundamental for the proof is that the objective function is separable convex as we showed in Proposition 2. Additionally, this algorithm can be implemented efficiently using the max heap data structure [START_REF] Cormen | Introduction to Algorithms[END_REF] with a computational complexity O(K +BN log K), and memory complexity O(K).

We give an example of optimal content distribution in Fig. 2. This numerical result corresponds to a mobility model with exponentially distributed residual times with rate λ, a Zipf law for content popularity, and two different cases of the patience time. For the case when T i is also assumed to follow Zipf, 1the figure shows that the replicas are not proportional to their popularity. Less popular contents (contents 8 and 9 in the figure) have the highest number of replicas in the network. The reason for the different number of replicas is that it is not just the popularity that plays a role in selecting the items, but there are other factors. On one hand, the cost functions are convex decreasing in the number of replicas, so adding one replica does not reduce the cost as much as it did a previous replica of the same content. On the other hand, the cost also depends on the patience time (T i , for i = 1, . . . , K) by the user for the contents. In the first example, the patience time by users for the more popular contents is assumed to be larger than for the other contents which explains why less replicas are sufficient. In the case when T i is constant (T i = 0.0067 in this case), the first 4 contents have the same optimal number of replicas (50 in this case) and then that number starts decreasing again with the popularity. Note that contents that are not shown in the figure (from 23 up to 10000) are not cached in the network (i.e., n * i = 0).

IV. FAIRNESS

In the previous section we derived the optimal number of copies for each content to minimize the global content retrieval cost for the user. Any placement of this number of copies (as far as the placement ensures that the same cache does not contain more than one replica of the same content) guarantees the minimal cost. This is due to the implicit assumption of homogeneous user mobility, so that residual times between the user and each cache are identically distributed. However, different placements could lead to the user experiencing different quality-of-service at different caches. For example, this would happen for constant patience times, if one cache would store only copies of the most popular contents and another only copies of the least popular contents. We denote this problem as geographical unfairness. The global optimization problem considered until now provides minimum global cost but not geographical fairness.

It appears natural to formulate then a different optimization problem, where we try to optimally allocate the replicas among the caches, so that the expected cost at each cache is almost the same. To argue formally, let us quantify the performance of our proposed WiFi caches network. Since the aim of this network is to reduce the costs for users, we can define the Fig. 3. Variability in caches' utilities shows that, even with optimal number of contents, the placement of these copies in the network can cause variable geographic performances (N = 50 caches randomly deployed, K = 10000 files, B = 10 buffer size, q i and T i are zipf distribution with parameter s = 1, λ = 0.04 is constant).

gain P of the network as the expected reduction of retrieval cost per request due to in-network caching. Since γ i is the reduction in the cost when content i is downloaded from the caches and F Y (i) (T i ) is the probability of finding this content, then

P = K i=1 q i γ i F Y (i) (T i ).
The utility of a cache would then be its contribution to the performance,

V j = contents i in cache j q i γ i F Y (i) (T i ) n * i ,
where we divided by n * i because when a request for a content i is satisfied, it is equally likely to be satisfied at any of the n * i caches having a replica of this content. Therefore, the utility (or value) U i of a replica of content i is given by

U i = q i γ i F Y (i) (Ti) n * i if n * i > 0, 0 if n * i = 0. (7)
For exponentially distributed residual times it is

U i = γ i q i 1-exp(-n * i λTi) n * i if n * i > 0, 0 if n * i = 0, (8)
so U i can be considered as a measure for the importance of a replica of a content.

A numerical example for the potential high variability of caches' utilities is given in Fig. 3. The figure shows that the optimal number of copies does not necessarily provide geographic fairness. In this specific example, the caches having the most popular contents happen to be placed at the edges of the geographic area, and thus users moving close to these areas would get higher satisfaction compared to users moving near the regions where caches have less popular contents. The question that arises here is the following: how can we preserve optimal performance (having always the optimal number of copies in the caches) and at the same time get geographic fairness (same geographic performance everywhere in the area)?

So our goal is to achieve lower performance variability across the area. Since there are different ways of distributing the replicas in the network, and each replica has a different performance value, for the purpose of fairness, it is desirable that this distribution of replicas minimizes the variability of values in different caches. Let us define

V max = max j∈V V j ,
which is the maximum utility of a caches. Therefore, to achieve fairness, we should distribute the contents in the caches in a way that minimizes V max while still guaranteeing the minimum cost. This problem can be formulated as the following Mixed Integer Linear Program (MILP):

minimize y subject to K i=1 U i x ij ≤ y, for j = 1, . . . , N, K i=1 x ij = B, for j = 1, . . . , N, N j=1 x ij = n * i , for i = 1, . . . , K, x ij ∈ {0, 1}, y ∈ R, (9)
where the variable x ij indicates whether a content i is placed in cache j or not. The newly introduced real variable y together with the first set of constraints guarantee that V max is minimized. The second set of constraints guarantees that each router is not exceeding its buffer capacity. And the third set of constraints is to ensure that the number of replicas for a content does not exceed its optimal value. However, this optimization turns out to be computationally hard, Proposition 4. The fairness problem (9) is NP-complete.

Proof: The proof is due to a reduction from the 3-partition problem (known to be NP-complete). The 3-partition problem is to decide whether a given set of S integers m 1 , . . . , m S can be partitioned into triples that all have the same sum. We show that given any 3-partition problem, it can be solved by the formulation (9) with a particular choices T i , q i , γ i , B, N , and K. There are different possibilities for these choices, one of them is to take B = 3, N = S/3, K = S to be the number of contents, and U i = m i (obtained by taking T i = ∞ so that n * i = 1 for all i = 1, . . . , K, γ i = i m i for all i, and q i = m i / i m i ). Thus by solving this problem by ( 9) and checking, in polynomial time, if the values of the caches are equal, we can solve the original 3-partition problem. Therefore, the optimal placement of contents in the caches is at least as hard as the 3-partition problem and this ends the proof.

Given the difficulty to determine the static content distribution which maximizes fairness, an alternative approach could be to achieve it on the long run by dynamically changing content distribution among the caches. A naive idea is to periodically reshuffle the copies among the caches, by calculating in a centralized way a random distribution of the copies and transmit them to each cache accordingly. This solution guarantees that all the caches have the same timeaverage utility on the long run (because of the strong law of large numbers), but has significant communication cost, given that N B copies should be re-transmitted at every iteration. In order to improve on top of this idea by reducing the amount of communication required, we consider shuffling algorithms based only on local exchanges among caches as in gossiping protocols [START_REF] Bénézit | Order-optimal consensus through randomized path averaging[END_REF] [START_REF] Boyd | Randomized gossip algorithms[END_REF], without the need to download B copies at each cache.

The movement of contents between caches is done by exchanging files through the communication network connecting these caches. We consider an asynchronous operation, where pairs of caches can autonomously decide to exchange some contents taking advantage of their mutual connectivity, which, as we mentioned before, may also be intermittent. We model then the network as discrete-time system where at every time slot two caches exchange files, then k exchanges occur in the network by slot k. Let G(k) be the graph representing the interaction at time k: this graph has the N caches as vertices and only one edge connecting the two caches that exchange contents at time k. We consider that the interaction among caches can be modeled probabilistically as follows:

Assumption 2. At a given communication iteration, the probability that caches i and j exchange contents is given by the following expression,

Prob((ij) ∈ G(k)) = σ ij > 0 if (i, j) ∈ E 0 else;
The probabilities σ ij take into account both the probability there is a communication opportunity between caches i and j (as we said the network could be only intermittently connected) and the decision of the caches to actually exchange some contents. Given this model, the union of all possible graphs G(k) is the connected graph G with probability one.

Once the connection between two caches is established, these caches should choose which files (contents) to exchange. We propose here two exchange rules for the files:

• Exchange Rule 1: Common files are not exchanged, only different files are randomly divided between the two caches. • Exchange Rule 2: Common files are not exchanged. The connected caches exchange one file (if it exists) that gives the maximum decrease in the local objective function V

(ij) max = max{V i , V j }. Let V j (k) be the value of the cache j after the k-th exchange such that V j (0) = V 0 j are the values of the caches due to an initial distribution of the optimal copies in the network (the initial values vector is denoted by V 0 ).

For both the exchange rules the total utility of the caches ( j V j (k)) does not vary neither with k nor with the initial placement, i.e.,

N j=1 V j (k) = K i=1 n * i U i , ∀k, ∀V 0 .
This is evident if we think that the total number of copies per content does not change or if we consider that at every iteration, local utilities change as follows:

V i (k) + V j (k) = V i (k -1) + V j (k -1) if (i, j) ∈ G(k) V i (k) = V i (k -1) else.
(10) Utilities of caches evolve randomly under both the exchange rules. For exchange rule 2 the randomness is only due to the random selection of the caches which interact. For exchange rule 1 the randomness is also due to the random reallocation of the files among the caches. For this reason it is useful to consider the following quantities

x j (k) = V 0 j if k = 0 E(V j (k)) if k > 1,
where the expected value is calculated over the possible sequence of pairs of interacting caches for exchange rule 2 and over the sequence of pairs of interacting caches and possible re-allocations of the files for exchange rule 1. We can now give the following definition:

Definition 1. We say that the caches achieve asymptotic geographic fairness under a given exchange rule if for any initial optimal distribution of contents (i.e., any

V 0 ) 1) lim k→∞ x j (k) = x for j = 1, . . . , N , 2) x = 1 N K i=1 n * i U i .
The first condition says that all caches would have on average the same value, i.e., it is indifferent to be in an area or in another because the caches have the same value. The other condition says that the distribution of contents remains optimal and its value is independent from the initial placement of contents (actually this corresponds to averaging the initial utilities of caches in the network).

Although the similarities illustrated above for the two exchange rules, there are some evident differences. For exchange rule 2 we have

|V i (k) -V j (k)| ≤ |V i (k + 1) -V j (k + 1)|
for every iteration k. Given that the sum of the utilities is constant, this inequality implies that the maximum over all the utilities (max i V i (k)) is a non-increasing sequence and then converges to a given value. This value is a local minimum that the greedy local exchange rule 2 cannot improve. On the other side, utilities in exchange rule 1 do not converge in general, because at any time slot some contents are randomly reallocated. Nevertheless, exchange rule 1 has the following desirable property: Proposition 5. Under Assumption 2 and the Exchange Rule 1, caches achieve asymptotic fairness.

Proof: The proof is provided in the Appendix.

V. PERFORMANCE STUDY A. Comparison of Optimal Placement with LRU Caches

An alternative approach to content placement is to use LRU (Least Recently Used) replacement policy in a noncoordinated way where each cache applies its own policy as a function of the traffic it sees. For our particular scenario where caches see the same traffic distribution, they would have the same probability distribution of having a content (mainly most popular contents are found everywhere). To compare with our optimal placement, we assume that when a user does not find a requested file in an LRU cache, it continues searching in nearby caches (according to the mobility model) and does not wait for the cache to download the content (even if the cache will eventually replace one of its existing file with the recently missed file using 3G synchronization, we don't include this O p t i mal D i st r i b u t i on T i = c on s t a n t L RU C h e Ap p r ox . T i = c on s t a n t Fig. 4. Cost comparison between having LRU replacement policy in caches (C LRU ) and the cost due to optimal replicas from problem (6) (C * ) while varying the number of caches in the network from 1 to 100 (α = 0, γ i = 1, K = 10000 files, B = 10 buffer size, q i and T i are zipf distribution with parameter s = 1, λ = 5 is constant). latter cost for the sake of comparison). Having said that, similar to the derivation of the expected cost in Proposition 1, we can write the expected cost per request C LRU for using LRU caches with exponential mobility model of rate λ as follows:

C LRU = α + K i=1 γ i q i E[exp(-Z i λT i )], (11)
where Z i is a random variable for the number of caches having content i. Since each of the caches are considered as independent LRU cache, then we have E[Z i ] = h i N where h i is the probability that a file i is present in one LRU cache (hit probability for file i). By applying Jensen's inequality for the convex functions φ i (Z i ) = exp(-Z i λT i ), we can derive a lower bound C LRU on C LRU : 2

C LRU ≥ α + K i=1 γ i q i exp(-h i N λT i ). (12)
Notice that by comparing equations ( 5) and ( 12), n i is now replaced by h i N , which is no other than the mean number of replicas of content i across the network of LRU caches. To find h i , we apply "Che approximation" [START_REF] Che | Hierarchical Web caching systems: modeling, design and experimental results[END_REF], which is proved to be very accurate [START_REF] Fricker | A Versatile and Accurate Approximation for LRU Cache Performance[END_REF], given as follows:

h i ≈ 1 -exp(-q i t C ) for i = 1, . . . , K, (13)
where t C is the unique root of the following equation

K i=1 (1 -exp(-q i t)) = B.
Fig. 4 shows that the average cost per request when caches run optimization problem [START_REF] Neves Dos Santos | CEDO: Content-centric Dissemination Algorithm for Delay-tolerant Networks[END_REF] is less than the one in the case caches use LRU replacement policy. Moreover, when the number of caches is small, adding one more cache in the 2 Alternatively, another approach can be used to derive the same lower bound, we observe that E[φ i (Z i )] is the probability of not finding content i in the network (i.e., probability of a miss) which can be calculated precisely by the following formula E[φ i (Z i )] = (h i exp(-λT i ) + (1 -h i )) N . Therefore, the mentioned lower bound can be also easily derived since network reduces quickly the cost for the optimal placement, which does not hold for the lower bound cost C LRU for LRU caches. At one point, it also seems from the figure that the difference between the two costs tends to be constant. It is important to recall that the cost given in the figure for the LRU does not include the cost of synchronization using the 3G (i.e., practically speaking we should add to the LRU cost curve in the figure the cost of downloading a new content after every miss). Having said that, the LRU has the advantage that it will automatically adapt to changes in the popularity distribution (q i 's) while in our case we should rerun our optimization problem and perform a new placement (that can be done incrementally starting from the existing one). We leave the study of dynamic traffic to our future research.

h i exp(-λT i ) + (1 -h i ) ≥ exp(-h i λT i ).

Number of iterations

B. Fairness Exchange Rules

As we have seen so far, fairness is not granted, and the caches should run an additional gossiping protocol to guarantee fairness while preserving optimal performance. The two exchange rules for gossiping the content provide two different (orthogonal) approaches to the problem. Fig. 5 shows the resulting utility of caches after applying 500 iterations of gossiping between caches. We consider a random geometric graph (RGG) 3 topology with n = 50 caches that communicate according to Assumption 2. We consider that all the edges in the network have the same rate σ ij = cte for all (ij) ∈ E. Thus, each edge can be selected with the same probability and the expression in equation ( 16) for the weight matrix W = E[W (k)] reduces to:

W = I + 1 2m (A -D),
where I is the identity matrix, A is the adjacency matrix of the graph (having entries a ij = 1 if (ij) ∈ E and 0 otherwise), D is a diagonal matrix having the degrees of the nodes, and m is the number of edges in the graph (i.e., m = |E|). We just remind the reader that in our fairness model, the lower V max achieved, the fairer is the content distribution. Considering Exchange Rule 1, the expected utility curve (in red) in the figure shows the expected value (or utility) of a given cache in the network, that can be calculated using the expected weight matrix W . The dotted (light blue) curve shows instead the evolution of the actual utility of a given node. If we imagine to repeat an infinite number of times the process, and average all the value curves for the same node, then we would get exactly the red one. The solid line curve (dark blue) shows the maximum value in the network (V max ). This curve shows that Exchange Rule 1 does not provide fairness at any time slot because V max always has high values. At the same time each cache is changing its contents regularly over time, and thus experiences different utilities such that, on average, all caches are equivalent. Fairness is only guaranteed by Exchange Rule 1 in this sense.

On the contrary, for Exchange Rule 2 the maximum value V max (green curve) keeps decreasing. In fact, as we have explained in the previous section, caches are performing a sort of distributed local search, by exchanging contents only when the total objective function V max can be reduced. This gives better fairness as time goes on in comparison to the initial placement (reduced from V max = 1.6 at iteration 0 till V max = 1.05 at iteration 500). Nevertheless, the final value to which V max converges is just a local optimum. Notice that the static placement that optimizes V max should give a value that lies between the red and the green curve (i.e., the optimal value in this case is V * max ∈ [0.7, 1.05]).

VI. CONCLUSION

We have presented in this paper an optimization problem for the optimal caching in mobile networks. Caching is supposed to be done at WiFi hotspots that are spread across a city. Caching is meant to offload traffic of the 3G/4G infrastructure, and hence to reduce costs for both the operator and the users. We showed that by an optimal replication and placement of contents in caches accounting for their popularity, the average cost per requested content is minimized. Replicas are then dynamically shuffled between caches to provide a balanced performance for users in different geographic places while preserving optimality for global network performance. We leveraged the theory of gossiping to realize content shuffling between caches and developed algorithms to ensure that the user cost is indifferent with respect to her/his place in the network. Up to our knowledge, this work is a first attempt to explore the problem of unfairness that might result from content placement, and that proposes shuffling of the content of caches to achieve fairness without compromising global performance. For future work, we will keep investigating this direction by focusing in particular on the effect of traffic dynamicity and the way it interacts with our developed algorithms.

APPENDIX PROOF OF PROPOSITION 5

To show the fairness property, let us study the dynamics and the convergence properties of the vector x(k) where x i (k) is its ith element. Notice that with every iteration (exchange), we have

x i (k + 1) = x j (k + 1) = xi(k)+xj (k) 2
if (i, j) ∈ G(k) x i (k + 1) = x i (k) else, (14) then we can write in matrix form

x(k + 1) = W (k)x(k), (15)
where W (k) is a matrix having off-diagonal elements equal to zero except for w ij = w ji = 0.5 such that (i, j) ∈ G(k), and all diagonal elements equal to 1 except for w ii = w jj = 0.5. Then the matrix W (k) satisfies the following properties for every k: W (k)1 = 1, W (k) = W (k) T , and W (k) 2 = W (k) (where 1 is the vector of all ones). These properties are important for studying the stability of gossip algorithms (see [START_REF] Bénézit | Order-optimal consensus through randomized path averaging[END_REF], [START_REF] Boyd | Randomized gossip algorithms[END_REF]). Notice that due to Assumption 2, W (0), W (1), . . . are all i.i.d. matrices. Let W = E[W (k)], then we have

W = (i,j)∈E Prob[(i, j) ∈ G(k)]W (k) = (i,j)∈E σ ij W (ij) , (16)
where W (ij) is a doubly stochastic matrix having off-diagonal elements equal to zero except for w ij = w ji = 0.5. Notice that W is a positive definite matrix (having dominant diagonal entries) and it is irreducible because we assume that the union graph G is connected. Thus we can say that the second largest eigenvalue (in magnitude) of W , λ 2 (W ), satisfies

0 ≤ λ 2 (W ) < 1, λ 2 (W ) = λ 1 (W -J), (17)
where λ 1 (.) is the largest eigenvalue (in magnitude) of a matrix, and J = 1 N 11 T . To show the convergence of equation [START_REF] Bénézit | Order-optimal consensus through randomized path averaging[END_REF]

E[α k ] ≤ Aλ k 2 , (18)
where A = ||y(0)|| 2 2 is just a constant. Using the Markovian inequality, we get that for any > 0,

Prob[α k > ] ≤ E[α k ] ≤ Aλ k 2 ,
and thus

∞ k=0 Prob[α k > ] ≤ A ∞ k=0 λ k 2 = A (1 -λ 2 ) < ∞,
then since the rate of convergence of E[α k ] is fast enough ( ∞ k=0 Prob[α k > ] < ∞), then α k converges almost surely to 0 [19, p. 138]. Therefore, x j (k) converges almost surely to x for j = 1, . . . , N and this ends the proof because all properties of geographic fairness with optimality in Definition 1 are satisfied.

Fig. 1 .

 1 Fig. 1. The mobility-network model for a user.

Fig. 2 .

 2 Fig.2. The optimal distribution of contents as function of their popularity given by Algorithm 1 for two different cases of the patience, T i being zipf and T i being constant (α = 0, γ i = 1, K = 10000 files, N = 50 caches, B = 10 buffer size, q i is zipf distribution with parameter s = 1, and λ = 5 is constant).

  i mal D i st r i b u t i on T i = z i p f L RU C h e Ap p r ox . T i = z i p f

Fig. 5 .

 5 Fig. 5. Utility of caches by considering different exchange rules.

  under Assumption 2 and the Exchange Rule, we define y(k) = x(k) -x1 such that x = 1 i U i . We will thus show the two properties for geographic fairness of Definition 1 by showing that y(k) converges to 0 almost surely. Letα k = ||y(k)|| 2 2 = y(k) T y(k), then we have E[α k+1 |y(k)] = E y(k) T W (k) T W (k)y(k)|y(k) = E y(k) T W (k)y(k)|y(k) = y(k) T W y(k) = y(k) T (W -J)y(k) ≤ λ 2 (W )y(k) T y(k). Since E[α k+1 ] = E [E[α k+1 |y(k)]],then by a simple recursive equation we get that

A Zipf distribution for T i may be justified if the user knows contents popularity. For a more popular content, the user would rather be more patient with the hope of finding it later in the network (because it is popular), but s/he would not be as patient for a less popular one as s/he knows that s/he is less likely to find it.

RGGs are random graphs where n nodes are placed uniformly at random in a 1 by 1 unit square area, and any two nodes within a connectivity threshold forms an edge in the graph.

ACKNOWLEDGMENT

The authors would like to thank Nicolas Nisse for his feedback regarding the NP-completeness of the MILP problem (9).


Mathieu Langlais

Jean-Paul Hugonin

Mondher Besbes

Philippe Ben-Abdallah

Cooperative electromagnetic interactions between nanoparticles for solar energy harvesting

The cooperative electromagnetic interactions between discrete resonators have been widely used to modify the optical properties of metamaterials. Here we propose a general approach for engineering these interactions both in the dipolar approximation and for any higher-order description. Finally we apply this strategy to design broadband absorbers in the visible range from simple n-ary arrays of metallic nanoparticles.

Introduction

Engineering light-matter interactions is a longstanding problem in physics and is of prime importance for numerous technological applications such as the photo and thermophovoltaic energy conversion [START_REF] Atwater | Plasmonics for improved photovoltaic devices[END_REF]2], the optical manipulation of nanoobjects [START_REF] Ashkin | Observation of a single-beam gradient force optical trap for dielectric particles[END_REF] or the quantum information treatment [START_REF] Poddubny | Spontaneous radiation of a finite-size dipole in hyperbolic media[END_REF]. Light interaction with resonant structures embedded inside a material is a natural way to modify its optical properties. To date, a large number of resonant structures have been developped following such a strategy. Among these, metamaterials based on metallo-dielectric structures have been proposed to operate at frequencies ranging from the microwave domain [START_REF] Landy | Perfect metamaterial absorber[END_REF] to the visible [START_REF] Aydin | Broadband polarization-independent resonant light absorption using ultrathin plasmonic super absorbers[END_REF].

The design of artificially constructed magnetodielectric resonators which strongly interact cooperatively is a very recent and promising way to generate metamaterials that highlight innovative physical [START_REF] Jenkins | Metamaterial transparency induced by cooperative electromagnetic interactions[END_REF][START_REF] Fedotov | Spectral collapse in ensembles of metamolecules[END_REF] and transport [START_REF] Ben-Abdallah | Heat superdiffusion in plasmonic nanostructure networks[END_REF]10] properties. However, so far, only heuristic approaches have been followed to identify the convenient meta-structures which display target functionalities. In this paper, we present a general theory to describe the multiple scattering interactions mechanisms in discret networks of resonators embedded in a host material and we propose a general method to identify the appropriate inner structure of networks that highlight a targeted optical property either by considering the interacting objects as simple (electric and magnetic) dipoles or as multipoles of arbitrary order. To illustrate the strong potential of cooperative interactions to tailor the optical properties of materials we design a broadband light absorber made with simple binary lattices of metallic nanoparticles immersed in a transparent host material.

Scattering by nanoparticle networks in the dipolar approximation

To start, let us consider a set of objects dispersed inside a host material as depicted in Fig. 1. Suppose this system is higlighted by an external harmonic electromagnetic field of wavelength much larger than the typical size of objects. In this condition, we can associate to each object an electric (E) and magnetic (H) dipole moment p m;A (A = E, H) (the higher orders contributions are discussed in the next paragraph). The local electromagnetic field A ext (r m ) at the dipoles location r m results from the superposition of external incident field, the field generated by the others dipoles and the auto-induced field which comes from the interactions with the interfaces. Therefore it takes the form [START_REF] Purcell | Scattering and absorption of light by nonspherical dielectric grains[END_REF][START_REF] Draine | Discrete-dipole approximation for periodic targets: theory and tests[END_REF][START_REF] Evlyukhin | Optical response features of Si-nanoparticle arrays[END_REF]]

A ext m = A inc m -iω ∑ B=E,H Γ AB (∆G AB mm p m;B + ∑ n =m G AB mn p n;B ), (1)
where

Γ EE Γ EH Γ HE Γ HH ≡ iω µ 0 ω/c ω/c -iωε 0
and G AB mn is the dyadic Green tensor in the host material which takes into account the presence of interfaces [START_REF] Tomas | Green function for multilayers: light scattering in planar cavities[END_REF] and gives the field A at the position r m given a B-dipole located in r n . ∆G AB defined as ∆G AB ≡ G AB -G AB 0 gives the contribution of interfaces only.

Here

G AB 0 (r m , r n ) = exp(ikr mn ) 4πr mn × 1 + ikr mn -1 k 2 r 2 mn 1 + 3-3ikr mn -k 2 r 2 mn k 2 r 2 mn r mn ⊗ r mn i f A = B ikr mn -1 kr mn L i f A = B
is the free space Green tensor in the host material defined with the unit vector r mn ≡ r mn /r mn . r mn denotes here the vector linking the center of dipoles m and n, while r i j =| r i j |, k is the wavector,1 the unit dyadic tensor and L =   0 -r mn,z rmn,y rmn,z 0 -r mn,x -r mn,y rmn,x 0   . Beside the dipoles location the auto-induced part of field does not exist anymore and it takes the simplified form

A ext (r) = A inc (r) -iω ∑ B=E,H Γ AB∑ j G AB (r -r j )p j;B . (2)
It immediately follows that, the dipolar moments associated to each object reads

p m;A = χ A ← → α m;A A ext m ( 3
)
where χ A represents either the vacuum permittivity ε 0 or the vacuum permeability µ 0 and ← → α i,A is the free polarizability tensor of m th object under the action of field A. By inserting the external contribution (1) of local field into relation [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF] we get the following system which relates all dipole moments

p m;A = χ A ← → α m;A [A inc m -iω ∑ n ∑ B=E,H G AB reg (r m , r n )p n;B ]. (4)
Here, where we have introduced the regularized Green tensor

G AB reg (r, r ) = Γ AB G AB (r, r ) i f r = r Γ AB ∆G AB (r, r ) i f r = r . ( 5
)
In the particular case of n-ary periodic lattices made with n arbitrary dipoles of free electric and magnetic polarizability ← → α m;A=E,H distributed in a unit cell we have, according to the periodicity, the supplementary relations for the incident fields A inc j β = A β exp(ik // .r j β ) and for the dipolar moments p j β ;A = p β ;A exp(ik // .r j β ). Here r j β is the position vector of the β th dipole inside the unit cell j of lattice and k // is the parallel component of wavector. Accordingly, Eq. ( 4) can be solve with respect to the incident field to give

p E p H = A -1 H E H . (6)
Here we have set

p A=E,H = p 1,A , • • • , p n,A t and A = A 1 , • • • , A n t
and we have define the block matrixes

H = diag(ε 0 ← → α 1;E , . . . , ε 0 ← → α n;E , µ 0 ← → α 1;H , . . . , µ 0 ← → α n;H ) (7)
and

A =                (1 + U EE 11 ) U EE 12 • • • U EE 1n U EH 11 • • • U EH 1n U EE 21 . . . . . . . . . . . . . . . . . . (1 + U EE nn ) U EE n-1,n U EH n1 • • • U EH nn U EE n1 • • • U EE n,n-1 (1 + V HH 11 ) V HH 12 • • • V HH 1n V HE 11 • • • V HE 1n V HH 21 . . . . . . . . . . . . . . . . . . . . . . . . V HH n-1,n V HE n1 • • • V HE nn V HH n1 • • • V HH n,n-1 (1 + V HH nn )                (8)
with

U EA lk = iε 0 ω ← → α l;E ∑ j G EA reg (r 0l , r jk )e ik // .(r jk -r 0l ) , (9)
V HA lk = iµ 0 ω ← → α l;H ∑ j G HA reg (r 0l , r jk )e ik // .(r jk -r 0l ) . (10)
These summations can calculated directly or using the Ewald's method [START_REF] Ewald | Die berechnung optischer und elektrostatischer gitterpotentiale[END_REF][START_REF] Capolino | Efficient computation of the 2-D Green's function for 1-D periodic structures using the Ewald method[END_REF] as in solid-state physics. Relation (6) defines the dressed polarizability tensor [START_REF] Castanie | Absorption by an optical dipole antenna in a structured environment[END_REF] Λ

≡ A -1 H = Λ EE Λ EH Λ HE Λ HH (11)
of resonators system within the unit cell of lattice. It takes into account both the intrinsic properties of isolated objects and their interactions with the environnement [START_REF] Ben-Abdallah | Many-body radiative heat transfer theory[END_REF] (particles and interfaces). The dipspersion relations of resonant modes inside the system of coupled resonators is then given by the eigenvalues of the dress polarizability tensor.

Beside the spectrum of nanoresonators network we can calculate the amount of energy which is dissipated by the electromagnetic field inside each resonator. According to the Pynting's theorem [START_REF] Jackson | Classical Electrodynamics[END_REF] the power dissipated at a frequency ω inside the m th resonator is given by the rate of doing work by the electric and mgnetic fields inside the resonator volume

V m P m (ω) = 1 2 ∑ A=E,H V m Re[j * m;A (r, ω).A(r, ω)]dr. ( 12
)
Here A denotes either the local electric or magnetic field E and H while j E and j H are the corresponding local current density. In the dipolar approximation

j m;A = -iωp m;A δ (r -r i ), expression (12)
can be recasted into the discrete form

P m (ω) = - ω 2 ∑ A=E,H {Im[p * m;A (ω).A ext m (ω)] - ω 3 µ 0 2 p * m;A Im[G AA 0 (r m , r m )]p m;A }. (13)
By inverting (1) after having replaced the dipole moments by their expression with respect to A ext m , we can express A ext m in term of A inc and explicitely calculate the power dissipated in each object under an external lighting.

For spherical particles of radius R the polarizability is straightforwardly derived from the Mie scattering theory [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF][START_REF] Evlyukhin | Collective resonances in metal nanoparticle arrays with dipole-quadrupole interactions[END_REF]. If those particles, of refractive index n m , are immersed inside a medium of index n h , we have

← → α A = α A 1 with α E -1 = k 3 0 n h 6π (C E -i), (14)
α H -1 = k 3 0 n 3 h 6π (C H -i). (15)
Here k 0 is the wavevector inside vacuum and

C E = ρ 2 m -ρ 2 h ρ 2 m ρ 2 h (Cosρ h + ρ h Sinρ h )(Sinρ m -ρ m Cosρ m ) + ρ m Cosρ h Cosρ m + ρ h Sinρ h Sinρ m ρ 2 h -ρ 2 m ρ 2 m ρ 2 h (Sinρ h -ρ h Cosρ h )(Sinρ m -ρ m Cosρ m ) -ρ m Sinρ h Cosρ m + ρ h Cosρ h Sinρ m , (16)
C H = -ρ 2 h Cosρ h (Sinρ m -ρ m Cosρ m ) + ρ 2 m Sinρ m (Cosρ h + ρ h Sinρ h ) ρ 2 h Sinρ h (Sinρ m -ρ m Cosρ m ) -ρ 2 m Sinρ m (Sinρ h -ρ h Cosρ h ) (17)
with ρ h = k 0 n h R and ρ m = k 0 n m R, n m being the refractive index of resonator. According to Eqs. ( 13), ( 16) and ( 17) it follows that the power dissipated in each particle can be expressed both in term of absorption cross-sections and of incident external field

P m (ω) = - ω 2 {ε 0 n h ω 3 6πc 3 Im[E ext * m (C E ← → α * E,m ← → α E,m )E ext m ] +µ 0 n 3 h ω 3 6πc 3 Im[H ext * m (C H ← → α * H,m ← → α H,m )H ext m ]} (18)

Generalization of scattering problem beyond the dipolar approximation

So far, we have only considered interactions between electric (resp. magnetic) dipoles. In this paragraph we describe how to take into account the multipolar interactions. The electromagntic field inside a medium of refractive index n h can be expressed in term of ingoing (-) and outgoing (+) vector spherical wave functions (which form a complete basis)

ψ ± pq = E ± pq H ± pq ( 19
)
where we have adopted the usual convention for the multipolar index (m, n) which are replaced by a single index p = n(n + 1) + m and where q set the polarization state (i.e. q = 1 for T E waves and q = 2 for TM waves).

The outgoing wave functions ψ + pq are solutions of Maxwell's equation (using the e -iωt convention)

∇ × E + pq = iω µH + pq + H S pq ∇ × H + pq = -iωεE + pq + E S pq . (20)
with the source term ψ S pq = E S pq H S pq (i.e. the multipole pq) defined by

E S pq = 0 H S pq = in 1/2 h r.D m n . ( 21
)
for the magnetic contributions and by

E S pq = -in -1/2 h r.D m n H S pq = 0 . ( 22
)
for the electric ones where D m n expresses in term of partial derivative of the Dirac distribution δ as

D m n = i (2k 0 n h ) n n! 8π(-1) m (2n + 1)(n -m)! n(n + 1)(n + m)! × z( ∂ ∂ x + i ∂ ∂ y ) -(x + iy) ∂ ∂ z (n+m) (- ∂ ∂ x + i ∂ ∂ y ) (n) δ . (23)
Let us consider an isolated particle of arbitrary shape immersed inside a host medium of index n h and highlighetd by an external electromagnetic field. By definition, this field can be decomposed on the complete basis of ψ ± pq as

A inc (r) = ∑ pq A pq inc ψ + pq (r) + ψ - pq (r) 2 . ( 24
)
Fig. 2. On the first column, absorption of simple and binary hexagonal lattices made with Ag and Au nanoparticles 30 nm radius immersed at h = 100nm from the surface in a transparent host medium of index n h = 1.5 with respect to the density in particles. On the second column, this absoption is compared with the absorption of single particles without multiple scattering interaction and, on the last column, with the results given by the effective medium theory with the same filling factor.

As for the diffracted field, which is an outgoing field, it reads

A di f f (r) = ∑ pq A pq di f f ψ + pq (r). (25)
To calculate the components A pq di f f of the diffracted field we first calculate, by reciprocity and using the Lorentz relations [START_REF] Jackson | Classical Electrodynamics[END_REF] , the action of the sources ψ S p q = E S p q H S p q

(with p = (n, -m)) on the incident field. Then using the fact that ψ ± pq functions form a complete basis they satisfy the following orthogonality relations

< ψ ± pq , ψ ± p q >= -4iδ pq,p q < ψ ± pq ψ ∓ p q >= 0 . ( 26
)
where the brackets < ., . > represents the scalar product in the Lorentz sens defined by

< ψ 1 pq , ψ 2 pq >= (E 1 × H 2 -E 2 × H 1 ).ndS. (27)
Here, integration in taken over an oriented surface (with the vector n) surrounding the particle. It follows by applying the Lorentz relation with the field ψ + p q generated by a source ψ S p q and the incident field A inc that < ψ + p q , A inc >= ψ S p q (r).A inc (r)dr ≡ I ψ S p q

[A inc ]. (28)
Note that I ψ [A] is the action of the distribution on the test function ψ [START_REF] Schwartz | Théorie des distributions[END_REF]. Then using the orthogonality relations [START_REF] Aubry | Plasmonic light-harvesting devices over the whole visible spectrum[END_REF] and according to (24)

A pq inc = i 2 I ψ S p q [A inc ]. (29)
Then, using the matrix T which relates the vectors A inc of components of incident field to the vector A di f f of diffracted field we have

A pq di f f = i 2 ∑ pq T pq,p q I ψ S p q [A inc ], (30)
Now, interactions between distinct particles dispersed inside a multilayer can be studied using a generalized form of the translation matrix as introduced by Stout et al. [START_REF] Stout | A transfer matrix approach to local field calculations in multiple scattering problems[END_REF] to express the field generated by a source inside a particle in term of components of incident field on another particle.

Broadband absorber design

Now the general theoretical framework needed to describe the cooperative electromagnetic interaction inside a network of optical resonators we discuss in this paragraph how to use it to design targeted optical properties. To start with this objective ans show the strong potential of cooperative interactions we first consider the simple geometric configurations as illustrated in Fig. 2, that is single and binary metallic [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF] particle arrays dispersed in regular hexagonal lattices of side length d and compare their absorption spectra with that ones of isolated particles and of homogeneous metallic film. All lattices are immersed in a transparent material of refractive index n h = 1.5 and are maintained at a distance h = 100nm from the surface. The results plotted in Fig. 2 clearly show that the resonance peaks in single particle lattices are essentialy centered at the resonance frequency of free particles. On the other hand the absorption spectrum of nanoparticles lattices is much broder and does not simply consist in a superposition of single particle spectra. Moreover, we see that the cooperative interactions allow increasing the absorption even in diluted lattices where the filling factor f is smaller than 3%. Finally, the comparison of the overall absorption of nanoparticle lattices with that of simple metallic films with a thickness defined, using the effective medium theory, from the nanoparticle filling factors points out the prime importance of cooperative effects to magnify the absorption level.

In binary lattices, new configurationnal resonances add up to the resonances of single lattices and naturally enlarge the absorption spectrum.

In light of these results we can introduce a rational design of cooperative electromagnetic interactions to optimze the optical properties of a composite structure made with a distribution of nanoparticles. For this purpose we present the inverse design of a broadband absorber in the visible range [START_REF] Aydin | Broadband polarization-independent resonant light absorption using ultrathin plasmonic super absorbers[END_REF][START_REF] Narimanov | Optical black hole: broadband omnidirectional ligh absorber[END_REF][START_REF] Aubry | Plasmonic light-harvesting devices over the whole visible spectrum[END_REF][START_REF] Sergeant | Design of wide-angle solar selective absorbers using aperiodic metal-dielectric stacks[END_REF] made with a n-ary array metallic spherical nanoparticles. A n-ary lattice is defined from a unit cell C of a two dimensional paving with a certain thickness (see Fig. 3). In the unit cell of a lattice we consider a set of n vectors r i and n positive reals R i that represent the location of particles center and the radius of particles, respectively. To avoid the particle interpenetration these vectors must satisfy to the supplemental constraint | r i -r j |> R i + R j . To design the n-ary lattice in order to maximize its overalll absorption we have to explore the large and complex space of all possible configurations. To do that we employ a genetic algorithm (GA) [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF] which is a stochastic global optimisation method that is based on natural selection rules in a similar way to the Darwin's theory of evolution. Evolutionary optimization has been yet successfully applied in numerous fields of optics [START_REF] Feichtner | Evolutionary optimization of optical antennas[END_REF][START_REF] Drevillon | Ab initio design of coherent thermal sources[END_REF].

Basically, a GA uses an initial population (Fig. 3) of typically few hundreds of structures also called individus which are randomly generated in size and position. For each individu we calculate the fitness parameter which is here the mean absorption (for a given polarization) A = (λ maxλ min ) -1 λ max λ min A(λ )dλ over the spectral range [λ min ; λ max ] where we want to increase the absorption. The monochromatic absorption A(λ ) at a given wavelength is simply given by the the sum of power dissipated inside the particles of the unit cell normalized by the incident flux φ inc (λ ) on its surface S that is

A(λ ) = ∑ m∈Cell P m (λ ) S φ inc (λ ) . (31)
The GA consists in maximizing the fitness function of structures (i.e. A → max). To do so, we select 90% of the highest fitness as future parents for the next generation of selecting process. Those parents are linear crossed and the new 'child' generation is completed by new individual structures (randomly generated) to keep the same total number of lattices for any generation. To By taking into account the multipolar interactions until the second order (i.e. quadrupolar interactions) we see that the level of aborption becomes close to one. The comparaison of these results with full electromagnetic simulations based on the finite element method shows that the higher order multipole moments do not contribute significantly to the overall absorption. The absorption enhancement can be understood by examinating the electromagnetic cooperative effects inside the structure. These effects are highlighted in Fig. 5 at two different wavelengths by plotting the local losses inside the gold (resp. silver) nanoparticles within the optimized binary lattice in presence or without silver (resp. gold) nanoparticles. At λ = 550nm, that is to say, at the resonance of gold particles (which corresponds to the region where we observe in Fig. 4 an important bump in the absorption spectrum when it is calculated in the dipolar approximation) we see that the presence of silver nanoparticles enhance by 20% the losses inside the gold particles. Similarly, at λ = 650nm, the gold particles enhance by a factor of 60% the dissipation inside the silver particles. However, because of the weakness of intrinsic losses inside isolated Ag particles this cooperative effect is not sufficiently important to increase the overall absorption of the structure. At low wavelength we have checked (not plotted in Fig. 5) that the high absorption levels as shown in Fig. 4 results from cooperative effects between the gold particles themselves. The silver particles do not play any role in the exaltation mechanism. Interestingly, the numerical simulations have shown also that the cooperative effects are not very sensitive to the presence of disorder. In Fig. 6 we show, by disturbing the optimal structure with a -random perturbation of particles locations by a maximum displacement of 20nm, that the discrepancy between the optimal structure and the perturbed ones, given by the mean square error ζ = [ λ max λ min (A(λ ) -A opt (λ )) 2 dλ ] 1/2 , remains small. For some realization a broadening of spectrum can be observed around 650 nm. This effect can attributed to the presence of new modes supported which give rise to new channels for dissipating light energy within the structure. However, the detailed study of random structures goes far beyond the scope of the present work and it will be carried out in a future work.

Conclusion

In conclusion, we have proposed a general method for engineering the cooperative electromagnetic interactions in resonators networks both in the dipolar approximation and for arbitrary multipolar orders. Our results have demonstrated the strong potential of these interactions to tailor the optical properties spectrum. We believe that this approach opens the way to a rational design of metamaterials and it could find broad applications in various fields of applied physics, as for instance, in the domain of photovoltaic energy conversion for the conception of more efficient solar cells, in optical information treatement for the design of quantum information systems and, according to the reciprocity principle [START_REF] Landau | Electromagnetics of Continuous Media[END_REF], in light extraction technologies to Fig. 6. Impact of disorder on the light spectrum at normal incidence in a binary Au-Ag lattice.The spatial location of particles is randomly perturbated by a displacement of 20nm. The red ciurve corresponds to the spectrum (in polarization TM at nomrla incidence) of the optimized structure and the dotted blue curve is the spectrum of a particular random realization (results in polarization TE, not plotted here are similar). The dashed area shows the maximum and minimm values of absorption spectrum of different random realizations. The histogram shows the discrepancy with the optimal fintess for different realizations of the structure. The red line on the histogram shows the mean error with respect to the number of realizations.The disorder is mimicked by using pseudoperiodic particle array with sufficiently large unit cells.

improve the performances of light emitting diodes.

Fig. 1 .

 1 Fig.1. Multiple light scattering interactions in a set of subwavelength plasmonic structures embeded in a transparent host material of refractive index n h . In the dipolar approximation each object is replaced by both a dipolar electric moment and a magnetic moment. The external field felt by each object decomposes into (1) the incident field, (2) the field radiated by the other objects and (3) the auto-induced field which comes from the interface after being emitted by the object itself. All dipoles radiate (4) in their surrounding.

Fig. 3 .

 3 Fig. 3. Evolutionary algorithm to optimize a n-ary lattice. (a) A random population of lattices (a physical view of an unit cell is plotted on the left) is randomly generated. (b) The best individus basd on the fitness function are selected as parents for the crossing over. (c) The next generation is created by linear crossing and completed by new individus (d) to keep the total population constant. (e) Mutations are aaplied on a few number of individus (typically 5%) in the current generation.

Fig. 4 .

 4 Fig. 4. Light absorption spectrum at normal incidence a binary Au-Ag lattice (red dashed curve) optimized by GA by taking into account all multipolar interactions until the second order (quadrupoles) and of a multilayer based on Au-Ag films of thickness defined with the filling factor in nanoparticles (i.e. effective medium theory). Circles curve shows the result obtained by solving the Maxwell's equations with a finite element method.

Fig. 5 .

 5 Fig. 5. Local losses at λ = 550nm in the particles of a gold lattice (a) with the same geometric parameters as in the optimized structure. Losses ℑ(ε)|E SG | 2 in the single particle lattice are normalized by the maximum loss. In (b) we show the normalized difference ℑ(ε)|E DG | 2 -ℑ(ε)|E SG | 2 of losses inside Au particles in presence and without Ag particles (white regions). Analogously, in (c) and (d) the cooperative effect induces by the presence of Au particles on the dissipation in the Ag particles is shown at λ = 650nm.

Acknowledgments

J.-P. H. acknowledges discussions with J.J. Greffet. P. B.-A. gratefully acknowledges the support of Total news energies.