<figure><img src="image_1.png" /></figure>

HAL
open science

# Locality approach to the bootstrap percolation paradox

Ivailo Hartarsky, Augusto Teixeira

► To cite this version:

Ivailo Hartarsky, Augusto Teixeira. Locality approach to the bootstrap percolation paradox. Physical Review Letters, 2025, 134 (11), pp.117102. 10.1103/PhysRevLett.134.117102. hal-05021998

HAL Id: hal-05021998

https://hal.science/hal-05021998v1

Submitted on 5 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Locality approach to the bootstrap percolation paradox

Ivailo Hartarsky*

Technische Universität Wien, Institut für Stochastik und Wirtschaftsmathematik,
Wiedner Hauptstraße 8-10, A-1040, Vienna, Austria

Augusto Teixeira†

IMPA, Estrada Dona Castorina, 110 - Rio de Janeiro, Brazil

(Dated: January 9, 2025)

We revisit the Bootstrap Percolation model, leveraging recent mathematical advances linking it with its local counterpart. This new perspective resolves, for the first time, historic discrepancies between Monte Carlo simulations and theoretical results: previously, those predictions disagreed even in the first-order asymptotics of the model. In contrast, our framework achieves excellent agreement between numerics and theory, which now match up to the third-order expansion, as the infection probability approaches zero. Our algorithm allows us to generate novel predictions for the model.

Monte Carlo simulations provide an accessible tool for probing physical systems and have been used extensively. Yet, the reliability of results based on such approaches is not guaranteed. Along with explosive percolation [1–3], bootstrap percolation (BP) is among the prime examples of how Monte Carlo simulations can yield misleading results. Indeed, for BP, going back to [4–6], this phenomenon is so common that it has been labelled the *BP paradox* [7] (also see [8–11]). BP is a paradigmatic statistical mechanics model of metastability with deep connections to the low temperature stochastic Ising model [12–15], kinetically constrained models of supercooled liquids [16–19] and the *k*-core model [20], as well as applications to areas as diverse as complexity theory [21], influence in social networks [22], sandpile formation [23] and many others. We direct the reader to the surveys [11, 24–27] for background on BP, simplistically schematised below.

In early years, physical systems for which BP it was relevant were identified. Not necessarily the ones suggested initially concerning magnetism, but rather in solid-state physics and rock geometry. Simulations were used to determine the large scale behaviour of the model. This gave rise to many predictions, attracting the application of heuristic arguments and theoretical techniques, as well as the attention of the probability community.

A much more striking aspect of history is the BP paradox. To summarise it, each prediction of asymptotic behaviour based on numerics was subsequently proved erroneous by rigorous results. The list is impressive: triviality of the phase transition [6, 28–30], critical scaling [9, 31–33], sharp threshold [9, 34–38], second order critical exponent [9–11, 39–42] and second order logarithmic corrections [43, 44] have all gone through this pattern over the years. Two decades ago [7, 39, 40], an important attempt at reconciling theory and numerics was made. Yet, this led to more discrepancies [11, 41] and ultimately only widened the gap between communities: “The striking conclusion is that not even the rigorous correction term can be captured reliably by the numerical

exact solution yet” [11]. Following our recent progress on the rigorous theory side [45] and, presently, the numerics side, we are now able to settle this recurrent controversy, explain its origin and, finally, show that it is possible to make accurate predictions for these models based on numerical results.

MODELS

We will discuss two closely related BP models on the square lattice $Z^2$. Start with a set $A_0 \subset Z^2$ of initially infected sites (also known as + spins, vacancies, etc. in other contexts) with density $p \in (0, 1)$, selected independently at random. At each discrete time step we define $A_t$ by adding to $A_{t-1}$ all sites $\mathbf{x} \in Z^2$ fulfilling the following condition.

* **MBP** - Modified two-neighbour BP: there is at least one infection in $A_{t-1}$ among both the horizontal neighbours $\{\mathbf{x} - \mathbf{e}_1, \mathbf{x} + \mathbf{e}_1\}$ and the vertical neighbours $\{\mathbf{x} - \mathbf{e}_2, \mathbf{x} + \mathbf{e}_2\}$ of $\mathbf{x}$ (e.g. $\mathbf{e}_1 \to \mathbf{e}_2$).

* **FBP** - Froböse BP: there is a four-cycle xyzw such that $\mathbf{y}, \mathbf{z}, \mathbf{w}$ are all infected in $A_{t-1}$ (e.g. $\mathbf{z} \to \mathbf{w}$).

We remark that the classical two-neighbour constraint can be treated like FBP [45], but is omitted for brevity.

Since, in infinite volume, BP exhibits a trivial phase transition with $p_c = 0$ [28], we focus on the low temperature regime $p \to 0$. As it is common in metastability, key quantities in BP can be expressed in terms of the frequency (density) $\rho$ of occurrence of certain critical droplets enabling nucleation [31]. In our models, a critical droplet is simply an axis-parallel square of side length $\Lambda = C \log(1/p)/p$ (with $C$ an appropriately large constant, for us $C = 2$ suffices) which is *internally filled*: the initial infections of $A_0$ inside the square are sufficient to infect it completely, see Fig. 1. Once $\rho$ is known, one can readily determine other quantities of interest like the
<figure><img src="image_3.png" /><figcaption>**FIG. 1:** Sample critical droplets conditioned to be internally filled. Colors represent infection times, with grey being the initial infections, while the bars to the left of each figure indicate the correspondence between colors and time of infection. The little windows contain a zoomed-in picture of the early stages of nucleation, with initial infections in lighter grey. Close inspection reveals that local FBP internally fills with the same initial condition, while local MBP does not. The dotted rectangles on the right show two moments where local MBP becomes confined.</figcaption></figure>

