Chao Xu

Vineeth S Varma

Min Sheng

Tony Q S Quek

Wireless Service Provider Selection and Bandwidth Resource Allocation in Multi-tier HCNs

In this work, the inter-linked problems of bandwidth allocation by wireless service providers (WSPs) and WSP selection by users in heterogeneous cellular networks (HCNs) are addressed. This work employs a multi-hierarchical game framework. While the interaction between users are modelled using evolutionary game theory, the interactions between competing WSP's are modelled as a non-cooperative spectrum bandwidth allocation game (N-BAG). Finally the interaction between the WSPs and the users is modelled as a multi-leader multi-follower Stackelberg game. Moreover, for the formulated evolutionary game, the existence and uniqueness of the evolutionary equilibrium (EE) was investigated. Conditioned on the obtained EE, the existence of a Nash equilibrium (NE) for the proposed N-BAG was proven and an offline algorithm to achieve the equilibrium state was proposed. Finally, simulation results verify the validity of the analysis and demonstrate that a unique NE would be achieved by the HCNs adopting the proposed scheme.

I. INTRODUCTION

Heterogeneous cellular networks (HCNs) are regarded as a promising approach to cater for the increased demand of high-data-rate services and requirement of ubiquitous access [START_REF] Ghosh | Heterogeneous cellular networks: From theory to practice[END_REF]- [START_REF] Xu | Distributed subchannel allocation for interference mitigation in OFDMA femtocells: A utility-based learning approach[END_REF]. In this light, HCNs is a novel network structure for the future. Particularly, in contrast to the traditional pure macro cellular network, HCNs would consist of high power macro base stations (MBSs) as well as various classes of low power nodes, which may include micro BSs (often called eNBs), femto BSs (FBSs), and possibly future radiating elements [START_REF] Andrews | Femtocells: Past, present, and future[END_REF].

Several works on game theory focused on HCN based scenarios can be found in literature [START_REF] Yang | Backhaul-constrained small cell networks: Refunding and qos provisioning[END_REF]- [START_REF] Jiang | Optimal pricing strategy for operators in cognitive femtocell networks[END_REF]. The major difference between the above mentioned works and ours the type of HCN considered in these works. Specifically, different tiers are deployed by different wireless service providers (WSPs) and each WSP only provides one service to users. Hence, their main concerns are spectrum leasing among different WSPs. This work however, jointly studies the WSP selection by users and spectrum allocation of WSPs.

In this work, the interactions between the users in an environment where HCNs are deployed by several WSPs is studied and the population share for each WSP is predicted using results from evolutionary game theory. The WSPs assume that the user population stabilizes at the EE and competes for a higher revenue. Each WSP determines the bandwidth allocation for each service while the users decide which WSP to subscribe to. Naturally, as every user and WSP are independent decision makers, this results in a situation best modelled through game theory.

II. SYSTEM MODEL

In the following sections, we consider an area of interest where a total of K different service would be provided to users (e.g., the macrocell service, picocell service and femtocell service) by N different WSPs having deployed N different HCNs. Additionally, for each HCN n, there are K n tiers, each of which has the BSs of a particular class, such as MBSs or FBSs. Note that the BSs in different tiers may be differ among each other due to differences in the transmit power, coverage area and spatial density [START_REF] Dhillon | Modeling and analysis of K-tier downlink heterogeneous cellular networks[END_REF], [START_REF] Soh | Energy efficient heterogeneous cellular networks[END_REF].

A. Network deployment model

For notational simplicity, we denote the set of WSPs as N = {1, 2, • • • , N } and meanwhile, for each WSP n, we denote the set of the K n tiers of networks as K n . Without loss of generality, we consider that K n ⊆ K = {1, 2 • • • K}, which means different WSPs may like to provide different services to users. For instance, some WSPs may provide macrocell service, picocell service and femtocell service, but some WSPs may only provide macrocell service.

The HCNs are operating on the orthogonal spectrum and each WSP n has to properly allocate the available spectrum bandwidth B n to different tiers, i.e., the WSP deployment framework is adopted [START_REF] Lin | Macro-femto heterogeneous network deployment and management: From business models to technical solutions[END_REF]. For the HCN operated by WSP n, we denote the spectrum bandwidth allocated to tier k by B n,k and the spectrum allocation profile by

B n = (B n,1 , B n,2 , • • • , B n,Kn ).
Therefore, there is no intertier interference, only intra-tier interference. Meanwhile, each user can independently decide to get the corresponding service from any WSP they prefer. For each user, let us denote the price charged by WSP n for service k as C n,k which is the fixed access fee having units $ per unit of time. For each service k, the WSPs supporting this service are given by N k .

Finally, we assume that the user population density is given by λ u k . Here, each service is associated to a certain user density. This is explained by the preferences of different kinds of users to different services. For example a highly mobile user population will prefer the macro-cell service in order to avoid regular hand-overs. Of course in practice, the same user can jump from one population set to another, but on an average the population of each set of user types can be assumed to be a constant. Each user in a certain population type, picks one out of the N WSPs to get his service. This choice is determined by the data rates and prices offered by the WSP, and results

B. Hierarchical game framework

For each WSP n, a key problem becomes how to accurately estimate the equilibrium behavior of users and thus make an informed decision. Here, we consider that all the WSPs and users are selfish, i.e., each individual user wants to maximize his payoff by selecting the best WSP and meanwhile, each WSP wants to attract more users and further improve his own revenue by properly allocating the spectrum resource to different tiers. To be more specific, after selecting the service provided by a particular WSP, each user will obtain some payoff determined by both the service price and expected transmission rate. Intuitively, the expected rate of each user will be determined by both the number of users choosing the same network (i.e., the effects of other users' actions) and the available spectrum bandwidth of the network (i.e., the effect of WSPs' actions). Similarly, to attract more users, each WSP also need to make his own decision on the spectrum partition by fully considering the behavior of other WSPs and the WSP selection of users. The scenario under investigation therefore allows two levels of competition and meanwhile one cyclic dependency as described above. 1 These competitions and cyclic dependency can be illustrated with the hierarchical framework shown Fig. 1, where three different game formulation are proposed to investigate the concerned problem. Specifically, in the lower level we formulate the WSP selection of users as an evolutionary game, with which the long-term equilibrium behavior of users can be captured. Meanwhile, in the upper level we have formulated a non-cooperative game to depict the competition among WSPs. In addition, the cyclic dependency interaction between the users and WSPs is modeled by applying the multi-leader multifollower Stackelberg game theory, where the WSPs are the 1 The definitions for both the user's payoff function and WSP's revenue will be sequentially given in the following section. leaders and the populations of users are the followers.

III. EQUILIBRIUM ANALYSIS

The formulated Stackelberg game can be solved with backward induction, i.e., the leaders (WSPs) can make their responses (i.e., B * ) by fully considering the followers' (users') equilibrium reactions (i.e., x * (B * )). The mapping x * (B * ) denotes the equilibrium behaviour of users, when given the WSPs' strategy profile B * . Here,

x = (x 1 , x 2 , • • • , x K ) where ∀k ∈ K , x k = (x a k ,k ; a k ∈ N k )
represents the population states. In this section, the evolutionary game modelling the user reaction to the WSP bandwidth allocation is studied first. Following which, the WSP's compete for a higher revenue assuming the user population to stabilize at the proposed EE based on their bandwidth allocation.

A. Evolutionary equilibrium analysis

For the WSP selection of users, the formulated evolutionary game can be formally depicted as follows:

• Population: The users that request for the various services. The density of these usrrs is denoted by λ u k .

• Strategy space: Each user in the k-th population can decide to get the service from which WSP. Therefore, his strategy space can be depicted with the set

N k = {n |n ∈ N , k ∈ K n } . (1)
• Population state: Given by x • Payoff function: The payoff is used to quantify the satisfaction level or fitness of a user adopting a strategy when given the population state. Mathematically speaking, the payoff function of an individual player choosing strategy a k can be considered as a mapping π a k : x → R, a k ∈ N k , ∀k. We define the payoff function of each user as:

π a k ,k = ln 1 + β a k ,k x a k ,k (2)
where

β a k ,k = Ra k ,k B a k ,k ln 2λ u k C a k ,k and R a k ,k
is the expected rate per bandwidth allocated for a user on WSP n accessing the network through tier k2 . We note that, from the economic perspective, the logarithmic payoff function π = ln (1 + y) can capture both the user's non-satiation, i.e., dπ dy > 0, and risk aversion, i.e., dπ 2 dy 2 < 0, properties. Hence, such a function has been widely used in wireless resource allocation to evaluate satisfaction degree of an agent [START_REF] Lin | Joint spectrum partition and user association in multitier heterogeneous networks[END_REF]. It is worth noticing that although only the logarithmic function is considered here, the following analysis also holds when other monotonically increasing functions are adopted. To depict the dynamically learning behavior of users, the replicator dynamic is adopted in this paper, which can be expressed as follows:

ẋa k = δx a k (π a k (x k ) -πk (x k )) (3) = δx a k   π a k - a k ∈A k x a k π a k (x k )   , ∀k ∈ K
where δ > 0 is the rate of strategy adaptation and π (x) is the average payoff in the population. Under this dynamic in each population, the percentage growth rate of the population share of each strategy is proportional to the excess of the strategy's payoff over the population's average payoff. This dynamic could be interpreted biologically as a model of natural selection, and economically as a model of imitation [START_REF] Cressman | Evolutionary dynamics and extensive form games[END_REF]. As shown in the above equation, the dynamics of service selection can be described with K first-order differential equations. Hence, in order to investigate the equilibrium behaviour of users, we consider the solution of this evolutionary dynamic as the fixed point of the differential equations, which is termed as the evolutionary equilibrium (EE) represented by x * . When the population state of the game is at EE, then no user would like to change his strategy.

Based on the support, the EE can be mainly divided into two classes, i.e., the boundary EE and interior EE. Particularly, let us consider a EE x * and denote its support as supp

(x * ) = (a k , k) a k ∈ A k , k ∈ K, x * a k > 0, x * a k ∈ x * . If supp (x * ) = (A k , k) k∈K we say x * is an interior EE.
Otherwise, we term it as a boundary EE. Generally, the class of the achieved EE for each replicator dynamic is determined by the initial state x(0), i.e., supp (x(0)) = supp (x * ). In the concerned problem, we note for the applied replicator dynamic the boundary EE is not stable since any small perturbation will make the system deviate from this boundary state. Hence, in the rest of this paper, we just focus on the interior EE. Furthermore, we have the following Theorem.

