The geometry of percolation fronts in two-dimensional lattices with spatially varying densities

Percolation theory is usually applied to lattices with a uniform probability p that a site is occupied or that a bond is closed. The more general case, where p is a function of the position x, has received less attention. Previous studies with long-range spatial variations in p(x) have only investigated cases where p has a finite, non-zero gradient at the critical point p_c. Here we extend the theory to two-dimensional cases in which the gradient can change from zero to infinity. We present scaling laws for the width and length of the hull (i.e. the boundary of the spanning cluster). We show that the scaling exponents for the width and the length depend on the shape of p(x), but they always have a constant ratio 4/3 so that the hull's fractal dimension D=7/4 is invariant. On this basis, we derive and verify numerically an asymptotic expression for the probability h(x) that a site at a given distance x from p_c is on the hull.

Site percolation on a triangular lattice where occupied sites (grey hexagons) are more abundant towards increasing x-coordinates. The hull (black curve) marks the boundary between the vacant spanning cluster on the left and the occupied spanning cluster on the right. The hull width w, defined in (1), is the standard deviation in the hull's x-coordinates.

Introduction
Percolation theory is a powerful tool in describing systems of randomly located objects that are placed on the sites or bonds of a network [1]. The sites or bonds are independently occupied with a probability p or remain vacant with probability 1 − p. The task of percolation theory is to identify the properties of the clusters formed by the connected components. In most theoretical studies, p is assumed to be equal everywhere (i.e. the system is assumed to be uniform). However, in many real-world applications, it is more natural to allow some variation in p among the sites or bonds. In spatial systems, p often depends on the position. For example, if occupied sites represent particles that diffuse from the right-hand edge of a two-dimensional lattice and are absorbed on the left-hand edge, then the concentration of particles exhibits a gradient [2] and the system looks qualitatively like the lattice in figure 1. Gradients also frequently arise in the spatial distribution of biological populations because the environmental conditions may change gradually in space (e.g. from lower to higher elevation along a hillside). If occupied sites represent the presence of a species, the population can exhibit a spatial transition from isolated individuals to a large connected component [3,4].
Given such a transition in connectivity, the focus of gradient percolation is on the hull [5], defined as the interface between adjacent occupied and vacant clusters that both fully span the lattice's y-direction (black curve in figure 1). ‡ The hull has been used in practice, for example, in studies of the surface structure of polymers [7], of the treeline dividing woodland from grassland [3,4], of the boundary of urban settlements [8] and as an example of a non-trapping self-avoiding random walk whose properties have been investigated theoretically [9,10] and with computer simulations [11,12]. Apart from its many applications, gradient percolation is also an active field of theoretical investigation. It allows highly accurate measurements of the smallest occupancy probability p c for which an infinite spanning cluster of occupied sites exists. The value of p c is equal for gradient and non-gradient (i.e. uniform) lattices. Introducing a gradient has the advantage that p does not need to be carefully fine-tuned to achieve excellent numerical results [13,14]. Furthermore, the geometry of the hull is particularly intriguing. Suppose the hull connects the sequence of coordinates (x 1 , y 1 = 0), . . . , (x l , y l = L y ) delineating the boundary of the spanning cluster, where L y is the system size in the y-direction. If p can be approximated near p c by a linear function p(x) = ax + p c , then the hull width w, defined as the standard deviation of the hull's x-coordinates, is known to scale as a −4/7 in the limit a → 0. The hull length l is proportional to L y a −3/7 and, on length scales below w, the hull's fractal dimension is D = 7/4 [2]. These non-trivial scaling laws -recently proved rigorously by Nolin [15,16] -have sparked a lot of interest (see [17] for a survey of the literature) and some researchers have generalized the results to gradients in more complicated lattice models [18][19][20][21][22][23]. So far, all studies have assumed that p(x) can be treated, at least locally, as a linear function (figure 2(a)). However, this is frequently not the case in reality. As a → 0, higher order terms in the Taylor expansion of p may become dominant so that w and l do not diverge, although this is predicted when interpreting w ∝ a −4/7 and l ∝ L y a −3/7 literally. Local linearization is not possible either if the gradient does not exist at all, for instance if p jumps discontinuously from one value to another.
Here we examine representative cases in which no linear approximation for p(x) exists at p c . We investigate the class of functions for b ≥ 0 and a > 0. If b = 0, p(x) makes a discontinuous jump over the percolation threshold p c ( figure 2(b)). § For 0 < b < 1, p(x) is continuous, but the gradient at p c this probability goes rapidly to zero as L y → ∞ [6], see [15] for a rigorous proof. We never encountered this situation in our numerical simulations with L y ≥ 4096. § For functions defined only on the discrete coordinates of a lattice, we cannot, strictly speaking, distinguish continuous from discontinuous behaviour. However, for p(x) defined in (2), we can apply the following criterion. Let x − and x + be the neighbouring left and right coordinates of x in the lattice. If and only if p(x) → p(x − ) and p(x) → p(x + ) in the limit a → 0 + , then the underlying function p of real-valued coordinates must be continuous at x. In this case, we also call the discretized function continuous.  (2)). (c) and (d): Continuous function with infinite or zero slope at x = 0 (0 < b < 1 or b > 1, respectively).
is infinite (figure 2(c)). If b > 1, the gradient is zero (figure 2(d)), so that |p(x) − p c | approaches zero faster than any linear function at x = 0. We have chosen the coefficient a in (2) to be equal on the positive and negative x-axis. This is a convenient choice because, as we show below, the probability to find the hull at a coordinate x is in this case symmetric about x = 0. The results described below can be easily generalized to cases where the coefficients differ in the positive and negative direction. We present numerical simulations for site percolation on triangular lattices, where p c = 1/2, and apply periodic boundary conditions in the y-direction (L y ≥ 4096). For a comparison, we also performed simulations for site and bond percolation on square lattices. We found the same scaling exponents, indicating that they are independent of the lattice type and presumably also valid in continuum models [24].