typical infection time of the origin in infinite volume, the critical initial infection density for a large box (or torus) to be internally filled with probability 1/2, the critical length of a box so that it becomes fully infected with probability 1/2 at given initial infection density, etc.

## BACKGROUND

It is useful to recall previous approaches for estimating the critical droplet density $\rho$. The most natural one consists in running a Monte Carlo simulation of a square of side $\Lambda$ with initial infection density $p$ and record the empirical proportion of samples infecting the entire square. The main drawback of this straightforward approach, which has been extensively used [6, 9, 29, 33–36], is that its running time is governed by $\rho^{-1}$ rather than $\Lambda$. In metastability settings like ours, this leads to an exponential complexity, since $\rho^{-1} \sim \exp(C'/p)$ [31], as we will see. This practically limits Monte Carlo simulations to moderate densities $p > 0.05$, see Fig. 2. Despite the significant computational power needed to obtain such results, extrapolating them to predict thermodynamic results is risky.

An approach to go beyond moderate densities was proposed in [7, 39, 40], based on a local version of BP. This simplified model [41] prohibiting polynucleation is defined as follows. Consider a growing sequence of nested rectangles $R_n$ as follows. The seed $R_0$ is a single initially infected site. At each step, $R_{n+1}$ is obtained from $R_n$ by increasing the height or width by 1, but we require that the additional line (row or column) contains at least one initial infection. For local MBP we also allow $R_{n+1}$ to be obtained from $R_n$ by adding both a row and a column to $R_n$, provided the corner at their intersection is initially

infected. The local density $\rho_\ell$ is the probability that, starting at $R_0 = \{0\}$, there exists such a sequence reaching a square of side $\Lambda$. The difference with the original BP is that we do not allow multiple non-trivial rectangles to first grow and then merge.

Local BP is much more tractable, because the probability of reaching a rectangle of given size can be computed recursively [40, 45]. In order to do this, we need to record some additional information on which sides of the current rectangle we have already attempted and failed to extend. Yet, there are only a finite number of such ‘frame states’ (seven for MBP and eight for FBP [45]), so these probabilities can be computed recursively. This yields a dynamic programming algorithm for computing $\rho_\ell$ with time and memory complexities of $\Lambda^2$ and $\Lambda$ respectively. As opposed to the exponential cost of Monte Carlo simulations, this local approach is polynomial in $1/p$. This allowed probing densities down to $p \approx 0.0035$ [39] for local MBP with minimal computational resources. Moreover, the output of this algorithm is deterministic and only subjected to numerical precision errors, which can be controlled to arbitrary precision.

Strikingly, local BP computations still led to incorrect predictions of the asymptotic scaling of $\rho_\ell$ for MBP as $p \to 0$ [41], even when informed of the rigorous result $\rho \sim \rho_\ell \sim \exp(-\pi^2/(3p))$ [37]. The disappointing results of this very promising method essentially extinguished the hope of numerical and rigorous results reaching agreement [11]. Consequently, the numerical side of the area has seen no activity in the last decade, since the last incorrect predictions of [43], based on a version of local BP, later refuted in [44]. During the same period, the rigorous theory has been flourishing with a number of very precise results for specific models [45–50] and the development of a unified universality theory for BP and the
closely related kinetically constrained models [51–58].

## RIGOROUS LOCALITY APPROACH

The present letter is largely inspired and complemented by our recent mathematical work [45]. It is our goal to elucidate the new perspective it offers along with its implications on the rigorous and numerical side of the BP problem. Let us start by stating one of our main results [45, Theorem 1.1, Lemma 3.5], the locality theorem. For FBP, the non-locality $log ρ - log ρ_{ℓ}$ is at most polylogarithmic in $1/p$ [59]. This establishes in a very strong sense that $ρ \sim ρ_{ℓ}$. We also provide compelling numerical support for this fact in Figs. 1a and 2a. While previously local FBP was introduced as a simplification of FBP, it is now clear that these two models are essentially equivalent. The locality theorem is important to both theory and numerics, because the local model is much more accessible. Owing to locality, we could prove [45, Theorem 1.2] that, for FBP,

$$ \log \rho \sim \log \rho_{\ell} \sim \frac{-\pi^2}{3p} + \frac{2\pi\sqrt{2+\sqrt{2}}}{\sqrt{p}}. \quad (1) $$

At present, we do not seek to discuss the proof of (1) (see [45]), but rather the implications of the locality theorem to numerics, keeping (1) as a verification reference.

Since the discrepancy between the local and non-local infection times is much smaller than the precision currently feasible for $log ρ$, the numerical approach to local BP in fact yields precise results also for the non-local model. Thus, it is all the more important to understand the reasons for its previous failure to predict the critical exponent $1/2$ of the second term in (1), even given the first order term. In order to do this, we have performed a contemporary implementation of the dynamical programming algorithm. While this is not a difficult task, some care is needed to obtain an efficient implementation.

There are two key novel features in our simulations. First, we take advantage of the parallel nature of the computation through the use of a GPU. Also, we do not store the desired probabilities directly and take care of error propagation when operating over them. This is necessary, since the probabilities $ρ_{ℓ}$ that we attain are as small as $10^{-2} \cdot 10^5$ and some off-diagonal values can be as small as $p^{-\Lambda} \approx 10^{-2700000}$. So a trade-off is needed between numerical precision and efficiency. See more details in the supplemental material. We alert the reader that, even though such extraordinarily small numbers have little meaning in physical context, the results they will be used to obtain do have tangible implications in regimes of physical interest.

Running the algorithm discussed above and available at the repository https://github.com/

augustoteixeira/bootstrap (see folder taichi), we are able to reach $p \approx 8 \cdot 10^{-6}$. The resulting values of $p \log ρ_{ℓ}$ are presented in Fig. 2 (see the supplemental material for the numerical values). For comparison, we have also implemented the straightforward Monte Carlo strategy (see supplemental material), displaying one sample internally filled critical droplet in Fig. 1 and data for $p^2 ρ$ in Fig. 2. A detailed explanation for the corrective factor $p^2$ is provided in the supplemental material, but it is necessary for roughly compensating the volume $\Lambda^2$ of the droplet and its contribution is only relevant for the large and moderate values of $p$ accessible to Monte Carlo simulations.

In the words of [40], “If we are to interpret quantitative data reliably, there is now a pressing need for advancement of the whole theoretical agenda.” Indeed, efficient data analysis is a crucial and novel point of our work. Namely, we develop a reliable scheme for extracting extremely precise asymptotics and determining their correct functional form. The details are provided in the supplemental material, since we expect such treatment to be of use well beyond the scope of BP. But let us present our results obtained by following these guiding principles: allow general functional forms rather than artificially imposing power laws; use a wide range of data to roughly estimate the correct functional form of the subsequent error term before estimating the main one; use only the latest data for precise estimates, given the functional form.

As an illustration, consider estimating the limit $π^2/3 \approx 3.28987$ from the data in Fig. 2a. What one should not do is to take the value of the last data point. In our case, this gives 3.2586 for FBP (and even 3.2122 for MBP). Instead, we use the last few points to infer the functional form of the leading and second order terms. Once this is done, performing a 3-parameter fit of the form $p \log ρ_{ℓ}^{-1} \sim a + bp^{-c}$, using the last four data points, gives $a \approx 3.290(3)$: four correct digits instead of one or two. Proceeding similarly, we are able to infer the functional form in (1) and the constants appearing, as well as the third order term, which is not known rigorously. For the second order constant $2\pi\sqrt{2+\sqrt{2}}$, we also obtain three correct digits. We estimate the third order exponent for FBP to be approximately 0.19(6) based only on the last two data points. Assuming it is 1/5 and using only the last data point, we obtain the remarkable 1-parameter fit depicted in Fig. 2a. The reason not to consider further corrective terms is that, strikingly, the error $|log ρ_{ℓ}^{-1} - π^2/(3p) + 2\pi\sqrt{2+\sqrt{2}}/\sqrt{p} - 10.153/\sqrt{p}|$ is less than 2 throughout the entire range $log p^{-1} \in [1, 12]$ and actually decays in the second half of this interval. In this interval, $log ρ_{ℓ}^{-1}$ increases from 2 to $5 \cdot 10^5$. This suggests that this might be the final divergent term in the asymptotics of $log ρ_{ℓ}^{-1}$ for FBP. This is supported, for instance by results in simpler but somewhat similar
<figure><img src="image_4.png" /><figcaption>(a) FBP: The solid orange line, the dashed red curve and the dotted magenta curve represent the first, second and third order asymptotics respectively. The first two are rigorously established.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) MBP: The dotted magenta line and the alternating dashed and dotted green curve represent the third and fourth order asymptotics respectively. Only the first order is known rigorously. The solid black box indicates the previously explored parameter range [40].</figcaption></figure>

