# Non-Markovian closure kinetics of flexible polymers with hydrodynamic interactions

N. Levernier, M. Dolgushev, O. Bénichou, A. Blumen, T. Guérin, R. Voituriez

► To cite this version:

N. Levernier, M. Dolgushev, O. Bénichou, A. Blumen, T. Guérin, et al. Non-Markovian closure kinetics of flexible polymers with hydrodynamic interactions. The Journal of Chemical Physics, 2015, 143 (20), pp.204108 (1-10). 10.1063/1.4935966. hal-01256795

HAL Id: hal-01256795

https://hal.science/hal-01256795

Submitted on 15 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>
# Non-Markovian closure kinetics of flexible polymers with hydrodynamic interactions

N. Levernier,<fnref n="1" /> M. Dolgushev,<fnref n="2" /> O. Bénichou,<fnref n="1" /> A. Blumen,<fnref n="2" /> T. Guérin,<fnref n="3" /> and R. Voituriez<fnref n="1" />,<fnref n="4" />

This paper presents a theoretical analysis of the closure kinetics of a polymer with hydrodynamic interactions. This analysis, which takes into account the non-Markovian dynamics of the end-to-end vector and relies on the preaveraging of the mobility tensor (Zimm dynamics), is shown to reproduce very accurately the results of numerical simulations of the complete nonlinear dynamics. It is found that Markovian treatments based on a Wilemski-Fixman approximation significantly overestimate cyclization times (up to a factor 2), showing the importance of memory effects in the dynamics. In addition, this analysis provides scaling laws of the mean first cyclization time (MFCT) with the polymer size *N* and capture radius *b*, which are identical in both Markovian and non-Markovian approaches. In particular, it is found that the scaling of the MFCT for large *N* is given by $T \sim N^{3/2} \ln(N/b^2)$, which differs from the case of the Rouse dynamics where $T \sim N^2$. The extension to the case of the reaction kinetics of a monomer of a Zimm polymer with an external target in a confined volume is also presented.

## I. INTRODUCTION

Reactions involving polymers are ubiquitous in nature. Among them, reactions of closure of linear chains are of particular interest since they are involved in a number of chemical and biological processes. Examples cover gene regulation by the formation of RNA-hairpins<sup>3,4</sup> or DNA-loops,<sup>5,6</sup> the folding of polypeptides,<sup>7,8</sup> as well as the appearance of cycles in synthetic polymers.<sup>9–11</sup> In the diffusion-controlled regime, the kinetics of contact formation strongly depends on the complex dynamics of the reactive monomers. As a result of the collective dynamics of all the monomers in the chain, the motion of a single monomer is often subdiffusive<sup>12</sup> and presents non-Markovian features<sup>13</sup> (i.e., memory effects), which lead to nontrivial reaction kinetics<sup>14</sup> even for the simplest models of polymers.

On the theoretical level, various approaches have been proposed to quantify the kinetics of polymer closure.<sup>1,2,15–26</sup> An important step has been provided by Wilemski and Fixman,<sup>15,16</sup> who made a local equilibrium assumption, thereby replacing the non-Markovian problem by an effective Markovian approach. Other theories include the so-called Szabo, Schulten, and Schulten (SSS) analytical approach,<sup>17</sup> which neglects important aspects of the polymer dynamics and the perturbative renormalization group theory,<sup>22</sup> which provides results at leading order in the parameter $\epsilon = 4 - d$, with *d* the spatial dimension. Recent works have improved these approaches in the case of the Rouse chain (i.e., flexible without hydrodynamic interactions) by introducing several methods: (i) a refined way to take into account the memory of the initial configurations,<sup>19</sup> (ii) an exact formal iterative

resolution scheme in one dimension,<sup>20</sup> (iii) a strong localized perturbation analysis,<sup>25</sup> and (iv) an approach based on the calculation of the distribution of chain configurations at the instant of cyclization,<sup>21,27–29</sup> which has a strong influence on the contact kinetics.

The Rouse model, however, provides an incorrect description of polymer chain dynamics in dilute solutions, where hydrodynamic interactions are long-ranged and deeply modify the motion of the polymer<sup>12,30</sup> as well as the closure kinetics.<sup>2,22</sup> Since hydrodynamic interactions are non-linear, most of available analytical treatments involve approximations. A standard approach consists in using a pre-averaged form introduced by Zimm,<sup>30</sup> which is known to be accurate.<sup>12</sup> In the context of closure kinetics, it was found in Ref. 1 by simulations that the Wilemski-Fixman treatment of the cyclization kinetics of Zimm chains systematically overestimates the reaction times, in a way that cannot be attributed to the simplified treatment of these interactions. This suggests that non-Markovian effects, which have not been described before in the presence of hydrodynamic interactions, play an important role in the cyclization kinetics.

This paper is devoted to the theoretical description of the cyclization kinetics of chains with hydrodynamic interactions in the diffusion controlled regime. This problem has already been discussed in the Wilemski-Fixman approach,<sup>1,2</sup> and the perturbative renormalization group theory.<sup>22</sup> However, a non-Markovian description of the closure kinetics has not been considered so far. The first goal consists in describing these effects by adapting a theory proposed recently,<sup>21,27,31</sup> which was so far restricted to the case of flexible and semi-flexible free draining chains. One could expect that the treatment

<footnote n="1">¹Laboratoire de Physique Théorique de la Matière Condensée, CNRS/UPMC, 4 Place Jussieu,
Paris 75005, France</footnote>

<footnote n="2">²Theoretical Polymer Physics, University of Freiburg, Hermann-Herder-St. 3, D-79104 Freiburg, Germany</footnote>

<footnote n="3">³Laboratoire Ondes et Matière d'Aquitaine, University of Bordeaux, Unité Mixte de Recherche 5798, CNRS,
F-33400 Talence, France</footnote>

<footnote n="4">⁴Laboratoire Jean Perrin, CNRS/UPMC, 4 Place Jussieu, Paris 75005, France</footnote>
of pre-averaged interactions is not valid for the problem of chain closure, since interactions averaged over equilibrium configurations could differ from the actual (non-preaveraged) interactions for closed chains, which could then play an important role for the contact kinetics. We will show however that this effect is relatively small: indeed, calculating the average shape of the polymer at the instant of cyclization for pre-averaged interactions yields a reaction time that is in quantitative agreement with the results of Brownian dynamic simulations performed without pre-averaging. The second goal of the paper is to derive new scaling laws for the mean cyclization time in the limit of small capture radius and long chains. In particular, we prove that for long chains, the mean cyclization time $T$ scales as $T \sim N^{3/2} \log(N/b^2)$, with $N$ the number of monomers and $b$ the capture radius, for both the Wilemski-Fixman approach and the non-Markovian theory, with however a different prefactor. In addition, we predict scaling laws for intermolecular reactions, by extending the formalism to the case of the search of an external target of size $b$ in a confined volume $V$. We show that in this case, the mean reaction time scales as $T \sim V \log(b/\sqrt{N})$ for long chains.

The outline of the paper is as follows. In Section II, we recall briefly some important results of the Zimm model. In Section III, we show how to adapt the recent non-Markovian theory<sup>21,27,31</sup> of cyclization kinetics to the Zimm model. We compare the numerical solution of the model to both Brownian dynamics simulations and to the results of a Wilemski-Fixman approach. Section IV is devoted to the derivation of scaling laws for the mean cyclization time with the chain length and the capture radius, for both the Markovian and the non-Markovian theories. Finally, we present in Section V the direct application of our calculations to the case of the reaction of one monomer with an external target in confinement and derive explicit asymptotic formulas for the reaction time for this problem.

## II. POLYMER DYNAMICS WITH HYDRODYNAMIC INTERACTIONS

We consider the dynamics in a three-dimensional (3d) space of a polymer represented by $N$ beads at positions $(\mathbf{x}_1, \dots, \mathbf{x}_N)$ (in this paper, quantities in bold represent vectors or tensors in the 3d space). The beads are linked by springs of stiffness $k$, so that the force $\mathbf{F}_j$ exerted on a bead $j$ by the neighboring beads is

$$ \mathbf{F}_j = -k(2\mathbf{x}_j - \mathbf{x}_{j+1} - \mathbf{x}_{j-1}) = -\sum_{k=1}^{N} M_{jk}\mathbf{x}_k, \quad (1) $$

where we used the convention $\mathbf{x}_0 = \mathbf{x}_1$ and $\mathbf{x}_{N+1} = \mathbf{x}_N$ and $M$ is proportional to the (Laplacian) tridiagonal connectivity matrix. This force $\mathbf{F}_j$ is balanced by the force due to the fluid solvent, and from Newton's third law, a force $\mathbf{F}_j$ is exerted on the fluid at the position of bead $j$. This force generates a velocity field in the fluid that influences the motion of all other beads. We introduce a non-isotropic mobility tensor $\mathbf{D}_{ij}$ to describe these hydrodynamic interactions, such that the average velocity of the fluid at the position of the bead $i$ is $\mathbf{D}_{ij}\mathbf{F}_j$. The dynamics in the presence of hydrodynamic interactions then follows the Langevin equation

in the overdamped limit

$$ \dot{\mathbf{x}}_i = \sum_{j=1}^{N} \mathbf{D}_{ij} \cdot \mathbf{F}_j + \zeta_i(t), \quad (2) $$

where $\zeta_i(t)$ is a stochastic Gaussian white noise term whose amplitude follows from the fluctuation-dissipation relation, $\langle \zeta_i(t) \otimes \zeta_j(t') \rangle = 2k_B T \mathbf{D}_{ij} \delta(t-t')$, with $k_B T$ the thermal energy. We also introduce $l_0 = \sqrt{k_B T/k}$ the typical bond length and $\tau_0 = 6\pi\eta a/k$ the single bond characteristic relaxation time, where $\eta$ is the fluid viscosity and $a$ the monomer radius.

Different choices of mobility matrix $\mathbf{D}_{ij}$ exist. The simplest one is the Oseen tensor, describing the fluid velocity induced by a punctual force in the fluid, and therefore valid for quasi-punctual monomers. However, for very small distances between beads, the Oseen tensor may not be positive definite, making dynamics (2) unphysical due to the appearance of negative relaxation times.<sup>32</sup> Instead, for numerical purposes, we use the so-called Rotne-Prager tensor, which does not have this problem. We denote by $\mathbf{r}_{ij} = \mathbf{x}_j - \mathbf{x}_i$ the vector from the particle $i$ to the particle $j$ and $r_{ij} = |\mathbf{r}_{ij}|$. The elements of the Rotne-Prager tensor are then given as follows:<sup>33</sup> when $i \neq j$, and when the beads $i, j$ do not overlap ($r_{ij} > 2a$),

$$ \mathbf{D}_{ij} = \frac{1}{8\pi\eta r_{ij}} \left[ \left( 1 + \frac{2a^2}{3r_{ij}^2} \right) \mathbf{I} + \frac{\mathbf{r}_{ij} \otimes \mathbf{r}_{ij}}{r_{ij}^2} \left( 1 - \frac{2a^2}{r_{ij}^2} \right) \right], \quad (3) $$

while for $i \neq j$ and $r_{ij} < 2a$,

$$ \mathbf{D}_{ij} = \frac{1}{6\pi\eta a} \left[ \left( 1 - \frac{9}{32} \frac{r_{ij}}{a} \right) \mathbf{I} + \frac{3}{32} \frac{\mathbf{r}_{ij} \otimes \mathbf{r}_{ij}}{r_{ij} a} \right], \quad (4) $$

and finally, when $i = j$,

$$ \mathbf{D}_{ii} = \frac{1}{6\pi\eta a} \mathbf{I}, \quad (5) $$

where $\mathbf{I}$ is the $3 \times 3$ identity tensor. Eqs. (1)–(5) define the polymer dynamics that will be analyzed in this paper. For numerical purposes, we will make use of the choice $a = 0.25 l_0$.<sup>34,35</sup>

The non-linear dependence of the mobility matrix $\mathbf{D}_{ij}$ on the positions $\mathbf{x}_j$ makes the Langevin equation (2) very difficult to solve. In particular it does not admit Gaussian solutions. In his pioneering work,<sup>30</sup> Zimm overcame this difficulty by replacing the mobility matrix $\mathbf{D}_{ij}$ by its average value over the equilibrium distribution, thereby making the equation linear. It has been shown that this approximation catches the main physics of hydrodynamic interactions.<sup>12,30</sup> In the model introduced above, it would be therefore relevant to use a pre-averaged version of the Rotne-Prager tensor rather than the Oseen tensor. However, it is shown in Appendix A that these tensors differ at most by a few percents for the parameters that we use. Hence, it is sufficient to use for analytical calculations the pre-averaged form of the Oseen tensor, given by<sup>12,30</sup>

$$ \tilde{\mathbf{D}}_{ij} \approx \frac{1}{6\pi\eta} \left[ \frac{\delta_{ij}}{a} + \frac{(1 - \delta_{ij})}{l_0} \sqrt{\frac{2}{\pi|i - j|}} \right] \mathbf{I}. \quad (6) $$

Note that the pre-averaged Oseen tensor is isotropic.

The effects of hydrodynamic interactions on polymer dynamics have already been studied in details. In particular,
it is well-known that any monomer of a long polymer chain performs a subdiffusive motion at intermediate time scales, which results from the collective dynamics of all the monomers. More precisely, for $N \gg 1$, we have three different regimes for the mean-square displacement of a monomer,³⁶

$$ \langle [\mathbf{r}_i(t) - \mathbf{r}_i(0)]^2 \rangle \sim \begin{cases} D_G t & \text{for } t \ll \tau_0 \\ At^{2/3} & \text{for } \tau_0 \ll t \ll \tau_0 N^{3/2}, \\ D_G t/\sqrt{N} & \text{for } \tau_0 N^{3/2} \ll t \end{cases} \quad (7) $$

where $D_G = k_B T/\eta l_0$, $A = (k_B T/\eta)^{2/3}$ (see Eq. (24)) and all the numerical constants have been omitted. We can already note that in the subdiffusion regime, the dimension of the random walk of a monomer (defined such that $\langle [r(t) - r(0)]^2 \rangle \sim t^{2/d_w}$) is given by $d_w = 3$. The process in 3d is therefore marginally recurrent, whereas it is transient for both short and long time scales (see Ref. 37). This will play a key role in the determination of the scaling laws that we present below.

The goal of this paper is to characterize the Mean First Cyclization Time (MFCT, denoted $T$), defined as the average first time for which the distance $|\mathbf{x}_N - \mathbf{x}_1|$ between the two chain ends becomes smaller than a capture radius $b$. Initially the polymer is in an equilibrium configuration, with the constraint that $|\mathbf{x}_N - \mathbf{x}_1| > b$. The effect of the complex monomer dynamics resulting from the hydrodynamic interactions will be analyzed with both numerical simulation methods (which we describe in Appendix D) and analytical means (see Section III).

### III. THEORIES OF CYCLIZATION KINETICS

We present here a theoretical approach that enables the determination of the mean cyclization time $T$, which can be seen as the mean first passage time of a non-Markovian problem. The approach consists in adapting a recent theory, developed so far for free-draining chains,²⁷,²⁸,³¹ in which the MFCT is expressed in terms of the conformational distribution of chains at the very instant of cyclization, whose moments are computed by solving a set of self-consistent equations.

The main steps of this approach can be summarized as follows. We introduce the joint probability density $f(\{\mathbf{x}\}, t)$ that the 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 positions $\{\mathbf{x}\} = (\mathbf{x}_1, \mathbf{x}_2, \dots)$. We partition the trajectories that lead to a configuration $\{\mathbf{x}\}$ (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{x}\}$ in a time $t - t'$. The mathematical formulation of this decomposition of events is

$$ P(\{\mathbf{x}\}, t|\text{ini}, 0) = \int_0^t dt' \int d\{\mathbf{x}'\} f(\{\mathbf{x}'\}, t') P(\{\mathbf{x}\}, t - t'|\{\mathbf{x}'\}), \quad (8) $$

