Bridges in the random-cluster model

The random-cluster model, a correlated bond percolation model, unifies a range of important models of statistical mechanics in one description, including independent bond percolation, the Potts model and uniform spanning trees. By introducing a classification of edges based on their relevance to the connectivity we study the stability of clusters in this model. We derive several exact relations for general graphs that allow us to derive unambiguously the finite-size scaling behavior of the density of bridges and non-bridges. For percolation, we are also able to characterize the point for which clusters become maximally fragile and show that it is connected to the concept of the bridge load. Combining our exact treatment with further results from conformal field theory, we uncover a surprising behavior of the variance of the number of (non-)bridges, showing that these diverge in two dimensions below the value $4\cos^2{(\pi/\sqrt{3})}=0.2315891\cdots$ of the cluster coupling $q$. Finally, it is shown that a partial or complete pruning of bridges from clusters enables estimates of the backbone fractal dimension that are much less encumbered by finite-size corrections than more conventional approaches.


Introduction
Percolation is probably the most widely discussed and arguably the simplest model of critical phenomena. Due to a combination of conceptual simplicity and wide applicability which is a signature of a problem of very general interest, its many incarnations including bond and site percolation on a lattice as well as continuum and non-equilibrium, directed variants, have been the subject of thousands of studies [1]. In statistical physics and mathematics alike, the quest to understand aspects of the percolation problem has led to developments of powerful and beautiful new techniques [2,3]. While the problem is well defined and interesting for any graph or lattice, both finite and infinite, the understanding of non-trivial cases is most advanced in two dimensions. There, the scaling limit of critical percolation can be related to the Coulomb gas [4] and conformal field theory [5], leading to exact results for most critical exponents and certain correlation functions. More recently, a rigorous approach to conformal field theory was pioneered by Schramm who used a mapping introduced by Löwner to construct a way of generating conformally invariant fractal random curves, the Stochastic (or Schramm) Löwner Evolutions (SLEs) [6], for which a range of properties, including fractal dimensions, can be calculated exactly. In this context, Smirnov and co-workers used the concept of discrete analyticity to establish rigorously that the scaling limits of critical percolation [7,8] and the Ising model [9] on the triangular lattice are indeed conformally invariant, and cluster boundaries in these models converge to certain classes of SLE traces.
The random-cluster (RC) model was suggested by Fortuin and Kasteleyn as a natural extension of the (bond) percolation problem, noting that there was a class of models fulfilling the series and parallel laws of electrical circuits that also included the Ising model [10]. Given a graph G = (V, E), it assigns to a spanning subgraph (V, A ⊆ E) a probability mass (in the following also referred to as RC measure) [11]  The presence of the cluster-weight factor q K(A) distinguishes (1.1) from the percolation problem and prevents P p,q,G [A] from being a Bernoulli product measure; only for q → 1, where the model reduces to the percolation problem and hence edges become independent, this property is restored. Although the cluster weight q can be any non-negative real number, integer values of q are particular in that for them the partition function is very closely related to that of the q-state Potts model [12] (see Eq. (2.3) below). For lattice graphs in at least two dimensions, the model undergoes a percolation phase transition at a critical value p c (q) of the bond probability where, for sufficiently large q the transition becomes discontinuous [11]. On the square lattice, self-duality allows to deduce the exact transition point p c (q) = p sd (q) = √ q/(1 + √ q) [13], and the location of the tricritical point is known to be q c = 4, beyond which the transition becomes of first order [14].
The crucial importance of the RC description for the understanding of critical phenomena is through its expression in purely geometrical terms. Hence understanding the geometric structure of the (correlated) percolation problem (1.1) provides a geometric route to the understanding of the thermal phase transition of the Potts model. Correspondingly, significant effort has been devoted in particular to investigations of the structure of the incipient percolating cluster. While initially it was assumed that it was a network of nodes connected by essentially one-dimensional links, results regarding the conductivity of the critical cluster implied that, instead, the structure is better described by the more elaborate 'linksnodes-blobs' picture [15]. If one fixes two distant points A and B on the cluster, those bonds that have independent, non-intersecting paths to both A and B form the backbone of the cluster. It is essentially the part that would carry current in case a voltage was applied between A and B. The remaining cluster mass is in dangling ends. The backbone itself consists of singly-connected or red bonds that destroy percolation or, equivalently, destroy conductivity between A and B if cut (the links which meet in the nodes) as well as bonds in cycles that are multiply connected (the blobs). These subsets of bonds have fractal scaling with associated exponents d F for the cluster mass itself, d BB for the backbone, and d RB for the red bonds. Obviously one has d RB ≤ d BB ≤ d F , and it turns out that strict inequality holds in the generic case and hence asymptotically almost all of the bonds of the incipient percolating cluster are in dangling ends, and almost all of the mass of the backbone is in the blobs. In two dimensions, exact expressions for d F and d RB are available from Coulomb gas arguments [16], but d BB is only known numerically [17,18].
As we discuss here, the distinction of singly-connected and multiply-connected bonds first suggested for understanding the backbone structure [15] is a more generally useful classification of bonds in a cluster configuration. Removing singly connected bonds, which we call bridges, 1 generates an additional connected component, while this is not the case for multiply connected bonds or non-bridges. Indeed, this separation allows us to control the effect of a local edge manipulation (open ↔ closed) on the weight of the configuration as in Eq. (1.1). We will exploit this below to derive expressions for the expected number of bridges and related classes of bonds as well as the scale of fluctuations in the former. A further partitioning of bridges was recently suggested for percolation in Ref. [18]: branches are bridges for which at least one of the connected clusters is a tree, whereas the remaining bridges are dubbed junctions. This distinction is related to but not congruent with the links-nodes-blobs picture: most of the red bonds will be junctions, but so are some bonds in dangling ends as well. Cutting the branches, on the other hand, will remove some part of the dangling ends (the treelike ones), and branches are not typically found in the backbone 2 . As a result of this incongruence, the number of junctions on the percolating cluster grows as L dF and not proportional to L dRB , where L is the linear size of the system. In contrast to the class of red bonds, therefore, both branches and junctions are extensive subsets of the bridges. For uncorrelated percolation Xu et al. [18] find numerically, however, that removing all bridges, i.e., both branches and junctions, leads to tightly connected remainders, corresponding to the blobs, whose size grows as L dBB . As we will show below in Sec. 6, it suffices to remove the junctions only to see the same behavior more generally for the RC model.
As, by definition, the removal of a bridge leads to the generation of a new component and hence the breakup of an existing cluster, understanding the properties of bridge bonds is crucial also to the understanding of fragmentation phenomena in the framework of a lattice model. This connection was investigated in our recent letter [20], where we studied the fragmentation rate and kernel, and related the associated scaling exponents to the more standard critical exponents and fractal dimensions. In the present paper, however, we follow a different line of study. In particular, we derive a linear relation between the expected density of bridges and the density of bonds; we provide the asymptotic densities of bridges in the thermodynamic limit for the square lattice and study the scaling corrections in finite systems in general. We employ Monte Carlo simulations based on the Swendsen-Wang-Chayes-Machta algorithm [21,22] and a recent new implementation [23] of Sweeny's algorithm [24] to confirm these results for the square lattice and study the behavior in three dimensions. Building on these results for the densities, i.e., the first moment of the bridge distribution, we move on to studying the fluctuations, corresponding to the second moment. We reveal an unexpected finite-size behavior leading to a nonspecific-heat variance singularity for q ≤ 4 cos (π/ √ 3) 2 , a direct consequence of an intriguing interplay between bridges and non-bridges.
The rest of this paper is organized as follows. In Section 2 we introduce the edge classification, derive the bridge-edge identity in general, consider some of its consequences and present a numerical study for the special case of the square lattice, together with a detailed analysis of finite-size corrections. In Section 3 we use symmetries of the RC model in order to derive connections to other quantities of interest in the percolation literature, such as the short range connectivity. In Section 4 we introduce the concept of bridge load which allows us to analytically characterize the point of maximal bridge density for independent bond percolation. In Section 5 we derive a second-order bridge-edge identity that allows us to study the variance of the number of bridges. In Section 6 we compare various new and old methods to extract the backbone fractal dimension, with an emphasis on the involved finite-size corrections. Section 7 contains our conclusions. [16,19]. While these are actually backbone bridges, we use the term in a more general sense here to denote singly-connected bonds anywhere in the cluster configuration.