<caption>FIG. 2: Exact numerical values of the local critical droplet density $\rho_{\ell}$ for $p = 2^{-k}$ with $k \in \{2, \dots, 17\}$ are given by the solid black dots. Empty blue circles represent Monte Carlo simulation data for $p^2\rho$ instead of $\rho_{\ell}$. The confidence intervals are smaller than the circles used in the figure.</caption>

settings [60], though the value of this exponent did not match our heuristic expectations.

A similar analysis allows us to detect the asymptotics

$$ \log \rho_{\ell} \sim \frac{-\pi^2}{3p} + \frac{\sqrt{2+\sqrt{2\log(1/p)}}}{\sqrt{p}} + \frac{6.54}{\sqrt{p}} \quad (2) $$

for local MBP, very clearly distinguishing the logarithmic correction (see Fig. 2b). Indeed, it is essential to consider richer functional expressions to obtain accurate results for this model [61]. In this case, we do expect at least one subsequent term with exponent around 0.19 to be neglected in (2). For reference, the first order term in (2) is known rigorously [37], the second one was recently conjectured [45, Section 8.2] based on theoretic considerations, while the third and subsequent order terms are new.

## WHAT WENT WRONG?

While the above results are rather satisfactory, for the sake of other models, it is important to reflect on what makes this method work, as opposed to the original local BP approach in [7, 39, 40]. Indeed, as far as the algorithm is concerned, disregarding implementation, the two are effectively identical. Yet, [40] mistakenly predicted a second term of order $p^{-2/3}$, while we achieve perfect agreement even to third order, based on the same kind of data.

So what went wrong twenty years ago? In fact, several things. Firstly, in the present approach, we do not artificially limit numerical computations to physically relevant scales. This may initially seem like a bad idea. However, we saw that even the last few data points, all corresponding to system sizes well beyond the size of the universe,

allow us to reliably extract very precise information, valid down to very small system sizes and consistent with rigorous asymptotic results previously accused of being irrelevant at physically meaningful scales [9]. For this reason it is essential that we take full advantage of the efficiency of the algorithm and its implementation. In other words, by considering large (and un-physical) scales, we effectively obtain more terms in the expansion of $\log \rho$ around infinity, which are crucial to understanding medium sized (and physically relevant) systems.

Secondly, in view of (1) and (2), one should allow some richness in the functional form to be fitted. Indeed, at our values of $p$, a logarithmic factor can easily be mistaken for a small power, if we do not take into account several data points to infer the correct shape. This is why, in [40], a second order exponent for MBP of 2/3 was estimated, instead of 1/2 with a logarithmic correction. It is thus important to consider fairly general expressions in order to test what kind of terms are actually needed, rather than artificially imposing a power law.

