Geometric universality and anomalous diffusion in frictional fingers

Frictional finger trees are patterns emerging from non-equilibrium processes in particle-fluid systems. Their formation share several properties with growth algorithms for minimum spanning trees (MSTs) in random energy landscapes. We propose that the frictional finger trees are indeed in the same geometric universality class as the MSTs, which is checked using updated numerical simulation algorithms for frictional fingers. We also propose a theoretical model for anomalous diffusion in these patterns, and discuss the role of diffusion as a tool to classify geometry.


Introduction
Frictional finger patterns are a result of flow instabilities in quasi-two-dimensional (2D) deformable media due to frictional and capillary forces [1,2]. Although these patterns have been studied for over a decade, the only means of characterizing their complex geometry has been their channel width. The fingers appear when liquid is withdrawn from a two-phase, particle-fluid system [1][2][3]. The particles are initially distributed throughout the system with an approximately uniform packing fraction f, before the moving fluid-air interface compactify the particle packing. This process leaves behind walls of particles while the invading air forms bifurcating fingers in a tree-like pattern (as illustrated in figures 2 and 3 ). The random geometry of the emerging patterns arises due to non-uniformity in the initial packing fraction. Figure 3 shows several patterns, displaying the range of sizes available. This set of figures is generated numerically, following the procedure outlined in appendix A. Figure 1 shows the 1D skeleton of the pattern, where the finger width has been contracted as in figure 2(c). It is the geometry of this skeleton tree we wish to understand.
The process that generates the frictional finger patterns is in many ways an optimal path finding process that happens in small bursts. The bursts take place along the existing interface where the force needed to overcome the compactified particle front is the smallest. This is very reminiscent of the formation of minimum spanning trees (MST). Here one assigns a weight e, often thought of as an energy, to every link in a graph or lattice. The MST is then the tree spanning all the vertices of the lattice (but not all bonds) such that the total energy is minimized [4]. Hence the MST is a geometry constructed on the basis of global optimization. For the frictional finger structure, a very similar thing happens although the process now is off-lattice. Both processes terminate when the structures are space-filling, i.e. they are both examples of random spanning trees.
The MST universality class is a famous one, to which many systems have been argued to belong [5]. Minimal paths on MST, minimal paths on invasion percolation clusters and watershed lines are examples of random planar curves with the same apparent fractal dimension of 1.22 [5]. Although the frictional finger trees and the MSTs follow similar dynamical construction rules it is not obvious that they share universality class. By numerically measuring various geometric exponents we will see that these differences seemingly does not significantly alter the resulting tree geometry and that the two are in the same universality class.
Once the geometric universality class is known we can predict other exponents by using existing scaling relations. Most interesting perhaps is the relation between the exponents that define the geometric universality class and the exponents of dynamical variables of a random walker. Random walks in random geometries, fractals and tree-like structures often display anomalous diffusion [6][7][8][9][10], whereby the mean-squared Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. displacement of the random walk has a nonlinear asymptotic scaling with time In open and uniform space the diffusion exponent α=1 for any dimension. By contrast, in complex geometries or under flow, the diffusion exponent α may depart from unity, being subdiffusive 0<α<1 or superdiffusive 1<α<2. This type of anomalous transport is typical in complex systems [11][12][13][14][15][16][17][18]. The most  famous example of subdiffusion is perhaps de Gennes 'la fourmi dans le labyrinthe' (the ant in the labyrinth), referring to a random walker on a 2D critical percolation cluster [9]. In this case the diffusion exponent α is determined by the critical percolation exponents [10]. Often α can be expressed in terms of a handfull of parameters reflecting the system geometry and boundary conditions. For example, in the case of a random comb model, the relevant geometric parameter is the tail index (scaling exponent) of the branch-length probability density [7]. In this way, α depends on the class of geometries constructed using a certain type of length distribution. Similarly, the universality class of MST specify such a class of random geometries in the case of spanning trees.
The rest of the paper is organized as follows. In section 2, we introduce the different scaling exponents and fractal dimensions in random tree-like geometries. We also discuss an effective model for anomalous diffusion in the finger geometry, relating the diffusion exponent to the systems geometry. Section 3 discusses statistical measures of branch ordering to classify different tree-like structures. This is then used to determine the Hack exponent of the system. Numerical measurements of scaling exponents are presented in section 4. Finally, concluding remarks are offered in section 5. Numerical details on the frictional finger labyrinths and pattern analysis are included in the two appendices.

