# An Exactly Soluble Hierarchical Clustering Model: Inverse Cascades, Self-Similarity, and Scaling

## I Introduction

Clustering and aggregation play an important role in many complex systems. In this paper, we present an inverse cascade model for the self-similar growth of clusters. Elements are introduced at the smallest scale, which then coalesce to form larger and larger clusters. The inverse cascade is terminated by the loss of the largest clusters. The system is thus in a quasi-steady state with the loss of elements in large clusters balanced by the introduction of new elements. The clustering process is recognized to be a branching network similar to a DLA cluster or a river network. Individual clusters are analogous to branches, and coalescence is equivalent to the joining of two branches.

There is a wide range of applications for this analysis. As a specific example, we consider the forest-fire model which has been said to exhibit self-organized criticality. In one version of the forest-fire model, a square grid of sites is considered. At each time step, a model tree or a model spark is dropped on a randomly chosen site. If the site is unoccupied when a tree is dropped, it is “planted.” The sparking frequency $`f`$ is the inverse number of attempted tree drops before a spark is dropped. If the spark is dropped on an empty site, nothing happens; if it is dropped on a tree, it ignites and “burns” all adjacent trees in a model forest fire. In this model, individual trees are introduced at the smallest scale, clusters of trees coalesce to form larger and larger clusters. Significant numbers of trees are lost only in the largest fires that terminate the inverse cascade. The noncumulative frequency-area distribution for the fires is well approximated by a power-law relation

$$N\frac{1}{A^\alpha }$$
(1)   
with $`\alpha 1`$. If the sparking frequency $`f`$ is relatively large, the largest fires are relatively small and the self-similar inverse cascade is valid only over a relatively small range of cluster sizes. If the sparking frequency $`f`$ is small, the fires that terminate the cascade are large and if $`f`$ is sufficiently small the fires will span the entire grid. The noncumulative frequency-area distribution of cluster sizes satisfies equation (1.1) with $`\alpha 2`$ and the cumulative distribution of clusters with area larger than $`A`$ satisfies equation (1.1) with $`\alpha 1`$. The behavior of the one-dimensional forest-fire model has been discussed in terms of a cascade by Paczuski and Bak. The inverse cascade analysis is also applicable to the sandpile model and the slider-block model. In the sandpile model the clusters are the metastable regions that participate in avalanches once they are triggered. In the slider-block model, the clusters are the metastable regions that participate in slip events once they are initiated.

One of the most striking patterns in biology is clusters or aggregations of animals. Examples range from bacteria to whales and include insects, fish, and birds. Bonabeau et al. showed that the frequency-number distribution of whales satisfy equation (1.1) with $`\alpha 1`$. The model we present here should also be applicable to these biological problems.

## II Hierarchical Clustering

We consider a system of stationary entities that we shall refer to as elements. In terms of the forest-fire model, the elements are the trees that are planted on a lattice. The system is growing due to the steady injection of new elements that are added to locations that are not already occupied by previously injected elements. We define connected sets of elements, i.e. groups of elements that are in contact, to be clusters. Note, however, that our model does not require that elements be confined to lattice points. Neighbors can be defined with any metric (e.g. distance) condition, or according to a defined graph structure (e.g. lattice). In the forest-fire model, clusters are the groups of adjacent trees that would burn in a fire if a spark dropped on one of the trees in the cluster. We construct rules for assigning rank to clusters in such a system, based in spirit on the Strahler classification that was originally developed for branching in river networks. In this classification system, a stream with no upstream tributaries is defined to be of rank one; when two rank-one streams combine, they form a stream of rank two, and so forth. However, when streams of different rank combine, the rank of the dominant stream prevails. Our model for the growth of clusters is an extension of a scheme developed earlier which only allowed for the coalescence of clusters of the same rank. The new model is much richer in that it accomodates the coalescence of clusters of all ranks and can, therefore, describe a much wider array of phenomena.

The rules for our cluster model are:

1. We define a single element that is added to a system to be a cluster of rank 1.
2. If a new element is added adjacent to an existing cluster, we say that it is added to the cluster without changing that cluster’s rank, unless the cluster is a single element. In that special case, we define the two elements as forming a cluster of rank 2.
3. If a new element connects two existing clusters of ranks $`i`$ and $`j`$, respectively, then the rank of this new cluster is defined as $`i+1`$ when $`i=j`$ and as $`\mathrm{max}\{i,j\}`$ when $`ij`$. In words, this is equivalent to saying that when two clusters of equal rank coalesce, then the rank increases by one; however, if the two clusters are not of equal rank, then the rank of the larger cluster prevails.
4. If a new element connects three or more clusters, then the rank of the new cluster is defined to be

* the maximal rank of these clusters, when one of the clusters has a rank exceeding that of all of the others, or
    * the maximal rank of these clusters plus one, when there are two or more clusters of the same maximal rank. 

