# Contact Kinetics in Fractal Macromolecules

Maxim Dolgushev, Thomas Guérin, Alexander Blumen, Olivier Bénichou,
Raphaël Voituriez

► To cite this version:

Maxim Dolgushev, Thomas Guérin, Alexander Blumen, Olivier Bénichou, Raphaël Voituriez. Contact Kinetics in Fractal Macromolecules. Physical Review Letters, 2015, 115 (20), pp.208301 (1-5). 10.1103/PhysRevLett.115.208301. hal-01251841

HAL Id: hal-01251841

https://hal.science/hal-01251841

Submitted on 6 Jan 2016

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_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
Contact Kinetics in Fractal Macromolecules

Maxim Dolgushev,¹ Thomas Guérin,² Alexander Blumen,¹ Olivier Bénichou,³ and Raphaël Voituriez³

¹Physikalisches Institut, Universität Freiburg, Hermann-Herder-Str. 3, 79104 Freiburg, Germany

²Université de Bordeaux and CNRS, Laboratoire Ondes et Matière d'Aquitaine (LOMA), UMR 5798, 33400 Talence, France

³Laboratoire de Physique Théorique de la Matière Condensée,
CNRS/UPMC, 4 Place Jussieu, 75005 Paris, France

(Dated: October 27, 2015)

We consider the kinetics of first contact between two monomers of the same macromolecule. Relying on a fractal description of the macromolecule, we develop an analytical method to compute the Mean First Contact Time (MFCT) for various molecular sizes. In our theoretical description, the non-Markovian feature of monomer motion, arising from the interactions with the other monomers, is captured by accounting for the non-equilibrium conformations of the macromolecule at the very instant of first contact. This analysis reveals a simple scaling relation for the MFCT between two monomers, which involves only their equilibrium distance and the spectral dimension of the macromolecule, independently of its microscopic details. Our theoretical predictions are in excellent agreement with numerical stochastic simulations.

PACS numbers: 82.35.Lr,05.40.-a,36.20.Ey,82.20.Uv

*Introduction.* Intramolecular reactions are ubiquitous in nature. Examples are provided by the formation of RNA hairpins [1] or DNA loops [2], the folding of polypeptides [3], as well as the appearance of cycles in synthetic polymers [4]. It is generally known that the structure of a macromolecule has a strong influence on the dynamics of its monomers [5, 6], and that this complex intramolecular dynamics implies nontrivial reaction kinetics [7] in the diffusion limited regime. Until now, most of the theoretical work devoted to the reaction times in macromolecules has been limited mainly to linear chains [7–16]. However, there are numerous examples of macromolecules which differ from linear polymer chains [17, 18]; their dynamic and static properties often suggest a fractal character [6, 18, 19]. In particular, the dynamics of fractal macromolecules is characterized by dynamical exponents that are different from those of linear chains [6, 19–21], leading presumably to distinct first contact kinetics that cannot be deduced from existing works on linear chains. It is important to note the significance of fractals: They provide typical models, e.g., for hyperbranched polymers [6], proteins [21, 22], sol-gel branched clusters [23], and colloidal aggregates [24].

The aim of this Letter is to propose a theoretical description of the mean first contact time (MFCT) between two monomers in a fractal macromolecule, described as a network of beads connected by springs. The cornerstone feature of such fractal models is the anomalous vibrational dynamics of the network [6, 20, 21]. It originates from the non-Debye density of states which is characterized through the so-called spectral dimension $d_s$ [25] (also known as "fracton" dimension [26]). In this Letter, we go beyond existing studies that focused on the specific case of linear chains only [7–16], and show on general grounds that $d_s$ (rather than the microscopic properties of the macromolecule) is the key parameter that controls intramolecular reaction kinetics. Indeed, the spectral dimension $d_s$ will be shown to leave its fingerprint in the

scaling behavior of MFCTs with the equilibrium distance between monomers, as will be confirmed by explicit computations on examples of fractal structures (see Fig. 1).

It is important to stress that in the case of macromolecules, the interactions between monomers lead to an effective non-Markovian motion, which is the hallmark of monomer dynamics [27]. In this Letter, we take into account such non-Markovian features and describe the contact kinetics for fractal structures. We show explicitly that the non-Markovian effects increase with the complexity and the degree of branching of the macromolecules.

*Model.* The macromolecular structures are represented by N beads located at positions $\mathbf{r}_i(t)$ in a 3-dimensional space and connected by springs of stiffness K. The

<figure><img src="image_2.png" /><figcaption>FIG. 1. (color online) Structure of fractal macromolecules investigated in this Letter: (a) Vicsek fractal, here of functionality (i.e. number of nearest-neighbors of the branching sites) $f = 4$ (VF4), (b) dual Sierpiński gasket (DSG), (c) T–fractal (TF). The reactive monomers for which we compute the MFCT are represented by red squares. These extended conformations show only the topology of the structures.</figcaption></figure>
free-draining dynamics of the structure is given by the Langevin equations [5, 6]:

$$ \zeta \frac{\partial}{\partial t} \mathbf{r}_i(t) + K \sum_{j=1}^{N} A_{ij} \mathbf{r}_j(t) = \mathbf{F}_i(t), \quad (1) $$

where $\mathbf{A} = (A_{kj})$ is the connectivity (Laplacian) matrix that describes the topology of the structure. The off-diagonal elements $A_{ij}$ are equal to $-1$ if beads $i$ and $j$ are connected, and 0 otherwise. For each bead $i$, the diagonal element $A_{ii}$ is equal to the number of bonds emanating from it. Also Eq. (A3) includes friction forces $-\zeta \partial_t \mathbf{r}_i$ and stochastic forces $\mathbf{F}_i(t)$ obeying white noise statistics with amplitude $\langle F_{i\alpha}(t) F_{j\beta}(t') \rangle = 2k_B T \zeta \delta(t-t') \delta_{ij} \delta_{\alpha\beta}$, where $k_B T$ represents the thermal energy and $\alpha, \beta$ are spatial coordinates $x, y, z$. It is natural to introduce the monomeric relaxation time $\tau_0 = \zeta/K$, and the characteristic microscopic length $l = (3k_B T/K)^{1/2}$ (in the case of structures without loops $l^2$ is the mean-squared bond length).

Here we consider the contact kinetics between two given monomers denoted "reactive monomers", whose indexes are called $q_1$ and $q_2$ (see Fig. 1). We introduce the vector $\mathbf{R}(t)$ that joins them:

$$ \mathbf{R}(t) = \mathbf{r}_{q_1}(t) - \mathbf{r}_{q_2}(t) = \sum_{i=1}^{N} h_i \mathbf{r}_i(t), \quad (2) $$

where $h$ is a $N$-dimensional vector defined by this equation, which has only 2 non-zero elements in positions $q_1$ and $q_2$. It is convenient to decompose the Gaussian vector $\mathbf{R}(t)$ as a sum of independent modes,

$$ \mathbf{R}(t) = \sum_{\lambda} b_{\lambda} \mathbf{a}_{\lambda}(t), \quad (3) $$

where $\lambda$ represents all the distinct nonvanishing eigenvalues of $\mathbf{A}$, $b_{\lambda}^2$ is the norm of the orthogonal projection of the vector $\mathbf{h}$ [defined in Eq. (A1)] on the eigensubspace associated to $\lambda$ [28]. The $\mathbf{a}_{\lambda}(t)$ evolve independently of each other with the correlation function:

