Robert Geist

Zachary H Jones

James Westall

Predicting Disk Scheduling Performance with Virtual Machines

A method for predicting the performance of disk scheduling algorithms on real machines using only their performance on virtual machines is suggested. The method uses a dynamically loaded kernel intercept probe (iprobe) to adjust low-level virtual device timing to match that of a simple model derived from the real device. An example is provided in which the performance of a newly proposed disk scheduling algorithm is compared with that of standard Linux algorithms. The advantage of the proposed method is that reasonable performance predictions may be made without dedicated measurement platforms and with only relatively limited knowledge of the performance characteristics of the targeted devices.

Introduction

In the last five years, the use of virtual computing systems has grown rapidly, from nearly non-existent to commonplace. Nevertheless, system virtualization has been of interest to the computing community since at least the mid-1960s, when IBM developed the CP/CMS (Control Program/Conversational Monitor System or Cambridge Monitor System) for the IBM 360/67 [START_REF] Creasy | The origin of the VM/370 time-sharing system[END_REF]. In this design, a low-level software system called a hypervisor or virtual machine monitor sits between the hardware and multiple guest operating systems, each of which runs unmodified. The hypervisor handles scheduling and memory management. Privileged instructions, those that trap if executed in user mode, are simulated by the hypervisor's trap handlers when executed by a guest OS.

The specific architecture of the host machine essentially determines the difficulty of constructing such a hypervisor. Popek and Goldberg [START_REF] Popek | Formal requirements for virtualizable third generation architectures[END_REF] characterize as sensitive those machine instructions that may modify or read resource configuration data. They show that an architecture is most readily virtualized if the sensitive instructions are a subset of the privileged instructions.

The principal roadblock to widespread deployment of virtual systems has been the basic Intel x86 architecture, the de facto standard, in which a relatively large collection of instructions are sensitive but not privileged. Thus a guest OS running at privilege level 3 may execute one of them without generating a trap that would allow the hypervisor to virtualize the effect of the instruction. A detailed analysis of the challenges to virtualization presented by these instructions is given by Robin and Irvine [START_REF] Robin | Analysis of the Intel R Pentium's TM ability to support a secure virtual machine monitor[END_REF].

The designers of VMWare provided the first solution to this problem by using a binary translation of guest OS code [START_REF] Vmware | Understanding full virtualization, paravirtualization, and hardware assist[END_REF]. Xen [START_REF] Barham | Xen and the art of virtualization[END_REF] provided an open-source virtualization of the x86 using para-virtualization, in which the hypervisor provides a virtual machine interface that is similar to the hardware interface but avoids the instructions whose virtualization would be problematic. Each guest OS must then be modified to run on the virtual machine interface.

The true catalysts to the widespread development and deployment of virtual machines appeared in 2005-2006 as extensions to the basic x86 architecture, the Intel VT-x and AMD-V. The extensions include a "guest" operating mode, which carries all the privilege levels of the normal operating mode, except that system software can request that certain instructions be trapped. The hardware state switch to/from guest mode includes control registers, segment registers, and instruction pointer. Exit from guest mode includes the cause of the exit [START_REF] Kivity | kvm: the Linux virtual machine monitor[END_REF]. These extensions have allowed the development of a full virtualization Xen, in which the guest operating systems can run unmodified, and the Kernel-based Virtual Machine (KVM), which uses a standard Linux kernel as hypervisor.

Nevertheless, VMWare, Xen, and KVM are principally aimed at facilitating user-level applications, and thus they export only a fairly generic view of the guest operating system(s), in which the devices to be managed are taken from a small, fixed collection of emulated components. This would seem to preclude the use of virtual machines in testing system-level performance, such as in a comparative study of scheduling algorithms. The devices of interest are unlikely to be among those emulated, and, even if they are, the emulated devices may have arbitrary implementation, thus exhibiting performance characteristics that bear little or no resemblance to those of the real devices. For example, an entire virtual disk may be cached in the main memory of a large NAS device supporting the virtual machine, thus providing a disk with constant (O(1)) service times.

In [START_REF] Geist | Virtualizing high-performance graphics cards for driver design and development[END_REF], we introduced the iprobe, an extremely light-weight extension to the Linux kernel that can be dynamically loaded and yet allows the interception and replacement of arbitrary kernel functions. The iprobe was seen to allow a straightforward implementation of PCI device emulators that could be dynamically loaded and yet were suitable for full, system-level driver design, development, and testing.

The goal of this effort is to extend the use of the iprobe to allow performance prediction for real devices using only virtual machines. We suggested this possibility in [START_REF] Geist | Virtualization of an advanced course in operating systems[END_REF]. In particular, we will predict the performance of a collection of disk scheduling algorithms, one of which is new, for a targeted file-server workload on a specific SCSI device. We then compare the results with measurements of the same algorithms on the real hardware. We will see that the advantage of our approach is that reasonable performance predictions may be made without dedicated measurement platforms and with relatively limited knowledge of the performance characteristics of the targeted devices.

The remainder of the paper is organized as follows. In the next section we provide background on both kernel probes and disk scheduling. In section 3, we propose a new disk scheduling algorithm that will serve as the focus of our tests.

In section 4, we show how the kernel probes may be used to implement a virtual performance throttle, the mechanism that allows performance prediction from virtual platforms. In section 5 we describe our virtual machine, the targeted (real) hardware, and a test workload. Section 6 contains results of algorithm performance prediction from the virtual machine and measurements of the same algorithms on the real hardware. Conclusions follow in section 7.

Background

Kernel Probes

The Linux kernel probe or kprobe utility was designed to facilitate kernel debugging [START_REF] Coffman | On the expected performance of scanning disks[END_REF]. All kprobes have the same basic operation. A kprobe structure is initialized, usually by a kernel module, i.e., a dynamically loaded kernel extension, to identify a target instruction and specify both pre-handler and post-handler functions. When the kprobe is registered, it saves the targeted instruction and replaces it with a breakpoint. When the breakpoint is hit, the pre-handler is executed, then the saved instruction is executed in single step mode, then the post-handler is executed. A return resumes execution after the breakpoint.

A variation on the kprobe, also supplied with Linux, is the jprobe, or jump probe, which is intended for probing function calls, rather than arbitrary kernel instructions. It is a kprobe with a two-stage pre-handler and an empty posthandler. On registration, it copies the first instruction of the registered function and replaces that with the breakpoint. When this breakpoint is hit, the firststage pre-handler, which is a fixed code sequence, is invoked. It copies both registers and stack, in addition to loading the saved instruction pointer with the address of the supplied, second-stage pre-handler. The second-stage pre-handler then sees the same register values and stack as the original function.

In [START_REF] Geist | Virtualizing high-performance graphics cards for driver design and development[END_REF], we introduced the intercept probe or iprobe, which is a modified jprobe. Our iprobe second-stage pre-handler decides whether or not to replace the original function. If it decides to do so, it makes a backup copy of the saved (function entry) instruction and then overwrites the saved instruction with a no-op. As is standard with a jprobe, the second-stage pre-handler then executes a jprobe return, which traps again to restore the original register values and stack. The saved instruction (which now could be a no-op) is then executed in single step mode. Next the post-handler runs. On a conventional jprobe, this is empty, but on the iprobe, the post-handler checks to see if replacement was called for by the second-stage pre-handler. If this is the case, the single-stepped instruction was a no-op. The registers and stack necessarily match those of the original function call. We simply load the instruction pointer with the address of the replacement function, restore the saved instruction from the backup copy (overwrite the noop) and return. With this method, we can intercept and dynamically replace any kernel function of our choice. Note that it is possible to have two calls to the same probed function, one that is to be intercepted and one that is not. A discussion of the handling of potential attendant race conditions in SMP systems may be found in [START_REF] Jones | A Framework for Virtual Device Driver Development and Virtual Device-Based Performance Modeling[END_REF].

Disk Scheduling

Scheduling algorithms that re-order pending requests for disks have been studied for at least four decades. Such scheduling algorithms represent a particularly attractive area for investigation in that the algorithms are not constrained to be work-conserving. Further, it is easy to dismiss naive (but commonly held) beliefs about such scheduling, in particular, that a greedy or shortest-access-time-first algorithm will deliver performance that is optimal with respect to any common performance measure, such as mean service time or mean response time. Consider a hypothetical system in which requests are identified by their starting blocks and service time between blocks is equal to distance. Suppose the read/write head is on block 100 and requests in queue are for blocks 20, 82, 120, and 200. The greedy schedule and the differing, optimal schedule are shown in Table 1. Table 1. Sub-optimal performance of the greedy algorithm.

Disk scheduling algorithms are well-known to be analytically intractable with respect to estimating response time moments. Early successes in this area, due to Coffman and Hofri [START_REF] Coffman | On the expected performance of scanning disks[END_REF] and Coffman and Gilbert [START_REF] Jr | Polling and greedy servers on a line[END_REF] were restricted to highly idealized, polling servers, in which the read/write head sweeps back and forth across all the cylinders, without regard to the extent of the requests that are actually queued.

Almost all knowledge of the performance of real schedulers is derived from simulation and measurement studies. Geist and Daniel described UNIX system measurements of the performance of a collection of "mixture" algorithms that blended scanning and greedy behavior [START_REF] Geist | A continuum of disk scheduling algorithms[END_REF]. Worthington, Ganger, and Patt [START_REF] Worthington | Scheduling algorithms for modern disk drives[END_REF] showed, in simulation, that scheduling with full knowledge of disk subsystem timing delays, including rotational delays and cache operations, could offer major performance improvements. They also concluded that knowledge of the cache operation was far more important than an accurate mapping of logical block to physical sector, a point which we will address.

More recently, Pratt and Heger [START_REF] Pratt | Workload dependent performance evaluation of the Linux 2.6 i/o schedulers[END_REF] provided a comparison of the four schedulers distributed with Linux 2.6 kernels. Until 2.6, Linux used a uni-directional or circular scan (CSCAN), in which requests are served in ascending order of logical block number until none remains, whereupon the read/write head sweeps back down to the lowest-numbered pending request. With recognition that the best scheduler is likely workload-dependent, Linux authors changed the 2.6 kernel to allow single-file, modular, drop-in schedulers that could be dynamically switched. Four schedulers were provided. The default is the completely fair queueing (cfq) algorithm, which has origins in network scheduling. Each process has its own logical queue, and requests at the front of each queue are batched, sorted and served.

The deadline scheduler was designed to limit response time variance. Each request sits in two queues, one sorted by CSCAN order, one FIFO, and each has a deadline. The CSCAN order is used, unless a deadline would be violated, and then FIFO is used. As with many algorithms, reads are separated from and given priority over writes because the requesting read process has usually suspended to await I/O completion, and so there are actually four (logical) queues. The anticipatory scheduler is no longer supported, and the noop scheduler is essentially FIFO, which delivers poor performance on almost all workloads, and thus these two will not be discussed.