Theory
For simplicity, we will make the assumption that the width of the fingers can be ignored. We therefore replace the 2D geometry of figure 2 with the one-dimensional (1D) tree shown in part (c) of the figure. This corresponds to studying the pattern on space and time scales that are much larger than the finger width. A much larger version of the 1D tree is shown in figure 1.

Non-Euclidean fractal measures
Let us consider a 1D tree  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 ,  Î there is a unique nonintersecting 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  and d(a, b) the geodesic distance between them. If j is the function that embeds  into the Euclidean plane, we will write x a a j = ( ) and 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: This is nothing but a circle in the geodesic metric. The average Euclidean distance to these points are then where P s ℓ | ( )| 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 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 which of course depends on the embedding j through the local average ... s á ñ .
The standard Euclidean fractal dimension d f  is defined by [10] m r r , 4 where m is the mass within Euclidean radius r. By assumption, our embedding produces d f =2 for the spacefilling frictional finger trees. We also introduce the scaling exponent of the connected mass m r r , 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  (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 The averaged mass B s, ℓ | ( )| measured using Euclidean length is nothing more that the connected mass in equation (5). Hence, we expect that d d d f c c 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 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.

Anomalous diffusion
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 where τ is the time step and r x x r ) is the projection onto the radial direction of the random walkers step. Expanding in small step size we find By performing an ensemble average we can read of the diffusion coefficient where ... ens á ñ is an ensemble average and D is defined so that the current takes the standard Fickian form j D D r r = - ¶ . 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] j x , Here μ is the mobility. The current should vanish in equilibrium, where the distribution takes on a Boltzmann form [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 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 R d m á ñ~D ℓ , where ... á ñ 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 Using this as the relevant conductivity, the Einstein relation gives the power-law scaling D r r d m -( ) . We can now solve the free diffusion problem. Using the continuity equation for our Fickian current we get the diffusion equation If we think of the diffusion problem with spatially depending diffusivity as a Langevin problem x g r t h = ( ) ( ) with δ-correlated white noise η and g 2 /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 r where d s =4/(2+ξ) is the so-called spectral dimension governing the scaling of the return probability . The distribution is normalized as In particular, the second moment scales as r t 2 á ñ~a, with α=2/(2+ξ). We know that ξ=d m , which implies the diffusion exponent 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].

