Geometric universality and anomalous diffusion in frictional fingers

Let us consider a 1D tree ${ \mathcal T }$ in which distances are measured by the intrinsic distance function d, given by the shortest-path or the geodesic distance between two points. Thus, our trees are metric spaces consisting of 1D curves that are topologically equivalent to line intervals. For any two points $a,b\in { \mathcal T }$ there is a unique non-intersecting curve connecting them, with a geodesic length d(a, b). This formally describes trees such as those in figures 1 and 2(c).

To convert from the geodesic distance to the Euclidean one, we need to embed our tree into the plane. The only requirement we put on our embedding is that the tree becomes space-filling, to mimic the frictional finger trees or the MST. The space-filling property is measured by the fractal dimension. Let us recall some well-known relations between various fractal dimensions. We will use the mass-length definition of fractal dimension, following the conventions of  [10].

Let a, b be points in ${ \mathcal T }$ and d(a, b) the geodesic distance between them. If φ is the function that embeds ${ \mathcal T }$ into the Euclidean plane, we will write ${{\boldsymbol{x}}}_{a}=\varphi (a)$ and ${{\boldsymbol{x}}}_{b}=\varphi (b)$ for the 2D vectors. The intrinsic distance d and the Euclidean distance are related by the minimum-path dimension d_m, typically introduced as [10]

$$\begin{eqnarray}&&d(a,b)\sim | {{\boldsymbol{x}}}_{b}-{{\boldsymbol{x}}}_{a}{| }^{{d}_{m}(\varphi )}.\end{eqnarray} \tag{ 2 }$$

We will use scalar variables r and ℓ for a generic Euclidean and geodesic distance respectively. To make a global estimate for the typical fractal dimension of shortest paths, we propose the following. Pick a point s inside the tree that is not a branching point or end point. Then consider the set of points a geodesic distance ℓ away from s:

$$\begin{eqnarray*}&&{P}_{s}(\ell )=\left\{p\in {T}| d(s,p)=\ell \right\}.\end{eqnarray*}$$

This is nothing but a circle in the geodesic metric. The average Euclidean distance to these points are then

$$\begin{eqnarray*}&&\langle r{\rangle }_{s}\equiv \displaystyle \frac{1}{| {P}_{s}({\ell })| }\displaystyle \sum _{p\in {P}_{s}({\ell })}| {{\boldsymbol{x}}}_{s}-{{\boldsymbol{x}}}_{p}| ,\end{eqnarray*}$$

where $| {P}_{s}({\ell })|$ are the number of points a geodesic distance ℓ  away from s. Based on this we could define a local estimate for the minimum path dimension. To go to a global estimate we pick a discrete set of sample points S along line segments and perform a weighted average

$$\begin{eqnarray*}&&\overline{r}({\ell })=\displaystyle \frac{{\displaystyle \sum }_{s\in S}{W}_{s}\langle r{\rangle }_{s}}{{\displaystyle \sum }_{s}{W}_{s}},\end{eqnarray*}$$

where the weights W_s are taken to be the length of the line segment containing s as measured between the two nearest branching points. In the remainder of this paper, the overline will signify such an average over sample points. We will then use the following definition for the global minimum path dimension

$$\begin{eqnarray}&&\overline{r}({\ell })\sim {{\ell }}^{1/{d}_{m}(\varphi )}\end{eqnarray} \tag{ 3 }$$

which of course depends on the embedding φ through the local average $\langle ...{\rangle }_{s}$.

The standard Euclidean fractal dimension d_f  is defined by [10]

$$\begin{eqnarray}&&m(r)\sim {r}^{{d}_{f}(\varphi )},\end{eqnarray} \tag{ 4 }$$

where m is the mass within Euclidean radius r. By assumption, our embedding produces d_f = 2 for the space-filling frictional finger trees. We also introduce the scaling exponent of the connected mass

$$\begin{eqnarray}&&{m}_{c}(r)\sim {r}^{{d}_{f}^{c}(\varphi )},\end{eqnarray} \tag{ 5 }$$

