Introduction

Percolation was first introduced by Flory to describe the gelation of polymers1 and later studied in the context of physics by Broadbent and Hammersley2. This model is considered the paradigm of connectivity and has been extensively applied in several different contexts, such as, conductor-insulator or superconductor-conductor transitions, flow through porous media, sol-gel transitions, random resistor network, epidemic spreading and resilience of network-like structures3,4,5,6,7,8,9,10. In the lattice version, lattice elements (either sites or bonds) are occupied with probability p and a continuous phase transition is observed at a critical probability pc, where for p < pc, as the correlation function decays exponentially, all clusters are of exponentially small size and for p > pc there is a spanning cluster. At pc, the spanning cluster is fractal11. In this article we focus on the shortest path, defined as the minimum number of lattice elements which belong to the spanning cluster and connect two opposite borders of the lattice12,13. The shortest path is related with the geometry of the spanning cluster12,14,15,16,17. Thus, studies of the shortest path resonate in several different fields. For example, the shortest path is used in models of hopping conductivity to compute the decay exponent for superlocalization in fractal objects18,19. It is also considered in the study of flow through porous media to estimate the breakthrough time in oil recovery20 and to compute the hydraulic path of flows through rock fractures21. The shortest path has even been analyzed in cold atoms experiments to study the breakdown of superfluidity22. However, despite its relevance, the fractal dimension of the shortest path is among the few critical exponents in two-dimensional percolation that are not known exactly23,24.

Let us consider critical site percolation on the triangular lattice, in a two-dimensional strip geometry of width Lx and height Ly (Ly > Lx), in units of lattice sites, see Fig. 1. Each site is occupied with probability p = pc. See Methods for details on the algorithm used to generate the curves. The largest cluster spans the lattice with non-zero probability and the average shortest path length 〈l〉, defined as the number of sites in the path, scales as , where dmin is the shortest path fractal dimension and its best estimation is dmin = 1.13077(2)24,25. There have been several attempts to compute exactly this fractal dimension26,27,28,29,30,31,32. Most tentatives were based on scaling relations, conformal invariance and Coulomb gas theory. But the existing conjectures have all been ruled out by precise numerical calculations. For example, Ziff computed the critical exponent g1 of the scaling function of the pair-connectiveness function in percolation using conformal invariance arguments33. g1 has been conjectured to be related to the fractal dimension of the shortest path31. In turn Deng et al. conjectured a relation between dmin and the Coulomb gas coupling for the random-cluster model32. Both conjectures were discarded by the latest numerical estimates of dmin24,34. Thence, as recognized by Schramm in his list of open problems, a solid theory for the shortest path is still considered one of the major unresolved questions in percolation35.

Figure 1
figure 1

A spanning cluster on the triangular lattice in a strip of vertical size Ly = 512.

The shortest path is in red and all the other sites belonging to the spanning cluster are in blue.

Full size image

Impressive progress has recently been made in the field of critical lattice models using the Schramm-Loewner Evolution theory (SLE). In SLE, random critical curves are parametrized by a single parameter κ, related to the diffusivity of Brownian motion. Let us consider the case of a non self-touching curve, like the shortest path, defined in the upper half plane , that starts at the origin and grows towards infinity. Under a proper choice of parameters, it is possible to define a unique conformal map gt from , i.e. the upper half-plane minus the curve γ[0, t], onto such that there exists a continuous real function ξt and gt satisfies the stochastic Loewner differential equation,

with g0(z) = z. The function ξt is called driving function. For details about the conformal map gt see Supplementary Information online. We define chordal SLEκ as the random collection of conformal maps in the upper-half plane that satisfy the Loewner equation with a driving function , where Bt is a one-dimensional Brownian motion.

With the value of κ, one can obtain exactly several probability distributions for the curve, allowing to compute, for example, crossing probabilities and critical exponents36,37,38. SLE has been shown to describe many conformally invariant scaling limits of interfaces of two-dimensional critical models. In particular, SLE6 has first been conjectured39 and later proved on the triangular lattice36 to describe the hull in critical percolation40. SLE has been successfully used to compute rigorously other critical exponents of percolation-related objects38,41 as, for example, the order parameter exponent β, the correlation length exponent ν and the susceptibility exponent γ38. More recently, the probability distributions of the hulls of the Ising model42,43,44,45 and of the Loop Erased Random Walks39,46,47 were computed exactly. Therefore, it is legitimate to ask if the SLE techniques can help solving the long standing problem of the fractal dimension of the shortest path.

Also, a possible description of the physical process through SLE gives interesting insights in new ways of generating the shortest path curves. Once SLEκ is established, the value of κ suffices to generate, from only a Brownian motion, curves having the same statistical properties as the shortest path48,49,50. This can be very useful in the case of problems involving optimization processes like the shortest path, watersheds51, or spin glass problems52,53,54, as traditional algorithms imply the exploration of large areas.

