For bond simulation, pc=0.5
p > pc, so it is possible.

     S = ( sum _{s = 1 to infinity}  n_s s² ) / ( sum _{s = 1 to infinity}  n_s s )

clusters in Figure 4-7:

cluster class      s      ns

1                      1      1
2                      4      1
3                      6      1
4                     46     1
 

substitute above {s, ns}  pairs into the formula.
 

For a simple cubic lattice, its critical exponents are:
 beta=0.41, gamma=1.80

P = 0                                         for p < pc
P = (p-pc)^beta = (p-pc)^0.41   for p > pc

S =  |p-pc|^-gamma =  |p-pc|^-1.8
 
 
 

The critical exponents for this lattice are:
beta=5/36, gamma=43/18, nu=4/3

S =  |p-pc|^-gamma =  |p-pc|^-43/18
xi  =  | p - pc | ^-nu =  | p - pc | ^-4/3

Exponent relations:

beta =  tau - 2/sigma
gamma =  3 - tau/sigma

beta = 0.48 - 2/2.31 =
gamma = 3 - 0.48/2.31 =
 
 

(1-p)² (p d/dp)²  sum _{s = 1 to infinity}  ( p^s ) =
(1-p)² sum _{s = 1 to infinity}  ( s²  p^s )

(1-p)² (p d/dp)( p/(1-p)2) =
(1-p)2 p [ 1/(1-p)2 + 2P/(1-p)3 ] =
p(1+p)/(1-p)
 

 sum _{s = 1 to infinity}  n_s s  = p
 

S = ( sum _{s = 1 to infinity}  n_s s² ) / ( sum _{s = 1 to infinity}  n_s s )

   = p(1+p)/(1-p)/p = (1+p)/(1-p)
 

p=exp(lnp)

and

-lnp ~ 1-p for p -> 1

and pc -> 1 (hence: -lnp ~ p_c-p)

we have:

g(r) = p^r = (exp(lnp))^r = exp(rlnp) = exp[-r(p_c-p)] =

        = exp[-r/(1/(p_c-p))]

Therefore:

 xi  ~ 1/(p_c - p)