A fundamental departure from greedy algorithms, scanning algorithms, and O(N ) mixtures thereof was offered by Geist and Ross [START_REF] Geist | Disk scheduling revisited: Can o(n 2 ) algorithms compete?[END_REF]. They observed that, over the preceding decades, CPU speeds had increased by several orders of magnitude while disk speeds remained essentially unchanged. They suggested that O(N 2 ) algorithms might be competitive and offered a statically optimal solution that was based on Bellman's dynamic programming [START_REF] Bellman | Dynamic Programming[END_REF], in which a table of size O(N 2 ) containing optimal completion sequences was constructed. Although their algorithm was shown to deliver excellent performance in tests on a real system, there were two easily-identifiable problems. It ignored the dynamics of the arrival process, and it ignored the effects of any on-board disk cache.

More recently, Geist, Steele, and Westall [START_REF] Geist | Enhancing webserver performance through the use of a drop-in, statically optimal, disk scheduler[END_REF] partially addressed the issue of arrival dynamics. Although at first glance counter-intuitive, it is often beneficial for disk schedulers with a non-empty queue of pending requests to do nothing at all [START_REF] Iyer | Anticipatory scheduling: a disk scheduling framework to overcome deceptive idleness in synchronous i/o[END_REF]. The process that issued the most recently serviced request is often likely to issue another request for a nearby sector, and that request could be served with almost no additional effort. They added a so-called busdriver delay, to mimic the actions of a bus driver who would wait at a stop for additional riders, to the table-based, dynamic programming algorithm of Geist and Ross, and showed that it delivered excellent performance, superior to any of the four schedulers distributed with Linux 2.6, for a fairly generic, web file-server workload designed by Barford and Crovella [START_REF] Barford | Generating representative web workloads for network and server performance evaluation[END_REF].

A New Scheduler

We now propose an extension of Geist-Steele-Westall algorithm to capture the effects of the on-board cache. On-board caches are common, although their effect on the performance of standard workloads is often minimal. All UNIX-derivative operating systems allocate a significant portion of main memory to I/O caching. In Linux this is called the page buffer cache. The file systems also issue readahead requests when they detect sequential reads of a file's logical blocks. Since the page buffer cache is usually an order of magnitude larger than any on-board disk cache, most of the benefits of caching are captured there. Nevertheless, users can force individual processes to avoid the page buffer cache, and so there are cases where the on-board disk cache could have significant effect.

We model the effects of the on-board cache with just three parameters, the number of segments, the number of sectors per segment, and the pre-fetch size, also in sectors. We assume the cache is fully associative with FIFO replacement, and, should a request exceed the segment size, we assume wrap-around. Although manufacturers are often secretive about the operations of on-board disk caches, the minimal information we require can usually be obtained from the SCSI mode page commands. The sdparm utility [START_REF] Gilbert | [END_REF] is a convenient tool for accessing such.

With this information we maintain a shadow cache within the scheduler. Upon entry to the scheduler, even if a previously computed optimal sequence is still valid and no O(N 2 ) table-building would be required, we check the entire arrival queue of pending requests against the shadow cache for predicted cache hits. If any request is found to be a predicted hit, it is scheduled immediately.

The shadow cache comprises pairs of integers denoting the start and end sectors of the span assumed to be contained within each segment. When a request is dispatched, the pair (start sector, end sector + pre-fetech) is written into the shadow cache at the current segment index, and the current segment index is advanced. To check for a predicted read hit, the requested span is compared against the spans contained in each segment.

Finally, we use a form of soft deadline to control response time variance. Reads are given priority over writes, and only reads use the table-building algorithm. Writes are served in CSCAN order. If the oldest pending read request exceeds a maximum READDELAY parameter, we forgo table-building and serve reads in CSCAN order as well, until the oldest reader no longer exceeds this age. The deadline is soft because it only guarantees service within the next sweep.

A Virtual Performance Throttle

We now describe how the real performance of these algorithms can be predicted using a virtual machine. We start by specifying a service time model for the targeted physical drive. A highly accurate model would require a detailed mapping of logical to physical blocks, which, for modern disks, is often serpentine in nature [START_REF] Qian | A Linux implementation validation of trackaligned extents and track-aligned raids[END_REF]. Within a single track, consecutive logical sectors usually map directly to consecutive physical sectors, but the mapping of logical tracks to physical tracks is considerably more complex because cylinder seek time is now less than head switching time. As a result, logical tracks are laid out in bands that allow multiple cylinder seeks per head switch when the disk is read in a logically sequential way. Nevertheless, as we will illustrate in the next section, at the macroscopic level, where seek distances are measured in millions of sectors, a linear model will suffice. Thus we assume we have a service time model of the form X r = R r /2 + S r (d r /D r ), where R r is rotation, S r is maximum seek time, and D r is maximum seek distance.

The idea of the virtual performance throttle (VPT) is to insert an iprobe into the SCSI path of the virtual system to force virtual service times that are proportional to real ones, with an identifiable constant of proportionality. If the virtual seek distance is d v with maximum D v , then we want the virtual system to deliver a service time of kX r , where d r /D r = d v /D v , for some system-wide constant, k. Instead, it will deliver a service time of X v . Clearly, we can achieve our goal by delaying the completed virtual request by kX r -X v , except that k is unknown to us. Further, due to factors outside our control, such as server load and network congestion on the hardware supporting the virtual machine, an appropriate k may change during the course of our tests. Thus we need a self-scaling system.

We specify an initial value of k and allow the VPT to constrain the flow of disk requests being served based on this value and the linear service time model. The VPT uses two kernel probes in the generic SCSI driver: a jprobe in the down path to record when requests leave for the virtual disk and an iprobe in the up path to intercept and delay completed requests upon return from the virtual disk. The jprobe calculates the target completion time of the request and passes it to the iprobe, which then determines how long the request should be delayed after its completion but before returning it to the requesting process. The iprobe's queue of delayed completions is checked periodically with a timer. Once the target completion time has passed, the request is injected back into the SCSI generic path.

As noted, the service time on the virtual disk is subject to change, and so it is possible that a request will arrive to the iprobe after its target completion time. This means that the scale factor, k, is too small. The iprobe increases it and passes the new value to the jprobe. Similarly, if the VPT develops a large queue of delayed requests, we know that the target performance is an underestimate, and the iprobe can decrease k to shorten the overall length of the simulation. The iprobe reports the current value of k on each change. In practice, we use preliminary tests with a dynamic k to find stable values and then fix k at a stable value for each real test.

Platform and Workload

The real test platform used in our study was a Linux (2.6.30) system with two, Intel Xeon 2.80GHz processors, 1 GB main memory, a Western Digital IDE system drive and two external Seagate Cheetah 15K.4 SCSI drives, each with its own Adaptec 39320A Ultra320 SCSI controller. Tests were restricted a single Cheetah drive. The disk is a model ST373454 with 4 recording surfaces and a formatted capacity of 73.4 GBytes. It rotates at 15,000 rpm yielding a rotation time of 4 ms. The disk has 50,864 tracks per recording surface and was formatted at 512 bytes per sector.

To approximate the linear service time at the macroscopic level, we disabled the on-board cache, opened /dev/sda using the O DIRECT mode, which forced a bypass of the page buffer cache, and read 100,000 randomly selected pages. The time required to read each page and the distance in sectors from the previously read page were captured. We sorted this data in order of increasing distance and plotted distance versus time. The result was a reasonably linear band of noise approximately 4 ms in width. The data was then smoothed using a filter that replaced each point with the average of the (up to) 1001 points centered at the point in question. The filtered data and the least squares approximation to it are plotted in Figure 1. The linear model is X r = 4.25 + 5.25(d r /D r ). The Cheetah manual [START_REF]Product Manual Cheetah 15K.4 SCSI. Pub[END_REF] indicates that 7,077KB is available for caching, which yields 221 512-byte sectors per segment. We assumed that a single span of requested sectors, plus the pre-fetch, would be cached in each segment. Single request spans larger than 157 sectors (221-64) were rare in our workloads, but for those cases we assumed a wrap-around with over-write within the segment.

The KVM-based virtual machine was hosted on an IBM 8853AC1 dual-Xeon blade. It was configured with a 73GB virtual SCSI disk, for which the available emulator was an LSI Logic / Symbios Logic 53c895. The virtual disk image was stored on a NetApp FAS960c and accessed via NFS.

We compared the performance of our cache-aware, table-scheduling (CATS) algorithm with cfq and deadline on both real and virtual platforms. The CATS busdriver delay was set to 7ms and the maximum READDELAY to 100ms.

We used two workloads in this study. Both were similar, at the process level, to that used by Geist, Steele, and Westall [START_REF] Geist | Enhancing webserver performance through the use of a drop-in, statically optimal, disk scheduler[END_REF], which was based on the approach used by Barford and Crovella [START_REF] Barford | Generating representative web workloads for network and server performance evaluation[END_REF] in building their Scalable URL Reference Generator (SURGE) tool. Each of 50 processes executed an ON/OFF infinite request loop, shown in pseudo-code in Figure 2.

The Pareto(α,k) distribution is a heavy-tailed distribution,

F X (x) = 1 -(k/x) α x ≥ k 0 elsewhere (1)
The discrete Zipf distribution is given by where k is a normalizing factor, specifically, the reciprocal of the N +1 st harmonic number. A continuous approximation,

p(i) = k/(i + 1), i = 0, 1, ..., N (2)
F X (x) = log(x + 1) log(N + 2) 0 ≤ x ≤ N + 1 (3)
suffices for our study. File count parameters were taken directly from the Barford and Crovella study, (α 1 , k 1 ) = (2.43,1.00). The shape parameter of the sleep interval, α 2 = 1.50, was also taken from this study, but we used a different scale parameter, k 2 = 2.0, because our sleep interval was milliseconds per block rather than seconds per file.

Both the virtual SCSI drive and the Cheetah drive were loaded with 1 million files in a two-level directory hierarchy where file sizes were randomly selected from a mixture distribution also suggested by Barford and Crovella. This mixture distribution is lognormal(9.357,1.318) below 133 KB and Pareto(1.1,133K) above, where the lognormal(µ,σ) distribution function is given by:

F Y (y) = y 0 e -(loget-µ) 2 2σ 2 /(tσ √ 2π)dt y > 0 (4)
To induce reasonable fragmentation, we erased a half million files, selected at random, and then added back a half million files with different, randomly selected sizes. Due to the relatively small capacity of these drives, we chose to truncate at 100 MB any files that would have exceeded that size.

For each algorithm, for each test run, we captured the arrival times, service initiation times, and service completion times of 50,000 requests. We captured these time stamps by directly instrumenting the kernel outside of the schedulers, and, during each test run, we stored the time stamps to a static kernel array. The time stamp data was extracted from the kernel array after the test run by using a custom system call.