Finally, let us mention one last important aspect of the problem, for which the theoretical understanding at the time of [40] was not yet ripe. The locality theorem stated above implies that the non-locality $\log \rho - \log \rho_{\ell}$ is very small for FBP. However, for MBP, we expect this quantity to be much larger, of order $p^{-1/2+o(1)}$ [45, Section 8.2], as also supported by Figs. 1b and 2b. Thus, going beyond the precision of (2) for local MBP will not lead to more accurate results for non-local MBP. Hence, for models other than FBP and two-neighbour BP, for which locality is proved in [45], it is crucial to assess the validity and precision of the local approximation rather than applying it directly, as done in [7, 43].
CONCLUSION

We proposed a new viewpoint on BP based on rigorous quantitative locality. Local BP enabled us to make progress both on the rigorous theory and on the numerics of the problem. For the former, we proved exact second order asymptotics for FBP and two-neighbour BP, thanks to the strong locality theorem for these models. Concerning numerical methods, we have revived the approach pioneered in [7] for MBP, extended it to other models and implemented it efficiently. This allowed us to obtain a simple functional form in remarkable agreement with rigorous asymptotics, Monte Carlo simulations, in the regime where they are feasible, and exact numerical computations of the local density of critical droplets. We identified the caveats associated with the BP problem and developed a reliable methodology for avoiding them. This enabled a detailed understanding of the shortcomings of previous methods.

While the study of FBP and two-neighbour BP appears to be near its end, following this progress, this is far from being the case for other BP and related models, particularly those relevant for understanding glassy dynamics. It is therefore important that we have developed a detailed understanding of a few paradigmatic models, which can be used for benchmarking and calibrating future theoretical and numerical investigations. We also believe the iterative approach to correctly extract asymptotic expressions from data to be of potential use in a much broader context, whenever slowly decaying lower order terms may be present.

In conclusion, let us extend an invitation to numerical physics to reopen the study of BP and related models, taking into account the present locality approach. We hope that this will allow matching the recent theoretical understanding of sharp thresholds in two dimensions, as well as BP and related models in higher dimensions, within the universality framework and beyond it. To mention a few possible directions, it would be interesting to investigate to what extent other models are close to their local counterparts, to determine the leading asymptotics of so-called unbalanced critical models (see [55]) or practically determine the rather implicit leading asymptotics of balanced models (see [49]), as well as to explore even the standard models studied here in three spatial dimensions (see [46]).

*Acknowledgments* We thank Persi Diaconis, Alexander Holroyd, Cristina Toninelli and the anonymous references for helpful remarks and Paul Thévenin for enlightening discussions regarding [60]. I.H. was supported by the Austrian Science Fund (FWF): P35428-N. A.T. was supported by grants “Projeto Universal” (406250/2016-2) and “Produtividade em Pesquisa” (304437/2018-2) from CNPq and “Jovem Cientista do Nosso Estado”, (202.716/2018) from FAPERJ.

* hartarsky@math.univ-lyon1.fr

† augusto@impa.br

[1] D. Achlioptas, R. M. D'Souza, and J. Spencer, Explosive percolation in random networks, Science 323, 1453 (2009).

[2] O. Riordan and L. Warnke, Explosive percolation is continuous, Science 333, 322 (2011).

[3] M. Li, J. Wang, and Y. Deng, Explosive percolation obeys standard finite-size scaling in an event-based ensemble, Phys. Rev. Lett. 130, 147101 (2023).

[4] M. Pollak and I. Riess, Application of percolation theory to 2D-3D Heisenberg ferromagnets, Phys. stat. sol. (b) 69, K15 (1975).

[5] J. Chalupa, P. L. Leath, and G. R. Reich, Bootstrap percolation on a Bethe lattice, J. Phys. C 12, L31 (1979).

[6] P. M. Kogut and P. L. Leath, Bootstrap percolation transitions on real lattices, J. Phys. C 14, 3187 (1981).

[7] P. de Gregorio, A. Lawlor, P. Bradley, and K. A. Dawson, Clarification of the bootstrap percolation paradox, Phys. Rev. Lett. 93, 025501, 4 (2004).

[8] A. C. D. van Enter, J. Adler, and J. A. M. S. Duarte, Finite-size effects for some bootstrap percolation models, J. Stat. Phys. 60, 323 (1990).

[9] J. Adler and U. Lev, Bootstrap percolation: visualizations and applications, Braz. J. Phys. 33, 641 (2003).

[10] J. Gravner and A. E. Holroyd, Slow convergence in bootstrap percolation, Ann. Appl. Probab. 18, 909 (2008).

[11] P. de Gregorio, A. Lawlor, and K. A. Dawson, Bootstrap percolation, in *Encyclopedia of complexity and systems science*, edited by R. A. Meyers (Springer, New York, NY, 2009) pp. 608–626.

[12] P. Dehghanpour and R. H. Schonmann, Metropolis dynamics relaxation via nucleation and growth, Comm. Math. Phys. 188, 89 (1997).

[13] L. R. Fontes, R. H. Schonmann, and V. Sidoravicius, Stretched exponential fixation in stochastic Ising models at zero temperature, Comm. Math. Phys. 228, 495 (2002).

[14] R. Morris, Zero-temperature Glauber dynamics on $Z^d$, Probab. Theory Related Fields 149, 417 (2011).

[15] R. Cerf and F. Manzo, Nucleation and growth for the Ising model in $d$ dimensions at very low temperatures, Ann. Probab. 41, 3697 (2013).

[16] F. Ritort and P. Sollich, Glassy dynamics of kinetically constrained models, Adv. Phys. 52, 219 (2003).

[17] N. Cancrini, F. Martinelli, C. Roberto, and C. Toninelli, Kinetically constrained spin models, Probab. Theory Related Fields 140, 459 (2008).