where m_c(r) is the mass of the connected part of the structure within radius r from a chosen reference point. That is, if a branch exists the disc of radius r and then enters again somewhere else, the disconnected part is ignored. On length scales where m_c(r)/m(r) is a constant, the two Euclidean fractal dimensions coincide. This is the scale where the system has a well-developed fractal behavior, but where still finite-size effects has not significantly entered. Note that for small radii the ratio m_c(r)/m(r) is close to 1 since all the mass is connected. In larger scales when the two masses deviate we must have that m_c(r) < m(r), so the graph of their ratio initially decrease with radius. However, since we are working with a finite system size the ratio will become unity again when the systems radius is hit. To avoid the finite size effects we therefore work in an intermediate range of radii where m_c(r)/m(r) has not yet began to increase towards 1. We will discuss this again for the frictional fingers in the numerical section.

Finally we introduce the connectivity dimension d_c as a fractal dimension that is measured using the metric of ${ \mathcal T }$ (geodesic distance d) and therefore does not depend on the embedding of the tree. For a ball of geodesic radius ℓcentered at a point s on the tree

$$\begin{eqnarray*}&&B(s,{\ell })=\left\{p\in { \mathcal T }| d(s,p)\lt {\ell }\right\},\end{eqnarray*}$$

we can define the connectivity dimension d_c by looking at how the mass within the ball scales with ℓ, i.e. $| B(s,{\ell })| \sim {{\ell }}^{{d}_{c}}$ [10, 19]. To get a global estimate we will again perform an average over a set of sample point, and we define

$$\begin{eqnarray}&&\overline{| B(s,{\ell })| }\sim {{\ell }}^{\overline{{d}_{c}}}.\end{eqnarray} \tag{ 6 }$$

Since we will only refer to the global connectivity dimension, we will write $\overline{{d}_{c}}\equiv {d}_{c}$ for simplicity. The three fractal measures, ${d}_{m},{d}_{f}$ and d_c, are related by the fact that the mass, i.e. total length, of the tree should be conserved under an embedding. Using equations (3) and (6) we see that

$$\begin{eqnarray*}&&\overline{| B(s,{\ell })| }\sim \overline{r}{\left({\ell }\right)}^{{d}_{c}{d}_{m}}.\end{eqnarray*}$$

The averaged mass $\overline{| B(s,{\ell })| }$ measured using Euclidean length is nothing more that the connected mass in equation (5). Hence, we expect that ${d}_{f}^{c}={d}_{c}{d}_{m}$. As mentioned above, there is a range of length scales where the connected mass and the total mass scales with the same dimension, so in this range we expect that ${d}_{f}^{c}={d}_{f}={d}_{c}{d}_{m}$. Hence we only need two of these three exponents. We already know that our trees have d_f = 2 by construction, so d_c and d_m are the interesting quantities.

To model anomalous diffusion in the space-filling frictional fingers we use an effective medium approach, where the tree structure is replaced with a homogeneous medium with a spatially varying diffusion coefficient.

We are interested in the transport in the radial direction. The current can be written

$$\begin{eqnarray*}&&{j}_{D}=[\rho (r)-\rho (r+\delta r)]\displaystyle \frac{\delta r}{\tau },\end{eqnarray*}$$

where τ is the time step and $\delta r=\hat{{\boldsymbol{r}}}\cdot ({{\boldsymbol{x}}}_{t+1}-{{\boldsymbol{x}}}_{t})$ is the projection onto the radial direction of the random walkers step. Expanding in small step size we find

$$\begin{eqnarray*}&&{j}_{D}=-\displaystyle \frac{\delta {r}^{2}}{\tau }{\partial }_{r}\rho .\end{eqnarray*}$$

By performing an ensemble average we can read of the diffusion coefficient

$$\begin{eqnarray*}&&D(r)={\left\langle \displaystyle \frac{\delta {r}^{2}}{\tau }\right\rangle }_{\mathrm{ens}},\end{eqnarray*}$$

where $\langle ...{\rangle }_{\mathrm{ens}}$ is an ensemble average and D is defined so that the current takes the standard Fickian form ${j}_{D}=-D{\partial }_{r}\rho$. Since the diffusion happens on the tree structure, there may be a non-trivial spatial dependence in the diffusion coefficient.