\[This is a rare event—akin to a four-body interaction—and it is ignored in the model equations given below.\]
5. We terminate the inverse cascade of elements from small to large clusters by eliminating clusters of a specified high rank.

In Figure 1, we illustrate how this model works.

We now wish to establish the dynamical equations governing the evolution of this system. Let us define $`N_i`$ to be the number of clusters with rank $`i`$, for $`i1`$. Let $`m_i`$ be the average mass—i.e., the number of elements—of a cluster of rank $`i`$. Then, the total mass $`M_i`$ of the clusters of rank $`i`$ is given by

$$M_i=N_im_i.$$
(2)   
For convenience, we will define the mass of a single element to be one, namely $`m_1=1`$. For example, in two dimensions, we can regard $`m_i`$ as the mean area $`A_i`$ of a cluster of rank $`i`$. This would be the case in the forest-fire model.

We now develop a mean-field approximation describing the dynamical evolution prescribed by the mapping rules given above. As indicated, we ignore the simultaneous coalescence of more than two clusters. We denote the instantaneous change in all quantities using the mapping symbol $``$. Accordingly, when two clusters of ranks $`i`$ and $`j`$ coalesce, the values of $`N_i`$ and $`M_i`$ are modified as follows. For $`i=j`$,

$$N_{i+1}N_{i+1}+1,N_iN_i2,$$
(3)   
$$M_{i+1}M_{i+1}+2m_i,M_iM_i2m_i,$$
(4)   
and for $`i<j`$,

$$N_iN_i1,N_jN_j,$$
(5)   
$$M_jM_j+m_i,M_iM_im_i,$$
(6)   
with equivalent expressions for $`j<i`$. In these equations for $`M_j`$, we have ignored the addition of an element that bridges or joins the two clusters. Since $`m_i`$ will be shown to increase in an essentially geometric progression with respect to the rank $`i`$, the omission of that solitary unit mass in the calculation does not influence the asymptotic properties as $`i\mathrm{}`$.

In our model, coalescence occurs when a new element connects two existing clusters. (We have already indicated that 4-body and higher order effects will be neglected.) Accordingly, in the mean field approximation, we assume that the rate $`r_{ij}`$ of coalescence between clusters of ranks $`i`$ and $`j`$ is proportional to the product of their total numbers, $`N_i`$ and $`N_j`$, and to the product of their boundary sizes, $`\mathrm{}_i`$ and $`\mathrm{}_j`$, and is naturally related to the joint probability of the new element connecting two pre-existing clusters. For example in two dimensions, $`\mathrm{}_i`$ refers to the effective length of the cluster boundary. Thus, we assume that

$$r_{ij}N_i\mathrm{}_iN_j\mathrm{}_j.$$
(7)   
This is an Euclidean approximation, and emerges in the spirit of classical kinetic theory, although the mechanics of this problem is entirely different. In sections IV and VII, this model will be modified to accommodate the possible fractal geometry of clusters.

We now define

$$L_i=N_i\mathrm{}_i$$
(8)   
to be the total size of the boundary associated with clusters of rank $`i`$. We select the normalization for our time-scale so that $`r_{ij}=L_iL_j`$. Accordingly, let $`C`$ be the injection rate of single elements, utilizing this time scale. The evolution of the system can be determined by appropriately adapting equations (2.2)–(2.5). From equations (2.2) and (2.4), we write

$$\dot{N}_1=C2L_1^2\underset{j=2}{\overset{\mathrm{}}{}}L_1L_j,$$
(9)   
$$\dot{N}_i=L_{i1}^22L_i^2\underset{j=i+1}{\overset{\mathrm{}}{}}L_iL_j,\mathrm{for}i>1.$$
(10)   
In equation (2.8), we observe that the rate of change in the number of clusters of rank 1 is equal to the injection rate minus the rate of coalescence of rank 1 clusters together with the rate of coalescence of rank 1 clusters with clusters of larger rank. The factor of 2 appears because two rank 1 clusters were lost in coalescing to form a rank 2 cluster. Meanwhile, in equation (2.9), we observe that the rate of change in the number of clusters of rank $`i`$ is equal to the rate of rank $`i`$ cluster formation from the coalescence of pairs of rank $`i1`$ clusters, minus the rate of coalescence of pairs of rank $`i`$ clusters, together with the rate of coalescence of rank $`i`$ clusters with clusters of larger rank $`j>i`$.

In a similar way, taking into account $`m_1=1`$, we can express the mass-balance in the system, derived from equations (2.3) and (2.5), according to

$$\dot{M}_1=C2L_1^2\underset{j=2}{\overset{\mathrm{}}{}}L_1L_j,$$
(11)   
$`\dot{M}_i=2L_{i1}^2m_{i1}`$ $`+`$ $`{\displaystyle \underset{k=1}{\overset{i1}{}}}L_iL_km_k2L_i^2m_i`$ (12)   
$``$ $`{\displaystyle \underset{j=i+1}{\overset{\mathrm{}}{}}}L_iL_jm_i,\mathrm{for}i>1.`$ (13)   
Note that equations (2.8) and (2.10) are identical since $`M_1=N_1`$.

