Fluctuations of ring polymers

We present an exact solution for the distribution of sample averaged monomer to monomer distance of ring polymers. For non-interacting and weakly-interacting models these distributions correspond to the distribution of the area under the reflected Bessel bridge and the Bessel excursion respectively, and are shown to be identical in dimension d greater or equal 2. A symmetry of the problem reveals that dimension d and 4 minus d are equivalent, thus the celebrated Airy distribution describing the areal distribution of the one dimensional Brownian excursion describes also a polymer in three dimensions. For a self-avoiding polymer in dimension d we find numerically that the fluctuations of the scaled averaged distance are nearly identical in dimensions 2 and 3, and are well described to a first approximation by the non-interacting excursion model in dimension 5.

The statistical mechanics of polymers has been well studied for many years due both to the numerous practical applications of polymers as well as their many interesting properties. One signal finding is that the overall size of a polymer of length N scales as N νp , where ν p is a dimension dependent critical exponent. This is reflected in the behavior of various observables, such as the average end-to-end distance, and the radius of gyration (R g ) [1,2]. The scaling exponent is known to be sensitive to the excluded volume interaction applied to part or all of the monomers. Other geometrical constraints applied to the chain, such as cyclization (where the first monomer is connected to the last one), leading to a "ring" polymer, affect only the prefactor for these quantities. Ring polymers are commonly found in many biological systems e.g., bacterial and mitochondrial genomes, as well as DNA plasmids used in many molecular biology experiments [3]. Recently, ring polymers were also studied in the context of a model for chromosome territories in the nucleus of eukaryotic cells [4].
The conformational fluctuations of some polymer models can be analyzed using the theory of random walks (RW) [1,2]. In particular, the fluctuations of the polymer size are of physical and biological importance. In the current paper we study the distribution of sample-averaged monomer-to-monomer distance of ring polymers, both for ideal, noninteracting polymers as well as for polymers with excluded volume. In particular, we exploit recent mathematical development on d dimensional constrained Brownian motion (defined below) [5] to find an exact expression for the distribution of sample averaged monomer to monomer distance for both ideal ring polymers and those with an additional excluded volume interaction applied at a single point. This observable yields insight on the sample averaged fluctuations of polymer sizes. An important ingredient of this calculation is the identification of the appropriate boundary conditions for the underlying equation, a variant of the Feynman-Kac equation, which depends on both the dimensionality and the interaction, and in turn yields a selection rule for the solution. The resulting distributions are then compared numerically to those measured in simulations for a ring polymer with full excluded volume constraints.
As we have noted, constrained random walks lie at the heart of our analysis. For rings, the primary constraint is that the path returns to the origin after N steps. Statistics of such constrained one dimensional Brownian paths have been the subject of much mathematical research [6][7][8][9]. These constrained paths have been given various names, depending on the additional constraints imposed. The basic case is that of a Brownian bridge where the return to the origin is the sole constraint. For a Brownian excursion, the path is also forbidden from reaching the origin in between. Majumdar and Comtet used Brownian excursions to determine statistical properties of the fluctuating Edwards Wilkinson interface in an interval [10,11]. The focus of most previous work has been on the constrained one dimensional Brownian paths which describes inherently non-interacting systems (note that the problems of non-intersecting Brownian excursions [12] or vicious random walkers [13] in higher dimensions are exceptions). For the case of the fluctuations of ring polymers, we need to extend the theory of Brownian excursions and bridges to other dimensions. We address the influence of different kinds of interactions on the polymer structure, both analytically (for a single point interaction) and numerically (for a polymer with excluded volume interactions). These models yield rich physical behaviors and open new questions.
Polymer Models. We consider three lattice models of ring polymers with N bonds, each of length b, in d dimensions. The simplest polymer model is an "ideal ring" -a closed chain without excluded volume, where different monomers can occupy the same lattice site. While such a polymer does not exist in nature, its global behavior is the same as that of a polymer at the Flory θ-temperature [1]. An ideal ring chain corresponds exactly to an unbiased RW in d dimensions which starts and ends at the origin, i.e. a d-dimensional bridge. In the second model, the "weakly interacting ring polymer", the first and last monomers are tied to the origin and no other monomer is allowed to occupy this lattice site. This case is equivalent to that of d-dimensional excursions. The third model is a ring polymer with excluded volume interactions, also called a self-avoiding walk (SAW). Further details on the polymer models and simulation methods are provided in the supplementary material (SM). We first consider the ideal and weakly interacting ring models, for which we can provide an analytic solution.
Bessel process. In the analogy between statistics of an ideal polymer and a RW, the position of the ith monomer, r i , corresponds to the position r i of the random walker after i time steps and N is proportional to the total observation time. The Bessel process [14,15] describes the dynamics of the distance r = |r| from the origin of a Brownian particle in d dimensions. This process is described by the following Langevin equation: where η(t) is Gaussian white noise satisfying η = 0 and η(t)η(t ) = δ(t − t ). One may map the polymer models to the Bessel process using r(t) 2 where R 2 is the mean square end-to-end distance of an ideal linear polymer chain without constraints. In what follows we choose b 2 = d and t = N . Bessel Excursions and Reflected Bridges. The process r(t) with the additional constraint of starting and ending at the origin, is called a reflected (since r ≥ 0) Bessel bridge. This process corresponds to an ideal (noninteracting) ring chain. Bessel excursions are paths still described by Eq. (1) however with the additional constraint that any path that reaches the origin (besides t = 0 and t = N ) is excluded. The Bessel excursion corresponds to what we have called the "weakly interacting" ring chain, where a multiple occupation of the origin is not allowed. The mapping of the polymer models to the Bessel process, allows us to extract statistical properties of the former with new tools developed in the stochastic community [2,5,11].
The Observable A. For a ring polymer, let r i be the position of the ith monomer in space where i = 0..N , and we place the origin at the position of the zeroth monomer, r 0 . For the weakly interacting chain this monomer is also the excluding one. We study a new measure, A, for the size of a ring polymer, defined by In the RW language, A is the area under a random process, and hence is a random variable itself. Clearlyl = A/N is the sample averaged distance of the monomers from the origin. Specifically, let the area under the random Bessel curve r(t) be denoted by A B = t 0 r(t)dt (the subscript B is for Bessel). More generally the mapping of the processes implies that in the limit of large N the distribution of A B / A B is identical to the distribution of A/ A (orl/ l ), with the corresponding constraints.
Numerical results. In Fig. 1 we plot the probability density function (PDF) of the scaled random variable A/ A . Both results of simulations and theory are presented, however at this stage let us focus on the main features as revealed in the simulations. For dimensions d = 2, . . . , 5, there is a clear trend of narrowing of the PDF for increasing d. This trend is explained by examining Eq. (1): As d increases, the noise term becomes negligible compared with the force term, resulting in smaller fluctuations and narrower tails. Against this expected trend are the results in d = 1 for the weakly interacting model. As we will show analytically, the weakly interacting model in dimensions one and three surprisingly have the same distribution even though d = 1 has a vanishing deterministic force term in Eq. (1) while for d = 3 the force is clearly not zero. In addition, we observe that for d ≥ 3 the shape of the distribution of weakly interacting and ideal ring chains coincide, indicating that weak interactions are negligible (when N → ∞). As we shall see, this is also observed in the theory. Indeed the theory discussed below suggests that these two distributions are already identical for d = 2. However, since this is a critical dimension, due to extremely slow convergence, we don't see this behavior in the simulations. This asymptotic convergence is logarithmic (see SM) and an ε expansion shows that it is reminiscent of critical slowing down. As for the SAW polymer, we see that the fluctuations are considerably reduced compared to the other models. This is due to the fact that the number of configurations of a SAW polymer is smaller than for the other models, hence fluctuations are smaller. A striking observation is that the two and three dimensional SAW results are identical, both being equal to the simulations of the d = 5 models. We now address these observations with theory.
Functionals of Constrained Bessel Processes. Our goal is to find the PDF P (A B , t) of the functional A B = t 0 r(t)dt of the Bessel process, constrained to start and end at the origin. We show that the difference between the weakly interacting model (the Bessel excursion) and the ideal polymer (the reflected Bessel bridge) enters through the boundary condition in the Feynman-Kac type of equations describing these functionals. The choice of boundary condition turns out to be non-trivial and controls the solution. Other aspects of the solution follow the steps in [5].
It is useful to find first the Laplace transform of P (A B , t), i.e.,P (s, t) = withG| t=0 = δ(r −r 0 ) and r 0 a cutoff which is eventually taken to zero. For d = 1, the second term on the right hand side vanishes, and we get the celebrated Feynman-Kac equation corresponding to Brownian functionals [17]. The third linear term −srG stems from the choice of our observable, namely our functional A B is linear in r [18].
Since we are describing a ring polymer, the Bessel process must start and end on the origin, and so, following [10], we need to calculatẽ The denominator gives the proper normalization condition. The first step in the calculation is to perform a similarity transformation: Using Eq. (3),K t (r, s|r 0 ) is the imaginary time propagator of a Schrödingier operator: with the effective Hamiltonian: The is constructed [5] from the eigenfunctions φ k ofĤ where λ k is the kth eigenvalue and the normalization condition is ∞ 0 φ 2 k (r)dr = 1. The subtle point in the analysis is the assignment of the appropriate boundary condition corresponding to the underlying polymer models we consider. The eigenfunctions at small r exhibit one of two behaviors: From the normalization condition, the φ − solution cannot be valid for d ≥ 4 and d ≤ 0. For the critical dimension d = 2 the two solutions are: We now solve the problem for the two boundary conditions and then show how to choose the relevant one for the physical models under investigation.
The distribution of A/ A . Following the Feynman-Kac formalism described above and performing the inverse Laplace transform [5], we find two solutions for the PDF of the scaled variable χ ≡ A/ A The solution is independent of N and valid in the limit of N → ∞. Here, |α| = |d − 2|/2, |ν| = 2|α|/3, and 2 F 2 (·) refers to the generalized hypergeometric functions.
The supplementary material provides a list of λ k and d k values for d = 1, . . . 4. For d = 1 the solution agrees with the known results [10,11,20,21], where the + solution is the celebrated Airy distribution [10,11]. The average of A is The + solution was previously presented in a slightly different form in [5] and here the question is how to choose the solution for the corresponding polymer models. Clearly, for d = 2, A + = A − , indicating that this is a critical dimension. Further A + in 1 and 3 dimensions are identical and so is A − as the result of the symmetry around d = 2 in Eq. (7). The scaling A ∝ N 3/2 is expected since r scales with the square root of N as for Brownian motion, so the integral over the random processes r scales like N 3/2 . We investigate the physical interpretation of the two possible boundary conditions. A mathematical classification of boundary conditions was provided in [15,19] and here we find the physical situations where these conditions apply. We examine the behavior of the probability current associated with the k th mode: The analysis is summarized in Table I. We see that in dimension two and higher, the current on the origin is either zero or positive. A positive current at the boundary means that probability is flowing into the system, which is an unphysical situation in our system. Hence we conclude that in dimension two and higher, the − solution is not relevant. This implies that statistics of excursion and reflected bridges (and equivalently, ideal and weakly interacting ring polymers) are identical for d ≥ 2 and correspond to the + solution.
In Fig. 1 we compare the results of the simulations of  the ideal and weakly interacting polymer models with our theoretical results for P + (A/ A ), as given in Eq. (10). As noted above, for d ≥ 3, we see that even for finite size chains the local interaction is not important, and that the theory and simulations perfectly match, while for d = 2 there are strong finite size effects in the weakly interacting case. Self-avoiding polymers. Extensive simulations of ring SAWs were performed on cubic lattices. As has already been pointed out, the global expansion of a polymer is characterized by the exponent ν p . Since A constitutes a measure of the overall size of a polymer, its behavior for large N should follow A ∼ N νp+1 . For the ideal and weakly-interacting chains ν p = 1/2, and A ∼ N 3/2 . For SAWs the exact value of the exponent depends on d, i.e., ν p = ν p (d). ν p = 1, 0.75, and 0.5 are known exactly for d = 1, 2, and 4 [22], respectively, while for d = 3, ν p , based on renormalization group considerations and Monte Carlo simulations, is ν p 0.588. These prediction were extensively tested numerically for the observable of interest A with a critical dimension of d = 4, characterized by a very slow convergence of the weakly-interacting model (see SM).
While the scaling behavior of the SAW model is different from that of the other two models (as reflected in ν p ), as we have noted, the scaled PDFs, P (A/ A ) are nevertheless similar. A striking observation is that the SAWs in d = 2 and 3 coincide to the precision of our measurements with the comparably narrow PDF of the d = 5 non-interacting model, (see Fig. 1). That these distribu-tion are narrower than the non-self-avoiding case can be qualitatively explained as follows: Since the interaction forbids many compact conformations, the fluctuations of the area become smaller. This is easily observed in the extreme case of a linear SAW in one dimension where only one conformation is allowed and the scaled PDF assumes the form of a δ-function.
Discussion. The mapping of ring polymer models to the reflected Bessel bridge and excursion is very promising since it implies that not only the observable A can be analytically computed, but also other measures of statistics of ring polymers. An example would be the maximal distance from one of the monomers to any other monomer, since that would relate to extreme value statistics of a correlated process. The famed Airy distribution describes both the one dimensional polymer, as well as the three dimensional one, due to the symmetry we have found in the underlying Hamiltonian. The case of d = 2 is critical in the sense that interaction on the origin becomes negligible, though for finite size chains it is still important. Boundary conditions of the Feynman-Kac equation were related to physical models, which allowed as to select the solutions relevant for physical models. The SAW polymer exhibits interesting behavior; the distribution of A/ A is identical (up to numerical precision) in dimension 2 and 3 and corresponds to the non-interacting models in dimension 5. Further work on this observation is required. In order to maintain the closure condition of the chain, we choose an array of length of N with N/2 components of (+1) and N/2 of (−1) in each of the directions (d such arrays), and then shuffle them for each direction separately. For each i, the components of our d dimensional step are the ith values of these arrays. The sum of all of the displacements in each direction is then naturally zero so that the last monomer is always positioned at the origin. For the ideal chain model we built 10 6 such conformations while for the weakly interacting we threw away all the conformations that crossed the origin prior to the final monomer. For each of the conformations we calculated A = N i=0 |r i − r 0 |, where r 0 = 0, and plotted the distribution of this parameter.
According to Eq. (11) in the paper, for N → ∞ we have A = c ± N 3/2 . By this we can check the convergence of the simulations to the theory as a function of N . At the critical dimension, d = 2 this convergence becomes very slow. In Fig. S1 we plot A /N 3/2 for different N values (in a logarithmic scale) in d = 2, where A = c + = 0.4922 is the theoretical value for N → ∞. One can observe the very slow logarithmic convergence to this value.

