The Olami–Feder–Christensen earthquake model in one dimension

We study the earthquake model by Olami, Feder and Christensen in one dimension. While the size distribution of earthquakes resembles a power law for systems with small sizes, it splits for systems with large sizes into two parts, one comprising small avalanches and showing a size-independent cutoff, and the other comprising avalanches of the order of the system size. We identify four different types of attractors of the dynamics of the system, which already exist for very small systems. For larger system sizes, these attractors contain large synchronized regions.


Introduction
The Olami-Feder-Christensen earthquake model [1] is probably the most studied nonconservative and supposedly self-organized critical model. Nevertheless, the origin of its power-law like avalanche-size distribution is still not clear. Apart from these power laws, the model shows a variety of other interesting and unusual features such as a marginal synchronization of neighbouring sites driven by the open boundary conditions [2], and the violation of finite-size scaling [3,4] together with a qualitative difference between system-wide earthquakes and smaller earthquakes [5]. Also, small changes in the model rules (like replacing open boundary conditions with periodic boundary conditions [6], or introducing frozen noise [7]), destroy the power laws. Recently, it was found that the results of computer simulations are strongly affected by the computing precision [8], and that the model exhibits sequences of foreshocks and aftershocks [9].
In order to better understand the model, we study here its one-dimensional version. The model is highly nontrivial even in one dimension, and some of its properties resemble those in two dimensions. Just as in two dimensions, we find large synchronized regions and a fundamental difference between the avalanches triggered at the boundaries and those triggered deep inside the system. We identify different types of attractors of the dynamics of the system and explain the features of the model in terms of the properties of these attractors. Our main finding is that the system in the stationary state can be separated into a boundary region, where all larger avalanches are triggered, and one (or two) synchronized inner region(s), the size of which can be varied without changing the behaviour of the boundary region.
The outline of this paper is as follows. In the next section, we introduce the model rules. In section 3, we focus on a system of up to four lattice sites and find its attractors. In section 4, we view the model from a dynamical systems' perspective and present a general analytical approach that allows us to classify the attractors into four different types. In section 5, we study larger systems. First, we investigate the approach to the stationary state as function of the system size and the model parameter. Then, we discuss the properties of the stationary state. Finally, we summarize and discuss our main findings in section 6.

The model
The Olami-Feder-Christensen model is a discretized and simplified version of the Burridge-Knopoff model of earthquakes [10]. In a one-dimensional system consisting of L sites, it is defined by the following rules: at each site i = 1, . . . , L, a continuous variable z i is defined that represents a local force. The force at all sites increases uniformly at a constant rate, which we set equal to 1. When the force z i exceeds the threshold value z c , which can be chosen to be z c = 1 without loss of generality, the force at this site is reset to zero, while the two nearest neighbours (or the only neighbour, if the toppling site is a boundary site) receive a force increment of αz i . The parameter α is the only parameter of the model, and it has a value in the interval (0, 0.5). If a neighbour is lifted above the threshold, the force on its neighbours is immediately increased according to the same rule, etc, until the 'avalanche' (the earthquake) is finished. The 'size' of an avalanche, s, is defined to be the number of toppling events during this avalanche. Such an earthquake is instantaneous on the time scale of driving. After the earthquake, the force is again increased with unit rate, until the next site reaches the threshold, triggering the next earthquake, and so on. The total force increment received by the system since the beginning defines the time. However, for the computer simulations it is often more useful to measure time as the number of topplings per site, which is the total number of topplings divided by the system size. We will use both definitions of time, and will always clearly indicate which one is used at a given place.
This model is deterministic, once the initial conditions are given. Usually, the initial conditions are chosen randomly from a uniform probability distribution for each site. Since the model is deterministic and dissipative, it has attractors of the dynamics.
From a dynamical systems' perspective, the model can be viewed as an L-dimensional map, which maps the state of the model after one avalanche (which may have size 1) on the state after the next avalanche. Due to the toppling, the map is noncontinuous.
If α > α c = ( √ 3 − 1)/2 0.366, a site that topples can in principle receive from its neighbours packages of a total size larger than 1, causing the first site to topple again. Throughout this paper, we assume that each site topples only once during an avalanche, and we limit therefore our numerical studies and analytical arguments to the case α < α c , except where indicated otherwise.
Let us briefly summarize a few known results that are relevant for our study of the onedimensional system. Firstly, it is found that for periodic boundary conditions the dynamics approach an attractor where every site gives to its neighbours only force packages of size α. This means that no force value exceeds the threshold value z c on the attractor. After a time 1 − 2α, each site has toppled once and has received two force packages of size α from its neighbours. This means that after time 1 − 2α the force on each site is again the same. Slightly increasing or decreasing force values gives again a periodic orbit with period 1 − 2α, as long as this change does not cause a toppling site to lift its neighbour above the threshold.
Secondly, the open boundary conditions are responsible for the occurrence of large avalanches and large synchronized regions, where neighbouring sites differ in force values only by a small amount. A nice explanation for this has been suggested by Middleton and Tang ten years ago [2]. They considered a system of two sites, where one site is driven at a slower rate than the other. This mimics the fact that sites close to the boundary receive on an average less force packages than those deep inside the system. The two sites settle on an attractor where the slower site always topples first and lifts the faster site above the threshold. The faster site therefore loses more force when it topples, and the slower site receives a larger force package. This compensates the different driving speed, and the two sites remain synchronized and always topple together.
Thirdly, the largest possible force package that a site can pass to its neighbour is α/(1 − α). This package size is reached if an avalanche passes through a region where all sites are at the threshold.