We observe that the equations above have the potential for self-similarity, since most of the sums are infinite in extent, and might be expected to be convergent. Intuitively, we expect that $`L_j`$ will diminish as $`j`$ increases; while the boundary size of individual clusters of rank $`j`$ increase, their absolute numbers will decrease even more rapidly so that the total boundary size in clusters of rank $`j`$ will be monotone decreasing. The finite sum, which appears in equation (2.11), is somewhat more involved. Nevertheless, it is reasonable to expect that the product of $`m_k`$ with $`L_k`$ will steadily diminish as $`k`$ becomes smaller and that negligible contributions emerge from low values of $`k`$. Finally, it is easy to see that all of the governing rate equations will quickly converge, in the sense of an inverse cascade from $`i=1`$ to some finite cut-off, as $`t\mathrm{}`$. As $`N_1`$ begins to grow, it provides a stimulus to the growth of $`N_2`$, and so on. Similarly, as the masses at each rank in the system grow, they will in turn cause the boundary size $`\mathrm{}_i`$ of each cluster of rank $`i`$ to grow, basically in proportion to some power in $`m_i`$. With this intuition in hand, we now obtain the steady-state solution for this system.

## III Steady State Solution: Cluster and Mass Scaling

We derive a steady state solution for an inverse cascade from equations (2.8) through (2.11). In our inverse cascade, single elements are introduced at the lowest level, and they coalesce to form larger and larger clusters. The inverse cascade is terminated by assuming that very large clusters are removed from the system. We assume that our system develops in a sufficiently large region, so that edge effects can be ignored over a long time. Otherwise, we will have a completely space-filling solution and percolation effects will govern. We can regard this (limited) steady-state solution to be an intermediate asymptotics for our system—our solution will describe the similitude that emerges before percolation and space-filling issues become significant. The steady state solution follows when the time derivatives in the left hand sides of equations (2.8)–(2.11) vanish with the result

$$C=2L_1^2+\underset{j=2}{\overset{\mathrm{}}{}}L_1L_j.$$
(14)   
$$L_{i1}^2=2L_i^2+\underset{j=i+1}{\overset{\mathrm{}}{}}L_iL_j,\mathrm{for}i>1.$$
(15)   
$`2L_{i1}^2m_{i1}`$ $`+`$ $`{\displaystyle \underset{k=1}{\overset{i1}{}}}L_iL_km_k=`$ (16)   
$`2L_i^2m_i`$ $`+`$ $`{\displaystyle \underset{j=i+1}{\overset{\mathrm{}}{}}}L_iL_jm_i,\mathrm{for}i>1.`$ (17)   
As noted earlier, equations (2.8) and (2.10) are equivalent.

Equation (3.2) has a self-similar solution, since that equation is invariant under $`ii+1`$, and depends only on $`L_j/L_i`$. Thus, we seek a solution having the form

$$L_i=ax^{i1}$$
(18)   
where $`0<x<1`$. The first of these constraints on $`x`$ corresponds to boundary sizes being positive, while the second is necessary for the summation to exist. We find that $`x`$ satisfies

$$2x^{2i2}+\underset{j=i+1}{\overset{\mathrm{}}{}}x^{i+j2}=x^{2i4}.$$
(19)   
Summing the infinite geometric series explicitly and dividing by $`x^{2i4}`$, we obtain

$$2x^2+\frac{x^3}{1x}=1,\mathrm{or}x^32x^2x+1=0.$$
(20)   
This equation has a single root in the range $`0<x<1`$, namely $`x=0.55495813\mathrm{}`$ . Given equations (3.1) and (3.6), we find that

$$C=a^2\left[2+x/\left(1x\right)\right]=a^2\mathrm{or}a=C^{1/2}.$$
(21)   
Substitution of these results into equation (3.4) gives

$$L_i=C^{1/2}\left(0.55495813\right)^{i1}.$$
(22)   
We now turn our attention to equation (3.3). We substitute equation (3.4) into equation (3.3), dividing by $`a^2x^{i3}`$ and taking into account equation (3.6). We then obtain

$`2x^{i1}m_{i1}+{\displaystyle \underset{k=1}{\overset{i1}{}}}x^{k+1}m_k`$ $`=`$ (23)   
$`2x^{i+1}m_i+{\displaystyle \frac{x^{i+2}}{1x}}m_i`$ $`=`$ $`x^{i1}m_i.`$ (24)   
This equation does not have an exactly self-similar solution, since it is not invariant under $`ii+1`$. Suppose that we make the substitution

$$x^{i1}m_i=y^{i1},$$
(25)   
assuming that $`y>1`$, whereupon we obtain from summing the finite series