Theorem 1: For the replicator dynamic in (3), the formulated WSP selection game has an unique interior EE state. Let us denote this EE by x * , and each element x * a k in x * can be expressed as

x * a k = η a k B a k ,k a k ∈A K η a k B a k ,k , ∀a k ∈ A k , ∀k ∈ K (4)
where

η a k ,k = λ a k ,k R0 a k ,k C a k ,k
. Intuitively, η a k ,k can be seen as the cost performance per unit bandwidth of the k-th service provided by WSP a k .

Proof: Substituting (2) into (3), we can equally transform the replicator dynamics of the formulated evolutionary service selection game as [START_REF] Andrews | Femtocells: Past, present, and future[END_REF]. Based on (5), we note that to get the EE x * we should make sure that ẋ *

a k = 0, ∀a k ∈ A k , ∀k ∈ K. Meanwhile, recalling that x * a k > 0 (∀a k ∈ A k , ∀k ∈ K) and a k ∈A k x a k = 1, we have x * a k = β a k ,k R0 a k ,k β a k ,k R0 a k ,k x * a k , ∀a k , a k ∈ A k (6)
and

a k ∈A k β a k ,k R0 a k ,k β a k ,k R0 a k ,k x * a k = 1 (7)
Further, we can get the expression of the EE as

x * a k = β a k ,k R0 a k ,k a k ∈A k β a k ,k R0 a k ,k = λ a k ,k B a k ,k C a k ,k R0 a k ,k a k ∈A k λ a k ,k B a k ,k C a k ,k R0 a k ,k (8)
where ∀a k ∈ A k , ∀k ∈ K.

Then, replacing the term

λ a k ,k B a k ,k C a k ,k
by η a k ,k we conclude the proof.

Remark 1: From an economic perspective, elements in x * denote the expected long-term equilibrium of market shares of different services. For instance, the expected market share of each WSP n can be represented as

k∈K a k ∈N k ,a k =n
x * a k . Therefore, with such statistical results, WSPs can get some insights about the effect of bandwidth allocation on users' WSP selection and further, appropriately deploy the HCNs.

B. Spectrum Bandwidth Allocation Game

In order to obtain a higher market share and revenue, a WSP must properly allocate the available spectrum. Since the spectrum allocation profile may not be frequently adjusted, each WSP should make his decision based on the equilibrium state of the service selection of users, i.e., the EE. 3 Hence, we define the utility function for the WSP as follows

U n (B n , B -n ) = k∈Kn U n,k = k∈Kn λ u k x * n,k C n,k (9)
= k∈Kn ϕ n,k B n,k m∈N k η m,k B m,k = k∈Kn ϕ n,k B n,k m∈N k ,m =n η m,k B m,k + η n,k B n,k
where

B -n = (B 1 , • • • , B n-1 , B n+1 , • • • , B N ) denotes bandwidth allocation profiles of WSPs other than WSP n. Where ϕ n,k = λ u k η n,k C n,k , ∀k ∈ K n , ∀n ∈ N .
We also impose the constraint B n,k ≥ b on each WSP's bandwidth allocation strategy, where b > 0 is a positive constant. 4As shown in ( 9), the utility of each WSP is determined by both its own bandwidth allocation profile and those of other WSPs. In contrast to the study in the previous section, we assume that each WSP is fully capable of acquiring all the information on HCNs deployed by its competitors (other operators) i.e., the expected average rate and pricing for each service. Conditioned on this, we can formulate the problem as

ẋa k = δx a k ln 1 + β a k ,k Ra k ,k,u x a k - a k ∈A k x a k ln 1 + β a k ,k Ra k ,k,u x a k , ∀a k ∈ A k , ∀k ∈ N (5)
a non-cooperative spectrum Bandwidth Allocation Game (N-BAG) and hereafter, the terms WSP and player would be used interchangeably in this section. Definition 1: N-BAG: This game can be represented by the tuple

G = Γ N , (B n ) n∈N , (U n ) n∈N . (10)
Here, N denotes the set of players which is identical to the WSP set. For each player n, its strategy space B n is defined as the set of available spectrum bandwidth allocation profiles

B n = B n B n = (B n,k ) k∈Kn , B n,k > b, k∈Kn B n,k = B n (11)
Based on ( 9), when given a strategy profile

(B n ) n∈N = (B 1 , B 2 , • • • , B N ) ∈ (B n ) n∈N , (12)
the utility function of each player n, n ∈ N , can be expressed as

U n (B n , B -n ) = k∈Kn ϕ n,k B n,k χ n,k (B -n ) + η n,k B n,k (13)
where the term

χ n,k (B -n ) = m∈N k ,m =n η m,k B m,k depicts
the competition caused to player n from other players. For a non-cooperative game, NE is the standard solution standing for the equilibrium state, under which no player can unilaterally improve its own utility by choosing a different strategy [START_REF] Myerson | Game theory: analysis of conflict[END_REF]. Accordingly, we have the following theorem about the existence of NE for N-BAG.

Theorem 2: For our formulated non-cooperative game N-BAG, there always admits at least one NE under which no WSP would like to unilaterally change his own bandwidth allocation strategy. Furthermore, when given other WSPs' strategy profile B -n , the best response of WSP n,

B * n = B * n,k k∈Kn
, is unique and can be expressed as

B * n,k = max b, λ u k C n,k χ n,k (B -n ) υ * n η n,k - χ n,k (B -n ) η n,k . ( 14
)
where

k∈Kn max b, λ u k C n,k χ n,k (B -n ) υ * n η n,k - χ n,k (B -n ) η n,k = B n . ( 15
)
Proof: First, we will prove the existence of NE for the formulated game N-BAG. After some algebraic computations, the utility of WSP n is transformed into:

U n (B n , B -n ) = τ n - k∈Kn λ u k C n,k χ n,k (B -n ) χ n,k (B -n ) + η n,k B n,k (16)
where

τ n = k∈Kn ϕ n,k η n,k = k∈Kn λ u k C n,k .
It is straightforward to verify that the n-th player's utility U n is concave over his strategy set B n for each fixed B -n and meanwhile, the strategy space (B n ) n∈N is a closed and bounded convex set [START_REF] Boyd | Convex optimization[END_REF]. Therefore, we can know that G is a concave game which always admits at least one NE [START_REF] Rosen | Existence and uniqueness of equilibrium points for concave N-person games[END_REF]. We note that this proposition can be proved through three steps. First, when given a strategy profile, the best response of each player should be compact and convex. This is true due to the continuity and concavity of the utility function. Then, it can be proved that that the best response function Ω : (B n ) n∈N → (B n ) n∈N is upper semi-continuous. Finally, by applying Kakutani fixed point theorem, it can be proved that there exists a fixed point

Ω((B * n ) n∈N ) = (B * n ) n∈N .
To be more specific, this fixed point (B * n ) n∈N is an NE for the considered concave game. Next, we will derive the expression of the best response of WSP n when given the strategy profile of other WSPs B -n . In this case, for each player n, selecting the strategy maximizing his own utility is equivalent to solving the following optimal spectrum allocation problem

P : min Bn k∈Kn λ u k C n,k χ n,k (B -n ) χ n,k (B -n ) + η n,k B n,k (17)
s.t. B n,k ≥ b, ∀k ∈ K n , (18) k∈Kn
B n,k = B n , (19)
which is a convex problem. For problem P, the Lagrangian can be written as

L n = k∈Kn λ u k C n,k χ n,k (B -n ) χ n,k (B -n ) + η n,k B n,k (20)
- k∈Kn θ n,k (B n,k -b) + υ n k∈Kn B n,k -B n
where θ n,k ≥ 0 and υ n are the Lagrange multipliers [START_REF] Boyd | Convex optimization[END_REF]. Moreover, the KKT conditions can be expressed as

∂Ln ∂B n,k = - λ u k C n,k χ n,k (B -n )η n,k (χ n,k (B -n ) + η n,k B n,k ) 2 -θ n,k + υ n = 0, B n,k ≥ b, k∈Kn B n,k = B n , θ n,k ≥ 0, θ n,k (B n,k -b) = 0, ∀k ∈ N k . (21)
Since problem P is a convex problem, the optimal variable B * n = B * n,k k∈Kn can be achieved by solving the KKT conditions (21), which could be further transformed as

υ n ≥ λ u k C n,k χ n,k (B-n)η n,k (χ n,k (B-n)+η n,k B n,k ) 2 , (B n,k -b) υ n - λ u k C n,k χ n,k (B-n)η n,k (χ n,k (B-n)+η n,k B n,k ) 2 = 0, B n,k -b ≥ 0, k∈Kn B n,k = B n , ∀k ∈ N k . (22)
To solve the above equations, two disjoint cases need to be considered. Firstly, if υ n <

λ u k C n,k χ n,k (B-n)η n,k (χ n,k (B-n)) 2 = λ u k C n,k η n,k
χ n,k (B-n) , the first three conditions in (22) only hold when

B n,k = λ u k C n,k χ n,k (B -n ) υ n η n,k - χ n,k (B -n ) η n,k ≥ b. (23)
or

B n,k = b, λ u k C n,k χ n,k (B -n ) υ n η n,k - χ n,k (B -n ) η n,k < b. ( 24
)
On the other hand, if

υ n ≥ λ u k C n,k η n,k
χ n,k (B-n) , then B n,k > b > 0 and the second condition in (22) can not be satisfied simultaneously, which further means B n,k can only equal to b. To this end, the optimal solution B * n = B * n,k k∈Kn can be expressed as

B * n,k = max b, λ u k C n,k χ n,k (B -n ) υ * n η n,k - χ n,k (B -n ) η n,k .
Meanwhile, recalling the constraint that k∈Kn B n,k = B n , we can get the optimal Lagrange multiplier υ * n by solving