The two workloads were identical at the process level but decidedly different at the drive level. For the first, each file was opened with mode flag O DIRECT, which, as noted earlier, forced the associated I/O to by-pass the main memory page buffer cache. The second workload differed from the first only in that the O DIRECT flag was not used. The page buffer cache, a standard feature of UNIX-derivative operating systems, is dynamic and can grow to become quite large. For tests described here, 65MB was often observed.

Finally, although the Cheetah drive supports tagged command queueing (TCQ), we disabled it for all tests. We found that, for all schedulers, allowing re-scheduling by the drive hardware decreased performance. We would have thought this to be an anomaly, but we have observed the same result for other SCSI drives on other Linux systems.

Results

The results for the first workload, using O DIRECT, are shown in Table 2 see that the virtual system uniformly predicted higher mean service, higher mean response, and lower throughput than was found from measurements of the real system. Nevertheless, on all three measures, the predicted performance rank of the three algorithms was correct: CATS performs better than deadline, which performs better than the Linux default, cfq. Thus algorithm selection could be made solely on the basis of the virtual system predictions. We also gauged the effectiveness of the shadow cache in predicting real cache hits by placing record markers within the captured time stamp trace of the CATS scheduler on all records that were predicted by the shadow cache to be hits. We examined the service time distribution of a large collection of singlesector requests, independent of any trace, and found a prominent initial spike at 250 microseconds. We then processed the CATS trace and marked any record with service time below 250 microseconds as an actual hit. We found that the shadow cache correctly predicted 97% of the 31,949 actual hits observed.

When I/O is staged through the main memory page buffer cache (the second workload), the results are decidedly different, as shown in Table 3. Again the virtual system uniformly overestimated mean service time and mean response time and underestimated throughput, but again it correctly predicted the performance rank of all three algorithms on all three measures.

Conclusions

We have suggested a method for predicting the performance of disk scheduling algorithms on real machines using only their measured performance on virtual machines. The method uses a dynamically loaded kernel intercept probe (iprobe) to adjust low-level virtual device timing to match that of a simple model derived from the real disk device. We used this method to predict the performance of three disk scheduling algorithms, one of which is new. Although the virtual system was seen to uniformly underestimate performance, it correctly predicted the relative performance of the three algorithms as measured on a real system under two workloads. It it fair to charge that we are simply using a virtual operating system as an elaborate simulator. Nevertheless, this simulator provides almost all the subtle nuances of a real operating system and yet requires almost no programming effort on our part. A low-level model of device service time performance and the drop-in iprobe are all that is required.

We are currently working on methods to achieve more accurate absolute predictions. As yet we have not accounted for measurement overhead inherent in our method. We believe that this is one of the causes of the uniformly underestimated performance. Another potential source of error is the model service time ascribed to a cache hit, which we fix at 250 microseconds, even though this only an observed maximum. Unfortunately, reducing this value, X r , will require an increase in scale factor, k, so that kX r -X V remains non-negative, and increasing k increases run-time. We may be faced with a trade-off between accuracy and (simulation) run-time. Nevertheless, we believe that accounting for these factors will result in more accurate predictions.

Fig. 1 .

 1 Fig. 1. Stochastic sector to sector costs

Table 2 .

 2 . We Performance on O DIRECT workload.

			real		virtual (k=8)
	algorithm	cats deadline	cfq	cats deadline	cfq
	mean service (ms)	1.96	2.71	1.39	2.58	3.24	2.36
	variance service	8.51	9.76	5.85	9.03	8.23	7.78
	mean response (ms)	37.35	59.87	124.70	53.79	78.27 117.13
	variance response	6961.50 561.15 839270.49 16641.07 633.28 28651.71
	throughput (sectors/ms)	8.19	6.08	2.19	6.15	5.06	3.38

Table 3 .

 3 Performance on non-O DIRECT workload.

			real		virtual (k=8)
	algorithm	cats deadline	cfq	cats deadline	cfq
	mean service (ms)	6.53	7.41	7.80	7.15	7.60	8.57
	variance service	11.13	8.80	17.62	6.22	6.05	10.30
	mean response (ms)	114.91 121.87 179.17 189.45 198.33 258.75
	variance response	8080.16 3296.87 35349.39 19292.52 6839.66 65796.33
	throughput (sectors/ms) 12.00	12.04	9.08	11.44	11.68	8.82


Vicente Casares-Giner
email: vcasares@upvnet.upv.es

On Lookahead Strategy for Movement-based Location Update. A General Formulation

Keywords:

Location management deals with the procedure to update the current location of a mobile terminal (MT) and with the procedure to deliver incoming calls to that called MT. Basically, the performance evaluation of location management procedures are dependent on the MT's mobility behavior, on the MT's cell residence time and on the call arrival process to the MT. In the open literature, the typical analysis for location management has been addressed under the assumptions of a random walk mobility model, the exponential cell sojourn time of the MT -some times relaxed to a general probability distribution functionand the exponential inter-arrival time distribution for incoming calls. However, the random walk model seems not to be valid as many mobile users follow some daily trajectories, such as from home to the working place, from the working place to the shopping center, etc. To reflect a more realistic movement pattern, we propose a directional oriented mobility model. And as a consequence of that, we also propose a lookahead procedure combined with the movement-based location update scheme, with the main idea of saving signaling traffic through the air interface. In the lookahead strategy we analytically derive closed form expressions for the mean number of location update (LU) messages triggered by the MT between two consecutive call arrivals and the paging probabilities to evaluate the paging cost under some selective paging strategies. The analysis has been carried out assuming a general cell residence time and a renewal point process for call arrivals to the MT.

Introduction

Mobility management is a key functionality that aims in providing ubiquitous telecommunication services in wireless mobile access networks. Mobility management procedures can be classified as handoff management and location management. The first procedure aims to guarantee the continuity of an on going call when the mobile terminal (MT) visits a new cell. Location management is defined as a set of procedures that allow a MT being located at anytime, anywhere, so that incoming calls may be delivered to that MT. This set consists of location update (LU) and call delivery (CD). The LU process consists of maintaining the MT location information updated in the system database. The database entry of an MT is updated whenever the MT triggers an LU message or an incoming call is delivered. The CD procedure is decomposed into two steps: interrogation and terminal paging (PG). Firstly, in the interrogation step, the system database is queried to obtain the Registration Area (RA) where the MT is supposed to be in. Afterwards, the PG procedure follows, i.e., the MT is searched by polling the set of cells of the RA.

There is a trade off between LU and PG procedures. As the number of LU messages increases, the LU cost becomes higher but the PG cost decreases as the MT position is known more accurately. On the other hand, the lower is the number of LUs, the lower is the LU cost; but the uncertainty of the MT position is higher and the PG cost increases.

The performance of a location management procedure is dependent on the mobility patterns of the MT. One of the most extended mobility models used in the evaluation of location tracking is the random walk mobility model, [START_REF] Camp | A survey of mobility models for ad hoc network research[END_REF] [START_REF] Bettstetter | Stochastic properties of the random waypoint mobility model[END_REF]. However, the random walk model is rather unrealistic since many MTs follow daily trajectories such as, as we have said before, from home to the working place, from home to the recreation center, etc. Therefore, in order to evaluate in a proper way the signaling load associated to mobility tracking, more realistic mobility models are required [START_REF] Martinez-Arrue | Location Management Based on the Mobility Patterns of Mobile Users[END_REF].

In this work we consider a mobility model with a certain directionality, and as a consequence of that information, an enhanced movement-based mobility tracking procedure is proposed. The proposal is based on the lookahead tracking principle and it is illustrated and evaluated in a 2-D scenario.

Movement-based Location Update

Scenario and mobility model

Although hexagonal and square cell layouts in a 2-D scenario can be considered, here we analyze our proposal in the hexagonal cell layout [START_REF] Escalle-García | Reducing location update and paging costs in a Personal Communications Services Network[END_REF]. Cells are located in a plane and they are enumerated according to Fig. 1. The mobility model is as follows. After a certain sojourn or cell dwell time in cell (i, j), characterized by a probability density function (pdf) f m;i,j (t), our MT will travel in the direction of angles {iπ/3} with probabilities {p i }, counterclockwise, for i = 0, 1, 2, 3, 4, 5, respectively. We will assume that all cells have the same size, and that f m;i,j (t) = f m (t). The following probabilities {p i } are considered for our single parameter mobility model

p 0 = β 3 D(β) ; p 1 = p 5 = β 2 D(β) ; p 2 = p 4 = β D(β) ; p 3 = 1 D(β) with D(β) = (β 3 + 2β 2 + 2β + 1) = (β + 1)(β 2 + β + 1);
where β is positive number, β ∈ (0, 1) or β ∈ (1, ∞), and D(β) is a normalization factor. Clearly, the parameter β provides a wide range of the directionality of the MT's movement. For instance, if β = 1, p i = 1/6, ∀i, we have the standard random walk model. For β > 1 we have p o > p 1 = p 5 > p 2 = p 4 > p 3 . Clearly when β → ∞ p 0 → 1 the MT will always travel in the same direction, i.e. according to Fig. 1, towards the right hand side. The higher is β the higher the directionality of our MT will be.

Lookahead strategy

Let us consider that at time t = 0 + our MT is located at cell (i, j) = (0, 0), see Fig. 1. We choose a set of target cells, the ∆ set, that our MT will likely find along its trajectory. For example, if β > 1 cells (1, 1), (2, 2), (3, 3), ... can be reached by our MT with higher probability than cells (-1, -1), (-2, -2), (-3, -3), ... . Then, each time the MT visits a cell that does not belong to the ∆ set, its movement counter is increased in one unit. When the movement counter reaches a certain threshold, say M , it will send a LU message to the network. But, if the MT visits any cell of the ∆ set before the threshold M is reached, it will reset its movement counter. Fig. 1 depicts the case of M = 4 and ∆ = {(3, 3), [START_REF] Martinez-Arrue | Location Management Based on the Mobility Patterns of Mobile Users[END_REF][START_REF] Escalle-García | Reducing location update and paging costs in a Personal Communications Services Network[END_REF], [START_REF] Escalle-García | Reducing location update and paging costs in a Personal Communications Services Network[END_REF][START_REF] Martinez-Arrue | Location Management Based on the Mobility Patterns of Mobile Users[END_REF], (4, 4)}. Obviously the distance d ∆ between the starting point, the cell (0, 0) and the nearest target cell, measured in terms of cells, must be such that 0 ≤ d ∆ ≤ M ; d ∆ = 3 in Fig. 1. The case d ∆ > M is equivalent to having no target cells since, in this case, our MT will trigger LU messages according to the standard movement-based scheme, i.e. no absorption into the ∆ set will happen. Notice that, each time the MT triggers a LU message the procedure starts anew at cell (0, 0), i.e. the scenario depicted in Fig. 1 is repeated.

Modeling and analysis

