Fabien Campillo
email: fabien.campillo@inria.fr

Rivo Rakotozafy
email: rrakotozafy@uni-fianar.mg

Vivien Rossi
email: vivien.rossi@cirad.fr

Bayesian numerical inference for markovian models. Application to tropical forest dynamics

Keywords: Markov chain Monte Carlo method, interacting chains, hidden Markov model, forest dynamic model

Bayesian modelling is fluently employed to assess natural ressources. It is associated with Monte Carlo Markov Chains (MCMC) to get an approximation of the distribution law of interest. Hence in such situations it is important to be able to propose N independent realizations of this distribution law. We propose a strategy for making N parallel Monte Carlo Markov Chains interact in order to get an approximation of an independent N-sample of a given target law. In this method each individual chain proposes candidates for all other chains. We prove that the set of interacting chains is itself a MCMC method for the product of N target measures. Compared to independent parallel chains this method is more time consuming, but we show through example that it possesses many advantages. This approach will be applied to a forest dynamic model.

Introduction

In many models arising in environment (ecology, population dynamics, renewable resource management etc.), measurements y 1 , . . . , y T are collected yearly or monthly so that the real-time constraint is not relevant even if the underlying law features a temporal structure. State-space modeling of these data consists in proposing a Markov process (x t , y t ) t=1•••T , where the state process x t is not observed and y t are the associated observation process. This process usually depends on some unknown parameter θ with given a priori law. The goal of the Bayesian inference is to determine the a posteriori law of (x 1:T , θ) given the measurements y 1:T . In the case of general state-space modeling, we need to utilize approximation procedures. The success of the Bayesian inference is mainly due to the development of efficient Monte Carlo approximation techniques [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. Among them, MCMC methods allow us to sample from almost any prescribed distribution law [START_REF] Gilks | Markov Chain Monte Carlo in practice[END_REF]. Still high dimensional or "tricky" distribution laws are barely tackled by these techniques and should be approached with realistic expectations. Together with numerical Bayesian inference, hidden Markov models for general state-space have been recently used in environment sciences and ecology, see e.g. [START_REF] Rivot | A Bayesian state-space modelling framework for fitting a salmon stage-structured population dynamic model to multiple time series of field data[END_REF] MCMC algorithms [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] allow us to draw samples from a probability distribution π(x) dx known up to a multiplicative constant. This consists in sequentially simulating a single Markov chain whose limit distribution is π(x) dx. There exist many techniques to speed up the convergence toward the target distribution by improving the mixing properties of the chain.

In practice one however can make use of several chains in parallel. It is then tempting to exchange information between these chains to improve mixing properties of the MCMC samplers. A general framework of "Population Monte Carlo" has been proposed in this context [START_REF] Cappé | Population Monte Carlo[END_REF]. In this paper we propose an interacting method between parallel chains which provides an independent sample from the target distribution. Contrary to papers previously cited, the proposal law in our work is given and does not adapt itself to the previous simulations. Hence, the problem of the choice of the proposal law still remains.

Parallel/interacting Metropolis within Gibbs (MwG) algorithm

Let π(x) be the probability density function of a target distribution defined on (IR n , B(IR n )). For ℓ = 1, . . . , n, we define the conditional laws:

π ℓ (x ℓ |x ¬ℓ ) def = π(x 1:n )/ π(x 1:n ) dx ¬ℓ . (1)
where ¬ℓ def = {m = 1 : n; m = ℓ}. When we know to sample from (1), we are able to use the Gibbs sampler. It is possible to adapt our interacting method to parallel Gibbs sampler. But very often we do not know how to sample from (1) and therefore we consider proposal conditional densities π prop ℓ (x ℓ ) defined for all ℓ. In this case, we use MwG algorithm. One iteration X → Z of the parallel/interacting MwG method consists in updating the components X ℓ successively for ℓ = 1, . . . , n, i.e. [X

1:n ] → [Z 1 X 2:n ] → [Z 1:2 X 3:n ] • • • [Z 1:n-1 X n ] → [Z 1:n ].
For each ℓ fixed, the subcomponents X i ℓ are updated sequentially for i = 1 : N in two steps:

1. Proposal step: We sample independently N candidates Y j ℓ ∈ IR according to:

Y j ℓ ∼ π ℓ,prop i,j (ξ| Z, X i ℓ , X i ℓ ) dξ , 1 ≤ j ≤ n where Z, ξ, X i ℓ def =          Z 1:ℓ-1 Z 1 ℓ . . . Z i-1 ℓ ξ X i+1 ℓ . . . X N ℓ X ℓ+1:n          .
We also use the following lighter notation: π ℓ,prop i,j (ξ|ξ ′ ) = π ℓ,prop i,j (ξ| Z, ξ ′ , X i ℓ ).

Selection step:

The subcomponent X i ℓ could be replaced by one of the N candidates Y 1:N ℓ or stay unchanged according to a multinomial sampling, the resulting value is called Z i ℓ ,

i.e.:

Z i ℓ ←            Y 1 ℓ with probability 1 N α i,1 ℓ (X i ℓ , Y 1 ℓ ) , . . . Y N ℓ with probability 1 N α i,N ℓ (X i ℓ , Y N ℓ ) , X i ℓ with probability ρi ℓ (X i ℓ , Y 1:N ℓ )
where:

α i,j ℓ (ξ, ξ ′ ) def = π ℓ (ξ ′ |X i ¬ℓ ) π ℓ (ξ|X i ¬ℓ ) π ℓ,prop i,j (ξ|ξ ′ ) π ℓ,prop i,j (ξ ′ |ξ) ∧ 1 , ρi ℓ (X i ℓ , Y 1:N ℓ ) def = 1 -1 N N j=1 α i,j ℓ (X i ℓ , Y j ℓ ) .
Proposition Let P (X, dZ) the Markov kernel on IR n×N associated with the MwG algorithm, the measure

Π(dX) = π(X 1 ) dX 1 • • • π(X N ) dX N is invariant for the kernel P , i.e. ΠP = Π.
The proof is detailed in [START_REF] Campillo | Parallel and interacting Markov chains Monte Carlo method[END_REF].

Numerical tests

Our aim is to show that making parallel samplers interact could speed up the convergence toward the stationary distribution.

An hidden Markov model. We apply the parallel/interacting MwG sampler to a toy problem where a good estimate π of the target distribution π is available. Consider

s ℓ+1 = a s ℓ + w ℓ , y ℓ = b s ℓ + v ℓ for ℓ = 1 • • • n, where s 1 ∼ N (s 1 , Q 1 )
, w 1:n and v 1:n are centered white Gaussian noises with variances σ 2 w and σ 2 v . Suppose that b is known and a = θ is unknown with a priori law N (µ θ , σ 2 θ ). We also suppose that w 1:n , v 1:n , s 1 and θ are mutually independent. The state variable is x 1:n+1 def = (s 1:n , θ) and the target law is π(s 1:n , ϑ) ds 1:n dϑ def = law(s 1:n , θ|y 1:n = y 1:n ). We will perform two algorithms: (i) N parallel/interacting MwG samplers and (ii) N parallel/independent MwG samplers. For methods (i) and (ii) we perform N = 50 parallel samplers so we compute ǫ (i ) and ǫ (ii ) . ǫ (•) is an indicator which must decrease and remain close to a small value when convergence toward the stationary distribution occurs. (ii) for the parallel/independent MwG sampler (--), and ǫ (i) for the parallel/interacting MH sampler (-). Right: Evolution of the indicator ǫ (ii) for the parallel/independent MwG sampler (--). After 5000 sec. CPU time, the convergence of this method is still unsatisfactory. ǫ (•) is based on the L 1 errors between π and the kernel density estimates of the target density based on the final values of method (i) or (ii) (see [START_REF] Campillo | Parallel and interacting Markov chains Monte Carlo method[END_REF] for details).

To compare fairly the parallel/independent MwG algorithm and the parallel/interacted MwG algorithm, we represent on Figure 1 the indicators ǫ (i) and ǫ (ii ) not as a function of the iteration number of the algorithm but as a function of the CPU time. In Figure 1 (left) we see that even if one iteration of algorithm (i) needs more CPU than one of (ii), still the first algorithm converges more rapidly than the second one. This shows the inefficiency of parallel/independent MwG on this simple model.

Real case study: a tropical forest. Forest dynamics models have emerged to assess the sustainability of the management of tropical forest. Among various types of models, matrix models rely on a description of the forest stand by a discrete diameter distribution. The state of a tree corresponds to its belonging to a diameter class and the shifts from one class to another occur by discrete time step. Matrix model of forest dynamics rely on four hypotheses (see [START_REF] Picard | Modelling diameter class distribution with a secondorder matrix model[END_REF] for details) of which Markov's hypothesis. So the parallel/interacted MwG algorithm can be used to calibrate a matrix model. An illustration on a data set obtained from experimental forest plots in a tropical rain forest of French Guiana will be presented.

Conclusion

This work showed that making parallel MCMC chains interact could improve their convergence properties. We presented the basic properties of the MCMC method, we did not prove that the proposed strategy speeds up the convergence. This difficult point is related to the problem of the rate of the convergence of the MCMC algorithms. Through simple examples we saw that the MwG strategy could be a poor strategy. In this situation our strategy improved the convergence properties.

Figure 1 :

 1 Figure 1: Left: Evolution of the indicators ǫ(ii) for the parallel/independent MwG sampler (--), and ǫ (i) for the parallel/interacting MH sampler (-). Right: Evolution of the indicator ǫ (ii) for the parallel/independent MwG sampler (--). After 5000 sec. CPU time, the convergence of this method is still unsatisfactory.


Perception of Surface Defects by Active Exploration with a Biomimetic Tactile Sensor Raphaël Candelier, Georges Debrégeas and Alexis Prevost

Abstract-We investigate the transduction of tactile information during active exploration of finely textured surfaces using a novel tactile sensor mimicking the human fingertip. The sensor has been designed by integrating a linear array of 10 micro-force sensors in an elastomer layer. We measure the sensors' response to the passage of elementary topographical features in the form of a small hole on a flat substrate. The response is found to strongly depend on the relative location of the sensor with respect to the substrate/skin contact zone. This result can be quantitatively interpreted within the scope of a linear model of mechanical transduction, taking into account both the intrinsic response of individual sensors and the context-dependent interfacial stress field within the contact zone. Consequences on robotics of touch are briefly discussed.

I. INTRODUCTION

T HE tactile sensitivity of hands and fingertips allows humans to dexterously manipulate objects and obtain a wealth of information such as their shape, weight or texture [START_REF] Hollins | Somesthetic senses[END_REF]. Specialized nerve endings (mechanoreceptors) located in the first layers of the skin trigger sequences of spikes traveling up the afferent fibers, that are further processed by the central nervous system to ultimately form a representation of the probed object. Biomimetic tactile sensing aims at artificially reproduce this mechanism by mimicking the fingerpad. Biomimetic sensors generally consist of a collection of force sensors covered with an elastic membrane [START_REF] Fearing | Tactile sensing mechanisms[END_REF]- [START_REF] Dahiya | Tactile sensingfrom humans to humanoids[END_REF]. Different technologies have been used to produce the artificial counterparts of the mechanoreceptors, such as MEMS [START_REF] Scheibert | The Role of Fingerprints in the Coding of Tactile Information Probed with a Biomimetic Sensor[END_REF]- [START_REF] Muhammad | Development of a biomimetic mems based capacitive tactile sensor[END_REF], strain-sensitive materials [START_REF] Howe | Dynamic tactile sensing: perception of fine surface features with stress rate sensing[END_REF]- [START_REF] Lacour | An elastically stretchable tft circuit[END_REF] or thin-film devices [START_REF] Maheshwari | High-Resolution Thin-Film Device to Sense Texture by Touch[END_REF]. More recently, the role of the skin itself in the transduction process has been studied. It has been shown in particular that the topography of the skin may significantly modify the sensors response [START_REF] Scheibert | The Role of Fingerprints in the Coding of Tactile Information Probed with a Biomimetic Sensor[END_REF], [START_REF] Rossi | Polymer based interfaces as bioinspired 'smart skins[END_REF].

Despite the rapid development of this field in recent years, replication of the full manipulative capabilities of the human hand is still years away. Beyond the technological challenge of designing realistic tactile devices, efforts are still needed to understand how one can efficiently and robustly extract physical information about the probed object from the sensors signal. One of the fundamental difficulties owes to the fact that tactile perception is an active process: the sensor's response not only depends on the intrinsic property of the sensor system but also on the exploratory conditions, i.e. the way the object is probed. This may explain why, in natural tactile perception, various types of exploratory procedures are spontaneously used depending on the type of information one seeks to extract (shape, texture, weight, etc.) [START_REF] Lederman | Tactile Roughness of Grooved Surfaces: The Touching Process and Effects of Macro and Microsurface Structure[END_REF]. How these various exploratory conditions shape the elicited mechanoreceptors response and thus effectively participate to the transduction and encoding of tactile information is to date mostly unknown.

In this Letter we focus on the role of one exploratory parameter, namely the location of the fingertip/surface contact zone with respect to the sub-cutaneous receptors. To that end, we investigate the response of a multi-sensor biomimetic device to the passage of elementary topological defects and examine how this response varies with the location of each sensor below the contact zone.

II. BIOMIMETIC FINGERTIP DESIGN AND CALIBRATION

A novel biomimetic tactile sensor has been designed with the aim to reproduce the functioning principle and main mechanical characteristics of the human fingertip. The force sensing device is a linear array of 10 MEMS force sensors 1 . Successive sensors are 1mm apart and aligned along the xaxis (see fig. 1-top), which is the scanning direction in dynamic sensing experiments. Sensors are numbered from left to right and a consistent color code is used in all figures, from blue (1) to red [START_REF] Rossi | Polymer based interfaces as bioinspired 'smart skins[END_REF]. The sensitive part of each force sensor consists of a vertical silicon cylinder of diameter 200µm and height 400µm attached to an horizontal suspended silicon membrane of radius 350µm (see inset in fig. 1-top). This membrane bears 4 pairs of piezoresistive gauges that measure its internal stress state, from which one can extract the three components of the force acting on the cylinder.

The array of force sensors is embedded in an elastic spherical cap mimicking the skin, with a radius of curvature R = 129.7mm and a maximum thickness h = 3.04mm. The apex of the elastic layer is positioned roughly above the middle of the array. This elastic cover is made of cross-linked Poly-DiMethylSiloxane (PDMS, Sylgard 184, Dow Corning) whose Young's modulus and Poisson ratio are respectively ∼ 3M P a and 0.5. Its surface is rendered rough at a micrometric scale by mechanically abrading the spherical mold used for its design. This fine texturing significantly reduces adhesive and frictional forces which may otherwise result in the deterioration of the sensor.

The biomimetic sensor is integrated in a mechanical setup which allows for controlled tactile stimulation. The probed substrate is attached to a three-axis motorized platform through a double cantilever system coupled with capacitive position sensors that enable one to simultaneous record the global normal force F z and tangential force F x .

Individual micro-force sensors are first calibrated using a standard quasi-punctual indentation protocol [START_REF] Howe | Dynamic tactile sensing: perception of fine surface features with stress rate sensing[END_REF], [START_REF] Scheibert | Experimental evidence of non-amontons behaviour at a multi-contact interface[END_REF] sketched in fig. 1-top. A 500µm in diameter metallic wire is indented normally to the skin at various positions under a controlled load F z in the range 0 -0.4N . Regardless of the indentor's location, the sensors' response varies linearly with the applied load. Figure 1-bottom displays this response for the ten sensors corresponding to a unit load (F z = 1N ) as a function of the position x of the indentor, for both the normal σ z and tangential σ x components.

These point-load response profiles are compared with the result of a finite element simulation using CAST3M2 : the stress field at the bottom of an elastic layer of finite thickness h rigidly attached to a solid substrate is calculated as its surface is normally indented by a 500µm flat punch. In the limit where the punch diameter is much smaller than the layer thickness, the calculated stress field is similar to that obtained for a punctual indentor (the so-called Green function) and is a function of h only. As shown in the bottom panels of figure 1, the measured sensors' response profiles and calculated stress fields can be correctly adjusted for both directions of the stress. This adjustment provides a Volt-to-Pascal calibration of the sensors, as well as an estimated value of the thickness of the layer at each sensor's location. The latter is consistent with the known geometry of the elastic cap: in particular, the layer thickness is found to vary continuously along the line of sensors as expected for a spherical cap, with a maximal thickness of the skin located close to sensor ♯6. This thickness variation can be clearly seen on fig. 1 through an inverse variation of the maxima of σ z (corresponding to the point where the indentor lays right above the sensor).

III. STRESS FIELD DURING EXPLORATION OF A SMOOTH

SUBSTRATE

We first investigate the response of the biomimetic sensor to a flat and smooth substrate, under both static and dynamic exploration. A Plexiglas plate is put into contact with the sensor under a constant load (F z = 0.8N ) forming a circular contact of diameter 4 -5mm. The interfacial contact stress field is in this case centro-symmetric with respect to the center of the contact zone. The corresponding normal and tangential stress σ z and σ x , measured with the micro-force sensors, display respectively symmetric and anti-symmetric profiles around the contact center, as shown in figure 2-middle. A slight asymmetry in the x-direction is however observed which can be accounted for by the development of a minute tangential force under normal loading associated with the presence of the cantilever [START_REF] Scheibert | Experimental evidence of non-amontons behaviour at a multi-contact interface[END_REF]. The substrate is then put into sliding motion along the x-axis at constant velocity v = 500µm.s -1 . The time-averaged profiles in steady sliding are shown in the bottom panel of figure 2. These profiles strongly differ Fig. 3. Effect of the finite thickness of the elastic layer. Left Sphere-on-plane contact pressure profile calculated for a semi-infinite elastic medium (dashed black line) and an elastic layer of thickness h (blue line) such that h/a = 1.3 where a is the contact radius. This ratio corresponds to the experimental conditions. Right Normal component gzz of the Green tensor in (x, y = 0, z = h): the profile for the semi-infinite layer (dashed line) can be rescaled (blue line) to approximate the profile obtained by finite element simulation for a layer of thickness h = 3.04mm (dotted line).

from their static counterpart as a result of the development of friction-induced tangential stress in the contact zone.

In order to interpret the sensors' response in both static and dynamic conditions, a simple mechanical model of linear elastic transduction is developed. Three ingredients are needed: the pressure field at the skin/substrate interface, a local law of friction to derive the tangential stress in the contact zone and a description of the stress propagation in an elastic layer of finite thickness.

• (i) Pressure field at the interface. The interfacial pressure field in a sphere-on-plane contact under normal load is given by Hertz calculation [START_REF] Johnson | Contact Mechanics[END_REF]. This profile however needs to be corrected to take into account the finite thickness of the elastic layer. Fretigny and Chateauminois have recently proposed a numerical solution to this problem [START_REF] Fretigny | Solution for the elastic field in a layered medium under axisymmetric contact loading[END_REF]. A comparison between the Hertz pressure profiles for an infinite and a finite thickness layer corresponding to our system geometry is shown in figure 3left. The modified pressure profile has a smaller contact radius and the pressure at the center is higher.

• (ii) Local law of friction. In steady sliding, we postulate a local Coulomb law everywhere in the contact zone, i.e. we assume a linear relationship between the tangential stress p x and the normal stress p z at the interface:

p x = µ d .p z , (1)
where the dynamical friction coefficient µ d ∼ 3 is the ratio of the global tangential and normal forces F x /F z . We also assume that this friction-induced tangential stress field does not affect the pressure field discussed above.

• (iii) Stress propagation through the elastic layer. Under linear elasticity and quasi-static hypothesis, the stress propagation through the elastic layer is entirely given by the Green tensor {g ij (x, y)} which characterizes the stress field at position (x, y, z = h) in response to a unit localized force applied at the surface of the layer at position (x = 0, y = 0). For a semi-infinite elastic medium, an analytical solution to this tensor exists [START_REF] Johnson | Contact Mechanics[END_REF]. Again, in the present configuration, this result needs to be corrected to take into account the finite thickness of the elastic layer and its attachment to a rigid base. The finite element analysis we have carried out for normal indentation indicates that a reasonable approximation for this correction can be obtained by rescaling the semiinfinite solution by a simple numerical factor k, while the spatial coordinate is rescaled by a factor k -1/2 to ensure force conservation. This is shown in fig. 3-right which compares the result of the finite-element calculation with the rescaled semi-infinite solution with a scale factor k = 1.48. The same factor yields a correct approximation to the tangential component of the stress (not shown). We make the assumption that the same correction holds for the other components of the Green tensor, i.e. for the response to a tangential localized force. With these three ingredients, an estimated stress profile along the sensor array axis (y i = 0) at depth z i = h can be derived by convoluting the interfacial stress field with the Green tensor:

σ z (x i ) = p x (x, y)g xz (x -x i , y, h)dxdy + p z (x, y)g zz (x -x i , y, h)dxdy, σ x (x i ) = p x (x, y)g xx (x -x i , y, h)dxdy + p z (x, y)g zx (x -x i , y, h)dxdy,
In spite of the numerous assumptions made in this model, fig. 2 demonstrates that the predicted stress profiles (black curves) consistently captures the experimental data (colored circles) in both static and dynamic conditions. Note that, for the static case, the observed residual tangential force at rest has been accounted for by postulating a tangential stress field p x proportional to the normal stress field p z .

IV. RESPONSE TO ELEMENTARY TOPOLOGICAL DEFECTS IN DYNAMIC EXPLORATION

We now turn to the response of the biomimetic tactile sensor to the presence of an elementary topological defect on the flat surface. The same Plexiglas plate is patterned with individual holes of diameter 500µm and depth 500µm sufficiently well separated in space so that only one hole at most can lay within the contact zone. In static contact, the presence of a hole within the contact zone cannot be detected by the sensor. However, in dynamic exploration, i.e. when the hole is scanned across the sensor along the MEMS array axis, the passage of the hole in the contact zone can be easily identified by subtracting the time signal of each sensor with its time-averaged value. We denote ς z (u) and ς x (u) the stress variation signals thus defined, with u being the position of the hole along the x direction of motion (see fig. 4-top). Fig. 4-left shows the profiles obtained by averaging over 14 passages of similar holes, for the ten sensors whose locations range from -4.5 to 4.5mm. Standard deviations are typically 0.6kP a for ς z and 0.4kP a for ς x . These dynamic response profiles, unlike the localized force responses shown in fig. 1, display very dissimilar shapes depending on the sensor's precise location with respect to the contact zone. The normal stress profiles ς z (u) evolve from a positive hump on the left-hand side sensors to a negative one on the right-hand side sensors via an intermediate antisymmetric profile (see sensors ♯6 and 7, in green). In contrast, the shear-stress profiles show a dominant hump on the left and right-hand side sensors, with an intermediate response exhibiting two comparable maxima. Note that the outmost sensors do respond to the defect's passage although they lay outside of the contact zone.

This result can be quantitatively interpreted based on the simple mechanical model presented before. The interfacial pressure field at the location of the hole is strictly zero since the elastic layer does not make contact with the substrate in this area. We further assume that the interfacial stress field in the rest of the contact is unchanged except for an overall uniform rescaling factor to ensure that the integral of the stress over the contact zone remains equal to the imposed overall load F z . Under this assumption, the stress modification induced by a hole located at (u,y = 0) can be approximated, for a sensor located in (x i , y i = 0, z i = h), as:

ς z (u) ≃ -[p x (u)g xz (u -x i ) + p z (u)g zz (u -x i )] S d ς x (u) ≃ -[p x (u)g xx (u -x i ) + p z (u)g zx (u -x i )] S d
where S d is the hole's surface area. These predictions are shown in the right panels of fig. 4 for ten 1mm-spaced sensors whose positions are depicted on the horizontal axis with colored circles, allowing for a direct comparison with the experimental measurements. Interestingly, the profiles predicted by the model are not only similar in shape to the measured ones but they are also quantitatively equivalent (compare for instance the profiles for sensor ♯5 in cyan, which can be almost superimposed). The major quantitative difference occurs for the sensors laying on the right-hand side of the contact (from ♯6 to 10): the measured stress variations amplitudes are up to 4 times larger than their predicted counterpart, both for the normal and tangential components. This discrepancy may reflect the need to refine our description of the interfacial stress field, in particular to account for the coupling between tangential and normal stresses: the friction-induced tangential stress field produces a torque on the elastic layer which is balanced by a dissymmetrization of the interfacial pressure profile [START_REF] Scheibert | Stress field at a sliding frictional contact: Experiments and calculations[END_REF]. This effect, ignored in our simple model, is expected to increase the signal experienced by the sensors sitting at the front of the contact.

V. CONCLUSION AND PERSPECTIVES

This study constitutes the first step towards the comprehension of mechanical transduction of roughness information in realistic exploratory conditions. It allowed us to exhibit individual sensors' response to the scanning of elementary topological defects. These dynamic response functions should enable one to predict the linear response to any number of such holes patterned on a smooth substrate. By comparing such predictions to actual measurements, one may in turn characterize possible corrections that would reflect the nonlinear properties of solid frictional interfaces.

A central issue in artificial touch perception is the implementation of an inversion method, i.e. a signal processing procedure allowing one to consistently extract a representation of the probed environment from the measured sensors' responses. A naive strategy would consist in considering each sensor independently and extract from their signal, based on their intrinsic characteristics, the physical information of the substrate that lay within their sensitive area. Our results show that such an approach is likely to fail: in realistic conditions of exploration, the contact zone between the surface and the sensor cannot be known with great precision but the response of each sensor will critically depends on their position with respect to the contact zone. This difficulty is reminiscent of the problem of face recognition in image processing for instance, where variability in the lighting conditions needs to be accounted for in the signal analysis. This analogy suggests that in tactile exploration as well, the signal processing should include a pre-treatment stage based on multi-sensors signal analysis to account for contact stress heterogeneities at the skin/substrate interface.

Fig. 1 .

 1 Fig. 1. Calibration of the biomimetic sensors. Top Sketch of the experimental system along with the indentation protocol: a line of 10 MEMS sensors are embedded in a spherical elastomer cap. It allows for the measurement of the subcutaneous stress induced by a local indentation with a wire of diameter 500µm. Inset Image as seen from above of one of the 10 MEMS sensors. Bottom Normal (right) and tangential (left) stress σz and σx measured by the 10 micro-force sensors as a function of the indentor's position x (y = 0) for a normal force Fz = 1N . The stress measured by sensor ♯5 (circles) is compared to the prediction of a finite element calculation (thick black curve).

Fig. 2 .

 2 Fig. 2. Static and dynamic average stress as measured by the sensors compared to the model prediction. Top Sketch of the experiment: a smooth substrate is pushed against the skin surface with a confining normal load Fz = 0.8N . In dynamic exploration, the substrate moves with a velocity v along the sensor array axis. Middle Measured average stress in the static regime for the 10 sensors (circles) compared to the model's prediction (black curve) for the normal component σz (left) and the tangential component along the direction of motion σx (right), as functions of the sensor's location x. Position x = 0 is the center of the contact. Error bars show the measured standard deviation and colors correspond to the locations depicted on the sketch. Bottom Same as above but in the steady sliding regime (v = 500µm.s -1 ).

Fig. 4 .

 4 Fig. 4. Response to an elementary topological defect scanned across the skin's surface. Top Sketch of the protocol: a smooth substrate is put into contact with the tactile sensor under a normal load Fz = 0.8N . A single defect (500µm hole) is scanned along the direction of the sensor's array at constant velocity v = 500µm.s -1 . Middle Normal stress fluctuations ςz measured by the sensors (left) and predicted by the model (right) as functions of the defect's position u. The 10 curves correspond to different probing locations with respect to the contact zone, shown with colored dots on the right plot. Bottom Same as above, but for the tangential component ςx along the direction of motion.

The MEMS sensors have been developed by Patrice Rey and collaborators at CEA-LETI in Grenoble, France.

CAST3M is an open source finite element simulation software. See http://www-cast3m.cea.fr/cast3m/index.jsp

ACKNOWLEDGMENT

The authors would like to thank A. Chateauminois for the computation of the modified Hertz stress field and J. Scheibert and J. Frelat for their help on the finite element simulation.

Dr. Raphaël Candelier is a post-doc fellow in the sys-bio team of the LPS-ENS, France. He graduated the ENS/EHESS master of Cognitive Sciences, studied psychoacoustics at the IML (Doshisha U., Kyoto) and psycholinguistics at the LSCP (ENS, Paris). He was a member of the Starflag project, an European experimental work aimed at understanding the rules of collective behavior. Nevertheless his major is Physics, which he studied at the ESPCI and during his PhD in the GIT-SPEC (CEA, Saclay), where he specialized in statistical physics.

Dr. Georges Debrégeas is a CNRS experimental physicist. He received his PhD in 1997, from UPMC. He then joined the James Frank Institute in Chicago as a Post-Doc. In 1999, he got a position as a CNRS researcher, first in Strasbourg (ICS) then in College de France in Paris (LFO), and finally in ENS in 2004. In the earlier years of his career, he studied the rheological properties of various complex systems, such as granular materials and liquid foams. More recently, he got interested in the mechanics of contact and friction and its application to tactile perception in both humans and rodents.

Dr. Alexis Prevost is a CNRS experimental physicist whose areas of expertise range from wetting phenomena and granular materials to friction physics. He received his PhD degree in 1999 from the Univ. of Paris XI, a work completed at the LPS-ENS in Paris, France. In 2002, he joined the Physics Department of Georgetown Univ. as a postdoctoral fellow to study granular media dynamics. Since 2002, he has been holding a full-time CNRS researcher position, at LFO-College de France in Paris until 2005, and at LPS-ENS ever since. He now studies the contact mechanics of elastomers on rigid surfaces. He has been developing contact imaging techniques and direct stress measurements using micro-force sensors, in particular to gain a better understanding of the mechanical transduction of the tactile information involved in the human tactile perception.


C Force

S Grévy

L Gaudefroy

O Sorlin

L Caceres

F Rotaru

J Mrazek

N L Achouri

J C Angélique

F Azaiez

B Bastin

R Borcea

A Buta

J M Daugas

Z Dlouhy

Zs Dombrádi

F De Oliveira

F Negoita

Y Penionzhkevich

M G Saint-Laurent

D Sohler

M Stanoiu

I Stefan

C Stodel

F Nowacki

Prolate-Spherical Shape Coexistence at N=28 in 44 S

Keywords: numbers: 23.20.Lv, 25.70.Mn; 27.40.+z, 29.30.Kv

The structure of 44 S has been studied using delayed γ and electron spectroscopy at ganil. The decay rates of the 0 + 2 isomeric state to the 2 + 1 and 0 + 1 states have been measured for the first time, leading to a reduced transition probability B(E2 : 2 + 1 →0 + 2 )= 8.4(26) e 2 fm 4 and a monopole strength ρ 2 (E0 : 0 + 2 →0 + 1 ) = 8.7(7)×10 -3 . Comparisons to shell model calculations point towards prolatespherical shape coexistence and a phenomenological two level mixing model is used to extract a weak mixing between the two configurations.

'Magic' nuclei exhibit large gaps between the occupied and valence orbits. They are cornerstones of the nuclear structure, being used (i) to test our understanding of the nuclear forces which form these gaps and (ii) to model more complicated systems having many valence nucleons. While nuclei having 8 and 20 protons (or neutrons) can be reproduced by modeling the atomic nucleus with an harmonic oscillator potential, a spin-orbit interaction must be added to describe heavier magic nuclei. This spin-orbit interaction strongly binds nucleons having their angular momenta ℓ aligned with their intrinsic spin value s, denoted as ℓ ↑ . This leads throughout the chart of nuclei to regular sequence of orbits ℓ ↑ , (ℓ -2) ↑ , (ℓ -2) ↓ , ℓ ↓ , with the so-called large spin-orbit gaps 14, 28,50,82 and 126 between the lowered ℓ ↑ orbit (ℓ=2, 3, 4, 5 and 6) and the others. Generally, in particular at the stability, these gaps are large enough to prevent excitations between occupied and valence orbits and these magic nuclei are spherical. However, as the orbits forming the gap are separated by two units of angular momentum, quadrupole excitations are likely to develop if for some reason the shell gap is reduced. In this hypothesis, the development of quadrupole excitations jeopardizes the rigidity of the spherical gap and conduct the nucleus to deform. Consequently the doubly magic nuclei which have proton and neutron spin-orbit shell closures could become vulnerable to quadrupole excitations, as both protons and neutrons could act co- * Electronic address: grevy@in2p3.fr herently to deform the nucleus. So far the prototypical deformed nucleus composed of such a double spin-orbit shell-closure is 42 14 Si 28 [1]. At N=28 a gradual development of deformation occurs between the spherical doubly magic 48 20 Ca 28 and the deformed 42 14 Si 28 . In between these two extremes, i.e. in 44 16 S 28 , competition between spherical and deformed shapes is expected to be present leading to shape coexistence [2][3][4]. Depending on the strength of the quadrupole correlations induced by the cross shell excitations either the spherical normal configuration, or the deformed one, becomes the ground state while the other configuration forms a low lying 0 + 2 state. Therefore the discovery and characterization of this 0 + 2 state in 44 S 28 represent crucial information for understanding the evolution of N=28 shell gap. The non spherical nature of the 44 S ground state was suggested by its short β halflife and weak neutron-delayed emission probability [7], by the low energy of the 2 + 1 state (1297(18) keV), and the enhanced reduced transition probability B(E2 : 2 + 1 →0 + 1 ) of 63(18) e 2 fm 4 [8]. However the 2 + 1 and B(E2) values are intermediate between a rigid rotor and a spherical nucleus. It suggests a possible mixing of spherical and deformed shapes which can be deduced by studying the properties of the 0 + 2 isomer at 1365(1) keV observed in [9]. Already the study of a 7/2 -isomer in 43 S [5,6] has shed light on shape coexistence in the N≈28 region. Other cases of shape coexistence around shell closures have been reported in [10,11].

The present letter reports on the determination of the monopole strength ρ 2 (E0 : 0 + 2 →0 + 1 ) and the reduced transition probability B(E2 : 2 + 1 →0 + 2 ) in 44 S, extracted from the measurement of the half-life and the branching ratio between the E0 and E2 decay of the isomeric 0 + 2 state. These pieces of information were obtained by using combined γ and electron delayed-spectroscopy and are used to demonstrate the shape coexistence in 44 S.

The experiment was carried out at the Grand Accélérateur National d'Ions Lourds (ganil) facility. A primary beam of 48 Ca at 60A•MeV (I∼2eµA) impinged onto a 138 mg/cm 2 Be target to produce neutron-rich fragments. They were separated by the lise3 spectrometer [12] using an achromatic 100 mg/cm 2 Be degrader. The magnetic rigidity was set to optimize the transmission of the 44 S nuclei, produced at a rate of 200 sec -1 , with a momentum acceptance of ±1.45%. Fragments were identified on an event by event basis by means of their energy loss and magnetic rigidity (Bρ) values. The Bρ was obtained from the position of the fragments at the dispersive focal plane given by a multi-wire proportional chamber (caviar) [START_REF] Perrot | 11 th Int. Conf. on heavy ion accelerator technology[END_REF]. The selected nuclei were implanted in a 125 µm kapton foil tilted at 20 degrees with respect to the beam axis. Before the foil, a stack of Si detectors, including a position-sensitive one, was used to adjust the implantation depth and to reconstruct the position of the ions in a plane perpendicular to the beam axis. A thick Si detector located downstream to the implantation foil was used as veto. The γ and electron decay events were registered up to 20 µsec after the implantation. Electrons were detected in four cooled 45*45 mm 2 , 4 mm thick Si(Li) detectors, placed 20 mm above and below the beam axis. The γ-rays were measured by two clover Ge detectors of the exogam array located on the side of the implantation foil, at a distance of 25 mm to the beam axis. The use of a parallel beam optics along 2 meters length enables to derive the ion implantation profile on the kapton foil from the position-sensitive Si detector. This ion profile, the geometry of the detectors and that of the chamber were used as ingredients in a geant4 simulation to derive the electron (ǫ e -) and γ (ǫ γ ) efficiencies. The simulated efficiencies compare well with the ones obtained with calibrated sources of 207 Bi and 152 Eu placed in calibration runs at 6 different positions on the implantation foil. Using these comparisons, ǫ γ =3.06(5)% and ǫ e -=13.3(6)% were adopted for a gamma-ray of 1329 keV and an electron of 1362.5 keV, respectively [14].

The decay of the 0 + 2 to the ground-state (E0) proceeds through the emission of an internal conversion electron (IC) and by internal pair formation (IPF). The electron spectrum, following the implantation of a 44 S nucleus, is shown in Fig. 1. A single peak is observed at 1362.5 (10) keV corresponding to an excitation energy of 1365(1) keV for the 0 + 2 state, after having corrected for the binding energy of the K electrons in the 44 S nucleus. The integral of the peak is I e - IC (E0)= 148(8)×10 3 . The low energy part of the spectrum is well accounted for by the pair formation (IPF) mechanism in which electrons and positrons share an energy of 1365-(2*511)=343 keV. The fit of the electron time distribution (insert of Fig. 1) leads to an half-life of 2.619(26) µsec, which agrees with the value of 2.3(5) µsec reported in [9]. The 0 + 2 →2 + 1 decay branch, occurs through a strongly converted E2 transition at 36(1) keV, an energy below the experimental threshold of the detection system. The energy of this unobserved transition is derived from the measured energies of the 0 + 2 and 2 + 1 states, the latter being obtained from the observation of a delayed 2 + 1 →0 + 1 transition at 1329.0(5) keV (half-life of 2.66 (23) µsec, in agreement with the value obtained from the electron spectrum) which follows the 0 + 2 →2 + 1 decay. The 1329 keV energy agrees with the value of 1297 (18) found in Ref. [8]. The yield of the 0 + 2 →2 + 1 transition, I γ (E2), has been extracted from the number of delayed 2 + 1 →0 + 1 γ-rays. As can be seen in the insert of Fig. 2, this transition is contaminated by the 1332.5 keV γ-ray of 60 Co arising from the activation of the last selection slits of the spectrometer, which also produce a 1173 keV γ-ray with the same intensity. The number of counts in the 1329 keV peak has been obtained by fitting the γ spectrum with two gaussians, the intensity of the 1332.5 keV transition being constrained by that of the 1173 keV γ-ray.

The resulting I γ (E2) is 56(3)×10 3 . The decay of the 0 + 2 state occurs through E2 and E0 transitions, the ratio of which is expressed as :

R = λ(E2) λ(E0) = I γ (E2) I e - IC (E0) 1 + α conv (2 + 1 → 0 + 1 ) 1 + ΩIP F ΩIC (1)
In this expression, the electronic factors for pair formation and internal conversion have been extrapolated for a nucleus with A=44 from Ref. [START_REF] Passoja | [END_REF][16][17] to be Ω IP F =1.495×10 7 sec -1 and Ω IC =1.1125×10 7 sec -1 , respectively. A value of 3.6×10 -5 has been taken for the conversion coefficient α conv (2 + 1 →0 + 1 ) [18]. Using the experimental values of electron I e - IC (E0) and γ-ray I γ (E2) yields derived above, the resulting branching ratio is R=0.163 [START_REF] Perrot | 11 th Int. Conf. on heavy ion accelerator technology[END_REF] values are obtained using the following equations :

ρ 2 (E0) = ln(2) T 1/2 (0 + 2 )(1 + R)(Ω IP F + Ω IC ) (2)
B(E2) = 5.65 × 10 -10 5E 5 γ T 1/2 (1 + 1 R )(1 + α conv (0 + 2 → 2 + 1 ) (3)
Using the measured branching ratio R, the halflife value T 1/2 (0 + 2 ) and α conv (0 + 2 →2 + 1 )=10.94(1) extrapolated from Ref. [18], the monopole strength ρ 2 (E0 : 0 + 2 →0 + 1 ) and the reduced transition probability B(E2 : 2 + 1 →0 + 2 ) have been determined to be 8.7(7)×10 -3 and 8.4(26) e 2 fm 4 , respectively.

The values of E(0 + 2 )=1365(1) keV and ρ 2 (E0)=8.7(7)×10 -3 are the smallest measured in this mass region, pointing to a weak mixing between the 0 + 1 ground state and the 0 + 2 isomer and therefore to shape coexistence. In case of a large mixing, these states would repel each other to exhibit a large energy spacing and a larger ρ 2 (E0) value. To obtain further understanding on the nature of the shape coexistence, data are compared to shell model calculations.

Shell model (SM) calculations have been performed for 44 16 S 28 using the antoine code [19] and the up-to-date sdpf-u interaction that accounts remarkably well for nuclear structure in this mass region [20]. The full sd (f p) valence space has been considered for protons (neutrons) using standard effective charges e π = 1.35 e (e ν = 0.35 e).

The results gathered in Table I show a good agreement with the experimental values, the only exception is a somewhat larger calculated B(E2 : 2 + 1 →0 + 2 ) value than measured. Nevertheless, both experiment and calculation agree with the fact that the 2 + 1 state connects much strongly with the 0 + 1 state than with the 0 + 2 one. Indeed, the experimental B(E2 : 2 + 1 →0 + 1 )/B(E2 : 2 + 1 →0 + 2 ) ratio is 7.5 whereas the calculated one is 3.2. Calculated excited states connected to these two 0 + states are presented in Fig. 3 with their intrinsic quadrupole moments

Expt FIG. 3: 44 S level scheme calculated within the present SM approach (left), compared with available experimental data (right). E2 transition probabilities (in e 2 fm 4 ) are reported on top of black arrows and intrinsic quadrupole moments (in ef m 2 ) are shown in light gray on the right side of calculated levels. The ground state of the nucleus is head of a rotational band (β ≃ 0.25) and coexists with the rather spherical low-lying 0 + 2 isomer. Calculated values of the N=28 gap and correlation energies (in MeV) are given for even-even N=28 isotones.

Q 0 . For sake of clarity only the states of present interest are shown in this picture. Remarkable is the presence of 2 + 1 , 4 + 2 and 6 + 2 states on top of the 0 + 1 ground state connected by large B(E2) values. These states present equal Q 0 values of about 60 ef m 2 . These two features characterize the presence of a rotational band from an axially deformed nucleus with β ≃ 0.25. The 2 + 2 state at 2.14 MeV has a smaller intrinsic quadrupole moment Q 0 =-0.3 ef m 2 compatible with a spherical shape. A candidate for the 2 + 2 state is proposed at 2335(39) keV by placing the previously reported 988 keV transition [21] on top of the 0 + 2 or 2 + 1 state. Hence SM calculations suggest a prolate-spherical shape coexistence in 44 S.

A detailed analysis of the components contributing to the total energy of the 0 + states has been performed in order to deepen our understanding on the evolution of the collectivity from 48 20 Ca to 42 14 Si. Within the SM framework, the total Hamiltonian can be separated into its monopole (i.e. spherical mean-field contribution) and multipole (i.e. correlations mainly of pairing and quadrupole type) parts [22]. As can be seen from the values reported in Fig. 3, correlations strongly increase from the dou-bly magic 48 Ca (≃2 MeV) down to the exotic deformed 42 Si (≃18 MeV), while the size of the N=28 shell gap gets slightly reduced [23]. This increase of correlations is favored on one hand by neutron quadrupole excitations across the N=28 gap between the f 7/2 and p 3/2 orbits [23], and on the other hand, by the degeneracy of proton s 1/2 and d 3/2 orbits and excitations from the d 5/2 shell [1,[24][25][26]. In both cases, quadrupole correlations are favored by the fact that occupied and valence states are separated by two units of angular momentum. Without considering multipole contributions to the 0 + 1 and 0 + 2 states in 44 S, both levels are found to be quasi-degenerate in energy, and the ground state of 44 S is spherical. A gain of 1.5 MeV from the multipole energy brings the deformed configuration at the minimum of binding energy, while the spherical configuration corresponds to the excited state. Similar multipole effects energetically favor the oblate 0 + state in 42 Si which is predicted to coexist with a prolate 0 + state [20] at 1293 keV. The shell model calculation uses an Harmonic Oscillator basis for the description of the atomic nucleus. From the definition of the E0 operator, the calculated E0 transition between states of the same harmonic oscillator shells, as for protons in the sd shells and neutron in the f p shells, is strictly zero. Therefore, in order to shed light on the amount of mixing between the 0 + 1,2 states and to deduce their shape before mixing, we use a phenomenological two interacting levels model. We assume two spherical (S) and deformed (D) states before mixing which interact to produce 0 + 1 and 0 + 2 states defined as :

|0 1 = cosθ|0 D + sinθ|0 S (4) |0 2 = -sinθ|0 D + cosθ|0 S ( 5
)
where θ is the mixing angle. The E2 transition between the 2 + 1 and 0 + 2 (or 0 + 1 ) states being mainly due to the D-component of these 0 + states, it follows that B(E2 : 2 + 1 →0 + 2 )/B(E2 : 2 + 1 →0 + 1 )∼tan 2 (θ) (eq. 2 of [27]). A mixing amplitude tan 2 (θ)=0.13 is deduced from the experimental B(E2) values whereas the shell model gives a somewhat larger value of 0.24, both being smaller than the case of a maximum mixing (tan 2 (θ)=1). Therefore, the shape coexistence is found to be more pronounced experimentally than calculated by the SM. The magnitude of the monopole matrix element can be written as a function of the mixing amplitude and of the difference of shapes, β S and β D , between the two configurations before mixing [28], ρ 2 (E0)=(3Ze/4π) 2 sin 2 θcos 2 θ(β 2 Dβ 2 S ) 2 . Using the experimental mixing amplitude value of tan 2 (θ) in this equation, the experimental monopole strength is reproduced only when deformations β D ≃ 0.274 and β S =0 are assumed. The deformation parameter β D is in close agreement with the values obtained after mixing from Coulomb excitation experiment [8], β = 0.258(36), and from the shell model calculations, β=0. 25. Altogether these values again point towards a deformed-spherical shape coexistence in 44 S.

To summarize, electron and γ delayed-spectroscopies have been used to determine the monopole strength ρ 2 (E0 : 0 + 2 →0 + 1 )=8.7(7)×10 -3 and the reduced transition probability B(E2 : 2 + 1 →0 + 2 )=8.4(26) e 2 fm 4 in the 44 16 S 28 nucleus. Using these values, the earlier measured B(E2 : 2 + 1 →0 + 1 ) [8], shell model calculations and a two level mixing model, it is found that 44 S exhibit a shape coexistence between a prolate ground state (β ≃ 0.25) and a rather spherical 0 + 2 state. This establishes how the onset of collectivity progressively develops between the spherical 48 20 Ca and the deformed 42 14 Si nuclei. This study completes uniquely the understanding of the shellbreaking mechanism at the spin-orbit closed-shell N=28, which is as well of importance for the evolution of other shell gaps having the same origin.

FIG. 1 :

 1 FIG. 1: Electron energy spectrum obtained from the Si(Li) detectors following the implantation of 44 S nuclei. The peak at 1362.5(1.0) keV corresponds to the 0 + 2 →0 + 1 E0 transition. The low energy part is due to pair creation. Insert: Time distribution of the 1362.5 keV electron transition from which a half-life of 2.619(26) µsec is extracted.

TABLE I :

 I Experimental and shell model values for the excitation energies, in MeV, and reduced transition probabilities B(E2), in e 2 fm 4 , of 44 S.

	E/B(E2)	2 + 1		0 + 2		2 + 2		2 + 1 →0 + 1	2 + 1 →0 + 2
	exp.		1.329(1) 1.365(1) 2.335(39) 63(18) 8.4(26)
	SM		1.172		1.137	2.140		75	19
		5	6	+ 2		61		Ca 48	Ar 46	S 44	Si 42
		4		120			Gap	4.6		4.1	3.6	3.1
	E (MeV)	3	4	+ 2		58		Corr.	2.3 10.6 16.6 17.8
		2		110		2 + 2	-0.3
		1	2	+ 1	75	56	0 + 2	55 19		2	+ 1	63
		0	0 + 1						0 + 1
						SM

We are grateful to the GANIL staff and the LISE team for support. We acknowledge P. Van Isaker for fruitful discussions. This work was supported by CNCSIS-UEFISCSU, proj. numb. PNII-IDEI 933/2007, Academy of Sciences of Czech Rep. and by the European Community through OTKA K68801.


Philippe Jacquet

Salman Malik

Networks and Telecommunications

Keywords: single-hop network, multi-hop network, capacity, throughput, slotted ALOHA, TDMA, Grid based TDMA, square, hexagonal, honeycomb

We present our analysis and results that allow us to conjecture that maximum capacity in wireless networks can be achieved if nodes transmitting simultaneously are positioned in a hexagonal grid pattern. But obviously, it is very difficult to realize such a protocol which ensures that active transmitters in the network are positioned in any specific grid pattern. We compare the optimal capacity in networks with grid positioned transmitters with the capacity of wireless networks where nodes are dispatched according to uniform distribution and use very simple ALOHA-based protocol for channel access. We will also extend this analysis to multi-hop case and characterize the maximum throughput achievable in wireless networks with ALOHA-based protocols and TDMA and grid based TDMA protocols.

Protocoles optimaux de transmission dans les réseaux sans fil et comparaisons avec des protocoles basés sur ALOHA Résumé : Nous présentons une analyse et des résultats qui nous permettent de conjecturer que la capacité maximale dans les réseaux sans fil est obtenue quand les noeuds transmettant simultanément sont placés dans une configuration de grille hexagonale. Naturellement, il est très difficile de réaliser en pratique un protocole qui s'assure que des émetteurs actifs dans le réseau soient placés dans une topologie spécifique en grille. L'objet de ce rapport est de comparer la capacité optimale dans les réseaux avec des émetteurs placés en grille avec la capacité de réseaux sans fil où des noeuds sont placés selon la distribution uniforme et utilisent un protocole ALOHA très simple pour l'accès au medium. Nous étendons l'analyse au cas multi-sauts et nous caractérisons le capacité maximale utilisable dans les réseaux sans fil avec des protocoles ALOHA, des protocoles TDMA et des protocoles de TDMA organisś en grille.

Mots-clés : multi-saut réseau, capacité, protocoles TDMA, protocoles de TDMA organisés en grille

Introduction

Performance in multi-hop networks is in general poor compared to nominal bandwidth. This is mainly due to the interference management which can effect the spatial reuse of wireless channel over the multi-hop path. In future, wireless networks are expected to be deployed widely and the question is how these networks can satisfy the increasing demand on throughput under these constraints. Many of the proposed protocols for optimizing throughput in wireless networks use channel reservation or node coloring schemes, for example [START_REF] Bao | A new approach to channel access scheduling for ad hoc networks[END_REF][START_REF] Gandham | Link scheduling in sensor networks: Distributed edge coloring revisited[END_REF][START_REF] Rhee | Drand: distributed randomized tdma scheduling for wireless ad-hoc networks[END_REF][START_REF] Wu | D-lsma: Distributed link scheduling multiple access protocol for qos in ad-hoc networks[END_REF]. Such protocols can be very complex to implement in real world.

In this report we will be characterizing the throughput of wireless networks under different protocols which can be a challenging task especially under multihop conditions. We will investigate how TDMA based protocols compare with much simpler ALOHA-based protocols. We will also evaluate the grid based TDMA protocol which ensures that simultaneous transmitters are positioned in a specific grid pattern to give maximum throughput in a wireless network. Wireless networks of regular grid topologies are studied in, for example, [START_REF] Hong | Throughput analysis of large wireless networks with regular topologies[END_REF][START_REF] Liu | Throughput analysis of fading sensor networks with regular and random topologies[END_REF] and compared to networks with randomly dispatched nodes. In our analysis of grid based TDMA, only the transmitters accessing the channel simultaneously form a regular grid pattern. Our aim is to lay the groundwork for the design of simpler protocols for multi-hop adhoc networks and the main result from this report is that complex reservation protocols like TDMA or grid based TDMA do not significantly improve the throughput in wireless networks as compared to the simple slotted ALOHA protocol.

This report is organized as follows. In the next section, we present the comparison of Poisson distributed networks and grid-based networks in terms of average information rate that any arbitrarily placed node can receive. We limit our analysis to the case where a node can receive at most one packet at a time provided its signal-to-noise ratio (SNR) condition is fulfilled. We extend our analysis to multi-hop case and evaluate maximum throughput of different protocols in section 3. In the last section 4, we analyze the optimality of grid arranged field of transmitters.

Capacity with randomly and grid positioned transmitters

The authors of [START_REF] Gupta | The capacity of wireless networks[END_REF] have shown that information rate received by every node is finite and is irrespective of the network density. This analysis was done under the hypothesis that a node can only receive from one transmitter at a time. An accurate assessment of the capacity of a wireless network shall take into account:

1. Geometry for positioning of the nodes 2. Physics for wave propagation and attenuation in the medium and 3. Information theory for extraction of information from received signal.

Article [START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF] evaluates the capacity of a wireless network in a realistic model. From explicit formula of the Laplace transform of received signal distribution, it derives an explicit formula for the average information rate received by a

RR n°7339

random node under the assumption that transmitters are dispatched over an infinite plane according to the Poisson distribution. This model involves the above mentioned three aspects in the following manner:

1. Nodes are distributed over the network plane according to uniform distribution and use any nominal transmit power 2. Attenuation is related to distance: 1 r α and random fading and 3. Information is extracted from parallel superimposed signals in a similar way as in multiple-input-multiple-output technology (MIMO).

It can be argued that the optimal positioning of transmitters in a network shall be the grid formation. This optimal positioning can be achieved under a hypothetical medium access protocol but designing such a protocol is very difficult because of the limitations introduced by wave propagation characteristics and node distribution and is beyond the scope of this report. In this section we analyze the networks with transmitters optimally positioned in square, hexagonal and honeycomb grids. We compare these results with [START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF] and limit to the case when at most one packet, satisfying the SNR condition, is received successfully. We are not considering the MIMO situation where a node can receive several packets simultaneously. Our comparison is in terms of the average information rate received by an arbitrarily placed Access Point (AP) in the network. An important result from this analysis is that information rate received by the AP is maximum when transmitters are positioned in hexagonal grid and is comparable albeit slightly higher than the average information rate received if the transmitters are distributed randomly.

In section 2.1 we present the main results from [START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF]. Evaluation of grid-based networks is discussed in section 2.2 and results from our numerical simulations are presented in the section 2.3.

We assume that fading is zero and we restrict to uniform unit nominal transmit power.

Model for protocols based on randomly distributed transmitters

Consider an infinite network on 2D plane with nodes randomly distributed.

We consider the following model for this network:

1. time is slotted 2. nodes are synchronized and use slotted ALOHA protocol 3. in each slot, nodes which simultaneously transmit are given by a uniform Poisson distribution of mean λ transmitters per unit square area.

This model has been investigated in detail in [START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF]. In order to make this report self contained, we will summarize some of the analytical results from [START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF].

We assume that the attenuation coefficient (α) is greater than 2.

The signal level of a transmission received at distance r is:

1 r α .
If S is the set of the locations of all transmitters, and z i ∈ S is the transmitter whose SNR at location z should be at least equal to K then:

|z -z i | -α ≥ K j =i |z -z j | -α (1)
where |z| represents the Euclidean norm of a vector z = (x, y)

: |z| = x 2 + y 2 . or W (z, {z i }) ≥ KW (z, S -{z i }) where W (z, S) = zj ∈S |z -z j | -α .

Distribution of signal levels

Figure 1: Signal Levels (in dB's) for a random network with attenuation coefficient α = 2.5

Figure 1 shows the function W (z, S) for z varying in the plane with S an arbitrary set of Poisson distributed transmitters. Figure 1 uses α = 2.5. It is clear that closer the receiver is to the transmitter, larger is the SNR. For each value of K we can draw an area, around each transmitter, where its signal can be received with SNR greater or equal to K. Figure 2 shows reception areas for the same set S, as in Figure 1, for various values of K. As can be seen, the reception areas do not overlap for K > 1 since there is only one dominant signal. For each value of K we can draw, around each transmitter, the area where its signal is received with SNR greater or equal to K. The aim is to find the average size of this area and how it is a function of λ, K and α.

W (z, S) depends on S and hence is also a random variable. The random variable W (z, S) has a distribution which is invariant by translation and therefore does not depend on z. Let w(S) be its density function. If S is given by a 2D Poisson process with intensity λ transmitters per slot per unit square area, Laplace transform of w(S), w(θ, λ), can be computed exactly.

The Laplace transform, w(θ, λ) = exp( (e -θr -α -1)rdr, satisfies the identity: From the above formula and by applying inverse Laplace transformation:

w(θ, λ) = exp(πλΓ(1 - 2 α )θ 2 α ) RR n°7339
P (W (λ) < x) = 1 2iπ +i∞ -i∞ w(θ, λ) θ e θx dθ Expanding w(θ, λ) = n≥0 (-Cλ) n n! θ nγ P (W (λ) < x) = 1 2iπ n≥0 (-Cλ) n n! +i∞ -i∞ θ nγ-1 e θx dθ
Then by bending the integration path towards negative axis:

1 2iπ +i∞ -i∞ θ nγ-1 e θx dθ = e iπnγ -e -iπnγ 2iπ ∞ 0 θ nγ-1 e -θx dθ = sin(πnγ) π Γ(nγ)x -nγ
we get,

P (W (λ) < x) = n≥0 (-Cλ) n n! sin(πnγ) π Γ(nγ)x -nγ (2)
where γ = 2 α and Γ(.) is the Gamma function.

Reception areas

Let p(λ, r, K, α) be the probability to receive a signal sent at distance r with SNR at least equal to K. Then: p(λ, r, K, α) = P (W (λ) < r -α K ) and average

RR n°7339

size of the reception area around an arbitrary transmitter with SNR at least equal to K is: σ(λ, K, α) = 2π p(λ, r, K, α)rdr.

The average size of the reception area around an arbitrary transmitter i with SNR at least equal to K satisfies the identity:

σ i (λ, K, α) = 1 λ sin( 2 α π) 2 α π K -2 α ( 3
)
The reception area σ i (λ, K, α) is inversely proportional to the density of transmitters λ and the product λσ i (λ, K, α) is a function of K and α.

We notice that when α → ∞, σ i (λ, K, ∞) → 1 λ . This is due to the fact that when α is very large, closest node acting as source of interference gives the far largest estimate and consequently the area of reception turns to be Voronoi cell around each transmitter. This holds for all values of K. The average size of Voronoi cell being equal to the inverse density of the transmitters, 1 λ , we get the asymptotic result. Note that when K grows as exp(O(α)), we have σ i (λ, K, α) ≈ 1 λ exp(-2 α log(K)), which suggests that the typical SNR as α → ∞ is of the order of exp(O(α)).

Secondly, as α → 2, we have σ i (λ, K, 2) → 0 because sin( 2 α π) → 0. Indeed, the contribution of remote nodes tends to diverge and makes the SNR approach to zero. This explains why σ i (λ, K, 2) → 0 for any fixed value of K.

Capacity of wireless networks with Poisson distributed transmitters

Consider an AP at an arbitrary location z and let N (z, λ, K, α) denote the number of K-reception areas it belongs to. The quantity N (z, λ, K, α) is sometimes called the K-hand-over number of the AP. Following theorem has been proved in [START_REF] Baccelli | On a Coverage Process Ranging from the Boolean Model to the Poisson Voronoi Tessellation With Applications to Wireless Communications[END_REF][START_REF] Jacquet | Shannon capacity in poisson wireless network model[END_REF]:

E(N (z, λ, K, α)) = λσ i (λ, K, α) = σ i (1, K, α) (4)
E(N (z, λ, K, α)) represents the average number of associations of an arbitrarily placed AP in the network. Under the hypothesis that a node can only receive at most one packet at a time, the average information rate received by the AP is also equal to E(N (z, λ, K, α)).

Note that E(N (z, λ, K, α)) is independent of λ and

lim α→∞ E(N (z, λ, K, α)) → 1.

Model for protocols based on transmitters positioned in grid formation

Consider an infinite 2D plane with transmitters distributed like a sequence of points (z 1 , z 2 , z 3 , ...) and arranged in grid formation where z i = (x i , y i ) is the location of transmitter number i.

We consider the following model for this network:

1. time is slotted and 2. in each slot, all transmitters will transmit simultaneously. In our analysis, we have covered grid layouts of square, hexagonal and honeycomb as shown in figure 3. Grids are constructed from d which defines the minimum distance in-between neighboring transmitters and can be derived from hop-distance parameter of a typical TDMA-based protocol. The actual meaning of the value of d is beyond the scope of this section. The density of transmitters, λ, is the number of grid points (transmitters) per unit square area. λ will depend on the type of grid layout and will be inversely proportional to d 2 .

We are not aware of any analytical formula, as derived for the case of network with Poisson distributed transmitters, for the size of reception areas in grid networks. Therefore we have used a numerical method for computing the size of reception area around an arbitrary transmitter in the network.

Reception areas

Our aim is to find the area, the size of set of points A i (λ, K, α), where signal from transmitter i is received with SNR at least equal to K. Here noise is the sum of received signal power from all transmitters excluding i. Now there can be two approaches for computing the size of reception area. First approach uses the Monte Carlo simulation method while the second approach involves computing the size of reception area via numerical integration. While we have implemented both approaches, the problem with Monte Carlo simulation method is that we are not able to get sufficiently accurate results because of the sampling error. Therefore, herein we will only present the numerical integration approach.

The set A i (λ, K, α) contains the point z i since here the SNR is infinite.

Our aim is to compute

σ i (λ, K, α) = |A i (λ, K, α)|. If C(K, α)
is the closed curve that makes the boundary of σ i (λ, K, α) and z is a point on C(K, α), we have

σ i (λ, K, α) = 1 2 C(K,α) det(z -z i , dl) (5)
where det(a, b) is the determinant of vectors a and b. Here it gives the area of the parallelogram formed by the two vectors.

In equation 5, dl is the vector tangent to C(K, α) at point z. C(K, α) is assumed counter-clockwise such that area σ i (λ, K, α) is always on its left. Note From equation 1, SNR S i (z) of transmitter i at point z can be written as:

S i (z) = |z -z i | -α j =i |z -z j | -α (6)
We can assume that at point z, S i (z) = K.

On point z we can define the gradient of

S i (z), ∇S i (z) = ∂ ∂x S i (z) ∂ ∂y S i (z) .
Details on computing ∇S i (z) can be found in appendix.

The vector dl, tangent to

C(K, α) at point z, is co-linear with J ∇Si(z) ||∇Si(z)||
where J is the rotation matrix, J = 0 1 -1 0 and ||.|| is the magnitude of the vector (to align with the direction of C(K, α), ∇S i (z) is rotated clockwise). Therefore, we can fix dl = J ∇Si(z) ||∇Si(z)|| δt and in equation 5,

det(z -z i , dl) = -(z -z i ). ∇S i (z) ||∇S i (z)|| δt = -(x -x i ) ∂ ∂x S i (z) ||∇S i (z)|| δt -(y -y i ) ∂ ∂y S i (z) ||∇S i (z)|| δt
where δt → 0. The sequence of points z(k), computed as

RR n°7339 z(0) = z z(1) = z(0) + J ∇S i (z(0)) ||∇S i (z(0))|| δt z(k + 1) = z(k) + J ∇S i (z(k)) ||∇S i (z(k))|| δt describes C(K, α) when δt → 0.
Therefore the expression 5 reduces to,

σ i (λ, K, α) = 1 2 C(K,α) det(z -z i , dl) (7)
≈ -

1 2 k (z(k) -z i ). ∇S i (z(k)) ||∇S i (z(k))|| δt
assuming that we stop the sequence z(k) when it loops back on or close to the point z. Details on how to find the first point z(0) = z can be found in Appendix. Note that the negative sign in equation 7 is automatically negated by the dot product of vectors z(k)z i and ∇S i (z(k)).

Capacity of wireless networks with transmitters positioned in grid patterns

As we do not have a closed form expression for the reception area σ i (λ, K, α),

N (z, λ, K, α) = λσ i (λ, K, α)
is computed numerically according to the method described above.

Numerical simulations

In this section we will compare grid networks where transmitters are arranged in square, hexagonal and honeycomb patterns respectively with network where transmitters are Poisson distributed. The comparison is in terms of

E(N (z, λ, K, α)) = λσ i (λ, K, α).
For grid arranged transmitters, we performed numerical simulations in a very large network spread over 2D square area with length of each side equal to 5000 meters. Transmitters are spread over this area in square, hexagonal or honeycomb pattern. To avoid edge affects, the transmitter i, whose reception area we will compute, is located in the center at origin, z i = (x i , y i ) = (0, 0). The network area is large enough so that the transmitters on the boundary have almost negligible effect on the reception area of transmitter i. We set d equal to 30 meters although it will have no effect on the validity of our conclusions as λσ i (λ, K, α) is independent of λ. In our numerical simulations, we set δt = 0.01. λ depends on the type of grid and is computed from the total number of transmitters spread over the network area of 5000 x 5000 square meters.

In case of randomly distributed transmitters, λσ i (λ, K, α) = σ i (1, K, α) is computed from the analytic expression 3. x-axis (m) The shape of reception areas, in all three grid layouts, is influenced most by the closest source of interference but as the desired minimum K increases, influence from more distant transmitters, spread around the origin, increases and shape of reception areas approaches the shape of a small circular disk with transmitter in the center. Figure 8 shows the comparison of λσ i (λ, K, α) for Poisson distributed transmitters and grid arranged transmitters with varying K and fixed α = 4.0.

K = 1 K = 2 K = 6 K = 10 K = 100

Comparison with varying attenuation coefficient (α)

Figures 9, 10 and 11 show the variation in reception areas with increasing α.

Note that when α → ∞, the area of correct reception of a transmitter tends to be the Voronoi cell around this transmitter. This is because reception area is influenced most by the interference from the nearest transmitter as compared to any other transmitter. The average area of Voronoi cell equals 1 λ . As α → ∞, reception area in square grid layout approaches the shape of a square with area d 2 for whatever is the K. Similarly, for hexagonal grid formation, as α → ∞, reception area approaches the shape of a hexagon with area x-axis (m) We numerically computed λσ i (λ, K, α) with increasing α up to 100 and from the results, in figure 12, we can say that asymptotic result of Poisson distributed transmitters and grid arranged transmitters is same. where N is the number of nodes in the network. In this section we will extend our analysis to the multi-hop case. Extending the model of previous section to evaluate the capacity in multi-hop networks is difficult. Therefore, we will develop an analytical model for computing throughput of a wireless network under multi-hop configuration.

K = 1 K = 2 K = 6 K = 10 K = 100
In one of the first papers on throughput in multi-hop networks, [START_REF] Nelson | The spatial capacity of slotted aloha multihop packet radio network with capture[END_REF] analyzed slotted ALOHA protocol and showed that the critical performance parameter is p -probability of transmission in any time-slot.

The authors of [START_REF] Baccelli | An aloha protocol for multihop mobile wireless networks[END_REF] developed an analytical model of slotted ALOHA in multihop networks and showed that transport capacity can be optimized by selecting an appropriate value of p.

[9] investigated the single-hop and multi-hop throughput using slotted ALOHA protocol and analyzed the Performance Accident when multi-hop throughput reaches the single-hop throughput and optimal probability of channel access -p -shifts dramatically from 1 N to 1 log N . We will use our analytical model with simulator to compare maximum throughput that can be achieved in a network under Constrained TDMA, Constrained Grid TDMA and simple ALOHA-based protocols. We will describe the model of Constrained TDMA and Constrained Grid TDMA protocols later in section 3.1.2.

This section is organized as follows. In section 3.1 we present the details of the network model we have used. Analytical model for computing throughput in a network is discussed in the section 3.2. Section 3.3 gives detail of our simulation methodology and results from comparison of different protocols can be found in the section 3.4. We performed experiments for:

• Comparison of slotted ALOHA and Constrained TDMA

• Comparison of Constrained TDMA and Constrained Grid TDMA

System model

Following are the major assumptions of our model.

Physical model

We have used the same physical model as was used in section 2.

Details of this model are:

• All transmitting nodes have the same unit nominal transmission power.

• The signal received at distance r from the transmitter is 1 r α . We have assumed α = 4.0.

• A packet is correctly received if its SNR is greater or at least equal to K where noise is the sum of the signal powers received from all others RR n°7339 transmitters, acting as sources of interference, in the same time-slot. This is identical to the SNR model used in section 2 in equation 1:

|z -z i | -α j =i |z -z j | -α ≥ K 3.1.

Protocol models

• The time is slotted.

• In slotted ALOHA access protocol every node transmits a packet on any given slot with probability p, independently of the other nodes.

• In our analysis we do not consider any particular TDMA-based protocol. We only assume that same time-slot shall not be used within a certain neighborhood and refer to the protocols following this reservation rule as constrained TDMA protocols. The radius of this neighborhood is defined in terms of two parameters: m representing the number of hops and r representing the distance of 1-hop in euclidean space (such that all nodes lying at distance less than r are at 1-hop). Therefore a time slot cannot be shared within a distance of less than m.r or, in other words, nodes transmitting in the same slot shall be located at a distance greater or equal to m.r from each other.

In this study we consider the following two types of constrained TDMA protocols:

1. Constrained TDMA In Constrained TDMA protocol, the list of simultaneously transmitting nodes is built by random draws while following the above specified reservation rule for channel access. This is similar to the TDMA model use in [START_REF] Banaouas | Performance evaluation of tdma versus csma-based protocols in sinr models[END_REF]. The set of transmitters S which are transmitting in the same time slot is built as follows:

Constrained Grid TDMA

In Constrained Grid TDMA protocol, the nodes which transmit simultaneously in a time slot form a specific grid pattern. The dimension d of this grid pattern is related to the parameters, m and r, as:

d = m.r
Figure 3 shows the realization of square, hexagonal and honeycomb grid patterns from dimension d.

In section 2 we have shown that hexagonal grid pattern gives the maximum local capacity i.e. the information rate received by an arbitrarily placed node. Using the hypothesis that maximum local capacity will also result in maximum multi-hop capacity1 , in this section we will evaluate the hexagonal grid pattern only. The set of transmitters S is built as follows.

(a) Initialize: N = {s 1 , s 2 , ..., s N } and S = ∅. (b) Randomly select a node s i from N and add it to the set S: S := S ∪ {s i }. (c) Construct a virtual hexagonal grid, from the position of s i , with dimension d. (d) All nodes which overlap with the points of virtual hexagonal grid will be added to the set S.

Traffic model

• The network throughput is computed under the assumption that traffic model is uniform and every node sends equal traffic to every other node in the network.

• We assume that every node has an infinite transmit buffer filled with packets (since we want to compute the ultimate throughput).

Network map

The network map is a square area A = 100 × 100m 2 with N nodes distributed over it.

Throughput computation

Here we will present the analytical model for computing throughput of the network in terms of packets per time-slot. The goal of the simulation run with N nodes can be limited to the measurement of the Average Reception Matrix. The average reception matrix C is the matrix whose coefficient: c ij represents the proportion of transmissions from node i which are successfully received by node j.

The average number of retransmissions, including the first transmission, required for delivering a packet from node i to node j will be 1 cij . Let us call the matrix 1 cij the transmission cost matrix D ij where D ij = 1 cij (we set D ii = 0 for all i). Therefore the average number of retransmissions needed to deliver a packet from node i to node j is D ij .

Each node i computes its best route to any destination j by performing Dijkstra algorithm on the matrix D. The best route is the route (i 1 , i 2 , ..., i k ) such that i 1 = i and i k = j and D i l D i l+1 is minimal. Let M ij be this minimal value and it will represent the minimum optimal number of retransmissions from node i to node j. We use the (min, +) algorithm to compute the min hop matrix M . We set M = D and iterate M ← M * D until M is stationary, with (M * D) ij = min l (M il + D lj ). There could be only a maximum of log N iterations.

Let R i be the traffic rate of node i and it is equal to the proportion of slots when node i is active. In case of slotted ALOHA, R i should be equal to the probability of medium access, p.

During T slots, there are on average T R i packet transmission attempts in the network. If ρ is the throughput of the system in packets per time-slot, the number of successfully delivered packets is ρT . Since each node sends equal traffic to every other node and there are N (N -1) pairs of source and destination, the number of packets delivered from source i to destination j during T slots is ρT N (N -1) .

The transmission cost from source i to destination j is M ij , therefore ρ i,j Mij N (N -1) T is equal to the average number of transmission attempts in the network during T slots and should be equal to T R i . Therefore, throughput of the network is:

ρ = N (N -1) R i i,j M ij
In case of slotted ALOHA, R i = pN and net throughput expression reduces to: ρ = N 2 p(N -1) i,j

Mij

which is similar to the analytic expression, derived in [START_REF] Jacquet | The Performance Accident in wireless multi-hop networks[END_REF], for throughput with slotted ALOHA protocol.

Simulation technique

We have used two different schemes for distribution of N nodes in the network area A. In the first scheme, N nodes are dispatched in the network area according to Poisson distribution and throughput of this network is computed with slotted ALOHA and Constrained TDMA medium access protocols.

Poisson distributed nodes may not give the maximum throughput with Constrained Grid TDMA protocol especially at low network densities. At lower network densities, many points of virtual grid, constructed from any randomly selected node, may not overlap with any other node and hence the number of active transmitters in the network, in any given slot, will be lower than the optimal number possible within the constraints of Constrained Grid TDMA protocol. Since we are interested in finding the ultimate throughput, we have devised another scheme for dispatching nodes in the network area A. Consider a hexagonal grid pattern spread over an infinite area also containing the network area A. Depending on d, let's say maximum n number of nodes will lie within A. We may call this hexagonal grid pattern, with n nodes in A, a grid frame.

There shall be at least N n such grid frames, with each frame contributing a maximum of n nodes to the network, in order to have N nodes in A.

All these grid frames shall be superimposed over each other such that each frame is displaced along x-and y-axes. This displacement is uniformly selected from the interval [-2d, +2d].

Figure 13 shows the distribution of 326 and 969 nodes in the network area of 100 x 100m 2 with d = m.r = 30m (m = 2, d = 15m). Note that at lower network densities, there will be fewer grid frames and clustering of nodes can be observed but as the network density increases, nodes appear well spread over the network area.

The average reception matrix C is computed via simulation of packet transmissions. Each node, when transmitting, sends out a broadcast packet in the network and, from SNR model, we can compute the coefficients c ij of the matrix C. Note that in case of Constrained Grid TDMA protocol, the received signal and interference levels, for the source and destination pair (i, j), remains constant and c ij is computed numerically.

Results

The transmission probability and number of nodes in the network are varied in order to understand the nature of throughput w.r.t. these parameters. Note that in case of constrained TDMA based protocols, transmission probability implicitly depends on m. Therefore for constrained TDMA based protocols we will vary m and investigate its effect on the network throughput.

Slotted ALOHA and Constrained TDMA protocols

We know that the single-hop throughput in a network with slotted ALOHA protocol is ρ s ≥ pN (1p) N -1 . Therefore an optimal value of p is ≈ 1 N which

RR n°7339

would give ρ s ≈ 1 e . Note that single-hop throughput is slightly higher than 1 e because of the spatial reuse in the network when more than one simultaneous transmissions can be successful. Since single-hop communication is a special case of multi-hop communication, multi-hop throughput ρ m ≥ ρ s . Gupta and Kumar in [START_REF] Gupta | The capacity of wireless networks[END_REF] placed an upper bound on the throughput of multi-hop adhoc networks in the order of N log N . To achieve this asymptotic behavior, we expect to get maximum throughput with p of the order of 1 log N when N is large. Figure 14 shows the throughput achieved with slotted ALOHA protocol in network with N Poisson distributed nodes. It can be observed that throughput increases initially, followed by a decrease and then an increase with rate in the order of √ N . The first local maxima observed at p ≈ 1 N occurs because of the single-hop nature of the transmissions. This accounts for the initial increase and subsequent decrease in throughput. As the number of nodes in the network increase, a second maxima starts to develop. This second maxima increases with increasing number of nodes and occurs at p ≈ c log N where c ≈ 0.25. Flat nature of the curvature at second maxima makes it difficult to get the accurate value of c. The dramatic shift of optimal value of p from 1 N to c log N is called performance accident and it can be observed easily in figure 15 when optimal value of p shifts dramatically for N ≈ 113. Performance accident is studied in detail in [START_REF] Jacquet | The Performance Accident in wireless multi-hop networks[END_REF]. Figure 16 shows the maximum throughput achieved by slotted ALOHA protocol with optimal value of p. Note that before the occurrence of performance accident, maximum throughput remains constant at a value slightly higher than 1 e . This almost constant maximum throughput can be explained by the singlehop nature of the transmissions. As the number of nodes in the network increase, the multi-hop throughput out-performs the single-hop throughput.

For Constrained TDMA protocol, we have set r = 15m and m is varied from 2 to 4. Under these conditions, the nodes transmitting simultaneously will be at least at a distance of d = m.r from each other. Figure 17 shows throughput in the same network with Constrained TDMA protocol. It can be observed that at very low network densities, throughput with m = 3 or m = 4 is higher than the throughput with m = 2. The reason can be the sparsely located nodes and higher number of active transmitters with m = 2. This may reduce the number of available intermediate nodes required for multi-hop communication from source to destination but as the network density increases, throughput with m = 2 out-performs throughput with higher values of m.

Constrained TDMA and Constrained Grid TDMA protocols

For comparing Constrained TDMA protocol with Constrained Grid TDMA protocol we have used the grid based distribution of nodes. We have set r = 15m and m = 2. Each grid frame is therefore constructed with the dimension d = m.r = 30m. Note that as the node density increases, throughput of Constrained Grid TDMA protocol is observed to be slightly higher than the typical Constrained TDMA protocol. This can be explained by the corresponding increase in the average density of successful receivers in the coverage area of any transmitter resulting in higher local capacity, λσ(λ, K, α), in grid based networks. This result also verifies our hypothesis that higher local capacity in grid based networks will result in higher multi-hop capacity.

Field optimization for wireless capacity

Consider an infinite set of points S = {z 1 , z 2 , . . . , z n , . . .} on the plane. We denote

x i (z) = |z-zi| -α j |z-zj | -α . Let α > 2.
In order to simplify notations, we will remove the reference to z when no ambiguity is possible. We define a function f (x i ) which can be continuous or integrable. For instance we will use f (x i ) = 1 xi>K for some given K.

We suppose that the set density ν(S) has a limit:

ν(S) = lim R→∞ 1 πR 2 i 1 |zi|≤R .

RR n°7339

We denote h(z) = i f (x i ) and we want to optimize E(h(z)) defined by the limit:

E(h(z)) = lim R→∞ 1 πR 2 |z|≤R h(z)dz 2 .
We denote σ i = f (x i )dz 2 and we have

E(h(z)) = lim R→∞ 1 πR 2 i 1 |zi|≤R σ i = ν(S)E(σ i ) with E(σ i ) = lim n→∞ 1 n i≤n σ i .
We denote ∇ i the operator of differentiation with respect to z i .

We have for i = j:

∇ i x j = αx i x j z -z i |z -z i | 2 and ∇ i x i = α(x 2 i -x i ) z -z i |z -z i | 2 . Therefore, ∇ i h(z) = αx i z -z i |z -z i | 2   -f ′ (x i ) + j x j f ′ (x j )   .
Although h(z)dz 2 = ∞, we nevertheless have a finite ∇ i h(z)dz 2 . In other words, the sum j ∇ i σ j converges for all i.

We also know that for all j in S, i ∇ i σ j = 0. Indeed this would be the differentiation of σ j when all points in S are translated by the same vector. Similarly

i ∇ i h(z)dz 2 = 0.
If S is a grid, then: ∇ i h(z)dz 2 = j ∇ i σ j and therefore would be null

(since i ∇ i h(z)dz 2 = 0).
We could conclude erroneously that:

• all grid sets are optimal;

• all grid sets give the same E(h(z)).

In fact this is wrong: we could also conclude that E(σ i ) does not vary but this will contradict that ν(S) must vary. The reason of this error is that a grid set cannot be modified into another grid set with uniformly bounded transformation, unless the two grid sets are just simply translated by a simple vector.

However we prove that the grid sets are locally optimal within sets that can be uniformly transformed between each other.

In order to cope with uniform transformation and to be able to transform a grid set to another grid set we introduce the linear group transformation. We only consider the grid situation, therefore:

E(h(z)) = σ 0 = f (x 0 )dz 2 .

RR n°7339

Linear group transformation

Here we assume that the points in z are modified according to a continuous linear transform M (t) where M(t) is a matrix with M(0) = I, for example M(t) = I + tA with A a matrix. Under this assumption we have

∂ ∂t σ 0 = i (Az i .∇ i σ 0 )
Therefore the derivative of σ 0 with respect to matrix A is

i Az i .∇ i σ 0 = tr(A( i z i ⊗ ∇ i σ 0 )).
In other words the derivative with respect to matrix is

exactly D = i z i ⊗ ∇ i σ 0 . D = D xx D xy D yx D yy
Let assume that M(t) = (1 + t)I, i.e. the linear transform is homothetic. In this case we have σ 0 (t) = (1 + t) 2 σ 0 and A = I. As a first property we have tr(D) = 2σ 0 since the derivative of σ 0 with respect to identity matrix I is exactly 2σ 0 (since σ ′ 0 (0) = 2σ 0 ). Less obvious is the fact that D is a symmetric matrix. The easiest proof of this property is to consider the derivative with respect to matrix J:

J = 0 -1 1 0
which is zero since J is the initial derivative for a rotation. Since tr(JD) = D yx -D xy , this implies that D is symmetric. Let T be defined as:

T = dz 2 i (z -z i ) ⊗ ∇ i f (x 0 ) .
The sum i (zz i ) ⊗ ∇ i f (x 0 ) leads to a symmetric matrix since

T = α f ′ (x 0 ) ×   x 2 0 -x 0 z -z 0 2 (z -z 0 ) ⊗ (z -z 0 ) + i =0 x 0 x i z -z i 2 (z -z i ) ⊗ (z -z i )   dz 2
and the left hand side is made of (zz i ) ⊗ (zz i ) that are symmetric matrices.

Since,

D = i z i ⊗ ∇ i σ 0 = dz 2 i z i ⊗ ∇ i f (x 0 )
we have

D + T = dz 2 i z ⊗ ∇ i f (x 0 ) RR n°7339 It remains to prove that i z ⊗ ∇ i f (x 0 )dz 2 is symmetric too. Since i ∇ i f (x 0 ) = -∇f (x 0 ) we have i z ⊗ ∇ i f (x 0 )dz 2 = - x ∂ ∂x f (x 0 )dxdy x ∂ ∂y f (x 0 )dxdy y ∂ ∂x f (x 0 )dxdy y ∂ ∂y f (x 0 )dxdy = σ 0 0 0 σ 0 = σ 0 I which is symmetric.
We know that T is symmetric and D + T = σ 0 I. Thus tr(T) = 0 When grid is optimal, we must have T = 0. In any case, the matrix T must be invariant with respect to isometric symmetries of the grid. For the square grid, we have all symmetries with respect to any horizontal or vertical axes of the grid and in particular with the rotation of π 2 represented by J. Therefore the Eigen system must be invariant by rotation of π 2 . This implies that the eigenvalues are the same and therefore null since tr(T) = 0.

Same argument for the hexagonal grid, with the invariance for π 3 rotation and for the honeycomb pattern with invariance in π 6 rotation. In both cases T = 0.

Of course this does not exclude other grids or patterns with T = 02 . Let r i = |zz i | and r j = |zz j |. The gradients of f and g will be:

Appendix

∇f = -α|z -z i | -α-1 z -z i |z -z i | = -α 1 r α+1 i z -z i r i = -α z -z i r α+2 i ∇g = j =i -α z -z j r α+2 j
Using identity ∇ f g = 1 g 2 (g∇ff ∇g), we get The SNR S i (z) at point z should be greater or at least equal to K. We assume that it is equal to K. Then using equation 6

∇S i (z) = ∂ ∂x S i (z)
|z-zi|

j =i
|z-zj | = K Our aim is to find the coordinates of point z which satisfy above relation. To simplify computation of point z on the closed curve, bounding reception area of transmitter i located at z i = (x i , y i ), its y coordinate can be fixed such that z = (x, y) = (x, y i ). This reduces the SNR equation to

|x -x i | j =i |z -z j | -K = 0 ( 8
)
Equation 8 is a function of variable x and can be solved using Newton's Method.

Remark: Newton's Method: Given a function f (x) and its derivative f ′ (x) begin with a first guess x 0 . Provided the function is reasonably well-behaved a better approximation x 1 is := x 0 -f (x0) f ′ (x0) . The process is repeated until a sufficiently accurate value is reached: The first derivative of the function f (x) is,

x n+1 = x n - f (x n ) f ′ (x n ) (9)
f ′ (x) = 1 h 2 (h d dx g -g d dx h) = 1 j =i r -α j   ( j =i r -α j )(-α x -x i r α+2 i ) -(r -α i )(-α z -z j r α+2 j )  
Newton's Method also requires first approximation of the root, x 0 . An approximate value, closer to the actual root, can significantly reduce the number of iterations in Newton's Method. The location of transmitters i and j and point z in the plane form three corners of a triangle with angle (θ) equal to 0 in case of square and hexagonal grid and 30 degrees is case of honeycomb grid layout. Figure 19 shows the location of transmitters z i and z j , point z and distances r i , r j and d. Using above relation between r i and r j and Law of Cosines, we get the solution of r i as

r i = -B± √ B 2 -4AC 2A where A = 1 -K 2 α B = -2.d.cos(θ) and C = d 2 .
d is the distance between transmitters i and j and is known parameter of the grid layout.

Remark: Select positive value of r i as the solution of the above quadratic equation. Using x 0 = x i + r i as the first approximate solution in Newton's Method (equation 9), and after a few iterations, we can get a sufficiently accurate value x n+1 which will be the x coordinate of the point z.

Figure 2 :

 2 Figure 2: Distribution of reception areas for various value of SNR. K = 1, 4, 10 for situation of figure 1

Figure 3 :

 3 Figure 3: Square, Hexagonal and Honey-Comb grid layouts

Figure 4 :

 4 Figure 4: Computation of size of reception area of transmitter i

Figure 5 :

 5 Figure 5: Reception Area in Square Grid network with varying K -d = 30m, α = 4.0

√ 3 2 d 2 3 4 d 2 .

 3232 and for honeycomb grid reception area around transmitter approaches the shape of an equilateral triangle with area 3 √ Figure12shows the comparison of grid arranged transmitters with the case of transmitters distributed according to Poisson distribution.For Poisson distributed transmitters, equation 3 implies that:

Figure 6 :Figure 7 :

 67 Figure 6: Reception Area in Hexagonal Grid network with varying K -d = 30m, α = 4.0

Figure 8 :

 8 Figure 8: Comparison of Poisson distributed transmitters and Grid arranged transmitters in terms of λσ i (λ, K, 4.0) with fixed α = 4.0

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Reception Area in Square Grid network with varying α -d = 30m, K = 10.0

Figure 12 :3

 12 Figure 12: Comparison of Poisson distributed transmitters and Grid arranged transmitters in terms of λσ i (λ, 10.0, α) with fixed K = 10.0

  (a) Initialize: N = {s 1 , s 2 , ..., s N } and S = ∅. (b) Randomly select a node s i from N and add it to the set S: S := S ∪ {s i }. (c) Remove s i from the set N . (d) Remove all nodes from the set N which are at distance less than m.r from s i . (e) If set N is non-empty, repeat from step (c).

Figure 13 :

 13 Figure 13: Location distribution of N nodes on 100 x 100m 2 area using grid based distribution scheme

Figure 14 :Figure 15 :

 1415 Figure 14: Throughput with slotted ALOHA protocol. Nodes are dispatched in network area of 100 x 100m 2 according to Poisson distribution.

Figure 16 :

 16 Figure 16: Throughput at optimal p with slotted ALOHA protocol. Nodes are dispatched in network area of 100 x 100m 2 according to Poisson distribution.

Figure 17 :

 17 Figure 17: Throughput with Constrained TDMA protocol with varying m. Parameters: r = 15m and d = m.r. Nodes are dispatched in network area of 100 x 100m 2 according to Poisson distribution

Figure 18 :

 18 Figure 18: Throughput comparison of Constrained TDMA and Constrained Grid TDMA protocols. Parameters: m = 2, r = 15m and d = m.r = 30m. Nodes are spread over an area of 100.100m 2 based on grid based distribution

( a )

 a Computation of ∇S i (z) From equation 6 S i (z) = |z-zi| -α j =i |z-zj |-α = f g where, f = |zz i | -α and g = j =i |zz i | -α .

  Locating the starting point z on closed curve bounding reception area

From 8 ,

 8 f (x) = |x-xi| -α j =i |z-zj | -α -K Assume, g = r -α i and h = j =i r -α j where, r i = |xx i | and r j = |zz j |.Newton's Method requires the derivative of f (x) which is computed as below. d dx g = -α x-

  In case of square and hexagonal grid networks (b) In case of honeycomb grid network

Figure 19 :

 19 Figure 19: Geometric representation of finding the first approximation of point z = (x 0 , y i )

RR n°7339

This can be verified using the simulation model we have developed in this section.RR n°7339

This remains to be proven.RR n°7339

The coordinates of point z will be: (x n+1 , y i ).

RR n°7339