k∈Kn max b, λ u k C n,k χ n,k (B -n ) υ * n η n,k - χ n,k (B -n ) η n,k = B n .
Finally, we will prove the uniqueness of B * n when given B -n . It should be noted that the left side of the above equation is a piece-wise-linear increasing function of 1 υ * n , hence the solution of this equation is unique when given B -n , which further means that the solution shown in ( 14) is also unique in the same case, i.e., given B -n . Based on the above, an iterative algorithm has been developed to obtain the NE of N-BAG, which is shown in Algorithm 1. Particularly, at the beginning of this algorithm, the starting point will be initialized based on the available bandwidth of each WSP, i.e., k∈Kn B n,k = B n , ∀n ∈ N . After that, the algorithm goes into a loop. At each iteration t, each WSP n calculates the optimal Lagrange multiplier υ * n with (15) and then, updates the best response B n (t) by applying [START_REF] Cressman | Evolutionary dynamics and extensive form games[END_REF]. Then, the loop will stop when the relative difference of the two optimal solutions obtained after the consecutive iterations (i.e., (B n (t -1)) n∈N and (B n (t)) n∈N ) is small enough. Mathematically speaking, the stopping criterion of Algorithm 1 can be expressed as Set t = t + 1 and B n (t) = B n (t -1) .

(B n (t)) n∈N -(B n (t -1)) n∈N (B n (t -1)) n∈N < ε (25)

4:

for n = 1 to N users do 5:

Calculate the optimal Lagrange multiplier υ * n with (15)and then, update the best response B n (t) by applying [START_REF] Cressman | Evolutionary dynamics and extensive form games[END_REF]. and the corresponding utility

U n,k (t) to 0, ∀k ∈ K n , ∀n ∈ N . 9: Set B * n = B n (t) and U * n = n∈N U n,k (t), ∀n ∈ N .
where • denotes a proper vector norm and ε is a small positive constant. Remark 2: It should be noted that when the initial point changes, the finally achieved solution may be different unless the NE of N-BAG is unique. However based on our extensive simulation results, we conjecture that the NE is unique and hence, the algorithm will converge to the same point irrespective of the starting allocation chosen.

Remark 3: From a mathematical perspective, note that this game and the corresponding IWF algorithm proposed is indeed quite similar to the well known IWFA for power allocation [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF]. However, a significant difference that makes the problem harder to prove uniqueness of the NE is that, the water level in each service k for this work are different unlike the uniform water level over channels seen in [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF]. Note that the developed algorithm is an offline algorithm. In addition, since we assume that each WSP can acquire the related information of the co-players, this algorithm can be independently implemented by each WSP and will converge to an NE. Since we observe from simulations that the NE is unique, all WSPs will naturally operate at the unique NE by adopting this algorithm.

IV. SIMULATION RESULTS

We have tested the uniqueness of the NE for an exhaustive set of parameters using Monte-Carlo simulations. In this section, we would like to study the efficiency of the NE. For this purpose, we compare it with:

• Global Optimum : In a hypothetical case where all WSPs are given instructions from a centralized agent, the WSPs choose the bandwidth allocation maximizing their sum utilities.

We note that this scheme is similar to the one proposed in [START_REF] Lin | Joint spectrum partition and user association in multitier heterogeneous networks[END_REF], which focuses on the single WSP scenario. Here, we consider there are two WSPs, each of which deploys a two-tier HCN 5 to provide services to users with the population density as (λ u 1 , λ u 2 ) = (20, 60) per km 2 . In addition, the available bandwidth of the three WSPs are set as (B 1 , B 2 ) = (10, 5) MHz, which means that the bandwidth obtained by each WPS is not the same and may sometimes happen in some practical cases. Based on the different deployment of HCNs, we consider that the price charged by WSPs are also different: (C 1,1 , C 1,2 ) = (2, 0.8) and (C 2,1 , C 2,2 ) = (0.8, 0.1). Also, we take R n,k ≈ 1.49nats/sec/Hz, ∀n ∈ {1, 2}, ∀k ∈ {1, 2}. The performance of these schemes in terms of the utility achieved by WSPs are compared in Fig. 2.

V. CONCLUSIONS

The bandwidth allocation of WSPs in a multi-tier HCN has been formulated as a non-cooperative game and the existence of Nash equilibrium (NE) has been proven. An algorithm used to achieve the NE has also been developed. Simulation results have shown the NE is unique. The resulting NE has also been compared to the global optimum and other strategies. These results can be used by competing WSPs to achieve the NE without having to run an online algorithm which could potentially take a long time to stabilize.

There are many interesting extensions for this work. One of the possible extensions is to investigate infrastructure sharing among WSPs, i.e., the WSPs co-operate to maximize their revenue while minimizing the cost of infrastructure.

Fig. 1 .

 1 Fig. 1. Illustration of the hierarchical game framework.

Algorithm 1 :

 1 Computation of NE (B * n ) n∈N . 1: Initialization: For ∀n ∈ N , initialize the starting point B n,k (0) satisfying k∈Kn B n,k = B n and set t = 0.

7 :

 7 until the stop criterion shown in (25) is satisfied. 8: For B n,k (t) being equal to b set its value (i.e., B n,k (t))

Fig. 2 .

 2 Fig. 2. The utility obtained by WSP 1 and WSP 2 when different schemes are adopted. The point indicated by "Our Alg. 1" converges is the NE.

R n,k is determined by the infrastructure, transmission power, noise and many such factors but for the purpose of this work, is treated as a constant which can be calculated using a stochastic geometry approach.

This consideration is mainly because in practice, each operator usually makes a long-term decision instead of changing his decision frequently[START_REF] Duan | Economic analysis of 4G upgrade timing[END_REF].

This b can be interpreted as a minimum bandwidth allocation required in order to deploy a certain service. b = 0 leads to a potentially undefined utility function when B n,k = 0, ∀n for any k. So from a mathematical point of view this constraint keeps the function continuous and easier to handle.

Note that considering a three or higher tier HCN with more WSPs is easily possible, however for the numerical part we consider the case of two for better presentation of the results.


Member IEEE Antoine Jay
email: antoine.jay@isae.fr

Member IEEE Mélanie Raine
email: melanie.raine@cea.fr

Member IEEE Nicolas Richard

Normand Mousseau
email: normand.mousseau@umontreal.ca

Member IEEE Vincent Goiffon

Anne Hémeryck
email: anne.hemeryck@laas.fr

Member, IEEE Pierre Magnan

Simulation of Single Particle Displacement Damage in Silicon -Part II: Generation and Long-Time Relaxation of Damage Structure

Keywords: Index terms-Silicon, Displacement Damage, Molecular Dynamics, kinetic Activation Relaxation Technique, Defects, Dark current, Random Telegraph signal. I

A statistical study of displacement cascades induced by silicon Primary Knock-on Atoms (PKA) in bulk silicon is performed by running a large number of molecular dynamics (MD) simulations. The choice of the PKA species and energy varying from 1 to 100 keV comes from a previous particle-matter simulation [1]. The electronic stopping power missing in standard MD simulations is here taken into account using the Two Temperature Model (TTM). This prevents from overestimating the number of created defects. The damaged atomic structures obtained after one nanosecond of MD simulation are not representative of what is observed in image sensors for example after several minutes. For this reason, the kinetic Activation Relaxation Technique (k-ART) is used in a second step, allowing to access longer simulation times of up to second. The obtained damaged structures can then be compared with experimental observations. Analysis reveals two possible links between the simulated structures and the measurements in solid-state image sensors. First, the cluster size distribution exhibits a shape similar to the measured exponential distribution of Dark Current (DC). Second, the temporal evolution of metastable atomic configurations resembles experimental DC-Random-Telegraph-Signals.

INTRODUCTION

hile displacement damage (DD) effects have been studied for a long time on various types of devices (diodes, bipolar transistors, image sensors...), a fundamental understanding of the link between the changes at the atomic structure level and the resulting electrical effect measured on devices remains to be achieved [START_REF] Srour | Displacement Damage Effects in Irradiated Semiconductor Devices[END_REF]. One of the key barriers limiting this understanding is the difficulty to replicate numerically the evolution of damaged structures on realistic observable timescales compatible with available experimental data. Indeed, molecular dynamical (MD) approaches often used to generate damage structures at the atomic level only treat timescales up to nanoseconds. This is far from the timescale of experimental data generally spanning from milliseconds to days. Simulations using such ab initio molecular dynamics have also been performed previously (see for example [START_REF] Beck | Disorder-recrystallization effects in low-energy beam-solid interactions[END_REF][START_REF] Beck | Quantum mechanical description of displacement damage formation[END_REF]) but this kind of simulation is too limited in time and size to describe realistic displacement cascades.

To solve this issue, a comprehensive simulation project is developed, combining different techniques, in order to simulate the succession of phenomena, from the interaction between the incident particle and the target material, to the resulting electrical effect measured at the device level. This approach is described in details in a dedicated paper corresponding to the first part of this study [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF], and summarized in the schematic of Fig. 1. Because of the complexity of this project and because each step can lead to its own interesting conclusions, we choose to dedicate one paper to each simulation step. The first step simulating the incident particle interaction with the target material using Monte Carlo (MC) methods is described in [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF]. This paper is devoted to the second step, corresponding to the generation of damage by a Primary Knock-On Atom (PKA) in the crystalline structure and the temporal evolution of the resulting defect structure. To achieve this, two different simulation tools are used.

First, MD simulations are performed to generate the damage structure and to observe its evolution up to nanoseconds. The novelty in this paper lies in the simulation of rather large PKA energies, compared to previous papers (see for example [START_REF] Otto | Defect characterization of low energy recoil events in silicon using classical molecular dynamics simulation[END_REF]), corresponding to realistic cases generated by the primary interaction and determined in the MC step in [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF]. Moreover, special care is given to the consideration of electronic stopping power in the MD simulations, using the Two-Temperature Model (TTM) introduced by Duffy et al. [START_REF] Duffy | Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations[END_REF]. Finally, these simulations are performed on a large number of PKA directions in order to be as realistic as possible and to construct a statistically significant database of damaged structures.

The second step detailed in this paper uses a new algorithm called the kinetic Activation Relaxation Technique (k-ART) [START_REF] Béland | Kinetic activation-relaxation technique[END_REF]. This technique allows the evolution of the damage structure at a fixed temperature on much longer timescales than usual studies, i.e. up to seconds and more after the damage generation. This approach thus avoids the use of thermally activated annealing, which can generate misleading reaction barriers and structures because of the artificially accelerated W Fig. 1 : Schematic of the global simulation project. The steps inside the scope of this paper are identified by red squares. structural relaxation. Here the k-ART algorithm has been used at 300 K, which is the temperature at the end of the MD step.