In this article, we will show that the numerical results are consistent with SLE predictions with κ = 1.04 ± 0.02. SLEκ curves have a fractal dimension df related to κ by 55. From the estimate of the fractal dimension of the shortest path, one deduces the value of the diffusion coefficient κ corresponding to an SLE curve of same fractal dimension; κfract = 1.0462 ± 0.0002. In what follows, we compute three different estimates of κ using different analyses and compare them to κfract. In particular we consider the variance of the winding angle39,56,57, the left passage probability58 and the statistics of the driving function53,59. All estimates are in agreement with the one predicted from the fractal dimension and therefore constitute a strong numerical evidence for the possibility of an SLE description of the shortest path.

Results

Winding angle

The first result related to SLE deals with the winding angle. For each shortest path curve we have a discrete set of points zi, called edges, on the lattice. The winding angle θi at each point zi can be computed iteratively as θi+1 = θi + αi, where αi is the turning angle between the two consecutive points zi and zi+1. Duplantier and Saleur computed the probability distribution of the winding angle for random curves using conformal invariance and Coulomb gas techniques56. According to their result39, for SLEκ, the winding angle along all the edges of the curve exhibits a Gaussian distribution of variance

where b is a constant and Ly is the vertical lattice size57. Therefore, κ/4 corresponds to the slope of 〈θ2〉 against ln(Ly). Figure 2 shows the results for the winding angle of the shortest path. The distribution is a Gaussian with a variance consistent with Eq. (2). The estimate κwinding = 1.046 ± 0.004 that we get from fitting the data with Eq. (2) is in agreement with the value deduced from the fractal dimension.

Figure 2
figure 2

Variance of the winding angle against the lattice size Ly.

The analysis has been done for Ly ranging from 16 to 16384. The statistics are computed over 104 samples. The error bars are smaller than the symbol size. By fitting the results with Eq. (2), one gets κwinding = 1.046 ± 0.004. In the inset, the probability distribution of the winding angle along the curve is compared to the predicted Gaussian distribution, drawn in green, of variance with κ = 1.046 and Ly = 16384.

Full size image

Left passage probability

In the following, we work with chordal SLE. Therefore, one has to conformally map the original curves into the upper half plane. This is done using an inverse Schwarz-Christoffel transformation (see Supplementary Information online).

The shortest path splits the domain into two parts: the left and the right parts of the curve. The curve is said to pass at the left of a given point if this point belongs to the right side of the curve, see inset of Fig. 3b. For chordal SLEκ curves, Schramm has computed the probability of a curve to go to the left of a given point z = Reiϕ, where R and ϕ are the polar coordinates of z58. For a chordal SLEκ curve in , the probability Pκ(ϕ) that it passes to the left of Reiϕ depends only on ϕ and is given by Schramm's formula,

where Γ is the Gamma function and 2F1 is the Gauss hypergeometric function. We define a set of sample points S in for which we numerically compute the probability P(z) that the curve passes to the left of these points. To estimate κ, we minimize the weighted mean square deviation Q(κ) defined as,

where |S| is the cardinality of the set S and ΔP(z)2 is defined as , where Ns is the number of samples60.

Figure 3
figure 3

Left passage probability test.

(a) Weighted mean square deviation Q(κ) as a function of κ, for Ly = 16384. The vertical blue line corresponds the minimum of Q(κ) and the green vertical line is a guide to the eye at κ = κfract. The minimum of the mean square deviation is at κLPP = 1.038 ± 0.019. The light blue area corresponds to the error bar on the value of κLPP. We define the error bar ΔQ for the minimum of Q(κ) using the fourth moment of the binomial distribution. The error Δκ is defined such that Q(κ ± Δκ) − ΔQ = Q(κ) + ΔQ. We considered 400 points, regularly spaced in [−0.1Lx, 0.1Lx] × [0.15Ly, 0.35Ly] which are then mapped through the inverse Schwarz-Christoffel mapping into 72. (b) Computed left passage probability as a function of ϕ/π for R ∈ [0.70, 0.75] and κ = 1.038. The blue line is a guide to the eye of Schramm's formula (3) for κ = 1.038.

Full size image

For a lattice size of Ly = 16384, the minimum of the mean square deviation is observed for κLPP = 1.04 ± 0.02 as shown in Fig. 3. This value is in agreement with the estimate of κ obtained from the fractal dimension and the winding angle.

Direct SLE

The winding angle and left passage analyses are indirect measurements of κ. Therefore we also test the properties of the driving function directly in order to see if it corresponds to a Brownian motion with the expected value of κ.

As for the left passage probability, we consider the chordal curves in the upper half plane, starting at the origin and growing towards infinity. We want to compute the driving function ξt underlying the process. For that, we numerically solve Eq. (1) by considering the driving function to be constant within a small time interval δt, thus one obtains the slit map equation48,61,

We start with ξt = 0 at t = 0 and the initial points of the curve and map recursively all the points , i > 0, of the curve to the points through the map , sending to the real axis by setting and in Eq. (5). Re{} and Im{} are respectively the real and imaginary parts. In the case of SLEκ the extracted driving function gives a Brownian motion of variance κ. The direct SLE test consists in verifying that the driving function is a Brownian motion and compute its variance to obtain the value of κ. The variance should behave as .