Whenever the particles move in an external potential V, the equations take the altered forms [20]

$$\begin{eqnarray}&&{{\rm{\partial }}}_{t}\rho =-{\rm{\nabla }}\cdot j(x),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&j(x)=-\mu (x)\rho {\rm{\nabla }}V-D{\rm{\nabla }}\rho .\end{eqnarray} \tag{ 8 }$$

Here μ is the mobility. The current should vanish in equilibrium, where the distribution takes on a Boltzmann form $\rho ={Z}^{-1}{{\rm{e}}}^{-V(x)/{k}_{{\rm{B}}}T}$. Using this to calculate the gradient ${\rm{\nabla }}\rho$ we find

$$\begin{eqnarray*}&&j(x)=[-\mu (x)+D(x)/{k}_{{\rm{B}}}T]({\rm{\nabla }}V)\rho .\end{eqnarray*}$$

For this current to vanish, the Einstein relation $\mu (x)=D(x)/{k}_{{\rm{B}}}T$ must hold locally [20].

The mapping between the frictional fingers and the effective medium is made through the Einstein relation for conductivity. We will demand that the effective medium satisfies the same Einstein relation as in the tree structure. In the presence of an electric field the mobility reads $\mu (x)=\sigma (x)/{{nq}}^{2}$, where σ is the conductivity, n the equilibrium particle mass density and q the particle charge [10]. At large times the particle density n is uniform throughout the space-filling tree, so it has no interesting scaling. The Einstein relation then implies the scaling D(r) ∼ σ(r), where we assumed a radial dependence only.

To extract the spatial scaling of the diffusion coefficient consider two large concentric circles in the frictional finger tree, centered on the initial position of the random walker. Let Δ R denote the radial distance separating the two circles. Since the tree has a statistical self-similarity we expect that if we remove all the shortest branches and increase ΔR the statistics of the paths connecting the two circles remains the same. In particular, the number of paths remains the same. A random walker starting at the inner circle will pick one of the paths on its journey to the outer circle, and the typical geodesic distance of the path is $\langle {\ell }\rangle \sim {\rm{\Delta }}{R}^{{d}_{m}}$, where $\langle ...\rangle$ is some average over the fixed number of paths. The conductivity of a single path scales with its inverse length, so the effective conductivity should behave a leading order scaling behavior

$$\begin{eqnarray*}&&\langle \sigma \rangle \sim \left\langle \displaystyle \frac{1}{{\ell }}\right\rangle =\displaystyle \frac{1}{\langle {\ell }\rangle }+...\sim {\rm{\Delta }}{R}^{-{d}_{m}}.\end{eqnarray*}$$

Using this as the relevant conductivity, the Einstein relation gives the power-law scaling $D(r)\sim {r}^{-{d}_{m}}$.

We can now solve the free diffusion problem. Using the continuity equation for our Fickian current ${j}_{D}=-D(r){\rm{\nabla }}\rho$ we get the diffusion equation

$$\begin{eqnarray}&&{{\rm{\partial }}}_{t}\rho -\displaystyle \frac{1}{r}{{\rm{\partial }}}_{r}\left[rD(r){{\rm{\partial }}}_{r}\rho \right]=0.\end{eqnarray} \tag{ 9 }$$

If we think of the diffusion problem with spatially depending diffusivity as a Langevin problem $\dot{x}=g(r)\eta (t)$ with δ-correlated white noise η and g²/2 = D an interpretation of the stochastic integrals is needed to derive the Fokker–Planck/diffusion equation. Equation (9) corresponds to the Hänggi interpretation, in contrast to the Ito and Stratonovich interpretations which have an additional drift term proportional to ∂_r D(r) [20, 21].

The solution of equation (9) for a power-law diffusivity $D(r)={D}_{0}{r}^{-\xi }$ is known to be [22]

$$\begin{eqnarray*}&&\rho (r,t)=\displaystyle \frac{2+\xi }{2\pi {\rm{\Gamma }}({d}_{s}/2)}{\left[\displaystyle \frac{1}{{D}_{0}{\left(2+\xi \right)}^{2}t}\right]}^{\tfrac{{d}_{s}}{2}}\exp \left(-\displaystyle \frac{{r}^{2+\xi }}{{D}_{0}{\left(2+\xi \right)}^{2}t}\right),\end{eqnarray*}$$

