D J Durian

H Bideaud

P Duringer

A Schröder

C M Marques

The shape and erosion of pebbles

Keywords: numbers: 45.70.-n, 83.80.Nb, 91.60.-x, 02.60.Jh

The shapes of flat pebbles may be characterized in terms of the statistical distribution of curvatures measured along their contours. We illustrate this new method for clay pebbles eroded in a controlled laboratory apparatus, and also for naturally-occurring rip-up clasts formed and eroded in the Mont St.-Michel bay. We find that the curvature distribution allows finer discrimination than traditional measures of aspect ratios. Furthermore, it connects to the microscopic action of erosion processes that are typically faster at protruding regions of high curvature. We discuss in detail how the curvature may be reliable deduced from digital photographs.

I. INTRODUCTION

The roundness of pebbles on a beach has long been a source of wonder and astonishment for scientists in many fields [1,2]. Explanations for the pebble shapes were born from the simple pleasure of understanding nature but also from the hope that a pebble, or a collection of pebbles, might carry lithographically imprinted the signature of their erosion history. Reading that imprint would then, for instance, reveal if a pebble was eroded on a beach, a river or a glacier, or if it traveled a long distance down a stream. It even perhaps would reveal for how long the erosion forces have been at work on that object. Of obvious interest in Geology [3], a physical understanding of the formation of erosion shapes would also allow for a better control of many industrial processes leading to rounded objects such as gem stone or clay bead grinding in tumblers or fruit and vegetable peeling in several mechanical devices. Diverse mathematical tools have been developed for geometrical shape analysis of crystalites, cell membranes, and other far from equilibrium systems [4][START_REF] Lipowsky | Structure and dynamics of membranes, Handbook of biological physics[END_REF][START_REF] Meakin | Fractals, scaling and growth far from equilibrium[END_REF][START_REF]Growth, dissolution and pattern formation in geosystems[END_REF]; however, these do not seem applicable to pebbles.

The evolution of a pebble shape under erosion can arguably be viewed as a succession of elementary cuts that act at the surface of the body to remove a given amount of material. This converts young, polyhedral-like shapes with a relatively small number of large sides and sharp vertices into more mature shapes with a high number of small sides and smooth vertices. The size and the shapes of each of these successive ablations, as well as the surface sites where the cutting happens, are determined both by the conditions under which erosion takes place and by the nature of the material being eroded. Exposure of a young, polyhedral-like shape to the rough tumbling of a steep stream slope will result in relatively large cuts of the angular sections, while exposure to the gentle erosion of wind or water is more likely to lead to small cuts almost parallel to the existing flat sides. Also, the same sequence of external forces acting on two identical original shapes of different materials will result into distinct forms due to weight, hardness or anisotropy differences. In spite of the diversity of factors at play in shape modification, the complete evolution of the pebble shape is fully determined by (i) the initial form described by some number of faces, edges and vertices and (ii) the position, size and orientation of the successive ablations.

Given that the erosion process evolves by a succession of localized events on the pebble surface, it is surprising that the majority of the precedent attempts to characterize the pebble shapes were restricted to the determination of global quantities such as the pebble mass or the lengths of its three main axes [3]. Clearly, in order to capture both the local nature of the erosion process and the statistical character of the successive elementary cuts, one needs to build a new detailed description of the pebble shapes based on quantities that are more microscopic and more closely connected to evolution processes. In Ref. [START_REF] Durian | [END_REF] we proposed curvature as a key microscopic variable, since, intuitively, protruding regions with large curvature erode faster than flatter regions of small curvature. We then proposed the distribution of curvature around a flat, two-dimensional, pebble as a new statistical tool for shape description. And finally we illustrated and tested these ideas by measuring and modeling the erosion of clay pebbles in a controlled laboratory apparatus.

In this paper we elaborate on our initial Letter [START_REF] Durian | [END_REF], and we apply our methods to naturally-occurring ripup clasts found in the tidal flats of the Mont St.-Michel bay. Section II begins with a survey of shape quantification for two-dimensional objects, in general, and recapitulates our new curvature-based method. Section III provides further details of the laboratory experiments on clay pebbles. Section IV presents a new field study of the Mont St.-Michel rip-up clasts. And finally, following the conclusion, two methods are presented in the Appendix for reliably extracting the local curvature from digital photographs.

II. 2D SHAPE QUANTIFICATION

The issue of rock shape is of long-standing interest in the field of sedimentology [9][10][11][12][13][14][15][16][17]. Two basic methods have become sufficiently well established as to be discussed in introductory textbooks [3]. The simplest is a visual chart for comparing a given rock against a standard sequence of rocks that vary in their sphericity and angularity. A rock has high "sphericity" if its three dimensions are nearly equal. And it is "very angular", independent of its sphericity, if the surface has cusps or sharp ridges; the opposite of "very angular" is "well rounded". While useful for exposition, such verbal distinctions are subjective and irreproducible. The second method is to form dimensionless shape indices based on the lengths of three orthogonal axes. From the ratios, and the ratios of differences, of the long to intermediate to short axes, one can readily distinguish rods from discs from spheres. And a given rock may be represented by a point on a triangular diagram according to the values of three such indices, with rod / disc / sphere attained at the corners. This practice is nearly half a century old [18]. Nevertheless, there is still much debate about which of the infinite number of possible shape indices are most useful [19][20][21][22][START_REF] Graham | Earth Surface Processes and Landforms[END_REF]. In any case, such indices cannot capture fine distinctions in shape, let alone the verbal distinctions of angular vs rounded. Furthermore, they provide no natural connection to the underlying physical process by which the rock was formed.

More recent methods of shape analysis employ Fourier [START_REF] Schwarcz | [END_REF][25][26][27][28], or even wavelet [29], transforms of the contour. This applies naturally to flat pebbles or grains, but also to flat images of three-dimensional objects. The advantage of Fourier analysis over shape indices is that, with enough terms in the series, the exact pebble contour can be reproduced. For simple shapes, the contour may be described in polar coordinates by radius (eg distance from center of mass) vs angle, r(θ), and the corresponding transform. However, this representation is not single-valued for complex shapes with pits or overhangs. Generally, the contour may be described by Cartesian coordinates vs arclength, {x(s), y(s)}, and the corresponding transforms. In any case, the relative amplitudes of different harmonics give an indication of shape in terms of roughness at different length scales. In a different area of science, Fourier representations have proven especially useful for analysis of fluctuations and instabilities of liquid interfaces, membranes, etc. [4][START_REF] Lipowsky | Structure and dynamics of membranes, Handbook of biological physics[END_REF][START_REF] Meakin | Fractals, scaling and growth far from equilibrium[END_REF]. In practice, for these systems, shape fluctuations are sampled during some time interval and then the average Fourier amplitudes extracted by averaging over many different realizations of the shape. Also, because these phenomena are linear, each Fourier component grows or shrinks at some amplitude-independent rate and the evolution is fully determined by a dispersion relation. Unfortunately these features do not hold for the erosion of pebbles. Because each pebble shape only provides one configuration, average quantities need to be built from a different pre-scription. Also, there is no a priori guaranty that the variables are Gaussian distributed, and one needs a direct space method to better assess the importance of nonlinear phenomena. Non-linearity is, we suspect, intrinsically embedded in the erosion mechanisms of pebbles. If one considers for instance a shape represented by a single harmonic in the r(θ) representation, it is clear that the peaks will wear more rapidly than the valleys. Therefore the erosion rate cannot be a function of the harmonic number only; it must either be a non-linear function of the amplitude itself or a function coupling many harmonics.

