What Happens When Threading is Suppressed in Blends of Ring and Linear Polymers?

Self-diffusivity of a large tracer ring polymer, Dr, immersed in a matrix of linear polymers with Nl monomers each shows unusual length dependence. Dr initially increases, and then decreases with increasing Nl. To understand the relationship between the nonmonotonic variation in Dr and threading by matrix chains, we perform equilibrium Monte Carlo simulations of ring-linear blends in which the uncrossability of ring and linear polymer contours is switched on (non-crossing), or artificially turned off (crossing). The Dr≈6.2×10−7Nl2/3 obtained from the crossing simulations, provides an upper bound for the Dr obtained for the regular, non-crossing simulations. The center-of-mass mean-squared displacement (g3(t)) curves for the crossing simulations are consistent with the Rouse model; we find g3(t)=6Drt. Analysis of the polymer structure indicates that the smaller matrix chains are able to infiltrate the space occupied by the ring probe more effectively, which is dynamically manifested as a larger frictional drag per ring monomer.

Uncontaminated and unknotted RPs are difficult to produce and isolate, but crucial; even a modest amount of contamination by linear polymers (LPs), often by-products of prior reactions, can drastically slow down the overall dynamics. For example, Kapnistos et al. found that intentional contamination of cyclic polystyrene rings (160 kDa) with less than 1% linear chains was sufficient to dramatically retard the linear viscoelastic response [15]. Robertson-Anderson and coworkers used fluorescence microscopy to study the dynamics of probe circular DNA, in different matrices, and found significant topological effects [16][17][18]. Recently, the use of RPs as probes of LP melts via neutron spin echo spectroscopy has been been pioneered, in which perturbation of the internal dynamics of the RP is used to glean insights about the matrix it is embedded in [19,20].
As a result of these findings, ring-linear blends (RLBs) have emerged as an important subject of scientific investigation in their own right. In this paper, we focus primarily on RLBs in which the concentration of the RP is small enough to regard these systems as ring probes diffusing in a LP matrix. In such probe or tracer RP systems, for sufficiently large molecular weight constituents, threading of RPs by LPs is implicated in arrested dynamics. This was first recognized in polystyrene tracer diffusion studies in the 1980s [21][22][23], and reestablished more recently with DNA tracer diffusion studies [16]. Computer simulations of flexible and semiflexible RLBs have yielded valuable insights D N m Figure 1. Diffusion of linear (blue) and ring (red) probes in linear matrices (gray). The diffusivity of the linear probe decreases monotonically as the molecular weight of the matrix polymers increases, before reaching a plateau. On the other hand, the diffusivity of ring probes varies non-monotonically. At large N m it decreases without reaching a plateau.
In this paper, we re-examine this speculation more carefully, by conducting fresh simulations in which threading between RPs and LPs is artificially suppressed, by letting ring and linear contours cross each other. An advantage of simulations is that we can carry out such "unnatural experiments", which would be hard or impossible to perform in experiments. They are performed with the intention of isolating the effect of particular interactions, while leaving other interactions untouched.

Materials and Methods
We use the bond-fluctuation model (BFM) [50], which is a lattice Monte Carlo method, to simulate the RLBs. In the past, we have successfully used the BFM to study RLBs [51][52][53][54], due to its ability to efficiently explore long-time scales, and relatively large system sizes. In the BFM, the C ∞ ≈ 1.2 at the typical melt density, compared to C ∞ = 1.74 in the Kremer-Grest off-lattice model [50]. Although the polymer represented in the BFM is flexible, it has been successfully mapped to solutions of dsDNA [18], by matching the number of Kuhn segments.