Hull width, length and fractal dimension
Let us first consider the hull width w in the discontinuous case (b = 0 in (2)). Because a > 0, the hull must navigate around those occupied clusters to the left that are not part of the occupied spanning cluster (figure 3). Similarly, on the right-hand side, excursions of the hull must steer clear of smaller vacant islands that have not merged with the vacant spanning cluster. The typical linear size of these smaller clusters is given by the correlation length ξ which scales in uniform systems ∝ |p − p c | −ν if p ≈ p c . Because the critical exponent ν is equal to 4/3 in two dimensions [1,26] and because the hull width scales in proportion to the cluster sizes, we expect w ∝ a −4/3 . Our numerical simulations confirm this intuition (black filled circles in figure 4a).
If b > 0, the situation is more complicated because p(x), and hence also ξ(p(x)), change gradually as a function of x. For b = 1, Sapoval et al. [2] have pointed out that we can still derive the correct scaling behaviour by assuming that w is proportional to the correlation length at w. The self-consistent relation w ∝ ξ(p(w)), where the local correlation length satisfies ξ(p(x)) ∝ |p(x) − p c | −ν and p(x) = p c + ax, yields the aforementioned w ∝ a −ν/(ν+1) = a −4/7 (black open squares in figure 4a). Carrying out the calculation more generally for p(x) given by (2) with arbitrary b and assuming a to be small, we obtain The same self-consistency requirement has been noticed and tested for the Ising quantum chain [28,29], where the exponent ν however has a value different from the exponent in percolation. In figure 4(a), we compare data from numerical simulations with the prediction of equation (3) based on the two-dimensional percolation exponent ν = 4/3. The data show excellent agreement. There is a similar scaling relation between the normalized hull length l/L y and a, where l is the number of steps in the hull in a lattice of vertical dimension L y . The scaling of l/L y cannot be derived as intuitively as the scaling of the hull width. However, previous numerical evidence in the case b = 1 has supported the hypothesis l/L y ∝ a −3/7 [2], now known to be correct [15]. For arbitrary b, our data (figure 4b) indicate the more general relation which includes the scaling law for b = 1 as a special case. The scaling of l and w is related to the fractal dimension D of the hull as follows. For a given value a = a 1 and system size L y,1 , let us denote the hull length by l 1 and the width by w 1 . If we coarse-grain the hull by measuring it with rulers whose length is A  times longer than the lattice constant, we need l 2 ∝ A −D l 1 rulers to cover the full length of the hull [30]. The width and the system size, however, are linear objects: it requires w 2 ∝ A −1 w 1 rulers to span across the hull's x-and L y,2 ∝ A −1 L y,1 rulers to cover its y-extent. Thus, measured in units that are increased by a factor A, the hull length, width and system size all appear smaller by either a factor A −D or A −1 . Now suppose that the hull geometry scales as w ∝ a −B and l ∝ L y a −C . We obtain the same w 2 and L y,2 as before if the lattice constant remains our fundamental length scale, but we first replace a 1 by a 2 ∝ A 1/B a 1 and then, from the original system of length L y,1 , we cut out a strip of size A −1 L y,1 . In this geometry Comparing this with our previous expression for l 2 , we conclude that Inserting B = ν/(νb + 1) and C = 1/(νb + 1) from equations (3) and (4) we find D = 1 + 1/ν = 7/4, independent of b. We also measured the fractal dimension directly using an independent numerical method (following [12,25] and outlined in the appendix) and again find good agreement with D = 7/4 (figure 4c). The fractal dimension of cluster hulls has recently received a lot of interest. In uniform percolation, where p = p c at all sites (i.e. a = 0 in (2)), D = 7/4 is now known to be exact, thanks to mathematical advances in stochastic Loewner evolution [26]. Recent work by Nolin and Werner [27] proves that the Hausdorff dimension remains 7/4 if p is near, but not exactly equal to p c . In constant gradients (b = 1 in (2)), D = 7/4 was already conjectured in early studies [2] and related to the scaling of w and l. Our results demonstrate that D = 7/4 holds much more generally. Although B = ν/(νb + 1) and C = 1/(νb + 1) depend on b and thus differ in the different scenarios depicted in figure 2, the ratio B/C = ν and hence the fractal dimension D = 1 + 1/ν are invariant properties of the hull geometry.