We extract the driving function ξt of the shortest path curves using the slit map, Eq. (5). Figure 4a shows the variance of the driving function as a function of the Loewner time t. We observe a linear scaling of the variance with t. The local slope κdSLE(t) is shown in the inset of Fig. 4a. In Fig. 4b, we plot the mean correlation function C(τ) = 〈C(t, τ)〉t of the increments δξt of the driving function, where the correlation function is defined as,

One sees that the correlation function vanishes after a few time steps. The initial decay is due to the finite lattice spacing, which introduces short range correlations. But in the continuum limit, the process is Markovian, with a correlation function dropping immediately to zero. In the inset of Fig. 4b, we show the probability distribution of the increments for different t. This distribution is well fitted by a Gaussian, in agreement with the hypothesis of a Brownian driving function. From this result and the estimates of the diffusion coefficient computed for several lattice sizes, we obtain κ = 0.9 ± 0.2.

Figure 4
figure 4

Driving function computed using the slit map algorithm.

(a) Mean square deviation of the driving function as a function of the Loewner time t. The diffusion coefficient κ is given by the slope of the curve. In the inset we see the local slope κdSLE(t). The thick green line is a guide to the eye corresponding to κdSLE = 0.92. (b) Plot of the correlation C(t, τ) given by Eq. (6) and averaged over 50 time steps. The averaged value is denoted C(τ). In the inset are shown the probability distributions of the driving function for three different Loewner times t1 = 1.2 × 10−3, t2 = 3.7 × 10−3 and t3 = 9.95 × 10−3. The solid lines are guides to the eye of the form , for i = 1, 2, 3.

Full size image

We note that the numerical results obtained with the direct SLE method are less precise than with the other analyses and, therefore, characterized by larger error bars, as is well known in the literature51,53,59,62,63,64. The result we have obtained for κ is in agreement with the ones obtained with the fractal dimension, winding angle and left-passage probability.

We also extracted the driving function of the curves in dipolar space, i.e. defining the curves as starting from the origin and growing in the strip (see Supplementary Information online). We also obtained a value of κ consistent with the fractal dimension.

Discussion

All tests are consistent with SLE predictions. The numerical results obtained with the winding angle, left-passage and direct SLE analyses are in agreement with the latest value of the fractal dimension. Being SLE implies that the shortest path fulfills two properties: conformal invariance and domain Markov property (DMP). Thus, the agreement with SLE predictions lends strong arguments in favor of conformal invariance and DMP of the shortest path.

The DMP is related to the evolution of the curve in the domain of definition. Let us consider the shortest path γ defined in a domain , starting in a and ending in b. We take a point c on the shortest path different from a and b. Then if the DMP holds, one would have that

where γ[c, b] is the shortest path starting in c and ending in b in the domain except the curve γ[a, c], denoted as and and are the probabilities in the domains and respectively. One can classify the models as the ones for which DMP holds already on the lattice and the ones for which it holds only in the scaling limit. Many classical models, like the percolation hulls, the LERW, or the Ising model65 for example, belong to the first case. But some two-dimensional spin glass models with quenched disorder52,53 are believed to only fulfill DMP in the scaling limit. Our numerical results suggest that, for the shortest path, DMP holds at least in the scaling limit. Further studies should be done to test the validity of DMP on finite lattices.

The second result we can expect if SLE is established for the shortest path is conformal invariance. Conformal invariance, being a powerful tool to compute critical exponents, is of interest for the study of the shortest path. Conformal invariance, associated to Coulomb gas theory for example, could be useful to develop a field theoretical approach of the shortest path. There is no proof of conformal invariance of the shortest path, but our numerical results give strong support to this hypothesis. For example, the expression of the winding angle is based on conformal invariance and agrees with the predictions based on the fractal dimension. Also the left passage probabilities and the direct SLE measurements have been performed on curves conformally mapped to the upper half plane and gave consistent results. In addition, we obtained the same estimate of κ by extracting the driving function in chordal and dipolar space. However, even if the scaling limit would not be conformally invariant, our results suggest that one could still apply SLE techniques to the study of this problem, as some SLE techniques have also been used to study off-critical and especially non conformal problems66,67,68,69,70.

Analyzing the shortest path in terms of an SLE process would give a deeper understanding of probability distributions of the shortest path, allowing to compute more quantities, like for example the hitting probability distribution of the shortest path on the upper boundary segment71.

Methods

We generate random site percolation configurations on a rectangular lattice Lx × Ly with triangular mesh, where Lx and Ly are respectively the horizontal and vertical lattice sizes, in units of lattice sites. The sites of the lattice are occupied randomly with the critical probability . If the configuration percolates, we obtain the spanning cluster and identify the shortest path between the top and bottom layers using a burning method5,13,16. In short, we burn the spanning cluster from the bottom sites, indexing the sites by the first time they have been reached and stop the burning when we reach for the first time the top line. We then start a second burning from the sites on the top line that have been reached by the first burning, burning only sites with lower index. With this procedure, we identify all shortest paths from the bottom line to the top one. We randomly choose with uniform probability one of these paths. The results presented in the paper are for Ly ranging from 16 to 16384 and an aspect ratio of Lx/Ly = 1/2. We generated 10000 samples and discarded the paths touching the vertical borders.