Bond-Fluctuation Model
In the BFM, 3D space is resolved into simple cubic grid of size L × L × L with periodic boundary conditions. We place n r non-concatenated RPs, each comprising N r beads or monomers, into a matrix containing n l LPs, each comprising N l monomers. To reproduce melt-like behavior, the total density of occupied lattice sites is maintained at ρ = (n r N r + n l N l )/L 3 = 0.5.
Trial moves, in which a randomly selected monomer is displaced to one of its 26 neighboring sites, are attempted. A trial move is accepted if three constraints are satisfied [50]: 1. Excluded Volume-the target lattice site is empty and available, 2. Finite Extensibility-none of the bond-lengths stretch beyond √ 3, and 3. Chain Uncrossability-mid-points of bonds do not intersect.
One Monte Carlo Step (MCS) corresponds to n r N r + n l N l trial moves. The units of distance and time are lattice spacing, and MCS, respectively. The protocol for equilibrating a RLB has been described previously [51][52][53]55]. We monitor the decorrelation of the end-to-end vector and the vector connecting beads 1 and N/2, for LPs and RPs, respectively. Equilibration is terminated once the correlation falls below the threshold value of 0.05. This brute-force protocol yields the correct internal bead distance distribution for situations where it is analytically known.
If the uncrossability constraint is relaxed, chains are allowed to cross and pass through each other. They still have to obey the excluded volume and finite extensibility conditions. Such simulations have been previously peformed on pure LP [50] and RP melts [56]. In this paper, we describe two types of simulations; in the regular or "non-crossing" (NC) simulations, the chain uncrossability constraint is strictly enforced. In the "crossing" (CX) simulations, the chain uncrossability constraint is selectively relaxed only for ring-linear interactions. Therefore, contours of RPs and LPs are allowed to pass through each other. However, uncrossability is strictly enforced for ring-ring, and linear-linear interactions. Thus, these CX simulations, allow us to isolate and explore the dynamical consequences of suppressing threading between RPs and LPs, while leaving other interactions untouched. The acceptance ratio of the Monte Carlo moves for large polymers is found to be independent of polymer architecture; for NC simulations it is about 0.18, while it is about 10% higher for CX simulations.

Self-Diffusion Coefficient
After equilibration, we perform production runs for simulation time τ sim . We monitor the mean-squared displacement (MSD) of the center of mass, where r C denotes the center-of-mass, and · denotes an average over all polymers and time-intervals t.
The self-diffusivity of the polymers can be obtained from g 3 (t), using the Einstein formula: The simulation time was chosen to be long enough to ensure that the polymers had diffused, on average, at least five times their radius of gyration. We used statistical bootstrap [57,58] to infer confidence intervals for the estimated self-diffusivities.

Systems Studied
We studied two sets of RLBs, (i) symmetric; and (ii) asymmetric. The details of these systems are presented in Table 1. In the bond-fluctuation model, at ρ = 0.5, the average number of monomers per entanglement segment is N e ≈ 30 [59][60][61].
• In symmetric blends, N r = N l = 300 was held fixed, while the linear fraction φ l = n l N l /(n l N l + n r N r ) was varied between 0 and 1. Note that for φ l = 0 (pure rings) and φ l =1 (pure linears), the crossing simulations and the non-crossing simulations were identical, since there are no ring-linear interactions to suppress in these pure systems. • In asymmetric blends, the concentration of the matrix LPs φ l = 0.9, and the number of RP monomers N r = 300, were held fixed, while the number of LP monomers was varied between N l = 10 − 300. To avoid ring-ring interactions in these probe systems, the concentration of the ring polymers was kept about 10 times lower than the overlap concentration [34]. Table 1. Description of the symmetric (N r = N l ) and asymmetric blends (N r = N l ) blends simulated in a periodic cubic box with L = 60, at a total density of ρ = 0.5. Simulation times for the non-crossing (NC) and crossing (CX) simulations are in units of 10 6 Monte Carlo Steps.

Results
In the following, we discuss the static and dynamic properties of the symmetric and asymmetric blends. The results of all the NC simulations have been previously reported, including the size [51], entanglement structure [53], free energy [54], and self-diffusion [55] of the symmetric N r = N l = 300 blends, and the size and diffusivity of the asymmetric probe ring blends [33,34].