Hull density profile
The hull width and length are simple measures of the hull geometry, but the distribution of the hull's x-coordinates contains yet more information. For fixed x, we define the hull density h(x) as the fraction of coordinates (x, y) in the lattice that are on the hull. We focus first on the case b = 0 where p(x) = p c ∓ a provided that a is sufficiently small. The ∓ sign is to be interpreted as a minus (plus) sign for negative (positive) x. Our simulations for b = 0 show that h(x) has a maximum at x = 0 and falls off exponentially and symmetrically to both sides (figure 5a), where E 0 and F are positive constants. The exponential decay of h(x) is reminiscent of the off-critical correlation function g c (r) ∝ exp(−|r|/ξ) in uniform percolation, where g c (r) denotes the probability that r and the origin are part of the same finite cluster. The decay rate (F a ν ) −1 has the same scaling as the correlation length ξ, namely ∝ |p − p c | −ν for both p < p c and p > p c . However, there is one important difference between (F a ν ) −1 and the correlation length. In uniform percolation, ξ has different amplitudes on the sub-and supercritical branch [31]. By contrast, we observe that F is equal for x < 0 and x ≥ 0.
Equation (6) has several interesting consequences. First, h(x) is a symmetric function (i.e. h(x) = h(−x)). We observed this symmetry in all investigated lattice types (triangular site percolation, square site and bond percolation). Second, (6) is consistent with the scaling relations (3) and (4) of the width and length for b = 0. This can be proved by inserting (6) after appropriate rescaling of h, the functionh(x) = a 1−ν h(a −νx ) is independent of a (see the data collapse in the inset of figure 5(a)).
The case b > 0, where p(x) changes continuously, requires a little more effort than b = 0, where p(x) is constant except for one discontinuity. Our goal is to approximate h(x) for b > 0 by discretizing p(x). To this end, let us consider the intermediate case where p has more than one discontinuity: p(x) = 0 for x < x 1 , p(x) = π 1 for x ∈ [x 1 , x 2 ), . . ., p(x) = π n−1 for x ∈ [x n−1 , x n ) and p(x) = 1 for x n ≤ x with 0 < π 1 < . . . π n−1 < 1. Let us assume that the length (x i − x i−1 ) of each plateau is larger than the exponential decay length (F |π i − p c | ν ) −1 . In this case, h consists of piecewise exponential functions ∝ exp(±F |π i −p c | ν x) that are continuously joined together at the breakpoints x 1 , . . . , x n (an example is shown in figure 5b). In terms of a differential equation, we can express this as where F is a constant. The plus sign holds for p(x) < p c ; otherwise the minus sign applies. We now make two approximations to obtain an analytic expression for h(x) if p(x) changes continuously. First we assume that the rate of change in p(x) is small compared to the characteristic decay length (F |p(x) − p c | ν ) −1 for all x. This allows us to insert (2) directly into (7). We obtain We know from (3) that the hull width is a vanishing fraction of the width of this interval for a → 0 + . Thus, our second approximation consists of extending (8) over the entire real line. The solution is where we have factored out a (ν−1)/(νb+1) from the integration constant so that, according to (4), the remaining factor E b does not depend on a.