$$2xy^{i2}+x^2\frac{y^{i1}1}{y1}=y^{i1}.$$
(26)   
We observe that the solution for $`y`$ in this equation depends upon $`i`$. However, for large $`i`$, equation (3.11) approximately implies, assuming that we can replace $`y^{i1}1`$ by $`y^{i1}`$, that $`2x+x^2y/(y1)=y`$ which we rewrite as

$$y^2(x+1)^2y+2x=0.$$
(27)   
This equation has a unique solution for $`y>1`$, namely $`y=1/x=1.8019377\mathrm{}`$ . Accordingly, for large $`i`$, we have asymptotic self-similarity with

$$m_i\alpha x^{1i}y^{i1}=\alpha c^{i1},$$
(28)   
where $`c=1/x^2=3.24697602\mathrm{}`$ . With $`m_1=1`$, we have

$$m_i\left(3.24697602\right)^{i1}.$$
(29)   
Before moving to issues dealing with fractals and branching, the solutions we have just obtained for $`L_i`$ and for $`m_i`$ can be immediately exploited. Since $`L_ix^i`$ and, approximately, $`m_ix^{2i}`$, we observe that $`L_i\sqrt{m_i}\mathrm{const}.`$ For example in two dimensions, recalling that $`L_iN_i\mathrm{}_i`$ and introducing the Euclidean relation that $`\mathrm{}_i\sqrt{m_i}`$, it follows that $`N_im_i\mathrm{const}.`$ or, equivalently, we find the number-mass or number-area relationships

$$N_i1/m_i1/A_i.$$
(30)   
This is equivalent to equation (1.1) with $`\alpha =1`$. The branch numbers $`N_i`$ are loosely equivalent to a logarithmic binning of cluster sizes. Logarithmic binning is equivalent to a cumulative distribution. Thus, the result given in equation (3.15) is in agreement with the distribution of cluster sizes obtained from the forest-fire model as discussed above. The concept of clusters can also be extended to both sandpile and slider-block models. In these cases, the clusters are the metastable regions that will avalanche or slip when an event is triggered. In both cases, the cumulative distribution of cluster sizes satisfy equation (1.1) with $`\alpha 1`$. These scaling relationships are archetypical of self-organized criticality. Remarkably, this scaling has been deduced using solely analytic means from our inverse-cascade hierarchical cluster model.

## IV Adaptation for Fractal Perimeter: Cluster and Mass Scaling

In the analysis given in the previous sections, we assumed that the rate of cluster coalescence $`r_{ij}`$ was proportional to the linear dimensions of the two clusters as given in Equation (2.6). We now generalize this dependence to account for the possibility of fractal clusters by introducing an “efficiency” factor $`ϵ<1`$, with an appropriate scaling such that

$$r_{ij}ϵ^{\left|ji\right|}N_i\mathrm{}_iN_j\mathrm{}_j=ϵ^{\left|ji\right|}L_iL_j.$$
(31)   
As before, $`r_{ij}`$ is the rate of coalescence between clusters of ranks $`i`$ and $`j`$. This modification can, for example, describe the increased efficiency with which a smaller cluster can coalesce with a larger one, since the smaller cluster can become attached inside one of the nooks and crannies that can characterize a fractal perimeter.

With this modification, we obtain analogs of equations (2.8)–(2.11)

$$\dot{N}_1=\dot{M}_1=C2L_1^2\underset{j=2}{\overset{\mathrm{}}{}}ϵ^{1j}L_1L_j,$$
(32)   
$$\dot{N}_i=L_{i1}^22L_i^2\underset{j=i+1}{\overset{\mathrm{}}{}}ϵ^{ij}L_iL_j,\mathrm{for}i>1,$$
(33)   
$`\dot{M}_i`$ $`=`$ $`2L_{i1}^2m_{i1}+{\displaystyle \underset{k=1}{\overset{i1}{}}}ϵ^{ki}L_iL_km_k`$ (34)   
$``$ $`2L_i^2m_i{\displaystyle \underset{j=i+1}{\overset{\mathrm{}}{}}}ϵ^{ij}L_iL_jm_i,\mathrm{for}i>1.`$ (35)   
In the steady state, we obtain analogs of equations (3.1)–(3.3)

$$C=2L_1^2+\underset{j=2}{\overset{\mathrm{}}{}}ϵ^{1j}L_1L_j.$$
(36)   
$$L_{i1}^2=2L_i^2+\underset{j=i+1}{\overset{\mathrm{}}{}}ϵ^{ij}L_iL_j,\mathrm{for}i>1.$$
(37)   
$`2L_{i1}^2m_{i1}`$ $`+`$ $`{\displaystyle \underset{k=1}{\overset{i1}{}}}ϵ^{ki}L_iL_km_k=`$ (38)   
$`2L_i^2m_i`$ $`+`$ $`{\displaystyle \underset{j=i+1}{\overset{\mathrm{}}{}}}ϵ^{ij}L_iL_jm_i,\mathrm{for}i>1.`$ (39)   
Substituting equation (3.4) into (4.6) we obtain an analog of (3.6)