MST universality class
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] 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 mst 2 . It is also known that the minimum path dimension takes the value d mst 1.22 0.01 m =  ( ) [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 mst . Using equation (10), which should hold for any space-filling tree structure, we also get an approximate diffusion exponent mst 0.62 a = ( ) . 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.

Horton-Strahler statistics
The ordering scheme due to Horton [24] and Strahler [25] is a way to classify topologically complex networks. Recall that our working definition of a tree is a connected set where for every two points there is a unique curve connecting them. We will need some more terminology for trees to proceed. A endpoint of the tree  is a point p such that by removing it, p ⧹ , we still have one connected component. A branching point is a point p  Î such that p ⧹ has at least three disconnected components. Similarly, removing point along a line segment will split the tree into two connected parts. The line segment connecting an endpoint to its closest branching point is called a leaf. By a root we will mean a designated endpoint of the tree, and the line segment connecting this point to a branching point is no longer considered a leaf.

Topological branch ordering
The pruning of a rooted tree is a transformation 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    = AE ( ) the tree has Horton-Strahler number  [26]. When we want to make the Horton-Strahler number of a tree explicit we will write   . An example of these ordering rules are shown in figure 4.
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 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 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 Figure 4. Example of a tree structure with highest order 3. This tree has seven order 1 branches, two order 2 branches and one order 3 branch. The root is denoted R.

The topological fractal dimension
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 . The total mass of the tree can be written as a sum over all orders m n L .
The topological fractal dimension d t can now be defined through this mass and the length of the highest order branch m L d t    ( ) [28]. This can be rewritten as 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 n r w .
Using equation (13) with L r L 2   =and assuming that r L /r B <1 we find as in [28] d r r 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  Î we denote by p  the subtree rooted at p containing all points further away from the main root. Let also m p  ℓ ( ) denote the geodesic length of the largest path containing p in such a subtree. The Hack exponent h is then defined through [30] m 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 m p p  ( ). For self-similar trees it is also expected that m p  ℓ ( )scales in the same way as the highest order branch in the subtree. In this case the Hack exponent should be h 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: 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.

Numerical results
In this section, we numerically calculate the various dimensions and exponent that we have discussed in the theory and Horton-Strahler sections. We also measure the diffusion exponent and compare with equation (18).
The frictional finger patterns used are generated numerically using the scheme presented in appendix A, and they are mapped into 1D trees using the pattern analysis method from appendix B. Table 1 summarizes the values for various exponents.

Connectivity and minimum paths
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 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.

Anomalous diffusion in frictional finger labyrinths
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.
The mean-squared displacement x x t t 0 2 á -ñ~a | | 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.
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][32][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 2 1 DFA a a = -( ) [33]. Hence we expect α DFA ≈1.32, signifying a time signal with positive correlations [31,32]. Table 1. Definitions and values for various exponents. The fourth column shows the value of various exponents based on direct measurements in the frictional finger trees. The last two columns shows the values for frictional finger trees and for MST based on expected scaling relations from equations (16)- (18). The values for MST are based on the value for d m given in [4], and the values for the frictional finger trees are based on the direct measurement of the Hack exponent. In both cases we assume d f =2 is known. The algorithms used for direct measurements are described in the text.

Topological scaling and Hack exponent
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 r 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. 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.
Using equations (16) 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. Figure 9. Graph showing the scaling of branch count n w (red squares) and mean branch length l w (blue circles) with branch order w. Black lines show fits of r B =4.1 and r L =2.15, ignoring the points corresponding to branch order 1 in each case. Uncertainties on n w are assumed equal to n w , and uncertainties on l w are estimated from the standard deviation. There is no uncertainty for l 8 or l 9 due to insufficient data to find a standard deviation, and these values are not used in the fit. The error bars on most other points are smaller than the symbols. For reference, lines with slope 6/10 and 7/10 has been included. We expect that the slopes have an uncertainty of the order±0.03.

Conclusion
We have argued that the frictional finger trees belong to the same universality class as the MSTs. Several geometric exponent were measured, both directly and indirectly, and compared with values for MSTs, which confirmed the hypothesis. The values of the geometric exponents also give a value for the diffusion exponent associated with the universality class, which agrees with random walk simulations. r n r r In addition to the front particle positions r i { }, the local front thickness needs to be stored and updated. When x r r nd i i i  + the front thickness must be updated simultaneously. This happens by the combined action of front stretching, which reduces L i , and mass accumulation, which increases L. The mass accumulation happens because a region of packing density f<1 becomes a f=1 region. The mass added to the front gives an addition f /(1−f ) dx to L. The stretching adds no mass to the front, so that this step conserves the area Ls (see figure A1), that is, d(Ls)=0, so that dL=−Lds/s. The two steps combined then gives  If the front folds back to meet itself, it is stopped, as is illustrated in figure A3. Randomness is introduced by adding a 10% white noise on the initial packing fraction f.

Appendix B. Pattern analysis methods
In order to measure branching statistics and fractal dimensions of our patterns, it is necessary to represent them as logical tree structures. Labyrinths are first rendered as binary images, and a binary closing operation is performed to remove small-scale structure on length scales below that of the frictional fingers. A skeletonizing algorithm then reduces all branches to single-pixel width. A custom algorithm (previously used in [34]) uses this skeletonized image to produce a form of the labyrinth expressed as a hierarchical tree of nodes, in which each node holds information on its parent and descendant nodes, on its position, and on the length and shape of its branch. Leaves with length below the length scale of the binary closing operation are pruned from the tree, as they are unlikely to represent real structural features.