$$ \langle a_{\lambda, \alpha}(t) a_{\lambda', \beta}(t') \rangle = l^2 \delta_{\alpha, \beta} \delta_{\lambda, \lambda'} e^{-\lambda|t-t'|/\tau_0} / (3\lambda). \quad (4) $$

In this picture of independent modes, the normalized temporal autocorrelation function of the Cartesian components ($R_x(t), R_y(t), R_z(t)$) of vector $\mathbf{R}(t)$ follows from Eqs. (3-4):

$$ \phi(t) = \frac{\langle R_{\alpha}(t) R_{\alpha}(0) \rangle}{\langle R_{\alpha}(0)^2 \rangle} = \left( \sum_{\lambda} \frac{b_{\lambda}^2}{\lambda} e^{-\lambda t / \tau_0} \right) / \sum_{\lambda} \frac{b_{\lambda}^2}{\lambda}. \quad (5) $$

Since $\mathbf{R}(t)$ is a Gaussian stochastic process, it is entirely characterized by the function $\phi(t)$. We also introduce $n = \langle [\mathbf{R}(t)/l]^2 \rangle = \sum_{\lambda} b_{\lambda}^2 / \lambda$, which for macromolecules devoid of loops can be shown to be simply the number of bonds connecting the reactive monomers.

Theories of first contact times. We now sketch briefly the method we use for the calculation of the MFCT,

<figure><img src="image_3.png" /><figcaption>FIG. 2. (color online) $T_{NM}$ (Eq. (8)) and $T_{WF}$ (Eq. (C1)) for Vicsek fractals of functionality $f = 4$ as a function of the capture radius $a$. The lines represent the results of Eq. (9).</figcaption></figure>

defined as the average time needed for the reactive monomers to be separated by a distance smaller than $a$ (called the capture radius), starting from an initial equilibrium configuration in which the reactive monomers are not in contact. We introduce the joint probability density $f(\{\mathbf{a}\}, t)$ that contact is made for the first time at time $t$ and that, at this first passage event, the macromolecule has a configuration described by the set of modes $\{\mathbf{a}\} = (\mathbf{a}_1, \mathbf{a}_2...)$ Let us partition the trajectories that lead to a configuration $\{\mathbf{a}\}$ (in which the contact condition is satisfied) into two steps, the first step consisting in reaching the target for the first time at $t'$, and the second step consisting in reaching the final configuration $\{\mathbf{a}\}$ in a time $t - t'$. The mathematical formulation that corresponds to this decomposition of events is

$$ p(\{\mathbf{a}\}, t) = \int_{0}^{t} dt' \int d\{\mathbf{a}'\} f(\{\mathbf{a}'\}, t') p(\{\mathbf{a}\}, t - t'|\{\mathbf{a}'\}), \quad (6) $$