[18] P. Garrahan, P. Sollich, and C. Toninelli, Kinetically constrained models, in *Dynamical heterogeneities in glasses, colloids and granular media and jamming transitions*, International series of monographs on physics 150, edited by L. Berthier, G. Biroli, J.-P. Bouchaud, L. Cipelletti, and W. van Saarloos (Oxford University Press, Oxford, 2011) pp. 341–369.

[19] F. Arceri, F. P. Landes, L. Berthier, and G. Biroli, A statistical mechanics perspective on glasses and aging, in *Encyclopedia of Complexity and Systems Science*, edited by R. A. Meyers (Springer, Berlin, Heidelberg, Berlin, Heidelberg, 2021) pp. 1–68.

[20] A. Guggiola and G. Semerjian, Minimal contagious sets
in random regular graphs, J. Stat. Phys. **158**, 300 (2015).

[21] E. Goles, P. Montealegre-Barba, and I. Todinca, The complexity of the bootstraping percolation and other problems, Theoret. Comput. Sci. **504**, 73 (2013).

[22] S. Banerjee, M. Jenamani, and D. K. Pratihar, A survey on influence maximization in a social network, Knowl. Inf. Syst. **62**, 3417 (2020).

[23] S. S. Manna, Abelian cascade dynamics in bootstrap percolation, Phys. A **261**, 351 (1998).

[24] J. Adler, Bootstrap percolation, Phys. A **171**, 453 (1991).

[25] R. Morris, Bootstrap percolation, and other automata, European J. Combin. **66**, 250 (2017).

[26] R. Morris, Monotone cellular automata, in *Surveys in combinatorics 2017*, London math. Soc. Lecture note ser., Vol. 440 (Cambridge University Press, Cambridge, 2017) pp. 312–371.

[27] J. Adler, S. Elfenbaum, and L. Sharir, Percolation in education and application in the 21st century, Eur. Phys. J. Spec. Top. **226**, 737 (2017).

[28] A. C. D. van Enter, Proof of Straley’s argument for bootstrap percolation, J. Stat. Phys. **48**, 943 (1987).

[29] J. Adler and A. Aharony, Diffusion percolation. I. Infinite time limit and bootstrap percolation, J. Phys. A **21**, 1387 (1988).

[30] R. H. Schonmann, On the behavior of some cellular automata related to bootstrap percolation, Ann. Probab. **20**, 174 (1992).

[31] M. Aizenman and J. L. Lebowitz, Metastability effects in bootstrap percolation, J. Phys. A **21**, 3801 (1988).

[32] R. Cerf and F. Manzo, The threshold regime of finite volume bootstrap percolation, Stochastic Process. Appl. **101**, 69 (2002).

[33] D. Kurtsiefer, Threshold value of three-dimensional bootstrap percolation, Int. J. Modern Phys. C **14**, 529 (2003).

[34] R. Lenormand and C. Zarcone, Growth of clusters during imbibition in a network of capillaries, in *Kinetics of aggregation and gelation*, edited by F. Family and D. P. Landau (Elsevier, Amsterdam, 1984) pp. 177–180.

[35] H. Nakanishi and H. Takano, Numerical study on the kinetic Ising model for glass transition, Phys. Lett. A **115**, 117 (1986).

[36] J. Adler, D. Stauffer, and A. Aharony, Comparison of bootstrap percolation models, J. Phys. A **22**, L297 (1989).

[37] A. E. Holroyd, Sharp metastability threshold for two-dimensional bootstrap percolation, Probab. Theory Related Fields **125**, 195 (2003).

[38] A. E. Holroyd, The metastability threshold for modified bootstrap percolation in $d$ dimensions, Electron. J. Probab. **11**, Paper No. 17, 418 (2006).

[39] P. de Gregorio, A. Lawlor, P. Bradley, and K. A. Dawson, Exact solution of a jamming transition: closed equations for a bootstrap percolation problem, Proc. Natl. Acad. Sci. USA **102**, 5669 (2005).

[40] P. de Gregorio, A. Lawlor, and K. A. Dawson, New approach to study mobility in the vicinity of dynamical arrest; exact application to a kinetically constrained model, EPL **74**, 287 (2006).

[41] J. Gravner and A. E. Holroyd, Local bootstrap percolation, Electron. J. Probab. **14**, Paper No. 14, 385 (2009).

[42] J. Gravner, A. E. Holroyd, and R. Morris, A sharper threshold for bootstrap percolation in two dimensions, Probab. Theory Related Fields **153**, 1 (2012).

[43] E. Teomy and Y. Shokef, Finite-density effects in

the Fredrickson-Andersen and Kob-Andersen kinetically-constrained models, J. Chem. Phys. **141**, 064110, 7 (2014).

[44] I. Hartarsky and R. Morris, The second term for two-neighbour bootstrap percolation in two dimensions, Trans. Amer. Math. Soc. **372**, 6465 (2019).

[45] I. Hartarsky and A. Teixeira, Bootstrap percolation is local, arXiv e-prints (2024), arXiv:2404.07903.

[46] J. Balogh, B. Bollobás, H. Duminil-Copin, and R. Morris, The sharp threshold for bootstrap percolation in all dimensions, Trans. Amer. Math. Soc. **364**, 2667 (2012).

[47] H. Duminil-Copin and A. C. D. van Enter, Sharp metastability threshold for an anisotropic bootstrap percolation model, Ann. Probab. **41**, 1218 (2013).

[48] B. Bollobás, H. Duminil-Copin, R. Morris, and P. Smith, The sharp threshold for the Duarte model, Ann. Probab. **45**, 4222 (2017).

[49] H. Duminil-Copin and I. Hartarsky, Sharp metastability transition for two-dimensional bootstrap percolation with symmetric isotropic threshold rules, Probab. Theory Related Fields **190**, 445 (2024).