Our aim is to provide an alternative measure of pebble shape that is well-defined, simple, and connects naturally to local properties involved in the evolution process. We restrict our attention to flat pebbles, where an obvious shape index is the aspect ratio of long to short axes. Since erosion processes generally act most strongly on the rough, pointed portions of a rock, we will focus on the local curvature of the pebble contour. Technically, curvature is a vector given by K = dT/ds, the derivative of the unit tangent vector with respect to arclength along the contour [START_REF] Weisstein | The CRC concise encyclopedia of mathematics[END_REF]. More intuitively, the magnitude of the curvature is the reciprocal of the radius of a circle that mimics the local behavior of the contour. Here we shall adopt the sign convention K > 0 where the contour is convex (as at the tip of a bump) and K < 0 where the contour is concave (as where a chip or bite has been removed from an otherwise round pebble). In the Appendix, we describe two means by which the curvature may be reliably measured at each point along the pebble contour. Note that the average curvature is simply related to the perimeter of the contour:

P = 2π/ K , (1)
which is obviously correct when the shape is a circle. To describe the shape of a pebble, a very natural quantity is the distribution of curvatures, ρ(K), defined such that ρ(K)dK is the probability that the curvature at some point along the contour lies between K and K +dK [START_REF] Durian | [END_REF]. In order to distinguish different distributions, as a practical matter, it is more reliable [START_REF] Press | Numerical Recipes in C[END_REF] to use the cumulative distribution of curvatures

f (K) = K 0 ρ(K ′ )dK ′ (2)
Literally, f (K) is the fraction of the perimeter with curvature less than K. Note that f (K) increases from 0 to 1 as K : 0 → ∞; the minimum curvature is where f (K) first rises above 0, the maximum curvature is where f (K) first reaches 1, and the median curvature is where f (K) = 1/2. Unlike for ρ(K), it is not necessary to bin the curvature data in order to deduce f (K). Instead, just sort the curvature data from smallest to largest and keep a running sum of the arclength segments, normalized by perimeter. Finally, so that the shapes of pebbles of different sizes may be compared, it is useful to remove the and (c) fraction of the perimeter f (K) with curvature less than K, vs K divided by the average curvature K , for a superellipse, oval, ellipse, and circle. The curve types match those for the shapes as shown in the inset. Note that, except for the circle, all have the same aspect ratio a = 3/2. The differences in shape are reflected in differences in the forms of f (K).

scale factor K , which is related to the total perimeter as noted above in Eq. (1). Altogether, we thus propose to quantify pebble shape by examining f (K) as a function of K/ K . To help build intuition, examples of f (K) are given in Fig. 1 for a few simple shapes. The simplest of all is a circle, where the curvature is the same at each point along the contour: K = K = 2π/P = 1/R. Thus f (K) = 0(1) for K < (>) K . The curvature distribution is the derivative of this step function, giving ρ(K) = δ(K -K ) as required. The other three shapes shown in Fig. 1 are all tangent at four points to the rectangular region {|x| < a = 3/2, |y| < 1}: an ellipse (x/a) 2 + y 2 =1; a superellipse (x/a) 4 + y 4 = 1; and an oval [x/a + (1 -1/a)y 2 /2] 2 + y 2 = 1. For a general ellipse, with long and short axes a and b respectively, one may compute

f (K) = E(sin -1 [1/ǫ(1-(2 √ 1 -ǫ 2 E(ǫ 2 )/(πK)) 2/3 ) 1/2 ], ǫ 2 )/E(ǫ 2 ) where ǫ = 1 -b 2 /a 2 is the ellipticity, E(x)
is the complete elliptical integral of the first kind and E(x, m) is the incomplete elliptic integral of the second kind [START_REF] Weisstein | The CRC concise encyclopedia of mathematics[END_REF].

While the ellipse, superellipse, and oval in Fig. 1 all have the same aspect ratio, a = 3/2, their shapes are obviously different. This emphasizes how a single number is insufficient to quantify shape. The shape differences do, however, show up nicely in the forms of f (K). The ellipse is closest to a circle, with a distribution of curvatures that is most narrowly distributed around the average and hence with an f (K) that is most like a step function. The superellipse is farthest from a circle, with four long nearly-flat sections and four high-curvature corners; its curvature distribution is broadest. The oval is intermediate.

III. LABORATORY EXPERIMENT

The examples given in Fig. 1 correspond to regular, highly symmetric shapes of two dimensional convex "pebbles". In practice, natural or artificial erosion processes lead to curvature functions with an important statistical component. In this section we elaborate on the laboratory experiments of Ref. [START_REF] Durian | [END_REF], designed to study both the statistical nature of the curvature distribution and the influence of the original shapes of the pebbles on the final output of a controlled erosion process.

Laboratory pebbles were formed from "chamotte" clay, a kind of clay made from Kaolin and purchased from Graphigro, France. The water content of the purchased clay was 22% in a state that could easily be kneaded. The clay was kept tightly packed before use in order to avoid water evaporation. Clay pebbles were produced using aluminium molds made in our laboratory. They consist of a polygonal well of 0.5 cm depth. Once the mould was filled with clay, it was left at rest for 24 hours, so that 98.5% of the water was removed by evaporation. All the experimental results presented here concern one day old pebbles. We noticed that pebbles older than 2 days were too fragile for our experimental configuration. The number of samples and the dimensions of the various pebbles are as follows: four squares of side 5 cm, five rectangles of sides 4 × 6 cm 2 , five regular pentagons of side 4.25 cm, one triangle with sides {7 cm, 7.5 cm, 9 cm}, one irregular polygon with 7 sides, one lozenge with acute angles of 45 • and four sides of 5 cm, and one circle of diameter 7 cm.

The wearing method that we chose relies on placing a pebble in the rotating apparatus sketched in Fig. 2. The apparatus is a square basin, of dimensions 30×30×7 cm 3 . The basin bottom is a 1 cm thick aluminium plate and the walls are made of 0.04 cm thick aluminium sheets. This rotating plate is fixed to a rod held by the jaws of a laboratory mixer, Heidolph RZR1. The mixer itself is fixed to a tripod, so that its inclination angle can be varied.

A typical trajectory of the pebble during the continuous rotation of the plate can be described as follows. First the pebble rotates with the basin until it reaches a high position. After the plate has rotated an angle between π/2 and π, the pebble begins to slide due to gravity, until it hits one of the walls in the bottom part of the container, and then rolls down along that wall as the basin keeps its rotation. After a short stop at a container corner the pebble starts a new cycle again. We performed preliminary tests in order to determine both ideal basin orientation and ideal rotation frequency for our experiments. As expected, above some maximum rotation frequency the pebble becomes immobilized in the basin: centrifugal forces maintain the pebble on a given position against the wall. Also, under some minimum inclination angle, no fall of the pebble is observed, while a high inclination doesn't allow the pebble to reach it's maximum altitude. Altogether, we found it suitable to operate at a basin angle of 45 • and a rotation frequency of one cycle per second. Using the latter experimental conditions, we observed that the in-plane dimension of a pebble decreased by around a factor 2 after 30 minutes. Thus, a significant wearing of a pebble could be observed after a few minutes rotation. In practice, each pebble was eroded under the described conditions during 30 minutes, while a picture of the pebble was taken after each 5 minutes wearing. Hence, for each of the pebbles studied, we obtained about 7 pictures, corresponding to the initial pebble and to six following states of the wearing pebble. For some of the pebbles 8 or 9 pictures at 5 minutes interval were taken. The images were then analyzed following the method described in the appendix.

