Threading of Unconcatenated Ring Polymers at High Concentrations: Double-Folded vs Time-Equilibrated Structures

Unconcatenated ring polymers in concentrated solutions and melt are remarkably well described as double-folded conformations on randomly branched primitive trees. This picture though contrasts recent evidence for extensive intermingling between close-by rings in the form of long-lived topological constraints or threadings. Here, we employ the concept of ring minimal surface to quantify the extent of threadings in polymer solutions of the double-folded rings vs rings in equilibrated molecular dynamics computer simulations. Our results show that the double-folded ring polymers are significantly less threaded compared to their counterparts at equilibrium. Second, threadings form through a slow process whose characteristic time-scale is of the same order of magnitude as that of the diffusion of the rings in solution. These findings are robust, being based on universal (model-independent) observables as the average fraction of threaded length or the total penetrations between close-by rings and the corresponding distribution functions.

Equilibrated conformations of branched polymers in melt were prepared according to the algorithm by Rosa and Everaers [1]. The algorithm consists of four steps: Step 1: Melt of lattice trees -We construct randomlybranching polymers on the 3d-cubic lattice by resorting to a modified version of the Monte Carlo (MC) "amoeba" algorithm by Seitz and Klein [2] with periodic boundary conditions. The length of the unit cell of the lattice is equal to the Kuhn length, K /σ = 10, of the final ring polymer model. Here, σ is the unit of length (=monomer diameter) in the final ring polymer model. The average Kuhn segments density per unit cell is ρ K 3 K = 5. Initially, polymers are randomly placed on the lattice as standard random-walks. At each MC step, one of the monomers with functionality = 1 is randomly selected and attached to any randomly chosen monomer of the same chain with functionality < 3. The move is accepted with probability acc(i → f ) = min 1, where: (1) n 1 (i) and n 3 (i) (respectively, n 1 (f ) and n 3 (f )) are the total numbers of 1-and 3-functional monomers in the initial (respectively, final) state. * jan.smrek@univie.ac.at † kremer@mpip-mainz.mpg.de ‡ anrosa@sissa.it (2) n K (i, site) (respectively, n K (i, site)) is the total number of Kuhn segments inside the unit cell centered at the corresponding lattice site in the initial (resp., final) state.
(3) µ br = −2.0 and v K = 4k B T are two phenomenological parameters. The first is chosen by imposing that the fraction of branching points in each tree is ≈ 40%. The second fixes the scale for the free energy penalty accounting for excluded volume between overlapping pairs of Kuhn segments within the same unit cell.
Step 2: Conversion of lattice trees into off-lattice beadspring branched polymers -To resolve the spatial overlap between polymer conformations, we replace each Kuhn segment in the original lattice model with 5 linearly connected beads of diameter = 2σ. This is to assure that in the step 3 below we can build non-overlapping ribbons inside the corresponding occupied volume. Overlaps are then removed by the gentle "push-off" Molecular Dynamics (MD) procedure described in [3,4] and implemented in LAMMPS [5]. The procedure ensures that chain monomers are displaced by spatial distances of the order of their own linear size.
The Kremer-Grest-like [3] model adopted for the offlattice branched conformations accounts for the connectivity, bending rigidity, excluded volume and local topology conservation of polymers: 1. Bonded interactions between pairs of monomers at spatial separation r are modeled by the finitely extensible nonlinear elastic (FENE) potential. It is made of two terms. The first, attractive, is given by: where k = 30 /σ 2 is the spring constant, = 1k B T fixes the energy scale and the maximum extension of the bond ∆ + R 0 is determined by choosing ∆/σ = 1 and R 0 /σ = 1.5. The second, repulsive, is described by the Lennard-Jones (LJ) expression: where r c /σ = 2 1/6 . 2. Bending rigidity between consecutive triplets of neighbouring beads along the ring with spatial coordinates ( r i−1 , r i , r i+1 ) are modeled according to: where θ i = cos −1 ( ri+1− ri)·( ri− ri−1) | ri+1− ri| | ri− ri−1| is the angle formed by the adjacent bonds along the chain and k θ = 30 k B T is the bending prefactor. The large value adopted for the stiffness is based on the need for preventing excessive polymer crumpling during simulations [1].
3. The push-off procedure consists of two stages. During the first stage, non-bonded excluded volume interactions are described by the non-diverging, softcore potential: with r soft c /σ = 2.4225 and A ramped from 0 to 100k B T in an MD run of a few τ MD ≡ σ m/ where m is the "conventional" monomer mass. During this run, overlapping monomers are pushed away from each other. This run is followed by a second one of same duration and with non-bonded interactions changed to the full LJ potential, Eq. (3).
Step 3: Generation of double-folded ring conformations -Non-overlapping branched polymers can now be transformed into bead-spring, crumpled ring conformations. To do so, we start from a (randomly-chosen) 1-functional monomer of each branched polymer and we keep placing monomers by moving along the branched backbone while remaining at a distance of σ/2 from it. During the process the distance between nearest neighbors monomers is bound to the range [0.8σ − 1.2σ] so to avoid unnaturally stretched chain bonds. At branching points (3-functional sites), we choose randomly amongst the two possible remaining directions. This construction ends when the proximities of the initial monomer are finally reached. We end the protocol by checking once again if nearest neighbor distances along the chain stay in the interval [0.8σ − 1.2σ]. If not, we correct for this wherever needed. At the end of this step, the total number of monomers of each system is four times those of corresponding systems of branched backbones. Finally, Kuhn segments density is ρ K 3 K = 10 [1,6].
Step 4: Density relaxation on the entanglement scale -Bonded and non-bonded interactions between monomers in ring polymers are described by the same Eqs. (2) and (3) with ∆ = 0. Bending stiffness is instead replaced by the expression: corresponding to the fiber Kuhn length, K , given by [4]: For the present model k θ = 5k B T and Eq. (7) implies that the fiber Kuhn length K /σ = 10. To achieve complete system relaxation, we perform a short (= 10τ MD ) MD run under the condition that monomers cannot move more than 0.05σ at each integration time step. Finally, density fluctuations up to the entanglement scale are levelled-off by performing standard MD runs up to 0.75 of the (estimated) entanglement time τ e ≈ 1600τ MD [6].
In this work, we have considered the same ring polymer solutions studied in Ref. [1]. A summary of the rings and systems sizes considered as well as further details concerning their numerical analysis are summarized in the IBP-model column of Table S1.