The final goal of this paper is to identify complex defect structures still present in the material after seconds and that could be associated with experimental measurements of Dark Current (DC) or Dark Current Random Telegraph Signal (DC-RTS) in solid-state image sensors [START_REF] Goiffon | Radiation Effects on CMOS Active Pixel UImage Sensors[END_REF].

The paper is organized as follows: in Sec. II, the two simulation steps and their computational details are described. In Sec. III, the MD-TTM results are given, the cascade properties and the generated defects are analyzed and discussed through a comparison with results coming from SRIM [9] and Geant4 [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF][START_REF] Allison | Geant4 Developments and Applications[END_REF] codes. Then, Sec. IV is devoted to the results obtained from the k-ART calculations. Results of the damage structure evolution after seconds are then analyzed, allowing the identification of different configurations of stable defects. Moreover, the analysis also reveals the presence of "flickering" defects, i.e. damage structures which switch between a limited number of configurations over time. These "flickers" are discussed in relation with the Dark Current phenomena measured in image sensors. Finally, Sec. V summarized the main conclusions and gives perspectives regarding the next step of our approach.

II. METHODOLOGY

A. Molecular Dynamics

MD simulations are performed with the LAMMPS code [START_REF] Plimton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF] to access the atomic scale of defect configurations generated by the PKA. The simulation supercell contains 1 million silicon atoms in their diamond-like crystalline form. They interact together via the Stillinger-Weber (SW) potential [START_REF] Stillinger | Computer simulation of local order in condensed phases of silicon[END_REF]. The simulated system evolves during 1 ns, with an increasing time step varying from 0.001 to 10 fs as the velocity of the atoms decreases. The structure does not evolve anymore during the last steps.

The Two-Temperature Model (TTM) [START_REF] Duffy | Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations[END_REF][START_REF]For documentation of the two-temperature model implemented in LAMMPS[END_REF] is used in order to take into account the electronic stopping power in addition to the force field governing the shocks between the atoms. This model provides a mechanism for energy exchange between radiation-excited carrier gas and atoms. When TTM is not used, as it is the case by default in standard MD, this leads to too high penetration ranges. This overestimation depends on the studied material and is approximately of 150% in bulk silicon. Indeed, if not included, all the kinetic energy of incident particles is lost through non-ionizing interactions, thus also overestimating the resulting damage. This TTM model has been mainly used for metals [START_REF] Duffy | Including the effects of electronic stopping and electron-ion interactions in radiation damage simulations[END_REF][START_REF] Chen | A semi-classical two-temperature model for ultrafast laser heating[END_REF] or insulators [START_REF] Toulemonde | Track creation in SiO 2 and BaFe 12 O 19 by swift heavy ions: a thermal spike description[END_REF][START_REF] Toulemonde | Transient thermal processes in heavy ion irradiation of crystalline inorganic insulators[END_REF]. It has also been used for semiconductors [START_REF] Chen | Numerical investigation of ultrashort laser damage in semiconductors[END_REF][START_REF] Osmani | Damage in crystalline silicon by swift heavy ion irradiation[END_REF] but, to our knowledge, never with temperature dependent parameters. The temperature dependency has been estimated through the electronic specific heat C e in eV/(e -.K) as follow:

𝐶 ! 𝑋 = 1.29.10 !! + -1.2899.10 !! -2.93276.10 !! 𝑋 + 2.29991.10 !! 𝑋 ! + 9.27036.10 !! 𝑋 ! 𝑒 ! !.!"!"#!$%! !
with X=T e /1000 where T e is the electronic temperature. C e is obtained from a fit of experimental data [START_REF] Glassbrenner | Thermal conductivity of silicon and germanium from 3k to melting point[END_REF] at low temperature and of a free electron gas The other parameters independent from temperature involved in the TTM model are:

• The electronic density is taken as ρ e =0.05 e -/Å 3 assuming that each atom can give one electron for the electronic thermal transport [START_REF] Toulemonde | Transient thermal processes in heavy ion irradiation of crystalline inorganic insulators[END_REF]. • The electronic thermal diffusivity: D e =20 000 Å 2 /ps [START_REF] Toulemonde | Transient thermal processes in heavy ion irradiation of crystalline inorganic insulators[END_REF] is the one of molten metal as the electronic temperatures obtained in the simulations are really high. • The electron-phonon friction coefficient γ p =1.741 g/(mol.ps) is the one optimized by Chettah et al. [START_REF] Chettah | Behavior of crystalline silicon under huge electronic excitations: A transient thermal spike description[END_REF] for TTM simulations.

• The electron stopping power coefficient: γ s =39. [START_REF] Lindhard | Energy loss in matter by fast particles of low charge[END_REF] g/(mol.ps) is obtained with SRIM tables in the Lindhard and Sharff model [START_REF] Lindhard | Energy loss in matter by fast particles of low charge[END_REF] using silicon PKA in bulk silicon. • The electronic stopping critical velocity v 0 =79.76Å/fs is taken as being two times the cohesive energy 4.63 eV [START_REF] Farid | Cohesive energies of crystals[END_REF]. • The grid of finite volume on which the equations of diffusion are solved in the electronic subsystem is 20*20*20.

At the beginning of the simulation, all the atoms have a random velocity corresponding to a temperature of 300 K. Then, to start the simulation, one atom, the PKA, is given an initial velocity corresponding to the wanted kinetic energy. The chosen energies (1 and 10 keV) and species (silicon atoms) come directly from a selection of relevant PKA in the outputs of the Monte Carlo (MC) calculations performed in [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF]. The use of this very large PKA energy for MD allows the simulation of realistic cases, representative of the kind of damages generated by incident neutrons or protons from hundreds of keV to hundreds of MeV. As detailed in the following, 10 keV PKA indeed leads to the creation of clusters of defects that are large enough to form amorphous pockets larger than 20 unit cells and that do not re-crystallize totally during the annealing step. Some 100 keV simulations are also performed to validate the choice of these representative energies, checking that larger incident PKA energies do not induce larger clusters.

A large number of simulations is performed (100) for each incident energy. Each of them can be distinguished by its PKA initial direction. The set of all these directions is equally distributed all over the space to optimize the sampling of the possible directions that can result from the particle-matter interaction. This permits to access good statistics on the resulting damaged structures.

Indeed, the penetration depth and the number of defects vary drastically with the PKA direction. As it is then impossible to draw conclusions from only one simulation, a statistic study based on numerous simulations is then mandatory. As an illustration, the penetration depths of the 100 simulations performed for a 10 keV PKA with different incident directions are given in Fig. 2. This histogram shows that each of the 100 incident directions produces a different penetration depth. For this reason, it is inaccurate to deduce the cascade properties from only one simulation. Furthermore, the standard deviation of this histogram increases with the PKA energy, justifying even more the need of statistical results for this kind of study. Hence, this MD-TTM approach has three main improvements: the use of TTM, of realistic PKA energies and of statistics on different incident directions. The result of this simulation step is a collection of 100 damaged atomic structures. Since this step only simulates the first nanosecond of the life of the structures, it is not relevant to directly compare them to the experimental observations that occur several minutes after the irradiation. The long time scale annealing needed to do this comparison can be done using the Activation Relaxation Technique in its kinetic format (k-ART).

B. k-ART

The kinetic Activation Relaxation Technique (k-ART) [START_REF] Béland | Kinetic activation-relaxation technique[END_REF] is a self-learning and off-lattice kinetic Monte Carlo algorithm based on the Activation-Relaxation Technique nouveau (ARTn) [START_REF] Malek | Dynamics of Lennard-Jones clusters: A characterization of the activation-relaxation technique[END_REF]. This method allows the exact treatment of elastic effects both for energy minima and saddle points in complex materials in order to follow the time evolution of large simulation boxes over timescales up to seconds or more. The principle is to go out of the potential well of the ongoing atomic configuration by searching an energetically more favorable one. The method is separated in two steps.

The first one consists in the exploration of all the atomic configurations of the energy landscape separated by only one maximum (i.e., the nearest valleys). An event is then defined through three atomic configurations: the starting one, the final one and the saddle point. The transition probability is calculated from their potential energy differences. Events are then classified based on the local topology which allows the study of disordered environments without any reference to a crystalline reference lattice. At each steps, all low energy barriers, associated with a probability of at least 0.01 per cent, are then reconstructed and fully relaxed to account for elastic deformation induced by the structural evolutions.

The second step follows the standard kinetic Monte Carlo algorithm, which computes the time step based on the rates of events for the specific structure and selects one move from the initial configuration to another. The catalog of possible configurations is then updated if necessary and the event lists reconstructed to account for the change in topology before applying the kinetic MC step again.

With its efficiency flicker handling algorithm and its use of local-force calculations, this algorithm has already proved its efficiency on the simulation of a number of complex materials including metallic alloys and relaxation of ion-implanted silicon with system sizes ranging from a few hundred to many tens of thousands of atoms [START_REF] Béland | Kinetic activation-relaxation technique[END_REF]. Similarly to the MD-TTM simulations, the result of these simulations is also a library of damaged atomic structures, but relaxed on timescales more similar to the experimental ones. Due to computational cost, only some of the 100 MD-TTM final configurations are annealed with the k-ART code (between 1 and 10 depending on the energies), chosen amongst the most common ones.

III. MD-TTM SIMULATIONS RESULTS

A. Steps of the cascade

As illustrated in Fig. 3, five main steps can be distinguished in the formation of the cascade. First (Fig. 3a), the PKA (in red in the figure) shocks the Secondary Knock-on Atoms (SKA, in blue) with high velocity. SKA are far from each other: on average 4 nm for a 10 keV PKA. Second (Fig. 3b), the SKA transmit their energy to the surrounding atoms still with a high velocity. These two steps correspond to the ballistic propagation and determine the width and depth of the cascade. Third (Fig. 3c), those atoms give their energy to create interstitial atoms and vacancies. They form amorphous pockets when close enough. Fourth (Fig. 3d), the elastic thermal wave propagates and moves the atoms with a displacement so small that they come back to their initial crystalline positions (this corresponds to transmitted energies lower than the 21 eV required to displace an atom in silicon). Fifth (Fig. 3e), with the residual energy, a large part of the isolated Frenkel pairs (~50%) recombine to reform the crystal.

B. Characterization of the defects along the cascade

Two main groups of defects are observed along the cascades: the point defects and the extended defects.