Symmetric Blends
When the non-crossing constraint is relaxed in the symmetric RLBs, the change in static properties is barely perceptible, while the self-diffusivity profiles change considerably. Symmetric blends provide a baseline from which the results of the more interesting asymmetric blends can be analyzed. Figure 2 plots the squared radius of gyration, R 2 , of the RPs and LPs in the blend. In the NC simulations, the mean radius of gyration of the LPs R NC l ≈ 11.6. Superscripts "NC" and "CX" are used to distinguish properties extracted from non-crossing and crossing simulations, and the subscripts "r" and "l" are used to represent RPs and LPs, respectively. In the crossing simulations, the mean R CX l increases from 11.6 ± 0.2 to 12.0 ± 0.4, as φ l decreases from 1.0 to 0.2. In previous comparisons of NC and CX simulations of pure LPs [50,62], no significant differences in polymer size were observed. However, in those simulations, the LPs were allowed to pass through each other, unlike our simulations in which LP-LP crossings are prohibited, and only LP-RP crossings are allowed. RPs in pure melts adopt highly compact conformations due to the non-catenation constraint between neighboring rings. Brown et al. [56] performed CX simulations of pure ring melts, in which the non-catenation constraint was relaxed by allowing RP-RP crossings. This caused the N r = 300 pure ring melts to swell from R NC r = 7.0 to R CX r = 8.3. In both our CX and NC calculations, RP-RP crossings are not permitted. However, a similar effect is observed by contaminating a pure ring melt with increasing levels of LPs. As the linear fraction φ l increases (Figure 2), LPs weaken the strength of the non-catenation constraint due to dilution of ring-ring interactions, causing the RP to expand. In RLBs that are LP-rich (large φ l ), the RPs adopt conformations consistent with Gaussian rings [51]. In the NC simulations, the size of the RP increases from 7.0 ± 0.1 to 8.1 ± 0.1, as φ l increases from 0 to 0.9. In the CX simulations, the RP size increases from 7.0 ± 0.1 to 8.3 ± 0.1.

Statics
The small enhancement in the LP size in RP-rich blends, and in RP size in LP-rich blends for the CX simulations relative to the corresponding NC simulations can be undestood in terms of a local "solvation effect". For example, a LP surrounded predominantly by "crossable" RPs in RP-rich environments feels that it is in a (partially) good solvent, prompting it to swell. The degree of relative swelling depends on the concentration of the opposite species in the blend. Figure 3 depicts the diffusivity of the LPs and RPs in the blend for NC and CX simulations. In the NC simulations, diffusivity of the LP, D l ≈ constant, although there is a shallow minima at large φ l , which is also observed experimentally [18]. On the other hand, the diffusivity of the RP decreases precipitously from the pure melt (φ l = 0) as the linear fraction increases. This is due to threading of the RPs by the LPs; the RPs are pinned down by LPs, and are effectively immobilized on the diffusion timescale (τ l ∼ R 2 l /D l ) of the threading LPs. When ring-linear threading is artifically switched off in the CX simulations, the RP is no longer constrained, and its diffusivity actually increases from D CX r = 1.4 × 10 −5 to 2.3 × 10 −5 as φ l increases from 0 to 0.8. The increase in linear in φ l and is well-described by the relation, D CX r (φ l ) = 1.4 × 10 −5 (1 + 0.8φ l ). The increase in D CX r with increasing φ l is due to the replacement of uncrossable ring neighbors with LPs, which a RP can cut through. This can again be thought of as a manifestation of the solvation effect, which was responsible for a small increase in R r ( Figure 2) with φ l . For the LPs, a corresponding solvation effect is manifested by the increase in D l from 2.2 × 10 −6 to 1.9 × 10 −5 as φ l decreases from 1.0 to 0.2. The data is well-described by D CX l ≈ 2.2 × 10 −6 exp(2.6(1 − φ l )). As we move away from a LP-rich to a RP-rich environment, the fraction of the medium offering topological resistance to a LP goes down. In the CX simulations, the RP or the LP sees the opposite species as a spatial correlated high-density solvent, offering some resistance. In other words, the presence of LPs as φ l → 1.0 slows RPs below their mobility in the absence of LPs. This is evident from the CX computations of Brown et al. [56] referenced earlier. In their work, the found that the diffusivity of N r = 300 RPs in a pure ring melt increased from 1.3 × 10 −5 to 4.1 × 10 −5 , when the RP-RP crossing was allowed. However, this diffusivity is still smaller than the diffusivity of single isolated non-crossing (D r = 1.4 × 10 −4 ) or crossing (D r = 1.6 × 10 −4 ) RPs.

Ring Probes in Linear Matrix
The results of the static and dynamic properties of symmetric RLBs seem intuitive. We now turn our attention to asymmetric blends, with N r = 300, in LP-rich matrices comprised of varying molecular weights N l .