In this paper we generalize the work [START_REF] Casares-Giner | On Movement-Based Location Update. A Lookahead Strategy[END_REF] in two aspects; first by allowing more than one cell in the ∆ set, secondly by considering a general renewal point process for modeling call arrivals to the MT. Then, for a predetermined number of movements m we are interested in the number of LU messages triggered by our MT between two consecutive call arrivals. We say that our MT is in state S i if it is roaming in cell i. Although a natural cell identification in a 2-D space is the use of two indexes, for the sake of simplicity in the mathematic notation we will use a single index. Since the RA where the MT is roaming is configured by a finite number of cells, a correspondence between the identification with two indexes and the identification with a single index can be established. Then, starting from state S 0 , cell (0, 0) in Fig. 1, if an absorption into the ∆ set happens at movement M or before, the MT will reset its movement counter. If the MT visits a cell that does not belong to the ∆ set it will increase in one unit its movement counter. When the movement counter reaches the threshold M without being absorbed into the ∆ set, the MT will trigger a LU message and the process starts anew at state S 0 .

The analysis can be formulated by building a suitable Markov chain and using some useful taboo probabilities [START_REF] Karlin | A first course in stochastic processes[END_REF]. We call H f (n) i,j the conditional probability that the Markov chain enters state S j for the first time at the nth step, having initially started from state S i and avoiding the set of states H.

With no loss of generality, in the sequel we assume that the ∆ set of absorbing states is configured by three cells, j, k and l, that is, ∆ = {j, k, l}. Also, we denote by ∆ s the complementary set of cells in ∆ with respect to cell s. In our case

∆ j = {k, l}, ∆ k = {j, l} and ∆ l = {j, k}.
Then, if initially the MT is in state S 0 the probability that in M movements no absorption into the ∆ set occurs, P nab(0,∆) (M ), is given by

P nab(0,∆) (M ) = 1 - M n=1 f (n) 0,∆ = 1 - M n=1 ∆j f (n) 0,j + ∆ k f (n) 0,k + ∆ l f (n) 0,l (1)
Also, if initially the MT is in state S s , for (s = j, k, l) the probability that in M movements no absorption into the ∆ set occurs, P nab(s,∆) (M ), is given by

P nab(s,∆) = 1 - M n=1 f (n) s,∆ = 1 - M n=1 ∆j f (n) s,j + ∆ k f (n) s,k + ∆ l f (n) s,l ; for s = j, k, l (2)
NOTE: Details on how the taboo probabilities

{ H f (n)
i,j } are evaluated are given in the internal report [START_REF] Casares-Giner | On lookahead strategy for movement-based location update. A general formulation[END_REF].

3 The number of location update messages , for s = j, k, l) denote the mean number of LU messages. Then, we have

M (m) 0 =    0; for m < M M n=1 f (n) 0,∆ M (m-n) ∆ + P nab(0,∆) (M )[1 + M (m-M ) 0 ]; for m ≥ M (3)
and

M (m) ∆ =    0; (a column vector of zeros); for m < M M n=1 f (n) ∆,∆ M (m-n) ∆ + P nab(∆,∆) (M )[1 + M (m-M ) 0 ]; for m ≥ M (4)
where (more details are given in [START_REF] Casares-Giner | On lookahead strategy for movement-based location update. A general formulation[END_REF])

M (m) ∆ = M (m) j , M (m) k , M (m) l T = ∞ i=0 iP (m) ∆ (i); with P (m) ∆ (i) = P (m) j (i), P (m) k (i), P (m) l (i)
T ; P nab(∆,∆) (M ) = P nab(j,∆) (M ), P nab(k,∆) (M ), P nab(l,∆) (M ) T ; f (n) 0,∆ = ∆j f (n) 0,j , ∆ k f (n) 0,k , ∆ l f (n) 0,l ; f (n) ∆,∆ =     ∆j f (n) j,j ∆ k f (n) j,k ∆ l f (n) j,l ∆j f (n) k,j ∆ k f (n) k,k ∆ l f (n) k,l ∆j f (n) l,j ∆ k f (n) l,k ∆ l f (n) l,l    
Therefore, the mean number of LU messages triggered by the MT between two consecutive call arrivals, # LU , can be expressed as

# LU = ∞ m=0 M (m) 0 α(m) = ∞ m=M M (m) 0 α(m) (5)
where α(m) is the probability of m cell boundary crossings between two consecutive call arrivals. α(m) is derived in the next section.

On the number of cell boundary crossing between two call arrivals

We assume that calls arrive to the MT following a renewal point process with inter-arrival time pdf given by f c (t). We denote f * c (s) the Laplace transform (LT) of f c (t). Then, α(k) is given by (see [START_REF] Casares-Giner | On lookahead strategy for movement-based location update. A general formulation[END_REF]),

α(k) =          1 2πi c+i∞ c-i∞ [1 -f * mr (s)] s f * c (-s)ds; for k = 0 1 2πi c+i∞ c-i∞ f * mr (s)[1 -f * m (s)][f * m (s)] k-1 s f * c (-s)ds; for k > 0 (6)
In ( 6), f * m (s) is the LT of f m (t); the pdf of the MT dwell time or residence time in a cell, with mean value 1/λ m = -f * ′ m (0). f * mr (s) is the LT of the residual cell residence time, f * mr (t). Following the residual life theorem [START_REF] Kleinrock | Queuing Theory[END_REF] we have

f mr (t) = λ m ∞ t f m (t)dt = λ m [1 -F m (t)]; f * mr (s) = λ m s [1 -f * m (s)] (7)

Mean number of LU between two consecutive call arrivals

In order to obtain a closed form expression for (5) first we define the following generating functions

M * 0 (z) = ∞ m=0 M (m) 0 z m ; M * ∆ (z) = ∞ k=0 M (k) ∆ z . (8)
F * 0,∆ (z) = M n=1 f (n) 0,∆ z n = F * 0,j (z), F * 0,k (z), F * 0,l (z) ; with F * 0,s (z) = M n=1 ∆s f (n) 0,s z n , for s ∈ ∆ ≡ {j, k, l} F * ∆,∆ (z) = M n=1 f (n) ∆,∆ z n =   F * j,j (z) F * j,k (z) F * j,l (z) F * k,j (z) F * k,k (z) F * k,l (z) F * l,j (z) F * l,k (z) F * l,l (z)   ; with F * r,s (z) = M n=1 ∆s f (n) r,s z n , for r, s ∈ ∆ ≡ {j, k, l}
Then, from (3) and (4) a linear relationship between M * 0 (z) and

M * ∆ (z) is obtained. That is M * 0 (z) = F * 0,∆ (z)M * ∆ (z) + P nab(0,∆) (M )z M 1 1 -z + M * 0 (z) M * ∆ (z) = F * ∆,∆ (z)M * ∆ (z) + P nab(∆,∆) (M )z M 1 1 -z + M * 0 (z) (9)
Therefore, solving (9) we get the expression for M * 0 (z)

M * 0 (z) = = 1 1 -z F * 0,∆ (z)[I -F * ∆,∆ (z)] -1 P nab(∆,∆) (M ) + P nab(0,∆) (M ) z M 1 -F * 0,∆ (z)[I -F * ∆,∆ (z)] -1 P nab(∆,∆) (M ) + P nab(0,∆) (M ) z M (10)
Secondly, from ( 3)-( 6) we can write, after some algebra, [START_REF] Casares-Giner | On lookahead strategy for movement-based location update. A general formulation[END_REF],

# LU = ∞ m=0 M (m) 0 α(m) = ∞ k=M M (k) 0 α(k) = = 1 2πi c+i∞ c-i∞ f * mr (s)[1 -f * m (s)] sf * m (s) M * 0 (f * m (s))f * c (-s)ds = - p∈σc Res s=p f * mr (s)[1 -f * m (s)] sf * m (s)] M * 0 (f * m (s))f * c (-s) (11)
The last equality in (11) is derived applying the Cauchy's Residue Theorem, i.e., σ c denotes the set of poles of f * c (-s) in the right half complex plane, and Res s=p denotes the residue at poles

s = p ∈ σ c . M * 0 (f * m (s)) is the expression (10) evaluated at z = f * m (s).

Terminal paging procedure

Blanket or non selective and selective or multi-step PG procedures are considered. For PG evaluation, in the first case we only need the number of cells that conforms the RA. Once the parameter M and the ∆ set are fixed, the total number of cells to poll can easily be obtained. For instance, when M = 4 and ∆ = {(3, 3), [START_REF] Martinez-Arrue | Location Management Based on the Mobility Patterns of Mobile Users[END_REF][START_REF] Escalle-García | Reducing location update and paging costs in a Personal Communications Services Network[END_REF], (4, 3), (4, 4)} Fig. 1 accounts for the set of cells, the Ω set, in which the MT can be in. The cardinality of Ω is C(Ω) = 73. Obviously, for selective PG we need the PG probabilities. They are derived in the next lines.

Paging probabilities

Let π

(m) i,t denotes the probability that starting at cell i the MT is located at cell t after m movements. Then, following the arguments used in section 3.1 we can write, for i = 0

π (m) 0,t =      p (m) 0,t ; for m < M M n=1 f (n) 0,∆ π (m-n) ∆,t + P nab(0,∆) (M )π (m-M ) 0,t ; for m ≥ M (12)
and for the ∆ set, cells i = s ∈ ∆ ≡ {j, k, l}

π (m) ∆,t =      p (m) ∆,t ; for m < M M n=1 f (n) ∆,∆ π (m-n) ∆,t + P nab(∆,∆) (M )π (m-M ) 0,t ; for m ≥ M (13) with π (m) ∆,t = π (m) j,t , π (m) k,t , π (m) l,t T ; p (m) ∆,t = p (m) j,t , p (m) k,t , p (m) l,t T ;
Therefore the set of probabilities {π i,t } can be expressed as

π i,t = ∞ m=0 π (m) i,t α(m) (14)
For i = 0, (14) gives the probability that the MT is located at cell j when a call arrival occurs. To evaluate the set {π 0,t } we follow a similar procedure to that in section 3. To that end, first, let us define the following generating functions, π * 0,t (z) and π * ∆,t (z)

π * 0,t (z) = ∞ m=0 π (m) 0,t z m ; π * ∆,t (z) = ∞ m=0 π (m) ∆,t z m ; (15)
Then, from equations ( 12) and (13) a linear relationships between π * 0,j (z) and π * ∆,t (z) is obtained. That is

[1 -z M P nab(0,∆) (M )]π * 0,t (z) -F * 0,∆ (z)π * ∆,t (z) = = M -1 n=0 p (n) 0,t z n - M -1 n=1 f (n) 0,∆ M -n-1 k=0 p (k) ∆,t z k z n = Υ (0, ∆, z, t); -z M P nab(∆,∆) (M )π * 0,t (z) + [I -F * ∆,∆ (z)]π * ∆,t (z) = = M -1 n=0 p (n) ∆,t z n - M -1 n=1 f (n) ∆,∆ M -n-1 k=0 p (k) ∆,t z k z n = Γ(∆, ∆, z, t); (16)
Therefore, solving (16) for π * 0,t (z) we have