Regarding point defects, i.e. isolated vacancies and interstitial atoms, the remaining vacancies (V) at the end of the simulation generally correspond to the original positions of the PKA and the SKA that move far away from their initial crystalline site. They are then located at the beginning of the cascade. On the contrary, the PKA and SKA become interstitial atoms (I) that are mainly located at the extremity of the cascade, i.e. in the amorphous pockets. For this reason, those pockets or clusters have generally an excess of atoms, whereas simple isolated defects are generally atom deficient (i.e. related to vacancies). An example of final structure is given in Fig. 4.

Those clusters form extended defects whose sizes and shapes vary. In the following, we define the size of a cluster as the longest edge of the smallest rectangular parallelepiped that contains the entire cluster and we consider that two defects are in the same cluster if the distance between them is lower than two times the crystalline interatomic distance, i.e. 2x2.35=4.7 Å in silicon. This allows distinguishing three kinds of clusters: the small ones, the medium ones and the large ones. Each of them will have its own behavior during the annealing step. Small clusters are the only ones simple enough to be easily categorized: they are formed by a bi, a tri or a quadri-vacancy (2V, 3V, 4V), or by a bi a tri or a quadriinterstitial atom (2I, 3I, 4I). Some of them are presented in Fig. 11 and discussed in [START_REF] Trochet | Diffusion of point defects in crystalline silicon using the kinetic activation-relaxation technique method[END_REF]. The medium clusters involve from tens to hundreds of defects; they are not dense enough to be called amorphous pockets. The large clusters are formed by more than one thousand defects and seem very similar to amorphous silicon. Depending on the PKA energy, general trends can be drawn from this first MD step.

For 1 keV PKA simulations, the final configurations are generally composed of two to ten isolated point defects and small clusters, and of one or two medium clusters involving 30 to 200 defects. An example of these final configurations is given in Fig 4.

For 10 keV PKA simulations, the final configurations are generally composed of ten to fifty point defects and small clusters, 3 to 10 medium clusters, and one or two large clusters.

For 100 keV PKA simulations, the results are relatively similar to 10 keV PKA simulations, except that the number of clusters of each type is larger. The fact that the large clusters have nearly equivalent sizes compared to those generated by 10 keV PKA simulations is not surprising. In fact, in these large energy simulations, the mean free path of the SKA is larger than the largest size of the amorphous pockets generated by 10 keV PKA. For this reason, the clusters of two SKA with an energy superior to 10 keV cannot overlap to form a larger cluster. Hence, 10 keV seems to be a good compromise, leading to amorphous pockets which sizes are also representative of the ones created by incident PKA with energies larger than 10 keV. This validates the choice made from analysis of neutronsilicon MC simulations in [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF].

The cluster size distribution at the end of MD simulations is given in Fig. 5 for 100 simulations of 10 keV PKA. It has an exponential shape that resembles the exponential shape of the dark current (DC) distribution measured in CMOS image sensors [START_REF] Belloir | Dark Current Spectroscopy in Proton, Carbon and Neutron Irradiated Pinned Photodiode CMOS Image Sensors[END_REF]. While the size of clusters is expected to decrease in the following k-ART simulation step, the shape of the curve should remain similar, only shifted to smaller cluster sizes. This resemblance could be evidence of the instinctive assumption that a small DC value is due to the presence of small clusters in the silicon photo-diode, whereas a large DC value is due to the presence of large amorphous areas in the irradiated photo-diode. This assumption still needs to be proven, for example by the calculation of the electronic quantum states in the clusters band gap. This is the last part of our global simulation project presented in Fig. 1, involving ab-initio calculations, and will be the subject of future papers.

C. Cascade properties -comparison with Geant4 results

The histogram of the number of secondary knocked-on atoms (SKA) generated by a 10 keV PKA is given in Fig. 6. It is compared to the results obtained with Geant4 (for details on the Geant4 calculations, see [START_REF] Raine | Simulation of single particle displacement damage in silicon -part I: global approach and primary interaction simulation[END_REF]). The average value obtained with MD-TTM is slightly lower than the one obtained with Geant4, but the order of magnitude is the same. In fact, the Binary Collision Approximation used in Geant4 is only considered valid beyond 1 keV. On the contrary, 100 keV is the upper value that limits the domain in which MD simulations are feasible due to the exponential increase of computational time. Hence, the comparison can only be done in a small range of value (1 to 100 keV), and the difference between the two models decreases with increasing the PKA energy increases.

D. Limits of Molecular Dynamics

The evolution of the number of defects as a function of time is given in Fig. 7. It is compared to a value coming from the SRIM code in full cascade mode, often used as reference for this kind of calculation. It shows that, even with the addition of the TTM, the number of defects obtained at the end of the MD-TTM run (1 ns) is larger than the SRIM value. In Fig. 7, we represent the time evolution up to 100 ps, because the number of defect doesn't change between 100 ps and 1 ns in MD. This difference comes not only from the shortness of the time scale accessible with the MD technique, but also from the fact that this simulation method makes impossible the access to other configurations. The use of the k-ART algorithm is a way to address this issue, as detailed in the next section. Fig. 7: Temporal evolution of the number of defects during MD-TTM for an initial PKA of 10 keV. The red point corresponds to the final number of vacancies given by SRIM. As the simulation begins with a perfect cristal, the number of interstitial atoms and of vacancies is the same and is here named number of defects.

IV. k-ART RESULTS

In this section, starting from atomic configurations obtained at the end of the MD simulation, we first describe the global annealing behavior of the damaged structure on long timescale accessible with k-ART. The annealing behavior of defects is then analyzed depending on their size. Finally we underline the particular defect structures that could be linked to DC and DC-RTS.

A. Global behavior

The evolution of damaged structure can be roughly split in four main parts, each of them having its own time scale.

First, unstable atomic configurations are stabilized very rapidly (~1 ns) locally during the first microseconds of the simulation. Indeed, the energy barriers that need to be overcome for their annealing do not exceed 0.2 eV. These configurations are part of the medium and large clusters.

Second, Frenkel pairs anneal also very rapidly (in ~100 ns).

Third, mono-vacancies (V), mono-interstitials (I) and bi-interstitial (2I) diffuse, because of their low energy barriers until they join into an amorphous pocket or meet another defect either to form a small cluster or to recrystallize.

Fourth, after ~1 s, only small clusters and residual large clusters are still present around the cascade and main evolutions occur over longer timescales (one event per minute or longer). This is due to the fact that the system reaches a single potential well surrounded by high energy barriers (> 0.8 eV), corresponding to a meta-stable configuration. These long-changing events involve the global displacement of a large number of atoms that is typically larger than 15. All these processes lead to a drastic diminution of the number of defects present in the structure, as shown in Fig. 8. This number is then in good agreement with the SRIM value, which was not the case after MD-TTM simulations. Examples of configurations reached after several seconds of simulation are presented in Fig. 11.

B. Behavior of point defects and small clusters

The three most simple defects, the mono-vacancy (V), the mono-interstitial (I) and the bi-interstitial (2I) have a different behavior than other small clusters. Those three defects are the only ones to diffuse rapidly in the crystal structure, because of their low diffusion energy barriers. Indeed, we find energy barriers of 0.51, 0.78 and 0.28 eV for V, I and 2I respectively compared to 1.00, 1.32, 1.52, 1.66 and 1.10 eV for 2V, 3V, 4V, 3I and 4I respectively. These values are in perfect agreement with [START_REF] Trochet | Diffusion of point defects in crystalline silicon using the kinetic activation-relaxation technique method[END_REF]. However, their migration is not random if other defects are close enough to them, so that they feel the lattice distortion induced by those surrounding defects. This distortion leads to an agglomeration of the defects. This is due to the potential energy of the so created cluster that is lower than the sum of the potential energies of the separated defects. If the agglomerated cluster is composed of an interstitial atom and a vacancy, i.e. a Frenkel pair, it re-crystallizes in a few microseconds. Otherwise, it forms a small cluster (2V, 3V, 4V, 3I or 4I) that does not diffuse anymore.

For this reason, all of these 3 simple defects formed during the cascade are not present at the end of the annealing (> 1 s). This behavior is described by the blue arrows in Fig. 9. On the contrary, the number of small clusters increases during annealing since simple defects agglomerate. This number cannot actually decrease because the migration energy barrier of small clusters is too high to allow their diffusion. Moreover, their V or I overage make their whole re-crystallization impossible. As shown in Fig. 9, this behavior is the same for all defects structures whatever the energy of their original PKA (1 keV or 10 keV).

C. Behavior of medium versus large clusters

The larger the PKA energy, the larger the size of the amorphous pockets or clusters located at the extremity of the cascade obtained by MD. There is a limit, located between 1 and 10 keV, which separates two kinds of phenomenon.

First, when these clusters generated by MD are small enough (medium clusters), i.e. typically composed of less than 200 non-crystalline atoms, they can totally recrystallize with time. This annealing is indicated by the green arrows in Fig. 9. If the number of vacancies (nV) and of interstitial atoms (nI) in one cluster is not the same, this total re-crystallization is not possible. In this case, the cluster quasi-re-crystallizes with time, leading to a small cluster, typically formed by a bi, a tri, or a quadriinterstitial atom or vacancy as a function of the difference nV-nI. As this number rarely exceeds 3 or 4, small clusters larger than 4V or 4I are very unusual. Nevertheless, they are sometimes formed by the agglomeration of a resulting 4V or 4I and of a diffusing point defect.

Second, when the clusters generated at the end of MD are large enough (large clusters), typically when they are composed of more than 1000 non-crystalline atoms, their annealing takes a longer time and is rarely total. This annealing is indicated by the orange arrow in Fig. 9. This difference of behavior between medium and large clusters is due to the fact that the re-crystallization of large amorphous pockets occurs simultaneously at different places. This leads to a non-uniform repartition of V and I forbidding their total annealing. Thus, a large cluster becomes an agglomeration of many closely spaced small clusters that are hardly distinguishable by experimentalists as their separating distance is lower than 4.7 Å. For this reason, what we name an "agglomeration of small clusters" that we are able to differentiate after annealing, is commonly named an amorphous pocket by experimentalists. Despite this, the size of the large clusters drastically decreases during annealing. They become smaller and more stable amorphous pockets than after MD-TTM. In the 10 keV PKA cascade annealing given as an example in Fig. 9, the large cluster is formed by 1500 defects at the last step of MD (1 ns) and by only 300 defects at the last step of the k-ART algorithm (1 min).