where d_s = 4/(2 + ξ) is the so-called spectral dimension governing the scaling of the return probability $\rho (0,t)\sim {t}^{-{d}_{s}/2}$. The distribution is normalized as ${\int }_{0}^{\infty }{\rm{d}}{r}(2\pi r)\rho =1$. The moments of the distribution reads

$$\begin{eqnarray*}&&\langle {r}^{m}\rangle =\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{m+2}{\xi +2}\right){\left[{D}_{0}{\left(\xi +2\right)}^{2}t\right]}^{\tfrac{m}{\xi +2}}}{{\rm{\Gamma }}\left(\tfrac{2}{\xi +2}\right)}.\end{eqnarray*}$$

In particular, the second moment scales as $\langle {r}^{2}\rangle \sim {t}^{\alpha }$, with α = 2/(2 + ξ). We know that ξ = d_m, which implies the diffusion exponent

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{2}{2+{d}_{m}}.\end{eqnarray} \tag{ 10 }$$

Had we not assumed that the particle density was uniform throughout a space-filling structure, the number 2 in the denominator would have been replaced with the fractal dimension, yielding the standard equation for diffusion exponent in trees [10].

An MST is an example of a spanning tree generated by finding minimal energy configurations on a lattice. If we write (i, j) for a link connecting sites i and j on the square lattice we assign some energy _ij to the links using some probability distribution P(). A subset S of the lattice now has the energy [4]

$$\begin{eqnarray*}&&E(S)=\displaystyle \sum _{(i,j)\in S\subset {{\mathbb{Z}}}^{2}}{\epsilon }_{{ij}}.\end{eqnarray*}$$

An MST is the tree on the lattice with lowest energy. This globally optimized structure can be obtained through algorithms based on a local optimization procedure [23]. Here one chooses an initial site in the lattice and grows a spanning tree by first choosing the bond connected to the initial site with the lowest energy, and the proceeds by adding minimal energy links connecting sites in the cluster to ones neighboring the cluster. One can only add bonds that does not form a loop in the cluster. The resulting geometry is independent of the initial site and is the unique MST on the lattice provided that the energies _ij of bonds are unique [23]. This can in practice be assumed to hold, since if it were not the case some infinitesimal perturbation of the energies can be applied to make them unique.

It is well known that there is some universality associated with the MST. In particular, the values of the link energies are irrelevant—only the ordering of energies matter [4]. Clearly, if we were to shift all energies  → f() in such a way that the order remains unchanged, the same links will be invaded at every time step. This set of transformations on the energy landscape leaves invariant the final geometry of the resulting MST. This freedom in choosing the energies by any order-preserving function f  is the a source of universality for the MST [4]. For example, it does not matter if the energies are distributed close to each other or with large fluctuations, as long as the energy hierarchy is the same.

Since it is a spanning tree, the MST has fractal dimension ${d}_{f}({\rm{mst}})=2$. It is also known that the minimum path dimension takes the value ${d}_{m}({\rm{mst}})=1.22\pm 0.01$ [4]. This is also the minimum path dimension of strands in invasion percolation, optimal path cracks and fractal watershed lines [5]. From equation (6) we see that the connectivity dimension for MSTs should be roughly ${d}_{c}({\rm{mst}})\approx 2/1.22\approx 1.64$. Using equation (10), which should hold for any space-filling tree structure, we also get an approximate diffusion exponent $\alpha ({\rm{mst}})=0.62$.

As we noted in the introduction, the formation of frictional finger trees follow similar rules as the MST. The interface in the fluid-particle system evolves by finding the region along the interface where the energy barrier is smallest, and then evolves until friction stops the growth. Hence the frictional finger trees are formed by local optimization rules. Also in this case it is the ordering that matters. If one were to partition the interface into boxes, each with an average energy barrier height, it is clear that the interface will evolve in the box with smallest energy. Once the interface has evolved it has become longer, and more boxes are needed for the partition. New energies corresponding to new boxes are then added to the hierarchy of energies. The evolution of the interface then proceeds again by finding the smallest energy. The distribution of energy barriers for the frictional finger case depends on the random initialization of the packing fraction. Although the frictional finger trees are constructed in the continuum while the MST on a lattice there is a lot of similarity in the two systems. Whether or not this analogy can be made into a statement of equivalence is one of the main questions asked in this paper.