L = 2
For L = 2, any state where the force difference between the two neighbours is larger than α is part of a cycle of period 1 − α. After a time 1 − α, each site has toppled once and has received a force package of size α due to the toppling of the neighbour. Furthermore, it has received a force increment 1 − α due to the uniform driving. The net change in force is therefore zero.

L = 3
For L = 3, we consider the state of the system whenever the middle site has toppled, i.e., when z 2 = 0. z 1 and z 3 then take values in [α, 1 + α + α 2 ] only. The lower value α is realized, if a site just toppled itself and lifted z 2 exactly to the threshold. The upper limit occurs only if all sites were at the threshold before.
Without loss of generality, we assume z 1 z 3 . Different regions in the z 1 -z 2 plane can be characterized by their sequence of topplings t i and of growth, until z 2 = 0 again. An example of this would be (t 3 , g, t 1 , t 2 ), for which first the right site topples, followed by a uniform growth, and then the toppling of the left site lifts z 2 over the threshold. There is a total of 14 different such regions, seven of which are marked by letters a to g in figure 1.
By investigating the transitions between these regions, one finds two attractive fixed point lines. One fixed point line is at z 1 = z * = α(1 + 2α)/(1 + α) with z 3 ∈ (z * , 1]. The other fixed point line is obtained by interchanging z 1 and z 3 . The corresponding attractor is a cycle of two avalanches, written as (g, t j , g, t i , t 2 ) in the above notation, with i being the site with z i = z * i , and j being the other neighbour of site 2.
A special case is the symmetric case z 1 = z 3 = z * . Sites 1 and 3 topple simultaneously, thereby lifting site 2 above the threshold, and receiving a package z * in turn.
To summarize, a system with L = 3 approaches a periodic attractor with two different avalanches.