D. Flickering structures

In most cases, a final configuration continues to evolve even after days of annealing. However, this evolution is not a total defect annealing, but a cyclic switch between a finite number of configurations. In fact, the final system is trapped into several potential wells, separated by small energy barriers. This wells form a basin that is surrounded by insurmountable high energy barriers as shown in Fig. 10a). As the system can thermally switch between all the configurations into the basin, we name this switching topology a "flicker". An average lifetime can be attributed to each of its atomic configurations as a function of the calculated energy barriers. It is then possible to evaluate their evolution as a function of time as shown in Fig. 10b). These evolutions seem very similar to the ones of the DC-RTS observed in CMOS Active Pixel Image Sensors (CMOS-APS) [START_REF] Béland | Kinetic activation-relaxation technique[END_REF][START_REF] Durnez | In-depth analysis on radiation induced multi-level dark current random telegraph signal in slicon solid state image sensors[END_REF]. The example of Fig. 10 is a small cluster: the 3I. The flicker switches between four configurations, two of them having tenth of microsecond lifetimes. This short lifetime is hardly observable experimentally. Thus this flicker is probably seen as a stable DC and not as a DC-RTS. This should be the same for all the small clusters. On the contrary, a cluster that still have hundreds of defects after annealing also forms a basin, containing more than thousands of metastable configurations. Among all of them, a dozen are surrounded by energy barrier superior to 0.8 eV, corresponding to an average lifetime of 10 s that is easily observable experimentally. This kind of structure could be the source of multi-levels DC-RTS.

E. Limits of k-ART

It is known from experimental data that transient annealing is very dependent on the excess carrier concentration present in bulk material or in a device [START_REF] Srour | Review of displacement damage effects in silicon devices[END_REF].The annealing time can also be decreased due to an excess of charge-states that may modify the defects mobility [START_REF] Gregory | Injection-stimulated vacancy reordering in ptype silicon at 76 K[END_REF].

To take these effects into account in the k-ART simulations and allow direct comparison with experimental data, nor our empirical potential (SW) neither those using a charge density (EAM) are adequate. In fact, the positions of the localized states in the gap, which result from a defect in the periodic crystalline structure, are directly affected by an excess or lack of charges. As a consequence, the defect's binding energies that govern its mobility are also modified. To evaluate them in a charged system, one then need to solve its Hamiltonian. This task requires ab initio calculations. This technique does not permit to treat sufficiently large systems due to the present computational cost. Dealing with interatomic potential is thus for the moment the best compromise that allows studying mesoscopic annealing including diffusion and displacement of a large number of atoms.

IV. CONCLUSION

To describe the displacement damage in silicon, we use an original approach coupling Molecular Dynamics and Two-Temperature Model followed by k-ART, using input data coming from the Monte Carlo particle-matter interaction code Geant4. The TTM allows to take into account the electronic stopping power in MD simulations. A large library of damage structures is then generated after 1 ns with MD-TTM. This large library of configurations is afterwards relaxed on long times with the k-ART algorithm, leading to a complete catalog of possible atomic moves with their appropriate attempt-frequencies and activation barriers. Different stable structures are identified, along with a possible source of DC-RTS, corresponding to the multi-levels flickers in amorphous pockets.

  (3/2*k b *T) above the Fermi temperature T F = E F /k b = 0.56 eV/k b =6 498 K where k b is the Boltzmann constant.

Fig. 2 :

 2 Fig. 2: Distribution of the penetration depths of the cascade (above) and of the PKA (below) for the 100 different incident directions for 10 keV PKA.

Fig. 3 :

 3 Fig. 3: The five characteristic steps occurring during the MD-TTM propagation of the cascade of defects. Red balls: positions of the PKA. Blue balls: positions of the SKA. Orange balls: positions of the atoms shocked with a high velocity but not by the PKA. Green balls: positions of the atoms with low velocities. Black balls: positions of the atoms involved in the thermal wave (movie available as supplementary material).

Fig. 4 :

 4 Fig. 4: Examples of clusters generated by a 1 keV PKA. Blue balls represent interstitial atoms, red balls vacancies and yellow balls silicon atoms that remain on a regular crystalline site.

Fig. 5 :

 5 Fig. 5: Distribution of the size of the clusters generated by 10 keV PKA.Histogram coming from the last step of the 100 MD-TTM simulations at 1 ns. Orange numbers give the average number of defects per cluster.

Fig. 6 :

 6 Fig. 6: Distribution of the number of SKA in a 10 keV PKA cascade formed by MD-TTM (in black) and by the binary collision approximation (Geant4 in red).

Fig. 8 :

 8 Fig.8: Temporal evolution of the number of interstitial atoms (= number of vacancies) during MD-TTM (black curve) for an initial PKA of 10 keV (same curve as in Fig.7) and during k-ART annealing (red curve). The red dot indicates the value given by the SRIM code.

Fig. 9

 9 Fig. 9 represents an example of the occurrence of each type of clusters at the end of the k-ART simulation (here after 1 minute) compared to the end of MD simulation (1 ns), for one direction of 1 and 10 keV PKA.

Fig. 9 :

 9 Fig. 9: Examples of the variation of the number of defects per clusters for a 1 keV PKA simulation (top) and for a 10 keV PKA simulation (bottom).Black histogram: last step of MD-TTM simulation (1 ns). Red histogram: last step of k-ART simulation (1 min). Blue arrows: agglomeration of point defects. Green arrows: medium clusters total annealing. Orange arrow: large cluster partial annealing.

Fig. 10 :

 10 Fig. 10: An example of simple flickering configuration: the tri-interstitial. a): the four configurations, their relative energy and the energy barrier. b): time spend in each configuration after 8 ms of annealing.

Fig. 11 :

 11 Fig. 11: Examples of point defects (I and V), of small clusters (2I, 2V, 3I, 3V) and of a large cluster obtained at the end of the annealing. All the clusters of defects have several possible configurations (flickers), and only the stable one are shown there.

ACKNOWLEDGMENT The authors would like to thank M. Trochet, M. Toutlemonde, J. Srour, P. Paillet, M. Gaillardin, L. Martin-Samos, S. Girard, JM. Belloir and C. Durnez for fruitful discussions. The computer simulations were performed using HPC resources from GENCI-CCRT (Grant 2016097474) and from CALMIP (Grant 1555).


Daniel Bonilla

Vineeth S Varma
email: vineethsvarma@gmail.com

Samson Lasaulce
email: lasaulce@lss.supelec.fr

Jamal Daafouz
email: jamal.daafouz@univ-lorraine.fr

Mounir Ghogho
email: m.ghogho@ieee.org

Daniel Bonilla Licea

Trajectory planning for energy-efficient robots with communications constraints

Keywords: Trajectory planning, vehicles, robots, communication constraints, energy-efficiency

published or not. The documents may come

INTRODUCTION

In wireless networks, it is generally assumed that the trajectory of a mobile terminal is chosen by the user. Under this assumption, the communication and motion aspects are obviously independent. However, there are more and more applications for which this classical paradigm needs to be questioned. For instance, when an unmanned aerial vehicle or a mobile robot (MR) has to collect data from a field of wireless sensors, it typically has to optimize its trajectory to receive the data correctly (see e.g., [START_REF] Wang | Wireless sensors in agriculture and food industry-Recent development and future perspective[END_REF] [START_REF] Murray | Recent research in cooperative control of multivehicle systems[END_REF] ). This situation has becomes more frequent and appears in the wide area of networked robotics [START_REF] Das | An efficient group communication protocol for mobile robots[END_REF], and more specifically for nano-networks [START_REF] Akyildiz | Nanonetworks: A new communication paradigm[END_REF]. The emergence of this new type of communications scenarios led the authors to state the problem described on this paper.

The problem under consideration is as follows. A mobile terminal has to move from a starting point to a target point while transmitting a certain amount of data along its trajectory. In addition this task has to be completed in a time T . The data has to be sent to a wireless base station (BS) which receives the signal with a signal-to-noise ratio (SNR) which depends on the distance between the mobile and the base. Therefore, the mobile has to choose a trajectory which allows the data to be uploaded successfully (which is made possible by having a sufficiently large SNR) and to manage the energy consumed for moving. To the authors' knowledge, one of the closest work to the work reported in the present paper is given by [START_REF] Ooi | 6. The figure shows how the buffer size and the vehicle velocity varies over time for the optimal control laws u * 2 (t)[END_REF]. Therein, a robot has to optimize its trajectory to minimize a cost which consists of the sum of the energy consumed by the transmitter and the motion energy. Another relevant paper is given by [START_REF] Bonilla Licea | Improving Radio Energy Harvesting in Robots Using Mobility Diversity[END_REF]; the goal of the robot is find a trajectory which allows wireless energy to be harvested. Then in [START_REF] Cruz | Towards optical wireless communications between micro unmanned aerial and ground systems[END_REF] the authors design a control law for a drone to follow a ground robot while maintaining a minimum data rate in an optical wireless communications link. And in [START_REF] Kantaros | Distributed communication-aware coverage control by mobile sensor networks[END_REF] the authors maximize the coverage area of a mobile sensor network while ensuring wireless communications between its members; in [START_REF] Meyer | Distributed Estimation with Information-Seeking Control in Agent Networks[END_REF] the authors consider a cooperative mobile sensor network and then design control laws so that at each iteration the sensor nodes gather a maximum amount of information. The setting of the present paper differs from these works by the fact that a different performance metric and different communication scenarios are considered but more importantly by the theoretical analysis to derive the optimal trajectories. In particular, the link with optics is completely new and opens for a large variety of extensions; this link is reminiscent to the work by where optimal packet routing through a wireless ad hoc network has been interpreted as the light propagation path [START_REF] Jacquet | Geometry of information propagation in massively dense ad hoc networks[END_REF].

Our task has two main objectives, transmit all the bits as soon as possible, and reaching the target spending as little energy as possible, corresponding to two separate optimality criteria. So the optimization of the trajectory for this task is a multi-objective optimization problem involving communications and robotics aspects. In this paper we will present different approaches for this multi-objective optimization problem, namely in sequential importance and by constructing a hybrid cost.

This paper is structured as follows. In section 2 we describe the communications and the vehicle's models; then in section 3 we derive the optimum trajectories for the problem considered in this paper. Then in section 4 we present simulation results to give some insight into the optimum trajectories obtained and then we conclude in section 5.