Molecular Dynamics simulations details
The static and kinetic properties of the polymers are studied using fixed-volume and constant-temperature MD simulations (NVT ensemble) with implicit solvent and periodic boundary conditions. MD simulations are performed using the LAMMPS engine [5]. The equations of motion are integrated using a velocity Verlet algorithm, in which all beads are weakly coupled to a Langevin heat bath with a local damping constant Γ = 0.5/τ MD . The integration time step is set to ∆τ = 0.006 τ MD for push-off procedures and ∆τ = 0.012 τ MD for relaxation runs.

B. Concentrated solutions and melts of ring polymers equilibrated by Molecular Dynamics computer simulations
In this work, we have considered two distinct setsnamed hereafter "EQ MD 1" and "EQ MD 2", respectively -of solutions of ring polymers equilibrated by large-scale Molecular Dynamics (MD) computer simulations. The two sets are described by two different microscopic polymer models, each model characterized by a specific local chain stiffness (or Kuhn length of the chain, K ) and overall monomer density (ρ).
To ensure a proper comparison between the two data sets, the measured observables have been discussed in terms of the polymer mass (L c ) measured in entanglement lengths, Z ≡ Lc Le , where L e is the entanglement length [8]. For polymer contour lengths larger than L e , the microscopic details of the employed polymer model become irrelevant and the predicted properties of the solution universal.
In general, the entanglement length, L e , is a complicated function of the parameters of the solution like chain stiffness and monomer density. For loosely entangled solutions as the ones considered in this work, L e can be calculated from a "primitive path" analysis [9] according to the interpolation formula: where ρ K = ρ/( K /σ) is the density of Kuhn segments. In the following, we provide more details on the two polymer models considered in this work. A concise summary of the systems analyzed and the statistics of the different polymer configurations is given in Table S1.
Data for Z = 14, 29, 57 were produced by long, standard MD computer simulations as described in [10]. Instead, since relaxation times for the largest rings with Z = 114 (L c = 3200σ) would be computationally prohibitive, we have adopted the alternative strategy of doubling the contour length of each ring conformation with Z = 57 (L c = 1600σ). However, bare isotropic inflation of the system by a factor of 2 followed by insertion of a bead in between each two consecutive beads would just create threading statistics biased by the smaller system. To avoid this, we adopted the following procedure: On each ring a segment of length ds = 25σ was selected at random, a bead was inserted in between each two beads of the segment and the ring re-bonded to incorporate the new beads. The insertion was performed by inflating the bead diameter of the newly introduced beads from 0.5σ to σ in 50 steps of duration 5τ M D each. After the insertion, the system was relaxed for (20ds 2 )τ M D which is long enough (i) to completely relax the structures on scales of ds, because this is below the entanglement length and (ii) to reproduce the conformational statistics as found empirically when doubling shorter systems. The insertion and relaxation was repeated 64 times to reach the new length of N = 3200. Let us stress that in each round the doubling segment was chosen uniformly randomly on the ring. This ensures that we are not biased by the threading present in the original system of rings. The doubling process was ran at constant pressure maintained by Berendsen barostat with τ P equal to the square root of twice the number of particles. Here the parameter τ P governs the relaxation of the box size by rescaling all coordinates by a factor of µ = 1 − βdt 3τ P (P 0 − P ) as r new = µ r old , where dt = 0.005τ M D is the time step used, β is the isothermal compressibility (set to unity) and P 0 − P is the difference between the target and the actual pressure (see [11] for details). After the ring mass of N = 3200 monomers was reached, the system was switched from N P T to N V T simulation by re-scaling all coordinates by a common factor in order to match target monomer density = 0.85σ −3 . The factor was typically less than 1.001 as the pressure of the system was almost perfectly equilibrated. The system was then run for another 2 × 10 6 τ M D . The total run time was about 2.8 × 10 6 τ M D , which is significantly shorter than the predicted diffusive relaxation time but it is about the same as the predicted conformational relaxation time for these polymer rings [10].
EQ MD 2 -The second data set is for MD simulations of rings whose initial states correspond to doublefolded polymers on branched primitive trees (IBP-model, Sec. SI A). In this case, K = 10σ and ρ = 0.1σ 3 imply (see Eq. (8)) L e = 4 K = 40σ. The total run time for each of these simulations was fixed being equal to 1.2 × 10 7 τ MD or ≈ 7500τ e .
As anticipated in the main text, ring dynamics is discussed in terms of two observables [10]: (1) Rings thermal diffusion can be used to monitor systems equilibration. Accordingly, we define the time mean-square displacement [3]:    tory and r m cm (t) ≡ 1 N N i=1 r m i (t) is the spatial position of the center of mass of the m-th ring of the system whose ith monomer has spatial coordinates r m i (t). g 3 (t) is then compared to the time evolution of the ring mean-square gyration radius: The ring relaxation time τ diff rel is then defined by the relationship: As shown in Fig. S1, our trajectories are long enough to achieve good equilibration for rings with Z = 5, 15, 38. Conversely, rings with Z = 115 are approximately equilibrated and corresponding data need to be taken with care. Specific data for τ diff rel = τ diff rel (Z) are summarized in Table S2.
(2) The second observable was introduced in Ref. [10] as a mean to characterize rings internal motion. For any ring in solution at time t we consider an arbitrarily chosen pair of vectors d 1 = d 1 (t) and d 2 = d 2 (t), each vector connecting two monomers of the same ring separated by a contour length Z/2. The two vectors are chosen so that the tails are separated by a contour length = Z/4. By taking the cross product c(t) = d 1 (t) × d 2 (t), the internal relaxation time τ int rmrel is defined by the integral over time of the corresponding correlation function: Numerical evaluations of Eq. (12) at each Z have been performed over the corresponding full MD trajectory, by stopping the integration whenever the relative error on the correlation function becomes > 20%. This choice appears reasonable as c(t)· c(0) c(0) 2 = O(10 −2 ) at the cut-off for all Z's. Our final results are summarized in Table S2: in agreement with Ref. [10], τ int rel (Z) < τ diff rel (Z).