Edwards-Sokal coupling
Although in this paper we focus on the RC model itself, its relation to the Potts model is of relevance for the physical interpretation of the (geometric) percolation transition in terms of a thermal phase transition. In particular, we make reference to how certain observables of the Potts model translate into quantities in the RC language. The Potts model describes interacting spins, and assigns one of the q spin values in {1, . . . , q} to each vertex, attributing a configurational probability or Boltzmann weight to a spin configuration σ = {σ 1 , . . . , σ |V | } ∈ Q = {1, · · · , q} |V | . Here, we write 1 {σx=σy} for the indicator function, which equals 1 if σ x = σ y and 0 otherwise. Analogous to Z RC (p, q, G), the Potts partition function Z P (β, q, G) corresponds to the sum of the unnormalized weights. Importantly, in the second equality, we need to identify p = 1 − e −β . The relation between Eqs. (1.1) and (2.1) is established by augmenting phase space to include both, spin variables σ and subgraph variables A, leading to the joint probability that is also known as Edwards-Sokal coupling [25]. It has the crucial property that the marginal distribution on the spins alone reduces to the Potts weight (2.1), while the marginal on the subgraphs gives the RC weight (1.1). Similarly, one can derive the conditional probabilities for the spins given the bonds or vice versa [11]. As was already shown by Fortuin and Kasteleyn [10], the corresponding partition functions coincide up to a (trivial) multiplicative factor if we relate the parameters β and p via p = 1 − e −β , This coupling is the basis for the Swendsen-Wang cluster algorithm [21] and its generalization to noninteger q due to Chayes and Machta [22]. The above coupling and the identity (2.3) allow in particular to relate expectation values in the two ensembles. For instance, the internal energy in the Potts model is up to constants identical to the expected number of edges in the RC language, That is two vertices can only have unequal spin when they are not connected in the bond configuration and were assigned two different (random) colors. Summing over the edges shows that the nearest-neighbor connectivity defined as the sum is essentially the energy [11,26,27] Clearly for hypercubic lattices |E|/|V | is bounded (more precisely equals d, where d is the dimensionality).

Edge classification
Intuitively, (1.1) indicates that for q > 1 (q < 1) the model favors configurations with a larger (smaller) number of components, compared to the case q = 1. In this paper we investigate how the cluster structure is influenced by the cluster weight q, in particular with respect to the fragility of clusters. The relevance of individual edges for the connectivity of A ⊆ E depends on whether they are pivotal or non-pivotal. According to our definition, pivotal edges change the number of components upon removal/insertion, i.e., for a pivotal edge e we have K(A e ) = K(A e ) − 1, where A e (A e ) is the configuration obtained from A by opening (closing) e (note that necessarily one of A e and A e is equal to A). Hence for a pivotal edge there is no alternative path connecting the vertices incident to it. As discussed above, for the purposes of this paper open pivotal edges are denoted as bridges, while we call pivotal closed edges candidate-bridges. Likewise, we denote open (closed) non-pivotal edges as non-bridges (candidate-non-bridges). Let B, C, B, C denote the set of bridges, non-bridges, candidate-bridges and candidate-non-bridges, respectively. We stress that these sets depend explicitly on the configuration A. We denote the corresponding densities as B = |B|/|E|, C = |C|/|E|, B = |B|/|E|, and C = |C|/|E|.

Derivation of the bridge-edge formula
In this sub-section we show how, for arbitrary (finite) graphs, the expected densities of all the above mentioned edge-types B, C, B, and C can be related to the expectation of the edge density. We establish these identities by utilizing a differential equation known as Russo-Margulis formula [28,29], which applies to expectations with respect to product measures such as the probability density of the bond percolation model, corresponding to the q = 1 RC model. The formula intuitively quantifies the "response" to an infinitesimal change in the bond density. Using this approach to quantify the response in the number of connected components K allows us to study the bridge density. We then proceed and apply the idea to the RC model with generic q > 0, by exploiting the fact that the partition function of the model can be interpreted as a particular expectation in the percolation model.
First consider the percolation model, i.e., the measure (1.1) with q = 1, for which we denote expectations as E p,G [·]. In this setting, the Russo-Margulis formula reads d dp where X(A) : Ω → R is an arbitrary observable and the quantity δ e X is called the influence of e on X (or the discrete e-derivative of X) and is given by The graph specific geometric properties of the model under study are "encoded" in the influences. In order to simplify the subsequent analysis, we show in Appendix A the following bijection identity: (2.8) We now apply this identity to the cluster-number observable K, for which K(A) − K(A e ) = −1 if e is a bridge, e ∈ B(A), and 0 otherwise. Application of the Russo-Margulis formula (2.7) and Eq. (2.8) to X = K then yields d dp (2.9) This differential equation allows us to determine the cluster numbers E p,G [K] once the p-dependence of E p,G [B] is known. We remark that this relation was derived before in [30], however with a different target application in mind. As already mentioned, we can extend the above to the RC model and obtain more generally the bridge density E p,q,G [B] of the RC model itself. To start with, we note that the RC model partition function is Z RC (p, q, G) = E p,G [q K ]. In other words, the RC model partition function is the expected value of q K in the bond percolation model with parameter p. Hence we can apply the Russo-Margulis formalism to q K . Using the same reasoning as that for E p, On the other hand, direct differentiation of the partition function Z RC (p, q, G) yields Comparing both expressions (2.10) and (2.11), we derive the following relationship for the RC model (referred to as bridge-edge formula), which is of central importance for the following analysis (see [31] for some related identities developed along different lines). We remark again that this holds for any graph and for any p ∈ (0, 1) and q ∈ (0, ∞)\ 1 (the cases q → 0, 1 are discussed below). Due to the obvious identity E p,q, one immediately arrives at an analogous equation for the density of non-bridges, We remark that in principle one can use the Russo-Margulis formalism to study the (partial) p-derivative of expectations of other observables in the RC model: which shows the relevance of bridges in general for the RC model. Furthermore, let us consider a few direct consequences and applications of (2.12): 1. Firstly, notice that (2.12) recovers the correct results for the limit of uniform spanning trees, where all open edges are bridges. This tree model can be obtained in the RC model in the limit q, p → 0 for q/p → 0, which for the square lattice includes the critical line p sd (q) [11].
2. It is possible to establish an upper bound for E p,q,G [B] for q ≥ 1 and any p on any graph. Using general comparison inequalities for the RC measure [32], it is possible to show that wherep(p, q) ≤ p for q ≥ 1. Together with (2.12) this yields the desired upper bound for E p,q,G [B] As is shown below for the Ising case q = 2 with an exact solution this bound is tight in the limit p → 0, but clearly not so for p p sd .
3. Another interesting consequence follows when one recasts (2.12) to express the edge density in terms of the bridge density (2.14) This emphasizes the importance of bridges for q = 1 as the source of finite-size corrections in , a manifestation of the introduction of correlations (deviation from the product measure). Moreover it also explains, in the setting of the RC model, how thermal fluctuations in the number of edges vanish in the percolation limit q → 1, namely by a vanishing factor (q − 1). This has been a plausible assumption in the literature so far, but to our knowledge, relation (2.14) is the first exact statement, see for instance [27,33].

