Critical fragmentation properties of random drilling: How many random holes need to be drilled to collapse a wooden cube?

5.5 MB·

by  K. J. Schrenk

Critical fragmentation properties of random drilling: How many random holes need to
be drilled to collapse a wooden cube?
K. J. Schrenk,1,2, M. R. Hil´ario,3,4, V. Sidoravicius,5,6,7, N. A. M. Arau´jo,8,
∗ † ‡ §
H. J. Herrmann,1,9, M. Thielmann,10, and A. Teixeira11,
¶ ∗∗ ††
1Computational Physics for Engineering Materials, IfB, ETH Zurich,
Wolfgang-Pauli-Strasse 27, CH-8093 Zurich, Switzerland
2Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, U.K.
3Departamento de Matema´tica, Universidade Federal de Minas Gerais,
Av. Antonio Carlos, 6627 - PO Box 702 - 30161-970, Belo Horizonte, MG, Brazil
6 4Section de Math´ematiques, Universit´e de Gen`eve, 2-4 Rue du Li`evre, 1211, Gen`eve, Switzerland
1 5Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012, US
0 6New York University – Shanghai, 1555 Century Avenue, Pudong New Area, Shanghai, 200122, China
2 7CEMADEN, Avenida Doutor Altino Bondensan,

n 500, Sa˜o Jos´e dos Campos, SP, 12247-016, Brazil
a 8Departamento de F´ısica, Faculdade de Ciˆencias, Universidade de Lisboa,
J 1749-016 Lisboa, Portugal, and Centro de F´ısica Te´orica e Computacional,

4 Universidade de Lisboa, 1749-016 Lisboa, Portugal
1 9Departamento de F´ısica, Universidade Federal do Ceara´, 60451-970, Fortaleza, CE, Brazil
10Bayerisches Geoinstitut, University of Bayreuth,

] Universita¨tsstraße 30, 95440 Bayreuth, Germany
h 11Instituto Nacional de Matema´tica Pura e Aplicada,
c
Est. Dona Castorina, 110, 22460-320, Rio de Janeiro, RJ, Brazil
e
m
A solid wooden cube fragments into pieces as we sequentially drill holes through it randomly.
- Thisseeminglystraightforwardobservationencompassesdeepandnontrivialgeometricalandprob-
at abilistic behavior that is discussed here. Combining numerical simulations and rigorous results, we
t find off-critical scale-free behavior and a continuous transition at a critical density of holes that
s
significantly differs from classical percolation.
.
t
a
m PACSnumbers: 02.50.Cw,89.75.Da,64.60.ah
-
d
The connectivity of a solid block of material strongly fraction 1 p square-cells on each face is randomly se-
n −
o depends on the density of defects. To systematically lected and all sites along the line perpendicular to that
c study this dependence one must find an experimental face are removed (second row of Fig. 1). Thirty years
[ way to create defects inside the solid. For example, in ago,Y.Kantor[3]numericallystudiedthismodelonlat-

1 2D one can simply punch holes in a sheet and measure tices of up to 106 sites and concluded that the critical
v the physical properties of the remaining material. But fragmentation properties of this model are in the same
4 in 3D, inducing localized defects is not simple. One con- universality class as random percolation [4–6]. Here, we
3
ventionalsolutionconsistsinperforatingthematerialby combine rigorous results and large-scale numerical sim-
5
drilling holes or laser ablation from the surface [1, 2]. ulations, considering lattices three orders of magnitude
3
0 In a table-top experiment, we start with a solid cube larger in size, to show that this is not the case. Re-
. of wood and plot on each face a square-lattice mesh of L moving entire rows at once induces strong long-range di-
1
by L cells. Initially, the cube has no holes. Sequentially, rectional correlations and the critical behavior departs
0
6 for each one of three perpendicular faces, we randomly fromrandompercolation. Alsoremarkably,whileinran-
1 choose one square-cell and drill a hole having a radius dom fragmentation power-law scaling is solely observed
: of 1/√2 cell lengths to the other side of the cube. We around the critical threshold, here we find it in an en-
v
i repeat this process iteratively until the entire structure tire off-critical region. These findings suggest that long-
X
collapses into small pieces and the bottom and top part range directional correlations lead to a rich spectrum of
r of the cube are no longer connected. The first row of critical phenomena which need to be understood. Pos-
a
Fig. 1(a)-(d) shows the result of the drilling process of a sible implications for other complex percolation models
real cube with edge length 6 cm (manufactured from 2 are discussed in the conclusions.
cm thick plates of medium-density fiberboard (MDF)), Threshold. The average total number of drilled holes
whereholesweredrilledwithadiameterof1cm. Asthe is 3(1 p)L2 and the asymptotic probability that a site
−
drilling proceeds, pieces get disconnected and eventually in the bulk is not removed is p3. We first measure the
the entire structure collapses. thresholdp atwhichthecubecollapsesfordifferentlat-
c
Numerically, we start with a three-dimensional cubic tice sizes, up to L = 1024, using different estimators
lattice of L3 sites and fix three perpendicular faces. A of the transition point, as discussed in the Supplemental
2
(a)
2
10
101 Slope: 0.52±0.04
ν
/β
L
∞
P 100 128
256
(a) (b) (c) (d)
512
1024
-1
10
FIG. 1. (color online) Cube drilling. The upper row shows -1 0 1 2
10 10 10 10
photos of the experimental setup; the lower panels are the
1/ν
corresponding numerical results. The faces of the cube are (p−p )L
c
divided with a square-lattice mesh of linear size L=6, such
thateachfacecanbedrilledL2 =36times. Fromlefttoright,
(a)–(d),thenumberofdrilledholesperfaceare(a)0,(b)1, (b)
(c) 6, and (d) 8. From our numerical results for the position
-1
10
ofthetransitioninthethermodynamiclimit(seemaintext),
one estimates that around 13 holes need to be drilled for the
cube to disconnect.
-3
ν 10 Slope: -2.3±0.1
/γ
-
L
Material[7]. ExtrapolatingthedatatothelimitL→∞, M’2 10-5 128
gives pc = 0.6339 0.0005, consistent with the value 256
±
estimated by Kantor using Monte Carlo renormaliza- 512
1024
tiongrouptechniques[3](seeSupplementalMaterial[7]). -7
10
Thisthresholdislargerthanthetwo-dimensionalsquare- -1 0 1 2
10 10 10 10
latticepercolationthreshold(p )[8,9]andsmallerthan
2D
1/ν
thecubicrootoftheoneforthethree-dimensionalsimple |p−p |L
c
cubic lattice [10].
Static exponents. We consider the fraction P of
∞
sites in the largest cluster of connected sites (see Sup-
FIG. 2. (color online) (a) Double-logarithmic plot of the
plemental Material for more data of P (p) and of other rescaled order parameter P Lβ/ν as function of the scaling
observables [7]). P is the standard∞order parameter variable (p p )L1/ν for diff∞erent lattice sizes L. The linear
c
∞ −
in percolation identifying the transition from a discon- part has a slope of β = 0.52 0.04, consistent with β/ν be-
±
nected to a globally connected state. For the drilling ing the same as for three-dimensional percolation, but with
model, the situation will turn out to be more compli- a different exponent 1/ν = 0.92. The present value of β is
different from the two-dimensional one β = 5/36 0.139
cated. Figure 2(a) shows a double-logarithmic plot of ≈
[5, 11] and the one in three dimensions, β 0.417 [12].
the order parameter, rescaled by a power of the lattice ≈
(b) Double-logarithmic plot of the rescaled second moment
spize Pp∞c LL1β//νν.aBsafsuendctoinonfionfitteh-seizdeisstcaanlcinegtoanthaleytsrisan[5s]i,tiwone Mab2l(cid:48)eL−pγ/νp, wLi1th/νγ,/wνit=h12/.0ν4=520,.a9s2,fufonrctdiioffneroefntthleatstciacleinsgizevsarLi-.
|fin−d th|at the critical exponent of the order parameter is Thes|ol−idbc|lacklineisaguidetotheeyewithslope 2.3. The
−
β =0.52 0.04, and the inverse of the correlation length twosetsofdatapointscorrespondtothesub-critical(p<p )
c
±
exponent is 1/ν = 0.92 0.01 (see Supplemental Ma- and super-critical (p>pc) regions.
±
terial). We note that both β and ν are different from
the corresponding values for 2D and 3D classical per-
colation. However, somehow surprisingly, the exponent largestcluster. AsshowninFig.2(b),thefinite-sizescal-
ratio β/ν = 0.50 0.04 is within error bars the same inganalysisgivesthatthesusceptibilitycriticalexponent
±
as for 3D percolation. Thus, while the fractal dimension is γ = 2.3 0.1. We note that our results for the static
±
of the largest cluster (given by d = d β/ν) is consis- critical exponents are within error bars consistent with
f
−
tent with the one for 3D percolation, the larger value of the scaling relation 2β+γ =dν (for d=3).
β (compared to 3D percolation) implies that the transi- Dynamical exponents. The transport properties of the
tion from the connected to the disconnected state is less largestclusteratthecriticalthreshold,p ,areintimately
c
abrupt (see also Supplemental Material). We consider related to dynamical critical exponents and they can be
next the behavior of the second moment of the cluster measured by quantifying three sets of sites in the largest
size distribution M , excluding the contribution of the cluster [13]. First, we consider the so-called red sites.
2(cid:48)
3
A site is considered a red site if its removal would lead 106
to the collapse of the largest cluster [14]. The red sites Red sites 2.12±0.08
form a fractal set of fractal dimension d =0.92 0.07 5 Shortest path
RS ± s 10 Backbone
(see Fig. 3), which is compatible with the inverse of the e
t
correlation length exponent ν that we obtained from the si 4 1.30±0.05
f 10
finite-size scaling analysis in Fig. 2, d = 1/ν, as pre- o
RS
dicted by Coniglio for classical percolation [15]. How- er 3
b 10
ever the value of dRS = 1/ν for the drilling transition m 0.92±0.07
is very different from the classical 3D percolation re- u
sult 1/ν = 1.1437 0.0006 [12]. Figure 3 also shows N 102
± p=0.6339
that the shortest path connecting the top and bottom
1
sides of the largest cluster is a fractal of fractal dimen- 10
1 2 3
sion d = 1.30 0.05. Finally, the backbone of the 10 10 10
SP
±
largest cluster between its bottom and top ends is de-
Lattice size L
fined as the set of sites that would carry current if a
potential difference is applied between the cluster ends
FIG. 3. (color online) Number of sites in the backbone
(also known as bi-connected component). The backbone
of the spanning cluster, length of its shortest path (chem-
fractal dimension is determined as d = 2.12 0.08,
BB ± ical distance), and red sites in the backbone, measured at
which is larger than in classical 3D percolation, where
p = p = 0.6339, as function of the lattice size L. Con-
c
dBB =1.875 0.003 [13, 16, 17]. Qualitatively, an in- sidering the local slopes of the data, we obtain the follow-
±
crease in the backbone fractal dimension is compatible ingfractaldimensions: d =2.12 0.08,d =1.30 0.05,
BB SP
± ±
with a simultaneous decrease in the shortest path frac- andd =0.92 0.07. Thefractaldimensionoftheredsites
RS
±
tal dimension, since both correspond to a more compact is, within error bars, compatible with the value 1/ν 0.915,
≈
found from the finite size scaling behavior of the p estima-
backbone, similar to what is observed in long-range cor- c
tors (see Fig. S1), and the relation d = 1/ν [15]. Here,
related percolation [18, 19]. Thus, although the fractal RS
d is larger than in classical three-dimensional percolation,
dimension of the largest cluster is similar in both classi- BB
where d =1.875 0.003 [13, 16, 17], while d is smaller
BB SP
cal percolation and drilling, the internal structure of the thantheclassicalva±lued =1.3756 0.0006[20]. Thesolid
SP
largestclusterissignificantlydifferent. Thisimpliesthat linesareguidestotheeye. Toextract±thefractaldimensions,
transportandmechanicalpropertiesofthelargestcluster weanalyzedthelocalslopesasproposedinRef.[21]. Results
follow a different scaling. are averages over at least 5 103 samples.
×
Cluster shape. Given the highly directional nature of
the drilling process, we analyze the symmetry of the dif-
the correlation length for two-dimensional percolation
ferent clusters. In particular, we consider them as rigid
with parameter p. The event (respectively ) in-
bodies, consisting of occupied sites at fixed relative po- x y
R R
dicates the existence of a path crossing R (respectively
sitions, and look at the eigenvalues and eigenvectors of x
R ) from bottom to top that has not been drilled in the
their inertia tensors [22, 23]. The numerical results show y
y (respectively x) direction. By our choice of c , and
that,whencomparedtoclassicalpercolationclusters,the 1 x
R
havepositiveprobability(uniformlyoverk),i.e.there
drillingtransitionclustersaremoreanisotropic,theirori- y
R
exists a δ >0 such that P( ) δ. In addition, if (A)
entations being mainly aligned along the direction of the x
R ≥ S
cube edges (see Supplemental Material for quantitative occurs, then there exists a cluster spanning A [0,ec1k]
×
frombottomtotopandwhoseprojectionintothe(x,y)-
details [7]).
planedoesnotextendbeyondA. Thustheprobabilityof
We now give a rigorous argument for the existence of
asymmetric clusters in drilling percolation. Fix some finding a cluster of radius k and height ec1k is bounded
from below by the probability that there exists a square
p (p ,p ), where p 0.5927 < p is the criti-
2D c 2D c
∈ ≈ A along the diagonal x = y for which (A)
cal threshold for 2D site percolation. Consider a lattice X Y
S ∩R ∩R
size L L L with L = L = L and L = eL occurs, which is greater than
X Y Z X Y Z
× ×
and take a square domain A in its base with side length
k = c log(L) where c is a positive constant that 1 (1 δ2Lc0log[p(1−p)])L(c0logL)−1/2
0 0 − −
is sma(cid:112)ller than −[log(p(1−p))]−1. Say that the event ≥1−exp(−c−01/2δ2L1+c0log[p(1−p)]log(L)−1/2) ,
(A) occurred if, along the z-direction, no point inside
S
A is drilled but all points on its boundary are. For whichconvergestounityasLincreases. Thisshowsthat
large L, this event happens with probability at least one expects to have clusters extremely aligned along the
Lc0log(p(1−p)). ConsideralsotworectanglesRxandRy in z axis, asnumericallyobserved(seeforexampleFig.S15
the (x,z) and (y,z)-planes, respectively, aligned with A. of the Supplemental Material [7]). In fact, the same ar-
Theserectangleshavebaselengthkandheightexp c k , gumentcanbestraightforwardlyextendedtoexplainthe
1
{ }
where c is an arbitrary positive constant smaller than alignment along the x and y directions, as also observed.
1
4
100 not been drilled in the y-direction (x-direction). As dis-
cussedindetailintheSupplementalMaterial[7],forwell
3)
0 =0.610-2 chosen values of α and n the probability of the event C
3) 10 ∏p( LL=X=186 is bounded from below by a constant not depending on
=0.6 10-1 10-4 LXX=32 LthXei.rFjuorinthteorcmcourrereBncaendimCplaireesitnhdeepexenisdteenntceevoefnats.clSuisntceer
p 1 2 3 4
( Aspect ratio r includingsitesfromthebottomandthetopofthelattice,
∏
-2
y 10 we conclude that
abilit 10-3 r=1 Π(p)≥P[B∩C]=P[B]P[C]
Prob 10-4 rr==24 ≥exp{−c2αnlog(LX)}c3 ≥LX−θ ,
101 102 wherec2,c3,andθ arepositiveconstantsthatdependon
p.
Lattice size L
Z The above argument also shows the existence of
anisotropic clusters, sharpening the numerical results
FpI=G.0.46.3<(coplco,raosnfluinnec)tiMonaionfptlhoet:laStptiacnensiinzgepLrZob,afobrilidtiyffΠer,enatt spirmesileanrtleydtobewfohraet. hFaoprppen<s pfo2rDu,nocnoerrhealastΠed(pr)an∼doem−c4pLeXr-,
aspect ratios r. Solid black lines are guides to the eye with
colation, where c also depends on p. This is due to the
slopes, 0.26, 0.58, and 1.22, for r=1, 2, and 4. The 4
inset sh−ows the−same prob−ability as function of the aspect fact that the projection of a path spanning the lattice
ratio r, for different fixed lattice sizes L . Results are based into at least one of the coordinate planes is a path that
X
on at least 107 samples. spans the corresponding face, which has exponentially
small probability in L , due to the classical exponential
X
decay of connectivity in the subcritical phase [26, 27].
Spanning probability. To understand the properties of Conclusion. We find unexpected critical behavior
the drilling transition in terms of global connectivity, we when sequentially drilling holes through a solid cube un-
consider the spanning probability Π(p), defined as the til it is completely fragmented. At the critical density
probability to have at least one cluster including sites of drilled holes, a continuous transition is observed in a
from the top and bottom of the lattice, at a given value differentuniversalityclassthantheoneofrandomperco-
of the control parameter p. Figure 4 shows the spanning lation. We also numerically observe off-critical scale-free
probabilitybelowthethresholdΠ(p=0.63),fordifferent behaviorthatwecanjustifyforawiderangeofdensities
lattice aspect ratios. The lattice size is L L L ofholesusingrigorousarguments. Thismodelisarepre-
X Y Z
× ×
with L =L and L =rL and the spanning prob- sentativeofmorecomplexpercolationmodelswheresites
X Y Z X
ability is measured in the z-direction. At the drilling areremovedinastronglycorrelatedmanner[28–30]. Ex-
transition, p = p , the spanning probability approaches amplesaremodelswherethesetofremovedsitesisgiven
c
a constant for large lattice sizes (see Supplemental Ma- by randomized trajectories, such as the so called Pac-
terial, Fig. S9), similar to what is observed for classical man and interlacement percolation models proposed to
percolation [24, 25]. By contrast, for values of p between study the relaxation at the glass transition [31], enzyme
p and p , the numerical results suggest a power-law gel degradation [32] and corrosion [33, 34], as well as
2D
∗
decay of the spanning probability with L , where the percolation models for distributed computation [35, 36].
Z
exponent increases with the aspect ratio r. For fixed Other examples are percolation models explicitly intro-
L , it decays exponentially with r (see inset of Fig. 4). duce strong directional correlations as in the removal of
X
It is possible to establish rigorously the off-critical cylinders [37] and different variants of the four-vertex
power-lawdecayofΠ(p),modifyingtheargumentforthe model[38]. Itwouldbeinterestingtoexploreuptowhich
existenceofanisotropicclusterspresentedabove. Specif- degree these models are in the same universality class or
ically,wecanshowthatΠ(p)≥L−Xθ,whereθ =θ(p,r)> share common features.
0, for any fixed r > 0 and p (p ,p ). For that, let B While the fractal dimension of the largest fragment is
2D c
∈
bethediagonalband (x,y); x L /2 αnlog(L /n), consistent with the one of random percolation, all the
X X
{ | − |≤
x y <2n in the center of the (x,y)-face of the cube, other critical exponents are different. This has practical
| − | }
whereαandnareconstantssettingthelengthandwidth, implicationsastheconnectivityandtransportproperties
respectively. Let us say that the event occurred if B do change considerably close to the threshold of connec-
B
is free of holes in the (x,y) plane i.e. if no sites in B tivity. For example, we find the exponent of the order
are drilled in the z-direction. Also say that the event parameter to be substantially larger than for usual per-
C
occurred if there exists a path σ starting at height z =0 colation which implies that the drilling transition is less
and finishing at height z =L whose projection into the abrupt. Since sites are removed along a line, it is nec-
Z
(x,y)-plane is contained in B and whose projection into essary to remove more sites to produce the same effect
the (x,z)-plane ((y,z)-plane) consists of sites that have in the largest fragment. We also find that, compared to
usual percolation, the fractal dimension of the backbone [18] S. Prakash, S. Havlin, M. Schwartz, and H. E. Stanley,
is larger and the one of the shortest path is smaller, cor- Phys. Rev. A 46, R1724 (1992).
responding to a more compact backbone and therefore [19] K. J. Schrenk, N. Pos´e, J. J. Kranz, L. V. M. van
Kessenich,N.A.M.Arau´jo, andH.J.Herrmann,Phys.
enhanced conductivity properties.
Rev. E 88, 052102 (2013).
We acknowledge financial support from the ETH Risk
[20] Z. Zhou, J. Yang, Y. Deng, and R. M. Ziff, Phys. Rev.
Center, the Brazilian institute INCT-SC, and grant
E 86, 061101 (2012).
number FP7-319968 of the European Research Coun- [21] P. Grassberger, Physica A 262, 251 (1999).
cil. KJS acknowledges support by the Swiss Na- [22] J. Rudnick and G. Gaspari, Science 237, 384 (1987).
tional Science Foundation under Grants nos. P2EZP2- [23] M.L.MansfieldandJ.F.Douglas,J.Chem.Phys.139,
152188 and P300P2-161078. MH was supported by 044901 (2013).
[24] J. L. Cardy, J. Phys. A 25, L201 (1992).
CNPq grant 248718/2013-4 and by ERC AG “COM-
[25] S. Smirnov, C. R. Acad. Sci. Paris I 333, 239 (2001).
PASP”.NAacknowledgesfinancialsupportfromthePor-
[26] M. V. Menshikov, Dokl. Akad. Nauk SSSR 288, 1308
tuguese Foundation for Science and Technology (FCT)
(1986).
under Contracts nos. EXCL/FIS-NAN/0083/2012, [27] M.AizenmanandD.J.Barsky,Comm.Math.Phys.108,
UID/FIS/00618/2013, and IF/00255/2013. AT is 489 (1987).
grateful for the financial support from CNPq, grants [28] N. A. M. Arau´jo, P. Grassberger, B. Kahng, K. J.
306348/2012-8 and 478577/2012-5. Part of the rigorous Schrenk, and R. M. Ziff, Eur. Phys. J. Spec. Top. 223,
2307 (2014).
arguments presented in this work were based on math-
[29] M. Hila´rio and V. Sidoravicius, arXiv:1509.06204.
ematical results developed in the PhD Thesis of one of
[30] X.-W. Liu, Y. Deng, and J. L. Jacobsen, Phys. Rev. E
the authors, MRH [39].
92, 010103(R) (2015).
[31] R. Pastore, M. P. Ciamarra, and A. Coniglio, Fractals
21, 1350021 (2013).
[32] T. Abete, A. de Candia, D. Lairez, and A. Coniglio,
Phys. Rev. Lett. 93, 228301 (2004).
∗ [e-mail] [33] A.-S. Sznitman, Ann. Probab. 37, 1715 (2009).
† [e-mail] [34] A.-S. Sznitman, Ann. of Math. 171, 2039 (2010).
‡ [e-mail] [35] D. Coppersmith, P. Tetali, and P. Winkler, SIAM J.
§ [e-mail] Discrete Math. 6, 363 (1993).
¶ [e-mail] [36] P. Winkler, Random Struct. Alg. 16, 58 (2000).
∗∗ [e-mail] [37] J. Tykesson and D. Windisch, Probab. Theory Relat.
†† [e-mail] Fields 154, 165 (2012).
[1] P.SimonandJ.Ihlemann,Appl.Phys.A63,505(1996). [38] G. Pete, Ann. Probab. 36, 1711 (2008).
[2] S. Nolte, C. Momma, H. Jacobs, A. Tu¨nnermann, B. N. [39] M.R.Hila´rio,CoordinatepercolationonZ3,Ph.D.thesis,
Chichkov, B. Wellegehausen, and H. Welling, J. Opt. IMPA, Rio de Janeiro (2011).
Soc. Am. B 14, 2716 (1997).
[3] Y. Kantor, Phys. Rev. B 33, 3522 (1986).
[4] D. Stauffer, Phys. Rep. 54, 1 (1979).
[5] D. Stauffer and A. Aharony, Introduction to Percolation
Theory, 2nd ed. (Taylor and Francis, London, 1994).
[6] G. R. Grimmett, Percolation, Grundlehren der mathe-
matischen Wissenschaften, Vol. 321 (Springer, 1999).
[7] See Supplemental Material at.
[8] R. M. Ziff, Phys. Procedia 15, 106 (2011).
[9] J. L. Jacobsen, J. Phys. A 47, 135001 (2014).
[10] J. Wang, Z. Zhou, W. Zhang, T. M. Garoni, and
Y. Deng, Phys. Rev. E 87, 052107 (2013).
[11] S. Smirnov and W. Werner, Math. Res. Lett. 8, 729
(2001).
[12] Y. Deng and H. W. J. Blo¨te, Phys. Rev. E 72, 016126
(2005).
[13] H.J.Herrmann,D.C.Hong, andH.E.Stanley,J.Phys.
A 17, L261 (1984).
[14] Numerically, we consider those sites that if removed will
disconnect all possible paths between a unit site in the
topandoneinthebottomfacesofthecube.Thesitesare
selectedsuchthattheirEuclideandistanceismaximized.
[15] A. Coniglio, Phys. Rev. Lett. 62, 3054 (1989).
[16] M. D. Rintoul and H. Nakanishi, J. Phys. A 27, 5445
(1994).
[17] Y. Deng and H. W. J. Blo¨te, Phys. Rev. E 70, 046106
(2004).
Supplemental Material: Critical fragmentation properties of random drilling: How
many random holes need to be drilled to collapse a wooden cube?
K. J. Schrenk, M. R. Hil´ario, V. Sidoravicius, N. A. M. Arau´jo, H. J. Herrmann, M. Thielmann, and A. Teixeira
I. LOCATING THE DRILLING TRANSITION
0.66
Todeterminethepositionofthedrillingtransition,we
0.65
considerinFig.S1thefinite-sizebehavioroffivedifferent
estimators: )L 0.76 0.64
(
• χti∞onroeffetrhsetloartgheestmcaluxsimteurmsizien. tThheisstiasntdhaerdpodseitviioan- or pc 00..7724 χMΠ∞2’i 0.603.0000 0.0125 0.0250
t
of the peaks in the curves in Fig. S6. a Π
m 0.70 J
i
• Mme20ntis,etxhcelupdoisnigtitohneocfontthreibpuetaioknoofftthheelsaercgoensdtcmluos-- ld est 00..6668 KKK123
o
ter, as plotted in Fig. S7. h
s 0.64
e
r
• Πi(L) is the value of p where the curves of Π(L,p) Th 0.62
and Π(L/2,p), i.e. the spanning probabilities seen 0.00 0.05 0.10 0.15
in Fig. S2, intersect [1].
−0.915
Inverse lattice size L
Π is the average value of p where the first cluster
•
spanning the lattice in z-direction appears [2].
FIG.S1. Estimatesforp asafunctionoftheinverselattice
c
• J is the point where the largest change in size of sizeL−0.915. Extrapolatingtothethermodynamiclimityields
the largest cluster occurs, for every given sample, p = 0.6339 0.0005. The data points are averages over at
c
±
averaged over all samples [3]. least 7700 samples (details in the text). Solid black lines are
guides to the eye.
K1, K2, and K3 are the values obtained by Y.
•
Kantor using Monte Carlo renormalization tech-
niques [4]. p = p = p = p . A system is considered to be con-
x y z
nected if there exists a path of occupied nearest neigh-
bors connecting two opposing faces of the cube, say in z
direction. LetusdenotebyΠ(p)theprobabilitythatthe
II. DETAILED MODEL DESCRIPTION system is connected. Then, in the thermodynamic limit,
Π(p) is a step function with Π(1)=1 and Π(0)=0 and
Consider a simple cubic lattice of linear size L with withatransitionatp=pc [15]. InFig.S2weseethecon-
free boundary conditions. To explore different fractions nectionprobabilityasafunctionofthecontrolparameter
of drilled holes, we use the following setup. In the initial p for different linear system sizes L. It can be observed
configuration, there are no holes, such that all N = L3 that the curves intersect around pc = 0.6339 0.0005
±
unit cells (sites) are present (occupied). (see also Fig.S1). Given its error bars, this value of p is
higher than the site percolation threshold for the square
Weconsiderthefollowingprocess. Oneachfaceofthe
lattice p =0.59274621 00000013 [1, 16, 17] and lower
cube i = x,y,z a site x is selected uniformly at ran- c
0,i ±
than the cubic root of the site percolation threshold of
dom among the unselected sites, the corresponding hole
the simple-cubic lattice: [p (s.c.)]1/3 0.67795 [6].
is drilled, and all occupied sites in the direction perpen- c
≈
dicularto theface andoriginating fromx are removed
0,i
(labeled as unoccupied). Therefore, in this step all oc-
cupied sites with coordinates x = x + ke are re- III. OBSERVABLES AS FUNCTION OF p
k,i 0,i i
moved, where e is the unit vector perpendicular to the
i
cube face and k =0,...,L 1, such that at most L sites The connection probability Π(p) is measured as the
−
can be drilled in one step. This step is executed 3L2 ratio of samples that are connected at parameter value
times until all holes have been drilled and consequently p and the total number of samples. For example, if all
all sites in the system are gone. If we denote by (1 p ) holes have been drilled, we have p = 0 and Π(0) = 0,
i
−
the fraction of drilled holes in face i, we can character- and likewise if no hole has been drilled, Π(1) = 1 (see
ize the state of the process with the control parameter Fig. S2).
2
TABLEI. Criticalexponents. Thefirstcolumnshowsthesymboloftheexponent. Thesecondcolumnshowstheresultsforthe
drillingtransition). Forcomparison,literaturevaluesfortwo-andthree-dimensionalclassicalpercolationaregivenincolumns
threeandfour[1,5]. Here,β/ν istheexponentrelatedtotheorderparameter,seeEq.(S1). γ/ν istheexponentrelatedtothe
second moment [see Eq. (S7)] of the cluster size distribution and to the susceptibility [see Eq. (S3)]. 1/ν is the inverse of the
correlation length critical exponent, as determined from the finite size scaling of threshold estimators (see Fig. S1) [1, 6]. β is
the order parameter critical exponent, determined from the behavior P (p p )β [see Fig. 2(a)] [1]. γ is the susceptibility
c
critical exponent, determined from the behavior M20 ∼|p−pc|−γ (see∞Fig∼. 2)−[1, 6]. dRS is the fractal dimension of the red
sites in the backbone of the largest cluster [6–9]. d is the fractal dimension of the shortest path in the largest cluster [10].
SP
d is the fractal dimension of the backbone of the largest cluster [11–14]. τ is the exponent of the cluster size distribution
BB
(see Fig. S11).
Exp. Drilling Classic 2D Classic 3D
β/ν 0.50 0.04 5/48 0.1042 0.4774 0.0001
± ≈ ±
γ/ν 2.04 0.05 43/24 1.7917 2.0452 0.0002
± ≈ ±
1/ν 0.915 0.010 3/4=0.75 1.1450 0.0007
± ±
β 0.52 0.04 5/36 0.1389 0.4169 0.0004
± ≈ ±
γ 2.3 0.1 43/18 2.3889 1.7862 0.0013
± ≈ ±
d 0.92 0.05 3/4=0.75 1.1450 0.0007
RS
± ±
d 1.30 0.05 1.13077 0.00002 1.3756 0.0006
SP
± ± ±
d 2.14 0.08 1.6432 0.0008 1.875 0.003
BB
± ± ±
τ 2.1 0.1 187/91 2.0549 2.18925 0.00005
± ≈ ±
1.0 0.06
8
0.8
16
32
0.04
0.6 64
128
Π 256 0.3 C
0.4 512
1024 0.02
0.2
0.2
C=0.0252(3)
0.1
0.633 0.634 0.635
0.0 0.00
0.5 0.6 0.7 0.8 0.0 0.2 0.4 0.6 0.8 1.0
p p
FIG. S2. Connection probability Π as a function of the FIG. S3. Number of clusters per site C as a function of the
control parameter p. The systems consist of N = L3 sites controlparameterp. Latticesizesandnumbersofsamplesare
where the largest considered linear system size is L = 1024 thesameasinFig.S2. Thecurvesseemtoapproachalimiting
and the smallest one is L = 8. Results have been averaged curve with increasing lattice size. The data extrapolates to
over at least 107 samples for the smallest lattice size and 104 C(0.6339)=0.0252 0.0002inthethermodynamiclimit. The
±
samples for the largest one. Π is measured as the fraction of solid black straight lines are guides to the eye intersecting in
connectedsamples. Thecurvesforlargelatticescrossaround C(0.6339) = 0.0252. The considered lattice sizes and colors
p =0.6339 0.0005andΠ(p )=0.21 0.03. Forclarity,the of the curves are the same as in Fig. S2.
c c
± ±
inset only shows the data for the four largest lattice sizes.
ofsitesN =L3. Forcriticalpercolation,thisisknownto
In the simple-cubic lattice, every site in the bulk has
be a lattice-dependent constant [21–24], the same seems
sixnearestneighborsites. Iftwooccupiedsitesarenear-
to hold for the drilling transition, Fig. S3.
est neighbors they are said to be connected. A set of
connected occupied sites is called a cluster. Two dis- ThelargestclustersizepersiteP isdefinedasthera-
tinctclustersarethereforeseparatedbyunoccupiedsites. tioofthemaximumclustersizeam∞ongallclustersinthe
The number of sites forming the cluster is called its size. system and the total number of sites N = L3. For clas-
A single occupied site which has six unoccupied nearest sical percolation, P acts as order parameter [1] which
neighbor sites is considered to be a cluster of unit size. is zero in the disor∞dered state and nonzero in the or-
The number of clusters per site C is defined as the ratio deredone. Thedataforthedrillingtransitionisshownin
of the number of distinct clusters and the total number Fig. S4. For percolation close to the percolation thresh-
3
1.0
0.08
2
0.8 ν 0.06
0.6 /βPL∞1 0.06 3-/γνL0.03
∞ 0 2∞
P -1 0 1 2 ∞ 0.04 χ
0.4 (p-p )L1/ν χ 0.00
c -1 1 3
0.2 0.02 (p-pc)L1/ν
0.0
0.00
0.4 0.6 0.8 1.0
0.0 0.2 0.4 0.6 0.8 1.0
p p
FIG.S4. LargestclustersizepersiteP asafunctionofthe
FIG. S6. Standard deviation of the largest cluster size per
∞
controlparameterp. Intheinsets,oneseestherescaleddata
site χ as a function of the control parameter p. The inset
for the four largest lattice sizes around pc. We obtain a data shows∞the data for the four largest lattice sizes rescaled with
collapse for p = 0.6339 and 1/ν = 0.915 and β/ν = 0.4774,
c p = 0.6339, 1/ν = 0.915, and γ/ν = 2.0452. The ratio
c
whereβ/νissettothevalueforthree-dimensionalpercolation
γ/ν is the one of three-dimensional percolation [6]. 1/ν =
[6].
0.915 0.010 is consistent with the scaling behavior of P
± ∞
(see Fig. S4) and the estimators for p (see Fig. S1), but it
c
differs from the value for percolation, 1/ν =1.1450 0.0007
∞∞00..55 ±
PP [6, 18–20].
r r 64
ee
stst 00..44 128
uu 256 at lower fractions of removed cells in the case of classical
clcl Drilling
percolation, see Fig. S5.
estest 00..33 The standard deviation of the largest cluster size per
gg
rr site χ is defined as the square root of the variance of
lala 00..22 Classic P : ∞

in in 64 1−pc 1−(pc)3 ∞
ionion 00..11 128 χ∞ = hP∞2 iS −hP∞i2S, (S2)
ctct 256 q
rara 00..00 where h·iS indicates averaging over independent realiza-
FF 00..55 00..66 00..77 00..88 tions (see Fig. S6). For percolation close to the percola-
FFrraaccttiioonn ddrriilllleedd ssiitteess tion threshold, one expects the following scaling behav-
ior:
FIG. S5. Fraction of sites in the largest cluster P as func-
tionofthefractionofsitesthathavebeendrilled,f∞ordrilling χ2∞(p,L)=L−d+γ/νFχ[(p−pc)L1/ν], (S3)
percolation (upper curves), and classical percolation (lower
where d is the spatial dimension of the lattice, γ is the
curves). For both models, data is shown for three different
critical exponent related to the susceptibility, and F is
latticesizesL. Thearrowsindicate1 p ,withp theclassi- χ
c c
calsitepercolationthreshold,and1 −(p )3,respectively. For a scaling function.
c
− The second moment of the cluster size distribution is
both models, this corresponds to the fraction of drilled sites
at which the percolation transition occurs, in the thermody- defined as
namic limit. Results are averaged over 104 samples.
1
M = s2, (S4)
2 N k
k
X
old, one expects the following scaling behavior:
where the sum runs over all clusters in the system and
P (p,L)=L−β/νFP[(p pc)L1/ν], (S1) sk is the size (number of sites) of cluster k. Excluding
∞ − the contribution of the largest cluster size, we arrive at
where β is the critical exponent related to the order pa- the following definition:
rameter, ν is the critical exponent related to the correla-
M =M s2 /N, (S5)
tion length, pc is the percolation threshold, and FP is a 20 2− max
scaling function [1]. If we compare classical percolation
where s is the largest cluster size. In Fig. S7, we see
max
anddrillingwithrespecttoP againstthefractionofre-
moved cells, we observe that∞the cube gets disconnected M /N =M /N s2 /N2. (S6)
20 2 − max
4
0.005 10-1
0.04
0.004
-2
ν 10
N 0.003 γ-/M’L20.02 −0.50±0.04
M’/2 0.002 0.00-2 0 2 10-3 P∞
(p-p*)L1/ν M2’/N
0.001 10-4 −0.96±0.05
0.000
-5
0.0 0.2 0.4 0.6 0.8 1.0 10
1 2 3
10 10 10
p
L
FIG. S7. Second moment of the cluster size distribution
withoutthecontributionofthelargestclusterpersiteM20/N FIG. S8. Lattice size dependence of the largest cluster size
asafunctionofthecontrolparameterp. Theinsetshowsthe persiteP andthesecondmomentoftheclustersizedistri-
data for the largest four lattice sizes rescaled using the same butionwit∞houtthecontributionofthelargestclusterpersite
values of γ/ν, 1/ν, and pc as in Fig. S6. M20/N atp=pc =0.6339. Theasymptoticslopesofthelines
give β/ν = 0.50 0.04 and γ/ν = 2.04 0.05, in agreement
± ±
with the reported values for three-dimensional percolation,
Close to the percolation threshold, one expects the fol- β/ν =0.4774 0.0001 and γ/ν =2.0452 0.0002 [6].
lowing scaling behavior: ± ±
M (p,L)=Lγ/νF [(p p )L1/ν], (S7)
2 M2 − c
whereF isascalingfunction. Accordingtothescaling
M2
forms in Eq. (S1) and (S7), one expects to observe Asafurthertest,wealsoconsideredalternativebound-
aryconditions,obtainingresultsconsistentwiththeones
P (pc,L) L−β/ν (S8) obtained for free boundary conditions. In particular, we
∞ ∼
employed a single-cluster growth method, similar to the
and
Leath algorithm for classical percolation [32, 33]. This
M (p ,L)/N L d+γ/ν (S9) method works at a fixed occupation probability p. We
20 c ∼ − begin with a lattice of length L equal to a small multi-
at the threshold p = p . The corresponding measure- ple of three. The holes in this initial cube are present
c
ments are shown in Fig. S8. with probability 1 p and absent with probability p. If
−
thesiteinthecenterofthelatticeisoccupied(thisisthe
casewithprobabilityp3),aclusterisgrownstartingfrom
IV. ALGORITHMS AND BOUNDARY there, connecting to all occupied neighbors. In case the
CONDITIONS cluster touches a site on the boundary of the lattice, the
linearlatticesizeisincreasedbyafactorofthree, thein-
To simulate the drilling model with an algorithm of formation on the drilled holes is accordingly propagated
complexity linear in the number of sites N = L3, we to the new lattice and the cluster growth continues. At
adopt the following strategy. We consider initially a sys- a certain cutoff lattice size L, the cluster growth is ter-
tem with all sites occupied. Starting from this setup, minated. Thismeasurementgivesthenumberofclusters
theholesaredrilledaccordingtotheprocessesdescribed as function of the size of the cluster at the origin, as
above. Foreachsite, thefractionofdrilledholes(1 p ) shown in Fig. S10, at p = p . Because the probability
x c
−
at which it is removed is recorded. Once all holes have of a randomly picked site to be in a cluster of size s is
been drilled, we proceed in the opposite direction, and proportionaltos,thismeasuresasizedistributionwhich
analyze the following percolation problem. Suppose we behavesassp(s). Thecorrespondingpower-lawexponent
start from a lattice where all sites are unoccupied (p = is measured as τ 1 = 1.165 0.01, see Fig. S11. The
x
− ±
0). Now, p is incremented and the sites of the lattice influence of varying p on the cluster size distribution is
x
becomeoccupiedinthe inverseorder oftheonein which shown in Fig. S12. Finally we verified that the cutoff of
they have been drilled. This procedure allows to keep theclustersizedistribution,s ,scalesasexpectedfrom
cut
track of the properties of the clusters of occupied sites thefractaldimensionofthelargestcluster: Wemeasured
as function of the control parameter p [2, 16, 17]. Ran- the fraction f(L) of clusters grown from the origin that
x
dom numbers have been generated with the algorithms reach the boundary of a box of length L and found that
discussed in Refs. [29–31]. f(L) L 0.48 0.10 (not shown). One expects f to be
− ±
∼
5
(a) 105
9) 100 104 9
3 s 3 27
3 r 10 81
0.6 10-1 uste 102 243
p= cl 101
∏( 10-2 of 100
y er 10-1
t b
li -3 m 10-2
bi 10 r=1 u -3
a r=2 N 10
b
o -4 r=4 10-4
r 10
P 1 2 10-5
10 10 0 1 2 3 4 5 6 7 8
10 10 10 10 10 10 10 10 10
Lattice size L Size of cluster at origin
Z
(b)
0 FIG. S10. Number of clusters of size s versus size s of
10
) the cluster at the origin, for p=p and different values of
9 c
3 themaximumconsideredlatticesize. Thehistogramsareob-
63 -1 tained from 105 clusters.
. 10
0
= 6
p 10
∏( 10-2 4
s 10
obability 10-3 LLLLXXX===X=1366248 f cluster 110002
L =128 o
Pr 10-4 X r 10-2
e
1 2 3 4 mb 10-4
u
Aspect ratio r N 10-6 Data
Slope: −1.165
-8
FIG.S9. (a)showstheconnectionprobabilityΠasfunction 10
0 1 2 3 4 5 6 7
of the lattice length L , with p=0.6339, for different lattice 10 10 10 10 10 10 10 10
Z
aspect ratios r. Π saturates to a constant for large LZ. As Size of cluster at origin
can be seen in (b), for large aspect ratios r, the connection
probability decays exponentially in r, similarly to classical
percolation [25–28]. The lattice size is L L L with FIG. S11. Number of clusters of size s versus size s of the
X Y Z
LX =LY and LZ =rLX and the connectio×n pro×bability is clusterattheorigin,forp=pcandmaximumlatticesize729.
measured in z-direction. Results are based on at least 108 Thedataisobtainedbyaveragingover102 histogramsof105
samples. clusters. The solid blue line is a guide to the eye with slope
1.165.
−
V. CLUSTER SHAPE AND POWER-LAW
related to the cluster size distribution p(s) s τ by DECAY
−
∼
To study the geometry of the drilling transition clus-
ters, we considered the eigenvalues and eigenvectors of
f(L)= ∞ p(s)sds Ldf−d, (S10) the cluster inertia tensors [37–39], see Fig. S13 to S15.
∼ The results show that the clusters of the drilling transi-
Zscut(L)
tionaremoreanisotropicandalignedwiththecubeedges
as the ones of classical 3D percolation.
Wenowpresentthedetailsinvolvedintheproofofthe
assuming scut Ldf [34, 35] and the validity of the scal- polynomial decay for the spanning probability Π(p) as a
∼
ing relation τ = 1+d/d [1, 36], in analogy to classical function of the lattice length L (recall that L = L
f X Y X
percolation. InFig.S8,wemeasured d= 0.50 0.04, and L =rL with a fixed aspect ratio r). We fix a p
f Z X
− − ± ∈
in agreement with the behavior of f. (p ,p ). Weconsiderthediagonalbandcontainedinthe
2D c