L = 4
For the system size L = 4, all attractors of the dynamics are periodic. We find a variety of different attractors for a given value of α. Figure 2 shows the period of the attractors found doing a scan of the same 128 4 different initial conditions for each value of α and for two different precisions. The period is measured once in terms of the number of topplings, # t , and once in terms of the number of avalanches, # a . One can discern the following features 1. Degeneracy: attractors with different numbers of avalanches per period have the same total number of topplings per period. We found that different attractors can have a different toppling sequence, while the force packages that each site receives from its neighbours are identical. One explanation for this is that there are mechanisms in the model that create a degeneracy of different sites, for instance when two sites have been the end points of the same avalanche (after which they both had zero force) and remain synchronized until they reach the threshold again. The toppling sequence then depends on the exact implementation of the algorithm and on rounding errors due to finite precision; one or the other version of the attractor can be obtained depending on the initial conditions. The total number of topplings and the package sizes can be identical on the two versions of the attractor if there is a site between the two synchronized sites that topples only after having received a package from both sides. 2. Persistence: attractors exist over a certain interval of α values. If α is changed slowly (such that the system can follow adiabatically), the attractor remains the same as long as the avalanches remain the same. Eventually, a point will be reached where an avalanche decays into two avalanches (because the size of a package is no more large enough to lift the neighbour over the threshold), or where two avalanches merge to form one avalanche (because the distance between two neighbours has become smaller than the package size). As we will see in the next section, the stability of the attractor can change at this point (but typically not before). If the state with the new toppling pattern is stable, we have a new attractor with the same period (if measured in number of topplings). Otherwise, if the new state is unstable, the system moves to a different attractor, and we obtain a step in figure 2. 3. Divergence of the periods at small α: at small α, the period of the attractors increases. This can be explained by considering a boundary site and its neighbour. During each time interval 1 − 2α, the boundary site receives a package of the order of α from its neighbour, while the neighbour receives two packages of the order of α. It therefore takes of the order of 1/α time intervals until the two sites have again roughly the same height. For this reason, the period diverges as 1/α. 4. Smallest possible period: for all values of α, the shortest attractor has four topplings. In configuration space, this corresponds to the state (z,z, 0, α + α 2 ), with an arbitrary force z ∈ [α + α 2 , 1], and with the toppling sequence g, t 1 , t 2 , (g), t 4 , t 3 . Whether this attractor is realized, depends on the implementation of the algorithm. For smaller α and for smaller precision, it occurs more often and eventually has the weight unity. 5. Vastly different periods: for a given α, there exist attractors with widely different periods.
The most prominent periods lie in two bands, which are clearly visible in the figure 2. For certain values of α, very large periods occur with a considerable weight. Attractors with these large periods typically have toppling sequences that are most of the time periodic with a much shorter period than that of the attractor, but the force values do not have this short period. 6. Sensitivity to the computing precision: attractors with larger periods occur (for the same value of α) less often when the computing precision is smaller. The reason for this is that on longer attractors, there are more states, and states that are close to each other are more likely to occur. In a simulation with smaller numerical precision, such states can become identical, and the period of the attractor becomes shorter.

Analytical approach and classification of attractors
We describe the state of the system as a difference vector The uniform increase in force does not change this vector, but toppling sites do change it. The force value of the first site toppling in an avalanche is decreased by one, while its two neighbours receive a package α. This process can be decribed by adding a vector to x (with the four nonzero elements at the appropriate places). In contrast, a subsequent toppling event can be described by applying to x t a matrix that is identical to the unity matrix everywhere except for a 4 × 4 block on the diagonal. There are two different matrices, corresponding to avalanches propagating to the right and to the left. We assume that site i just toppled. Then, the next toppling event is represented by the matrix if the avalanche is moving to the left. The nontrivial column (number i − 1) describes the toppling of site (i − 1), which was lifted above the threshold by the prior toppling of site i.
If the avalanche is moving to the right, the corresponding matrix reads with the nontrivial column being number i. If site i − 1 and site i + 1 both topple simultaneously, the matrix contains both nontrivial columns. If the toppling site is a boundary site or the site next to it, surplus columns and rows of the nontrivial block of M i l or M i r are to be removed. We now focus our interest on the difference δ x between two systems and assume that the two systems have the same toppling sequence. This will be the case as long as they are sufficiently close to each other. Then they will both be updated by adding the same vectors and multiplying the same matrices in the same order. Adding the same vector to the state x of both systems has no effect on the difference δ x. The difference between the states of the two systems will therefore evolve solely by multiplying matrices of the form M i l and M i r to it (except if the first toppling site lifts both neighbours above the threshold; in this case the first matrix associated with this avalanche contains two nontrivial columns, while the other matrices have the usual form, since the two branches of the avalanche commute after the first toppling).
Let one of the two systems be on a periodic orbit. Whether the other system will approach the orbit, depends on the largest eigenvalue of the product with P being the total number of matrices occurring during one period, i(p) being the nontrivial column index for matrix number p, and ν(p) being l or r depending on whether an avalanche is moving to the left or to the right. If the largest eigenvalue of S is larger than one, the orbit is unstable and cannot be an attractor of the dynamics. If a site is lifted exactly to the threshold, the next toppling can be either considered as being part of the avalanche, or it can be viewed as an isolated toppling that is not accompanied by a matrix M. For these two cases, there are two different matrices S, which may have different eigenvalues. This means that the stability with respect to perturbations that lift the site slightly above the threshold can be different from the stability with respect to perturbations that lift the site slightly below the threshold.
Below, we describe four types of attractors that we have found in the model. All these attractors occur already in a small system of L = 4, but are also seen in large systems. We obtained these attractors and their properties using a combination of different methods. The matrix S is obtained by starting with the unit matrix and multiplying the appropropriate matrices M i l or M i r one after the other. The effect of M i l (or M i r ) on a matrix is that row number i (or row number i − 1), multiplied by a certain factor ±α or (1 + α), is added to three neighbouring rows and is itself multiplied by −α. For very small systems, we calculated S and its eigenvalues analytically. For larger systems, we evaluated during the computer simulations the position of the 1s in the matrix S. This means that we performed an expansion in powers of α and kept only the elements of order α 0 . In order α 0 , the matrices M i l and M i r simply add one row to a neighbouring one and replace the original row with all zeros. If a boundary site is caused to topple by its neighbour, all 1s that have been in the boundary row are flushed out of the system. Starting with a unit matrix, all 1s that remain in the system are in the same row after multiplying a sufficient number of M i ν matrices. We have never seen an attractor where this does not happen. Any given site of the system is reached by an avalanche that starts near the boundary, and therefore there cannot be 1s left in different rows. In order to estimate the eigenvalues for large systems, we observed the approach to the stationary state starting from a neighbouring configuration, and plotting the force value of a given site after each period.