The square lattice
For the infinite square lattice Z 2 , the self-dual point p sd (q) = √ q/(1 + √ q) is critical 3 . For the Potts model at criticality, a mapping to an exactly-solved six-vertex model gives extensive analytical results [12,14]; it is assumed that by analytic continuation these carry over to general values of q. In particular, the critical internal energy density in this case is u βc(q),q,Z 2 = 1 + 1/ √ q and hence from Eq. (2.4) one for any graph G (analogous for probabilities), where of course the critical line is given by p c (q) = p sd (q) for Z 2 . It follows from (2.12) and (2.13) that one has We note that Eqs. (2.15) can be written in terms of p sd (q) and E q, The above expressions allow us interpret p sd (q) as the expected fraction of non-bridges among all open edges. Clearly, the remaining fraction 1 − p sd (q) of open edges have to be bridges. We note that because of e ∈ B ⇒ e ∈ A we have P p,q,G [e ∈ B|e ∈ A] = P p,q,G [e ∈ B]/P p,q,G [e ∈ A] and hence in particular P q,Z 2 [e ∈ B|e ∈ A] = 1 − p sd (q) and similarly P q,Z 2 [e ∈ C|e ∈ A] = p sd (q). Interestingly, an immediate consequence of (2.15) is that our RC result reduces in the percolation limit, , which recovers a result recently derived in Ref. [18].
Going beyond criticality, for the Ising model, corresponding to q = 2, we can use the available exact solution for the internal energy and the relation (2.4) to find the exact p-dependence of the densities of edges as well as bridges and non-bridges. This allows us to study, in an exact setting, the fragility as the bond-density/temperature is varied. As the corresponding expression for the internal energy is not very instructive we refrain from reproducing it here, see, e.g., Ref. [14] for details. Instead, we show in the upper panel of Figure 1 the exact asymptotic results for Z 2 . For comparison, we also show numerically estimated values of the density of bridges for other cluster weights (for details regarding the numerical procedure, see Sec. 2.6 below) in the lower panel of Fig. 1. As it is clear from Fig. 1, the density of bridges is not maximized at the critical point but for some p f (q) < p sd (q). We shall give an explanation for this phenomenon in terms of the relationship between nearest-neighbor connectivity and bridge density for the percolation case q = 1 in Sec. 2.5. Regarding the results for the Ising model, we note that it is also possible to produce exact expressions for the internal energy on finite lattices [34], resulting in corresponding exact results for the bridge densities on Z 2 L , the L × L square lattice with periodic boundary conditions. It is apparent from the lower panel of Fig. 1 that a characteristic of the critical point is the extremal . Moreover, the relation (2.12) suggests that the slope becomes singular due to a thermal singularity for q ≥ 2. Indeed, the occurrence of a singular slope can be explained by using (2.10) connecting E p,q,G [B] to the first p-derivative of the log (Z p,q ), which allows one to derive . (2.17) this latter relation together with (2.4) yields the asymptotic scaling [14] and L α/ν for all other values in q ∈ (0, 4], with multiplicative logarithmic correction proportional to log −3/2 (L) for the tricritical case q = 4 [26,[35][36][37]. Here, f is the free-energy density. We note that Coulomb gas arguments yield expressions for α/ν for q ∈ (0, 4] that predict α/ν > 0 whenever q > 2 [4], implying a finite slope for q < 2 even in the limit L → ∞. We refer the reader to Section 2.7 for more details about size dependent effects and the universality of the above arguments.

The percolation case
Returning to the case of general graphs, we see that the bridge-density formula (2.12) is singular for q = 1. To deduce the correct result in the percolation limit q → 1, we use L'Hôpital's rule to find Here, we have used that E p,q,G [N ] → p as q → 1. The l.h.s. being a density, which in particular must be non-negative, shows that the covariance between the number of components K and the density of edges N is non-positive. This is plausible, because adding edges can never increase the number of components.
In fact on more mathematical grounds, for the monotone case, q ≥ 1, this can even be proven, as it follows from the Fortuin-Kasteleyn-Ginibre (FKG) inequality [11], applicable to the RC model whenever q ≥ 1. Applied to K and N , the FKG inequality yields E p,q, , which in turn particularly applied to q = 1 shows that Cov p, Note that the covariance in (2.18) was analyzed for the percolation model in Ref. [38], however without establishing a connection to the bridge-density. There, the authors discuss the size dependence of Cov pc,Z d Z d L is the d dimensional hypercubic lattice with linear dimension L and periodic boundary conditions. The authors find that the quantity has a leading finite-size correction that scales as L 1/ν−d . In Sec. 2.7 we provide an explanation for this. in terms of the bridge density.
Before we proceed to discussing the behavior on finite lattices in the RC setting, we derive a symmetry relation for the density of bridges, valid for Z 2 , which for the special case p = 1/2 implies the previously shown result that E 1/2,Z 2 [B] = 1/4 for critical percolation on Z 2 . Let us naturally view the vertex set of Z 2 , denoted by V (Z 2 ), as the set of (ordered) tuples, whose entries are integers, and the set of edges of Z 2 , denoted by E(Z 2 ), is the set of pairs of vertices which differ in precisely one coordinate by ±1. Then we define a sub-box of Z 2 , denoted by n . This allows us to consider a sequence of sub-boxes (in n), that "approaches" Z 2 . It is not hard to see that each of the G n 's is planar, hence we can use Euler's formula for planar graphs, stating that for any A ⊆ . Red (hollow) and blue (solid) edges correspond to bridges and non-bridges, respectively. The dashed edges in the right panel are back-edges. Note that the depth-first search tree together with its back-edges can be decomposed into two chains, corresponding to the two back-edges, and any non-bridge is in at least one chain.
in any planar embedding of the graph (V n , A). Taking expectations on both sites and dividing by |E n |, we obtain Furthermore it can be verified that for any planar graph the number of components in the primal graph (V, A) equals the number of faces induced by the dual configuration A in the dual graph (V , A ), in other words we have applied to this case, F Gn (A) = K G n (A ), see e.g. [32]. Thus, we can re-write the above which uses the duality of bond percolation, i.e., P p,G (A) = P 1−p,G (A ), valid for any planar graph G.
We can now differentiate both sides with respect to p and apply the Russo-Margulis formula (2.9) to obtain and exploiting self-duality of Z 2 in the form of the plausible assumption that lim n→∞ Finally, for the particular choice of p = 1/2 we recover the (asymptotic) result E

Numerical analysis for finite lattices
In order to confirm the asymptotic results (2.15) for the square lattice and to study the finite-size corrections for finite lattices, we performed Monte Carlo simulations of the RC model on the L × L square lattice with periodic boundary conditions, in the range of continuous phase transitions 0 < q ≤ 4. We used a recent implementation [23] of (the Metropolis variant of) Sweeny's algorithm [24] for q < 1 and the Chayes-Machta-Swendsen-Wang algorithm [22] for q ≥ 1. We determined the number of bridges in a given configuration by means of the algorithm introduced in Ref. [39]. In contrast to the approach used in [18,40], this is a linear time algorithm applicable to any finite graph. It does not depend on a medial graph analysis, which facilitates the study of higher dimensional and non-planar systems. Here we briefly describe the algorithm and refer the reader to Ref. [39] for more details. The main idea is to construct a depth-first search (DFS) tree for any component in the spanning subgraph (V, A). Any edge in A that is not in a DFS tree is called a back-edge (there are precisely |A| − |V | + K(A) back-edges). Because a back-edge connects vertices already connected in a DFS tree, it follows that any back-edge is a non-bridge in (V, A); see Fig. 2 for an illustration of this situation. In order to be able to determine whether a given edge in a DFS tree is a bridge in (V, A), the algorithm introduces a chain decomposition of the graph, where for any back-edge there is exactly one chain (ignoring chains consisting only of one vertex). Now, if a given tree edge is part of such a chain, then it is in at least one cycle and hence a non-bridge. Finally, a careful construction that avoids the iteration over overlapping chains ensures that the algorithm classifies all edges into bridges and non-bridges in linear time.

Finite-size corrections
The bridge-edge relation (2.12) also allows us to understand the finite-size corrections to the bridge density. It relates E p,q,G [B] to the density of edges which, in turn, is related to the internal energy via (2.4). Finite-size corrections to the energy density u β,q,G close to a point of a second order phase transition, on the other hand, are widely studied and well understood from the theory of finite-size scaling. To formulate it, we need to associate a linear size L to G, such that |V | = L d , as is the case for a lattice graph with spatial dimension d. Then, the standard ansatz for the singular part f s of the free-energy density f = − log [Z P (β, q, G)]/(βL d ) as a function of the reduced temperature t = (T − T c )/T c , the external field h, and the leading irrelevant field v is given by [41] The internal energy u βc,q,Z d L can be expressed in terms of the first derivative of f with respect to t, evaluated at t = h = 0. Furthermore it is plausible to assume [41,42] that the non-singular part of f has no size-dependence and directly yields the value obtained in the infinite-volume limit. We hence expect the following scaling of the bridge density: (2.20) Here, we introduced the L-subscript to denote the expectation for the model on any lattice graph with linear dimension L, corresponding to the universality class of the q-state Potts model on Z d L . For the case of Z 2 L , the asymptotic density E q,∞ [B] is the one of Z 2 , as given by Eq. (2.15). For other 2D lattices such as the triangular and honeycomb lattice E q,∞ [B] can be derived in a similar way using well-known expressions for the critical internal energy there, see, e.g., Ref. [12]. The scaling form (2.20) itself, however, is universally valid.
We now turn to a more detailed analysis of the corrections for the case of the square lattice. We used the above bridge detection algorithm to numerically determine the size dependence of the bridge density for the critical RC model on Regarding the prefactor b of the leading term L 1/ν−d extracted from the fits, we find that it is negative for the bridge density, hence the asymptotic result 1/[2(1 + √ q)] is approached from below.
The corresponding constant for E q,Z 2 L [C] is positive such that the limit is approached from above. This holds for all values analyzed in q ∈ [0, 4]. Thus in illustrative words, the configurations on Z 2 L still "feel" the "extra" edges closing the lattice to form a torus and hence they are typically less "fragile" than corresponding configurations on Z 2 . Interestingly, when one considers the finite-size dependence of E q,Z 2 L [N ] one finds an excess of edges (relative to the asymptotic result for Z 2 ) for q > 1 and a shortfall for q < 1. Now, due to the fact that E p,q,G [N ] = E p,q,G [B] + E p,q,G [C], we have that in particular the amplitudes of the finite-size corrections of E q,Z 2 L [B] and E q,Z 2 L [C] cancel each other for q = 1, whereas the modulus of the amplitude for E q,Z 2 L [B] is larger (smaller) than the one for E q,Z 2 L [C] for q < 1 (q > 1), respectively. We omit the details on the size dependence of E q,Z 2 L [C], as no new mechanism appears and the above results on E q,Z 2 L [B] can be easily adapted to this case. Following Xu et al. [18], we extend our study of finite-size corrections on the torus Z 2 L from bridges to the larger class of "type-1 edges", that is edges that that have both of their associated loop arcs in the same loop in the associated loop configuration. While on a planar lattice all such bonds are bridges, this is not the case for Z 2 L , where also non-bridges can be of type 1. Such edges are hence called pseudobridges [18,40]. If we denote by L 1 (A) and L 2 (A) the set of (open) edges in (V, A) that have both of their associated loop arcs in the same loop (type-1 edges) or in two different loops (type-2 edges), respectively, we have that (2.21) The second equality follows from B ⊆ L 1 , that is any bridge has necessarily both loop arcs in the same loop (as it cannot enclose any face). Equivalently, we have E p,q,Z 2 , where P(A) ≡ |e ∈ E : e ∈ L 1 (A), e ∈ C(A)| /m, i.e., the density of pseudo-bridges in A, and 1 denotes the fraction of edges that are of type 1. In order to understand the general size dependence of E q,Z 2 L [ 1 ], we 12 determined the number of type-1 edges in our numerical simulations and fitted a finite scaling ansatz E q,Z 2 L [ 1 ] = a + bL − 1 + cL − 2 to our Monte-Carlo data. The resulting estimates are consistent with a = 1/[2(1 + √ q)] and 1 = 2 − 1/ν as well as 2 = x 2 . In a second step, we also performed fits with c fixed to its exact value and 1 and 2 fixed to the values implied by the Coulomb gas for the identifications 1 = 2 − 1/ν and 2 = x 2 , such that b and c were the only free parameters in a now linear fit. As it is apparent from the parameters collected in Table 1, the quality of these fits is very good, even including the smallest system sizes L = 4, 8, 16, 32, strongly suggesting the following asymptotic form for the square lattice where again E q,Z 2 [B] = 1/[2(1 + √ q)]. Figure 5 shows our numerical data for the two cluster weights 0.5 and 2, together with the corresponding best fits.
It is interesting to compare the general form (2.22) to the results found in Ref. [18] for the percolation case, corresponding to q → 1. There, finite-size corrections to the bridge density on the torus Z 2 L were found to decay with exponent −x 2 = −5/4. While this is numerically consistent with our findings since 1/ν − d = −x 2 for percolation [17,[44][45][46], it is clear from the present analysis that the relevant exponent describing finite-size corrections to the bridge density is 1/ν − d. Indeed, the values of −x 2 and 1/ν − d strongly differ for q = 1, cf. Fig. 4. Moreover, for the square lattice with periodic boundary conditions Xu et al. [18] showed rigorously that P 1/2,1,Z 2 L [e ∈ L 1 ] = E 1/2,1,Z 2 L [ 1 ] = 1/4, independent of L. This is consistent with our general form (2.22) for the percolation case with 1/ν − d = −x 2 if the amplitudes b and c cancel. Indeed, this is what we find from the fit data in Table 1, which clearly support the statement that b = −c as q → 1, as required by the rigorous finding that E 1/2,1,Z 2 L [ 1 ] = 1/4. Hence, it is a particular cancellation of finite-size effects that leads to the size-independence observed for percolation. Further, it is known that for q → 0 with p = p sd (q), one recovers the uniform spanning tree model, for which clearly no pseudo-bridges exist, that is one can easily verify that, on Z 2 L , E q,Z 2 L [ 1 ] → 1 2 − 1 2 L −2 for q → 0, which is in agreement with (2.22), due to 1/ν → 0 and the numerical observation (see Table  1) that a → −1/2 and b → 0 for q → 0. We remark that the fitting estimates for q close to one need to −0.1410(7) 0.504 (7)  Q denotes the quality of fit [43] or confidence level, that is the probability that χ 2 would exceed the empirical value, under the assumption that the imposed statistical model is correct. In what follows we refer to fits with confidence level ≥ 10% simply as "good" fits.
be treated with caution, as here the two exponents 1/ν − 2 and −x 2 come very close to each other, and hence it is numerically very difficult to distinguish the two contributions and the corresponding constants a, b. However, outside a suitable "safety"-window around q = 1, it appears that both a, b are increasing with q.
We close the analysis of size-dependent effects by establishing a relationship for percolation between pseudo-bridges and a particular arm event which shows the observed L −x2 decay. A pseudo-bridge is a non-bridge that resides on a cross cluster that winds simultaneously along both directions on the torus [18]. Furthermore, any pseudo-bridge has both associated loop-arcs in the medial graph in the same loop. In other words, a pseudo-bridge is a non-bridge on a cross cluster that is pivotal to the existence of the cross cluster, i.e., if A is a configuration that contains a cross cluster, then removing a pseudobridge e implies that A e has no remaining cross cluster. Therefore pseudo-bridges are in a certain sense the "bridges" of the cross cluster. We can translate this into a (polychromatic) four-arm event [3,47] as follows: Any edge, say e, that is pivotal to the cross cluster has the property that starting from e two paths of open edges extend to a distance L/2 (here we use the toroidal geometry and in particular the absence of boundaries). Since e is pivotal for the cross configuration, there must also be two paths of closed (dual) edges starting at the dual vertices associated to e and extending to a distance L/2, which precisely ensure that there can be no alternative, e-avoiding, path that would preserve the cross configuration property upon removal of e. Thus what we constructed is a particular instance of an polychromatic four-arm event, of alternating primal and dual "color", in the annulus centered around the edge e with outer radius L/2. The probability of such an event decays for large L as L −x (P) 4 , where x (P) 4 is the poly-chromatic four-arm exponent, which equals 5/4 for critical percolation 4 . Using the fact that x (P) 4 = x 2 , where the latter is the two-arm exponent as defined for instance in [46] and referred to by Xu et al. in [18] we have the appearance of the exponent x 2 .
We note that the notion of pseudo-bridges and the corresponding scaling analysis given above is specific to the case of locally planar lattices. The concentration on the square lattice was for practical reasons only, however, and we expect the same results, only with different amplitudes, for the cases of other lattices such as the triangular or honeycomb cases with periodic boundary conditions.

Other types of pivotal edges
So far our analysis has focused on open edges, that is bridges and non-bridges. Surprisingly, the analysis of closed edges turns out to be very interesting too, in that it allows us to link several seemingly unrelated quantities, studied in the literature in different contexts, to the density of bridges. Moreover the study of closed edges provides a further probe for the cluster structure, for instance a large density of candidate non-bridges suggests that the clusters are more likely to self-entangle than to overlap with other clusters. In order to analyze the corresponding densities consider a given closed edge e = (x, y). Clearly we have that when x y then e is a candidate bridge (and vice versa), i.e., a candidate bridge corresponds to a pair of nearest-neighbor vertices that are not connected. This observation allows us to state the following almost trivial but useful identity that holds for arbitrary e = (x, y) ∈ E: In fact, it is possible to relate the r.h.s. above to P p,q,G [e ∈ B], the probability that e is a bridge. To see this, note that for fixed e = (x, y) any configuration A ⊆ E that has the property x y can be related one-to-one to the configuration A + e, in which e is a bridge. Furthermore, as the insertion of e reduces the number of components by one, it follows that P p,q,G [A] = q(1 − p)P p,q,G [A e ]/p which in turn implies (3.1) Hence we have P p,q,G [e ∈ B] = q(1−p) p P p,q,G [e ∈ B], and therefore using the bridge formula (2.12) and the definition (2.6) of the nearest-neighbor connectivity E we obtain the general result The corresponding asymptotic values for the square lattice follow immediately: As a result of the relation to the edge density, we again conclude that the leading finite-size corrections are given by L 1/ν−d . These findings are in line with those of Hu et al. in Ref. [27], who find for percolation that for L → ∞ on Z 2 L one has E p,Z 2 L [E] → 3/4. On the other hand, it was shown in Ref. [18] that the asymptotic bridge density for critical percolation on Z 2 L is 1/4 (which also follows from our result (2.12)). These two results are hence clearly consistent with the relation (3.1). Finally from (3.2) and (3.3) one concludes that which is consistent with the results of Ref. [27].

Maximal bridge density for percolation
It is apparent from Fig. 1 that the expected bridge density E p,q,Z 2 L [B] is maximal for p f (q) < p sd (q) for any q ∈ (0, 4] [20]. The reason for this is a priori not clear and here we focus on the percolation case q = 1 and obtain an expression for d dp E p,G [B], that is for general graphs G, to determine its zero(s) one of which must correspond to the location of the sought maximum. In principle we could again apply the Russo-Margulis formula directly to |B| to obtain the required derivative. However, there is an alternative, more geometric, approach. In Ref. [44], Coniglio showed that the p-derivative of P p,G [x ↔ y] for two vertices x, y ∈ V (not necessarily neighbors) satisfies the following differential equation: where λ x,y (A) is the number of bridges on any 5 self-avoiding path in A connecting x and y. We emphasize that for a given edge (x, y) ∈ E the quantity λ x,y can be non-zero, even if (x, y) / ∈ B(A). This is because λ x,y counts the number of bridges in the spanning subgraph (V, A) that would disconnect x from y upon removal. As an example, consider vertices b and e in the left panel in Figure 2. Even with (b, e) not being part of the considered spanning subgraph, we have λ b,e = 2, because edges (in fact bridges) (b, c) and (c, f ) are essential for the connectivity of b and e. Note that bridge (a, b) is not essential for the connectivity of b and e as its removal does not disconnect b and e. Note that relationship (4.1) also follows straightforwardly from the Russo-Margulis formalism applied to the indicator function of the event x ↔ y. As a result of the bijection relation (3.1) we can transform (4.1) into a differential equation for the bridge density: We emphasize that now λ x,y is only evaluated for nearest-neighbor pairs, that is for  Figure 2, here edges (a, b), (b, c), (c, f ) have bridge loads 2, 3, 3, respectively. It is evident that in a certain sense ρ e generalizes the bridge density and encodes how relevant a given bridge is for the connectivity at large. We emphasize that one only has to count connected neighbor pairs. We finally obtain: where we definedρ(A) ≡ 1 |E| e∈E ρ e (A). Therefore we find that the l.h.s. in (4.4) vanishes for p = p f , where p f is the solution of the following equation (4.5) We remark that p f is a graph dependent quantity and we suppress the explicit dependence. In Fig. 6 we illustrate this condition for a small system size and the case G = Z 2 L . Moreover, under translational invariance (which is given, e.g., for Z 2 L ), (4.5) is equivalent to for arbitrary e ∈ E. This has a nice intuitive interpretation: E p,G [ρ e ] counts the typical number of nearest-neighbor pairs at the interface between two clusters that are glued together only by the presence of bridge e. There are two effects at work here: increasing p typically increases the size of the clusters and therefore possibly also the number of nearest neighbors on the interface. Additionally, however, increasing p will eventually create additional links between the two clusters, in which case e ceases to be a bridge, such that this effect decreases P p,G [e ∈ B]. Now, the maximal fragility (maximum of P p,G [e ∈ B]) is attained for the bond density p f , where the expected number of closed edges between two clusters glued together by e equals precisely the probability that e is a bridge. Increasing p further strengthens the cohesion between the clusters, and hence the probability of e being a bridge starts to decrease afterwards.
In what follows we restrict the analysis to the case G = Z 2 L in order to extract the asymptotic value of p f for Z 2 . Consider Figure 6, where we numerically confirm condition (4.5) and estimate p f = 0.4056 (5)  : Numerically estimated density of bridges and bridge load for percolation on Z 2 L for L = 16 and L = 128, corresponding to the smaller and larger markers, respectively. Additionally, the inset shows our estimates for p f for system sizes from L = 8 up to L = 128. The estimates for the location of the crossing "saturate" to a value 0.4056 (5), and therefore we conclude that p f = 0.4056(5) for Z 2 , according to (4.5). We point out that we used a naive depth-first search traversal based algorithm for determining the bridge load, so we could not study systems with L > 128 in practise. The lines in the figure are obtained by interpolating the data points using splines.
for Z 2 . We note that, albeit the maximum of P p,Z 2 [e ∈ B] is attained for p = p f , the bridge load E p,Z 2 [ρ e ] increases further and is apparently maximized for p = p sd (1) = p c = 1/2. This nicely reflects the intricate cluster structure at criticality. One might wonder how the bridge density can decrease but the overlap continues to increase. This can happen when in most (probable) configurations e is not a bridge, but for configurations where e is a bridge one has typically a very large self-entanglement of the clusters attached to the two ends of e, outweighing the effective decrease of P p,Z 2 L [e ∈ B]. This effect is most drastic for p = p c . We note that Coulomb gas arguments predict that the l.h.s. of (4.4) remains finite at p = p sd (1) = 1/2, which clearly implies that E p,Z 2 [ρ e ] remains finite at p = 1/2. More precisely we can show that the peak of E q,Z 2 L [ρ e ] scales as L 2−2x2 for q ∈ [0, 4], where in particular x 2 = 5/4 for percolation [see the discussion in the next section and in particular (5.3)]. Unfortunately, so far we have not been able to establish an exact value for p f (1) let alone in the correlated setting q = 1 for Z 2 . This is a challenging problem we wish to address in a forthcoming study.

Bridge fluctuations
As we have seen, the bridge-edge formula (2.12) together with the theory of finite-size scaling, provides a rather complete understanding of the critical bridge density. Higher moments of the bridge distribution can be discussed with similar techniques as we will show now for the example of the variance. Naively, one might expect that this variance Var p,q,G [|B|] only depends on the fluctuations of N . In this case the critical variance would follow Var q,L [|B|]/L d ≈ L α/ν , which, in two dimensions, in turn would imply a divergence with L for q ≥ 2 and a "saturation" for q < 2 [46]. As we will show however, the story is not quite as simple.
To work out the fluctuations of the bridge density, we apply the Russo-Margulis formula to the second derivative of the partition function Z p,q,G ≡ Z RC (p, q, G), which we then in turn equate to the expression one obtains by explicit differentiation. To start with, we have, using (2.11), Let us now focus on the second term: We can split the inner sum into two sums corresponding to e ∈ B(A) and e ∈ C(A), of which we first consider the former A few comments are in order: If e = f we clearly have that e = f / ∈ B(A e ) hence we recover the expected bridge density. For e = f it is important to observe that removing a bridge e cannot influence the pivotality of an occupied edge f , in particular when f ∈ B(A) then also f ∈ B(A e ), whenever e is a bridge in A. The above can therefore be re-written as For the non-bridge sum of (5.1) note that only summands with e = f contribute, because otherwise if e ∈ C(A) then trivially f = e / ∈ B(A), B(A e ). Furthermore, if e = f such that e ∈ C(A) and f ∈ B(A) then we have also a vanishing contribution because deleting a non-bridge cannot change the fact that f is a bridge. Thus, there can only be a contribution for a configuration A such that e ∈ C(A) and f ∈ C(A) as well as both edges are in one cycle and deleting e will destroy the cycle and hence cast f into a bridge in A e . Moreover, both edges e and f must only be in one cycle that has no edge overlap with other cycles (imagine two clusters glued together in parallel by the edges e and f , hence in particular there are no additional links between the two clusters). Write e A ⇔ f for the above event involving edge e and f . We obtain for the second term: We therefore obtain eventually for ∂ 2 p Z p,q : On the other hand one can show by explicit differentiation and after some straightforward but tedious algebra that one has: Equating both expressions for ∂ 2 p Z p,q,G with subsequent rearranging yields Using the bridge edge identity (2.12) we can simplify the above to We emphasize that (5.2) is an exact result valid for any p and q as well as any graph G, thus it has the same status as the bridge-edge identity (2.12). Furthermore, we remark that for e = (u, v) and f = (x, y) we have by a bijection argument, similar to the one used in Sec. 3 to derive the relationship between the bridge density and nearest-neighbor connectivity: Here P p,q,G [(u, v) (x, y)] for (u, v), (x, y) ∈ E is the probability that the two nearest-neighbor pairs belong to two different clusters, such that the two distinct clusters each contain one vertex from {u, v} and one vertex from {x, y} (see Fig. 1(c) in Ref. [48]).
We pause for a moment, and remark that completely analogous arguments can be used to obtain an expression for the p-derivative of E p,q,G [B], studied in Section 4, in terms of the probabilities of events e ⇔ f . A consequence of this is that we can obtain an explicit expression for E p,q,G [ρ]: In particular, this shows explicitly that the bridge load can increase further beyond the point p f (G), for which E p,G [B] is maximal. Moreover, as we show below, we find that at criticality, in two dimensions, , which together with x 2 = 5/4 for the case q = 1 in two dimensions, shows that the peak at p = 1/2 in Figure 6 stays bounded as L → ∞.
Coming back to the study of fluctuations of |B|, we focus in what follows on the self-dual line for the RC model on Z 2 L and study the continuum limit. As shown by Vasseur et al. [48] one has for two pairs of neighboring vertices (u, v) and (x, y) at distance r the following asymptotic for large r This follows from the construction of a four-leg watermelon event, due to the four hulls propagating from the neighborhood of (u, v) to the neighborhood of (x, y), which are associated to the two clusters 20 involved. We hence expect the following asymptotic behavior where α(q) and β(q) are two q dependent constants.
Inspecting the form (5.2), we hence see that the fluctuation of the bridge density has two contributions, one related to the variance of the edge density that scales as L α/ν at criticality, and another one related to the above mentioned watermelon event with scaling proportional to L 2−2x2 . In order to decide which contribution constitutes the dominant leading scaling, it is useful to recall some relevant exact results from the Coulomb gas approach to two dimensional critical phenomena [4]. One finds that α/ν = 4 − 12/g and x 2 = 2 − (4 − g)(4 + 3g)/8g, where for critical 2D RC models one has to relate the Coulomb gas coupling g to the cluster weight q via q = 2 + 2 cos (gπ/2) [46]. A small calculation shows then that L α/ν constitutes the leading contribution to the system size scaling of the variance whenever q ≥ 1. In particular, for q > 2, where α/ν > 0, this translates into a divergence of Var q,Z 2 L [|B|]/L 2 with L, resembling the specificheat singularity in the same cluster weight regime. On the other hand, for q < 1 the leading term is proportional to L 2−2x2 . Interestingly (and somewhat unexpected), one can explicitly verify that for q <q = 4 cos 2 (π/ √ 3) = 0.2315891 · · · the variance of |B| diverges with the exponent 2 − 2x 2 > 0. This particular valueq follows from the condition 2 − 2x 2 (q) = 0 ⇔ x 2 (q) = 1, which, using the above expressions, yields a quadratic equation forg = g(q) (orq): 8g = (4 −g)(4 + 3g), with the two solutions g = ±4/ √ 3. Taking only the positive solution we find thatq = 2 + 2 cos (2π/ √ 3) = 4 cos (π/ √ 3) 2 . The situation is summarized in the dotted and solid lines of Fig. 7, showing the Coulomb gas values of α/ν and 2 − 2x 2 , respectively. We also analyzed the variance of the bridge density numerically. The fitting functions and resulting parameters are summarized in Table 2, and the corresponding parameter estimates are indicated by the symbols in Fig. 7. Clearly, we find excellent agreement with the expectations from Eq. (5.2) discussed above. For the marginal valueq = 0.2315891 · · · we expect a logarithmic divergence, and we indeed find the corresponding form to yield the best fit to our simulation data.
It remains to discuss the percolation case q = 1 where the bridge-edge identity (2.12) becomes singular and the above derivation needs to be revisited. Here, we derive the singular behavior based on another bijection argument that allows us to harness recent results on logarithmic observables emerging from a careful analysis of the appropriate logarithmic conformal field theory description of critical percolation [27,48]. Before we do so, note that in order to extract the asymptotic scaling of the variance it suffices to study the (co-)variance η e,f ≡ P p,q,G [e ∈ B, f ∈ B] − P p,q,G [e ∈ B]P p,q,G [f ∈ B], which relates to Var p,q,G [|B|] via the well known identity To start with, recall that we consider critical bond percolation on Z 2 L , i.e. q = 1 and p = 1/2 and write P[·] for P 1/2,1,Z 2 L [·]. Furthermore, note that Z 2 L is a transitive graph, and hence none of the following events depend on the explicit edge or vertex, used in the arguments. Now, fix two edges e = (x 1 , y 1 ), f = (x 2 , y 2 ) that are distance r L apart. Consider the event {e ∈ B ∧ f ∈ B}. All configurations contributing to this event can be further sub-divided into two events, depending on whether e and f belong to the same connected component in (V, A) or not. Denote the two events by Ω 1 and Ω 2 , respectively. Choose  Table 2.
a configuration A that belongs to Ω 2 , i.e. the two edges e and f are bridges in (V, A) and belong to two different connected components. The crucial point is that we can relate A one-to-one to A − {e, f }, a configuration where x 1 , y 1 , x 2 , y 2 belong to four different clusters. Denote all configurations in which the four vertices belong to four different components byΩ 2 (this is a event). Due to the choice of q = 1 and p = 1/2 we have that P[Ω 2 ] = P Ω 2 . Let us now consider the event Ω 1 , i.e. the set of configurations for which e and f belong to the same component and both e and f are bridges. Now, because both edges are pivotal we can relate any such configuration A ∈ Ω 1 one-to-one to a configuration where x 1 and y 1 as well as x 2 and y 2 are disconnected, and for which x 1 , y 1 , x 2 , y 2 belong to three different clusters, of which one cluster contains one vertex of {x 1 , y 1 } and one of {x 2 , y 2 }. Denote the corresponding eventΩ 1 . Note that any configuration A ∈ Ω 1 , where e and f are bridges belonging to the same cluster, must yield 3 disconnected clusters in A = A − {e, f }. This is because the alternative case of 2 disconnected clusters in A would imply that e and f are in a cycle in A, which is obviously a contradiction. As before we have P[Ω 1 ] = P Ω 1 . Finally, the probabilities P Ω 1 and P Ω 2 were studied 6 in [48] in the framework of a corresponding logarithmic conformal field theory. We note that in the continuum limit both probabilities  [48]: where x 2 = 5/4 is the two-arm exponent for critical percolation and a, b, a , b and c are constants. Our numerical analysis confirmed this scaling, as shown in Table 2 and Figure 7. Thus considering the variance of bridges for critical percolation yields yet another manifestation of the underlying logarithmic conformal field theory [48]. We remark that the same logarithmic multiplicative corrections are expected  [|B|]. The two rightmost columns show the exact values obtained from the Coulomb gas mapping. The exponent values in bold face indicate the asymptotically dominant behavior. For the cluster weights q = 1.5 and q = 0.5 we also performed a fit to the form a + bL κ , which yielded slightly worse results, due to the proximity of the two (negative) exponents. In all cases the quality-of-fit Q was at least 5%.
in higher dimensions [49]. More precisely, one expects even the same scaling form, that is where we used the well known fact that x 2 = d − 1/ν holds for critical percolation [44]. It would be interesting to to verify this theoretical prediction numerically.
Before we proceed to the next section we briefly discuss the finite-size scaling of Var q,Z 2 L [|C|]/|E|. This can be worked out by observing that |C| = |A| − |B| which in turn implies the general identity Additionally, one can verify that (see [11] or by explicit differentiation of Z p,q ) An immediate consequence of these observations is that Var q,Z d L [|C|] is governed by the same finite-size scaling as Var q,Z d L [|B|] and we therefore omit further details.

Bridge-free clusters
In concluding the present study, we return to the discussion of the structure of the incipient percolating cluster outlined above in the introduction. In the 'links-nodes-blobs' picture of Ref. [15] the attention is focused on the backbone of the cluster as the skeleton ensuring long-range connectivity. Bridges on the backbone are singled out as the bonds carrying full current for the case of an external potential difference being applied, hence their name 'red bonds'. Due to the importance of these structures for a number of applications of percolation theory, significant effort has been invested in the determination of the associated fractal dimensions d BB and d RB of the backbone and the red bonds, respectively [50]. Grassberger (Bonds) Figure 8: Effective backbone scaling exponents according to Eq. (6.1) for the q = 1 RC model extracted from the scaling of the backbone sites, the backbone bonds (data from Ref. [50]), the junction-free clusters and the bridge-free clusters (our data), respectively.
known exactly from Coulomb gas and further arguments [4], there is no exact expression for the fractal dimension of the backbone.
Clearly, the notions of backbone and red bonds appear to depend on the concept of long-range conductivity, and it is not obvious how they could be defined in terms of the more local connectivity structure of bridges and non-bridges. The concepts of junctions and bridges introduced in Ref. [18] are related, but not commensurate to the idea of separating bridges in the backbone (red bonds) and bridges in the blobs. We expect that, away from q → 0, most bridges in the backbone are junctions as the 'linksnodes-blobs' picture tells us that the backbone has cycles. Still, there are junctions in the dangling ends and there are branches in the backbone. Also, as Xu et al. established for percolation in Ref. [18], branches and junctions are asymptotically finite fractions of the bridge set and hence of the edge set. It is therefore interesting to study the fractal dimensions of the related edge classes. Removal of branches from critical configurations will typically only shave off the last part of dangling ends, and we do not expect this to alter the fractal dimension of the incipient percolating cluster. This is indeed what Xu et al. report [18] and we find the same to be true also for the RC model. On the other hand, removing the junctions means deleting most of the red bonds and hence breaks down the percolating cluster into individual blobs. Their size scales with the backbone fractal dimension as ∝ L dBB . The same holds true for the case of bridge-free clusters where both, branches and junctions, have been removed. Studying these sets hence provides an alternative route to the determination of d BB .
To investigate the utility of this approach for the determination of d BB , we determined the scaling of the junction-free and bridge-free clusters for the percolation case q = 1. If we consider the effective, system-size dependent fractal dimensions [50] d BB,eff ≡ log[N (2L)/N (L/2)] log 4 , (6.1) where N (L) denotes the number of sites in the junction-free or bridge-free clusters at size L, respectively, we can easily compare corrections to scaling for the different approaches. This is illustrated in Fig. 8, comparing the exponents extracted from the scaling of the junction-free and bridge-free clusters to the effective exponents found from the scaling of the backbone itself [50]. Clearly, the scaling of bridgefree clusters is significantly less affected by scaling corrections than the scaling of the backbone itself. Compared to this, junction-free clusters show slightly increased corrections, however interestingly with a correction amplitude of the opposite sign. We remark that our estimate for percolation d BB = 1.643(1) is consistent with existing literature values, such as d BB = 1.64336(1), 1.6434(2), 1.6432 (8) in [17,18,50], respectively. As this is a small-scale study not using the particularly efficient algorithms available for the uncorrelated percolation case, the reduced scaling corrections do not immediately lead to an improved estimate. Instead, we focus on the estimation of d BB for a wide range of q values for which no estimates have been previously reported [17]. Whereas for the determination of d BB directly from the backbone corrections to scaling need to be carefully taken into account, it turns out that the bridge-free clusters provide a scaling route less encumbered with such problems. As is evident from the results collected in Table 3, over the full range of q values we achieve excellent fits without the inclusion of scaling corrections. : Estimated backbone fractal dimension using the method described in [18]. For comparison also values from [17] are shown. The disagreement for q = 4 is most likely due to known logarithmic corrections.

Discussion and outlook
As we have seen, the natural classification of edges into bridges and non-bridges provides an understanding of how the cluster weight influences the connectivity in the random-cluster model. For any q = 1 we established that both bridges and non-bridges have expected densities that are linearly related to the overall density of open edges, such that there is always a non-zero fraction of both edge types. The derivation is based on an application of the Russo-Margulis formula to the partition function of the model, a connection that allows to express analytical derivatives in terms of combinatorial and geometric quantities. For the percolation limit q → 1 the linear relation does not hold and rather transforms into an intuitive identity relating the bridge density to the covariance between the number of connected components and open edges. The latter identity explains previously observed finite-size corrections for this covariance in the percolation model [38]. Thus, the established bridge-edge formula allows us to connect the finite-size corrections of the densities of bridges and non-bridges to corrections in the energy density governed by the exponent d − 1/ν. The lack of finite-size effects in the density of the so-called type-1 edges for the special case of a 2D lattice with periodic boundary conditions, corresponding to a torus graph, can be understood from a cancellation of correction terms, and we clarified the character and origin of correction terms with exponents d − 1/ν and −x 2 , respectively [18].
By studying the dependence of the bridge density on the bond occupation p we observe that its maximum is attained for p f (q) strictly below the corresponding critical point p sd (q). We explained this observation for the percolation case by studying the newly introduced bridge load of an edge. For percolation on Z 2 L we estimate p f = 0.4056 (5). Our numerical study of this quantity was limited by the computational cost of determining the bridge load for each edge, which we were not able to calculate with linear (worst-case) computational effort for all edges. However, we remark that it is possible to express a condition for p f in terms of easily computable covariances and moments (in linear time), which however requires to estimate the difference of quantities of order O(L d ) to a precision of order O(1), which implies a statistical inefficiency. The construction of a more efficient algorithm for the determination of the bridge load constitutes an interesting open problem.
Having established the expectation values of the bridge and non-bridge densities, we turned to a study of the fluctuations. Again using the Russo-Margulis formalism we established a second bridge-edge identity relating the variances Var p,q,G [|B|] and Var p,q,G [|C|] to the variance of the bond density and an extra term relating to a specific four-leg watermelon event [48]. For the special case of the two-dimensional model, we use results on the scaling of this event derived from conformal field theory [48] to predict a singularity in both variances for q ≤q ≡ 4 cos 2 (π/ √ 3) = 0.2315891 · · · , an effect completely absent in the fluctuations of the overall number of open edges. Numerical simulations confirm these predictions to a high precision. The limiting percolation case required a separate analysis, and here we are able to show that the finite-size corrections to the density of bridges for critical percolation have logarithmic multiplicative corrections, a direct consequence of the underlying logarithmic conformal field theory. It would be interesting to study these fluctuations for higher dimensional systems to see whether such a singularity appears generically for some q < 1, and how this threshold value depends on dimensionality. It follows from hyperscaling, νd = 2−α, and the relation 1/ν = d R = d−x 2 [16] valid only for percolation (q = 1) that in this case α/ν = d − 2x 2 , irrespective of dimension d. Hence, q = 1 is a marginal case in all dimensions below the upper critical one, yielding a multiplicative logarithm in the variance of the number of bridges. It would be interesting to see whether d − 2x 2 > α/ν for q < 1 in 3D also (see [46] and [49] for two definitions of x 2 that can be generalized to d ≥ 3), thus hinting at an interesting physical/geometrical difference between the (anti-monotone) q < 1 and (monotone) q > 1 regimes.
Removing bridges leads to fragmentation of clusters as discussed recently in Ref. [20]. Removing all bridges decomposes the percolating cluster into blobs, and it turns out that this approach allows for a very precise determination of the backbone fractal dimension, much less affected by finite-size corrections than the more traditional approach of studying the backbone directly [50]. We provide estimates of this important dimension, which is one of the few exponents of the random-cluster model in two dimensions for which no exact expression is known, for a wide range of values in 0 < q < 4, many of which have been studied here for the first time.
An even more detailed understanding of the role of bridges could be expected from a generalization of the random-cluster model giving different weights to bridge-and non-bridge bonds. This problem of a generalized bridge percolation [19] is the subject of a forthcoming study. The O(n) loop model is another general class of frequently studied systems [51]. As these are graph polynomials as well (namely over a restricted set of spanning subgraphs, i.e., all loop configurations or Eulerian subgraphs [52]), it is therefore tempting to investigate in how far a similar Russo-Margulis approach can be applied there or, more generally, to Eulerian subgraph models [52], possibly yielding novel identities.

Acknowledgement
We would like to thank Timothy M. Garoni and Youjin Deng for valuable discussions. E.M.E. is grateful to the members of the School of Mathematical Sciences at Monash University for their warm hospitality during his stay in summer 2014, where part of this work was done. The authors acknowledge funding from the EC FP7 Programme (PIRSES-GA-2013-612707).

A. Influence of an edge
In order to prove (2.8), we first decompose the sum over subgraphs A ⊆ E into two contributions, corresponding to configurations for which e ∈ A and the complementary set of configurations with e / ∈ A.

27
Then observe that for A such that e ∈ A we clearly have A e = A and similarly for A with e / ∈ A we have A e = A. Furthermore we can also relate the probabilities of A and the modified set A e or A e once we know whether e ∈ A. Using these rules, we arrive at We remark that in the last sum one can safely discard the condition e ∈ A, because for A ⊂ E with e / ∈ A one trivially has A e = A and hence X(A) − X(A e ) = 0.

B. Free energy derivative
In order to evaluate E p,G [q K ], note that K only changes on removing a bond e if it is a bridge, i.e., We can now apply (2.8) to obtain: