Analytical expressions for the closure probability of a stiff wormlike chain for finite capture radius

Estimating the probability that two monomers of the same polymer chain are close together is a key ingredient to characterize intramolecular reactions and polymer looping. In the case of stiff wormlike polymers (rigid fluctuating elastic rods), for which end-to-end encounters are rare events, we derive an explicit analytical formula for the probability η(rc) that the distance between the chain extremities is smaller than some capture radius rc. The formula is asymptotically exact in the limit of stiff chains, and it leads to the identification of two distinct scaling regimes for the closure factor, originating from a strong variation of the fluctuations of the chain orientation at closure. Our theory is compatible with existing analytical results from the literature that cover the cases of a vanishing capture radius and of nearly fully extended chains.


I. INTRODUCTION
Contact formation and reactions between distant sites of a macromolecule are ubiquitous in nature.For example, DNA looping is a key process in the regulation of gene expression [1][2][3][4][5][6], and loops are also involved in the folding pathway of polypeptide chains [7] as well as in the higher-order structure of proteins [8] and chromatin in eukariots [9].Since stiffer chains have rarer looping events, chain closure experiments can also be used to probe the elastic properties of macromolecules as diverse as carbon nanotubes [10], wormlike micelles [11], or DNA [12][13][14][15][16][17][18][19].
An important step in characterizing contact formation consists of quantifying the probability that two given monomers, say the end monomers, are closer than some capture radius r c at equilibrium, possibly under orientational constraints at the chain extremities [3,20,21].This closure probability gives access to the looping time in the reaction controlled regime, when many end-to-end encounters are necessary to actually form a link; it is also an important quantity in all studies of first contact kinetics in polymers [19,[22][23][24][25][26].In general, the looping probability can depend on geometrical constraints such as twist alignment restriction (for example due to DNA helical repeat [13,27]) or constraints on the opening angle (appearing in particular in the situation of protein-mediated contact in DNA [16,[28][29][30][31][32][33]).It is also interesting to understand looping in the absence of orientational constraints.In this case, one aims to answer the following simple question: what is the probability η(r c ) to find the end monomers closer than the capture radius r c ?The quantity η(r c ), on which we focus in this work, is relevant, for example, when closure occurs through the pairing of (flexible) single-stranded DNA [15,18] or via a protein bridge [3,34,35] in the limit of a very soft protein; η(r c ) can also be used to derive upper bounds for the looping probability with orientational constraints [15,18].
We consider here the calculation of η(r c ) in polymers represented as thin elastic rods (wormlike chains [36]).In this simple model, which realistically describes a number of polymers such as DNA [37][38][39][40][41], nonlinearities make the calculation of the looping probability nontrivial.Analytical formulas for η (or, equivalently, for the distribution of endto-end vectors) have been proposed for very long chains (see, e.g., Refs.[27,41,42]) or chains that are only weakly bent [43] (a situation that is not relevant to the closure problem).Analytical results also exist in the limit of rigid chains, where noise-induced end-to-end contacts become rare events, since they can occur only after overcoming a large barrier of bending energy.In this regime, Shimada and Yamakawa (SY) [27] provided analytical and asymptotically exact expressions for the looping probability that are valid in the limit of an infinitely small capture radius.Other theoretical approaches have been proposed to calculate the looping probability in more general situations: exact numerical results have been obtained by using Monte Carlo [17,34,[44][45][46][47][48] or Brownian dynamics [49] simulations, transfer-matrix methods [50,51], or by exploiting the analogy in Fourier space with a quantum-mechanical problem [35,52,53], with use of infinite continued fractions [54][55][56][57] or recursive relations [58].Other approaches use weak-noise approximations to obtain the closure likelihood in discrete [31,59,60] or continuous models via path-integral techniques (in two dimensions) [61], or numerical evaluation of functional determinants [62].Approximate extrapolation formulas of SY's expressions have also been proposed [35,63].However, for stiff chains, explicit and asymptotically exact formulas for the dependence of the encounter probability η(r c ) with the capture radius in three dimensions, with an identification of scalings, are not provided by these approaches, and they are the subject of the present work.
In this paper, we present an analytical calculation of the end-to-end encounter probability η(r c ) for stiff wormlike chains, and for the distribution p(r) of end-to-end vectors r [see Eqs.(3) and ( 9)].Our analytical formulas are asymptotically exact at leading order in the limit of stiff chains, when the end-to-end distance is not too close to the chain contour length.The formulas are in excellent agreement with the numerical results obtained in Ref. [55] for a broad range of parameters, and they lead to the identification of two distinct scalings.The presence of these distinct scaling regimes can be associated with a strong variation of the fluctuations of the chain orientation at closure.The identification of these distinct scalings, together with the explicit analytical expressions for η(r c ) and p(r), are the main results of this work.
The outline of this paper is as follows.In Sec.II, we define the model and summarize the analytical formulas derived in the paper, including their physical interpretation.In Sec.III, we briefly recall the equations describing the looping configuration of minimal bending energy.We calculate p(r) in the case of finite r in Sec.IV.The case of a smaller capture radius is considered in Sec.V.In Sec.VI, we provide explicit expressions for the probability η(r c ) to observe the end monomers separated by a distance less than r c , and the associated radius-dependent closure factor.In Sec.VII, we consider a model in which torsional stiffness is also included, and we check that η(r c ) is not modified by this additional parameter.Concluding remarks are presented in Sec.VIII.

A. Model and notations
We consider a polymer represented as an inextensible continuous curve of contour length L in a three-dimensional space.We call û(s) the local unit tangent vector, s being the curvilinear coordinate along the chain (with 0 s L).The energy of a configuration is assumed to be that of a thin elastic rod resisting bending, Here, the prime denotes differentiation with respect to s, and κ is the bending elastic modulus, related to the (bend) persistence length l p by κ = l p k B T , with k B the Boltzmann constant and T the temperature.We use torque-free conditions at the ends by imposing that all allowed configurations satisfy (∂ s û) s={0,L} = 0. Note that we do not include a torsional stiffness term in the expression (1), but we check in Sec.VII that torsion does not influence the value of the distribution of end-to-end vectors, and it is thus ignored here.We focus on the calculation of the probability distribution function (PDF) p(r) of the end-to-end vector r = (r x ,r y ,r z ) = L 0 û(s)ds.We also introduce the probability η(r c ) that the distance between the end monomers is smaller than a capture radius r c , with η(r c ) = |r| r c p(r)dr.To be clear with the notations, we specify that dr = dr x dr y dr z and that p(r)dr is the probability that each spatial coordinate of r is observed in the interval ]r i ,r i + dr i [.Note that, for symmetry reasons, p(r) does not depend on the orientation of r and is thus related to the PDF q(r) of the end-to-end distance r ≡ |r| by p(r) = q(r)/(4πr 2 ).

B. Results for the distribution of end-to-end vectors of stiff chains
The aim of the present paper is to generalize the formula derived by Shimada and Yamakawa (SY) [27] for stiff chains (L l p ) and vanishing end-to-end distances: L 5 e −14.055 l p /L+0.246L/ l p . (2) Here, p is exponentially small with 14.055l p /L, which represents the minimal bending energy (in units of k B T ) required to form an end-to-end contact, and the term in front of the exponential, called the preexponential factor, comes from the integration over fluctuations.The SY expression (2) and the ring closure probability derived in the same paper [27] have been used to characterize the closure probability of  3) of p(r) (continuous lines).We also show asymptotic values of h,E * (dashed lines) showing the compatibility of the results with those of Refs.[27,43].
The main result of this work is the following formula, which generalizes Eq. ( 2) for nonvanishing end-to-end distances: where E * is the dimensionless minimal bending energy to form a loop of size r ≡ |r|, λ 0 = −(E * ) (0) 21.55, and h is a dimensionless function describing the variation of the entropic prefactor with r ≡ r/L (explicit expressions of E * and h are given below).The entropic prefactor in the above expression is asymptotically exact at leading order in the stiff limit (l p L), when the ratio r/L does not tend to unity at the same time (more precisely, in both cases when r/L or rl p /L 2 is kept constant).The above expression is therefore relevant to describe the end-to-end encounter probabilities in the closure problem of stiff chains, for which r is typically not close to the contour length.Its validity can be controlled by looking at Fig. 2, which shows that our analytical expression accurately predicts the values of p for parameters where it varies by eight orders of magnitude.On the other hand, our formula requires that L − r L 2 /l p to be valid, and it fails to describe p(r) for stiff chains near full extension, for which the reduction in the number of available configurations (which even vanishes for r L) must be taken into account.This regime is described by the analytical theory of Ref. [43], and we show below that our Eq.(3) shares one validity regime with this theory (see Sec. II D).Note also that in (3) we have included a next-to-leading-order correction term e αL/ l p with α = 0.246 as in the SY formula; the rigorous calculation of these corrections for any capture radius is beyond the scope of this work, but α was numerically found to remain in the interval [0.2; 0.25] for r < 0.8L and does not significantly vary with r.
The functions E * and h are represented in Fig. 1 and can be expressed explicitly in terms of a parameter m defined by where K(m) and E(m) are, respectively, the complete elliptic integrals of the first and second kind, with elliptic parameter m.For finite r, this parameter m is related to the most probable angle between r and u(0) at closure, denoted π − ϕ 0 * (see Fig. 3 FIG.2. Distribution p(r) of end-to-end vectors r calculated in this work [Eq.(3), black continuous lines], compared to its exact value obtained numerically in Ref. [55] (green symbols).We also show the results of Shimada and Yamakawa (SY) [27] for vanishing r (black squares) and those of Wilhelm and Frey (WF) [43] valid near full extension (dashed red lines).Inset: exact p(r) (for the same values of L/ l p , black continuous lines) compared to a naive estimate that assumes that the variations of p come only from the variations of the minimal bending energy to form a loop of size r, E b (r) ≡ κE * /L (dashed blue lines).
for angular convention), by and explicit expressions of E * and h are where The values of the optimal angle (5) and of the rescaled lowest bending energy (6) are well known [34,35].The main contribution of the present paper is the identification of the leading-order behavior of the preexponential factor for a finite capture radius in Eqs. ( 3) and (7).Our formula reduces exactly to the SY result for vanishing end-to-end distances since the equality 2λh = 112.04holds for r = 0 (Fig. 1).We also show that the probability η(r c ) that the end-to-end distance is smaller than r c is well approximated by

C. Scalings and geometric interpretation
The expression for p can be recast in the form FIG. 3. Examples of configurations that minimize the bending energy when the end-to-end vector is fixed to r.In (a) r is finite, and such configurations are defined up to a single rotation of angle α, while in (b) r vanishes, and the looping configurations are defined up to three angles, denoted α,μ,ω.The initial angle ϕ 0 * is also indicated.
defining a preexponential factor K * .A remarkable prediction of Eq. ( 3) is the existence of two distinct scaling regimes for K * , namely The existence of these distinct scaling regimes can be expected by considering the following geometric argument, illustrated in Fig. 3. Consider the configurations of minimal bending energy that satisfy the constraint that the three-dimensional end-to-end vector is fixed to a value r.Then, for vanishing r, one realizes that these configurations are defined up to three angular degrees of freedom [Fig.3(b)] that describe the chain's orientation at closure.The situation is different when r is finite, in which case the lowest energy configurations that have a fixed value of r are defined up to one rotational degree of freedom [Fig.3(a)].Consequently, there is a strong variation with the end-to-end radius of the "number" of lowest energy configurations with an end-to-end vector inside a fixed infinitesimal volume dV centered around r.One can thus expect that the entropic prefactor K * behaves differently in the regimes of small or finite capture radius r.The calculation presented below enables us to identify the crossover length scale between these regimes to be L 2 /l p and to quantify the value of K * .This crossover between two regimes also appears in the expression of η: since G(u → 0) ∼ u 2 and G(∞) = 1, we see that the preexponential factor in Eq. ( 9) scales as r 3 c l 2 p /L 5 for vanishing r c , whereas for larger r c the prefactor scales as r c /L.Thus, the distinct scaling regimes for the prefactor, present in the expression of p, also appear in the expressions of η and are due to the variation of the chain orientational fluctuations at closure.Equation (12) means that the entropic prefactor K * strongly varies with the end-to-end radius.To illustrate this effect, we represent in the inset of Fig. 2 the exact value of p(r) computed in Ref. [55] compared to the naive guess p = p SY e −[E b (r)−E b (0)]/k B T , which takes only into account the variation of the bending energy barrier to form a loop [called E b (r)] and not the variation with r of the entropic prefactor.Contrary to our formula (3), this estimate is clearly inaccurate for finite r by several orders of magnitude (Fig. 2, inset), demonstrating the importance of taking into account entropic effects in the calculation of p.

D. Link with existing results for nearly extended chains
We observe in Fig. 2 that our theory is not valid anymore when r is close to L: the reduction of the number of available configurations occurring when r → L is not taken into account by our approach.The consequence of this effect is a strong variation of p when L − r is comparable to ∼L 2 /l p .Our theory is not valid at these length scales, and in particular it cannot describe the value of p at its maximum.This regime is described by a result obtained by Wilhelm and Frey (WF) [43], who performed an expansion around fully extended configurations and showed that with Here the factor N is used to impose that p WF is a normalized PDF.We present in Appendix C a derivation of this result with our formalism, and we show there that the Laplace transform of H takes a simple form, and that the formula is valid up to corrections of order O(e −L/ l p ), which are smaller than 0.5% for l p L/2.The WF result ( 13) can be explicitly linked to our theory by noting that in the joint limit (L − r)l p L 2 and l p /L 1 it reads and we see that this expression coincides with our expression (3) when r → L, since one can check that h(1) = π/2, and The above expression is valid when one can neglect the terms proportional to the second derivative of E * in Eq. ( 3), leading to the validity regime (L/ l p ) 1/2 (L − r)/L L/ l p .Our expression and the WF result have therefore one validity regime in common, as illustrated in Fig. 2, providing a supplementary consistency test of our theoretical approach.

III. LOOPING CONFIGURATIONS OF MINIMAL BENDING ENERGY
The configurations of minimal bending energy that satisfy the constraints on the end-to-end vector r = (r x ,r y ,r z ) = a êx are well known; see, e.g., Ref. [65] for a = 0 and Refs.[34,35] for a > 0.Here we briefly describe these configurations.We introduce the spherical angles θ (s),ϕ(s) describing the orientation of û(s), such that û = (sin θ cos ϕ, sin θ sin ϕ, cos θ) (16) in a fixed Cartesian reference frame (ê x ,ê y ,ê z ).The energy of the chain [Eq.( 1)] expressed in terms of these angles reads The torque-free conditions at the ends impose that all allowed configurations satisfy θ | s=0,L = ϕ | s=0,L = 0. From now on, we set the units of length and energy such that L = 1 and k B T = 1 and we consider the stiff limit κ → ∞, equivalent to a weak-noise limit.
To describe the configurations that minimize E under the constraint that r = a êx , we introduce a Lagrange multiplier λ associated with the constraint r x = a, and we define Note that we could add two supplementary Lagrange multipliers associated with the constraints r y = r z = 0, but one can check that they vanish and therefore do need to be included.The lowest energy looping configurations are planar; for simplicity, we describe the optimal configurations lying in the horizontal plane, with θ = π/2,ϕ = ϕ * (s).The condition δF/δϕ = 0 leads to with ϕ * (0) = ϕ * (1) = 0.This equation can be integrated once, leading to the first-order nonlinear differential equation where ϕ 0 * = ϕ * (0) is the angle that the polymer extremity makes with the direction êx in the optimal configuration (see Fig. 3).

IV. DISTRIBUTION OF END-TO-END VECTORS FOR FINITE CONTACT RADIUS
In this section, we consider the distribution of end-to-end vectors obtained in the stiff limit, l p /L → ∞, by keeping the parameter r/L constant (that is, in our units, the limit κ → ∞ when r = a is constant).When a is finite, the lowest energy configurations are degenerate and defined up to one rotational degree of freedom [see Fig. 3(a)].Selecting one optimal looping configuration requires us to specify another variable; here we consider the quantity z, representing a weighted average vertical extension of the chain, defined as where we use the notation f |g = 1 0 ds f (s)g(s) for any functions f,g.At this stage, the choice of the weight function sin ϕ * seems arbitrary, but it will become clearer below; other choices would lead to more complicated calculations (with the same final result).We consider the joint probability distribution function (PDF) of r and z, denoted Q(r,z), and we will see later how p can be deduced from Q.By definition, (25) where Z is a normalization constant that will be calculated below, and D[cos θ ] D[ϕ] is understood as being proportional to the product n i=1 sin θ i dθ i dϕ i when one decomposes [0; 1] into a large number of n subintervals.Here, the notation δ 3 is used for the three-dimensional Dirac delta function, All chains that satisfy both conditions r = a êx and z = 0 are necessarily around one of the two lowest energy configurations lying in the horizontal plane.Therefore, the integral ( 25) is dominated by the contribution of the configurations of the form θ = π/2 + θ 1 and ϕ = ϕ * + ϕ 1 , where θ 1 and ϕ 1 represent small deviations.Replacing the energy E by F [Eq. ( 18)] expanded at quadratic order, and expanding at linear order the terms inside the δ functions in Eq. ( 25), we obtain × δ( sin ϕ * |ϕ 1 )δ( cos ϕ * |ϕ 1 )δ( 1|θ 1 ) with The factor 2 in Eq. ( 26) comes from the fact that the fluctuations around −ϕ * (s) or +ϕ * (s) give the same contribution to Q.The operators −∂ 2 s + R θ and −∂ 2 s + R ϕ can be called fluctuation operators [34,67].It is important to note that the operator −∂ 2 s + R θ has one vanishing eigenvalue, associated with the eigenfunction sin ϕ * (s) [which is why we used this particular function in the definition (24) of z].Physically, the presence of the vanishing eigenvalue is linked to the fact that rotating the chain by a small angle δα around the axis êx changes neither the bending energy nor the end-to-end vector, and this transformation is Evaluating (26) now amounts to performing integrals involving multivariate Gaussian distributions.In Appendix A, we recall some formulas that give the results of such integrations, including the case of vanishing eigenvalues.Using these formulas, Q can be expressed in terms of Green's functions and functional determinants, namely where det (L) represents the determinant after extraction of the zero mode (thus the product of all nonvanishing eigenvalues of L).Note that we use the notation for any functions f,g and symmetric A. In Eq. ( 28), G ϕ (s,s ) is the Green's function of the Sturm-Liouville operator −∂ 2 s + R ϕ ; it can be expressed as where y ϕ (s) is the solution of The function y ϕ is therefore a homogeneous solution of the operator [−∂ 2 s + R ϕ (s)], but it satisfies only one-half of the Neumann boundary conditions.
We now calculate the normalization constant Z.We consider the average angles φ = 1 0 ds ϕ(s) = 1|ϕ and θ = 1|θ , and we denote as q( θ, φ) their joint PDF for free chains.In the stiff limit, polymers in solution are in nearly straight configurations where the angles θ,ϕ are almost equal to their average value θ, φ.We can therefore pose ϕ φ + ϕ 1 and θ = θ + θ 1 , and we expand the energy (18) at quadratic order to get Noting that the uniform function is an eigenfunction of the operator −∂ 2 s with vanishing eigenvalue, we apply again the formulas of Gaussian integration (see Appendix A) and find Isotropy imposes that we also have Comparing the above formulas (for θ = π/2) leads to the following expression for the normalization constant: Next, the calculation of functional determinants such as those appearing in Eqs. ( 28) and ( 35) has a long history [68,69], which is not discussed here.Applying the general formulas of Refs.[70,71] where y ϕ is defined in Eq. (31).Note that the factor 2π/κ in the above equation differs from formula (39) of Ref. [70] because of a slightly different definition of the determinant without the zero mode.
A geometrical argument can be used to relate Q and p.We call p z(0|r = aê x )dz the probability that z ∈ [0,dz] given that r = a êx .Bayes' theorem leads to We call α the angle between an optimal looping configuration and the horizontal plane; see Fig. 3(a).For such an optimal configuration, we note from elementary trigonometry that cos θ = û • êz = sin α sin ϕ * , and thus z = sin α 1 0 ds(sin ϕ * ) 2 .Since α is uniformly distributed on [0,2π[, we have where the factor 2 in the numerator comes from the fact that z = 0 is obtained either when α = 0 or when α = π .The above relations imply that Collecting all the above results [Eqs.( 28), ( 35), (36), and (37)] and inserting them into Eq.( 40), we finally obtain where we have defined the dimensionless function h as Reestablishing (temporarily) homogeneity, we obtain We have thus identified how the preexponential factor depends on the properties of the lowest energy configuration (the initial angle ϕ 0 * and the dimensionless force λ = −∂ r E * ), and on the fluctuation operators linearized around it (the function y φ ).
We now show that this result can be simplified by noting that the function y ϕ can be calculated explicitly.Let us consider the functions It can be checked that y 1 and y 2 are independent solutions of the second-order differential equation −y + R ϕ (s)y = 0. We can thus express y ϕ as a linear combination of y 1 ,y 2 : where a 1 ,a 2 are identified by using the boundary conditions y ϕ (0) = 1,y ϕ (0) = 0.After a few algebraic manipulations, described in Appendix B, we obtain analytical expressions for a 1 ,a 2 and therefore G ϕ , and the formula for h finally simplifies to Let us now discuss the validity regime of Eq. ( 43).In its derivation, we assumed that all chains that satisfy the conditions z = 0 and r = aê x are close to the configuration of minimal energy satisfying these two conditions.This requires κ 1.Additionally, in Appendix D, we show that the conditional variance of ϕ 1 (s) (when r,z are fixed) varies as 1/(κr) for small r.The assumed hypothesis of small fluctuations is therefore correct if r 1/κ (or, equivalently, r L 2 /l p with dimensions).Equation ( 43) is therefore not valid for infinitely small end-to-end distance, since in this limit the fluctuations around the minimal energy configuration diverge, leading to an (apparent) divergence of the prefactor in (43), which is similar to the divergence within the WKB (Wentzel-Kramers-Brillouin) approximation of the probability amplitude (in quantum mechanics) or the light intensity (in optics) near caustics [69].This divergence is due to the presence of an infinity of classical paths (or optical rays) at caustics, and it can therefore be compared to the explosion of the number of lowest energy looping configurations when the capture radius r becomes small (Fig. 3).A calculation of p(r) valid for vanishing r is presented in Sec.V.

V. DISTRIBUTION OF END-TO-END VECTORS FOR SMALL CONTACT RADIUS
In this section, we consider the distribution of end-to-end vectors for small r (in the sense that r 1/κ 1/2 ; see below).The distribution of end-to-end vectors for arbitrarily small r is not given by Eq. ( 43), since the prefactor in this expression diverges for r → 0. The major change arising in this limit is the degeneracy of the optimal configurations, which are now defined up to three rotational angles denoted = {α,μ,ω} [see Fig. 3(b)] instead of only one.We call 0 the set of angles such that the looping configuration is in the horizontal plane, described by θ = π/2,ϕ = ϕ * (s), as represented by the blue curve in Fig. 3

(b). Let us set λ
Interpreting κλ 0 as an effective force, we see that there is an energy reduction κλ 0 r x upon opening the loop by a distance r x in the direction êx when is fixed to 0 .Therefore, the probability density of r given that 0 takes the form where C 0 does not depend on r.When writing the above expression, we implicitly assumed that the variation of energy with r is well approximated by using the first derivative of E * only.This approximation holds if e κ(E * ) (0)r 2 x /2 1, which is realized when r 1/ √ κ.In this section, we assume that this condition is satisfied.
Let us call q(r| = 0 ) the PDF of the radius r = r given that = 0 , which can be obtained by integrating Eq. ( 47) over the orientation of r, Now, due to isotropy, q(r| = 0 ) does not depend on 0 and is simply equal to the radial distribution q(r).Performing the integral in the above expression leads to The distribution of end-to-end vectors p(r) does not depend on the orientation of r, and it can be deduced from the radial distribution: The constant C 0 will be identified by requiring that the solution (50) matches the solution (41) obtained for nonvanishing r.If we write Eq. ( 41) in the limit r 1/ √ κ, we obtain We can also consider Eq. ( 50) in the regime r 1/κ: The above expressions ( 51) and ( 52) should be identical since each of them is valid for κ −1 r κ −1/2 .They can be matched if we set the value of C 0 to Inserting this result into Eq.( 50), we thus find that the distribution of end-to-end vectors in the limit of vanishing radius (in the sense that r Reestablishing homogeneity, we obtain The above expression is valid when r L 3/2 /l 1/2 p .It has in common with Eq. ( 43) the validity regime L 2 /l p r L 3/2 /l 1/2 p in which both expressions reduce to Combining Eqs. ( 43) and ( 55), one can build the solution (3), which provides a correct estimate of the distribution of end-to-end vectors in all the regimes.The main result of this section is the functional form (55), which describes the transition between the regimes in which the optimal configurations of fixed r are defined up to, respectively, one or three angular degrees of freedom.We can check this function by representing f = K * (r)/K * (0) [where K * is the preexponential factor, defined by p(r) = K * (r)e −κE * (r) ] as a function of u ≡ λ 0 rl p /L 2 .Our theory predicts that the obtained curves should converge in the stiff limit to give This prediction is consistent with what is observed in Fig. 4.
Here it is interesting to note that it has already been observed in Ref. [63] that the preexponential factor should differ from a constant for small values of r; there it was proposed that K * (r) [see Eq. ( 19) in Ref. [63], in which we have used λ 0 r [E * (−r) − E * (r)]/2, valid for small r/L].Here I 0 is a modified Bessel function of the first kind.We see in Fig. 4 that the above functional form captures the behavior of K * (r) for small rκ, but it differs for larger values, indicating that it (1 e -2u )/(2u) FIG. 4. Value of the ratio of preexponential factors f = K * (r)/K * (0) [defined by p(r) = K * (r)e −κE * (r) ] as a function of u ≡ λ 0 l p r/L 2 .Symbols are data from Ref. [55], dashed lines with the same colors are the theoretical predictions of (3).The theoretical small noise limit f = (1 − e −2u )/(2u) is represented [thick black line, see Eq. ( 55)], as well as the form f = I 0 (u)e −u proposed in Ref. [63] (red dash-dotted line).
is not an exact expression.We conclude that the behavior of p(r) for r ∼ L 2 /l p is asymptotically given by ( 55) in the limit of stiff chains.

VI. RADIUS-DEPENDENT CLOSURE FACTOR AND CUMULATIVE END-TO-END DISTRIBUTION FUNCTION
We consider here the probability η(r c ) that the end-to-end distance r is smaller than a capture radius r c , [where the last equality holds because p(r) does not depend on the orientation of r].We also derive formulas for the so-called J factor [72], defined as the ratio of equilibrium reaction constants of cyclization and dimerization.J can also be interpreted as an effective concentration of end monomers with respect to the monomers at the other extremity.If looping occurs with equal probability for all r satisfying r < r c , then [35] Here we give explicit formulas for η(r c ); corresponding expressions for J are straightforwardly obtained from Eq. ( 60).First, for a small capture radius (meaning that r c L 3/2 /l 1/2 p ), using Eq. ( 55) we get where we have introduced which has the properties Note that, for r c L 2 /l p , one has Now, in the regime r c L 2 /l p , the values of r that most contribute to the integral (59) are those that are close to r c , hence we can perform this integral by using the value of p given by Eq. ( 43), where we use a linear approximation of E * for r r c : where we recall that λ = −∂ r E * (r).Performing the above integral leads to which is valid when r c L 2 /l p , still in the stiff limit l p L. The expressions ( 61) and ( 67) have in common the validity regime L 2 /l p r c L 3/2 /l 1/2 p .Combining the above formulas, and multiplying by the additional factor exp(0.246L/ l p ) (which takes approximately into account the corrections at next order in L/ l p ), we obtain where ), Eq. (63).The above expression is asymptotically exact at leading order in the limit l p /L → ∞ when either rl p /L 2 or r/L is kept constant, and the next-order corrections are approximately taken into account.
We see that the preexponential factors in Eqs. ( 65) and ( 67) are different.This change of scaling comes from the difference of the magnitude of the orientational fluctuations of the chain at closure.In Fig. 5, we check that our approximate expression (68) can accurately reproduce the results of simulations performed in Ref. [48].

VII. EFFECT OF TORSIONAL STIFFNESS ON CLOSURE PROBABILITY
In this section, we consider the effect of torsional stiffness C on the value of p(r).We associate with each coordinate s a local orthonormal basis (ê 1 (s),ê 2 (s),ê 3 (s)), with ê3 = û, and a rotation vector w(s) such that ê i = w × êi .For simplicity, we restrict the discussion to the model called KP-1 in Refs.[73,74], where the elastic energy is given by where C is the torsional rigidity.We denote by (ê r ,ê θ ,ê ϕ ) the unit basis vectors of the spherical coordinate system in an external fixed reference frame.Let ψ be the angle between êθ and ê1 .The elastic energy can now be written [73,74] Simulations Theory [Eq.( 9)] FIG. 5. Probability η(r c ) of observing the end monomers closer than r c : theoretical expression (9) (continuous lines), simulations of Ref. [48] (black circles).We also show the result η (4πr 3 c /3)p SY (dotted lines) obtained by assuming that the end-to-end encounter probability density p(r) does not depend on r.Here, l p = 50 nm.
where E is the bending energy.The suffix T in E T (and other quantities) indicates that torsion is included.The torque-free boundary conditions for the supplementary field ψ are Neumann conditions ∂ s ψ| s={0,L} = 0. Consider now the functional derivatives We see that imposing ∂ s ψ + cos θ∂ s ϕ = 0 removes all the terms proportional to C in Eqs. ( 71) and ( 72) and also enables us to satisfy δE T /δψ = 0.The lower-energy configurations are thus not modified by the torsional stiffness.
Next, let us evaluate the energy functional at quadratic order around a (horizontal) planar optimal configuration with θ where ) is the quadratic expansion of the energy in the absence of torsion.In the presence of torsion, the joint PDF Q T of r and z becomes where Z T is the partition function of free chains with torsional stiffness, For free chains, at quadratic order we get 2 and therefore Using this expression and the quadratic expansion (74), we obtain by setting The joint PDF Q(r; z) is therefore not modified by the presence of torsional stiffness.Since Eqs. ( 38) and (39), which relate Q and p(r), are not modified by the presence of torsion, we conclude that the value of the torsional stiffness has no influence on the closure probability, at least in the limit of stiff chains (l p L) with finite (or vanishing) r/L.This generalizes the same result obtained for vanishing end-to-end distances in Ref. [27].

VIII. CONCLUDING REMARKS
In conclusion, we have performed, in the limit of stiff chains, a calculation of the equilibrium probability density p(r) of the end-to-end vector r of a stiff wormlike chain.Up to now, this quantity was analytically known only for vanishing endto-end distances [27] or for nearly extended chains [43].Here, the closure factor is analytically calculated for any value of the end-to-end distance that is not too close to the contour length, as relevant in the problem of DNA closure assisted by flexible protein bridges or through the pairing of single-stranded DNA.The main result is the existence of two distinct scaling regimes for the preexponential factor, reflecting the fact that the fluctuations of the chain orientation at closure are highly dependent on the value of r.Thus, the strong variation of the preexponential factor is physically associated with a change in the "number" of lowest energy looping configurations occurring when r increases.The apparent divergence of the preexponential factor of p for small r in Eq. ( 43) is similar to the divergence within the WKB approximation of the probability amplitude (in quantum mechanics) or the light intensity (in optics) near caustics [69].This divergence is due to the presence of an infinity of classical paths (or optical rays) at caustics, and it can therefore be compared to the explosion of the number of lowest energy looping configurations when the capture radius r becomes small.
The presence of two distinct scaling regimes for the preexponential factor is not limited to the simple version of the wormlike chain model considered here; since it originates from geometrical considerations, it should also appear in more complex models of semiflexible polymers.The present approach could also be adapted to more realistic model geometries, for example in the case when one of the monomers is attached to a surface, as frequently occurs in experiments.
In such a situation, the surface modifies the "number" of available lowest energy looping configurations, and it would be interesting to see if our calculations could be adapted to quantify the modification of the closure factor that was recently observed in numerical simulations for polymers near a surface [47].It has also been noted recently [77] that the model proposed in Ref. [78], which includes a twist-bend coupling term in the chain elastic energy, could provide a better description of the results of recent single molecule experiments.Here we have shown that torsional stiffness does not modify the end-to-end PDF in the weak noise limit, due to the decoupling between twist and bend angles.Such decoupling could disappear in the presence of twist bend coupling, and it would be interesting to see how to adapt our calculations of the looping probability to this case.and it is useful to generalize it in the presence of vanishing eigenvalues.Consider now the case in which c is an eigenvector of L, associated with an eigenvalue μ, which tends to 0. In this case, L −1 • c = c/μ, and using Eq.(A5) for n = 1 in the limit μ → 0 gives where we denote det (L) the determinant after extraction of the zero mode (i.e., the product of all nonvanishing eigenvalues of L).Next, considering a second vector b orthogonal to c (which is still associated with an eigenvalue μ → 0), we obtain from Eq. (A5) for n = 2 where f is the vector that satisfies L • f = b with the condition t c • f = 0.In the main text, we apply these formulas in the continuous case, with the inverse matrices and the discrete determinants replaced by Green's functions and functional determinants, respectively.

APPENDIX B: DETAILS ON THE CALCULATION OF THE FUNCTION h
Here we describe the derivation of Eq. ( 46).First, the coefficients a 1 ,a 2 appearing in Eq. ( 45) are identified by using the boundary conditions y ϕ (0) = 1,y ϕ (0) = 0. Noting that y 1 (0) = 0 and y 2 (0) = −1/y 1 (0), we obtain and The value of y 2 (0) is found by inserting the known expression (21) of ϕ * into (44), leading to We now consider the constant M defined by Using the value (30) of G ϕ , taking into account the symmetry properties y 1 (1 − s) = y 1 (s) and y 2 (1 − s) = −y 2 (s), and noting that sin ϕ * = −y 1 /λ, we obtain We write M as where we calculate the terms K ij separately.First, where we have used the fact that y 1 vanishes at the boundaries.Using the same reasoning, we have where we have used the definition (44) of y 2 in the last equality.
The above expression can be simplified by integration by parts: Next, we consider Integrating by parts yields Using again the definition (44) of y 2 , one finds and one can again use integration by parts to find The last term is Integrating the last term yields Comparing with (B11), we find All the terms a i ,K ij have now been calculated.Inserting them into (B6) and using

APPENDIX C: EXPANSION NEAR FULLY EXTENDED CONFIGURATIONS
Here we briefly show how the results of Ref. [43] can be obtained with our formalism.The unit of length is again such that L = 1.We investigate the behavior of p near fully extended configurations: we set r = 1 − and define (C1) For small enough , the configurations that contribute to P ( ) are nearly extended ones, with θ(s) = π/2 + θ 1 (s) and ϕ(s) = ϕ 1 (s), θ 1 ,ϕ 1 being small deviations.At leading order, the coordinates of the end-to-end vector r = (1 − ,0,0) are 0 = r y = Note that these expansions are valid up to quadratic order in θ 1 ,ϕ 1 , which is crucial for r x , in which linear terms vanish.Using the above approximations for r x ,r y ,r z , and using a quadratic expansion of F near fully extended configurations, P ( ) reads Applying the formulas of Gaussian integration (A5), and replacing Z by its value (35), we obtain det −κ∂ 2 s (2π ) det −κ∂ 2 s + q (2π ) where G is the Green's function of −∂ 2 s + q/κ.The quantities appearing in this formula can be calculated if one introduces the function y(s) defined as −y (s) + (q/κ)y(s) = 0, y(0) = 1, y (0) = 0, (C8) whose solution is y(s) = cosh(s q/κ). (C9) The functional determinants and the Green's function appearing in (C7) are calculated by using Eqs.( 37) and (30) (where y ϕ is replaced by y), leading to P (q) = √ q/κ 4π sinh( √ q/κ) .(C10) The Laplace transform of P ( ) thus has a simple expression.This Laplace transform can be inverted by using the Bromwich-Mellin formula.With a residue calculation, we find p[(1 − )e x ] = P ( ) κ N H ( κ), (C11) with H (x) = ∞ n=1 (−1) n+1 n 2 e −xn 2 π 2 π/2, which is the result of Ref. [43].Here a normalization factor N has been introduced to enforce the normalization of (C11).The expression of N is thus where we have used the variable u = κ.Replacing the upper bound of the integral 1/κ by infinity is equivalent to neglecting terms of order e −1/κ .With this approximation, using the smallq expansion of (C10) yields The difference between the above expression and the numerical evaluation of Eq. (C12) is less than 0.5% for l p L/2.We finally note that we can obtain the behavior of P for → 0 by noting that, for large q, H (q) e − √ q √ q/(2π), which can be inverted to H (u) e −1/(4u) (1 − 2u)/(8π 2 u 5/2 ), (C14) which is the first term of the series appearing in Eq. (4) of Ref. [43].

APPENDIX D: FLUCTUATIONS NEAR A MINIMAL ENERGY CONFIGURATION
In the derivation of Eq. ( 41), one assumes that, when r and z are fixed, the fluctuations of the chain around the configuration of minimal energy are small.Here we compute the magnitude of these fluctuations.Let us call ([ϕ 1 ,θ 1 ]) the distribution of ϕ 1 ,θ 1 conditional to z = 0 and r = 0, which factorizes into × δ( sin ϕ * |ϕ 1 )δ( cos ϕ * |ϕ 1 ).(D1) We call ϕ 1 (s) 2 the variance of ϕ 1 over the distribution ϕ .
Using formulas for the conditional covariances of multivariate Gaussian distributions [76], one finds

FIG. 1 .
FIG.1.Dimensionless functions E * (a) and h (b) appearing in the expression (3) of p(r) (continuous lines).We also show asymptotic values of h,E * (dashed lines) showing the compatibility of the results with those of Refs.[27,43].