π * 0,t (z) = = F * 0,∆ (z)[I -F * ∆,∆ (z)] -1 Γ(∆, ∆, z, t) + Υ (0, ∆, z, t) 1 -F * 0,∆ (z)[I -F * ∆,∆ (z)] -1 P nab(∆,∆) (M ) + P nab(0,∆) (M ) z M (17)
Secondly, from (12)-( 13) and ( 6) the paging probabilities, (14) for i = 0, can be written as

π 0,t = ∞ k=0 π (k) 0,t α(k) = π (0) 0,t α(0) + ∞ k=1 π (k) 0,t α(k) = = 1 2πi c+i∞ c-i∞ [f * m (s) -f * mr (s)] sf * m (s) p (0) 0,t + f * mr (s)[1 -f * m (s)] sf * m (s) π * 0,t (f * m (s)) f * c (-s)ds = - p∈σc Res s=p [f * m (s) -f * mr (s)] sf * m (s) p (0) 0,t + f * mr (s)[1 -f * m (s)] sf * m (s) π * 0,t (f * m (s)) f * c (-s)
(18) As in (11), Cauchy's Residue Theorem is used to reach the last equality in (18).

Paging procedure

We have evaluated the PG cost for the conventional blanket or single step paging, η = 1 and for two selective paging schemes with delays η = 2 and η = 3. In the illustrative example we have used a single absorbing state for the ∆ set. Also, we have implemented the line-paging philosophy suggested in [START_REF] Hwang | A direction-based location update scheme with a line-paging strategy for PCS networks[END_REF]. To that end we have defined three sets of disjoint cells, A, B and C, that is,

A ∩ B = ⊘, A ∩ C = ⊘, B ∩ C = ⊘,
C{B} = 2C{A} + 4; if ∆ < M -1 2C{A} + 1; if M -1 ≤ ∆ ≤ M , C = {The remaining set of cells: Ω -A ∪ B}, C{C} = C{Ω} -C{A ∪ B} = C{Ω} -C{A} -C{B} (19)
and C{Ω} = 3(M -1) 2 + 3(M -1) + 1 + ∆(2M -1).

Numerical Analysis

The following assumptions have been established. A Gamma distributed for the cell residence time, with pdf and LT given by with mean value 1/λ m and with the squared coefficient of variation C 2 m = 1/γ, [START_REF] Kleinrock | Queuing Theory[END_REF]. For the call inter-arrival time, on one hand, the Erlangian-Geo-type f c,e (t), with LT

f m (t) = (γλ m ) γ t γ-1 Γ (γ) e -γλmt , f * m (s) = γλ m s + γλ m γ (20)
f * c,e (s) = n-1 i=0 λ i s + λ i , λ i = λ c a i 1 -a n 1 -a (21)
with mean value 1/λ c and

C 2 c,e = (1 + a n )(1 -a) (1 -a n )(1 + a)
≤ 1 for 0 < a < 1 and n any natural number. On the other hand, the hyper-exponential distribution H 2 , f c,h (t), with LT

f * c,h (s) = a λ 1 s + λ 1 + (1 -a) λ 2 s + λ 2 , λ 2 = nλ 1 , 1 λ c = a λ 1 + (1 -a) nλ 1 (22)
with mean value equal to 1/λ c and

C 2 c,h = 1 + 2an(n -1) -a 2 (n -1) 2 1 + 2a(n -1) + a 2 (n -1) 2 ≥ 1 for 0 < a < 1, equivalently 1 n ≤ λ 1 λ c ≤ 1,
and n any natural number. It is easy to verify that the maximum C 2 h is achieved for a = 1/(n + 1). All plots have been obtained for a Call-to-Mobility ratio (CMR) equals to λ c /λ m = 0.10. From Fig. 3 we observe that for C 2 m ≤ 1 the mean number of LU messages, # LU in [START_REF] Casares-Giner | On Movement-Based Location Update. A Lookahead Strategy[END_REF], reduces when C 2 c increases. However, when C 2 m > 1 the mean number of LU messages, # LU is rather insensitive to C 2 c . A significant saving is achieved when selective paging is implemented as it can be seen in Fig. 4 and Fig. 5. The saving is quite noticeable from blanket polling (not shown in the paper, since it is quite trivial) to the two step PG or polling scheme. Also a noticeable saving is obtained from two step polling to three step polling. As additional observation we realize that, for a fixed number of PG steps, and for a fixed C 2 m , the paging cost is no sensitive to the interarrival process. Also, for a fixed number of paging steps, for a fixed value of C 2 c and for C 2 m ≤ 1 (C 2 m > 1) the PG cost is rather insensitive (quite sensitive) to the inter-arrival time distribution.

Conclusions

In this paper we have proposed and studied within a general framework, a lookahead location update procedure with selective paging. We have seen that the information about the directionality of the mobile terminal's movement can be used to save signaling cost associated to movement-based location management procedures. With a lookahead scheme the size of the registration area increases and a selective paging scheme is recommended to avoid a noticeable increment in the paging cost. The net effect is a significant saving in the total signaling traffic on the common air interface. We believe that this scheme could be used by mobile operators where mobile terminals follow certain daily or customary trajectory.

Fig. 1 .

 1 Fig. 1. Illustration of the lookahead strategy for movement based location update, with movement-counter threshold M = 4 and ∆ = {(3, 3), (3, 4), (4, 3), (4, 4)}.

3. 1

 1 Mean number of LU messages for a fixed number of movementsLet P , for s = j, k, l) denote the probability of having i LU messages triggered in m movements, under the assumption that the MT starts its trajectory in cell 0 (s = j, k, l). Let M

  and A ∪ B ∪ C = Ω, being ⊘ the empty set and Ω the full set of cells where the MT can be in. Denoting by C{X} the cardinality of X, we have chosen A = {(0, 0), (1, 1), ... , (min{∆ + M -1, 2∆}, min{∆ + M -1, 2∆})}, C{A} = min(∆ + M -1, 2∆) + 1; B = {The first ring of cells around the set A},

Fig. 2 .

 2 Fig. 2. Paging schemes for M = 5 and ∆ = {(3, 3)}. Cells are grouped into three sets and they are sorted for paging procedure. Blanket polling: set Ω = A ∪ B ∪ C. Two step paging: first A, secondly B ∪ C. Three step paging: first A, secondly B, thirdly C.

βFig. 3 .

 3 Fig. 3. Average number of LU messages between two consecutive call arrivals in terms of the mobility model parameter β with the target ∆ cell as varying parameter (∆ = i means the target cell (i, i) in Fig. 2). Inter-arrival call process: Erlangian-Geo with n= 10 stages and C 2 c,e = 0.109 (low-variance) on the left and hyper-exponential with n=20 and C 2 c,h = 10.025 (high-variance) on the right. Cell residence time, from top to bottom: Gamma distribution with C 2 m = 0.01 (low-variance) , C 2 m = 1 (exponential) and C 2 m = 100 (hig-variance), respectively.

ββββββFig. 4 .

 4 Fig. 4. Two steps paging. Average number of cells polled per call arrival in terms of the mobility model parameter β with the target ∆ cell as varying parameter (∆ = i means the target cell (i, i) in Fig. 2). Inter-arrival call process: Erlangian-Geo with n= 10 stages and C 2 c,e = 0.109 (low-variance) on the left and hyper-exponential with n=20 and C 2 c,h = 10.025 (high-variance) on the right. Cell residence time, from top to bottom: Gamma distribution with C 2 m = 0.01 (low-variance) , C 2 m = 1 (exponential) and C 2 m = 100 (hig-variance), respectively.

Acknowledgments

This research has been financed by the Spanish Ministerio de Ciencia e Innovación project TSI2007-66869-C02-02 and by the Universidad Politécnica de Valencia under the grant PAID-05-08. Also, thanks for the support received from the European Network of Excellence Euro NF and to its research community for many scientific and fruitful interactions with its partners.


Axelle Amon
email: axelle.amon@univ-rennes1.fr

Baptiste Blanc

Jean-Christophe Géminard

Avalanche precursors in a frictional model

Keywords: numbers: 45.70.Ht, 45.70.-n, 46.55.+d, 45.05.+x

Avalanche precursors in a frictional model

I. INTRODUCTION

Identifying precursors to avalanches has inspired numerous works with the hope of being able to detect catastrophic events before they occur. A physicist approach consists in simplifying the system in order to single out the fundamental mechanisms underlying the phenomenon. For example, models based on cellular automaton have been proposed to study avalanches [1]. From an experimental point of view, a model set-up to study the behavior of a granular material before an avalanche takes place consists in the progressive inclination of a box filled with a granular material [2][3][4][START_REF] Duranteau | AIP Conference Proceedings[END_REF][START_REF] Gravish | [END_REF]. In such experiments, it has been evidenced that the response of the system consists in the superposition of two different behaviors: on the one hand, small rearrangements implying only a few number of grains occuring without any obvious synchronisation across the system and, on the other hand, large correlated events implying a large fraction of the system. Those last events have been shown to emerge from an angle at about half the avalanche angle [4] and then to occur at regular angle increments as the inclination increases until the destabilization of the pile [2][3][4]. To our knowledge, no explanation of the regularity of those pre-avalanche events exists in the litterature, nor any numerical observation of the phenomenon [7][8][9]. Yet, such micro-envents have been observed in different loading configurations [2,10,11] and are reminiscent of precursors observed in studies of the onset of frictional sliding [12].

Here, we present a one-dimensional frictional model that reproduces most of the features reported in the experiments. In particular, the regular micro-ruptures are reproduced and an interpretation of the origin of the regularity of the phenomenon as a stick-slip response is given. Synchronization of this periodic response in large, disordered, systems is exhibited by incorporating a global coupling in the model. The physical origin of this coupling is discussed. The simplicity of our model makes possible to provide master equations to describe the evolution state of the system as the inclination is increased. This analytical approach reveals an internal angle signing an intensification of the process.

The article is structured as follow: part II is devoted to the description of the numerical model. Part III presents the results obtained for three typical sets of parameters: a small system without global coupling, a large system without global coupling and finally a large system with global coupling. Part IV details a statistical approach to rationalize our observation. The last part discusses the understanding of the pre-avalanche behavior that can be obtained from our approach and compare our results to previous studies from the litterature.

The system under study is very close in spirit to the one-dimensional uniform Burridge-Knopoff model [13]. This seminal work has inspired numerous studies aiming at understanding the dynamics of earthquake faults using spring-blocks models [13][14][15]. Similar models are also studied to understand the onset of frictional sliding between two interfaces in tribological studies [16,17]. In our case we build on a model which has been previously proposed to study the effects of minute temperature changes on the stability of granular materials [18,19]. The main features of the present study compared to previous ones are the following: (i) we consider an homogeneous gravitational loading due to the progressive inclination of the blocks; (ii) we are interested in the evolution of the system from an initial preparation to the critical state which immediately preceed the avalanche; (iii) we propose an implementation of some long-range coupling in the system. We discuss the physical origin of this global coupling in section V.