An example for the shape evolution produced by this method is given in Fig. 3, with photographs shown every five minutes. The corresponding cumulative curvature distributions f (K) are given in Fig. 4, where the inset shows the extracted contours. Here the initial shape is square, with four long nearly-flat regions and four short high-curvature regions. Thus the initial f (K) rises steeply around K = 0 and extends with relatively little weight out to K ≫ K . At first, the action of erosion is most rapid at the high-curvature corners, with the flat regions in between relatively unaffected. Thus the high-K tail of f (K) at first is suppressed, and weight builds up across (0.5 -2) K . Next the rounded corners erode further and gradually extend across the flat sections. Thus weight in f (K) is gradually concentrated more and more toward K . After about 15-20 minutes, when the flat sections are nearly gone, the form of f (K) fluctuates slightly but ceases to change in any systematic manner. In otherwords, the shape of the pebble has reached a final limiting form. Further erosion will affect pebble size, but not pebble shape!

To test the universality of the final shape, we repeat the same experiment both for other squares as well as for a variety of other initial shapes such as rectangles, triangles, and circles. A number of different examples showing both the initial and final shapes are shown in Fig. 5. In all cases, the cumulative curvature distribution f (K) shows a systematic evolution at short times and slight fluctuations about some average shape at later times, just as in Fig. 4. The more angular or oblong the initial shape, the more erosion is needed to reach a stationary final shape. The average final f (K) is shown for the various initial shapes in Fig. 6. Evidently, these all display the same quantitative form independent of the initial shape. Even f (K) for an initially-circular pebble broadens from a step function to the same form as all the

K), for the evolving pebble depicted in the inset (and pictured in Fig. 3). As the pebble becomes progressively rounder, the curvature distribution narrows and approaches a final shape. The time interval between successive contours is 5 minutes; the early contours of evolving shapes are shown as different blue dashes, while the later contours of stationary shape are shown as solid red. The same color and curve types are used in the main plot. For the inset, the axes are given in pixel units, equal to 0.132 mm. For contrast, a circle is shown by points.

others.

The final f (K) for all initial shapes can thus be averaged together for a more accurate description of the stationary shape produced by the laboratory erosion machine. The result is shown by the open circles in the same plot, Fig. 6. Differentiating, we obtain the actual curvature distribution, ρ(K), shown on the right axis. It is fairly broad, with a full-width at half-maximum equal to about 1.6 K . The actual shape is not quite symmetrical, skewed toward higher curvatures. The closest simple analytic form would be a Gaussian, exp[-(K -K ) 2 /(2σ 2 )]. The actual distribution is slightly skewed toward higher curvatures, but the best fit gives a standard deviation of σ = 0.70 K , as shown in Fig. 6. It is easy to imagine that the width of this distribution could be set by the strength of the erosion process. For example, if the angle of the rotating pan were lowered, then the erosion would be more gradual and more like polishing; in which case a rounder stationary shape may be attained with a narrower distribution of curvatures. The form of f (K), as well as its width, could also be affected. These types of questions can be addressed, both in laboratory and field studies, now that we have an incisive tool like f (K) for quantifying shape.

To further study the erosion produced by our laboratory apparatus, we now consider how the perimeter of the pebble decreases with time, P (t). Since the initial behavior depends on the specific initial shape, we focus on subsequent erosion once the (universal) stationary shape is achieved. If the final stationary shape of the curvature distribution is reached at time t 0 , then the quantity of interest is really P (t)/P (t 0 ) vs t -t 0 . The results, averaged over all laboratory pebbles, are shown in Fig. 7. Though the dynamic range is not great, the data are consistent with an exponential decrease, P (t) = P (t 0 ) exp[-(t -t 0 )/τ ]. The best fit to this form is shown by a solid curve; it gives a decay constant of τ = 44 min. Exponential erosion is, in fact, observed in field and laboratory studies [12]. It is to be expected whenever the strength of the erosion is proportional to the pebble size, as in our lab experiments where the impulse upon collision is proportional to the pebble's weight.

IV. FIELD STUDY

As a first field test of our method of analysis, we collected mud pebbles in the Mont St.-Michel bay, France. The littoral environment located at the inner part of the Norman-Breton Gulf is characterized by a macro-tidal dynamics. This location exhibits the second largest tide in the world after the Bay of Fundy, Canada. During the spring tide periods the upper part of the tidal flats collects a muddy sediment. This mud dries up during the following neap tide period where the sediments are exposed to the air. In certain areas, between the large equinoxial neap tides, the exposure of mud sediments to air may last for several months. During this period, this mud layer will develop a vast network of desiccation cracks. This network then leads to fragmented plates of a polygonal shape with 20 to 40 cm size. During the next spring tide period, the plates in the erosional area can be eroded by tide currents, thus re-incorporated into the sedimentary cycle. During the tide process these clasts are progressively eroded over many months. Thus the mud cohesion allows enough observation time for the life of a clast to be observed within a distance of order of one hectometer, from the original erosional area down to the latest stages of abrasion. Our approach is similar in spirit to field studies of river pebbles, where downstream distance serves as a surrogate for time. The mud clasts here have the advantage of remaining in a smaller area and of eroding on a human timescale.

We have analyzed the shapes of three classes of ripup clast photographed at three distinct locations on the tidal flats near Mont St.-Michel. The first is large subangular cobble found near the site of formation. Ten samples were examined; for these immature pebbles, the average perimeter is 550 mm. The second class is medium sub-mature pebbles found further "downstream". As a result of erosion, these pebbles are smaller and smoother than the cobble. Thirty five samples were examined; for these, the average perimeter is 180 mm. The third class is rounded, mature pebbles found on a nearby sand bar. The relation of these pebbles to the other two classes is not clear. Seventeen samples were examined; for these, the average perimeter is 220 mm.

Typical photographs for each of these classes are shown in Fig. 8. The average of the cumulative curvature distribution for all samples in each class is shown in Fig. 9. The angularity of the large cobble is reflected in the breadth of the curvature distribution. Roughly a quarter of the perimeter has negative curvature, and roughly a tenth has curvature five times greater than the average. For the other two classes, the curvature distribution is progressively more narrow. The relative steepness of the f (K)'s shows that all of these shapes are less round than the final pebbles produced by the laboratory erosion machine.

The width of the curvature distribution can be specified quantitatively by the standard deviation, σ. Results are normalized by the average curvature, and are shown for the field and laboratory pebbles in Table I. The pebbles with steeper f (K) indeed have smaller widths. For example, the width for the immature field cobble is about 3-4 times that of the average laboratory pebble. While the dimensionless width of the distribution, σ/ K , is a useful number for comparisons, it does not distinguish between curvature distributions of different shape. The actual functional form of the curvature distribution can be specified to some extent by comparing its moments with that of a Gaussian. In particular, the "skewness" and "kurtosis" are dimensionless numbers defined by the third and fourth moments, respectively, in such a way as to vanish for a perfect Gaussian. The results in Table I show that the four classes of pebbles have curvature distributions of four distinct forms. Of these, the laboratory pebbles are closest to a Gaussian. Note that the curvature distribution data show no evidence of a stationary shape, in which the clasts erode away without changing shape. Rather, the three classes of clasts all have different sizes and shapes. Thus it would not be fruitful to compare or fit the clast erosion to the cutting model we introduced in Ref. [START_REF] Durian | [END_REF]. However, since the sub-mature and mature clasts are smaller than the immature clasts but are not circular, we can rule out the polishing model.

Class