Marginally stable attractors
Marginally stable attractors occur if the largest eigenvalue of S is identical to 1. We observed attractors regularly, where a column i of S is identical to the unit vector e i . This means that e i is an eigenvector of S, and adding a small multiple e i to the periodic orbit gives again a periodic orbit with the same toppling sequence. In terms of the matrix S this means that all 1s are in row i, where they stay forever. In terms of forces z j this means that increasing or decreasing all force values z j with j i by a small amount results again in a periodic orbit. The product (3) contains no matrix M i l or M i+1 r . Sites i + 1 and i never cause each other to topple. There is no avalanche that includes simultaneously site i + 1 and site i. We say that there is a 'barrier' between sites i + 1 and i. We found that the total size of the force packages that site i + 1 gives to site i during one period is identical to the total size of the force packages that site i gives to site i + 1. For L = 4, the barrier is always in the middle of the system. For larger sizes, it need not be in the middle, but it is often found at the centre of a system, since the synchronization proceeds at constant speed from the boundaries (see below).
Sites i and i + 1 to the right and left of the barrier must topple equally often during one period. If this was not the case and if site i toppled more often, there would be an instance where site i topples twice without site i + 1 toppling in between. After site i has toppled for the first time, its force is zero, and that of site i + 1 is at least as large as α. In order for site i to reach the threshold before site i + 1, it must receive a package from its left neighbour that is larger than α, while site i + 1 receives no package. The largest possible package size is α/(1 − α), and therefore site i + 1 has at least the force 1 − α 2 /(1 − α) > 1 − α at the moment where site i reaches the threshold for the second time. This means that site i + 1 is lifted above the threshold by the toppling of site i, in contradiction to our assumption that there is a barrier between the two sites. Therefore, the two sites must topple equally often.
There exists no attractor with two or more barriers. We show this in two steps. First, let us assume that all force packages passed over the barriers are of size α. Then the region between the two barriers is like a system with periodic boundary conditions, and no package passed on within this region is larger than α. The two sites immediately outside the barriers must not be lifted above the threshold by their neighbours. Otherwise, they would pass packages larger than α over the barrier. Furthermore, the two sites immediately outside the barriers must not topple more often than the sites between the barriers. Therefore, they cannot receive packages larger than α from their outward neighbours. This means that the sites immediately outside the barriers, in fact, also belong to the domain of sites that are never lifted above the threshold and that always receive packages of size α. By repeating this argument, we find that no site in the system can be lifted above the threshold. However, this situation cannot be realized with open boundary conditions. It occurs for periodic boundary conditions. Now, since we have ruled out the possibility that the region between the two barriers receives only force packages of size α from outside, let us assume next that they receive on average packages larger than α. We simulated the region between the two barriers by inserting packages of size larger than α at its boundaries immediately after the boundary sites have toppled. This leads to attractors where avalanches are triggered at the centre of the system and are running outwards. The attractors and therefore the number of topplings per unit time of the boundary site are determined by the size of the region and the size of the packages received from outside. On the other hand, this number of topplings must be identical to the number of topplings of the site on the other side of the barrier in the original system. However, there is no free continuous parameter left to match this condition, and therefore it usually cannot be satisfied. This problem does not arise in the case of a single barrier, because the state of the system can be symmetric about the barrier, thus satisfying the matching conditions.
Finally, let us consider a periodic orbit at the boundary of the basin of attraction of the marginally stable attractors. In order to obtain this orbit, we increase or decrease all force values z j with j i by an amount such that there is a moment in time where site i (or i + 1) is lifted by site i + 1 (or i) exactly to the threshold. The metastable orbit has now become degenerate with an orbit where site i (or i + 1) is lifted by site i + 1 (or i) infinitesimally above the threshold. This orbit has no barrier, and the matrix S corresponding to this periodic orbit is different from the one corresponding to the metastable orbit. Its largest eigenvalue will therefore be different from 1. We have seen realizations of the interesting case that the largest eigenvalue becomes smaller than 1. This means that the periodic orbit that is at the boundary of the basin of attraction of the marginally stable attractors can itself be an attractor that is reached from a nonvanishing set of initial conditions.