Self-Avoiding Ring Polymers
Monte Carlo (MC) simulations have been applied to self-avoiding ring polymer models [1] on square, simple cubic, and d = 4 hyper-cubic lattices. The polymer consists of N monomers (and N bonds), where the first monomer is attached to the origin of the coordinate system on the lattice, and the N th monomer is the nearest-neighbor to the origin. At step j of the MC process, monomer k (1 ≤ k ≤ N − 1) is selected at random (i.e., with probability 1/(N − 1)) and the segment of m monomers following k (i.e., k + 1, k + 2, , k + m) become subject to change in the MC process; the rest of the chain (i.e., monomers 1 to k and k + m + 1 to N ) is held fixed (notice that if k is at the end of the chain, N − m + 1 < k ≤ N − 1, m decreases correspondingly from m − 1 to 1). Thus, this current segment is temporarily removed and a scanning procedure is used to calculate all the possible segment configurations of m monomers satisfying the excluded volume interaction and the loop closure condition (i.e., the segment of m monomers should start at k and its last monomer, k + m is a nearest neighbor to monomer k + m + 1; notice that the initial segment configuration is generated as well). The segment configuration for step j is chosen at random out of the set of L configurations generated by the scanning procedure and the MC process continues.
This process starts from a given ring configuration, whose transient influence is eliminated by a long initial simulation, which leads to typical equilibrium chain configurations. Then, every certain constant amount of MC steps the current ring configuration is stored in a file to create a final sample of n rings from which the averages and fluctuations of the physical properties of interest are calculated. The segment sizes used are m = 10 for the square lattice, m = 6 for the simple cubic lattice, and m = 4 for d = 4. For each lattice several chain lengths, N are studied.
We calculated the averages of R 2 g , where and A = N −1 i=0 |r i − r 0 |. For large N , R g = R 2 g increases as ∼ N νp where ν p is a critical exponent and as discussed in the text, ν A = ν p + 1. However, we are mainly interested in the fluctuations of A, i.e., in the shape of the scaled distribution of A/ A for different N . To check the reliability of the simulations we provide in Table S1 the results obtained for R g . The ν calculated results for d = 2 and d = 3 are equal within the error bars to those of ν predicted while for d = 4 ν calculated is too large due to a logarithmic correction to scaling, which would become insignificant only for much larger N . In fact, considering this correction in the analysis has led indeed to ν calculated 0.5. The same quality of results for ν A has been obtained for the observable A (where ν A 1.5 for d = 4).
The d = 4 SAW case is a critical one, since the critical exponent, ν p for lower dimensions significantly differs from the ν p = 1/2 of the non-interacting models, and for d ≥ 4 the interactions become unimportant for N → ∞. Hence, we expect P (A/ A ) of the d = 4 SAW to coincide with that of the non-interacting model. However, for finite N the