C. Construction of minimal surfaces for ring polymers
To obtain the minimal surfaces spanned on MDequilibrated and branched ring polymers in nonprohibitive computing time, we follow a slightly modified version of the numerical procedure described in [12]. Comparison of the minimal surface from the original minimization algorithm presented in Ref. [12] (left) and the one obtained from the coarser procedure adopted in this work (right) for a single ring with Lc = 800σ (Z = 29) extracted from data set EQ MD 1 (Sec. SI B).
Firstly, we initialize the surface as a union of L c triangles, where each triangle is spanned on two successive monomer positions and the center of the mass of the ring. Next, we refine the surface once, by subdividing each triangle into four smaller ones. The surfaces, represented by 4L c triangles, are subsequently minimized by mean curvature flow with restructuring of the mesh towards Delaunay triangulation (equiangulation), weeding out of small triangles and vertex averaging, as discussed in [12]. After the initial minimization (that takes about 20000 steps of the mean curvature flow) the relative surface area change is being measured during the run. Surface evolution is stopped when the relative surface area does not change by more than 0.1% over 240 steps with equiangulations and, in case of decreasing time step, together with vertex averaging. We have found empirically that the above procedure produces surfaces close enough to a (local) minimum that is characterized by a vast majority of vertices (> 99%) with the absolute value of the mean curvature smaller than 0.01. As shown in [12] the use of several different protocols for simulated annealing ensures that the procedure leads sufficiently close to the minimal surface. The minimal surfaces presented here are slightly coarser than the ones presented in [12], mainly because of using only one mesh refinement instead of two. Yet, the past and the present methodologies give remarkably close results, as illustrated in Fig. S2. The distribution of the minimal areas is similar in both cases and the mean is higher by only about 4% for coarser surfaces. As a final check, we have also verified that the threading statistics is not sensitive to these small changes (Fig. S3). The following script (in Surface Evolver language) has been used to minimize the surfaces. For few surfaces, this procedure was unable to converge due to internal surface evolver error (when trying to free a loop edge). In those few cases the w commands were set to 0.1, 0.001 or 0.0001 which solved the issue.