Strongly stable attractors
There is the possibility that all the 1s are flushed out of the system by avalanches that extend from inside the system to the boundary. In this case, the largest eigenvalue of S is of the order α, and the attractor is quickly approached. We have seen many examples of such strongly stable attractors.

Weakly stable attractors
If not all 1s are flushed out of the system, and if there is no barrier in the system, the largest eigenvalue of S belonging to an attractor is 1 − O(α n ) with some power n of α. This corresponds to the situation, where the row containing the 1s remains in the system and is moved around by avalanches coming from both directions.
If α is small or n is large, attractors are approached very slowly. For α = 0.2, we have seen attractors with n = 2 for L = 4 and an estimated n = 11 for L = 20. However, we did not attempt a systematic survey of relaxation times towards the attractors as a function of α and L.

Complex attractors
There exist attractors with strikingly long periods of the order of many thousands for L = 4 or much longer periods for larger L (if period is measured in total force increment per site). Typically, these attractors contain long quasi-periodic sections where the sequence of avalanches remains the same but the force values change slowly, just as one would expect close to a weakly stable or weakly unstable periodic orbit. This quasi-periodic sequence is eventually interrupted by an intermittent phase containing other avalanches, until the quasi-periodic phase is entered again. One can understand the origin of such complex attractors in the following way. Imagine a weakly stable attractor for a certain value of α. Now change α slowly and let the system follow adiabatically. The largest eigenvalue of S on the resulting attractor will have the same coefficients if expanded in powers of α, as long as the avalanches remain the same. Eventually, a value of α will be reached where two avalanches merge or an avalanche splits, changing the product of M matrices, which now can have an eigenvalue larger than 1. Nevertheless, there may still be a region nearby in state space where the old sequence of avalanches can be maintained for a long time if the largest eigenvalue of S of the old attractor was only slightly smaller than 1 (i.e., if the old attractor was weakly stable).