Discussion
Straightforward integration proves that (9) is consistent with the scaling of the hull width in (3). Because (9) also satisfies the scaling of the length (4), it implies the  Figure 6. The transformation from h to f b according to (10) shows that the hull density decays as exp(−F a ν x νb+1 /(νb + 1)) with the same factor F ≈ 4.3 for all b.
correct fractal dimension 7/4 according to (5). In addition to the scaling laws, (9) also provides an analytic expression for the functional form of the hull density profile: h(x) is, in this approximation, proportional to a generalized normal distribution [32]. For b = 1/ν, we predict a Gaussian profile. Another special case is b = 0, where we retrieve the Laplace distribution of (6). For 0 ≤ b < 1/ν the distribution is leptokurtic (i.e. it has a narrower peak and fatter tails than a Gaussian); for b > 1/ν it is platykurtic (broader peak, thinner tails).
For any particular b in (9), the rescaled hull densityh(x) = a (1−ν)/(νb+1) h(a −ν/(νb+1)x ) is independent of a (insets in figures 5(a) and 5(c)). Rescaling can even go further. Above we have assumed that the constant F in (7) is independent of b so that the functions ought to be linear functions with the same slope −F for all b. Figure 6 demonstrates that this is indeed true forx 1, which is valid if x w. From linear regression to f b , we estimate F = 4.3 ± 0.1. For b ≪ 1, a simple line in figure 6 fits the data very well for allx > 0 and this is presumably exact for b = 0. However, as b becomes larger, there are visible departures from the asymptotic line forx ≪ 1. In this regime, p(x) varies significantly over length scales comparable to the correlation length, thus violating the assumption behind (8). If b ≫ 1, we are effectively dealing with a system at the critical point confined to the strip −(p c /a) 1/b < x < ((1 − p c )/a) 1/b . Some properties of critical cluster geometry in a strip are known exactly [33][34][35]. From the known results for the "one-pinch point density" in rectangular domains [36], it is plausible that h(x) is analytic and has a Taylor expansion around x = 0 with leading terms h(x) = h(0) + h ′′ (0)x 2 /2 + O(x 4 ) and h ′′ (0) < 0. As a consequence, as x → 0 we expect log(h) to be approximately quadratic (i.e. log(h(x)) ≈ G 0 − H 0 x 2 ), in contrast to log(h(x)) ≈ G ∞ − H ∞ |x| νb+1 for |x| → ∞ F depends on the lattice type, unlike B, C and D. The stated value is for triangular site percolation. For square site percolation, we find F = 3.9 ± 0.1, for square bond percolation F = 5.4 ± 0.1.  Figure 7. If the hull density is locally approximated as a function proportional to exp(−H|x| J ), the best fitting exponent J increases as νb + 1 in the tail (i.e. as x → ∞), but is close to 2 for b > 3/4 in the core (i.e. as x → 0). as predicted by (9).
We test this hypothesis with nonlinear least-squares regressions of the general form G − H|x| J to fit either the core or the tail of log(h(x)). The best fitting exponents for |x| → ∞ are indeed consistent with J = νb + 1 (figure 7) for all b. Performing the regression around x = 0, on the other hand, we find a crossover from J = νb + 1 to J = 2 at b ≈ 1/ν = 3/4. This suggests that for b < 1/ν, (9) approximates the hull density well, but for b > 1/ν we should replace it with h(x) ≈ E b a (ν−1)/(νb+1) exp − F νb + 1 a ν |x| νb+1 − K b a 2ν/(νb+1) x 2 .
This approximation has the observed asymptotic behaviour for both x = 0 and |x| → ∞: h(x) is Gaussian in the center, but drops more rapidly in the tails. In the inset of figure 5(c) we have used F = 4.3 from above and fitted E b and K b to the hull density for b = 4. The regression curve indeed fits the core of the distribution better than (9) and also provides an excellent approximation in the tails.

Conclusion
In summary, we have investigated a generalization of gradient percolation where the occupancy probability p(x) is nonlinear at the critical point p c . If, to lowest-order, |p(x) − p c | = a|x| b for a > 0 and b ≥ 0, then the hull width and length scale as w ∝ a −ν/(νb+1) and l ∝ a −1/(νb+1) . The hull's fractal dimension D = 7/4 is independent of b. The hull density is symmetric and drops for x far away from p c in proportion to exp(−F a ν |x| νb+1 /(νb + 1)), where F is a lattice-dependent constant independent of a and b. Closer to p c , the hull density becomes approximately Gaussian if b > 1/ν.