[50] I. Hartarsky, F. Martinelli, and C. Toninelli, Sharp threshold for the FA-2f kinetically constrained model, Probab. Theory Related Fields **185**, 993 (2023).

[51] F. Martinelli, R. Morris, and C. Toninelli, Universality results for kinetically constrained spin models in two dimensions, Comm. Math. Phys. **369**, 761 (2019).

[52] L. Marêché, F. Martinelli, and C. Toninelli, Exact asymptotics for Duarte and supercritical rooted kinetically constrained models, Ann. Probab. **48**, 317 (2020).

[53] P. Balister, B. Bollobás, R. Morris, and P. Smith, Universality for monotone cellular automata, arXiv e-prints (2022), arXiv:2203.13806.

[54] I. Hartarsky and L. Marêché, Refined universality for critical KCM: lower bounds, Combin. Probab. Comput. **31**, 879 (2022).

[55] B. Bollobás, H. Duminil-Copin, R. Morris, and P. Smith, Universality for two-dimensional critical cellular automata, Proc. Lond. Math. Soc. (3) **126**, 620 (2023).

[56] P. Balister, B. Bollobás, R. Morris, and P. Smith, Subcritical monotone cellular automata, Random Structures Algorithms **64**, 38 (2024).

[57] I. Hartarsky, Refined universality for critical KCM: upper bounds, Comm. Math. Phys. **405**, Paper No. 13, 68 (2024).

[58] P. Balister, B. Bollobás, R. Morris, and P. Smith, The critical length for growing a droplet, Mem. Amer. Math. Soc. (To appear).

[59] In [45], the polylogarithmic bound on non-locality features the large exponent 20. This is only for the sake of simplifying the proof. We do not expect this exponent to be large. Moreover, in view of Fig. 2a, it is plausible that non-locality may be $2\log(1/p) + o(\log(1/p))$.

[60] A. Elvey Price, W. Fang, and M. Wallner, Compacted binary trees admit a stretched exponential, J. Combin. Theory Ser. A **177**, Paper No. 105306, 40 (2021).

[61] I. Hartarsky, Sensitive bootstrap percolation second term, Electron. Commun. Probab. **28**, Paper No. 29, 7 (2023).
# Supplemental Material: Locality approach to the bootstrap percolation paradox

This supplemental material documents the data generation and analysis discussed in the letter.

## EXACT COMPUTATION FOR LOCAL DENSITIES

In order to study the local infection model, one can make use of dynamic programming techniques to obtain the exact densities of critical droplets. The precise transition rates used in the algorithm can be found in [45, Table 1] for FBP, while for MBP they can be read off the source code file bootstrap/taichi/modified.py, available at the repository https://github.com/augustoteixeira bootstrap. The resulting data is gathered in Tab. I. Besides the efficient algorithm we have employed, two technical ingredients were necessary for us to reach such low values of $p$. First, we made use of a GPU to take advantage of the high parallelism present in the calculations. We ran this code on a Intel Xeon E3-1240 v3, with a single NVIDIA GeForce RTX 4060 8GB during 24 hours. Secondly, in order to deal with the very low numbers that appear in the calculations without losing precision, we store the logarithms of each probability ($a = \log(r)$), using the expression $\log(r + s) = \log(e^a + e^b) = \max\{a, b\} + \log(1 + e^{-|a-b|})$, in order to obtain the log of the sum or two small probabilities without losing precision. Based on this data, we progressively build the expressions in (1) and (2) and subsequent terms, assuming no prior knowledge for a start, see the scaling discussions below. The code for the various fits below can be found at https://github.com/augustoteixeira bootstrap in fitting.py.

<table><caption>**TABLE I:** Numerical data for local critical droplet densities.</caption><thead><tr><th>log<sub>2</sub> p<sup>-1</sup></th><th>log ρ<sub>ℓ</sub><sup>-1</sup>(p), FBP</th><th>log ρ<sub>ℓ</sub><sup>-1</sup>(p), MBP</th></tr></thead><tbody><tr><td>2</td><td>3.6462939089044335</td><td>3.025003004824336</td></tr><tr><td>3</td><td>9.48534315586599</td><td>6.778614767734161</td></tr><tr><td>4</td><td>24.785862065200995</td><td>17.63216670792452</td></tr><tr><td>5</td><td>61.09464730696058</td><td>45.66021724467772</td></tr><tr><td>6</td><td>142.44209408918184</td><td>112.51140378895116</td></tr><tr><td>7</td><td>318.20988111558466</td><td>263.82432820233373</td></tr><tr><td>8</td><td>689.051877876013</td><td>594.4366647140112</td></tr><tr><td>9</td><td>1458.978748960122</td><td>1299.3999937139088</td></tr><tr><td>10</td><td>3039.8354477597804</td><td>2776.8561233741584</td></tr><tr><td>11</td><td>6260.721269989636</td><td>5835.068010227808</td></tr><tr><td>12</td><td>12787.496049133362</td><td>12108.37054514238</td></tr><tr><td>13</td><td>25962.723826269754</td><td>24891.836236292278</td></tr><tr><td>14</td><td>52486.88654620641</td><td>50814.67137292014</td></tr><tr><td>15</td><td>105782.68428205681</td><td>103192.960947481</td></tr><tr><td>16</td><td>212726.28468173486</td><td>208743.86987754496</td></tr><tr><td>17</td><td>427113.5701763713</td><td>421026.7811793454</td></tr></tbody></table>

## MONTE CARLO SIMULATION OF NON-LOCAL DENSITIES