We consider N identical frictional sliders of mass m lying on a rigid incline making the angle α with the horizontal (Fig. 1). The sliders are connected to one another by linear springs of stiffness k. In addition, each slider is connected to a rigid upper plate by a cantilever spring of stiffness k c . The upper plate insures a global coupling between the sliders. The springs as well as the upper plate are massless. In this configuration, the sliders are subjected to the elastic forces due to the springs, to their own weight mg, with g the acceleration due to gravity, and to the reaction force from the incline which includes a frictional force. The notations used in the following are shown in Fig. 1.

α k x 1 2 N-1 N g l * x x 1 x x 2 x x N 1 (c) 2 (c) (c) N c k FIG.
1. Sketch of the system under study.

One of the fundamental differences between the system under study and the usual Burridge-Knopoff (BK) models is found in the loading process. Indeed, we focus on the response of the system to a progressive increase of the inclination angle α, whereas the system is usually loaded by pulling it at one extremity. We prepare the system at horizontal, and then progressively increase the tilt in a quasi-static way until all the sliders descend the slope, which corresponds to the final avalanche. We are interested in the temporal and spatial distribution of the rearrangements that lead to the avalanche. In particular, we want to study how the distribution of the static frictional forces exerted on the sliders are modified from an initially uniform preparation because of the progressive loading.

Another difference is found in the modeling of the frictional contact. Whereas, in most of the BK-models velocity weakening of the friction force is introduced, we here characterize the frictional contact between the incline and the sliders by static and dynamic frictional coefficients. However, we consider that, due to the heterogeneity of the local properties of the incline surface, the static coefficient µ s,n takes a different value for each of the individual sliders, indexed by n (from 1 to N , Fig. 1). By constrast, because the sliders in motion see average properties of the incline surface, we consider a single value µ d of the dynamical frictional coefficient for all the sliders. Note that, in addition, we consider that µ d < µ s,n (∀n). Thus, in summary, we write that the slider n starts moving if

|f n+1→n + f n-1→n + f C→n + mg sin α| > µ s,n mg cos α (1)
where f n+1→n and f n-1→n are the elastic forces due to the neighbor sliders n -1 and n + 1, and f C→n the elastic force due to the cantilever spring that connects it to the upper plate. When the slider n is in motion with the velocity ẋn , the frictional force exerted by the incline is

f d,n = -µ d mg S( ẋn ) cos α (2)
where S denotes the sign function [S(u) = 1 if u > 0 and

S(u) = -1 if u < 0].

B. Dimensionless set of equations

The dynamics of the system is characterized by the timescale τ dyn = m/k. We can then introduce the dimensionless time t = t/τ dyn and position x = x/gτ 2 dyn . We then denote l the dimensionless natural length of the springs that link the sliders to one another and define ξ ≡ k c /k, the ratio of the stiffnesses of the two kinds of springs introduced previously.

Using these latter dimensionless variables, we can write the equations governing the dynamics of a slider in motion. We have (Eq. 2):

ẍ1 = -[x 1 -x2 + l] -ξ[x 1 - x(c) 1 ] -µ d S( ẋ1 ) cos α + sin α (n = 1) (3) ẍn = -[2x n -(x n+1 + xn-1 )] -ξ[x n -x(c) n ] -µ d S( ẋn ) cos α + sin α (n = 1, N ) (4) ẍN = -[x N -xN-1 -l] -ξ[x N - x(c) N ] -µ d S( ẋN ) cos α + sin α (n = N ) (5)
where

x(c) n = 1 N N i=1 xi + n - N + 1 2 l * (6)
with l * the dimensionless, constant, distance ( l * ≡ l * /gτ 2 dyn ) between the attachment positions x (c) n of the cantilever springs to the upper plate (Fig. 1, refer to next section addressing the preparation of the system). Eq. ( 6) reflects the fact that the upper plate moves with the barycenter of the sliders.

Depending on its position indexed by n, a slider, initially at rest, starts moving if the following condition is fulfilled (Eq. 1):

x2 -x1 -l -ξ x1 - x(c) 1 + sin α| > µ s,1 cos α (n = 1) (7) xn+1 + xn-1 -2x n -ξ xn -x(c) n + sin α| > µ s,n cos α (n = 1, N ) (8) xN-1 -xN + l -ξ xN - x(c) N + sin α| > µ s,N cos α (n = N ) (9)
C. Numerical method

In this section, we detail the numerical methods, starting with the crucial preparation of the initial state of the system.

Preparation of the system

First, the initial set of static frictional coefficients µ (0) s,n accounting for the contact betweent the slider n and the incline (angle α = 0, initially) are drawn at random from a Gaussian distribution:

p(µ s ) = 1 2πσ 2 µ exp - (µ s -μs ) 2 2σ 2 µ ( 10
)
where μs is the mean value and σ µ the width of the distribution. However, if a random value µ d exceeds µ s , we draw it again at random until we have µ d < µ s .

Once the initial set of static frictional coefficients µ