D. Threading detection
We define the penetration of the minimal surface of ring A by ring B iff a line segment between any two consecutive monomers of ring B intersect the internal area of any triangle from the triangulation of the minimal surface of ring A. As we label the line segments of B, it is straightforward to compute the (threading) length between two successive penetrations. We observe that a ring can also penetrate its own surface, however, in this work, we focus only on inter-ring threadings and we do not take self-threadings into account.
Determining the relative position of a line segment and a triangle is a standard problem in computer graphics that can be addressed in various ways (see e.g. one pedagogical approach in the appendix of work [13]).
In this work, we follow the algorithm presented in [14] that uses Plücker coordinates and the side operator: • Plücker coordinates are a six-dimensional representation of a line segment a defined by two boundary points in 3d space p = (p x , p y , p z ) and q = (q x , q y , q z ). Each Plücker coordinate is one determinant of the 2 × 2 minor of the matrix Explicitly the coordinates of the line segment a are a = (p x q y − q x p y , p x q z − q x p z , p x − q x , p y q z − q y p z , p z − q z , q y − p y ). Sometimes another representation with clearer geometrical meaning is used a = ( p − q, p × ( p − q)), which differs only by permutation and signs.
• The side operator between two lines a and b is a permuted inner product in the six-dimensional space defined as s(a, b) = a 0 b 1 +a 1 b 2 +a 2 b 3 +a 3 b 4 +a 4 b 5 + a 5 b 0 .
It can be shown that if all side operators of the line segment with the three triangle edges have the same sign, the line intersects the inner area of the triangle. From the values of the side operators one can extract more information of the relative positions of the line and the triangle such as the intersection at a vertex of the triangle (see [14] for details).