Perimeter (mm) σ/ K Skewness Kurtosis immature 550±100 2.4±0.3 1.4±0.3 3.8±1.4 sub-mature 180±60 2.0±0.6 0.0±0.9 3.3±2.9 mature 220±70 1.4±0.3 -0.2±0.5 1.3±1.0 lab-final 122±25 0.8±0.1 0.1±0.5 1.0±1.5 TABLE I: Characteristics of curvature distribution for the three classes of field pebbles. Final laboratory pebble shape is added for comparison.

V. CONCLUSION

We have studied the formation of two-dimensional pebble shapes. As in Ref. [START_REF] Durian | [END_REF] we introduced a local description of the erosion process, based on the distribution function of the curvature, measured along the pebble contours. This description captures both the local character of the erosion events, and the statistical nature of the erosion process.

For pebbles generated in the laboratory, we have shown that the curvature distribution has two important properties. First, the erosion drives the distribution towards a stationary form. When this stationary state is reached, the pebble contour still changes but, within small fluctuations, its curvature distribution remains the same, provided that the curvature is normalized by its average value. Secondly, we have found that the final stationary form of the distribution is independent of the original pebble shapes. This not only shows that the curvature distribution is a property of the erosion process itself, but it also opens the interesting possibility of establishing a classification of different erosion processes according to the type of curvature distribution they generate.

For pebbles collected in the field, we have made a first attempt to study a special class of rip-up clasts from the St.-Michel bay. These mud pebbles can be collected at very different erosion stages within a relatively small area of the tidal flats. We showed that the curvature distribution sharpens with the wearing degree, without getting however as sharp as the distribution obtained in the laboratory experiments.

The results presented in this experimental paper suggest a number of directions for modeling the formation of flat pebbles. Of central importance is the intrinsic statistical nature of the erosion process itself. As first hinted in Ref. [START_REF] Durian | [END_REF], a sequence of cuts of a noiseless, deterministic nature typically leads to a trivial curvature distribution like that of a circle. We also demonstrated in Ref. [START_REF] Durian | [END_REF] that a "cutting" simulation, with an appropriate distribution of cutting lengths, acting most strongly on regions of high curvature in accord with Aristotle's intuition [1], can reproduce the curvature distribution from the laboratory experiments. We will address these and other questions relevant for the theoretical modeling of pebble formation in a forthcoming paper.

APPENDIX: CURVATURE ANALYSIS

The goal of this section is to provide a detailed, practical description of two means to measure the local curvature at each point along the pebble contour. In both cases the starting point is an image of the pebble. For our work we use a digital camera Canon Power Shot G1 with a resolution of 1024 × 768 pixels. It should be just as effective to scan conventional photographs, or even to scan pebbles themselves. To determine the {x,y} coordinates of the contour, we import the images into NIH Image [START_REF]NIH Image[END_REF], which includes routines for finding edges and for skeletonizing the result. An example of the digitized pebble contour, and the smooth reconstructions to be discussed below, is given in Fig. 10. We show the contour in pixel units, since the actual calibration is needed only to determine size, not shape. Note from the inset that the contour points are indeed pixelized and skeletonized, with each point having only two neighbors located at either ±1 or 0 units away in the x-and y-directions. If the digitization process is faithful, then the uncertainty in each pixel is about ±0.5 units in each direction. This is not small compared to the distance between neighboring pixels, so a smoothing or fitting routine is necessary to reconstruct the actual contour and thereby extract the curvature distribution.

To illustrate the difficulties of extracting curvature, let us begin with two methods that are, in fact, unsatisfactory. It is perhaps tempting to simply smooth the data, replacing each point with a weighted average of neighbors lying within some window. Weights could be cleverly chosen to de-emphasize points at the edge of the window, for example. This fails, however, since it's far from obvious how to choose a suitable window size. For instance, the pixelized representation of the straight section given by y = 0.1x for 0 < x < 10 is a step function y = 0(1) for x < (>)5. A large window would be needed to even approximately reconstruct the original line. However, such a sufficiently large smoothing window would erase fine features if applied elsewhere along the contour. Since smoothing filters provide no feedback on quality, visual inspection of the result would be necessary to choose an optimal window size at each point along the contour. This is not only subjective, but rather impractical. As an alternative, it is perhaps tempting to implement an automated version of Wentworth's curvature gauge [9]. This is a device with circular notches of various diameters into which portions of a pebble may be pressed. The computational analogue would be to find the best nonlinear-least-squares fit to a circle at each point along the pebble contour. As with smoothing, one difficulty is to find the optimal window over which to do the fitting. A compounding difficulty is that essentially nowhere is the pebble exactly circular, so even with 'the' optimal window there is substantial disagreement between data and fit. A spectacular example of this problem is at an inflection point, where the curvature changes from positive to negative.

To overcome such difficulties, we propose to fit digitized contour data to a cubic polynomial at each point along the contour. A cubic is the lowest order polynomial needed in order to avoid systematic error when the curvature varies gradually across the fitting window, which is the usual case. In order to avoid having to rotate the coordinate system to ensure that the contour y(x) is a single-valued function, we instead convert to polar coordinates. Thus we define the origin by the centerof-mass of the contour and perform fits to r(θ) where r = x 2 + y 2 and θ = tan -1 (y/x). Once a satisfactory fit is achieved, the curvature may be deduced from the value and derivatives of the cubic polynomial by

K = r 2 + 2r 2 θ -rr θθ (r 2 + r 2 θ ) 3/2 , (A.1)
where r θ = dr/dθ and r θθ = d 2 r/dθ 2 [START_REF] Weisstein | The CRC concise encyclopedia of mathematics[END_REF]. Two tricks seem necessary to achieve satisfactory results. The first is to weight the data most heavily near the center of the window. We use a Gaussian weighting function with a standard deviation equal to 1/4 of the width of the window. This ensures that points at the edges have essentially no influence. Therefore, the fitting results do not vary rapidly as the window is slid along the contour. This guarantees that the reconstructed curve and its first two derivatives are continuous, which is a crucial requirement for measuring the curvature.

The second trick is to choose the window size appropriately. This is actually the most difficult and subtle aspect of the whole problem. If the window is too small, then the fit will reproduce the bumps and wiggles of the pixelization process; usually the curvature will be overestimated. If the window is too large, then the fit will significantly deviate from the data; usually the curvature will be underestimated. And while the curvature tends to decrease systematically with window size, there is in general, unfortunately, no plateau between these two extremes where the curvature is relatively independent of window size and hence clearly represents the true value. To pick the window appropriately requires careful understanding of the numerical fitting procedure and the feedback it provides. Since the fitting function is a polynomial, the minimization of the χ 2 total square deviation from the data reduces to solving a set of linear equations. This in turn reduces to inverting a matrix. If the window is too small then the fit will be 'ambiguous' in the sense that χ 2 is small but the error in the fitting parameters is large. Mathematically, the matrix to be inverted is essentially singular. A good strategy is therefore to start with a small window and increase its size until the matrix is no longer singular. This can be accomplished using a linear least-squares fitting routine based on singular-value decomposition [START_REF] Press | Numerical Recipes in C[END_REF]. However, the uncertainty in fitting parameters for the first suitable window is generally too large (nearly 100%). So we increase the window size two pixels at a time until the error in curvature has been reduced by a factor of ten. This defines the largest suitable window, beyond which systematic errors due to incorrect functional form begin to appear. We have not been able to define the largest suitable window based on the value of χ 2 . For the final result, we take a weighted average of the cubic fit parameters over all suitable windows, where the weights are set by the uncertainties in fitting parameters as returned by the fitting routine. An example of a reconstructed contour from this procedure is shown by the solid curve in Fig. 10. The recontruction is satisfactorily smooth; also, it clearly avoids the pixel noise without smoothing over significant small-scale features in the contour.