The pruning of a rooted tree is a transformation

$$\begin{eqnarray*}&&{ \mathcal P }:{ \mathcal T }\to { \mathcal P }({ \mathcal T })\end{eqnarray*}$$

that gives a new tree obtained by removing all leaves from the original tree [26]. The order of a line segment in the tree is now defined as the number of pruning transformations needed to remove it. The union of a collection of connected line segments with the same order is called a branch. The Horton–Strahler number of the tree is defined as the number of pruning transformations needed to eliminate the tree in its entirety, i.e. if ${{ \mathcal P }}^{{ \mathcal H }}({ \mathcal T })=\varnothing$ the tree has Horton–Strahler number ${ \mathcal H }$ [26]. When we want to make the Horton–Strahler number of a tree explicit we will write ${{ \mathcal T }}_{{ \mathcal H }}$. An example of these ordering rules are shown in figure 4.

Zoom In Zoom Out Reset image size

Note that the Horton–Strahler number is a topological invariant of the tree—there is no reference to a metric when defining it. It is also a measure of the size and complexity of the tree. Another interesting topological invariant for trees is the bifurcation ratio. Let n_w be the number of branches with order w. Following [27], the bifurcation ratio is defined as

$$\begin{eqnarray}&&{r}_{B}(w)=\displaystyle \frac{{n}_{w}}{{n}_{w+1}}.\end{eqnarray} \tag{ 11 }$$

This quantity contains information regarding the self-similarity of the tree. The type of self-similarity is a topological one because it only relies on the counting of branches- if the bifurcation ratio is independent of branch order r_B(w) = r_B the structure has r_B more branches at order w than at order w + 1, and r_B more branches at order w + 1 than at order w + 2 et cetera. The termination of this process is dictated by the Horton–Strahler number, i.e. when $w={ \mathcal H }-1$.

Analogously to the bifurcation ratio we can define a length scaling ratio. Let L_w be the average internal length of a branch of order w. Then the length scaling ratio is defined as

$$\begin{eqnarray}&&{r}_{L}(w)=\displaystyle \frac{{L}_{w}}{{L}_{w-1}}.\end{eqnarray} \tag{ 12 }$$

With the Horton–Strahler parameters defined we can analyze another generalized dimension of our system. This analysis closely resembles that in [28, 29], but we review it here for the sake of completeness.

Consider a tree with given Horton–Strahler parameters r_B(w), r_L(w) and ${ \mathcal H }$. The total mass of the tree can be written as a sum over all orders

$$\begin{eqnarray*}&&m({{ \mathcal T }}_{{ \mathcal H }})=\displaystyle \sum _{w=1}^{{ \mathcal H }}{n}_{w}{L}_{w}.\end{eqnarray*}$$

The topological fractal dimension d_t can now be defined through this mass and the length of the highest order branch $m({{ \mathcal T }}_{{ \mathcal H }})\sim {L}_{{ \mathcal H }}^{{d}_{t}}$ [28]. This can be rewritten as

$$\begin{eqnarray}&&{d}_{t}({{ \mathcal T }}_{{ \mathcal H }})=\mathop{\mathrm{lim}}\limits_{{ \mathcal H }\to \infty }\displaystyle \frac{\mathrm{log}m({{ \mathcal T }}_{{ \mathcal H }})}{\mathrm{log}{L}_{{ \mathcal H }}}.\end{eqnarray} \tag{ 13 }$$

The number of branches of a given order n_w can be written as a product of bifurcation ratios. Using the definition in equation (11), we have

$$\begin{eqnarray*}&&{n}_{w}=\displaystyle \prod _{w^{\prime} =w}^{{ \mathcal H }}{r}_{B}(w^{\prime} ).\end{eqnarray*}$$

Similarly, using equation (12) the average lengths can be written