E. More details on threading statistics
From the presented statistics, it is not clear what fraction of its neighbors each ring threads. To quantify it, we have constructed histograms of the number of threaded neighbors p(n t ) of a ring (Fig. S4) for the two sets of MD-equilibrated systems. As for the average number of times a ring penetrates the minimal surface of any other single ring (n p , Fig. 2 in the main text), we considered a neighbor to be threaded, only if at least one threading length L t that contributes to L sep (Eq. (1) in the main text) is longer than the entanglement length L e . Notice that if threadings of any length were considered, as it was done in [12], the mean values of n t agree with the data reported there (see Fig. 2a of that paper). In this case though, p(n t ) for the EQ MD 1 data set would not agree with the EQ MD 2 data set because of natural differences in the microscopic details of the two models. With the new definition adopted here, we see that the two data sets agree very well (Fig. S4, see for instance the overlap between the two results for Z = 114 (EQ MD 1) and Z = 115 (EQ MD 2)).
In general, the n t distributions appear unimodal for all Z. Only the very short rings (Z ≤ 15) have nonvanishing value at n t = 0. This means that every ring threads at least one neighbor, but often many neighbors. Surprisingly, the longest rings (Z = 114) have still higher averagen t than for Z = 57, although in the asymptotic limit Z → ∞, one would expect this number to saturate due to the compact regime of the ring size (R g ∼ N ν with ν = 1/3). Interestingly, the variance of the n t distribution grows with Z, which could have an important effect on the dynamics. We do not evaluate the fraction of neighbors that are threaded as this should depend on how we define a neighbor. In [7] the neighbor number K(Z) is defined as the rings with center of mass within distance R g . There the K(57) 16. Here, we see that when neglecting the short threadings, every ring threads fewer rings than K. On the other hand when threadings of any length are considered, the n t is larger than K (not shown), because also short threadings for distances larger than R g contribute.

II. ROLE OF NUMBER OF CHAINS INnp STATISTICS
We investigated briefly what is the effect of small number of chains on the threading statistics. We used the polymer model of EQ MD 2 and run MD to equilibrate two more sets of systems: (i) EQ MD 3: the starting configurations are (Hilbert or Moore) space-filling curves (SFC) as detailed in [1] with M = 8 chains each and for Z = 14, 38 and 115 and (ii) EQ MD 4: the starting configurations are IBP as in EQ MD 2, but the number of chains is M = 16 for Z = 38 and M = 8 for Z = 115. Fig. S5 summarizes our findings: 1. The equilibrium values of EQ MD 1 is below the equilibrium value of EQ MD 3 for both, the Z = 15, and for Z = 115. We speculate that this is because in this case the low chain number M = 8 causes that each chain interacts with itself through periodic boundary conditions and this could cause a "more open" conformation which is thus more easily threaded by other chains.
2. Then p of late (t > 10 4 τ e ) EQ MD 4 simulation for Z = 38 with M = 16 is below the EQ MD 3 with M = 8 and above EQ MD 2 with M = 256, which supports the hypothesis above.