Since the cubic fitting method is rather involved, and since the choice of window sizes is still slightly subjective, we have developed an alternative method. The starting point is the fact that the actual digital representation of the contour depends on the location and orientation of the pebble with respect to the grid of pixels. If the pebble were shifted or rotated, then the pixelized representation would be slightly different. For example, imagine the pixelization of a line making various angles with the grid. Perhaps the ideal experimental measurement procedure would be to systematically reposition the pebble, pixelize, then compute the average of all such representations. However this procedure does not lend itself to automatization, and would be impractically timeconsuming. Instead, we propose to do more or less the same thing numerically. The idea is to take the current best guess for the contour, pixelize it with respect to a random grid position and orientation, then use the new representation to update the best guess. When iterated, this procedure converges to a satisfactory reconstruction of the actual contour with two provisos. First, at each step, we locally smooth the trial pixelization by replacing each point by its average with its two immediate neighbors. Second, we keep only every fourth or fifth point in the original pixelized data and perform all operations on this subset. When done, we compute the curvature literally by the change of slope with respect to arclength using the straight segments between adjacent points.

The cumulative curvature distribution given by this iteration scheme is shown in Fig. 11 as a function of the number of points kept. When too many are kept, the reconstructed curve follows the bumps and wiggles of the original pixelization too closely; the curvature distribution is too broad. When too few are kept, the reconstructed curve incorrectly smooths over small-scale fea- tures; the curvature distribution is too narrow. When only every fourth or fifth point are kept, the distributions are nearly equal; furthermore, they are indistinguishable from that given by the cubic polynomial fitting. The actual reconstructed contour is also shown in Fig. 10. The plateau in the curvature distribution vs number of points kept, and the good agreement with the other method, both give confidence in this new iterative reconstruction scheme.

[1] Aristotle, Minor Works, Mechanical Problems (Harvard University Press, Cambridge, 2000), translation by W.S. Hett. Question 15 asks: "Why are the stones on the seashore which are called pebbles round, when they are originally made from long stones and shells? Surely it is because in movement what is further from the middle moves more rapidly. For the middle is the center, and the distance from this is the radius. And from an equal movement the greater radius describes a greater circle. But that which travels a greater distance in an equal time describes a greater circle. Things travelling with a greater velocity over a greater distance strike harder, and things which strike harder are themselves struck harder. So that the parts further from the middle must always get worn down. As this happens to them they become round. In the case of pebbles, owing to the movement of the sea and the fact that they are moving with the sea, they are perpetually in motion and are liable to friction as they roll. But this must occur most of all at their extremities.".

  FIG.1:(Color online) (a) Radius, (b) normalized curvature and (c) fraction of the perimeter f (K) with curvature less than K, vs K divided by the average curvature K , for a superellipse, oval, ellipse, and circle. The curve types match those for the shapes as shown in the inset. Note that, except for the circle, all have the same aspect ratio a = 3/2. The differences in shape are reflected in differences in the forms of f (K).

FIG. 2 :

 2 FIG. 2:The wearing apparatus used for the laboratory experiments. The rotating metal tray is 30 × 30 × 7 cm 3 .

FIG. 3 :

 3 FIG. 3: Shape evolution of a 5 × 5 cm square pebble eroded in the laboratory, by the method explained in the text. The images were taken at 5 min intervals, and are shown at the same magnification.

  FIG.4: (Color online) Cumulative curvature distribution, f (K), for the evolving pebble depicted in the inset (and pictured in Fig.3). As the pebble becomes progressively rounder, the curvature distribution narrows and approaches a final shape. The time interval between successive contours is 5 minutes; the early contours of evolving shapes are shown as different blue dashes, while the later contours of stationary shape are shown as solid red. The same color and curve types are used in the main plot. For the inset, the axes are given in pixel units, equal to 0.132 mm. For contrast, a circle is shown by points.

FIG. 5 :

 5 FIG. 5: The initial and final forms of different shapes eroded in our experiment, all shown with the same magnification. The erosion times for the final shapes are 40 min for square, 30 min for rectangle, 30 min for triangle, 35 min for polygon, 35 min for circle, and 30 min for lozenge. Note that the final shapes are all roughly circular.

FIG. 6 :FIG. 7 :

 67 FIG. 6: (Color online) Integrated curvature distribution, f (K), for the final shapes of laboratory pebbles with various initial shapes, as labeled. These are indistinguishable to within measurement accuracy; their average is shown by the open circles. The corresponding average curvature distribution, ρ(K) is obtained by differentiation and is shown on the right axis along with a fit to a Gaussian shape.

FIG. 8 :

 8 FIG. 8: Typical shapes of Mont St.-Michel rip-up clasts. The immature pebbles in the top row were collected close to their origin; the sub-mature pebbles in the middle row were collected further downstream; the smooth, mature pebbles in the bottom row were collected on a nearby sandbar. As these pebbles eroded, their shapes became rounder, an effect quantified in the next figure. The bars indicate, in each row, a length of 2 cm.

FIG. 9 :

 9 FIG. 9: (Color online) Cumulative curvature distribution, f (K), for the average shapes of Mont St.-Michel rip-up clasts.Even the roundest shapes remained less circular than the final shape in the laboratory study.

FIG. 10 :

 10 FIG. 10: (Color online) Reconstruction of a smooth pebble contour from the pixelized digital representation. The solid curve is based on fitting to cubic polynomials at each point. The dashed curve is based on the iteration scheme of Fig. 11.

FIG. 11 :

 11 FIG.11:(Color online) Trial integrated curvature distributions from the iterative smoothing scheme, vs the number of points kept. Keeping either 4 or 5 points seems optimal: the resulting distributions are identical and they agree with that based on cubic fits.

  [2] Lord Rayleigh, F.R.S., Nature 3901, 169 (1944). [3] S. Boggs, Jr., Principles of sedimentology and stratigraphy (Prentice Hall, N.J., 2001), 3rd ed. [4] P. M. Chaikin and T. C. Lubensky, Principles of Con-

Acknowledgments

We acknowledge insightful discussions with F. Thalmann and experimental insight by P. Boltenhagen. This work was supported by the Chemistry Department of the CNRS, under AIP "Soutien aux Jeunes Equipes" (CM). It was also supported by the National Science Foundation under Grant DMR-0514705 (DJD).


Dora Boviz
email: dora.boviz@nokia-bell-labs.com

Chung Shue Chen
email: chungshue.chen@nokia-bell-labs.com

Sheng Yang
email: sheng.yang@centralesupelec.fr

Cost-Aware Fronthaul Rate Allocation to Maximize Benefit of Multi-User Reception in C-RAN

Aiming to throughput enhancement in future mobile networks, dense deployment of wireless access points (APs) is intended. Spectral efficiency on the uplink can be improved via joint processing of users located in areas covered by several APs. Thanks to centralized processing of such users enabled by Cloud Radio Access Network (C-RAN) architecture, low-latency multi-cell cooperation becomes possible. However, the price to pay to benefit from its advantages is the cost of transferring data from distributed access points to the data center through limited-capacity fronthaul links. A mobile network operator using C-RAN would have the objective to maximize data transmission rates with lowest possible fronthaul rate. We consider in this work the optimal tradeoff between the amount of fronthaul allocated and the sum rate of each user group in a multiple access scheme, where the wireless resource is reused among the user groups. By performing low-complexity joint reception of signals quantized according to attributed fronthaul rate, we can maximize the benefit of uplink transmission in C-RAN despite constrained fronthaul links. With optimal fronthaul allocation, the net benefit of the uplink transmission improves by 10% in a usual configuration.