$`2x^2+{\displaystyle \frac{x^3}{ϵx}}=1,`$ $`\mathrm{or}`$ $`x^32ϵx^2x+ϵ=0,`$ (40)   
$`\mathrm{or}`$ $`ϵ={\displaystyle \frac{xx^3}{12x^2}}.`$ (41)   
As the $`L_i`$ are positive, $`x`$ must be positive from its definition in (3.4). Suppose that $`ϵ=x`$. Then $`x2x^3=xx^3`$, giving $`x=0`$, a contradiction. Accordingly, for positive $`x`$, the sign of $`ϵx=x^3/\left(12x^2\right)`$ changes only at $`x=\sqrt{1/2}`$ where $`ϵ`$ changes sign as it passes through infinity (due to the denominator). It is easy to see that the sign of both $`ϵx`$ and $`ϵ`$ is positive for $`0<x<\sqrt{1/2}`$, and negative for $`x>\sqrt{1/2}`$. As $`x<ϵ`$ is nesessary for the summation of the geometric series to exist, this implies that $`x<\sqrt{1/2}`$. In addition, the condition $`ϵ<1`$ requires that $`x<0.55495813\mathrm{}`$. For example, $`x=0.5`$ corresponds to $`ϵ=3/4`$. From equations (4.5) and (4.8), we obtain that $`a=xC^{1/2}`$, for any $`ϵ`$.

Let us turn now to the mass balance equation (4.7). Substituting (3.4) and assuming $`x^{i1}m_i=y^{i1}`$, we obtain an analog of equation (3.11)

$$2xy^{i2}+x^2ϵ^{1i}\frac{(ϵy)^{i1}1}{ϵy1}=y^{i1}.$$
(42)   
We assume that $`ϵy>1`$. Precisely as in equation (3.11), we observe that the solution for $`y`$ in this equation depends upon $`i`$. However, for large $`i`$, equation (4.9) approximately implies that $`2x+x^2y/(ϵy1)=y`$ which we rewrite as

$$ϵy^2(x^2+2ϵx+1)y+2x=0.$$
(43)   
Due to equation (4.8), equation (4.10) has a solution $`y=1/x`$, for any $`ϵ`$. Note that condition $`ϵy>1`$ is satisfied for $`y=1/x`$. Accordingly, for large $`i`$, we have asymptotic self-similarity with $`m_i\alpha c^{i1}`$, where $`c=1/x^2`$, as in equation (3.13). For example, when $`ϵ=3/4`$, we have $`c=4`$.

It is important to remember that $`ϵ`$ describes the perimetric fractal scaling for the clusters. The relationship between perimetric and areal scaling remains a controversial topic. However, assuming that one can identify an appropriate link between the two, for example in the context of forest fire or other models, then the preceding discussion makes it possible to identify the frequency-area relationship for fractal clusters, in analogy to the $`N1/A`$ relationship we identified previously for Euclidian clusters.

## V Branching Numbers

In the analogy between clustering and river networks that we have discussed above, we can write for our clusters

$$\frac{N_{i+1}}{N_i}=x^2$$
(44)   
which is known as the bifurcation ratio for river networks. Also, we have

$$\frac{\mathrm{}_i}{\mathrm{}_{i+1}}=x$$
(45)   
which is known as the length-order ratio for river networks. For river networks, the fact that these two ratios are almost constant is known as Horton’s laws.

A major step forward in classifying river networks was made by Tokunaga. He extended the Strahler ordering system to include side branching. A first-order branch joining another first-order branch is denoted by the subscript “11” and the number of such branches is $`N_{11}`$, a first-order branch joining a second-order branch is subscripted “12” and the number of such branches is $`N_{12}`$; a second-order branch joining a second-order branch is subscripted “22” and the number of such branches is $`N_{22}`$.

In order to apply the concept of side branching to the coalescence of clusters, let us suppose that we have a coalescence of two clusters, of ranks $`i`$ and $`j`$. In the case $`i<j`$, the cluster of rank $`i`$ becomes a branch of the cluster of rank $`j`$. Note that, if the smaller cluster has its own branches, these branches are not counted as branches of the larger cluster. However, these branches, together with all of their branches, etc. are counted as subclusters of the larger cluster. In analogy to river networks, branches are to tributaries as clusters are to drainage basins. A branch formed by the cluster of rank $`i`$ is considered to be a subcluster too, and is assigned the rank $`i`$. Any other subcluster is assigned the rank of a cluster from which it first formed as a branch. In analogy to river networks, subclusters of a cluster correspond with the streams in a drainage basin. The case $`i>j`$ is treated similarly. In the case $`i=j`$, both clusters of rank $`i`$ become branches of rank $`i`$ of the new cluster of rank $`i+1`$. Subclusters and their ranks are defined the same way as above.