In order to contrast the results for the local model with those of the original Bootstrap Percolation dynamics, we have also performed Monte Carlo simulations of both FBP and MBP for moderately small values of $p$. In Figure 2, these results have been represented with blue circles. The code for such simulations can be found in the folder bootstrap/rust_mc/ of the repository https://github.com/augustoteixeira bootstrap and it does not involve any refined optimisation, besides using a low level programming language and making use of various cores in parallel. The results can be found in Tab. II.

## FIRST ORDER SCALING

We start by determining the power $\alpha$ such that $\log \rho_{\ell} \sim p^{-\alpha}$. While Fig. 2a clearly suggests that $\alpha \approx 1$, we would like to be more quantitative. As a first approximation, the last three discrete derivatives of $\log \log \rho_{\ell}^{-1}$ against $\log p^{-1}$ are 1.011, 1.008 and 1.006. The fact that they are close to 1, but varying suggests that an additional corrective term is needed. We attempt a four-parameter fit of the form $\log \log \rho_{\ell}^{-1} \sim \alpha \log p^{-1} + a + bp^{c}$, $c > 0$. We would like to use only large values of $1/p$, but need at least four of them. Using only the last five data points, for FBP, we obtain $\alpha \approx 0.9999(0)$. Moreover, this value is stable if we decide to use a different set (or number) of data points instead. This indicates that we have found the correct functional form for the leading term. Moreover, even though our guess for a pure power

<table><caption>**TABLE II:** Monte Carlo data for non-local critical droplet densities.</caption><thead><tr><th>log 5/2 p<sup>-1</sup></th><th>log ρ<sup>-1</sup>(p), FBP</th><th>log ρ<sup>-1</sup>(p), MBP</th></tr></thead><tbody><tr><td>10</td><td>1.66297868</td><td>0.88591602</td></tr><tr><td>11</td><td>1.90880149</td><td>0.93774197</td></tr><tr><td>12</td><td>2.31055326</td><td>1.20306221</td></tr><tr><td>13</td><td>2.9839127</td><td>1.38709304</td></tr><tr><td>14</td><td>3.7385026</td><td>1.50706354</td></tr><tr><td>15</td><td>5.12426155</td><td>1.96544661</td></tr><tr><td>16</td><td>6.8547289</td><td>2.71492423</td></tr><tr><td>17</td><td>8.76065061</td><td>3.726459</td></tr><tr><td>18</td><td>11.40863243</td><td>5.16653646</td></tr><tr><td>19</td><td>14.59271272</td><td>7.1772021</td></tr><tr><td>20</td><td>18.43952811</td><td>9.56943828</td></tr><tr><td>21</td><td></td><td>12.94289428</td></tr></tbody></table>
second term is wrong for MBP (see (2)), the above four-parameter fit still gives $\alpha \approx 0.999(7)$, as compared to the last discrete derivative of 1.012.

## FIRST ORDER CONSTANT

We have determined that $\alpha = 1$. We repeat the same fit with the remaining three parameters $\log \log \rho_{\ell}^{-1} \sim \alpha \log p^{-1} + a + bp^c$, $c > 0$ to find $a$. This leads to $a \approx 3.290(3)$ and this is also stable with respect to the choice of data points. Based on this and the rigorous lower bound on $\rho_{\ell}$ from [31] it is natural to guess that $a = \pi^2/3 \approx 3.28987$, which could have been done without the matching rigorous lower bound of [37]. The same approach gives $a \approx 3.29(2)$ for MBP.

## SECOND ORDER SCALING

Whether we guessed that the leading asymptotic term is $\pi^2/(3p)$ as above or we knew it from [37], our next task is to identify the functional form of the second order term. We set $\tilde{\rho} = \log \rho_{\ell}^{-1} - \pi^2/(3p)$ for convenience and we restart our analysis with $\tilde{\rho}$ instead of $\log \rho_{\ell}^{-1}$. We start by examining the discrete derivatives of $\log \tilde{\rho}$. For FBP, the last three are 0.514, 0.511, 0.509, while the middle one is 0.536. In comparison, the corresponding values for MBP are 0.579, 0.574 and 0.570, and the middle one is 0.609. This suggests a corrective term once again. However, in view of [61], it is important to allow a more general form for it. We perform a five parameter fit of the form $\log \tilde{\rho} \sim a \log p^{-1} + b \log \log p^{-1} + c + dp^e$, $e > 0$. With so many parameters, we face a problem: there are many local minima corresponding to spurious fits and local optimisers either get trapped in one of them or fail to converge, depending on initialisation. We solve this issue by first using more data points when fitting, say the last 10, to get a rough estimate of the relevant fit. We then narrow down the possible bounds for the five parameters and use only the last six data points to obtain an accurate prediction of their values. Finally, we verify that the resulting parameter values do produce a very good fit (mean squared errors are of order $10^{-13}$).

For FBP, this yields $a \approx 0.497$ and $b \approx 0.06$, which supports the correct conclusion that $a = 1/2$ and $b = 0$. Instead, for MBP, we obtain $a \approx 0.496$ and $b \approx 1.05$, confirming that $a = 1/2$ and $b = 1$. It should be noted that in this case optimisation is particularly hard, due to the presence of a slightly worse fit around $a \approx 0.508$, $b \approx 0.71$. Its mean squared error is only about 3 times larger, still of order $10^{-13}$ but it fluctuates more with the choice of data points, leading us to prefer the correct parameter choice. As we will see, this difficulty is due to our incorrect choice of third order correction.

We further remark that, if we had allowed a similar five-parameter fit for the first order term, the additional logarithmic correction is rejected by choosing a parameter close to 0 for it, as in the case of the second order term of FBP.

## SECOND ORDER CONSTANT

To determine the correct constant $c$ above, we simply fix the values $a = 1/2$ and $b = 0$ for FBP respectively and perform the fit on the remaining three variables. This leads to $\sqrt{\tilde{p}}\tilde{\rho} \sim 11.58(9)$ to be compared with the rigorous value $2\pi\sqrt{2+\sqrt{2}} \approx 11.6098$ from (1).