SYSTEM MODEL

We denote by p(t) the position of the vehicle at time t and the signal received by the base station (BS) by the vehicle is:

y(t, p(t)) = x(t) p(t) -p B α/2 2 + n(t) (1)
where x(t) is the signal transmitted by the vehicle with fixed power E[|x(t)| 2 ] = P ; p B is the position of the closest BS; α is the path loss coefficient and n(t) ∼ CN (0, σ 2 n ) is the zeromean additive white Gaussian noise at the BS with power σ 2 n . So the SNR observed by the BS when the vehicle is located at p(t) is:

SNR(p(t)) = P p(t) -p B α 2 σ 2 n ( 2
)
The content of the vehicle's buffer at time instant t is:

b(t) = N - t 0 R b (p(τ ))dτ + ( 3
)
where N is the number of bits initially stored in the buffer;

[a] + = a if a > 0 and [a] + = 0 if not; R b (•)
is the bit transmission rate at position p(t), which increases as the SNR increases. Although strictly speaking b(t) is a discrete variable (since it represent the number of bits) we will consider it to be a continuous variable for mathematical simplicity. In addition, we select p B = 0 for mathematical convenience without any loss in generality. Now, if the SNR is below a minimum value SNR min determined by the sensitivity of the BS then the vehicle will not transmit and so the bit rate will be zero. But if the SNR is at least SNR min or higher then we will consider two different policies for determining the bit rate:

1. Discrete bit rate. In this case there is a finite number of different bit rates available. This can be the case when M-QAM modulation is used and we control the size of the constellation to determine the bit rate. Assume that there are M possible rates achieved based on the distance (and consequently the SNR) and we select

R b (p(t)) = R m if ||p t || 2 ∈ [r m , r m+1 ) for all m ∈ {1, 2, . . . , M } where r M +1 = 0 and r 0 → ∞ re- sulting in R 0 = 0. In addition R 1 < R 2 < • • • < R M .
We will also define A m as the region in which the bit rate R m is used. So A 0 is the area that is not covered by the BS.

Continuous bit rate

In this case we will assume the bit rate R b (SN R(p(t))) to be a continuous increasing function of the SNR.

Finally, the position of the robot is given by:

p(t) = s + t 0 u(τ )dτ (4)
where u(τ ) is the velocity of the robot at time τ , and is constrained by

u(t) 2 ≤ v max (5)
where v max is the maximum velocity of the robot and the energy consumed for motion from time t a to t b is:

E m (t a , t b ) = t b ta u(τ ) 2 dτ. (6)

OPTIMUM TRAJECTORY PLANNING

The task to be completed in a time T by the vehicle has two objectives: transmit all the content of the buffer and then reach the goal point g. We would like the buffer to be emptied quickly and to reach the goal point g using little mechanical energy. In this section we will present two different approaches to this problem.

Sequential importance

The simplest way to present this multi-objective optimization is by first giving full importance to the optimality criteria associated with the transmission of all the N bits. And then, once all the bits are transmitted, give full importance to the optimality criteria associated with spending a minimum amount of energy reaching goal point g. Let us first define t b (u) as the time in which the buffer is emptied when using the control law u:

t b (u) = min{t ∈ [0, T ] | b(t) = 0}. (7)
Now, the multi-objective optimization problem mentioned above can be then written as the following coupled pair of optimization problems:

minimize u t b (u), s.t. p(t) = s + t 0 u(τ )dτ, u(t) 2 ≤ v max , b(t) = N - t 0 R b (p(τ ))dτ + , p(T ) = g. (8)
which formulates the first optimization problem. Now, the first two constraints in [START_REF] Kantaros | Distributed communication-aware coverage control by mobile sensor networks[END_REF] describe the dynamics of the vehicle and the maximum speed of the vehicle while the third constraint describes the dynamics of the buffer. The fourth constraint simply states that the vehicle must reach the goal point g at time T . The set of all control laws that solve ( 8) is denoted by U b , i.e., the set of all control laws that minimize the time taken to empty the buffer. Now, the second optimization problem considers the minimization of the energy used to reach the goal point g once the buffer is emptied:

minimize u(t)∈U b T 0 u(τ ) 2 dτ s.t. p(t) = s + t 0 u(τ )dτ, u(t) 2 ≤ v max . (9)
Note that ( 9) is minimized over U b so coupling this optimization problem with [START_REF] Kantaros | Distributed communication-aware coverage control by mobile sensor networks[END_REF]. We have to mention that in order to solve the multi-objective optimization problem we need to find an optimum control law that simultaneously solves both ( 8) and ( 9).

Before we try to solve both ( 8) and ( 9) we need to see if there exist a solution for the value of T . One way to do this is by solving a feasibility problem. In this case the feasibility problem consist in minimizing the time t a (u) in which the vehicle reaches the goal point with its buffer empty. Now, the feasibility problem is:

minimize u t a (u), s.t., p(t) = s + t 0 u(τ )dτ, u(t) 2 ≤ v max , b(t) = N - t 0 R b (p(τ ))dτ + , t a (u) = min{t ∈ R + |p(t a ) = g, b(t a ) = 0}. ( 10
)
Denote by u * a (t) the control law that solves [START_REF] Jacquet | Geometry of information propagation in massively dense ad hoc networks[END_REF] and t * a the minimum value of t a (u * a ). If t * a ≤ T then the multi-objective optimization problem defined by ( 8) and ( 9) has a solution. In addition, it is easy to see that if t * a = T then both the feasibility and the multi-objective optimization problems are equivalent. Now, let us define the time set T * where

T * = {t ∈ T |R b (p(t)) = max{R b (p(t))}} , (11)
We also define t 1 = min{T * }, t 2 = max{T * }. Then we have the following result.

Proposition 1. The trajectory satisfying [START_REF] Jacquet | Geometry of information propagation in massively dense ad hoc networks[END_REF] is such that its segments (s, p(t 1 )) and (p(t 2 ), g) are straight lines and the vehicle uses v max in these segments. When the rates are continuous, p(t 1 ) = p(t 2 ). In the discrete rate case, the robot moves in a straight line from p(t 1 ) to p(t 2 ) where

||p(t 1 )|| 2 = ||p(t 2 )|| 2 .
Therefore solving [START_REF] Jacquet | Geometry of information propagation in massively dense ad hoc networks[END_REF] is equivalent to finding the optimum points p(t 1 ) and p(t 2 ) (t 2 , t 1 will be automatically determined based on v max and N ). And as mentioned before once we solve [START_REF] Jacquet | Geometry of information propagation in massively dense ad hoc networks[END_REF] we can determine if the multi-objective optimization problem (8), ( 9) is feasible. If this criteria is satisfied, then we can use the following proposition to find the optimal path satisfying ( 8) and ( 9) by using the same T * and t 1 ,t 2 defined earlier.

Corollary 1. The paths satisfying ( 8) and ( 9) is such that its segments (s, p(t 1 )) and (p(t 2 ), g) are straight lines and the vehicle uses v max in the segment (s, p(t 1 )). Additionally, ||p(t 1 )|| 2 , ||p(t 1 )|| 2 satisfy the same properties as in proposition 1.

Special case: When T is considerably larger than t * a . It can be easily observed that to minimize t b (u(t)) in [START_REF] Kantaros | Distributed communication-aware coverage control by mobile sensor networks[END_REF] we must increase the bit rate R b (SN R(p(t))) as quick as possible. This is achieved by heading directly towards the BS. Therefore, in the continuous rate case, for t ∈ [0, t 1 ], the vehicle must go from s towards the BS in straight line at maximum velocity v max . Then the optimum control law for t ∈ [0, t 1 ] case is:

u * (t) = - s s v max (12)
while the robot moves in a straight line from t 1 to t 2 (for the discrete case), but with a velocity determined by N such that b(t 2 ) = 0 . It can also be verified that for optimizing [START_REF] Meyer | Distributed Estimation with Information-Seeking Control in Agent Networks[END_REF], for t > t * b , the control law is.

u * (t) = g -p(t 2 ) g -p(t 2 ) 2 g -p(t * b ) T -t * b . (13)
• If t * b < t 2 then the robot moves in a straight line from p(t 1 ) to g.

• Otherwise p(t 2 ) is such that, given t 2 , ||p(t 2 ) -g|| is minimized.
So now that we have solved the multi-objective optimization problem considered in this section we will consider in the next section another approach to optimize the same task.

Combined importance

In the previous section first we gave full importance to the optimality criteria associated with the transmission of the N bits and once this objective was fulfilled we gave full importance to the optimality criteria associated with reaching the goal point g. An alternative approach is to simultaneously consider both optimality criteria in the design of the trajectory and also give similar importance to both terms. So we the control law that produces such optimum trajectory can be obtained from solving:

minimize u T 0 u(τ ) 2 2 v 2 max dτ + T 0 b(τ ) N dτ s.t. p(t) = s + t 0 u(τ )dτ, u(t) 2 ≤ v max , p(T ) = g, b(t) = N - t 0 R b (SN R(p(τ )))dτ + , b(T ) = 0 (14)
The first term of the optimization target in (14) represents the mechanical energy used during the whole task while the second is a measure of how fast the buffer is emptied. The normalization of each term (by v 2 max for the first term and by N for the second term) is done to make both terms to lie in the same order of magnitude and consequently having similar importance in the optimization target in (14). Now, we will solve (14) for the discreet case. But we first define u m (t) as any control law that takes the vehicle through the regions {A k } m k=0 . The set of all control laws u m (t) will be denoted as U m and U = ∪ M m=0 U m . The control law u(t) ∈ U m that solves (14) will be denoted by u * m (t) while the control law u(t) ∈ U that solves (14) will be denoted by u * (t). To find u * (t) we will first find all the u * m (t) and then compare them and select u * (t) as the u * m (t) that minimizes over m the optimization target in (14).

It is not hard to see that the optimum control law u * m (t) must make the vehicle enter and exit the convex hull of each region {A k } m k=0 at most once. We will refer to these input and output points to the convex hull of the area A m as i m and o m respectively.

We also define the following set of points

P m = {s, i 1 , i 2 , • • • , i m , o m , o m-1 , • • • , o 1 , g}.
In addition we will index its elements as follows:

p m 0 = s, p m k = i k , for k = 1, 2, • • • , m, p m k = o 2m+1-k , for k = m + 1, m + 2, • • • , 2m, p m 2m+1 = g. ( 15
)
The time instant when the vehicle is located at p m k is t k and τ k = t k -t k-1 . Let us also write these points in polar coordinates as