Let $`t_{ij}`$ be the average number of branches of rank $`i`$ in a cluster of rank $`j`$, for $`i<j`$, and let $`n_{ij}`$ be the total number of sub-clusters of rank $`i`$ in a cluster of rank $`j`$. For $`i=j`$, we define $`t_{ii}=n_{ii}=1`$. By definition, for $`i<j`$ we have

$$n_{ij}=\underset{k=i}{\overset{j1}{}}n_{ik}t_{kj}.$$
(46)   
Moreover, let $`N_{ij}=N_jn_{ij}`$ be the total number of sub-clusters of rank $`i`$ for all clusters of rank $`j`$, and let $`T_{ij}=N_jt_{ij}`$ be the total number of branches. This classification scheme is illustrated in Figure 2. In (a), we have a cluster of rank “1” which corresponds to a single tree in the forest-fire model. In (b), two clusters of rank “1” have coalesced to form a cluster of rank “2.” This cluster has been joined by a cluster of rank “1.” In the forest-fire model, two trees on adjacent grid points have been joined by a third tree. In (c) and (d), clusters of rank “3” and “4” are illustrated. For this example, we have $`n_{12}=n_{23}=n_{34}=3`$, $`n_{13}=n_{24}=11`$, $`n_{14}=43`$, $`t_{12}=t_{23}=t_{34}=3`$, $`t_{13}=t_{24}=2`$, and $`t_{14}=4`$.

As before, we regard the coalescence of more than two clusters as being exceedingly rare and neglect them in our treatment. When two clusters, of ranks $`i`$ and $`j`$ coalesce, we prescribe the mappings for $`N_{ki},N_{kj},T_{ki}`$, and $`T_{ij}`$ as described below. When $`i=j`$,

$$N_{k,i+1}N_{k,i+1}+2n_{ki},N_{ki}N_{ki}2n_{ki},\mathrm{for}k<i,$$
(47)   
$$T_{i,i+1}T_{i,i+1}+2,T_{k,i}T_{ki}2t_{ki},\mathrm{for}ki;$$
(48)   
and when $`i<j`$,

$$N_{kj}N_{kj}+n_{ki},N_{ki}N_{ki}n_{ki},\mathrm{for}ki,$$
(49)   
$$T_{ij}T_{ij}+1,T_{ki}T_{ki}t_{ki},\mathrm{for}ki.$$
(50)   
Given the rate of coalescence $`r_{ij}=L_iL_j`$, we describe the time evolution of the branching process by the following equations

$`\dot{N}_{kj}=2L_{j1}^2n_{k,j1}`$ $`+`$ $`{\displaystyle \underset{i=k}{\overset{j1}{}}}L_iL_jn_{ki}2L_j^2n_{kj}`$ (51)   
$``$ $`{\displaystyle \underset{i=j+1}{\overset{\mathrm{}}{}}}L_iL_jn_{kj},\mathrm{for}k<j,`$ (52)   
from equations (5.4) and (5.6), and

$`\dot{T}_{j1,j}`$ $`=`$ $`2L_{j1}^2+L_{j1}L_j2L_j^2t_{j1,j}`$ (53)   
$``$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}L_kL_jt_{j1,j},\mathrm{for}j>1,`$ (54)   
$$\dot{T}_{ij}=L_iL_j2L_j^2t_{ij}\underset{k=j+1}{\overset{\mathrm{}}{}}L_kL_jt_{ij},\mathrm{for}i<j1,$$
(55)   
from equations (5.5) and (5.7). As before, we turn our focus to the steady state solution of equations (5.8) through (5.10).

## VI Steady State: Branching Numbers

We begin for the steady state case by setting the time derivatives in the left hand sides of equations (5.8)–(5.10) to zero. We obtain

$`2L_{j1}^2n_{k,j1}`$ $`+`$ $`{\displaystyle \underset{i=k}{\overset{j1}{}}}L_iL_jn_{ki}=`$ (56)   
$`2L_j^2n_{kj}`$ $`+`$ $`{\displaystyle \underset{i=j+1}{\overset{\mathrm{}}{}}}L_iL_jn_{kj},\mathrm{for}k<j.`$ (57)   
$`2L_{j1}^2+L_{j1}L_j`$ $`=`$ (58)   
$`2L_j^2t_{j1,j}`$ $`+`$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}L_kL_jt_{j1,j},\mathrm{for}j>1.`$ (59)   
$$L_iL_j=2L_j^2t_{ij}+\underset{k=j+1}{\overset{\mathrm{}}{}}L_kL_jt_{ij},\mathrm{for}i<j1.$$
(60)   
We observe that, due to the finite summation present in equation (6.1), it is not invariant under $`jkjk+1`$ and its solution is not exactly self-similar in $`jk`$. However, we now employ the same methodology used in §III and obtain asymptotically valid approximate solution. In particular, we substitute (3.4) into equation (6.1) and divide by $`a^2x^{j+k4}`$, and we obtain