For MBP, proceeding similarly for $b = 1$ tends to prioritise values of $e$ close to 0. Indeed, as it is visible from Fig. 2b, we rather expect a third term of the form $p^{-1/2+o(1)}$, which means $e = 0$. This suggests that the four parameter fit $\sqrt{\tilde{p}}\tilde{\rho} \sim a \log p^{-1}+b+cp^d$, $d > 0$ is more appropriate. This gives $a \approx 1.848(4)$, again in agreement with the conjectured value $\sqrt{2+\sqrt{2}} \approx 1.84776$ from (2).

## THIRD ORDER SCALING

In order to investigate the third order term, it is convenient to define $\tilde{\rho} = \tilde{\rho} - 2\pi\sqrt{2+\sqrt{2}}/\sqrt{\tilde{p}}$ for FBP and $\tilde{\rho} = \sqrt{2+\sqrt{2}\frac{\log p^{-1}}{\sqrt{\tilde{p}}}} - \tilde{\rho}$ for MBP instead. For FBP, examining the discrete derivatives of $\log \tilde{\rho}$ reveals that they go down from $0.224$ at $p = 2^{-2}$, to $0.188$ at $p = 2^{-11}$ and then up to $0.195$ at $p = 2^{-16}$. Moreover, adding logarithmic corrections is unsuccessful, so we expect a pure power third order term with an exponent close to $1/5$.

For MBP, as seen from Fig. 2b, a constant with power law correction gives a good fit for $\tilde{\rho}\sqrt{\tilde{p}}$. On the other hand, proceeding as above, we could verify that more exotic possibilities such as polylogarithmic or iterated logarithmic leading terms deteriorate the result. Consequently, we expect that $\tilde{\rho}\sqrt{\tilde{p}}$ converges to a constant.

## THIRD ORDER CONSTANT

For FBP, already the straightforward third order constant obtained as $p^{1/5}\tilde{\rho}$ evaluated at the last data point yields an exceptionally good fit throughout the entire range we have access to (see Fig. 2a). We therefore focus on MBP and perform a three parameter fit $\tilde{\rho}\sqrt{\tilde{p}} \sim a+bp^c$, $c > 0$. Using the last four data points (again, the result is not sensitive to this choice). This gives $a \approx 6.539$, $b \approx -7.4(4)$ and $c \approx 0.31$. This also leads to a convincing fit for all available values of $p$ (see Fig. 2b). Interestingly, the value of the exponent $c$ suggests a fourth order term in MBP with exponent very close to the one of the third order term of FBP.
NON-LOCALITY

We now explain in more detail the reason for compar-
ing $\rho_\ell$ with $p^2\rho$ and assess non-locality based on the com-
parison between local and non-local simulations. Firstly,
the samples displayed in Fig. 1 are not specially selected,
but simply the first ones found. The one in Fig. 1a hap-
pens to be local in the sense that there exists an initially
infected starting $1 \times 1$ square for which the local dynam-
ics yields the final $93 \times 93$ square. This suggests that
many internally filled critical droplets are actually local.

Let us assume for a moment that all internally filled
critical droplets are local and use this assumption to ar-
gue that $p^2\rho \sim \rho_\ell$. In this case, the event that the droplet
is internally filled is the union over starting points of the
probability that the local dynamics started there pro-
duces the final square. Thinking of a typical starting
point, this probability is more or less $\rho_\ell$, itself essentially
equal to the probability of growing to infinity from this
starting point. If the events corresponding to different
starting points were disjoint and had probability exactly
$\rho_\ell$, we would have $\rho = \Lambda^2\rho_\ell$.

This is only approximately true for multiple reasons.
Firstly, starting points close to the boundary are less
likely to succeed, because their growth is constraint to go
in a particular direction from early on. Secondly, events
are not quite disjoint. Indeed, there are always at least
three valid starting points contained in a 2×2 square (see
Fig. 1a). Sometimes there may be more possible starting
points, but rarely many and they are always clustered

together (it is not hard to prove along the lines of [45]
that distant successful starting points are asymptotically
unlikely). Recall that $\Lambda^2 = 4 \log^2(1/p)/p^2$ is a bit big-
ger than $1/p^2$. We neglect the remaining polylogarithmic
factor in order to balance out the two effects above, which
are difficult to quantify, but certainly at most polyloga-
rithmic as well. This leads to the desired $\rho_\ell \sim p^2\rho$, up
to logarithmic corrections, which have negligible effect
already at the largest scales accessible to Monte Carlo
simulations.

The perfect fit already at moderate values of *p* be-
tween $\rho_{\ell}$ and $p^{2}\rho$ visible in Fig. 2a further confirms the
correctness of this reasoning and suggests that in FBP,
locality is even stronger than what is proved in the local-
ity theorem. For comparison, we have applied the same
procedure to MBP. Here, already Fig. 1b displays two
non-local moves necessary for internally filling the criti-
cal droplet: for creating a 6 × 9 rectangle from a 4 × 7
one and then a 10 × 13 rectangle from a 8 × 11 one. If,
despite this warning of non-locality, we compare $p^{2}\rho$ and
$\rho_{\ell}$ for MBP as in Fig. 2b, we observe that the two are
not getting closer as $p$ decreases, but on the contrary.
Moreover, the plot suggests that $(\log \rho_{\ell} - \log \rho)\sqrt{p}$ does
not decay to 0. While the data seems insufficient to ad-
vance more precise predictions, let us point out that this
is in agreement with the heuristics of [45, Section 8.2].
In any case, the evidence is clear that locality for MBP
is much weaker and only arises once droplets reach some
intermediate scale.