Transient stage
We now present simulation results for larger systems. We first study the transition from a random initial state to a stationary state. Figure 3 shows the force values throughout a system of size L = 1000 for different times. One can see that starting from the boundary more and more sites become synchronized, until in the stationary state all sites apart from a few sites at the boundary have almost the same force z i . We have evaluated the transient time using two different measures of the degree of synchronization: (i) the standard deviation as a function of α and L, where the bar denotes the average taken over all sites i. (ii) Nearestneighbour deviation where again the bar denotes the average over all i and z NN = 1 Time was measured as the number of topplings per site. We decided to take the minimum value of both measures over the time rather than the mean value, since they are very fast oscillating functions, and the minimum value gives smoother data. Since we always begin with random initial conditions, the minima decrease with time as long as the system is not yet in the stationary state. For other, more correlated initial conditions, one might have to measure the closeness to the stationary state in a different way. Figure 4 shows our results for these two synchronization measures, averaged over 1000 random initial configurations. The two figures at the top show the behaviour of σ NN (left) and σ (right) as a function of time and of the system size L. The lower plateau indicates clearly the stationary state. One finds that the transient time is proportional to the system size. This is due to the inward proceeding synchronization, which takes place at a constant rate, if L is sufficiently large. For small L, the boundary layer takes a large part of the system, and there is therefore little synchronization. Apart from the transition to the stationary state with a small value of σ and σ NN , one can also distinguish an earlier transition, where the two measures leave the value 1/12 corresponding to a random initial configuration. We interpret this transition as the onset of the formation of synchronized blocks, after the boundary layer has been set up. The characteristic shape of the curves between these two transitions is strikingly different. While the nearestneighbour deviation decreases rather fast once the synchronization starts, the standard deviation remains on a second plateau until its final decrease. This comes from the two synchronized blocks, which usually have a different force value, as one can see in the second picture of figure 3.
A similar behaviour is found for the dependence on α as shown in the lower two figures for a fixed system size L = 200. For values of α below the critical value α c , the nearest-neighbour deviation depends only weakly on α. The sharp transition at the end of the high plateau of σ NN is linear in α. The bulky part above α c must be due to the fact that a site can now topple twice during the same avalanche. While the onset of synchronization is more clearly visible in the data for σ NN , which decay very rapidly, the transition to the stationary state is much more clearly visible in the data for σ, which remain close to the initial value for a longer time. A more detailed investigation of the transition time to the stationary state reveals the following: (i) over a wide range of α values this transition time depends exponentially on α, as shown in figure 5, where the time is plotted that is needed to reach σ = 0.01 and (ii) for small values of α, the data show a power   law with an exponent around −2.84 ( figure 6). Since the synchronization proceeds very slowly, we did not measure the time to the stationary state, but the time to show for the first time an avalanche larger than 21.
The following analytical arguments suggest that the transient time should indeed diverge at least as fast as α −2 with α. Let us define a time unit as the time during which a force 1 − 2α is added to the system. A site at the centre receives two packages of size α from its neighbours per unit time and topples on average once per unit time. A boundary site receives only one package and has on an average approximately 1 − α topplings per unit time. A site in the synchronized block topples on average y times per unit time, with y being intermediate between these two limit cases, 1 − α < y < 1. Initially, the force difference between the synchronized block and the site that will be synchronized next is of the order of 1. In order to decrease this force difference to a value of the order of α, the difference in the total number of topplings between the block and its neighbour must be of the order of 1/α, which is achieved after a time of the order of 1/α(1 − y). This increases with decreasing α at least as fast as 1/α 2 .

Stationary state
We now turn to systems already in the stationary state and give an overview of their statistical behaviour.
The most striking feature of the stationary state are the large synchronized blocks. Sites within the blocks topple the same number of times, while sites closer to the boundaries topple less often. From time to time, a large avalanche that begins outside the block runs through the entire block. Between the large avalanches, the sites within a block topple mostly one by one, lifting each other almost exactly to the threshold. Let us consider a small region within such a synchronized block and let us show that the sites in this region must have approximately the same height, given the dynamics just described. When the sites topple one by one, their height differences are exactly the same as before, after each site has toppled once. When an avalanche enters the region from outside and extends several sites beyond it, the change in height differences is calculated by multiplying the state vector x with the appropriate product S av of M i ν matrices. and · · · · · ·, analytical result forᾱ = 0.2.ḡ was choosen to be 0.6168576.
If the avalanche passes our region from the right to the left, the elements of S av in a row i belonging to our region are with j ini + 1 being the site that triggered the avalanche. If we denote with x i and x i the values of the force differences within our region before and after the avalanche, we have for an avalanche passing through the region from the right to the left. The asymptotic values of the x i after many avalanches satisfy x i = α(x i+1 + x i−1 ) within the synchronized region. This condition can only be satisfied with all x i being zero or with x i decreasing by a factor of the order α from one value of i to the next. Deep inside the synchronized region, the x i become therefore very small. Directly related to these blocks is the behaviour of t i , the mean number of topplings at a site i (normalized by dividing by the total number of topplings in the system), which we observed as a function of i, averaged over a long time and over many systems, as shown in figure 7.
The fact that sites in the synchronized regions topple the same number of times, while those at the boundaries topple less often (due to the missing neighbours, or due to neighbours toppling less often), can also be explained analytically in the following way: at every site, a local balance equation has to be fulfilled. Letḡ be the mean force increment per unit time, andᾱ the mean package size, which we assume to be the same at each site. The balance equation then is for each site i, with the boundary condition t 0 = t L+1 = 0. Equation (5) can be written in matrix form as where Γ is tridiagonal and given by and d L is a vector with constant entries 1. The numerical solution of equation (5) is also shown in figure 7. An analytical solution can be found by making a continuum approximation to equation (5), The solution that satisfies the boundary conditions is This mean-field result predicts that the thickness of the boundary layer is proportional to α, which agrees with our previous numerical finding that the time until the onset of synchronization is proportional to α. We also considered s i , the mean size of all the avalanches triggered at site i. The result for α = 0.2 and L = 100, averaged over 10 000 synchronized systems, is shown in figure 8. Almost all of the large avalanches are triggered near the boundaries. Also shown in figure 8 is the relative number of avalanches triggered at site i, which also shows narrow peaks at the boundaries of the system, but also a broad peak in the centre.
Combining the two data sets, we arrive at the following scenario: in the stationary state, most of the avalanches are single topplings. All large avalanches are triggered near the boundaries and extended far into the synchronized block. If they do not reach the end of the synchronized block, the rest of the block topples in a series of smaller avalanches, mostly of size 1. These small avalanches cause the broad peak at the centre of figure 8. The structure of the peaks at the boundary of the curves in figure 8 depends on α, and results from averaging over many different stationary states.
We now turn to the size distribution of the avalanches. For two-dimensional systems this is believed to obey a power law (see for example [4,5]). However, the one-dimensional systems show this feature only for short times or for small system sizes. In the stationary state, the avalanche-size distribution looks like the upper left curve of figure 9. It was obtained by averaging over 2140 different systems and 10 9 topplings, where we neglected the first 10 9 transient topplings of each system. For small avalanche sizes, a power law is visible, but only for a single decade and up to the sharp cutoff. We see a large gap, followed by peaks centred at system size and half the system size. The shape of n(s) for small s depends only on the value of α and on the precision used in the simulations, but not on L. We checked this by inserting more sites into the synchronized blocks and comparing the resulting avalanche-size distributions, which only differed for the large avalanches. (See also the bottom right curve of figure 9.) For larger precision, more avalanches are found. (See the upper right curves of figure 9.) The reason for this is that the period of the stationary state is longer, as stated already before. Smaller values of α result also in smaller avalanches, because a higher precision is needed in order to resolve force differences (which scale with powers of α). Note also the nonvanishing weight for avalanches of size 2L for α > α c in the lower left figure.
The results for n(s) confirm the picture that the system is composed of a boundary layer that controls the dynamics and determines the stationary state, and a synchronized block of sites that topple the same number of times and that can be made larger without modifying the boundaries.