$$\begin{eqnarray*}&&{L}_{w}{/L}_{1}=\displaystyle \prod _{{w}^{{\rm{{\prime} }}}=2}^{w}{r}_{L}({w}^{{\rm{{\prime} }}}).\end{eqnarray*}$$

We will set the leaf lengths L₁ = 1, for simplicity. Assuming that we can approximate our tree with a self similar tree we get ${n}_{w}={r}_{B}^{{ \mathcal H }-w}$ and ${L}_{w}={r}_{L}^{w-2}$. The mass is then in the form of a geometric series

$$\begin{eqnarray*}&&m({{ \mathcal T }}_{{ \mathcal H }})=\displaystyle \frac{{r}_{B}^{{ \mathcal H }}}{{r}_{L}^{2}}\displaystyle \sum _{w=1}^{{ \mathcal H }}{\left(\displaystyle \frac{{r}_{L}}{{r}_{B}}\right)}^{w}=\displaystyle \frac{{r}_{B}^{{ \mathcal H }}}{{r}_{L}({r}_{B}-{r}_{L})}\left[1-{\left(\displaystyle \frac{{r}_{L}}{{r}_{B}}\right)}^{{ \mathcal H }}\right].\end{eqnarray*}$$

Using equation (13) with ${L}_{{ \mathcal H }}={r}_{L}^{{ \mathcal H }-2}$ and assuming that r_L/r_B < 1 we find as in [28]

$$\begin{eqnarray}&&{d}_{t}({{ \mathcal T }}_{\infty })=\displaystyle \frac{\mathrm{log}{r}_{B}}{\mathrm{log}{r}_{L}}.\end{eqnarray} \tag{ 14 }$$

This topological fractal dimension is expected to satisfy the same relation as the connectivity dimension, i.e. d_f = d_t d_m [30], and hence we expect that d_c = d_t. This will be checked numerically in the next section.

The topological fractal dimension is closely related to the so-called Hack exponent. For any point $p\in { \mathcal T }$ we denote by ${{ \mathcal T }}_{p}$ the subtree rooted at p containing all points further away from the main root. Let also ${{\ell }}_{m}({{ \mathcal T }}_{p})$ denote the geodesic length of the largest path containing p in such a subtree. The Hack exponent h is then defined through [30]

$$\begin{eqnarray}&&{m({{ \mathcal T }}_{p})}^{h}\sim {{\ell }}_{m}({{ \mathcal T }}_{p}).\end{eqnarray} \tag{ 15 }$$

In the study of river network topology, this exponent is typically defined through the relation between the drainage basin area connected to p, sometimes denoted a_p, and the maximal geodesic length inside it, but for space-filling structures we expect ${a}_{p}\sim m({{ \mathcal Y }}_{p})$. For self-similar trees it is also expected that ${{\ell }}_{m}({{ \mathcal T }}_{p})$ scales in the same way as the highest order branch in the subtree. In this case the Hack exponent should be $h=1/{d}_{t}$. It is known that the Hack exponent indeed satisfies the inverse of equation (14) [30]. We will measure the Hack exponent independently in addition to the Horton–Strahler ratios in the next section.

Before we turn to the numerical section we want to point out that for space-filling systems with d_f = 2 every other geometric exponent discussed in this paper can be expressed in terms of the Hack exponent by using the discussed scaling relations. In summary:

$$\begin{eqnarray}&&{d}_{c}={d}_{t}=1/h,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{d}_{m}=2h,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\alpha ={\left(1+h\right)}^{-1},\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{r}_{L}={r}_{B}^{h}.\end{eqnarray} \tag{ 19 }$$

The last of these equations can be seen as a consistency condition between the Horton–Strahler ratios and the Hack exponent for self-similar systems in the thermodynamic limit.

The connectivity dimension d_c and the minimum-path dimension d_m are defined by equations (6) and (3) respectively. A set of sampling points a is chosen such that it contains one random point along the length of each line segment of the tree. Then for a series of lengths l along the branches, we measure both the total branch length within l from a, $| B(a,l)|$, and the mean Euclidean distance D from a of the set of points exactly l  from a. For each l, we take a weighted average by branch length of both $| B(a,l)|$ and D across all points a to derive mean values for the whole pattern. To avoid influence from the pattern's edge, only points a at least a geodesic distance l from the labyrinth perimeter are considered.