$`2x^{jk}n_{k,j1}`$ $`+`$ $`{\displaystyle \underset{i=k}{\overset{j1}{}}}x^{ik+2}n_{ki}=`$ (61)   
$`2x^{jk+2}n_{kj}`$ $`+`$ $`{\displaystyle \frac{x^{jk+3}}{1x}}n_{kj}=x^{jk}n_{kj}.`$ (62)   
Based on our result obtained using equation (3.10), we introduce

$$x^{jk}n_{kj}=z^{jk}$$
(63)   
assuming $`z>1`$, and we obtain from summing the finite series in equation (6.4)

$$2xz^{jk1}+x^2\frac{z^{jk}1}{z1}=z^{jk}.$$
(64)   
Approximating $`z^{jk}1`$ by $`z^{jk}`$ in the asymptotic limit $`jk`$, equation (6.6) approximately implies that $`2x+x^2z/(z1)=z`$, or

$$z^2(x+1)^2z+2x=0.$$
(65)   
This latter equation is identical to equation (3.12), and has a unique solution $`z>1`$, namely $`z=1/x=1.8019377\mathrm{}`$ and, thereby, demonstrates that the branching network description preserves the same structural character. Accordingly, for $`jk`$, we have

$$n_{kj}\beta x^{kj}z^{jk}=\beta c^{jk},$$
(66)   
where $`c=1/x^2=3.24697602\mathrm{}`$ as before. Thus, we have approximately

$$n_{kj}\left(3.24697602\right)^{jk}$$
(67)   
in the limit $`jk`$. For the deterministic example given in Fig. 1, we have $`n_{kj}4^{j1}`$ for $`jk`$. Substituting (3.4) into equation (6.3) and dividing by $`a^2x^{2j4}`$, we obtain

$$x^{ij+2}=2x^2t_{ij}+\frac{x^3}{1x}t_{ij}=t_{ij},\mathrm{for}i<j1$$
(68)   
which establishes that

$$t_{ij}=x^{ij+2},\mathrm{for}i<j1.$$
(69)   
\[For the special case that $`i=j1`$, we have from equation (6.2) that $`2+x=t_{j1,j}`$.\] This, now, is functionally equivalent to the similitude relationship assumed by Tokunaga, namely

$$t_{ij}=t_{ji}=ax^{ij}.$$
(70)   
Importantly, the behavior that Tokunaga assumed to be valid emerges in a completely natural way from the underlying mathematics of our inverse cascade. Since $`x=0.55495813`$, we have for our inverse cascade

$$t_{ij}=\left(0.55495813\right)^{ij+1}.$$
(71)   
For the deterministic example given in Fig. 1, we have $`t_{ij}=\left(1/2\right)^{ij+1}`$.

Finally, the connection between our treatment of branching and our earlier treatment of clustering needs to be established. In particular, we observe that $`m_j`$ turns out to be equivalent to $`n_{1j}`$ and that both scale as $`c^{j1}`$ where, as we have already seen, $`c=1/x^2`$.

## VII Adaptation for Fractal Perimeter: Branching Numbers

The branching analysis given in the previous section is easily modified to include the fractal perimeter dependence introduced in equation (4.1). Introducing this relation into equations (5.8)–(5.10), we obtain

$`\dot{N}_{kj}`$ $`=`$ $`2L_{j1}^2n_{k,j1}+{\displaystyle \underset{i=k}{\overset{j1}{}}}ϵ^{ij}L_iL_jn_{ki}`$ (72)   
$``$ $`2L_j^2n_{kj}{\displaystyle \underset{i=j+1}{\overset{\mathrm{}}{}}}ϵ^{ji}L_iL_jn_{kj},\mathrm{for}k<j,`$ (73)   
$`\dot{T}_{j1,j}`$ $`=`$ $`2L_{j1}^2+ϵ^1L_{j1}L_j2L_j^2t_{j1,j}`$ (74)   
$``$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}ϵ^{jk}L_kL_jt_{j1,j},\mathrm{for}j>1,`$ (75)   
$`\dot{T}_{ij}`$ $`=`$ $`ϵ^{ij}L_iL_j2L_j^2t_{ij}`$ (76)   
$``$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}ϵ^{jk}L_kL_jt_{ij},\mathrm{for}i<j1.`$ (77)   
In the steady state, we obtain analogs of equations (6.1)–(6.3)