Discussion
The investigation of the one-dimensional version of the self-organized critical earthquake models has revealed many intriguing features. When viewed as a dynamical system, the system shows four different types of attractors, all of them being periodic. In contrast to a two-site version of the model, where the variables can only change continuously with time [11], and to a many-site version, where the reset rule is z i → z i − z c [12], we do not find chaotic attractors. In contrast to the one-dimensional Zhang-model, which is conservative and has a stochastic force input [13], the phase space volume does not necessarily shrink for systems that have the same sequence of topplings and avalanches. and · · · · · ·, precision 2 −28 . Bottom left: L = 1000; --, α = 0.1; ---, α = 0.2; · · · · · ·, α = 0.3; and -·-·-·-·, α = 0.4. Bottom right: α = 0.15; --, L = 100; ---, L = 500; and · · · · · ·, L = 1000.
In the stationary state, the model consists of two boundary layers, the thickness of which is larger for longer attractors, and an inner part consisting of one or two synchronized blocks, where all sites have approximately the same force value. The synchronized blocks can be made larger without changing the dynamics of the boundary layer or the period of the attractor. Large avalanches are always triggered near the boundary. These features are clearly reflected in the avalanche-size distribution, where the small avalanches are independent of the system size for sufficiently large systems, while the large avalanches are proportional to it.
Several features of the one-dimensional model are very similar to the two-dimensional model, while others differ. In both versions, large avalanches are only triggered at the boundaries [5], and synchronization proceeds inwards according to a power law in time [2,14]. The inner part is dominated by avalanches of size 1 [3,8,15] even in the stationary state. The computing precision affects the avalanche-size distribution [8]. However, while there are fewer large avalanches for smaller computing precision in one dimension, there are more large avalanches in two dimensions. That the inner part can be made larger without changing the dynamics of the boundary region, must be a special property of the one-dimensional system due to the fact that the boundary of a synchronized block is merely a point and that avalanches can propagate only along lines. Nevertheless, the inner part is to some extent slaved to the boundary region also in two dimensions. The exact interplay between the two still needs to be clarified. However, we want to conjecture that in two dimensions the avalanche-size distribution also separates into two parts, when the system size is only made large enough. We expect the first part to become essentially independent of the system size, while the second part becomes proportional to it. However, in order to verify (or falsify) this conjecture, larger and faster computer simulations are needed than those that have been performed up to now.