p m k = r m k [cos(φ m k ) sin(φ m k )] T .
From the definition of i m and o m we know that they lie in a circle of radius r m which is known. Therefore we know {r m k } 2m k=1 and as a consequence the only unknowns to uniquely determine P j are the angles 1 {φ m k } 2m k=1 , where the φ m k is the angle of the point respect to the BS. Now, the optimum control law u * m (t) takes the vehicle from p m 0 up to p m 2m+1 in ascending order through each point 1 Since φ 0 and φ 2j+1 are the angles of s and g they are also known.

in P m . We also can see that the second term in the optimization target of (14) depends only the durations τ k and not on the shape of the particular path followed by the vehicle nor by its velocity profile. So velocity profile and the path must be selected to minimize the second term in (14). To do so the vehicle must go from p m k-1 to p m k in a time τ k (to be determined) using minimum energy. Using calculus of variations we can show that this is achieved by:

u m (t) = p m k -p m k-1 τ k ∀ t ∈ [t k-1 , k). ( 16
)
Therefore if we optimize the durations {τ k } 2m+1 k=1 and the angles {φ m k } 2m k=1 we obtain the candidate for the optimum control law u * m (t). To achieve this we need to solve ( 14) with the additional constraint u(t) ∈ U m . By adding this constraint and taking into account (16) the optimization problem (14) becomes:

min {τ k } 2m+1 k=1 ,{φ m k } 2j k=1 2m k=1 p m k -p m k-1 2 v 2 max τ k + T 0 b(t) N dt s.t. 2m+1 k=1 τ k = T, τ m+1 R m + m k=0 (τ k+1 + τ 2m+1-k )R k ≥ N, p m k -p m k-1 τ k ≤ v max (17)
where the first constraint ensures that the sum of all durations is equal to the duration T of the whole task. The second constraint ensures that all the N bits are transmitted. Note that the inequality allows for the possibility of transmitting all the N bits before the vehicle exits the coverage area. The next constraint ensures that the vehicle does not surpass the maximum velocity. It is not difficult to see that b(t), for n ≤ m, is given by:

b(t) = [N -S(n -1) -(t -t n )R n ] + ∀ t ∈ [t n , t n+1 ] (18
) and for m < n ≤ 2m + 1 we have:

b(t) = [N -S(m) -Z(n) -(t -t n )R n ] + ∀ t ∈ [t n , t n+1 ], (19)
where

S(n) = n k=0 τ k+1 R k and Z(n) = n-1 k=m+1 τ k+1 R 2m-k .
Proposition 2. Given a set of durations {τ k } 2m+1 k=1 the optimum set of angles {φ m k } 2m k=1 is a solution to the following set of nonlinear equations: To optimize {φ m k } 2m k=1 given {τ k } 2m+1 k=1 we find numerically all the solutions to (20) and then compare them in the optimization target (17). For lack of space we will not give more details on this optimization process but it will be detailed in future journal extension. Very interestingly, one can note the similarity of the equations (20) with the law of refraction in optics, showing the potential of generalizing the considered problem to a large variety of scenarios such as the case where the data rate varies continuously with the distance and the multiple base stations case.

τ n+1 r n-1 sin(φ n -φ n-1 ) = τ n r n+1 sin(φ n+1 -φ n ), for n = 1, 2, • • • , 2m. (20)
Therefore given {τ k } 2m+1 k=1 we can calculate the angles {φ m k } 2m k=1 as a consequence the optimization problem (17) now depends only on {τ k } 2j+1 k=1 . We solve (17) by using a modified version of the simulated annealing algorithm that takes into account all the restrictions on the variables {τ k } 2m+1 k=1 and calculates the angles {φ m k } 2m k=1 as function of the durations {τ k } 2m+1 k=1 .

SIMULATIONS

Now, to illustrate the trajectory optimization and gain more insight into the technique described in this paper we will now present some simulation results. For the simulations we will consider two different bit rates R 1 = 1 Mbps and R 2 = 2 Mbps. In addition the radii of the areas A 1 and A 2 will be r 1 = 1km and r 2 = 500m. Finally, the maximum velocity of the vehicle will be v max = 10m/s. We remind the reader that the task consist in transmitting the whole content of the buffer and then reaching the goal point g in a time T . The starting and goal points are arbitrarily located so that if the vehicle goes from s to g in straight line it never crosses a coverage area and therefore the task cannot be completed.

We will consider three different scenarios: (I) N = 700 Mb and T = 2000 s, this represents a case in which the amount of data to be transmitted is reasonable for the time to complete the task; (II) N = 150 Mb and T = 2000 s, this Fig. 2. The starting (resp. target) point of the vehicle is marked by a red (resp. black) point and is on the left (resp. right). The disk delineated by the pink circles correspond to geographical zones where the data rates are different. Each figure depicts, for u * 1 (t) (left) and u * 2 (t) (right), the vehicle trajectory which aims at uploading its data to the base station and minimizing the cost of motion to reach the target. The similarity with the refraction law in optics is well illustrated. represents a case in which the amount of data to be transmitted is considerably low respect to the time to complete the full task; and (III) N = 700 Mb and T = 1400 s, this represents the case in which the time to complete the full task is tight respect to the amount of data to be transmitted.

In Figs. 234we observe the path obtained by each of the two optimum control law candidates (i.e., u * 1 (t) and u * 2 (t)) for the three scenarios. Then in Fig. 5 and Fig. 6 we observe the temporal evolution of the buffer b(t) with the respective velocity profile of the optimum control law candidates. Then in tables 1 and 2 we observe, for both u * 1 (t) and u * 2 (t) , the value of the optimization target in (17) and the amount of mechanical energy used for the three cases considered. First of all we have to note that the shape of the path for u * m (t) depends on both the number N of bits to be transmitted and on the amount of time T to complete the task. This is because these parameters will determine the distribution of the durations {τ k } 2m+1 k=1 and the coefficients of nonlinear equations (20) which determine the path depend on such distribution.

Mechanical energy (J) Optimization target

In case I the optimum control law u * (t) is u * 1 (t), as can be observed from comparing tables 1 and 2. This is because the number N of bits to transmit is relatively small so reaching A 2 is unnecessary and only deviates more the vehicle from g therefore costing more energy than simply staying in A 1 .

Then in case II the amount of bits to be transmitted in the same amount of time is significantly larger. In this case the benefits of reaching A 2 to use the rate R 2 for the transmission are considerably and so the optimum control law is u * 2 (t). In case III the optimum control law also is u * 2 (t) but we can observe that, as opposed to case II, the length of the path segment in A 2 is longer. This is because in this case the vehicle is short of time but it needs to stay some time in A 2 to finish the transmission. So in order to save time and energy while it is in A 2 it tries to get as near as possible to the closest point of A 2 to g which lies in the line from the BS to g.

It is also worth noticing in Figs. [START_REF] Ooi | 6. The figure shows how the buffer size and the vehicle velocity varies over time for the optimal control laws u * 2 (t)[END_REF] and Figs. 6 the velocity profile. In the cases presented, and various other with are not shown due to lack of space, we observed some interesting patterns in the shape of the velocity profile. First of all, the velocity used between the starting point s and i 1 is always v max but the optimum location of i 1 is not, in general, the one that minimizes the between s and A 2 . The use of the maximum velocity to reach the coverage area comes from wanting to initiate the transmission quickly but the fact that i 1 does not minimize its distance with s and is slightly biased towards the goal point g comes from wanting to reach g without wasting much energy. This is a fair tradeoff between the two optimality criteria.

Another interesting aspect in the velocity profile is the fact that for u * m (t) the velocity used in the region A m is always the lowest. This is because the vehicle wants to take advantage of the highest rate in the trajectory R m to transmit as much bits as possible in that time. But it does not stop to keep reducing its distance with the goal point g.

CONCLUSIONS

We have proposed a series of techniques to design optimum trajectories that allow a vehicle (or an autonomous vehicle) to efficiently reach a target while simultaneously transmit data as quick as possible when operating under temporal constraints in a workspace sparsely covered by a network. This type of trajectories could be of interest when time critical data needs to be transmitted by vehicles which are energy constrained. In addition, sparsely covered networks assisted by vehicles are of interest since they could replace costly infrastructure by cheap autonomous vehicles. So this type of interdisciplinary algorithms that consider both the communi-Fig. 5. The figure shows how the buffer size and the vehicle velocity varies over time for the optimal control laws u * 1 (t) under different conditions. cations aspects and the energetic aspects of the vehicles will promote the development of such networks. The problem presented in this paper is still at an early stage of development and further promising research needs to be done. In particular: the generalization of the derived counterpart of the refraction law to the case of continuous variations of the data rate might be addressed; the case of multiple base-stations and multiple vehicles is also of high interest e.g., for 5G communications-driven car networks; the effect of shadowing in the wireless channel needs to be taken into account as well as considering the effect of the vehicle's inertia.

Fig. 1 .

 1 Fig. 1. Illustration of the angle differences considered in (20).

Fig. 3 .

 3 Fig. 3. The starting (resp. target) point of the vehicle is marked by a red (resp. black) point and is on the left (resp. right). The disk delineated by the pink circles correspond to geographical zones where the data rates are different. Each figure depicts, for u * 1 (t) (left) and u * 2 (t) (right), the vehicle trajectory which aims at uploading its data to the base station and minimizing the cost of motion to reach the target. The similarity with the refraction law in optics is well illustrated.

Fig. 4 .

 4 Fig. 4. The starting (resp. target) point of the vehicle is marked by a red (resp. black) point and is on the left (resp. right). The disk delineated by the pink circles correspond to geographical zones where the data rates are different. Each figure depicts, for u * 1 (t) (left) and u * 2 (t) (right), the vehicle trajectory which aims at uploading its data to the base station and minimizing the cost of motion to reach the target. The similarity with the refraction law in optics is well illustrated.

Table 1 .

 1 Results for u * 1 (t) under different scenarios.

	Case I	35721	635
	Case II	35607	909
	Case III	44328	1002
		Mechanical energy (J) Optimization target
	Case I	37165	651
	Case II	40206	847
	Case III	48569	926

Table 2 .

 2 Results for u * 2 (t) under different scenarios.

The author acknowledges the funding of CONACYT (México).