Statics
The results of the NC simulations have been reported in a recent publication [34]. To summarize, the radius of gyration of the LPs R NC l ∼ N 1/2 l , while that of the probe RP is essentially a constant over a large range of matrix molecular weights. As N l falls below N * l = √ N r , the matrix LPs act as a solvent causing the ring probe to swell by about 10% (Figure 4).
R CX l of the LPs in the CX simulations essentially overlaps with the R NC l from the NC simulations. This is anticipated from Figure 2. At φ l = 0.9, a LP is essentially surrounded by other non-crossable LPs. The presence of a few crossable probe RPs in the melt induces the LPs to expand; however, their concentration is too low to produce significant observable differences.
On the other hand, the size of the ring probes is noticeably different in the CX and NC simulations. At φ l = 0.9, the neighborhood of a RP consists predominantly of crossable LPs, which allow the RP to swell. Indeed the plateau value of R 2 r for N l ≥ 75 increases from 65.0 ± 3.4 to 75.4 ± 3.8 as we move from NC to CX calculations (Figure 4). The upturn observed in R r for N l < N * l , is muted in the CX simulations due to the weakened impact of the solvation effect.  Figure 5 depicts the self-diffusion constant of the CX and NC probe ring systems. The NC simulations were previously reported [34]. The diffusivity of the LPs in the blend D l is a monotonically decreasing function of N l . At φ l = 0.9, D l tracks the diffusivity of pure LPs, which for N l > 75, varies as D l ∼ N −2.4 l [34]. The diffusivity of the RP D NC r , on the other hand, exhibits a surprising non-monotonic behavior, which was alluded to previously. Primitive path analysis suggested that for RP probes, there were two regimes: when N l was small, the number of ring-linear entanglements decreased with N l , eventually crossing over to a plateau as the size of the LP matrix chains became comparable with that of the RP. It was speculated that the increased degree of entanglement more than offset the effect of the increasing mobility of the matrix chains, in the small N l regime, leading to a non-monotonic variation in D r .

Dynamics
The dashed blue lines and symbols in Figure 5 depict D l and D r in the CX simulations. At φ l = 0.9, a LP sees a LP-rich environment. Since LP-LP contour crossing events are prohibited, D l in the CX simulations closely tracks the D l observed in the NC simulations. The presence of the crossable RPs causes the D CX l > D NC l , as the LPs are marginally more mobile. This effect is stronger at N l = 300, than at N l = 10, and as N l decreases, the enhancement in mobility weakens systematically. When N l falls below the entanglement threshold N e ≈ 30, the strength of the non-crossing topological constraint is reduced even in the NC simulations, and the difference between the CX and NC simulations becomes less important.
This reduction in the significance of ring-linear threading as N l decreases is also evident in the diffusivity of the RP, D r . In Figure 5, the values of D r in both the CX and NC simulations merge at small N l , as one would expect. Indeed, the D CX r in the CX simulations provides an upper-bound for D NC r . For N l < N e , the LPs do not effectively constrain the RPs. As N l increases to 75 (about 2N e ) and above, threading of the RPs by the matrix chains begins controlling the long-time dynamics of the RP, and the D r in the CX and NC simulations start to diverge. In the NC simulations, the D r decreases as N l increases. In the CX simulations, for N l ≥ 30 the power-law dependence D CX r ≈ 6.2 × 10 −7 N 2/3 l describes the data quite well.