The fractal dimensions d_c and d_m can be found by plotting $| B(a,l)|$ and D respectively against l, as is shown in figure 6, using data from the largest labyrinth with Horton–Strahler number 9. Numerical values are obtained of d_c = 1.67 ± 0.05 and d_m = 1.25 ± 0.03. Note that ${d}_{c}{d}_{m}={d}_{f}=2.09\pm 0.08$, which, in theory, should correspond to the Euclidean fractal dimension d_f.

The range corresponding to the unshaded region in figure 6 is taken from figure 5, which shows the ratio m_c(r)/m(r). We see that as expected the graph interpolated between 1 at small r and 1 at large r, and in between stabilizes in a range of radii. This range depends on the size of the system. The data in figure 6 corresponds to the largest system.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Diffusion in frictional fingers is studied by Brownian random walkers. Figure 7 shows the schematic setup of the simulations. A discrete random walk is released inside the frictional fingers, with a lattice spacing that is smaller than the finger width by one order of magnitude. For the sake of simplicity, we use hard-wall boundary conditions, i.e. when the particle hits the walls that step is discarded and a new step is taken.

Zoom In Zoom Out Reset image size

The mean-squared displacement $\langle | {x}_{t}-{x}_{0}{| }^{2}\rangle \sim {t}^{\alpha }$ is then calculated for labyrinths of different sizes. Figure 8 shows the simulation results in systems where the sizes differ by a factor of 16. The largest system corresponds to figure 1. We see that the diffusion exponent decreases with system size, and for the biggest system we have α ≈ 0.64.

Zoom In Zoom Out Reset image size

In figure 8 the slopes are found by the best fit of the data points. The diffusion exponent could also be found through detrended fluctuation analysis (DFA). In general, DFA is a tool that can be used to study correlations or scaling for long time series [31–33]. By applying the methods of DFA to the random walkers position x_t, one can through the scaling exponent α_DFA of the DFA fluctuation function find the anomalous diffusion exponent through $\alpha =2({\alpha }_{\mathrm{DFA}}-1)$ [33]. Hence we expect α_DFA ≈ 1.32, signifying a time signal with positive correlations [31, 32].

From the simplified 1D tree structure, it is possible to find the Horton–Strahler order of each branch, and to report the number and mean length of branches of each order. These trees are not perfectly self-similar, but to a good approximation we can use an average value of the Horton–Strahler ratios when doing our calculations. Figure 9 shows the number of branches of a given order and their average geodesic length in a logarithmic scale. Note that in a self-similar tree we should have ${n}_{w}={r}_{B}^{{ \mathcal H }-w}$, so that $\mathrm{log}{n}_{w}\sim -w\mathrm{log}({r}_{B})$. Hence the slope in this plot determines the bifurcation ratio for a self-similar structure, giving r_B ≈ 4.1 ± 0.15. Similarly we find r_L = 2.15 ± 0.10 from the slope of the blue (circular) data points in figure 9.

Zoom In Zoom Out Reset image size

Figure 10 shows the mass of subtrees versus maximal geodesic length, which gives us the Hack exponent h = 0.60 ± 0.015. Using equation (19) one can show that the measured value for h and the values for r_B and r_L indeed are consistent, within the uncertainties. Using the numerical values of the ratios and equation (19), solving for h gives the value 0.54 ± 0.06. However, there are larger sources of error in the measurements of the ratios, so the direct measurement of the Hack exponent is more reliable.

Zoom In Zoom Out Reset image size

Using equations (16)–(18) we have from the measured Hack exponent the values

$$\begin{eqnarray*}&&{d}_{t}=1.67\pm 0.04;{d}_{m}=1.20\pm 0.03;\alpha =0.63\pm 0.01.\end{eqnarray*}$$

We see this value of d_t agrees very well with the direct measurement of the connectivity dimension d_c as expected. The diffusion exponent also agrees with simulations on frictional fingers. Within the uncertainties these values all agree with those of the MST universality class.