where $d\{\mathbf{x}\} = d\mathbf{x}_1 d\mathbf{x}_2 \dots d\mathbf{x}_N$, $P(\{\mathbf{x}\}, t|\{\mathbf{x}'\})$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from $\{\mathbf{x}'\}$ at $t = 0$ while $P(\{\mathbf{x}\}, t|\text{ini}, 0)$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from the initial conditions

at $t = 0$ (in which the chain is at equilibrium, with the condition that the reactive monomers are not necessarily in contact). Next, taking the Laplace transform of (8) and expanding for small values of the Laplace variable, we obtain²⁷

$$ T P_{\text{stat}}(\{\mathbf{x}\}) = \int_{0}^{\infty} dt \left[ \int d\Omega P(\{\mathbf{x}\}, t|\pi_{\Omega}) - P(\{\mathbf{x}\}, t|\text{ini}, 0) \right]. \quad (9) $$

In this equation, $\Omega$ represents the angular directions parametrized by $\theta$ and $\varphi$ in spherical coordinates, with the normalization $\int d\Omega = 1$; $\pi_{\Omega}(\{\mathbf{x}\})$ represents the probability distribution of configurations at the very instant of cyclization, given that the angular direction of the end-to-end vector at first contact is $\Omega$, and $P(\{\mathbf{x}\}, t|\pi_{\Omega})$ is the probability of $\{\mathbf{x}\}$ at $t$ starting from the distribution $\pi_{\Omega}$ at initial time. Equation (9) is exact and does not depend on the particular hypotheses of chain dynamics. It is derived in details in Ref. 27. However, this equation cannot be solved explicitly to the best of our knowledge and approximations have to be introduced.

The first simplifying step is to approximate the dynamics by a Gaussian dynamics in order to be able to evaluate the propagators appearing in (9); hence, we consider only the Zimm dynamics with the pre-averaged mobility tensor. Next, the simplest approach is to make a Markovian approximation, which consists in neglecting any memory effect by assuming that the distribution $\pi_{\Omega}$ is the equilibrium distribution conditional to $\mathbf{x}_N - \mathbf{x}_1 = b\hat{\mathbf{u}}_r(\Omega)$ (with $\hat{\mathbf{u}}_r(\Omega)$ the unit vector pointing in the direction $\Omega$). This corresponds to the so-called Wilemski-Fixman (WF) approximation.¹⁵,¹⁸,²⁷ Introducing this approximation into Eq. (8), integrating over all configurations, and taking the long time limit lead to the estimate $T_{\text{WF}}$ of the MFCT,²⁷

$$ T_{\text{WF}} = \int_0^\infty \frac{dt}{[1 - \phi(t)^2]^{3/2}} \left\{ e^{-b^2 \phi(t)^2 / [2\psi(t)]} - \frac{Z(b, \psi(t))}{Z(b, L^2)} \right\}, \quad (10) $$

where $\phi(t)$, $L^2$, $\psi(t)$ characterize the dynamics of the end-to-end vector $\mathbf{r}_{ee} = \mathbf{x}_N - \mathbf{x}_1$, which is assumed to be Gaussian, and $Z(b, x) = \int_b^\infty dR_0 e^{-R_0^2/2x}$. The function $\phi$ is the normalized temporal auto-correlation function of any coordinate of the end-to-end vector $\mathbf{r}_{ee}$,

$$ \phi(\tau) = \frac{\langle x_{ee}(t + \tau)x_{ee}(t) \rangle}{\langle x_{ee}(t)x_{ee}(t) \rangle}, \quad (11) $$

where the spatial coordinates of $\mathbf{r}_{ee}$ are $(x_{ee}, y_{ee}, z_{ee})$, $L^2$ is the equilibrium mean-square length

$$ L^2 = \langle x_{ee}(t)^2 \rangle = (N - 1)t_0^2, \quad (12) $$

and $\psi(t)$ is the mean-square-displacement of any coordinate of $\mathbf{r}_{ee}$ when the initial value of $\mathbf{r}_{ee}$ is fixed and is easily shown to be related to $\phi$ by

$$ \psi(t) = \text{Var}(x_{ee}(t)|x_{ee}(0) = x_{ee}^0) = L^2[1 - \phi(t)^2], \quad (13) $$

where we denote $\text{Var}(y|B)$ the variance of the variable $y$ given the event $B$ is realized. Note that the original formula of Wilemski and Fixman is expressed in terms of the sink correlation function $C(t) = \iint d\mathbf{r} d\mathbf{r}_0 S(\mathbf{r}) S(\mathbf{r}_0) P(\mathbf{r}, t; \mathbf{r}_0, 0)$, where the sink function $S(\mathbf{r})$ is proportional to the absorption
probability per unit time for chains with end to end vector **r**. Eq. (10) can be recovered by setting the first sink function to $S(\mathbf{r}) = \delta(\mathbf{r})$ and the second one to $S(\mathbf{r}_0) = \delta(b - r_0)$, and by additionally taking into account the finite probability to start inside the target. It therefore corresponds to an imbalanced delta sink approximation (see Refs. 18, 27, and 29).

Going beyond the Wilemski-Fixman approximation requires a more precise description of the distribution $\pi_{\Omega}$. Here, our key hypothesis is to assume that $\pi_{\Omega}$ is a multivariate Gaussian distribution. For symmetry reasons, the average of $\mathbf{x}_i$ over $\pi_{\Omega}$ is along $\hat{\mathbf{u}}_r(\Omega)$. We denote this quantity $E_{\pi_{\Omega}}(\mathbf{x}_i) = m_i^{\pi} \hat{\mathbf{u}}_r(\Omega)$. We make the additional assumption that the covariance matrix of $\pi_{\Omega}$ is the same covariance matrix that characterizes the equilibrium distribution of chains with the constraint $\mathbf{r}_{ee} = b \hat{\mathbf{u}}_r(\Omega)$. Such an approximation is not necessary but largely simplifies the calculations. Previous studies on free-draining chains$^{21,28}$ have revealed that releasing this “stationary covariance approximation” only slightly improves the estimate of the cyclization time at an important calculation cost. A set of self-consistent equations that defines the $m_i^{\pi}$’s is found by multiplying (9) by $\mathbf{x}_i \delta(\mathbf{x}_N - \mathbf{x}_1 - \mathbf{r}^*)$ (for a fixed $\mathbf{r}^*$ satisfying $|\mathbf{r}^*| < b$) and integrating over all configurations. Adapting existing calculations for free-draining chains,$^{27}$ the resulting equation is (see the complete derivation in Appendix C)

$$ \int_0^\infty \frac{dt}{\psi^{3/2}} \left\{ e^{-\frac{R_\pi^2}{2\psi}} \left[ \frac{R_\pi}{3\psi} \left( \mu_i^\pi - \frac{\beta_i}{\psi} R_\pi \right) + \frac{\beta_i}{\psi} - \frac{\alpha_i}{L^2} \right] + \left( \frac{\beta_i}{\psi} - \frac{\alpha_i}{L^2} \right) \frac{b^3 e^{-b^2/2\psi}}{3Z(b, L^2)} \right\} = 0. \quad (14) $$

Here, we have denoted by $\mu_i^\pi(t)$ the average of $\mathbf{x}_i$ at time $t$ after the instant of cyclization in the direction $\hat{\mathbf{u}}(\Omega)$, $R_\pi = \mu_N^\pi - \mu_1^\pi$, $\beta_i(t)$ is the covariance between $x_i(t)$ and $x_{ee}(t)$ when $x_{ee}(0)$ is fixed, while $\alpha_i$ is the covariance between $x_i$ and $x_{ee}$ at equilibrium. In other words, $\alpha_i$ and $\beta_i$ characterize the dynamics and the equilibrium of the motion starting from a constrained equilibrium state,

$$ \beta_i(t) = \mathrm{Cov}(x_i(t), x_{ee}(t) | x_{ee}(0) = x_{ee}^0), \quad (15) $$

$$ \alpha_i = \mathrm{Cov}_{\mathrm{stat}}(x_i, x_{ee}) \quad (16) $$

(note that $\beta_i$ does not depend on $x_{ee}^0$), and $\mu_i^\pi$, $R_\pi$ characterize the motion of the chain in the future of the first contact,

$$ \mu_i^\pi(t) = \left\langle \mathbf{x}_i(t + t^*) \cdot \frac{\mathbf{r}_{ee}(t^*)}{|\mathbf{r}_{ee}(t^*)|} \right\rangle, \quad (17) $$

$$ R_{\pi}(t) = \left\langle \mathbf{r}_{ee}(t + t^*) \cdot \frac{\mathbf{r}_{ee}(t^*)}{|\mathbf{r}_{ee}(t^*)|} \right\rangle = \mu_{N}^{\pi}(t) - \mu_{1}^{\pi}(t), \quad (18) $$

where $t^*$ is the first cyclization time. The time evolution of $\mu_i^\pi$ follows from the Langevin equation (2),

$$ \partial_t \mu_i^\pi(t) = - \sum_{j,k=1}^{N} D_{ij} M_{jk} \mu_k^\pi(t), \quad \mu_i^\pi(0) = m_i^\pi. \quad (19) $$

Here, the unknowns $m_i^\pi$ are contained in the $\mu_i^\pi$ as initial conditions of dynamical system (19), while $\alpha_i$, $\beta_i$, $\phi$, $\psi$ characterize the dynamics of the chain and are analytically calculated in Appendix B. We stress that Equation (14) is fully general for a 3d isotropic Gaussian non-Markovian process

and does not depend on the particular structure of the mobility matrix $D$ or the connectivity matrix $M$. Then, the expression of the mean reaction time is obtained by multiplying (9) by $\delta(\mathbf{x}_N - \mathbf{x}_1)$ and integrating over all configurations, leading to

$$ \frac{T}{L^3} = \int_0^\infty \frac{dt}{\psi^{3/2}} \left[ e^{-R_\pi^2/2\psi} - \frac{Z(b, \psi)}{Z(b, L^2)} \right]. \quad (20) $$

To conclude, Eq. (19) provides a system of $N$ equations for the unknowns $\mu_i^\pi(t)$. Note that this system is of rank $N-2$, because the polymer center of mass can be set arbitrarily, and the constraint $m_N^\pi - m_1^\pi = b$ must hold. Solving this system allows us to compute the function $R_\pi(t)$, and then to calculate the mean first-passage time $T$ from (20).

We have solved numerically these equations for different values of the number of monomers $N$ and capture radius $b$. The results are shown in Fig. 1. We clearly see that the non-Markovian theory accurately predicts the mean-cyclization time, whereas the Wilemski-Fixman approach quickly fails when the capture radius is large enough. More precisely, the theory is in quantitative agreement with the simulations performed with preaveraging. Remarkably, the theory is also in quantitative agreement with the simulations performed with the exact (non-pre-averaged) interactions, in which the full non-linear stochastic dynamics is taken into account when the capture radius is not small.

The average positions $m_i^\pi$ of the monomers at the instant of cyclization are shown in Fig. 2, which shows that the Zimm chain is significantly more elongated in the direction of the end-to-end vector at the reaction than in an equilibrium looped configuration. This allows the chain to perform cyclization more rapidly than if it had to equilibrate. The monomers neighboring the reactive ones are on average outside the reactive region. Note that the small shift between the theoretical prediction and the simulation indicates that the theory is not exact. The shape of the curve of the $m_i^\pi$ is however significantly closer to the simulations for the

<figure><img src="image_2.png" /><figcaption>FIG. 1. Mean cyclization time of flexible chains with hydrodynamic interactions MFCT as a function of the number of monomers $N$ for different capture radii ($b = l_0$ in black, $b = 4l_0$ in blue and $b = 10l_0$ in red). Plain lines stand for the non-Markovian approach, dashed lines for the Markovian approach. Dots stand for simulations with the pre-averaged Rotne-Prager tensor, and triangles for simulations without any pre-averaging. The asymptotic scaling $T \sim N^{3/2}$ obtained analytically (green thick line) is well reproduced.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Average position of the monomers $m_i^\pi$ at the instant of cyclization in the direction of the end-to-end vector at this instant. We represent $m_i^\pi = \langle \mathbf{r}_i(t^*) \cdot \mathbf{u}(t^*) \rangle$, with $t^*$ the first cyclization time and $\mathbf{u} = \mathbf{r}_{ee}/|\mathbf{r}_{ee}|$ for Zimm chains in the WF approximation (red line), the non-Markovian theory (green line), simulations with pre-averaging (blue diamonds) and without (diamonds in cyan). The result of the non-Markovian theory for Rouse chains is also shown for comparison (black line). Parameters $N = 80$, $b = 4l_0$.</figcaption></figure>

non-Markovian theory than in the case of the Wilemski-Fixman approximation.

For small capture radius, the cyclization times obtained numerically without pre-averaging are slightly larger than those obtained numerically with pre-averaging, as well as those predicted by the non-Markovian theory. This can be understood from the fact that, for very small capture radius, the monomer motion at short time scales plays a key role. For small $t$, the mean square displacement (MSD) of the end-to-end distance reads $\psi(t) \approx 4(D_{11} + D_{1N})t$, and the motion is diffusive. With the pre-averaging, $D_{1N} \sim 1/\sqrt{N}$ is negligible compared to $D_{11}$. However, for chains that are close to form a loop, the actual mobility tensor reads $D_{1N} \sim D_{11}$ since the distance between the end-beads is small. Hence, the effective diffusion constant at small time scales is not correctly estimated with the pre-averaging procedure, and we can expect discrepancies for the cyclization times due to the pre-averaging for small capture radius, as observed in Fig. 1.

# IV. ASYMPTOTIC BEHAVIOR OF THE CYCLIZATION TIME

## A. Scalings of the MFCT in the Wilemski-Fixman theory

Let us now derive the different scaling behaviors of the MFCT with $N$ and the capture radius $b$, first in the Wilemski-Fixman approximation. The MFCT are then computed with Eq. (10). We first consider the limit of small capture radius, $b \to 0$ at fixed chain length $N$. If $b = 0$, integral (10) is divergent due to the linear behavior of the MSD function $\psi$ at short times. Introducing the effective short time diffusion coefficient $D_s$ (such that $\psi \approx_{t \to 0} 2D_s t$), and considering that integral (10) is governed by the short time regime, we get

$$T \approx L^3 \int_0^\infty dt \frac{e^{-b^2/4D_s t}}{(2D_s t)^{3/2}} = \frac{\sqrt{\pi} l_0^3 (N-1)^{3/2}}{\sqrt{2} D_s b}. \quad (21)$$

As discussed above, the actual value of $D_s$ could be quantitatively underestimated by the pre-averaging procedure, which computes it for equilibrium rather than looped configurations. Nevertheless, $D_s$ remains of the order of $k_B T/(6\pi\eta l_0)$. Scaling (21) is similar to that appearing for free-draining flexible (Rouse) chains.18,21,25,27 It is known that such scaling does not depend on the structure of the chain since it is the same for semi-flexible chains31 and for hyperbranched structures;40 here it is clear that it also appears in the presence of hydrodynamic interactions. Eq. (21) means that in the small capture radius regime, the MFCT is (up to a prefactor) the time needed for a diffusive particle to find a target of size $b$ in a confining volume $L^3$37,38 and does not result from the collective dynamics of the monomers.

We now consider the scaling of the MFCT in the limit of long chains. In this limit, we use the commonly used dominant diagonal approximation,2,12 in which the orthogonal matrix $Q$ that diagonalizes the Laplacian matrix $M$ is assumed to diagonalize also the product $DM$. One can check in Fig. 3 that the MSD function $\psi$ calculated with this approximation is very close to its exact value even for moderately large $N$. In this approximation, the correlation function in the continuous limit can be shown2,12 to be

$$\phi(t) \approx \sum_{p \text{ odd}} \frac{8}{p^2 \pi^2} e^{-p^{3/2}t/\tau_1}, \quad (22)$$

where $\tau_1$ is the slowest relaxation time scale of the chain,

$$\tau_1 = \frac{3\eta l_0^3}{\sqrt{\pi k_B T}} N^{3/2}. \quad (23)$$

In the limit $t/\tau_1 \ll 1$, sum (22) can be replaced by a continuous integral, leading to the identification of the short time subdiffusive behavior,

$$\psi(t) = \langle [x_{ee}(t) - x_{ee}(0)]^2 \rangle \approx \frac{8\Gamma(1/3)}{3^{2/3}\pi^{5/3}} \left( \frac{k_B T t}{\eta} \right)^{2/3}, \quad (24)$$

which shows no dependence on the bond size $l_0$ or the chain length $N$. Here $\Gamma(\cdot)$ represents the Gamma function. Defining $\tau = t/\tau_1$, and introducing the rescaled functions $\phi(t) = \Phi(t/\tau_1) = \Phi(\tau)$ and $\psi(t) = l_0^2 N \Psi(\tau)$, with $\Psi$ and $\Phi$

<figure><img src="image_4.png" /><figcaption>FIG. 3. Validity of the dominant diagonal approximation. The red and blue curves show the function $\psi(t)$ with and without the dominant diagonal approximation for $N = 100$, respectively. Inset: same figure in double logarithmic scales.</figcaption></figure>
independent on $N$ in the continuous limit, a simple change of variable in Eq. (10) leads to

$$T = \tau_1 \int_0^\infty d\tau \left[ \frac{e^{-\tilde{b}^2 \Phi(\tau)^2 / 2\Psi(\tau)}}{\Psi(\tau)^{3/2}} - 1 \right] \equiv \tau_1 f(\tilde{b}), \quad (25)$$

with $\tilde{b} = b/(l_0\sqrt{N})$ and $f$ a dimensionless function. Equation (25) finally gives the behavior of the mean cyclization time for continuous Zimm chains.

We focus now on the asymptotics of $f$ for small rescaled capture radius $\tilde{b}$, which is highly dependent of the short time behavior of $\Psi \approx \kappa\tau^{2/3}$, where $\kappa$ is deduced from Eqs. (23) and (24) and reads

$$\kappa = 8\Gamma(1/3)/\pi^2. \quad (26)$$

Because of this subdiffusive behavior, replacing $\tilde{b}$ by zero in (25) leads to a divergent integral, meaning that $f$ diverges with $\tilde{b}$ for small capture radii. The reason for this divergence is that, when $\tau \sim b^3$, the term $\tilde{b}^2/\Psi$ becomes of order 1 and cannot be replaced by 0. This suggests to introduce an intermediate “time” scale $\varepsilon$ such that $\tilde{b}^3 \ll \varepsilon \ll 1$, and to split the integral (25) into two contributions, leading to

$$f(\tilde{b}) \approx \int_0^{\frac{\varepsilon}{\tilde{b}^3}} du \frac{e^{-\frac{1}{2\kappa u^{2/3}}}}{\kappa^{3/2} u} + \int_{\varepsilon}^{\infty} d\tau \left( \frac{1}{[1 - \Phi^2]^{3/2}} - 1 \right) \quad (27)$$

in which we have used the short time expression (24) for $\tau < \varepsilon$, and simply set $\tilde{b} = 0$ in the contribution coming from the large times. In the joint limit $\varepsilon \to 0$ and $\varepsilon/\tilde{b}^3 \to \infty$, the integral (27) can be recast under the form

$$f(\tilde{b}) \approx -\frac{\ln(\tilde{b}^3)}{\kappa^{3/2}} + \int_0^\infty du \left[ \frac{e^{-\frac{1}{2\kappa u^{2/3}}}}{\kappa^{3/2} u} - \frac{\theta(u-1)}{\kappa^{3/2} u} \right] + \int_0^\infty d\tau \left( \frac{1}{[1 - \Phi(\tau)^2]^{3/2}} - 1 - \frac{\theta(1-\tau)}{\kappa^{3/2} \tau} \right) \quad (28)$$

with $\theta(\cdot)$ the Heaviside step function. All integrals appearing in this equation are convergent, and their numerical evaluation leads to the scaling form

$$f(\tilde{b}) \approx \frac{3}{\kappa^{3/2}}(-\ln \tilde{b} + 0.721 \dots). \quad (29)$$

In other words, we have identified the asymptotic scaling for the mean cyclization time of long Zimm chains,

$$T \approx \frac{9\pi^{5/2}}{[8\Gamma(1/3)]^{3/2}} \frac{\eta l_0^3 N^{3/2}}{k_B T} \ln \left( \frac{2.06 l_0 \sqrt{N}}{b} \right). \quad (30)$$

We note that the weak logarithmic dependence on the size of the target is due to the fact that the motion of a monomer is a marginally compact process in this regime (the dimension of the walk satisfies $d_w = d = 3$).$^{37,38}$ The scaling $T \sim N^{3/2} \log N$ is in particular consistent with the renormalization group approach of Ref. 22. It refines the earlier analysis of Ref. 2, based on a similar Wilemski-Fixman approach, which however did not predict the logarithmic correction. The effect of hydrodynamic interactions is clearly visible, since this scaling of $T$ with $N$ is very different from the case of flexible chains without hydrodynamic interactions, where $T$ for long chains scales as $N^2$ and is independent of the capture radius. Finally, let us note that a non-Markovian

analysis reproduces the scaling of the MFCT of Eq. (30), with however a different numerical prefactor.

## V. CASE OF AN EXTERNAL TARGET

The formalism presented above can be adapted to the case of an external target, which models an inter-molecular reaction between one monomer at the end of the chain and a fixed target of radius $b$ in a confined volume $V$. This case has already been studied for a Rouse chain,$^{28}$ and we here briefly describe how the case of a Zimm chain can be treated. Assuming that the target is at the origin, the quantity of interest is $\mathbf{R} = \mathbf{R}_1$ if we assume that the reactive monomer is the first one. We present below the scaling of the mean reaction-time $T$ with the target, for a chain confined in a large volume $V$ in Wilemski-Fixman approximation. Eq. (9) remains valid, and by integrating it over all configurations such that $\mathbf{x}_1 = \mathbf{0}$, one gets

$$\begin{aligned} T P_{\text{stat}}(\mathbf{x}_1 = \mathbf{0}) &= \int_0^\infty dt [P(\mathbf{x}_1 = \mathbf{0}|\pi, 0) - P(\mathbf{x}_1 = \mathbf{0}|\text{ini}, 0)] . \end{aligned} \quad (31)$$

The propagators appearing above are a priori the propagators in a confined volume. In the limit of large volume, we argue that $P_{\text{stat}}(\mathbf{x}_1 = \mathbf{0})$ is simply equal to $1/V$. The other propagators are well defined in infinite volume. The Wilemski-Fixman approximation then consists in assuming that the reactive conformations given by $\pi$ are simply the equilibrium configurations such that $|\mathbf{x}_1| = b$. Hence, one has

$$\frac{T}{V} = \int_{0}^{\infty} dt \left[ \frac{e^{-b^2/[2\psi(t)]}}{[2\pi\psi(t)]^{3/2}} - \frac{e^{-R_0^2/[2\psi(t)]}}{[2\pi\psi(t)]^{3/2}} \right], \quad (32)$$

where $R_0$ is the initial position of the reactive monomer and $\psi(t) = \langle [x_1(t) - x_1(0)]^2 \rangle$ the MSD of the reactive monomer, whose asymptotic behavior is

$$\psi(t) \approx \begin{cases} \frac{\alpha' t}{4\Gamma(1/3)} \left( \frac{k_B T t}{3\eta} \right)^{2/3} & \text{for } t \ll \tau_0 \\ 2D_{\text{cm}} t & \text{for } \tau_0 \ll t \ll \tau_1, \\ 0 & \text{for } \tau_1 \ll t \end{cases} \quad (33)$$

where $D_{\text{cm}}$ is the large time center-of-mass diffusion coefficient. For large $N$, the following expressions are found:$^{36}$

$$\alpha' = \frac{4 k_B T}{3\sqrt{\pi} \pi \eta l_0}, \quad (34)$$

$$D_{\text{cm}} = \frac{4\sqrt{2}}{9\pi^{3/2}} \frac{k_B T}{\eta l_0 \sqrt{N}}. \quad (35)$$

The expression of the mean reaction time $T$ in the Wilemski-Fixman approximation can then be deduced from Eq. (32) and yields in the limit of chains starting far from the target ($R_0 \to \infty$),

$$T \sim \begin{cases} V/(4\pi D_{\text{cm}}b) & \text{for } \sqrt{N}l_0 \ll b \\ \nu \frac{V\eta}{k_B T} \log\left(\frac{\sqrt{N}l_0}{b}\right) & \text{for } l_0 \ll b \ll \sqrt{N}l_0, \\ \frac{V}{(2\pi\alpha'b)} & \text{for } b \ll l_0 \end{cases} \quad (36)$$
where v is a numerical coefficient that reads

$$v = \frac{9\pi}{8^{3/2}\Gamma(1/3)^{3/2}}. \qquad (37)$$

As in the case of cyclization dynamics, the scaling of the mean reaction time with N is different from the case of a Rouse chain. In particular, for intermediate capture radii b, the very weak (logarithmic) dependence on N leads to reaction times that can be significantly shorter than for a Rouse polymer, for which$^{28}$ $T \propto V\sqrt{N}$. The non-Markovian analysis of the reaction time with an external target closely follows the steps developed above for the determination of the MFCT and is not described here. As in the case of the cyclization dynamics, the scaling laws with b and N would be the same as predicted by the Wilemski-Fixman treatment (up to numerical prefactors), given in Eq. (36).

## VI. CONCLUSION

In conclusion, we have presented in this paper a theoretical analysis of the closure kinetics of a polymer with hydrodynamic interactions. We have provided both a Markovian (Wilemski-Fixman) and a non-Markovian analytical approach, based on a recent method introduced in Refs. 21 and 27–29 Although our theory relies on the preaveraging of the mobility tensor (Zimm dynamics), it is shown to reproduce very accurately the results of numerical simulations of the complete nonlinear dynamics. It is found that the Markovian treatment significantly overestimates cyclization times (up to a factor 2), showing the importance of memory effects in the dynamics. Such non-Markovian effects can be understood by analyzing the distribution of the polymer conformations at the instant of reaction, which is found to significantly depart from the equilibrium distribution.

In addition, we derived asymptotic expressions of the mean cyclization time with the polymer size N and capture radius b, which are identical in both Markovian and non-Markovian approaches, but with different prefactors. We computed the precise values of the prefactors in the case of the Wilemski-Fixman approach. In particular, it is found that the scaling of the MFCT for large N is given by $T \sim N^{3/2} \ln(N/b^2)$, whereas for the Rouse chain, one has $T \sim N^2$ (see Ref. 29 for review). Hydrodynamic interactions therefore change both the dependence on N and b. This difference comes from the fact that subdiffusive exponent that characterizes the monomer dynamics at intermediate length scales is different in both models.

The present work demonstrates that the physics of cyclization kinetics in realistic models of polymers can be described by taking into account non-Markovian effects, which turn out to be much more important than the errors due to the approximate treatment of hydrodynamic interactions.

## ACKNOWLEDGMENTS

N.L., T.G., O.B., and R.V. acknowledge the support of the Campus France (Project No. 28252XE) and of the European Research Council starting Grant No. FPTOpt-277998. M.D. and A.B. acknowledge the support of the DAAD through the

PROCOPE program (Project No. 55853833) and of the DFG through Grant No. B1 142/11-1 and through IRTG Soft Matter Science (Grant No. GRK 1642/1).

## APPENDIX A: PRE-AVERAGED ROTNE-PRAGER TENSOR

We give in this section some details about the pre-averaging of the Rotne-Prager tensor, in order to quantify the differences obtained with the pre-averaged Oseen tensor. Let us first remind that the hydrodynamic interactions do not modify the equilibrium state of the polymer, which is that of a flexible Gaussian chain, hence the equilibrium probability density of the vector $\mathbf{r}_{ij}$ is

$$P_{\text{eq}}(\mathbf{r}_{ij}) = (2\pi\sigma_{ij})^{-3/2} e^{-r_{ij}^2/(2\sigma_{ij})}, \qquad (\text{A1})$$

where $\sigma_{ij} = l_0^2|i-j|$ is the variance of each spatial coordinate of $\mathbf{r}_{ij}$ at equilibrium. The pre-averaged Rotne-Prager tensor is defined as

$$\bar{\mathbf{D}}_{ij} = \int d\mathbf{r}_{ij} P_{\text{eq}}(\mathbf{r}_{ij}) \mathbf{D}_{ij}(\mathbf{r}_{ij}). \qquad (\text{A2})$$

We note that the average over all orientations of the non-isotropic tensor appearing in the expression of $\mathbf{D}_{ij}$ is

$$\left\langle \mathbf{r}_{ij} \otimes \mathbf{r}_{ij} \right\rangle = \frac{1}{3} \mathbf{I}. \qquad (\text{A3})$$

Using this relation, the integral over rotational degrees of freedom can be performed in Eq. (A2); using Eqs. (3) and (4), we obtain an integral over radial components only,

$$\bar{\mathbf{D}}_{ij} = \mathbf{I} \int_{0}^{\infty} d\rho \frac{4\pi\rho^2 e^{-\frac{\rho^2}{2\sigma_{ij}^2}}}{(2\pi\sigma_{ij})^{3/2}} \left[ \frac{H(\rho - 2a)}{6\pi\eta\rho} + \frac{H(2a - \rho)}{6\pi\eta a} \left(1 - \frac{\rho}{4a}\right) \right], \qquad (\text{A4})$$

with H the Heaviside step function. The result of this integral in the limit of small a is

$$\bar{\mathbf{D}}_{ij} \approx \frac{1}{6\pi\eta} \sqrt{\frac{2}{\pi\sigma_{ij}}} \left(1 - \frac{a^2}{3\sigma_{ij}}\right) \mathbf{I}. \qquad (\text{A5})$$

Note that for $a = 0$, one recovers pre-averaged Oseen tensor (6). The corrections are of order $a^2/(3l_0^2|i-j|)$; with our choice of parameter $a = l_0/4$ these corrections are always less than ≈2%, which justifies to chose to take the averaged Oseen tensor instead of the averaged Rotne-Prager tensor in the analysis.

## APPENDIX B: EXPLICIT EXPRESSIONS OF THE DYNAMIC QUANTITIES $\alpha_i, \beta_i, \phi, \psi$

We consider the Fokker-Planck equation describing the chain dynamics with the pre-averaged mobility tensor. Since this tensor is isotropic, we consider the dynamics of a single spatial component, say $x_i$, of the vector positions $\mathbf{x}_i$. To the Langevin equation (2), we can associate the Fokker-Planck
equation

$$
\frac{\partial P}{\partial t} = \sum_{i,j=1}^{N} \frac{\partial}{\partial x_i} D_{ij} \left( k_B T \frac{\partial P}{\partial x_j} + \sum_{k=1}^{N} M_{jk} x_k P \right). \quad (B1)
$$

By replacing the mobility tensor by its pre-averaged version, Gaussian solutions do exist. The evolution of the mean vector $\mu = (\mu_1, \dots, \mu_N)$ and the covariance matrix $\Gamma_{ij}$ of $x_i, x_j$ with time can be found by multiplying (B1) by $x_i$ or $x_m x_n$, integrate over all $x_k$, and use the divergence theorem, resulting in Ref. 39,

$$
\frac{d\mu}{dt} = -K\mu, \qquad (B2)
$$

$$
\frac{d\Gamma}{dt} = 2k_B T D - K\Gamma - \Gamma'^{\prime}K, \quad (B3)
$$

where $K = DM$ is the product of the pre-averaged mobility matrix $D$ by the connectivity matrix $M$. To solve these equations, it is convenient to introduce modes, which diagonalize $K$ (which ones coincide with Rouse modes if $D_{ij} \propto \delta_{ij}$). We define an invertible matrix $P$ such that $PKP^{-1} = \text{Diag}(\nu_1, \dots, \nu_N)$ is diagonal, with $0 = \nu_1 < \nu_2 < \nu_3 < \dots$. The vanishing eigenvalue is associated to the translational motion of the chain, while all other eigenvalues describe the internal conformational degrees of freedom. We consider the amplitudes of the Zimm modes, defined by

$$
a_i = \sum_{j=1}^{N} P_{ij} x_j \qquad (B4)
$$

and the associated mean vector $\mu_i = \langle a_i \rangle$ and covariance matrix $\kappa_{ij} = \text{Cov}(a_i, a_j)$ read

$$
u = P\mu, \quad \kappa = P \Gamma'^{\prime} P, \quad B = 2k_B T P D \Gamma^{-1} P. \quad (B5)
$$

Note that $P$ is not an orthogonal matrix. The evolution of $u$
and $\kappa$ reads

$$
\frac{d u_i}{dt} = -v_i u_i, \qquad (B6)
$$

$$
\frac{d\kappa_{ij}}{dt} = -(\nu_i + \nu_j)\kappa_{ij} + B_{ij}. \quad (B7)
$$

The solutions are

$$
u_i(t) = u_i(0) e^{-\nu_i t}, \quad (B8)
$$

$$
\kappa_{ij}(t) = \left( \kappa_{ij}(0) - \frac{B_{ij}}{\nu_i + \nu_j} \right) e^{-(\nu_i + \nu_j)t} + \frac{B_{ij}}{\nu_i + \nu_j}. \quad (B9)
$$

At large times, the covariance matrix reaches the stationary
value (for i, j ≥ 2),

$$
\kappa_{ij}^{s} = \frac{B_{ij}}{\nu_i + \nu_j}. \qquad (B10)
$$

We introduce the set of coefficients $\tilde{c}$,

$$
\tilde{c}_k = (P^{-1})_{kN} - (P^{-1})_{k1} = P^{-1}h, \quad (B11)
$$

where we remind that $h = (-1, 0, \dots, 0, 1)^t$. The vector $\tilde{c}$ is
such that

$$
x_{ee} = \sum_{i=1}^{N} \tilde{c}_i a_i . \qquad (B12)
$$

Note that $\tilde{c}_1 = 0$, meaning that the motion of the drift center is not involved in the evolution of the internal variable $x_{ee}$. The

equilibrium end-to-end distance is

$$
L^2 = \lim_{t \to \infty} {}^t h \Gamma h = \sum_{m,n=2}^{N} \tilde{c}_m \kappa_{mn}^s \tilde{c}_n. \quad (B13)
$$

Similarly, we obtain for $\alpha_i$

$$
\alpha_i = \sum_{j=2}^{N} \sum_{k=1}^{N} (P^{-1})_{ij} \kappa_{jk}^{s} \tilde{c}_k. \quad (B14)
$$

Consider the correlation matrix $C_{ij}(t) = \langle a_i(t)a_j(0) \rangle$, starting from an equilibrium configuration. This matrix is a solution of

$$
\partial_t C_{ij} = -v_i C_{ij}; \quad C_{ij}(0) = \kappa_{ij}^s. \qquad (B15)
$$

The solution is straightforward: $C_{ij} = \kappa_{ij}^s e^{-\nu_i t}$. Given that $\langle x_{ee}(t)x_{ee}(0) \rangle = \sum_{ij} \tilde{c}_i \tilde{c}_j C_{ij}$, we obtain for the function $\phi$ [defined in Eq. (11)],

$$
\phi(t) = \frac{1}{L^2} \sum_{m,n=2}^{N} \tilde{c}_n \kappa_{mn}^{s} \tilde{c}_m e^{-\nu_m t}. \quad (B16)
$$

Next, we denote $\kappa^*$ the covariance matrix of $a_i, a_j$ at equilibrium with the constraint of fixed $x_{ee}$, which from (C3) reads

$$
\kappa_{ij}^* = \kappa_{ij}^s - \sum_{m,n=1}^{N} \frac{\kappa_{im}^s \kappa_{jn}^s \tilde{c}_m \tilde{c}_n}{L^2}. \quad (B17)
$$

Consider now

$$
\kappa_{ij}(t) = \mathrm{Cov}(a_i(t), a_j(t) | x_{ee}(0) = 0) \quad (B18)
$$

which is related to $\beta_i$ by

$$
\beta_i = \sum_{j,k=1}^{N} (P^{-1})_{ij} \kappa_{jk} \tilde{c}_k. \quad (B19)
$$

Taking covariance matrix (B17) as an initial condition for
dynamics (B7), we obtain

$$
\kappa_{ij}(t) = \kappa_{ij}^{s} - \frac{\sum_{n} \kappa_{in}^{s} \tilde{b}_{n} \sum_{m} \kappa_{jm}^{s} \tilde{b}_{m}}{L^2} e^{-(\nu_{i} + \nu_{j})t}. \quad (B20)
$$

Using (B16) and (B19), we obtain for $\beta_i$

$$
\beta_i = \alpha_i - \phi(t) \sum_j \sum_k P_{ij}^{-1} \kappa_{jk}^s e^{-\nu_j t} \tilde{b}_k. \quad (B21)
$$

Finally, the function $\psi = \beta_N - \beta_1$ reads

$$
\psi(t) = L^2(1 - \phi^2(t)), \quad (B22)
$$

where we have again used (B16). Hence, the expressions of all the dynamical quantities $\alpha_i, \beta_i, \phi, \psi$ are given explicitly in this section.

APPENDIX C: DERIVATION OF THE
SELF-CONSISTENT EQUATION (14)

We present here a derivation of the set of equations (14)
for the moments $m_i^\pi$, which adapts the method used for Rouse
polymers.²⁷ First, we multiply (9) by $x_{iz}\delta(\mathbf{x}_N - \mathbf{x}_1 - \mathbf{R}_f\hat{\mathbf{e}}_z)$
(where $\hat{\mathbf{e}}_z$ is the unit vector in the vertical direction, and
$R_f$ is fixed and satisfies $0 < R_f < b$) and integrate over all configurations, to get

$$ 
\begin{aligned}
& T p_{\text{stat}}(R_f \hat{\mathbf{e}}_z) \mathbb{E}_{\text{stat}}(x_{iz} | R_f) \\
& = \int_0^\infty dt \int d\Omega [p(R_f \hat{\mathbf{e}}_z, t | \pi_\Omega, 0) \mathbb{E}(x_{iz}, t | R_f \hat{\mathbf{e}}_z, t; \pi_\Omega, 0) \\
& \quad - p(R_f \hat{\mathbf{e}}_z, t | P_{\text{ini}, \Omega}, 0) \mathbb{E}(x_{iz}, t | R_f \hat{\mathbf{e}}_z, t; P_{\text{ini}, \Omega}, 0)],
\end{aligned} 
\quad \text{(C1)}
$$

with $p$ representing the probability distribution function of the end-to-end vector, $p(R_f \hat{\mathbf{e}}_z, t | \pi_\Omega, 0)$ is the probability that $\mathbf{r}_{ee} = R_f \hat{\mathbf{e}}_z$ at $t$ with an initial distribution $\pi_\Omega$, and $\mathbb{E}(x_{iz}, t | R_f \hat{\mathbf{e}}_z, t; \pi_\Omega, 0)$ is the conditional average of the $z$ coordinate of $\mathbf{x}_i$ at $t$ given that $\mathbf{r}_{ee} = R_f \hat{\mathbf{e}}_z$ at the same time $t$ and starting from $\pi_\Omega$ initially; other notations are similar. The mention “$P_{\text{ini},\Omega}$” means the equilibrium distribution with $\mathbf{r}_{ee} = R_0 \mathbf{u}_r(\Omega)$ initially, we perform the average over initial end-to-end distances at the end of the calculation.

Let us remind the following formula for the conditional mean of a Gaussian variable $X$ given that a second Gaussian variable $Y$ takes the value $Y_0$:

$$ \mathbb{E}(X|Y = Y_0) = \mathbb{E}(X) - \frac{\text{Cov}(X, Y)}{\text{Cov}(Y, Y)} [\mathbb{E}(Y) - Y_0]. \quad \text{(C2)} $$

A similar general formula for the conditional covariances of Gaussian variables is

$$ 
\begin{aligned}
\text{Cov}(X_1, X_2 | Y = Y_0) \\
= \text{Cov}(X_1, X_2) - \frac{\text{Cov}(X_1, Y) \text{Cov}(X_2, Y)}{\text{Cov}(Y, Y)}.
\end{aligned} 
\quad \text{(C3)}
$$

Now, we consider a fixed angular direction $\Omega$ and define $\theta$ the angle with the vertical direction. We write

$$ x_{iz} = x_{ir} \cos \theta - x_{i\theta} \sin \theta, \quad \text{(C4)} $$

with $x_{ir}$ the component of $\mathbf{x}_i$ in the direction $\hat{\mathbf{u}}(\Omega)$. We note that conditioning the end-to-end vector to have the value $R_f \hat{\mathbf{e}}_z$ imposes that its component in the direction $\hat{\mathbf{u}}(\Omega)$ takes the value $R_f \cos \theta$. Then, applying (C2) and using the definitions of $R_\pi, \beta_i, \psi$, we get

$$ \mathbb{E}(x_{ir}, t | R_f \hat{\mathbf{e}}_z, t, \pi_\Omega, 0) = \mu_i^\pi - \frac{\beta_i}{\psi} (R_\pi - R_f \cos \theta). \quad \text{(C5)} $$

Applying the same reasoning in the direction normal to $\mathbf{u}(\Omega)$, we get

$$ \mathbb{E}(x_{i\theta}, t | R_f \hat{\mathbf{e}}_z, t, \pi_\Omega, 0) = -\frac{\beta_i}{\psi} R_f \sin \theta. \quad \text{(C6)} $$

Inserting these expressions into (C4) leads to

$$ \mathbb{E}(x_{iz}, t | R_f \hat{\mathbf{e}}_z, t, \pi_\Omega, 0) = \cos \theta \left( \mu_i^\pi - \frac{\beta_i}{\psi} R_\pi \right) + \frac{R_f \beta_i}{\psi}. \quad \text{(C7)} $$

Let us pose now

$$ \mu_i^{\text{stat}, R_0} = \mathbb{E}(x_{iz}(t) | x_{ee}(0) = R_0, 0), \quad \text{(C8)} $$

$$ R = \mathbb{E}(x_{ee}(t) | x_{ee}(0) = R_0, 0), \quad \text{(C9)} $$

where initial equilibrium conditions (apart from the constraint for $x_{ee}(0)$) are understood. A reasoning similar to that leading to Eq. (C7) gives

$$ 
\begin{aligned}
& \mathbb{E}(x_{iz}, t | R_f \hat{\mathbf{e}}_z, t; P_{\text{ini}, \Omega}, 0) \\
& = \cos \theta \left( \mu_i^{\text{stat}, R_0} - \frac{\beta_i}{\psi} R \right) + \frac{R_f \beta_i}{\psi}.
\end{aligned} 
\quad \text{(C10)}
$$

Similarly, from (C2) and (16), we obtain

$$ \mathbb{E}_{\text{stat}}(x_{iz} | R_f \hat{\mathbf{e}}_z) = \frac{\alpha_i}{L^2} R_f. \quad \text{(C11)} $$

Finally, the propagators for the end-to-end distance read

$$ p(R_f \hat{\mathbf{e}}_z, t | P_{\text{ini}, \Omega}, 0) = \frac{1}{(2\pi\psi)^{3/2}} e^{-\frac{(R_f \mathbf{u}_z - R_f \mathbf{u}_r)^2}{2\psi}}, \quad \text{(C12)} $$

$$ p(R_f \hat{\mathbf{e}}_z, t | \pi_\Omega, 0) = \frac{1}{(2\pi\psi)^{3/2}} e^{-\frac{(R_f \mathbf{u}_z - R_f \pi \mathbf{u}_r)^2}{2\psi}}, \quad \text{(C13)} $$

$$ p_{\text{stat}}(R_f \hat{\mathbf{e}}_z) = \frac{1}{(2\pi L^2)^{3/2}} \exp\left(-\frac{R_f^2}{2L^2}\right). \quad \text{(C14)} $$

All the terms appearing in (C1) have been evaluated, the self-consistent equation becomes

$$ 
\begin{aligned}
& \int_0^\infty dt \int_0^\pi d\theta \frac{\sin \theta}{2} \left[ \frac{1}{(2\pi\psi)^{3/2}} e^{-\frac{(R_f \mathbf{u}_z - R_f \mathbf{u}_r)^2}{2\psi}} \right. \\
& \qquad \left. \times \left( \cos \theta \mu_i^\pi(t) - \frac{\beta_i}{\psi} (\cos \theta R_\pi - R_f) \right) \right. \\
& \qquad \left. - \frac{1}{(2\pi\psi)^{3/2}} e^{-\frac{(R_f \mathbf{u}_z - R_f \pi \mathbf{u}_r)^2}{2\psi}} \right. \\
& \qquad \left. \times \left( \cos \theta \mu_i^{\text{stat}, R_0}(t) - \frac{\beta_i}{\psi} (\cos \theta R - R_f) \right) \right] \\
& = T R_f \frac{\alpha_i}{L^2} p_{\text{stat}}(R_f).
\end{aligned} 
\quad \text{(C15)}
$$

This equation should be verified for any $R_f$ between 0 and $b$. We choose to write it in the limit $R_f \to 0$ (more precisely, we develop both expressions at first order in $R_f$). Noting that $\hat{\mathbf{u}}_r \cdot \hat{\mathbf{e}}_z = \cos \theta$, the integration over $\theta$ can be performed, leading to

$$ 
\begin{aligned}
& \int_0^\infty dt \left\{ \frac{e^{-R_\pi^2/2\psi}}{(2\pi\psi)^{3/2}} \left[ \frac{\beta_i}{\psi} + \frac{R_\pi}{3\psi} \left( \mu_i^\pi - \frac{\beta_i}{\psi} R_\pi \right) \right] \right. \\
& \qquad \left. - \frac{e^{-R^2/2\psi}}{(2\pi\psi)^{3/2}} \left[ \frac{\beta_i}{\psi} + \frac{R}{3\psi} \left( \mu_i^{\text{stat}, R_0} - \frac{\beta_i}{\psi} R \right) \right] \right\} \\
& = T \frac{\alpha_i}{L^2} \frac{1}{(2\pi L^2)^{3/2}}.
\end{aligned} 
\quad \text{(C16)}
$$

We can eliminate $T$ by replacing its expression

$$ \frac{T}{L^3} = \int_0^\infty \frac{dt}{\psi^{3/2}} \left\{ e^{-R_\pi^2/2\psi} - e^{-R^2/2\psi} \right\}. \quad \text{(C17)} $$

Let us now find the dependence of $R$ and $\mu_i^{\text{stat}, R_0}$ with $R_0$. From (C8) and (C3), we deduce

$$ \mu_i^{\text{stat}, R_0} = R_0 \frac{\text{Cov}(x_i(t), x_{ee}(0))}{L^2}. \quad \text{(C18)} $$

Furthermore, applying the formula for the conditional covariances in (15), we get

$$ 
\begin{aligned}
\beta_i = \text{Cov}(x_i(t)x_{ee}(t)) & \\
& - \frac{\text{Cov}(x_i(t)x_{ee}(0))\text{Cov}(x_{ee}(t)x_{ee}(0))}{\text{Cov}(x_{ee}(0)x_{ee}(0))}.
\end{aligned} 
\quad \text{(C19)}
$$

From which we obtain

$$ \mu_i^{\text{stat}, R_0} = \frac{R_0(\alpha_i - \beta_i)}{L^2 \phi(t)} \quad \text{(C20)} $$
and

$$R = \mu_{N}^{\text{stat}, R_{0}} - \mu_{1}^{\text{stat}, R_{0}} = R_{0}\phi(t). \quad (\text{C21})$$

The last step consists in adapting the reasoning to the case of initial end-to-end distance that is distributed $p_r(R_0) = R_0^2 e^{-R_0^2/2L^2} / \int_b^\infty r^2 e^{-r^2/2L^2} dr$. This can be achieved by keeping an average over $R_0$ at all steps of the derivation. Averaging (C16) with respect to $R_0$ [by using (C20) and (C21)] leads to the final form of self-consistent equation (14) for the moments given in the main text.

## APPENDIX D: NUMERICAL METHODS

We performed the numerical integration of stochastic equation (2) by using either the full Rotne-Prager tensor or its pre-averaged form. At each time step of size $\Delta t$, the positions of the monomers evolve according to the algorithm of Ermak and McCammon,<fnref n="32" />

$$\mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + \sum_{j=1}^{N} \mathbf{D}_{ij}(t)\mathbf{F}_j(t)\Delta t + \xi_i(\Delta t), \quad (\text{D1})$$

with the same notations as before. $\xi_i(\Delta t)$ is a random Gaussian noise with zero mean and covariance $\langle \xi_i(\Delta t) \xi_j(\Delta t) \rangle = 2k_B T \mathbf{D}_{ij}(t)\Delta t$. The generation of the 3N random numbers $\xi_i$ requires to find a Cholesky decomposition of the mobility tensor, which can be done in $N^3$ operations. In the case of simulations using pre-averaged mobility tensor (6), this decomposition needs to be performed only once, whereas simulations using the full Rotne-Prager tensor require to perform a Cholesky decomposition at each time step, resulting in considerably longer computational times which prevented us to explore the same range of parameters. Initial configurations are generated from the equilibrium Gibbs-Boltzmann Gaussian distribution $P_{\text{eq}}(\{\mathbf{x}\}) \propto e^{-k \sum_{i=1}^{N-1} (\mathbf{x}_{i+1} - \mathbf{x}_i)^2 / 2k_B T}$ and are rejected if the condition $|\mathbf{x}_N - \mathbf{x}_1| > b$ is not satisfied. Once an equilibrium configuration is generated, the positions evolve through (D1) until reaching a configuration $|\mathbf{x}_N - \mathbf{x}_1| \le b$, the cyclization time $t_i$ for this run is recorded. The MFCT is finally found by ensemble averaging $t_i$ over many runs. The time step is chosen as suggested by Pastor et al.<fnref n="18" /> Noting $R$ the end-to-end distance,

$$\Delta t = \Delta_{\text{low}} + \Delta_{\text{high}} \sin \left( \frac{\pi}{6} (R^2 - b^2) \right) \quad (\text{D2})$$

if $R^2 < b^2 + 3l_0$, and

$$\Delta t = \Delta_{\text{low}} + \Delta_{\text{high}} \quad (\text{D3})$$

if not.

We chose $\Delta_{\text{high}} = \sqrt{N}10^{-5}$ and $\Delta_{\text{low}} = 10^{-5}$, and we controlled that convergence was reached. The results of simulations are represented in Fig. 1.

<footnote n="1">¹M. Ortiz-Repiso, J. J. Freire, and A. Rey, *Macromolecules* **31**, 8356–8362 (1998).</footnote>

<footnote n="2">²R. Chakrabarti, *Physica A* **391**, 4081 (2012).</footnote>

<footnote n="3">³G. Bonnet, O. Krichevsky, and A. Libchaber, *Proc. Natl. Acad. Sci. U. S. A.* **95**, 8602 (1998).</footnote>

<footnote n="4">⁴X. Wang and W. M. Nau, *J. Am. Chem. Soc.* **126**, 808 (2004).</footnote>

<footnote n="5">⁵M. I. Wallace, L. Ying, S. Balasubramanian, and D. Klenerman, *Proc. Natl. Acad. Sci. U. S. A.* **98**, 5584 (2001).</footnote>

<footnote n="6">⁶J.-F. Allemand, S. Cocco, N. Douarche, and G. Lia, *Eur. Phys. J. E* **19**, 293 (2006).</footnote>

<footnote n="7">⁷L. J. Lapidus, W. A. Eaton, and J. Hofrichter, *Proc. Natl. Acad. Sci. U. S. A.* **97**, 7220 (2000).</footnote>

<footnote n="8">⁸A. Möglich, K. Joder, and T. Kiefhaber, *Proc. Natl. Acad. Sci. U. S. A.* **103**, 12394 (2006).</footnote>

<footnote n="9">⁹J. K. Gooden, M. L. Gross, A. Mueller, A. D. Stefanescu, and K. L. Wooley, *J. Am. Chem. Soc.* **120**, 10180 (1998).</footnote>

<footnote n="10">¹⁰Y. Zheng, H. Cao, B. Newland, Y. Dong, A. Pandit, and W. Wang, *J. Am. Chem. Soc.* **133**, 13130 (2011).</footnote>

<footnote n="11">¹¹A. Burgath, A. Sunder, and H. Frey, *Macromol. Chem. Phys.* **201**, 782 (2000).</footnote>

<footnote n="12">¹²M. Doi and S. F. Edwards, *The Theory of Polymer Dynamics* (Oxford Science Publications, 1988).</footnote>

<footnote n="13">¹³D. Panja, *J. Stat. Mech.: Theory Exp.* **2010**, P06011.</footnote>

<footnote n="14">¹⁴P.-G. De Gennes, *J. Chem. Phys.* **76**, 3316 (1982).</footnote>

<footnote n="15">¹⁵G. Wilemski and M. Fixman, *J. Chem. Phys.* **60**, 866–877 (1974).</footnote>

<footnote n="16">¹⁶G. Wilemski and M. Fixman, *J. Chem. Phys.* **60**, 878 (1974).</footnote>

<footnote n="17">¹⁷A. Szabo, K. Schulten, and Z. Schulten, *J. Chem. Phys.* **72**, 4350 (1980).</footnote>

<footnote n="18">¹⁸R. W. Pastor, R. Zwanzig, and A. Szabo, *J. Chem. Phys.* **105**, 3878 (1996).</footnote>

<footnote n="19">¹⁹I. M. Sokolov, *Phys. Rev. Lett.* **90**, 080601 (2003).</footnote>

<footnote n="20">²⁰A. E. Likhman and C. M. Marques, *Europhys. Lett.* **75**, 971 (2006).</footnote>

<footnote n="21">²¹T. Guerin, O. Benichou, and R. Voituriez, *Nat. Chem.* **4**, 568 (2012).</footnote>

<footnote n="22">²²B. Friedman and B. O'Shaughnessy, *Phys. Rev. A* **40**, 5950 (1989).</footnote>

<footnote n="23">²³R. Chakrabarti, *Physica A* **391**, 5326 (2012).</footnote>

<footnote n="24">²⁴P. Bhattacharyya, R. Sharma, and B. J. Cherayil, *J. Chem. Phys.* **136**, 234903 (2012).</footnote>

<footnote n="25">²⁵A. Amitai, I. Kupka, and D. Holcman, *Phys. Rev. Lett.* **109**, 108302 (2012).</footnote>

<footnote n="26">²⁶N. M. Toan, G. Morrison, C. Hyeon, and D. Thirumalai, *J. Phys. Chem. B* **112**, 6094 (2008).</footnote>

<footnote n="27">²⁷T. Guerin, O. Benichou, and R. Voituriez, *J. Chem. Phys.* **138** (2013).</footnote>

<footnote n="28">²⁸T. Guerin, O. Benichou, and R. Voituriez, *Phys. Rev. E* **87** (2013).</footnote>

<footnote n="29">²⁹O. Bénichou, T. Guérin, and R. Voituriez, *J. Phys. A: Math. Theor.* **48**, 163001 (2015).</footnote>

<footnote n="30">³⁰B. H. Zimm, *J. Chem. Phys.* **24**, 269–278 (1956).</footnote>

<footnote n="31">³¹T. Guerin, M. Dolgushev, O. Benichou, R. Voituriez, and A. Blumen, *Phys. Rev. E* **90**, 052601 (2014).</footnote>

<footnote n="32">³²D. L. Ermak and J. A. McCammon, *J. Chem. Phys.* **69**, 1352 (1978).</footnote>

<footnote n="33">³³T. Geyer and U. Winter, *J. Chem. Phys.* **130**, 114905 (2009).</footnote>

<footnote n="34">³⁴K. Osaki, *Macromolecules* **5**, 141 (1972).</footnote>

<footnote n="35">³⁵K. Osaki, J. L. Schrag, and J. D. Ferry, *Macromolecules* **5**, 144 (1972).</footnote>

<footnote n="36">³⁶I. Teraoka, *Polymer Solutions* (Wiley, 2002).</footnote>

<footnote n="37">³⁷S. Condamín, O. Bénichou, V. Tejedor, R. Voituriez, and J. Klafter, *Nature* **450**, 77 (2007).</footnote>

<footnote n="38">³⁸O. Bénichou and R. Voituriez, *Phys. Rep.* **539**, 225 (2014).</footnote>

<footnote n="39">³⁹N. Van Kampen, *Stochastic Processes in Physics and Chemistry*, 3rd ed. (North-Holland Personal Library, Amsterdam, 1992).</footnote>

<footnote n="40">⁴⁰M. Dolgushev, T. Guérin, A. Blumen, O. Bénichou, and R. Voituriez, *Phys. Rev. Lett.* **115**, 208301 (2015).</footnote>