Discussion
When Shaffer [50] performed BFM simulations with pure LP melts by switching the uncrossability criterion on and off, he found Rouse-like scaling of the LP diffusivity in the crossing simulations The key question to address, to completely understand the non-monotonic varition of D NC r (N l ), is "why does D CX r increase monotonically with N l in Figure 5?" We saw from Figure 4, that the size of the LP in the CX simulations was relatively unchanged from the NC simulations, while the RP was somewhat more expanded (except at the smallest N l explored). We wanted to examine how this affects the local microenvironment of a RP; in particular, to quantitatively address the question, "how many neighboring polymers infiltrate the space occupied by a ring probe?" Therefore, we considered a probe RP and computed the number of polymers of either species (RP or LP) whose centers-of-mass were contained within the radius of gyration R r of the probe. Figure 6 shows that the number of RPs within this region is nearly zero for both the CX and NC simulations. This is expected, since at φ l = 0.9, the concentration of the RPs is significantly below its overlap concentration. The number of LPs within this region increases as N l decreases. As N l decreases, the LPs become smaller in size, and are better able to sneak into small voids available within a RP. The inset to the figure shows the number of LPs on a double logarithmic scale. The difference between the CX and NC simulations becomes smaller as N l decreases. For larger N l , the number of LPs in the CX simulations lies above that in the NC simulations, presumably due to the larger size of the probe rings in these simulations (see Figure 4). Since the number of LPs within R r (let us label it n R l ) varies with N l as n R l ∼ N −1 l from the inset to Figure 6, one can conclude that n R l N l ≈ constant. This suggests that the number of LP monomers contained within R r of a ring probe, n R l N l , is independent of the molecular weight of the LP. In other words, fraction of sites occupied by the LPs in the local microenvironment of a RP is roughly unchanged as a function of N l . For the NC simulations, there appears to be a stronger dependence of n R l on N l at smaller N l . This observation is consistent with previous primitive path simulations [34], which showed enhanced ring-linear entanglement in this regime. The shorter LPs are more mobile, and better dispersed within the pervaded volume (see Figure 7). They can thus offer better frictional resistance to the internal modes of motion of the RP. Figures 8 depicts the g 3 (t) for the crossing and non-crossing ring probes as N l is varied, over small to intermediate timescales. The g 3 (t) for the NC RPs are similar to the g 3 (t) curves obtained in previous tracer diffusion studies [34]. In tracer diffusion studies of a large LP immersed in a matrix of shorter LPs, the g 3 (t) curves of the tracer LP become independent of the matrix molecular weight at short timescales. In constrast, the g 3 (t) curves of a tracer or probe RP at short timescales increases with increasing N l -a characteristic which is also observed in simulations of pure ring melts of varying molar mass [34]. As N l decreases, the matrix chains are able to infiltrate the space occluded by a ring probe more effectively. Therefore the g 3 (t) of the RP feels the effect of the matrix chains at shorter times. At intermediate timescales, the effect of threading becomes visible, especially for the longest N l . Unlike the complicated structure of the g 3 (t) curves in the NC simulations, the g 3 (t) curves in the CX simulations are simple. The terminal diffusive regime, as attested by the transition to the t 1 scaling, appears to be attained at smaller timescales. It is clear that the ring probes are faster at all timescales in the N l = 300 matrix compared to the the N l = 10 matrix. The g 3 (t) curves are effectively parallel to each other, with the diffusion constant reported in Figure 5 serving as the scaling factor. Indeed, this is strongly reminiscent of the g 3 (t) curves for RPs expected from Rouse theory [63]. In the Rouse model, g 3 (t) = 6D r t, with D r = k B T/N r ζ, where k B is the Boltzmann constant, T is the absolute temperature, N r the number of RP monomers, and ζ is the drag per bead.
From Figure 8, the g 3 (t) curves for the probe RP depend on N l . Since D CX r ∼ N 2/3 l is an increasing function of N l , it implies that ζ ∼ N −2/3 l is a decreasing function of N l . In other words, the effective drag on the beads of the Rouse RP is large when the matrix LPs are small, and are able to enter and disperse into the volume pervaded by the RP more effectively. Dynamics of isolated RPs in a matrix of fixed obstacles at low obstacle density [64] also show qualitatively similar profiles for g 3 (t). As the obstacle density is increased from zero, the frictional drag per RP segment increases. The RP is more effectively slowed down, even as the shape of g 3 (t) remains essentially unchanged. Beyond a certain obstacle density, the shape of the curve changes.
Finally, Figure 9 compares the g 3 (t) curves between the NC and CX simulations from Figure 8, directly for three different values of N l . When the LP is small, N l ≈ 10, the g 3 (t) of the RP probes is nearly the same in the CX and NC simulations. As N l increases to the entanglement threshold and beyond, the g 3 (t) of the probes in the CX simulations show greater mobility than in the NC simulations. For N l = 75 in Figure 9, the g 3 (t) of the crossing and non-crossing curves are roughly parallel over the timescales reported. As N l increases further, to say N l = 300, and the threading by LPs becomes more persistent, the two lines are no longer parallel to each other. The signature of the threading by the LPs is evident in the decreased slope in the t > 10 5 time range.