(0) s,n
is chosen, we set the initial positions of the sliders, x(0) n , insuring the mechanical stability of the initial configuration in absence of global coupling (The cantilever springs are deconnected, i.e. ξ = 0). To do so, a set of random tangential forces are drawn according to a uniform distribution in ] -μs , μs [. The corresponding positions, x(0) n , of the sliders are then computed. When considering the behavior of the system in absence of coupling (ξ = 0), the obtained configuration is our initial condition. When the behavior of the system with coupling (ξ = 0) is considered, we alter the configuration in the following way. First, we set l * ≡ (x (0) N -

x(0) 1 )/(N -1), the average distance between the sliders. Then, the coupling is introduced by connecting the cantilever springs to the upper plate (by setting ξ to a non-zero value). In response, some sliders can loose stability. If so, the dynamical equations (3) to [START_REF] Duranteau | AIP Conference Proceedings[END_REF] are integrated until a mechanically-stable configuration is reached. The final positions, x(0) n , of the sliders are then our initial condition.

Numerical integration

From the initial horizontal state (α = 0), the tilt angle α is progressively increased. From the knowledge of the static frictional coefficients µ (0) s,n and positions x (0) n , we can determine the value of the angle α (1) leading to a first rearrangement using Eqs. (7) to (9). The dynamical equations are then solved to obtain the set of new steady positions x

(1) n of the sliders that moved. To do so, we set the angle α to α (1) and integrate, using a standard velocity Verlet integrator (with a timestep ∆t τ dyn ), the dynamical equations (3) to ( 5) for all the sliders that enter in motion. At each integration step, we check if any of the motionless sliders is destabilized by any displacement of its neighbors. In case of such event, the integration step is reduced until the initiation of the motion of the corresponding slider. In the same way, we check if any slider in motion comes to rest and we adapt the integration step accordingly. A new value of the static frictional coefficient is drawn at random, according to Eq. ( 10), for those sliders that come to rest. We consider that the rearrangement ends when all the sliders have come back to rest. The new state of the system is then accounted for by the sets of new steady positions x From those new initial conditions the procedure can be iterated, leading to a set of destabilization angles α (k) associated to positions x (k) n and static frictional coefficients µ (k) s,n . The procedure ends when the angle α reaches a critical value, α c , such that all the sliders enter in motion and continuously accelerate downwards. This final angle is, by definition, the avalanche angle.

III. NUMERICAL RESULTS

In order to get insights in the behavior of the model, we report numerical results obtained for three typical systems, small and large systems with or without coupling: In all the simulations we set the parameters to typical values, i.e. l = 10, µ d = 0.55, μs = 0.6, and σ µ = 0.01.

A. Small system, no global coupling

In Fig. 2 we report results obtained for a small system of N = 10 sliders in absence of coupling (ξ = 0).

We observe in Fig. 2a that, far from being uniformly distributed neither between the sliders nor during the inclination process, the displacements are gathered in bursts of displacements of adjacent sliders at given angles separated by quiescent inclination intervals. Moreover, the interval between the inclinations angles α that lead to rearrangments seems to be almost constant. These features can be understood by considering the evolution of the tangential component of the force, due to the springs and the weight, on each of the sliders as function of α (Fig. 2b): Ft = fn+1→n + fn-1→n + fC→n +sin α. While α is increased, the component of the weight along the incline increases as sin α (dash-dotted line). Starting from the horizontal (α = 0), we observe that the first rearrangment occurs for α 0.1, when the most unstable slider looses stability. The displacement of this slider leads to a drop of the force it is submitted to. Simultaneously, the tangential forces exerted on its two direct neighbors present a sudden increase. At about α 0.2, the same block again looses stability, and drags downwards one of its neighbors which then also looses stability and moves. A further increase of the inclination by a few degrees leads to another rearrangment that involves three sliders, thus a larger region. Sliders, that were unperturbed up-to-now, are now adjacent to sliders that move, and the tangential forces they are subjected to exhibit sudden increases, which places them closer to the instability threshold. Hence, step by step, the bursts of reorganizations imply more and more sliders. Finally, a critical angle, α c 0.5, is reached at which all the sliders destabilize and accelerate down the incline. Note in Fig. 2b that the tangential forces exerted on the sliders that rearranged at least once are bounded by two well-defined values B. Large system, no global coupling For a large system, we observe the same typical behavior of the system with the main difference that bursts of rearrangments are observed at different places in the system (Fig. 3a). The picture is more complex as the rearrangments that occur in different regions are not synchronized. In addition, the intervals between the inclinations angles α are, on average, smaller than previously observed in a smaller system. Note also that, upon increasing inclination, the typical size of the active regions increases, which leads to their coalescence. Reporting the tangential components of the forces (Fig. 3b), we observe that rearrangements initiate at different locations as the instability threshold is reached by several sliders early in the inclination process. Upon further increase of the inclination angle, each of the active regions evolves, and grows in size, as previously described for the smaller system, which leads, finally, to the avalanche of the whole system.

C. Large system with global coupling

From now on, it is of particular interest to observe the behavior of the exact same system when a slight global coupling is introduced. To do so, we set the ratio of the spring constants k and k c to the small, but non zero, value ξ = 0.001.

Comparing the behaviour of the system with (Fig. 4) and without (Fig. 3) coupling, we immediately notice the coupling leads to larger rearrangments, involving more sliders and larger displacements than in the uncoupled system. In addition, in Fig. 4, we observe two main active regions involving the two extremities of the system. From an angle of about α ∼ 0.4, the upper region extends over more than half the size of the system. We also note also that the two active regions synchronize α 0.42 even if they do not overlap until the avalanche of the whole system which occurs for α 0.51. In Fig. 4b, we observe that the tangential components of the forces behave as previously observed in absence of coupling (Fig. 3b). Those components, for the sliders that moved at least once, are bounded by two welldefined values. The effect of the global coupling is clearly visible in the evolution of the tangential component of the forces for α 0.42: for each rearrangment, the tangential component corresponding to the sliders that are not involved exhibit a large increase. This synchronization is also visible by focussing on the behavior of the forces between the two dashed lines in Fig. 3b which correspond to the sliders that moved.

IV. ANALYTICAL MODEL

We present in this section a simple approach which allows to understand the main features that have been observed in the simulation. First we will study the case of a single slider surrounded by two immobile blocks (Sec. IV A). In Sec. IV B, we will see how a statistical model can be built from the results obtained with a single slider. In this section no global coupling is ever considered (ξ = 0). The effect of the coupling, as well as the comparison with experiments, will be the subject of the discussion (Sec. V).

A. A single slider

Let us first consider a single slider bonded by two springs to two immobile walls as sketched in Fig. 5.

We consider that the walls are located at x * -and x * + . Initially, the system is horizontal (α = 0) and the slider at stable position x (0) . For the sake of simplicity, we assume that static frictional coefficient is constant and equal to µ s . Indeed, in practice σ µ 2µ s , such that the associated width of the distribution of the static frictional force can be neglected when compared to the variation of the tangential force associated to the choice of x (0) .

From Eq. ( 1), we can determine the angle α (1) corresponding to the limit of stability of the slider:

F (0) t + sin α (1) = µ s cos α (1) , (11)
where F (0)

t = x * + + x * --2x (0)
Once α (1) is reached by tilting the system, for this simple case the dynamical Eqs. ( 2) can be solved analytically, which leads to:

∆x (1) ≡ x(1) -x(0) = (µ s -µ d ) cos α (1) (12)
The total tangential force exerted on the block in its new position is

F (1) t = µ s cos α (1) + ∆ F (1) , (13)
where ∆ F (1) = -2(µ s -µ d ) cos α (1) In the same way, upon further inclination, the slider destabilizes for α (2) verifying:

F (1)
t -sin α (1) + sin α (2) = µ s cos α (2) (14

)
Assuming that ∆α (2) = α (2) -α (1) is small, we get to the first order:

∆α (2) = 2(µ s -µ d ) 1 + µ s tan α (1) (15)
which remains small provided that µ s -µ d 1. The reasoning can be reproduced for any of the rearrangments, leading to:

∆α (i+1) = 2(µ s -µ d ) 1 + µ s tan α (i) (16)
Thus, a generic behavior emerges as shown in Fig. 5. We observe, first, a loading phase, which depends on the initial position x(0) . The tangential component of the force Ft increases continuously until it reaches the instability threshold for α = α (1) and suddenly drops by ∆ F (1) t . Then, rearrangements occur for successive values α (i) of the inclination angle. They are marked by sudden drops ∆ F (i)

t
of the tangential component, followed by continuous increase towards the next threshold. In a first order approximation in α, the force drops ∆ F (i) do not depend on the initial position of the slider. In addition, for small enough angle α (i) , the interval ∆α (i+1) is almost constant and equal to ∆α 2 (µ s -µ d ). Such periodic behaviour is reminiscent of the well-known stickslip motion. Finally, note that a sliding event leading to a drop ∆ F (i) t of the force exerted to the slider, results in an increase by ∆ F (i) t /2 = (µ s -µ d ) cos α of the force exerted on each of the walls.

B. Statistical point of view

As will be seen in the discussion (Sec. V), the simple analysis of a single slider between two walls helps to understand most of the features observed for a small system. For large systems however, where the pattern that emerges is more complicated (Sec. III B), a statistical approach is more appropriate to model the system.

In the present section, we seek for master equations describing the evolution of the distribution of the tangential forces Ft in a system as a function of the inclination angle α. We will make strong approximations to be able to provide an analytical model, but we will compare the solutions we obtain to numerical simulations of large systems in absence of global coupling (ξ = 0).

Qualitative analysis

We report in Fig. 6 the probability distribution function of the tangential forces, Ft , exerted on the sliders at different values of the inclination angle α, in absence of global coupling (ξ = 0).

Initially, at α = 0 (dark blue curve in Fig. 6), the distribution is flat in the interval ] -µ s , µ s [, which corresponds to the initial, random, preparation of the system at horizontal (Sec. II C). Upon inclination, we observe that the distribution functions exhibits two main parts:

-A rather flat part for Ft 0.4 with a well-defined plateau at the constant initial value. The lowest values of Ft are slowly depleted.

-A peak, for 0.4 Ft μs , that grows rapidly.

The qualitative behaviour of the distribution function can be understood as follows. Upon inclination, the lowest values of Ft are depleted because of the increase of the projection of the weight along the incline. By this simple effect, the lower bound of the distribution increases according to -μ s cos α + sin α. Conversely, upon increasing inclination, sliders loose stability. The upper bound, which equals μs cos α according to (Eq. 1), slightly decreases, but this variation can be neglected at the first order. The value of Ft for the sliders that moved remains close to the upper bound. The peak grows as their number increases. This picture holds as long as the plateau and the peak are well separated, i.e. for α 0.25. Then, the contribution of the sliders that are involved in the rearrangements become predominent.

In order to get more insights in the contribution of the rearrangments, we focus on the distribution of the values of Ft for only the sliders that moved (Figure 7). Defining F * t ≡ (2µ d -μs ) cos α, we observe that, for small inclination, the distribution function is peaked around Ft / F * t = 1. This observation shows that the rearrangements initially involve single sliders, the width of the peak is directly linked to the width of the distribution of µ s . By contrast, for large inclination, the distribution function exhibits a peak close to the value Ft / F * t 1.1, which is an upper limit for all the distribution functions. We can understand the accumulation of the tangential forces at this limit value in the following way. First, let us mention that Ft / F * t 1.1 corresponds to Ft = µ d cos α. For rearrangements implying several interacting sliders, the movement becomes complex enough for the vanishing velocities not to be determined by a sinusoidal movement, but by the overall slowing down of the dynamics due to dissipation. At the limit of vanishing accelerations, the sum of elastic forces and weight on each block compensates the dynamical friction so that, when the block stops, the value of the resultant tangential force is typically less than µ d cos α, thus controlled by the dynamical frictional coefficient.

In conclusion, after a sliding event, for the sliders implied in the event, the tangential forces Ft drop by a quantity comprised between ∆µ cos α and 2∆µ cos α with ∆µ = μs -µ d . This is the main ingredient of the model we present in the next section IV B 2.

Model

Here we present a model based on the prior observations. We consider in Fig. 8 a schematic representation as a bar histogram of the distribution function reported in Fig. 6. To do so, we discretize the probability density function in bins of width 2∆µ with ∆µ = μs -µ d . In addition, in a first approach, we neglect the effect of the inclination on the upper bound, which we assume to equal μs . We denote ν the number of bins of width 2∆µ, starting from the upper bound, μs , and including the lower bound, sin α -μs . Thus, ν is a function of α that decreases with α as: ν(α) = 2μs-sin α 2∆µ . For convenience, the bins are numbered from large to small forces such that the bin No. 1 denotes the peak on the righthand-side, corresponding to the peak of the distribution. We denote P i the height of bin No. i. With these definitions, the probability for a slider to be submitted to a tangential force belonging to the bin No. i is given by 2∆µP i . As a result of an angular increment δα, all the tangential forces increase because of the increase of the tangential component of the weight which goes like sin α. Consequently, in the distribution, N cos αP i δα sliders shift from the bin i to the bin i -1, excepted for the bin No. 1. This process is sketched by the red solid arrows in Fig. 8. For the bin No. 1, N cos αP 1 δα sliders loose stability. From the qualitative description of the system (Sec. IV B 1), we know that the final value of Ft after the rearrangement lies mostly between µ d cos α and (µ s -2∆µ) cos α. Thus, after the rearrangments, the sliders remain in bin No. 1. This process is sketched by the blue dashed arrow in Fig. 8. This simple process can be put in equations in the form:

dP ν dα = - 1 2∆µ P ν dP i dα = 1 2∆µ (P i+1 -P i ) (i = 1, ν) (17)
dP 1 dα = 1 2∆µ P 2
Starting from the condition that the distribution is uni-form for α = 0, thus from P i (0) = 1 2μs , we get:

P ν (α) = 1 2μ s e -α 2∆µ P i (α) = 1 2μ s ν-i k=0 1 k! α 2∆µ k e -α 2∆µ (18) P 1 (α) = 1 2μ s    ν -1 - ν-2 k=0   k j=0 1 j! α 2∆µ j   e -α 2∆µ   
It is of particular interest to focus on the behavior of the solution far from the lower boundary. For realistic values of µ s and µ d , ν(α) remains large and, in the limit, ν -i 1, we have ν-i k=0

x k k! e x (∀x). Consequently, for small angle α so that ν(α) is still large (ν(0) = μs ∆µ 1 for realistic values of µ s and µ d ) and for i ν, i.e. far from the lower bound of the system, we have:

P i (α) 1 2μ s (i = 1) ∧ (i ν) P 1 (α) 1 2μ s 1 + α 2∆µ (19)
Note that, as is, P 1 (α) accounts for the amplitude of the peak at large Ft as function of the inclination whereas the constant P i (α) correspond to the plateau.

We can now compare our model (Eq. 18) and its limiting case (Eq. 19) to the results of the numerical simulations. To do so, we report in Fig. 9, the probability P i as function of the inclination angle α for the data already reported in Fig. 6.

First, we observe that Eq. ( 19) correctly describes the linear increase of the amplitude of the peak P 1 at large Ft in a large range of inclination angle α (0 to typically 0.38). Second, Eq. ( 19) also accounts for the constant value of P i in the same range of inclination, but only for the largest values of Ft . Slightly less crude, Eq. ( 18) predicts the decrease of the P i (i = 1) upon increasing α. However, reporting the probability of the smallest value of Ft , we observe the model underestimates the depletion of the smallest values of Ft . This point will be further commented in the discussion (Sec. V).

In spite of this reasonably good agreement between the model and the numerical results, a trend observed in Fig. 9, which is not described at all by the model, is the linear increase of P 2 starting at α 0.2, followed by a sharp increase for α 0.4. We observe that, for 0.2 α 0.4, the slope dP 2 /dα = 1 4µs∆µ , thus takes the same value that dP 1 /dα. The model fails in accounting for the increase of P 2 because we did not consider the decrease of the upper bound upon inclination (the maximum force value evolves like cos α).

As a matter of fact, due to the shift of the upper bound, the bin No. 2 takes the place of bin No. 1 as the inclination proceeds.

Up-to-now, we commented mainly on the behavior of the system for limited values of the inclination angle 9. Height Pi vs. angle α for the data of Fig. 6. The dashed straight lines correspond to prediction of the crude model (Eq. 19). The exponential decay of Pν (α) is given in Eq. ( 18). The dotted line is a guide for the eye. The vertical solid line is the threshold from Eq. ( 21 α 0.38. For α 0.38, a drastic change in the behavior of all the probabilies P i is observed. Indeed, the probability of the largest values, i.e. P 1 + P 2 , drastically increases whereas the probabilities of smaller values P i (i > 2) decreases. One interesting question is whether our model is able to predict the typical angle, α t 0.4, at the transition. To answer the question, consider the probability P 3 of the force just below the peak:

P 3 (α) = 1 2µ s   ν(α)-3 k=0 1 k! α 2∆µ k   e -α 2∆µ 1 2µ s 1 - 1 (ν(α) -2)! α 2∆µ ν(α)-2
We previously considered that ν(α) was large enough for the term on the right-hand-side in the parenthesis to be negligible. Let assume that this is not the case anymore and that this term is of 1/10, i.e. that:

1 (ν(α) -2)! α 2∆µ ν(α)-2 = 0.1. ( 20
)
Using the Stirling approximation and the definition of ν(α), we get:

α t 2µ s 1 + e 0.54 µ s (21)
which, we point out, does not significantly depend on the choice of the value 1/10 in Eq. (20). We thus get α t ∼ 0.32 for μs = 0.6, which underestimates the observed value. However, in regard of the crude approximations we made, we can consider that this last estimate is reasonable.

V. DISCUSSION

A. Discussion of the model

The very simple model of a single slider (Sec. IV A) provides a good understanding of the occurence of the, local, quasi-periodic dynamics since the beginning of the inclination process, i.e. for small tilt angle α (typically ≤ 0.2). Even in a large system, the first sliders that destabilize are isolated from one another and behave like the single slider between two walls, which consist of the two immobile neighbors. Upon inclination, rearrangments involving initially only one slider lead to an increase of the tangential force applied to its neighbors sufficient to destabilize at least one of them. The pattern is then more complex but the picture of a process consisting of a series of, quasiperiodic, localized destabilizations holds. The number of sliders involved in each rearrangment increases upon inclination. When several sliders are involved, the rearrangements repeat with the period ∆α ∼ (μ s -µ d ) rather than ∆α ∼ 2(μ s -µ d ) (Fig 2). Indeed, as discussed in Sec. IV B 1, when several sliders enter in motion, the progressive slowing down of the dynamics due to dissipation leads to the value µ d cos α [instead of (2µ d -μs ) cos α for a single slider] of the tangential force.

However, in large systems, in absence of global coupling, the previous picture holds only locally and the overall dynamics is not quasi-periodic (Fig. 3). The rearrangments occuring in different regions are not correlated. The growth in the size of the active regions upon inclination leads them to merge before the whole system looses stability.

The introduction of a global coupling, even small (ξ = 10 -3 ), leads to the synchronisation of the previously isolated local rearrangements and, thus, to the occurence of a quasi-periodic dynamics of the whole system (Fig. 4). This synchronization is reminiscent of the synchronization of oscillators, which is a well-known phenomenon [20].

A particularly interesting, and striking, result is that one can distinguish a clear change in the activity of the system, with or without global coupling, at an inclination angle α of about 0.25, thus typically half of the avalanche angle α c : for α 0.25, the system is rather quiescent with few rearrangements occuring locally whereas, for larger α, numerous large events of increasing size occur. The model shows that a change in the dynamics is indeed expected for α t 0.54 μs . For 0 ≤ α α t , the number of active sliders increases linearly with α. By contrast, for α α t , the number of active sliders drastically increases (Fig. 9). In practice, 0.54μ s is close to half the angle of avalanche angle, α c , so that we can consider that, to the level of approximation of our model, the transition is well accounted for.

B. Comparison to experiments

Finally, we underline strong similarities between features of our model and the experiments mentionned in the introduction [2][3][4].

Upon inclination of a container filled with a granular material, one observed, for small tilt angle, the occurence of small, uncorrelated, rearrangements randomly distributed in the system. At approximately half the avalanche angle, the activity drastically increases in the form of bursts of synchronized displacements.

The order of magnitude of the angular period observed experimentally is typically of a few degrees. In our model, we predict that the periodicity is linked to the difference between the static and the dynamic friction coefficients: ∆α ∼ (µ s -µ d ). For granular matter, those friction coefficients can be assimilated to the tangent of, respectively, the angle of avalanche, at which the granular surface looses stability, and of the angle of repose, at which corresponds the granular surface stops flowing. In granular matter (µ s -µ d ) ∼ 0.1 [21] leading thus to δα indeed of a few degrees.

In our model, we obtain synchronization only when some global coupling is introduced, which leads to the question of the physical significance of this global coupling in experiments. In the classical Burridge-Knopoff model a rearrangement affects only the immediate neighbors. Such local model is irrealistic for modeling faults as the bulk materials at each side of the fault are of finite stiffness and mediate long-range elastic forces. Taking into account such long-range elastic redistribution of the stress after a reaarrengement has been taken into account in Burridge-Knopoff models by coupling blocks elastically to a various number of distant neighbors [22]. Such long-range elastic coupling has also been discussed in tribology-related studies (see [23] and references therein). In granular matter, it has been shown that when a local, plastic rearrangement takes place, an elastic response with long-ranged redistribution of the stress is observed [24,25]. This justifies the introduction of a longrange coupling in the model. Another clue is provided by the observation, in experiments, that the pseudo-period of the bursts depends on the cohesion in the system [4,[START_REF] Duranteau | AIP Conference Proceedings[END_REF]. This is in accordance with a picture of a synchronization due to a long-range coupling.

It has to be underlined that our model is not in contradiction with previous works, in particular with studies based on cellular automata reported in Ref. [3]. The authors reported an exponential growth of the accumulated plastic activity upon inclination. This trend, also observed experimentally, has been considered as a test validating the model. Such an exponential evolution could seem, at first sight, to be in contradiction with the linear evolution of the number of active sliders which arises in our model. However, the contradiction is only apparent. Actually, reporting the cumulated displacements A(α) in the whole system from the beginning of the inclination process:

A(α) = ic i=0 N n=1 x(i+1) n -x(i) n , with
x(k) n the position of the n th slider after the i th rearrangement and i c the total number of rearrangements we also observe a typically exponential increase for our model (Fig. 10a). This is explained by an exponential increase of the typical displacement associated to the rearrangements with the inclination angle α (Fig. 10b).

VI. CONCLUSION

We have presented a one-dimensional frictional model consisting of a chain of elastically coupled sliders in frictional contact with an incline. The model reproduces most of the features observed experimentally when quasistatically tilting a box filled with a granular material. In particular, the system reproduces the quasi-periodic series of rearrangments observed in experiments.

We explain the regularity of the bursts of rearrangments as the result of a combination of stick-slip and synchronization due to a large scale coupling. The phenomenon is thus related to the problem of synchronization of oscillators. The coupling originates from the redistribution of the forces, at large scale, in the bulk of the material.

In addition, our toy model provides a statistical approach to describe the evolution of the state of the system which approaches the avalanche. In particular, we revealed an internal threshold before the avalanche occurs. This threshold delimits the transition, at typically half the avalanche angle, from an initial smooth increase of the number of active sliders to an accelerating regime where a dramatic increase of the activity takes place. In regard to the prediction of avalanches, the determination of such a threshold is precious as it delimits a rather regular and predictible evolution of the system and a rapid growth regime which announces the avalanche.

Our model remains simplistic and, in accordance, provides no more than a qualitative description of the phenomenon. In addition, it has to be noted that we considered a unique preparation, consisting of a uniform distribution of the initial tangential forces. Nevertheless, considering that the model is based on very general arguments (disorder, solid friction and coupling), we expect that it unveils the core mechanisms at play in different experimental configurations.

  A. -N = 10 sliders without coupling (ξ = 0). B. -N = 400 sliders without coupling (ξ = 0). C. -N = 400 sliders with coupling (ξ = 0).

FIG. 2 .

 2 FIG. 2. Behavior of a small system without global coupling upon inclination -(a) Hot map of the displacements ∆xi vs. angle α. Each horizontal line corresponds to one slider. The darker the points are, the larger is the displacement of the slider. The vertical dashed line corresponds to the avalanche angle and ∆µ = μs -µ d . (b) Tangential force Ft ( Ft = fn+1→n + fn-1→n + fC→n + sin α) vs. angle α (The slider number is indicated by the color bar). The upper dashed line is the stability limit: μs cos α. The lower dashed line is the stability limit (2µ d -μs) cos α. The dash-dotted line: sin α [N = 10, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

FIG. 3 .

 3 FIG. 3. Behavior of a large system without global coupling upon inclination -(a) Hot map of the displacements ∆xi vs. angle α. (b) Tangential force Ft vs. angle α (The slider number is indicated by the color bar). The upper dashed line is the stability limit: μs cos α. The lower dashed line is the arrest limit (2µ d -μs) cos α [N = 400, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

FIG. 4 .

 4 FIG. 4. Behavior of a large system with global coupling upon inclination -(a) Hot map of the displacements ∆xi vs. angle α. (b) Tangential force Ft vs. angle α (The slider number is indicated by the color bar). The upper dashed line is the stability limit: μs cos α. The lower dashed line is the arrest limit (2µ d -μs) cos α [N = 10, ξ = 0.001, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

FIG. 5 .

 5 FIG. 5. Tangential force Ft vs. angle α -The upper dashed line is the stability limit: μs cos α. The lower dashed line is the arrest limit (2µ d -μs) cos α. Inset: sketch of the system under study [µ d = 0.55 and µ d = 0.6].

FIG. 6 .

 6 FIG. 6. Probability distribution function of the tangential force Ft upon inclination -at different values of the inclination angle -Each distribution corresponds to an average, at a given α (indicated by the colorscale), over 20 numerical runs. Inset: Enlargment of the distribution for Ft < 0.4 [N = 100, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

FIG. 7 .

 7 FIG. 7. Probability distribution function of the ratio Ft/ F * t for the sliders that moved upon inclination -The data are from the 20 same numerical runs as in Fig. 6 [N = 100, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

FIG. 8 .

 8 FIG. 8. Schematic representation of the histogram equivalent to the probability distribution function of the tangential force, Ft at a given inclination angle α. The arrows indicate the main modification of the distribution upon further inclination by δα (see Sec. IV B 2).

  ) [N = 100, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].

  The tangential force applied to it increases until it reaches the threshold μs cos α [upper dashed line in Fig 2(b)] and then drops by 2(μ s -µ d ) cos α, then reaching the value (2µ d -μs ) cos α [lower dashed line in Fig 2(b)]. Each rearrangment leads to the increase of the tangential force applied to the neighbors by (μ s -µ d ) cos α.

FIG. 10 .

 10 FIG. 10. (a) Cumulated displacement A vs. angle α -Upon inclination, A increases exponentially in spite of the linear increase of the number of sliders involved in the rearrangments with α. (b) Average displacement during a rearrangement ∆x vs. angle α. [N = 100, ξ = 0, l = 10, µ d = 0.55, μs = 0.6, and σµ = 0.01].