$`2L_{j1}^2n_{k,j1}`$ $`+`$ $`{\displaystyle \underset{i=k}{\overset{j1}{}}}ϵ^{ij}L_iL_jn_{ki}=`$ (78)   
$`2L_j^2n_{kj}`$ $`+`$ $`{\displaystyle \underset{i=j+1}{\overset{\mathrm{}}{}}}ϵ^{ji}L_iL_jn_{kj},\mathrm{for}k<j.`$ (79)   
$`2L_{j1}^2+ϵ^1L_{j1}L_j`$ $`=`$ (80)   
$`2L_j^2t_{j1,j}`$ $`+`$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}ϵ^{jk}L_kL_jt_{j1,j},\mathrm{for}j>1.`$ (81)   
$`ϵ^{ij}L_iL_j`$ $`=`$ (82)   
$`2L_j^2t_{ij}`$ $`+`$ $`{\displaystyle \underset{k=j+1}{\overset{\mathrm{}}{}}}ϵ^{jk}L_kL_jt_{ij},\mathrm{for}i<j1.`$ (83)   
Substituting equation (3.4) into (7.4) and assuming that $`x^{jk}n_{kj}=z^{jk}`$, we obtain, due to (4.8), an analog of equation (6.6), namely

$$2xz^{jk1}+x^2ϵ^{kj}\frac{(ϵz)^{jk}1}{ϵz1}=z^{jk}.$$
(84)   
Assuming $`ϵz>1`$, we approximate $`(ϵz)^{jk}1`$ by $`(ϵz)^{jk}`$ in the asymptotic limit $`jk`$. In this case, equation (7.7) approximately implies that $`2x+x^2z/(ϵz1)=z`$, or

$$ϵz^2(x^2+2ϵx+1)z+2x=0.$$
(85)   
This latter equation is identical to equation (4.10), and has a solution $`z=1/x`$, for any $`ϵ`$. Accordingly, for $`jk`$, we have asymptotic self-similarity with $`n_{kj}\beta x^{kj}z^{jk}=\beta c^{jk}`$, where $`c=1/x^2`$ as before.

Substituting equation (3.4) into equation (7.6) and dividing by $`a^2x^{2j4}`$, we obtain

$$ϵ^{ij}x^{ij+2}=2x^2t_{ij}+\frac{x^3}{ϵx}t_{ij}=t_{ij},\mathrm{for}i<j1$$
(86)   
which establishes that

$$t_{ij}=ϵ^{ij}x^{ij+2}\mathrm{for}i<j1.$$
(87)   
\[For the special case that $`i=j1`$, we have from equation (7.5) that $`2+x/ϵ=t_{j1,j}`$.\] Thus, our modification of the Euclidean model to accommodate fractal perimetric behavior is complete, and the self-similar description of the branching process has been shown to follow in a completely analogous way.

## VIII Conclusions and Discussion

In this paper, we have presented an inverse cascade model for clustering. This model requires:

1. The addition of single elements at a prescribed small scale;
2. The consideration of the clustering process as a hierarchical tree with side branching;
3. The probability that a cluster of one order will coalesce with another cluster of the same or different order is proportional to the product of the number of trees of the two orders and the square root of their masses (or areas); and
4. Clusters are lost (destroyed) at a prescribed large scale.

Our inverse cascade model provides a general explanation for the behavior of several models that have been considered to exhibit behavior which has often been described as “self-organized criticality” and occurs in various settings including the “forest-fire” model. In this model, the planting of individual trees is the introduction of single elements, and coalescence occurs when a planted tree bridges the gap between two existing clusters. The model “fires” burn significant numbers of trees only in the largest clusters and this terminates the inverse cascade. Our model gives the number-mass (or area) distribution to be $`N1/A`$; this is also found to be the case for the forest-fire model. Our model is also applicable for the sandpile and slider-block models. In the sandpile model, the cluster is the region over which an avalanche will spread once it is initiated. In the slider-block model, the cluster is the region over which a slip event will spread once it is initiated. The initiation of an avalanche in the sandpile model and the initiation of a slip event in the slider-block model are equivalent to a spark being dropped on a tree. In both models the clusters grow by coalescence.

We conclude that these models, which are said to exhibit self-organized criticality, are neither critical nor self-organized. Instead, their behavior is associated with an inverse cascade which asymptotically approaches (so long as the largest scales are not involved) power-law (“fractal”) scaling. This behavior is related to the self-similar direct cascade associated with the inertial-range of fully-developed isotropic turbulence. This behavior qualifies as a form of “intermediate asymptotics”. It is interesting to note that earthquakes, landslides, and actual forest fires also have analogous power-law frequency-area distributions.

We have quantified our inverse cascade in terms of a branching tree hierarchy with side branching. We have adapted the taxonomy used for river networks to the growth of our clusters. The order of each cluster is specified and, in our mean-field approximation, the number of clusters of each order is obtained. We find that this distribution is identical to the self-similar side branching distribution introduced empirically by Tokunaga. This distribution has been found to be applicable for river networks, DLA clusters, and vein structures of leaves.

###### Acknowledgements.

We wish to acknowledge the support of NSF Grant EAR 9804859. We are also grateful to Gleb Morein for several useful discussions.