Proposed Experiments and Simulations
Although the BFM models a flexible polymer, the provocative results reported in this work have direct implications for semiflexible polymers. First, most practically realizable RP systems are composed of semiflexible polymers; therefore to experimentally validate the BFM observations, we need to figure out the appropriate experiments to perform. Fortunately, the BFM has previously been mapped to solutions of dsDNA polymers. In the BFM, at the densities studied, the relationship between the number of monomers, and the number of Kuhn steps N K is N K = 0.83 N. Using the blob theory, the number of Kuhn steps corresponding to a 45 kbp dsDNA in a 1 mg/mL solution is found to be N K = 263 [18], which corresponds to N = 263/0.83 ≈ 317 in the BFM. Thus, the N = 300 probe, studied in this and previous work [34], is quite similar in size to a 45 kbp dsDNA in a 1 mg/mL solution. Thus, the N = 10-150 matrix polymers correspond approximately to dsDNA of lengths between 1.5-22.5 kbp. Thus experiments in which the diffusivity of a 45 kbp RP is measured in LP matrices whose lengths are varied between 1.5-22.5 kbp at a total concentration of 1 mg/mL will be useful to validate the BFM results.
Furthermore, the effect of stiffness in both the probe and the matrix polymers is an interesting topic to explore, especially using an off-lattice model. The observed non-monotonic effect is due to non-intuitive interplay between the structure and the dynamics of the probe RP. While the semiflexibility of the polymers definitely affects the structure, whether it enhances or subdues the observed non-monotonicity is an open question.

Summary and Conclusions
Unlike LP probes, the diffusivity of RP probes D r in linear matrices, exhibits an unusual non-monotonic dependence on the matrix molecular weight. For N r = 300, the diffusivity initially rises for N l < N e , reaches maxima around N l ≈ 2 N e , and then decreases with N l due to threading by matrix LPs. Previous work on the entanglement structure suggested that the non-monotonicity in D r could result from a competition between the number of ring-linear threading events and their persistence.
In order to investigate this hypothesis directly, we performed equilibrium Monte Carlo simulations of ring-linear blends using the BFM. In the BFM, the uncrossability criterion, which ensures that contours of polymers do not cut through each other, can be selectively switched off. In our crossing simulations, we turned off the uncrossability condition between RPs and LPs, while keeping it on for ring-ring and linear-linear interactions. This surgically eliminates the role of threading in the CX simulations, thereby allowing us to appreciate its role in the regular NC simulations.
We performed simulations on two sets of systems: symmetric and asymmetric. For the symmetric simulations, the change in the static and dynamics properties in the CX simulations relative to the NC simulations was expected. The size and self-diffusivity both increased as the composition of the opposite species in the blend was increased. The increase in size was marginal, while the increase in self-diffusion coefficient was more dramatic. The changes could be interpreted by a "solvation effect": the replacement of neighbors by "crossable" polymers of the opposite species allows the polymers to expand, and move faster.
We then explored the size and diffusivity in LP-rich asymmetric blends. The linear fraction was maintained at φ l = 0.9, N r = 300 was held fixed, and N l was varied between 10 and 300. For LPs, there was no significant change in R 2 or D l between the NC and CX simulations. This was expected, since at φ l = 0.9, the environment of a typical LP consists of mostly other LPs. The RP probe expanded slightly due to the solvation effect mentioned above. In the CX simulations, D CX r increased monotonically with N l and provided and upper-bound for D NC r . For N l ≥ 30, it was found that D CX r ≈ 6.2 × 10 −7 N 2/3 l . Investigation of the center-of-mass mean-squared displacement showed that g CX 3 (t) = 6D CX r t), consistent with the Rouse model over all the timescales (Figure 8). Analysis of the polymer structure in the CX simulations indicated that the number of LPs that infiltrate the volume occupied by a RP scales as N −1 l . In the NC simulations, for large N l the number of infiltrating LPs also varied as N −1 l , while the slope was stronger at shorter N l . This is consistent with previously reported primitive path analysis on the NC systems [34], which showed an increase in the topological interactions in this regime. Together, they indicate that the space occupied by the RP is incrementally enriched with matrix polymers as N l decreases as shown in Figure 7. This infiltration is manifested as a larger frictional drag per RP bead in the g CX 3 (t) curves.