where $p(\{\mathbf{a}\}, t|\{\mathbf{a}'\})$ is the probability of $\{\mathbf{a}\}$ at $t$ starting from $\{\mathbf{a}'\}$ at $t = 0$ while $p(\{\mathbf{a}\}, t)$ is the probability of $\{\mathbf{a}\}$ at $t$ starting from the initial conditions. Although Eq. (6) is exact, it is in general very difficult to solve. A classical approximation, introduced by Wilemski and Fixman (WF) [8], assumes that $f(\{\mathbf{a}\}, t)$ is proportional to the equilibrium distribution of configurations that satisfy the constraint that a contact is formed. Introducing this approximation into Eq. (6), integrating over all configurations and taking the long time limit lead to the estimate $T_{WF}$ of the mean first contact time [15]

$$ T_{WF} = \int_{0}^{\infty} dt \left\{ \frac{e^{-a^2 \phi(t)^2 / [2\psi(t)]}}{[1 - \phi(t)^2]^{3/2}} - \frac{Z(a/\psi(t))^{1/2}}{Z(a\sqrt{3}/l\sqrt{n})} \right\}, \quad (7) $$

where $\psi(t) = nl^2[1 - \phi(t)^2]/3$ is the mean-square displacement of $R_x(t)$ and $Z(y) = \int_y^\infty dx \, x^2 e^{-x^2/2}$. While this approximation takes into account some aspects of the complex dynamics of the macromolecule through the correlation function $\phi$, it neglects non-Markovian (NM)
effects which can be quantitatively important. These NM effects can be described by considering the distribution of configuration at the first contact event $\pi(\{\mathbf{a}\}) = \int_0^\infty dt f(\{\mathbf{a}\}, t)$. The analytic expression for $\pi(\{\mathbf{a}\})$ is unknown; following the case of linear polymers [12], we assume that it is well approximated by a Gaussian distribution, which is therefore characterized by its first and second moments. Its first moments, denoted $m_\lambda$, are the average mode amplitudes $\mathbf{a}_\lambda$ at the first contact instant $t^*$ (in the direction of the vector $\mathbf{R}(t^*)$), while the second moments can be approximated by their equilibrium values. Then, a precise estimate of the MFCT is obtained by integrating (6) over all contact configurations, leading to

$$T_{\text{NM}} = \int_{0}^{\infty} dt \left\{ \frac{e^{-R_{\pi}(t)^2/[2\psi(t)]}}{[1 - \phi(t)^2]^{3/2}} - \frac{Z(a/\psi(t)^{1/2})}{Z(a\sqrt{3}/l\sqrt{n})} \right\}. \quad (8)$$

The difference between this expression and Eq. (C1) lies in the presence of the reactive trajectory $R_\pi(t)$, defined as the average of $\mathbf{R}$ at a time $t$ after the first contact $t^*$ in the direction of $\mathbf{R}(t^*)$, and thus reads $R_\pi(t) = \sum_\lambda b_\lambda m_\lambda e^{-\lambda t/\tau_0}$. The involved $m_\lambda$ are obtained from a set of self-consistent equations, see the Supplemental Material [28] for details.

**Cyclization for small capture radius.** In the limit of $a \to 0$ (while keeping a fixed number of monomers), only the dynamics at small time scales matters, where $\psi(t) \simeq 4D_0t$ with the local diffusion constant $D_0 = k_B T/\zeta$. Also, at short times, we can write $R_\pi(t) \simeq R_\pi(0) = a$. Introducing these approximations into Eqs. (C1,8), we deduce that both the WF and the NM theories predict a MFCT which reads

$$T \simeq \tau_0 \pi^{1/2} n^{3/2} l/(2\sqrt{6} a) \quad (\text{for } a \to 0). \quad (9)$$

Thus, in this regime the MFCT is independent of the particular polymeric structure, and has the same expression as in the case of linear chains [12, 15]. The convergence of both Eqs. (C1,8) to the scaling form (9) for small capture radius is demonstrated in Fig. 2.

**Scaling for large $n$.** When the size of the structure grows, Eq. (9) is not valid anymore and the MFCT reflects then the complex monomer dynamics, which becomes subdiffusive, $\langle[\mathbf{r}_q(t) - \mathbf{r}_q(0)]^2\rangle \sim t^\gamma$, with $\gamma$ a subdiffusive exponent related to the spectral dimension $d_s$ of the structure, $\gamma = 1 - d_s/2$ for $d_s < 2$ [20, 25, 26, 29]. We first present a simple qualitative argument to derive the behavior of the MFCT as a function of $n$, which we explicitly check for the WF theory and validate next for the NM theory by a numerical analysis of the analytical results of Eq. (8).

For a Markovian subdiffusive walker whose subdiffusive exponent is $\gamma$, the time needed to find a point-like target in a confining volume of size $L$ starting at a random position is proportional to $L^{2/\gamma}$ [30]. Given that in our case the typical length scale is $L \sim \sqrt{n}$, the scaling argument

$$T \sim \tau_0 n^{1/\gamma} \quad (10)$$

<figure><img src="image_4.png" /><figcaption>FIG. 3. (color online) (a) NM and (b) WF MFCT for different fractal structures as a function of $n^{1/\gamma}$. The parameter $\gamma$ is obtained from the known values [20, 32, 33] of the spectral dimension $d_s$: for Vicsek fractals of functionality $f$, $\gamma = \ln(3)/\ln(3f + 3)$, for dual Sierpiński gasket $\gamma = \ln(5/3)/\ln(5)$, and for T-fractal $\gamma = \ln(2)/\ln(6)$. The capture radius is $a = l$. (c) NM MFCT for the dual Sierpiński gasket for different values of $a$.</figcaption></figure>

follows. We do not expect a strong dependence on the capture radius since the dynamics of a monomer is compact (or recurrent) [7, 30]. In fact, the scaling (10) can be derived in the WF approximation by noting that for extremal monomers (i.e. monomers whose relative distance is maximal), the time $\tau_0 n^{1/\gamma}$ is also of the order of the maximal relaxation time $\tau_N$ [26, 31], such that the correlation function scales as $\phi(t) = \tau_N \Phi(t/\tau_N)$. Once this equality is reported into Eq. (C1) one obtains that $T \sim \tau_N$, which is the expected behavior (10), see the Supplemental Material [28] for details. Strikingly, the MFCT is found to depend on the polymer structure through $d_s$ only (and not on its microscopic details). In particular, this shows that the MFCT of highly branched structures ($\gamma \to 0$) differs significantly from that of linear chains ($\gamma = 1/2$).

Now, checking the scaling (10) requires the actual computation of the correlation function $\phi(t)$, which itself involves the diagonalization of $\mathbf{A}$. However, a naive diagonalization of this matrix (as usually done for linear chains, where the analytic diagonalization is possible) does not allow to deal with structures large enough. This difficulty can be overcome by exploiting the highly symmetric nature of the fractal macromolecules we are considering. First, we remark that the actual number of
<table><caption>TABLE I. MFCT for the structures of Fig. 1 computed with the WF ($T_{WF}$) and NM ($T_{NM}$) theories [with Eqs. (C1,8)], compared to the results $T_{simu}$ of the stochastic simulations of (A3). $1/\gamma$ indicates numerical values of the scaling exponent in (10) (in the case of linear chains it is equal to 2 [8, 9, 12]). All structures are of generation $g = 3$. The times are units of $\tau_0$.</caption><thead><tr><th>a/l</th><th>structure</th><th>1/γ</th><th>T<sub>WF</sub></th><th>T<sub>NM</sub></th><th>T<sub>simu</sub></th></tr></thead><tbody><tr><td rowspan="3">1</td><td>DSG</td><td>3.15</td><td>3.85</td><td>2.38</td><td>2.41 ± 0.02</td></tr><tr><td>TF</td><td>2.58</td><td>13.31</td><td>9.48</td><td>9.53 ± 0.07</td></tr><tr><td>VF3</td><td>2.26</td><td>166.8</td><td>111.5</td><td>113.6 ± 0.9</td></tr><tr><td rowspan="3">2.5</td><td>VF4</td><td>2.46</td><td>224.9</td><td>135.1</td><td>137.4 ± 1.4</td></tr><tr><td>VF3</td><td>2.26</td><td>87.81</td><td>37.52</td><td>37.47 ± 0.35</td></tr><tr><td>VF4</td><td>2.46</td><td>128.3</td><td>47.1</td><td>45.3 ± 0.5</td></tr></tbody></table>

variables needed to be taken into account is the number of *distinct* eigenvalues, which is much lower than the total number of beads $N$. Then, one can set up a decimation procedure inspired by that used to find iterative formulas for the eigenvalues [20, 34]. Adapting this decimation approach to the iterative computation of the coefficients $b_{\lambda}^2$ amounts to solving a linear algebra problem, as described in the Supplemental Material [28]. In practice, by using the WF and NM formalisms we were able to calculate the MFCT for macromolecules containing as many as 800,000 beads, which would not have been possible through direct diagonalization.

Based on the iterative computation of $b_{\lambda}^2$, we test the scaling behavior (Fig. 3) for different fractal structures (represented in Fig. 1) and therefore for several values of the subdiffusive exponent $\gamma$ (or $d_s$). As can be seen on Fig. 3, for all these structures, the scaling $T \sim n^{1/\gamma}$ is in good agreement with the predictions of both the NM theory and the WF theory. This confirms that the spectral dimension $d_s$ of the structure plays a fundamental role for MFCT. In particular, the functionality $f$, which determines $\gamma$, plays a crucial role and yields scaling behaviors that can differ significantly from those of linear chains. However, the presence of many loops in the dual Sierpiński gasket does not modify the scaling behavior of the MFCT: In Fig. 3(c) we show the results of the corresponding MFCT $T_{NM}$ for 3 different values of $a$. We observe that all curves scale in the same way for large $n$, independently of the capture radius, fact consistent with Eq. (10).

**Comparison with numerical simulations.** We checked the validity of our theoretical predictions by performing Brownian dynamics simulations of Eq. (A3), using the algorithm of Ref. [35] with fixed time steps. The results are presented in Table I for different fractal structures and capture radius $a$. We observe that the WF theory systematically overestimates the MFCT, whereas the NM theory describes the simulation data almost quantitatively, thereby validating the accuracy of our analysis.

**Average conformations at first contact.** Inspecting Fig. 3 reveals that the WF theory overestimates the

<figure><img src="image_5.png" /><figcaption>FIG. 4. (color online) Average spatial position $\langle z_i \rangle_{\pi}$ of the monomers (numbered through $i = 1...n$) connecting reactants at the instant cyclization in the direction of the reaction, for Vicsek fractals of different functionalities $f$ with generation $g = 3$, as predicted by the NM (symbols) and WF theories (solid line). Inset: same quantity determined from simulations, with the same color code.</figcaption></figure>

MFCT by a numerical factor which grows when $\gamma$ decreases. This fact is also clearly seen in Table I, and confirmed by the numerical simulations. This means that, with decreasing subdiffusive exponent $\gamma$, the macromolecular conformations at the instant of first contact differ more and more from equilibrium ones. To illustrate this fact, we present on Fig. 4 the average spatial positions of the monomers in Vicsek fractals at the instant of first contact $\langle z_i \rangle_{\pi} \equiv \langle \mathbf{r}_i(t^*) \cdot \hat{\mathbf{u}}(t^*) \rangle$, where $t^*$ represents the first contact instant, $\hat{\mathbf{u}}(t^*) = \mathbf{R}(t^*)/|\mathbf{R}(t^*)|$ gives the direction between the reactive monomers at $t^*$ and $i$ are the indices of the beads that lie between the reactive monomers (on the chain). In the equilibrium WF theory, all monomers between the reactants (on the chain) are on average also between the reactants in space at the instant of first contact, irrespectively of the functionality of the macromolecule. The NM theory, in turn, predicts that those monomers which are close to the reactants are in average outside the capture radius at $t^*$ (Fig. 4, symbols). Moreover, this effect is more pronounced for Vicsek fractals of higher functionalities, meaning that NM effects increase with the degree of hyperbranching of the macromolecule. As shown in the inset of Fig. 4, the simulations confirm these conclusions, although the value of $\langle z_i \rangle_{\pi}$ is slightly overestimated in the NM theory.

**Conclusion.** Summarizing, we have studied the kinetics of first contact between two monomers belonging to the same fractal macromolecule. We identified two regimes: (i) for very small capture radius, the MFCT becomes independent of the macromolecular structure and originates essentially from the microscopic diffusive motion of the monomers, whereas (ii) for larger capture radius, the MFCT scales as a power-law with the mean-square distance $n$ between the monomers, with an exponent related to the spectral dimension $d_s$ and indepen-
dent of microscopic details. We confirmed this scaling law for a wide variety of structures (with and without loops). The non-Markovian effects are included by calculating the average equilibrium conformation of the whole macromolecule at the instant of first contact, and are found to be more important when the degree of hyperbranching of the structures increases. Finally, it would be interesting to explore how the MFCT varies with the location of the reactants in the structure or to generalize the theory to include the effect of hydrodynamic interactions, which typically would lead to new scalings [31].

## ACKNOWLEDGMENTS

*Acknowledgments.* M.D. and A.B. acknowledge the support of the DAAD through the PROCOPE program (Project No. 55853833), of the DFG through Grant No. Bl 142/11-1 and through IRTG Soft Matter Science (GRK 1642/1). T.G., O.B., and R.V. acknowledge the support of Campus France (project No. 28252XE) and of the European Research Council starting Grant FPTOpt-277998.

## Supplementary information

## Appendix A: Determination of the autocorrelation functions

Here we present general expressions for the determination of the temporal autocorrelation function of the vector $\mathbf{R}(t)$ connecting the reactive monomers $q_1$ and $q_2$, defined as

$$ \mathbf{R}(t) = \mathbf{r}_{q_1}(t) - \mathbf{r}_{q_2}(t) \equiv \sum_{i=1}^{N} h_i \mathbf{r}_i(t). \quad (A1) $$

where $\mathbf{h} = (h_i)$ is the $N$-components vector

$$ h_i = \delta_{i,q_1} - \delta_{i,q_2}. \quad (A2) $$

We recall that the positions $\{\mathbf{r}_i(t)\}$ of the $N$ beads in the macromolecule obey the Langevin equation [Eq.(1) in the main text]:

$$ \zeta \frac{\partial}{\partial t} \mathbf{r}_i(t) + K \sum_{j=1}^{N} A_{ij} \mathbf{r}_j(t) = \mathbf{F}_i(t), \quad (A3) $$

where $\mathbf{A}$ is the connectivity matrix that describes the topology of the structure. The off-diagonal elements $A_{ij}$ are equal to $-1$ if beads $i$ and $j$ are connected, and to 0 otherwise. For each bead $i$, the diagonal element $A_{ii}$ is equal to the number of bonds emanating from it. Finally, in Eq. (A3), the stochastic forces $\mathbf{F}_i(t)$ follow the statistics

$$ \langle F_{i\alpha}(t) F_{j\beta}(t') \rangle = 2 k_B T \zeta \delta(t-t') \delta_{ij} \delta_{\alpha\beta}, \quad (A4) $$

where $k_B T$ represents the thermal energy and $\alpha, \beta$ are the spatial coordinates $x, y, z$. It is natural to introduce the microscopic time scale $\tau_0 = \zeta/K$.

We introduce the eigenvectors of the matrix $\mathbf{A}$ and denote by $u_i^{(\lambda,q)}$ the $i$th coordinate of the $q$th eigenvector associated with the eigenvalue $\lambda \neq 0$. If $\lambda$ has degeneracy $p_\lambda$, then $q$ is bounded by $1 \le q \le p_\lambda$. Since the matrix $\mathbf{A}$ is symmetric and real, we can assume that these eigenvectors are orthogonal and normalized. It is then standard to decompose the motion of each bead in the structure in a sum of eigenmode amplitudes,

$$ \mathbf{r}_i(t) = \sum_{\lambda} \sum_{q=1}^{p_{\lambda}} u_i^{\lambda,q} \tilde{\mathbf{a}}_{\lambda,q}(t). \quad (A5) $$

Inserting (A5) into (A3) and taking into account the fact that the vectors $u_i^{\lambda,q}$ are orthogonal, we obtain that all the eigenmode amplitudes evolve independently of each other, with the correlation function

$$ \langle \tilde{a}_{\lambda,q,\alpha}(t) \tilde{a}_{\lambda',q',\beta}(t') \rangle = \frac{k_B T}{\lambda K} \delta_{\alpha,\beta} \delta_{q,q'} \delta_{\lambda,\lambda'} e^{-\lambda|t-t'|/\tau_0}. \quad (A6) $$

Let us introduce now, for each eigenvalue $\lambda$, the positive coefficient $b_\lambda$ and the stochastic variable $\mathbf{a}_\lambda(t)$ defined as

$$ b_{\lambda}^{2} = \sum_{q=1}^{p_{\lambda}} \left( \sum_{i=1}^{N} h_{i} u_{i}^{(\lambda,q)} \right)^{2}, \quad (A7) $$

$$ \mathbf{a}_{\lambda}(t) = \frac{1}{b_{\lambda}} \sum_{q=1}^{p_{\lambda}} \sum_{i=1}^{N} h_i u_i^{(\lambda,q)} \tilde{\mathbf{a}}_{\lambda,q}(t). \quad (A8) $$

Comparing (A8) and (A6), it is clear that the autocorrelation function of the stochastic variables $\mathbf{a}_\lambda$ is

$$ \langle a_{\lambda,\alpha}(t) a_{\lambda',\beta}(t') \rangle = \frac{k_B T}{K\lambda} \delta_{\alpha,\beta} \delta_{\lambda,\lambda'} e^{-\lambda|t-t'|/\tau_0}, \quad (A9) $$

which is exactly Eq. (4) of the main text. Furthermore, it follows from (A1,A5,A7) that

$$ \mathbf{R}(t) = \sum_{\lambda} \sum_{q=1}^{p_{\lambda}} \sum_{i=1}^{N} h_i u_i^{(\lambda,q)} \tilde{\mathbf{a}}_{\lambda,q}(t) = \sum_{\lambda} b_{\lambda} \mathbf{a}_{\lambda}(t) \quad (A10) $$

which is Eq. (3) of the main text. Finally, we note that (A7) can be written as

$$ b_{\lambda}^{2} = \sum_{q=1}^{p_{\lambda}} (\mathbf{h} | \mathbf{u}^{(\lambda,q)})^{2} = (\mathbf{h} | \hat{P}_{\lambda} | \mathbf{h}), \quad (A11) $$

where we have considered $\mathbf{h}$ and $\mathbf{u}^{(\lambda,q)}$ as $N$-components vectors, and where the projection operator $\hat{P}_\lambda$ is

$$ \hat{P}_{\lambda} = \sum_{q=1}^{p_{\lambda}} |\mathbf{u}^{(\lambda,q)})(\mathbf{u}^{(\lambda,q)}|. \quad (A12) $$
Thus, $b_{\lambda}^2$ is the norm of the orthogonal projection of the $N$-component vector $\mathbf{h}$ on the eigensubspace associated to $\lambda$, as claimed in the main text.

For the practical calculation of $b_{\lambda}^2$, it is useful to consider a basis of vectors $\{|\mathbf{w}^{(\lambda,q')}\}$ which span the orthogonal complementary eigensubspace associated to $\lambda$ (i.e. $(\mathbf{w}^{(\lambda,q')}|\mathbf{u}^{(\lambda,q)}) = 0$). The $N-p_{\lambda}$ vectors $\{|\mathbf{w}^{(\lambda,q')}\}$ are not necessarily orthogonal, but they are linearly independent.

Introducing the $(N-p_{\lambda}) \times N$ matrix $\mathbf{W}_{\lambda} \equiv (|\mathbf{w}^{(\lambda,(p_{\lambda}+1))}|, \dots, |\mathbf{w}^{(\lambda,N)})^T$ and its transpose $\mathbf{W}_{\lambda}^T$ we can construct the complementary projection operator $\hat{Q}_{\lambda} \equiv \hat{I} - \hat{P}_{\lambda}$ (see, e.g., Ref. [36]),

$$ \hat{Q}_{\lambda} = \mathbf{W}_{\lambda}^{T} (\mathbf{W}_{\lambda} \mathbf{W}_{\lambda}^{T})^{-1} \mathbf{W}_{\lambda}. \quad (A13) $$

Hence, based on Eq. (A13), Eq. (A11) transforms to

$$ \begin{aligned} b_{\lambda}^{2} &= (\mathbf{h}|\hat{I} - \hat{Q}_{\lambda}|\mathbf{h}) \\ &= 2 - (\mathbf{h}|\mathbf{W}_{\lambda}^{T}(\mathbf{W}_{\lambda}\mathbf{W}_{\lambda}^{T})^{-1}\mathbf{W}_{\lambda}|\mathbf{h}), \end{aligned} \quad (A14) $$

where we have used $(\mathbf{h}|\mathbf{h}) = 2$, see Eq. (A1).

## Appendix B: Derivation of the equations of the non-Markovian (NM) theory

Here we briefly describe the main steps that lead to the equations of the NM theory.

Subtracting the stationary probability distribution $p_s(\{\mathbf{a}\})$ from Eq. (6) in the main text, and integrating the result between $t=0$ and $t=\infty$, one obtains

$$ T p_s(\{\mathbf{a}\}) = \int_0^\infty dt \int d\{\mathbf{a}'\} \pi(\{\mathbf{a}'\}) [p(\{\mathbf{a}\}, t|\{\mathbf{a}'\}) - p(\{\mathbf{a}\}, t)], \quad (B1) $$

with the unknown quantities are the MFCT ($T$) and the (normalized) probability distribution $\pi(\{\mathbf{a}\})$ of the configurations at the instant of first contact. The probabilities $p(\{\mathbf{a}\}, t|\{\mathbf{a}'\})$, $p(\{\mathbf{a}\}, t)$ and $p_s(\{\mathbf{a}\})$ are Gaussian and admit analytical expressions. Consider now a unit vector $\hat{\mathbf{u}}$ and let $\pi_{\hat{\mathbf{u}}}$ be the distribution of modes at the first contact time $t^*$ with the constraint $\mathbf{R}(t^*) = a\hat{\mathbf{u}}$. In other words, $\pi_{\hat{\mathbf{u}}}(\{\mathbf{a}\}) \equiv C \pi(\mathbf{a})\delta(\sum_{\lambda'} a_{\lambda'}b_{\lambda'} - a\hat{\mathbf{u}})$, with $C$ a normalization constant.

In the NM theory, one assumes that $\pi_{\hat{\mathbf{u}}}$ is a multivariate Gaussian, with mean vector $m_{\lambda}\hat{\mathbf{u}}$, while the covariance matrix is approximated by that of the equilibrium configurations satisfying the constraint $\mathbf{R} = a\hat{\mathbf{u}}$. Then, the $\{m_{\lambda}\}$ are found by solving a set of equations that make the theory self-consistent. These equations are found by multiplying Eq. (B1) by $a_{\lambda}\delta(\sum_{\lambda'} b_{\lambda'}a_{\lambda'} - \mathbf{r}^*)$ for fixed $\lambda$ and $\mathbf{r}^*$ (satisfying $|\mathbf{r}^*| < a$) and integrating over all configurations. Such calculations involve the evaluation of a number of Gaussian integrals, as well as averages over angular directions $\hat{\mathbf{u}}$, and are done carefully in the appendix C of Ref. [15]. The main difference between

the calculations of Ref. [15] and ours is the fact that, here, the sums are carried out over the distinct values of the eigenvalues $\lambda$ instead of over all the indexes $i$. The resulting equation is

$$ \int_0^\infty \frac{dt}{\psi^{5/2}} \left\{ \left[ \frac{\mu_\lambda^{\pi,0} R_\pi}{3} + \frac{b_\lambda \phi(\phi - e^{-\frac{\lambda t}{\tau_0}})}{\lambda} \right] e^{-\frac{R_\pi^2}{2\psi}} - \frac{b_\lambda \phi(\phi - e^{-\frac{\lambda t}{\tau_0}})}{Z_0(a, nl^2/3)\lambda} \left[ Z_0(a, \psi) - \frac{G_0(a, \psi)}{3\psi} \right] \right\} = 0, \quad (B2) $$

where

$$ \mu_{\lambda}^{\pi,0} = m_{\lambda} e^{-\frac{\lambda t}{\tau_0}} - \frac{R_{\pi} b_{\lambda} (1 - \phi e^{-\frac{\lambda t}{\tau_0}})}{\lambda \psi} \quad (B3) $$

and

$$ Z_0(a, h) = \int_a^\infty dx \, x^2 e^{-\frac{x^2}{2h}}, \quad (B4) $$

$$ G_0(a, h) = \int_a^\infty dx \, x^4 e^{-\frac{x^2}{2h}}. \quad (B5) $$

By integration by parts, one obtains $Z_0(a, h) - G_0(a, h)/(3h) = -(a^3/3)e^{-a^2/(2h)}$ so that Eq. (B2) can be written in the compact form

$$ \int_0^\infty dt \left\{ \left[ m_\lambda e^{-\frac{\lambda t}{\tau_0}} - \frac{R_\pi b_\lambda (1 - \phi e^{-\frac{\lambda t}{\tau_0}})}{\lambda \psi} \right] \frac{R_\pi e^{-\frac{R_\pi^2}{2\psi}}}{3 \psi^{5/2}} + \frac{b_\lambda \phi(\phi - e^{-\frac{\lambda t}{\tau_0}})}{\lambda \psi^{5/2}} \left[ e^{-\frac{R_\pi^2}{2\psi}} - \frac{(a^3/3) e^{-a^2/2\psi}}{\int_a^\infty dr r^2 e^{-\frac{3r^2}{2nl^2}}} \right] \right\} = 0. \quad (B6) $$

There are as many equations in Eq. (B6) as the number of distinct nonvanishing eigenvalues $\lambda$, so that the values of $m_{\lambda}$ can be determined numerically, enabling the evaluation of the NM contact time with Eq. (8) of the main text. Given that for many fractal structures the eigenvalue spectra contain a lot of degenerate eigenvalues, this method allows to treat very large fractal structures, as shown in the main text.

## Appendix C: MFCT scaling for large fractal structures

Let us consider Eq. (7) of the main text,

$$ T = \int_0^\infty dt \left\{ \frac{e^{-a^2 \phi(t)^2 / [2\psi(t)]}}{[1 - \phi(t)^2]^{3/2}} - \frac{Z(a/\psi(t))^{1/2}}{Z(a\sqrt{3}/l\sqrt{n})} \right\}, \quad (C1) $$

where $Z(y) = \int_y^\infty dx \, x^2 e^{-x^2/2}$, in the case of large $n$. In our case the reactive monomers are at extremities such that $n^{1/2} \sim R_g/l$ ($R_g$ is the gyration radius of the whole structure). Thus, the largest time scale is given by the
largest relaxation time of the whole structure $\tau_N$. For this it is known [29, 31] that

$$ \tau_N \sim \tau_0 (R_g / l)^{4/(2-d_s)}. \quad (C2) $$

With $\gamma = 1 - d_s/2$ and with $n^{1/2} \sim R_g/l$ we get

$$ \tau_N \sim \tau_0 n^{1/\gamma}. \quad (C3) $$

If the monomers are separated (in physical space) by $n^{1/2} \sim R_g/l$, we can reasonably assume that

$$ \phi(t) = \Phi(t/\tau_N), \quad (C4) $$

where $\Phi$ is a dimensionless function independent on $N$, with $\Phi(\tau) \simeq 1$ at short times, and exponentially decaying for $\tau \gg 1$. Now, using that the function $Z(y)$ involved in Eq. (C1) goes for small $y$ as $Z(y) \simeq \sqrt{\pi/2} - y^3/3$, for a target of size $a$ which is much smaller than the typical distance $n$ between the reactants, Eq. (C1) reads

$$ T = \int_0^\infty dt \left\{ \frac{1}{[1 - \phi^2]^{3/2}} - 1 \right\}. \quad (C5) $$

Changing of variables in (C5) according to Eq. (C4), we obtain

$$ T = \tau_N \int_0^\infty d\tau \left\{ \frac{1}{[1 - \Phi^2]^{3/2}} - 1 \right\} \sim \tau_N \sim \tau_0 n^{1/\gamma} \quad (C6) $$

Note that the above integral converges because $\Phi(\tau)$ decays exponentially at large times, and behaves as $1 - \tau^\gamma$ for small times, with $\gamma < 2/3$.

Now, the expression for MFCT in non-Markovian case has a similar structure as Eq. (C1). The difference being in the first term, where $a^2\phi(t)^2$ in the exponent is replaced through a square of the reactive conformation $R_\pi(t) = \sum_\lambda b_\lambda m_\lambda e^{-\lambda t/\tau_0}$. Given that the analytic structure of the $m_\lambda$ involved is unknown ($m_\lambda$ are obtained from a set of self-consistent Eq. (B6)), we cannot copy steps of WF and hence numerical evaluations are needed. For this the iterative procedures described below are of much help.

## Appendix D: Calculation of the coefficients $b_\lambda$ by decimation of Vicsek fractals

In this subsection we show for Vicsek fractals how to construct the matrices $W_\lambda$ involved in Eq. (A14). We adapt a decimation procedure that was already used in Ref. [20] for the computation of the eigenvalues of the dynamical matrix. The topology of a Vicsek fractal is shown on Fig. 5, where we have used two symbols for the beads: circles and squares. It is clear that removing all the circles from the structure, one obtains another Vicsek fractal but of a lower generation, formed by squares in Fig. 5. Equivalently, one can say that the structure at the current generation is formed by the beads that were

<figure><img src="image_6.png" /><figcaption>FIG. 5. Vicsek fractal (VF) of functionality $f = 4$ and generation $g = 3$. The reactive beads are colored in red. The picture shows only the topology of the structure, particular VF conformations may come up in vastly different geometric forms.</figcaption></figure>

already part of the previous generation (the squares in Fig. 5, their indices are denoted by greek letters, say $\mu$) and by "new" beads (the circles on Fig. 5, their indices are denoted by latin letters, say by $k$).

Now, consider an eigenvector $\Phi^{(g)} = (\phi_1, ..., \phi_N)$ of the dynamical matrix of a Vicsek fractal of generation $g$ which is associated to the eigenvalue $\lambda^{(g)}$. It is shown in Ref. [20] that, if $\lambda^{(g)} \neq \{0, 1, f+1\}$, then necessarily the vector $\Phi^{(g-1)} = \{\phi_\mu\}$ is itself an eigenvector of a Vicsek fractal of the former generation ($g-1$), where $\mu$ represents the ensemble of indices of the beads present at the former generation. The eigenvalue $\lambda^{(g-1)}$ to which $\{\phi_\mu\}$ is associated is then

$$ \lambda^{(g-1)} = \lambda^{(g)}(\lambda^{(g)} - 3)(\lambda^{(g)} - f - 1). \quad (D1) $$

If the values of $\{\phi_\mu\}$ are known for the sites present at generation $g-1$, the values of $\phi_k$ for the "new" sites, are deduced from the equations $\sum_{i=1}^N (A_{ki} - \lambda^{(g)}\delta_{ki})\Phi_i^{(g)} = 0$ for all the new sites $k$. Hence, if one knows the matrix $W_{\lambda^{(g-1)}}$, then the matrix $W_{\lambda^{(g)}}$ at next generation is

$$ W_{\lambda^{(g)}}^{(g)} = \begin{pmatrix} W_{\lambda^{(g-1)}}^{(g-1)} & 0 \\ \text{Lines of } A^{(g)} - \lambda^{(g)}I & \text{for the "new" sites } k \end{pmatrix}. \quad (D2) $$

This relation allows the iterative construction of the matrices $W_\lambda$. Actually, since (D1) is a cubic equation, an eigenvalue at generation $g-1$ will generate three eigenvalues at generation $g$ and the 3 associated matrices $W_\lambda$.

The procedure needs to be initialized by expliciting the matrix $W_1^{(g)}$ associated to the eigenvalue $\lambda = 1$ at generation $g$. An eigenvector $\Phi = (\phi_1, ..., \phi_N)$ is associated
<figure><img src="image_7.png" /><figcaption>FIG. 6. Coefficients $b_{\lambda}^2$ for a Vicsek fractal of generation $g = 4$ and functionality $f = 4$. The reactive beads are at the extremities of the structure, as depicted on Fig. 5. The red crosses represent the results of a brute-force calculation (Eq. (A7)) while the blue circles are obtained with a decimation procedure (Eq. (A14)).</figcaption></figure>

to the eigenvalue $\lambda = 1$ if and only if $(\mathbf{A} - \mathbf{I})\mathbf{\Phi} = \mathbf{0}$. This leads to the equations for all the beads $\mu$ of functionality $f$:

$$ \phi_{\mu} = 0, \quad \sum_{i \in \text{NN of } \mu} \phi_i = 0, \qquad (\text{D3}) $$

where NN means nearest-neighbors. Also, for each pair of sites of functionality 2, we obtain the $(f+1)^{g-1}-1$ equations

$$ \phi_i - \phi_j = 0. \qquad (\text{D4}) $$

The equations (D3,D4) are linearly independent and determine the eigensubspace related to $\lambda = 1$. We can rewrite these equations under the matrix form $\mathbf{W}_1^{(g)}\mathbf{\Phi} = \mathbf{0}$, leading to the ready identification of the matrix $\mathbf{W}_1^{(g)}$.

We note, that we did not consider the eigenvalues that are generated iteratively from the eigenvalue $\lambda = f+1$, which are non-degenerate. These eigenvalues are associated [20, 37] to eigenvectors $\{\mathbf{\Phi}\}$ that all satisfy $\phi_{q_1} = \phi_{q_2}$, where $q_1$ and $q_2$ are the indices of the reactive beads (the red beads on Fig. 5), located at the extremities of the structure in this work. Hence, based on Eqs. (A2, A7) we see that the coefficients $b_{\lambda}^2$ vanish for these non-degenerate eigenvalues, which do not need to be considered anymore. We also note that the eigenvalue $\lambda = 0$, associated to the motion of center-of-mass, has no relevance in our problem concerned with intramolecular reactions.

Finally, the knowledge of all the matrices $\mathbf{W}_{\lambda}$ (which are very sparse) allows one to compute the coefficients $b_{\lambda}^2$ based on Eq. (A14). We checked our iterative procedure by comparing our results to the direct computation issued from brute force diagonalization and the use of Eq. (A7) (see Fig. 6)

## Appendix E: Calculation of the coefficients $b_{\lambda}$ for the dual Sierpinski gaskets by a decimation procedure.

We now describe the decimation procedure for the dual Sierpinski gasket. We represented a part of such fractals on Fig. 7a. As seen on this figure, the structure contains small triangles, and if one replaces each triangle by a new site, one obtains another dual Sierpinski gasket, at former generation (Fig. 7b). This remark is at the basis of the decimation procedure.

<figure><img src="image_8.png" /><figcaption>FIG. 7. (a) Topology of one part of a dual Sierpinski gasket. (b) decimated dual Sierpinski gasket.</figcaption></figure>

Let us consider an eigenvector $\mathbf{\Phi} = (\phi_1, \phi_2, ...)$ of the dynamical matrix $\mathbf{A}$ associated to the eigenvalue $\lambda$, which therefore satisfies

$$ \sum_{j=1}^{N} A_{ij} \phi_j = \lambda \phi_i. \qquad (\text{E1}) $$

Let us consider another vector $\mathbf{\Psi} = (\psi_{\alpha}, \psi_{\beta}, ...)$ (whose indices are the sites in the decimated structure), defined such as the coordinate $\psi_{\alpha'}$ is the sum of the values of $\phi_{i'}$ for the sites $i'$ in the triangle around $\alpha'$, and of the sites that are neighbors to these sites. For example, for the site $\alpha'$, we write (see the labeling on Fig.7a)

$$ \psi_{\alpha'} = \phi_{1'} + \phi_{2'} + \phi_{3'} + \phi_{4'} + \phi_{7'} + \phi_{5}, \qquad (\text{E2}) $$

For the sites at the boundary of the structure (i.e. of functionality 2), we modify the definition of $\psi$ by setting

$$ \psi_{\alpha} = 2\phi_{1} + \phi_{2} + \phi_{3} + \phi_{4} + \phi_{7}. \qquad (\text{E3}) $$

The relations (E1) (written for the sites $1', 2', ..., 9'$) and (E3,E2) (written for $\alpha, \beta, \gamma$) can be inverted to derive

$$ \phi_{i'} = \sum_{\mu = \{\alpha, \beta, \gamma\}} C_{i, \mu} \psi_{\mu'}, \qquad (\text{E4}) $$

where $C$ is a 9 × 3 matrix, of elements

$$ C_{1\alpha} = C_{5\beta} = C_{9\gamma} = L[-7 - \lambda(\lambda - 6)], \qquad (\text{E5}) $$

$$ C_{1\beta} = C_{1\gamma} = C_{5\alpha} = C_{5\gamma} = C_{9\alpha} = C_{9\beta} = L, \qquad (\text{E6}) $$

$$ C_{2\alpha} = C_{3\alpha} = C_{4\beta} = C_{6\beta} = C_{7\gamma} = C_{8\gamma} \\ = L(-4 + \lambda), \qquad (\text{E7}) $$

$$ C_{7\alpha} = C_{3\gamma} = C_{6\gamma} = C_{4\alpha} = C_{2\beta} = C_{8\beta} = -L. \qquad (\text{E8}) $$
All other elements of C are zero, and here L is the constant

$$L = [(\lambda - 5)(\lambda - 3)(\lambda - 2)]^{-1}. \qquad (E9)$$

Similarly, considering the 9 beads (1,2,...9) near the upper boundary of the structure, we obtain

$$\phi_i = \sum_{\mu = \{\alpha, \beta, \gamma\}} C_{i, \mu} \psi_{\mu}, \qquad (E10)$$

where the matrix C is the same as in (E4). Thus, we have obtained an inverse relation for the coordinates $\phi_i$ as a function of the vector defined on the decimated structure $\psi_\alpha$, which holds as soon as $\lambda \neq \{2, 3, 5\}$ (otherwise L is infinite in (E9)). Furthermore, we obtain the relations

$$[3 - (5 - \lambda)\lambda]\psi_{\alpha'} = \psi_{\beta'} + \psi_{\gamma'} + \psi_{\beta} \qquad (E11)$$

and similar relations hold for $\psi_\beta, \psi_\gamma$, while for the site at the boundary,

$$[2 - (5 - \lambda)\lambda]\psi_{\alpha} = \psi_{\beta} + \psi_{\gamma}. \qquad (E12)$$

Equations (E11,E12) mean that $\Psi$ is an eigenvector of the decimated structure, with the eigenvalue

$$\tilde{\lambda} = (5 - \lambda)\lambda. \qquad (E13)$$

Equation (E13) is an iterative relation between the eigenvalues at one generation ($\tilde{\lambda}$) and the eigenvalues at larger generation ($\lambda$), and is already known [32, 38]. Thus, an eigenvector $\Phi$ associated to an eigenvalue $\lambda \neq \{2, 3, 5\}$ can be expressed as a product $C_\lambda \Psi$, where $C_\lambda$ is the matrix defined in (E5-E8), and $\Psi$ is an eigenvector of the decimated structure (at generation $g-1$) associated to the eigenvalue $\lambda^{(g-1)} = (5-\lambda)\lambda$. Let us iterate the procedure and consider $s$ decimation steps, until one reaches an eigenvalue $\lambda^{(g-s)} = \{3, 5\}$: then an eigenvector $\Phi$ is expressed as $\Phi = B\Psi^{(g-s)}$, where $\Psi^{(g-s)}$ is an eigenvector at generation $g-s$, and $B = \Pi_{s'=0}^{s-1} C_\lambda^{(g-s')}$. Let $W_3$ be the matrix of vectors which span the subspace orthogonal to all the eigenvectors related to the eigenvalue $\lambda = 3$ at generation $g-s$. Then, in the basis related to the eigenvalue $\lambda = 3$ at generation $g-s$ Eq.(A14) transforms to

$$b_{\lambda}^{2} = (\mathbf{h}|\mathbf{B}(\mathbf{B}^{T}\mathbf{B})^{-1}\mathbf{B}^{T}|\mathbf{h}) - (\mathbf{h}|\mathbf{B}\mathbf{W}_{3}^{T}(\mathbf{W}_{3}\mathbf{B}^{T}\mathbf{B}\mathbf{W}_{3}^{T})^{-1}\mathbf{W}_{3}\mathbf{B}^{T}|\mathbf{h}). \qquad (E14)$$

The eigenvalue $\lambda = 3$ plays a fundamental role in the whole decimation procedure. The eigenvectors related to

this eigenvalue are the so-called symmetric modes [32]. This means that for each eigenvector $\Phi$ associated to $\lambda = 3$, the following relations hold [32]:

$$\phi_i + \phi_j + \phi_k = 0, \qquad (E15)$$

for each triplet $(i, j, k)$ of nearest-neighboring beads forming the smallest triangles of the structure (see Fig. 7). Furthermore, we have also [32]

$$\phi_n = \phi_m, \qquad (E16)$$

<figure><img src="image_9.png" /><figcaption>FIG. 8. Coefficients $b_\lambda^2$ for a dual Sierpiński gasket of generation $g=8$. The reactive beads are at extremities, as depicted on Fig.1 in the main text. Red crosses represent the results of the brute-force calculation (Eq. (A7)) and blue circles of the decimation procedure (Eq. (E14)).</figcaption></figure>

for each pair $(n, m)$ of nearest-neighboring beads connecting the smallest triangles. Equations (E15,E16) are linearly independent and determine the eigensubspace of $A$ associated to 3, they can therefore be written in matrix form $W_3\Phi = 0$, from which the (sparse) matrix $W_3$ reads.

On the contrary, the eigenvalue $\lambda = 5$, associated to antisymmetric modes [32], does not play any role in our problem: when the reactive beads are located at the extremities of the structure (Fig. 1b in the main text), the coefficients $b_\lambda^2$ vanish for symmetry reasons for $\lambda = 5$ and for any eigenvalue generated from the eigenvalue 5.

We implemented a recursive algorithm which calculates the coefficients $b_\lambda$ with our decimation procedure. The results were checked by comparing them to the output of the brute force diagonalization of the dynamical matrix, one sees on Fig. 8 that our procedure gives the correct results. We used our decimation approach to treat much larger structures and to generate Fig. 3 of the main text.

[1] G. Bonnet, O. Krichevsky, and A. Libchaber, Proc. Natl. Acad. Sci. USA **95**, 8602 (1998); X. Wang and W. M. Nau, J. Am. Chem. Soc. **126**, 808 (2004).

[2] M. I. Wallace, L. Ying, S. Balasubramanian, and D. Klenerman, Proc. Natl. Acad. Sci. USA **98**, 5584 (2001); J.-F. Allemand, S. Cocco, N. Douarche, and
G. Lia, Eur. Phys. J. E **19**, 293 (2006).

[3] L. J. Lapidus, W. A. Eaton, and J. Hofrichter, Proc. Natl. Acad. Sci. USA **97**, 7220 (2000); A. Möglich, K. Joder, and T. Kiefhaber, *ibid.* **103**, 12394 (2006).

[4] J. K. Gooden, M. L. Gross, A. Mueller, A. D. Stefanescu, and K. L. Wooley, J. Am. Chem. Soc. **120**, 10180 (1998); Y. Zheng, H. Cao, B. Newland, Y. Dong, A. Pandit, and W. Wang, *ibid.* **133**, 13130 (2011); A. Burgath, A. Sunder, and H. Frey, Macromol. Chem. Phys. **201**, 782 (2000).

[5] M. Doi and S. F. Edwards, *The Theory of Polymer Dynamics* (Clarendon Press, 1988).

[6] A. A. Gurtovenko and A. Blumen, Adv. Polym. Sci. **182**, 171 (2005).

[7] P.-G. De Gennes, J. Chem. Phys. **76**, 3316 (1982).

[8] G. Wilemski and M. Fixman, J. Chem. Phys. **60**, 866 (1974); **60**, 878 (1974).

[9] A. Szabo, K. Schulten, and Z. Schulten, J. Chem. Phys. **72**, 4350 (1980).

[10] B. Friedman and B. O'Shaughnessy, Macromolecules **26**, 5726 (1993).

[11] I. M. Sokolov, Phys. Rev. Lett. **90**, 080601 (2003).

[12] T. Guérin, O. Bénichou, and R. Voituriez, Nat. Chem. **4**, 568 (2012).

[13] A. Amitai, I. Kupka, and D. Holcman, Phys. Rev. Lett. **109**, 108302 (2012).

[14] R. Afra and B. A. Todd, J. Chem. Phys. **138**, 174908 (2013).

[15] T. Guérin, O. Bénichou, and R. Voituriez, J. Chem. Phys. **138**, 094908 (2013); Phys. Rev. E **87**, 032601 (2013); T. Guérin, M. Dolgushev, O. Bénichou, R. Voituriez, and A. Blumen, **90**, 052601 (2014); O. Bénichou, T. Guérin, and R. Voituriez, J. Phys. A **48**, 163001 (2015).

[16] J. Shin, A. G. Cherstvy, and R. Metzler, *Soft Matter* **11**, 472 (2015); ACS Macro Lett. **4**, 202 (2015).

[17] B. I. Voit and A. Lederer, Chem. Rev. **109**, 5924 (2009).

[18] W. Burchard, Adv. Polym. Sci. **143**, 113 (1999).

[19] I. M. Sokolov, J. Klafter, and A. Blumen, Phys. Today **55**, 48 (2002).

[20] A. Blumen, C. von Ferber, A. Jurjiu, and T. Koslowski, Macromolecules **37**, 638 (2004).

[21] S. Reuveni, R. Granek, and J. Klafter, Phys. Rev. Lett. **100**, 208101 (2008); Proc. Natl. Acad. Sci. USA **107**, 13696 (2010); S. Reuveni, J. Klafter, and R. Granek, Phys. Rev. Lett. **108**, 068101 (2012).

[22] R. Burioni, D. Cassi, F. Cecconi, and A. Vulpiani, Proteins **55**, 529 (2004); A. Banerji and I. Ghosh, Cell. Mol. Life Sci. **68**, 2711 (2011).

[23] A. G. Zilman and R. Granek, Phys. Rev. E **58**, R2725 (1998).

[24] C. M. Sorensen, Aerosol Sci. Technol. **35**, 648 (2001).

[25] T. Nakayama, K. Yakubo, and R. L. Orbach, Rev. Mod. Phys. **66**, 381 (1994).

[26] S. Alexander and R. Orbach, J. Phys. Lett. **43**, 625 (1982).

[27] D. Panja, J. Stat. Mech. **06**, 06011 (2010); S. Gupta, A. Rosso, and C. Texier, Phys. Rev. Lett. **111**, 210601 (2013).

[28] See Supplemental Material at [URL will be inserted by publisher] for the identification of $b_{\lambda}^2$ and the iterative calculation of correlation functions, as well as for the details of the MFCT theories, which includes Refs.[15, 20, 29, 31, 32, 36–38].

[29] J.-U. Sommer and A. Blumen, J. Phys. A: Math. Gen. **28**, 6669 (1995).

[30] S. Condamín, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, Nature **450**, 77 (2007).

[31] S. Reuveni, J. Klafter, and R. Granek, Phys. Rev. E **85**, 011906 (2012).

[32] M. G. Cosenza and R. Kapral, Phys. Rev. A **46**, 1850 (1992).

[33] E. Agliari, Phys. Rev. E **77**, 011128 (2008).

[34] B. Meyer, E. Agliari, O. Bénichou, and R. Voituriez, Phys. Rev. E **85**, 026113 (2012).

[35] R. W. Pastor and R. Zwanzig, A. Szabo, J. Chem. Phys. **105**, 3878 (1996).

[36] C. D. Meyer, *Matrix analysis and applied linear algebra* (Siam, 2000).

[37] F. Fürstenberg, M. Dolgushev, and A. Blumen, J. Chem. Phys. **138**, 034904 (2013).

[38] R. Rammal and G. Toulouse, J. Phys. Lett. **44**, 13 (1983).