I. INTRODUCTION

The idea of offloading base-band functions from several cell sites to a centralized data center arises the great opportunity of joint signal processing, but also the practical challenge of data transfer between remote Access Points (APs) and the Central Office (CO). Furthermore, in 5G architecture, centralization is identified as a key disruption [START_REF] Boccardi | Five disruptive technology directions for 5G[END_REF], since it enables user-centric processing which can greatly facilitate multi-cell cooperation. Finding an efficient and affordable data transfer solution over limited capacity fronthaul (FH) links between APs and the CO is crucial to realize Cloud Radio Access Network (C-RAN) deployments and fulfill 5G requirements [START_REF]Nokia, "5G use cases and requirements[END_REF].

Mobile network operators have to offer better Qualityof-Service (QoS) at lower cost, thus the efficient use of wireless spectrum, APs, and computational resources become more important then ever. In this paper, we will propose an optimized transmission scheme for C-RAN to increase uplink (UL) transmission efficiency from the User Equipments (UEs) to the end of physical layer (PHY) processing. To enhance spectral efficiency, low latency multi-cell cooperation is realized through partially non-orthogonal uplink scheduling. We aim to benefit from Multiple Input Multiple Output (MIMO) enabled APs and centralization of signal processing while keeping user coordination and receiver complexity low. To deal with limited capacity fronthaul links between the APs and the CO, it is important to optimize the fronthaul rate allocation in order to adapt quantization level of transferred signals according to the quality of the wireless transmission.

Constrained Non-Orthogonal Multiple-Access (NOMA) is proposed in [START_REF] Al-Imari | Uplink nonorthogonal multiple access for 5G wireless networks[END_REF] for uplinks where users are organized into several subsets and within each of them, the users transmit simultaneously over the same wireless channel, whereas orthogonal resources are allocated to each group in order to eliminate the interference between the groups. On the receiver side, one can use Minimum Mean Square Error (MMSE) detection to separate signals of users transmitting together. Thanks to the limited number of users in each group, the receiver complexity can be relatively low compared to that needed with fully nonorthogonal transmission of all users. Maximizing the sum rate of UL multi-user transmission through enhanced user association is described in [START_REF] Rui | Robust user pairing algorithm under channel estimation errors for uplink virtual multiple-input multiple-output systems[END_REF], where the authors consider channel estimation error and compare different methods for sum rate improvement. Joint fronthaul and power allocation in fully orthogonal multi-user model is studied in [START_REF] Liu | Joint power control and fronthaul rate allocation for throughput maximization in OFDMA-Based Cloud Radio Access Network[END_REF]. The authors provide practical quantization scheme and show that the optimization of fronthaul allocation improves significantly the overall throughput compared to uniform distribution of available fronthaul. The improvement is obtained by the novel scheme which takes into account different link gains of received signals that are forwarded on the same fronthaul link. Impact of signal compression in C-RAN architecture to satisfy fronthaul constraint is widely studied, see for example [START_REF] Kang | Joint signal and channel state information compression for uplink network MIMO systems[END_REF]. In our previous work [START_REF] Boviz | Optimal fronthaul capacity allocation and training for joint uplink receiver in C-RAN[END_REF], we have characterized the optimal fronthaul rate for a single MIMO channel. In this paper, we extend the optimization to a system of several user groups that transmit independently and share the fronhaul links. We show that optimally distributing available fronthaul rate between multi-user transmissions allows to improve the overall efficiency of the transmission. In an typical scenario, end-to-end benefit increases by 10% thanks to the proposed tradeoff between sum rate and fronthaul usage.

The rest of the paper is organized as follows. We describe the system model and the multi-user transmission scheme in Section II. Then, in Section III we develop the metric linking system performance to the fronthaul rate. The optimization problem of fronthaul allocation is formulated in Section IV. Section V provides numerical evaluation of the proposed fronthaul allocation scheme and observations. Finally, Section VI contains the conclusion.

II. SYSTEM MODEL

The system that we study is depicted in Figure 1. We consider a C-RAN system with M RRHs located at the cellsites with a antennas each, thus the total number of antennas is K = a • M . All of the RRHs are connected to the same CO where the major part of signal processing is realized for the cells. We assume that all the N users are located in the cell-edge region and are able to communicate with all the M RRHs. To avoid interference between them, they are benefiting from joint detection (JD) at the receiver. In the system that we study, user-PHY functions of all the N users are placed in the CO in order to enable multi-cell joint processing on the uplink [START_REF] Boviz | Fronthaul for 5G: low bit-rate design enabling joint transmission and reception[END_REF].

The users can be associated to each other randomly or in a deterministic way to create the user groups to which we allocate a wireless resource. Orthogonal resource blocks (RBs) are allocated to different groups, but users who transmit using the same RB will interfere with each other. The RRHs, after receiving the combination of uplink signals transmitted by every group, will forward the quantized version of this signal to the CO through limited-capacity fronthaul links. The total capacity of each fronthaul link is shared among all of the user groups.

We assume a block fading channel with a coherence time that is long enough to consider allocating fronthaul capacity to the groups using real-time channel realizations. In terms of user mobility, e.g., for pedestrian users moving at 5 km/h, the transmission characteristics change approximately every 720 ms, which is a time window large enough to update network parameters.

We use the following notational conventions: for random variables, upper case letters with bold and non-italic fonts, e.g., V, for vectors, and bold and sans serif fonts, e.g., M, for matrices. Deterministic quantities are denoted with italic letters, e.g., a scalar x, lowercase bold for a vector v v v, and uppercase bold letters for a matrix M M M . Logarithms are in base 2 and superscript (.) H denotes the conjugate transpose of a vector or a matrix.

We have N users uniformly distributed in the region covered by every RRH. The channel of each user n ∈ {1, ..., N } towards all antennas at both RRHs is denoted by h h h n which is a K dimensional array following the Gaussian distribution

N (0, R R R n ) with R R R n ∈ C K×K .
The number of user groups is denoted L, there are s l users in the group l ∈ {1, ..., L}. The multi-user channel of a group towards the whole set of antennas is the K × s l dimensional matrix H l . We assume that the channel is perfectly known at the receiver. Though, channel estimators that are generally used result in an estimation error, its impact on the fronthaul allocation is not significant.

The Gaussian channel noise vector is denoted n n n l ∼ N (0, σ 2 z I I I K ). Power of the input signal is normalized, so that noise covariance scales with σ 2 z = 1 SNR . The received signal for group l by the whole set of antennas is the K-dimensional vector y y y l and it is given by the superposition of the signals sent by all of the s l users in the group denoted by Π l = {π l 1 , ..., π l s l }.

y y y l = s l i=1 h h h π l i x π l i + n n n l (1)
The fronthaul capacity of the link between the m-th RRH and the CO allocated to group l is denoted by c (l) m with l ∈ {1, ..., L} and m ∈ {1, ..., M }. We can write this total fronthaul rate available using fronthaul rates c lk dedicated to the transmission of the signal from group l received by antenna k located at RRH m as

c (l) m = m•a k=(m-1)•a+1 c lk . (2)
We denote by c c c (l) = {c l1 , ...c lK } the set of capacity values attributed to group l and c m the total fronthaul capacity available between RRH m and the CO.

III. SUM RATE OF MULTI-USER TRANSMISSION WITH LIMITED FRONTHAUL

Once user groups are created and scheduled, users transmit following the scheduling decisions and RRHs receive the signals of every group -being the superposition of the signals transmitted by all the UEs in the group -and forward frequency-domain I/Q symbols to the CO. We perform receive Fast Fourier Transform in the RRH to decorrelate the subcarriers and enable to select the ones that need to be forwarded (since they are required for the JD). This decorrelation improves also the throughput by reducing quantization noise in the forwarded signal.

The channel gain depends on subcarrier frequency, so the fronthaul allocation can change with a different scheduling decision or a modification of user groups. In the following computations, we consider only one subcarrier per group, but the extension to several subcarriers is straightforward.

A. Quantization noise

We define in the following the equivalent quantization noise with limited capacity fronthaul for a group with s l users transmitting towards K receive antennas. The achievable rate with a fronthaul capacity c lk can be defined using the mutual information between the received signal and the forwarded one, given the distortion between them.

r lk ≥ min p ŷlk |y lk :D≤σ 2 d lk I(Y lk ; Ŷlk | H l ). (3)
We can derive the following relation between the received signal power and the variance of the distortion noise

σ 2 d lk ≤ σ 2 y lk |H l 2 -c lk (4)
where

σ 2 y lk |H l = s l i=1 (| h kπ l i | 2 ) + σ 2 z . (5)
We use a scaling factor α lk in order to adapt the power of forwarded signal to the fronthaul capacity used, i.e.,

α lk = σ 2 y lk |H l -σ 2 d lk σ 2 y lk |H l ∀k ∈ {1, ..., K}. (6)
Scaling factors for each antenna form the matrix A A A l = diag k={1,...,K} (α lk ). The distortion is then characterized by the following upper bound:

σ 2 d lk α lk = σ 2 d lk σ 2 y lk |H l σ 2 y lk |H l -σ 2 d lk ≤ σ 2 y lk |H l 2 -c lk 1 -2 -c lk , ∀k ∈ {1, ..., K}, ∀l ∈ {1, ..., L}. (7)

B. Achievable sum rate

We can compute the achievable sum rate for a given joint transmission group using the mutual information between the signal sent by all users in the group and the one received in the CO:

s l i=1 r i ≥ I(X l ; Ŷl | H l ) = h(X l ) -h(X l | Ŷl , H l ) (8)
We compute both entropy terms in order to find the lower bound of the sum rate. The first term describes the quantity of information sent by the users, and thus depends on the transmission power:

h(X l ) = log(det(2πeE[X l X H l ])) = log((2πe)) (9)
The second term comes from the loss of information between the UEs and the CO and can be computed using the linear MMSE covariance C e .

h(X l | Ŷl , H l ) ≤ log(det(2πeC e )). (10)
For this, we use the definition of received signal by the CO:

ŷ y y l = A A A l s l i=1 h h h π l i x x x i + n n n l (11)
where n n n = z z z + d d d is the equivalent noise containing Gaussian channel noise and the quantization noise. The covariance matrix of this equivalent noise is

C N = σ 2 z I K + diag k={1,...,K} σ 2 d lk α lk .
We can compute the linear MMSE covariance based on the channel matrix

C e = I I I s l -H H H H l (H H H l H H H H l + C N ) -1 H H H l . ( 12
)
Using the inversion lemma on the lower bound of the mutual information we get the following expressing of the achievable sum rate with a givel channel matrix

H H H l s l i=1 r i ≥ log det I s l + H H H H l C -1 N H H H l . (13)
From ( 7) and (13), we get the achievable sum rate of a group of N users transmitting towards K receive antennas, each of them using c j , j ∈ {1, ..K} bits of fronthaul.

s l i=1 r i ≥ log det I s l + H H H H l V V V -1 s l H H H l ( 14
)
with the equivalent noise covariance:

V V V s l = σ 2 z I I I K + diag k={1,...,K} σ 2 y lk |H l 2 -c lk 1 -2 -c lk . (15)
IV. FRONTHAUL ALLOCATION OPTIMIZATION After the CO receives forwarded signals from the RRHs, we can compute the sum rate (14) of every user group using channel realizations in order to evaluate the performance of the multi-user reception in C-RAN under the constraint of limited fronthaul capacity. Using the sum rate expression (14), our aim is to maximize the benefit of the transmission given a constraint on total rate available at each fronthaul link. We provide here a metric enabling the allocation of fronthaul capacity to user groups, which should allow to maximize the efficiency of the transmission. When the optimal fronthaul allocation for each fronthaul link is computed, quantization of forwarded signal is adapted at the RRHs, in order to improve the performance of the complete transmission.

A. Objective function: net benefit of the transmission

We use the upper bound of the sum rate (14) to formulate the objective function allowing to maximize the end-to-end benefit of the uplink transmission of N users forming L groups towards the M RRHs with a antennas each. The parameters of this function are the following:

• The Gaussian channel noise variance σ 2 z . • The average received signal power from group l at antenna k given the channel of the group l: σ 2 y lk |H l . • The fronthaul capacity c lk with k ∈ {1, ..., K} used to forward to the CO digital base-band I/Q symbols of group l received by antenna k. We denote c c c (l) = (c l1 , ...c lK ) T the set of capacity values for group l. • The cost λ k of the fronthaul capacity used for the transmission of the signal received on antenna k.

In this work we consider only linear fronthaul cost corresponding to exploitation costs such as transmission energy or cost renting the needed transmission rate from the owner of the fronthaul infrastructure. Note that any non-negative convex cost function can be used following the deployment scenario considered. The following function characterizes the net benefit of the transmission of group l towards the whole set of receive antennas when the fronthaul capacity allocated to the group is

c (l) m = m•a k=(m-1)•a+1
c lk for the fronthaul link between RRH m and the CO. Given the parameters σ 2 z , σ 2

y lk |H l , λ k ∀k ∈ {1, ..., K}, f (H H H l , c c c (l) ) = log det I s l + H H H H l V V V -1 s l H H H l - K k=1 λ k c lk (16)
where V V V s l is defined in (15). The first term of the function f (.) in ( 16) gives the instantaneous sum rate during the coherence time block where the channel matrix H H H l holds and the second term is the total cost of the fronthaul transmission for group l over the whole set of fronthaul links connecting the RRHs to the CO.

B. Limited fronthaul link rate

First, we aim to find the optimal capacity allocation {c c c (1) * , ..., c c c (L) * } that maximizes the overall benefit of the uplink transmission of the N users distributed in L groups while the fronthaul capacity used for all groups between each RRH m and the CO is not more than c m . We need to solve the optimization problem including this constraint when fronthaul links are physically limited to a given rate, for example in low capacity deployments using other transport solution than optical fibers or when the allocation of rates between various services exploiting the same transport links is fixed. 17) is concave, thus admits a unique solution that gives the optimal capacity allocation scheme.

Proof: The constraint in ( 17) is linear, as well as the cost term in (16), thus the concavity of the first term of f (.) is sufficient to show that the problem is concave. Notice that the second term of (16), i.e., K k=1 λ k c lk is linear w.r.t. c c c (l) , so it can be considered as concave. The function log det(A A A) is concave if and only if the matrix A A A is non-negative definite. The sum of two non-negative definite matrices is also non-negative definite. Since the identity matrix satisfies this condition, we only need to show that the second term of the argument of the log det(.) in ( 16) is non-negative definite. The equivalent SNR matrix V V V s l is diagonal with positive elements which are its eigenvalues, thus it is positive definite. This property stands also for its inverse.

If a positive definite matrix M M M is multiplied by another matrix and its hermitian as

B B B H M M MB B B, the result is also positive definite if B B B is full rank. This is true for H H H H l V V V -1 s l H H H l since the columns of H H H are independent, thus rank(H H H) = s l .
Consequently, the matrix being the argument of log det(.) is positive definite and also non-negative definite, thus the first term of (16) which implies with the above reasons that (17) is concave.

C. Unconstrained fronthaul links

In some deployment scenarios, very high capacity fronthaul links are available between RRHs and the CO. These can be considered in practice as unlimited links. For example, when optical fibers with Dense Wavelength Division Multiplexing are used, they can satisfy bandwidth requirements of mobile fronthaul, but installation and consequently usage cost may be high. Unconstrained optimization of the benefit of uplink transmission would allow to allocate fronthaul optimally for a given fronthaul cost.

Without constraint on the maximal fronthaul rate available at each link, to solve the optimization problem we need to maximize the same concave function as in (17). Intervals of possible fronthaul capacity values are specified according to the transmission scheme, see below (19). However, they do not constrain the problem in practice.

A. System parameters

We have evaluated the results of the above fronthaul allocation optimization (see equations ( 17) and ( 19)) with N = 40 users located at the edge of 2 cells with 1 RRH each. Note that the optimization problem can be solved efficiently using standard convex programming [START_REF] Boyd | Convex Optimization[END_REF]. We have K = 8 antennas equally distributed between the RRHs. Channel realizations are modeled using independent one-ring scatterer model for each user [START_REF] Shiu | Fading correlation and its effect on the capacity of multielement antenna systems[END_REF]. Users are transmitting jointly on the uplink by groups of 4 users each, each user group having its own RB.

B. Efficient transmission with constrained fronthaul

In a scenario where we aim to find the optimal fronthaul allocation scheme in presence of limited fronthaul links, one can be interested to evaluate if the available amount of fronthaul capacity allows an efficient transmission. Furthermore, while designing fronthaul infrastructures, the evaluation of required capacity to get optimal benefit from the transmission with expected price values (λ k ) can contribute to correct the dimensioning of the links to be deployed.

In Figure 2, we compare the proposed fronthaul allocation scheme that optimizes the net benefit of the transmission against uniform fronthaul allocation for different values of available fronthaul capacity. In case of uniform fronthaul allocation, available fronthaul capacity is equally distributed between all groups and all antennas. When available fronthaul is low, both uniform and optimized allocation result in similar efficiency, since the constraint does not allow to achieve higher sum rate. One can notice that optimized fronthaul allocation allows to achieve higher transmission benefit, since in case of sufficient fronthaul, group sum rates can be improved by allocating more fronthaul to the received signals with higher powers. In other words fronthaul allocation is adapted to the variations of channel gains for different users and antennas. Even in an ideal scenario with λ k = 0.1 and a relatively high available fronthaul rate, e.g., 300 bits/channel use, the proposed allocation scheme increases net benefit of the transmission by approximately 10%. Furthermore, by optimizing the fronthaul allocation, in the operation regime where the sum rate is limited by the transmission power, allocated fronthaul does not continue to increase so that the cost remains moderate. In uniform allocation cost continues to increase while the sum rate is limited, which result in the drop of the benefit that we can observe in Figure 2. Obviously, if the price of the fronthaul usage is higher, the transmission is less efficient for both allocation strategies, since the maximal sum rate does not increase while the cost does.

C. At what price is it still useful to transmit?

When fronthaul capacity has a fixed limitation (including practically infinite capacity available), the net benefit of the transmission can turn to negative values when uniform allocation scheme is used and cost is too high. By optimizing capacity allocation, we can avoid having higher fronthaul cost than the utility of the transmission, i.e., the sum rate. However, reducing allocated fronthaul rate when the price factor is high would result in a decrease of the sum rate. Depending on the transmission requirements, setting a lower bound on the sum rate that allows to achieve the required quality-of-service can be used to select whether it is useful to transmit with a given price or not.

We can see in Figure 3 that the more expensive the fronthaul is, the more we improve transmission benefit by fronthaul optimization since we can avoid allocating very high cost fronthaul rate. Consequently, optimized fronthaul allocation can reduce the sum rate of the transmission, as we can observe in Figure 4. If the sum rate is very low because of low FH rate that should be allocated, it can be more reasonable to not to transmit. One can also note that the optimal sum rate and net benefit values are almost equal for different fronthaul constraints. A gap can be observed only for low pricing and low available fronthaul rate (see the curve FH limit = 240), since we cannot achieve high benefit despite the optimal fronthaul allocation. We compare the amounts of allocated fronthaul rate per link with different constraints (including infinite) and varying fronthaul price in order to evaluate the difference between constrained and unconstrained optimization strategies. We can observe in Figure 5 that if available fronthaul rate is high enough (see the curve FH limit = 400) to enable to achieve maximal sum rate, the amount of optimally allocated fronthaul is almost the same with and without the constraint on available fronthaul for any value of the price factor. However, with less fronthaul rate available, obviously less of it is allocated per link even with a low price. For higher price values, we get similar results regardless of the amount of available FH rate. The reason is that under low fronthaul price, the sum rate term dominates the optimization, so that it is highly motivated to increase the sum rate in the optimization, however it will be limited by the total amount of available fronthaul. For higher price factor, fronthaul cost must be kept low in order to maximize the benefit of the transmission, so the allocated fronthaul will become similar when the price increases with and without constraint.

VI. CONCLUSION AND FUTURE DIRECTIONS

A main limitation of C-RAN architecture is the capacity limited and expensive fronthaul links. In this paper, we have analyzed a low complexity multi-user multi-cell transmission scheme where fronthaul rate allocation is performed to maximize uplink throughput. We provide a novel system model where constrained NOMA is used for several cells and multi-antenna RRHs are connected to a CO where UL JD is performed. To optimize the gain with respect to the limited fronthaul capacity and also price, we allocate the FH rate based on a performance metric exploiting uplink channel gains. We evaluate the net benefit of the transmission which allows a tradeoff between uplink sum rate including practical parameters and the cost of fronthaul usage.

A future work is to explore possible further improvement of the system performance by associating in each group users that can achieve higher sum rate together, i.e., user grouping optimization. This would be challenging due to the increased system complexity and also the fact that only limited channel state information is available before the radio transmission scheduling. However, we can expect further performance gains by adding the optimization of user association.

Fig. 1 :

 1 Fig. 1: System model with several user groups (example M = 2 RRHs), where JD stands for joint detection and BBU stands for base-band unit.

FindcProposition 1 .

 1 {c c c (1) * , ..., c c c (L) * } = argmax {c c c (1) ,...,c c c (L) } L l=1 f (H H H l , c c c (l) ) subject to L l=1 c (l) m ≤ c m ∀m ∈ {1, ..., M }. lk ≤ c m ∀m ∈ {1, ..., M } (18) The problem (

Find

  {c c c (1) * , ..., c c c (L) * } = argmax {c c c (1) ,...,c c c (L) } L l=1 f (H H H l , c c c (l) ) with 0 < c * lk < σ 2 y lk |H l ∀k ∈ {1, ..., K}, ∀l ∈ {1, ..., L}.

2 Fig. 2 :

 22 Fig. 2: Net benefit of uplink transmission with constrained fronthaul

Fig. 3 :

 3 Fig. 3: Net benefit of uplink transmission with various fronthaul constraints vs. fronthaul price factor

Fig. 4 :Fig. 5 :

 45 Fig. 4: sum rate of uplink transmission with various fronthaul constraints vs. fronthaul price factor