Ring polymers in two-dimensional melts double-fold around randomly branching “primitive shapes”

Abstract

Drawing inspiration from the concept of the “primitive path” of a linear chain in melt conditions, we introduce here a numerical protocol which allows us to detect, in an unambiguous manner, the “primitive shapes” of ring polymers in two-dimensional melts. Then, by analysing the conformational properties of these primitive shapes, we demonstrate that they conform to the statistics of two-dimensional branched polymers (or, trees) in the same melt conditions, in agreement with seminal theoretical work by Khokhlov, Nechaev and Rubinstein. Results for polymer dynamics in light of the branched nature of the rings are also presented and discussed.

---

1 Introduction

Edwards introduced the concept of a “primitive path” long ago¹ to explain the viscoelastic behavior of entangled polymer solutions and melts (see panel (a) of Fig. 1). This idea has since become a cornerstone of polymer physics. In a polymer melt the motion of each chain is constrained by its neighbors, which perform as uncrossable, topological obstacles caging the chain within a tube-like region.² The primitive path represents the central axis of this tube, and it is identified as the shortest path which the polymer chain can be contracted into, at fixed chain ends and without crossing the topological obstacles represented by the surrounding chains or each other.^3,4 This definition, when applied to computer-generated conformations of entangled linear polymer melts,⁵ allows researchers to map a melt of polymer chains into a melt of primitive paths. Quite remarkably, this mapping provides a universal framework⁶ for describing the macroscopic rheological properties of entangled polymer liquids, regardless of their specific chemical composition.

Fig. 1 (a) Classical illustration of the primitive path (dashed line) of a linear chain (solid line) in an entangled polymer melt. By representing the effects of the constraints exerted by neighboring chains as a network of uncrossable topological obstacles (×), the primitive path is identified as the shortest path which the polymer chain can be contracted into, at fixed chain ends and without crossing the obstacles. (b) In a melt of unknotted and non-concatenated rings, the equivalent topological barriers induce each ring to double-fold around a randomly branching path.

The conventional definition of a primitive path is intuitive and simple, but it relies on the presence of free chain ends. This raises a problem: can this concept be generalized to entangled polymers with no free ends?

In this work we explore this question by considering melts of unknotted and non-concatenated ring polymers, a particular class of systems that has challenged theorists,^7–17 experimentalists^18–23 and computational^24–31 physicists for decades. In particular, Khokhlov and Nechaev⁷ and Rubinstein and coworkers^8,10 suggested that, because of the constraint of non-concatenation, rings minimize their mutual overlap by double-folding around randomly branching, or tree-like, primitive paths (see panel (b) of Fig. 1).

The main goal of this work is to provide a direct numerical test of the Khokhlov–Nechaev–Rubinstein picture. Specifically, we introduce a numerical algorithm to straightforwardly and unambiguously extract the primitive path (to which, to distinguish it from the primitive path of linear chains, we prefer from now on the term “primitive shape”) of rings in two-dimensional computer-generated melts. Our results show that these primitive shapes exhibit complex, non-trivial behavior which is, indeed, strongly reminiscent of branched, tree-like conformations. To demonstrate that, we present a detailed characterization of their statistical properties using the set of observables and distribution functions previously introduced by one of us (A.R.) and R. Everaers^32–34 for studying two- and three-dimensional lattice trees in various dilution conditions.

The paper is organized as follows. In Section 2.1 we introduce the polymer model and the kinetic Monte Carlo algorithm used for the numerical simulations of ring melts, in Section 2.2 we provide details on the simulated systems, while in Section 2.3 we describe the algorithm used to reconstruct the primitive shape of each ring. Section 3 is devoted to recapitulate the relevant state-of-the-art theoretical ideas describing the structure of interacting random trees, focusing in particular on observables (Section 3.1) and distribution functions (Section 3.2). In Section 4 we present the main results, where we frame the extracted primitive shapes of rings in terms of the theoretical ideas for 2d melts of random trees. Then (Section 5), we discuss these results at the light of ring dynamics, highlighting in particular the fact that 2d rings represent a quite special situation with respect to 3d rings. Outlook and conclusions are presented in Section 6.

2 Polymer model, simulation methods, analysis of polymer conformations

2.1 Kinetic Monte Carlo algorithm

We adapt the efficient kinetic Monte Carlo (MC) algorithm employed in^31,35–37 to construct equilibrated dense melts of self-avoiding ring polymers on the two-dimensional triangular lattice (see Fig. 2). The lattice spacing, denoted by a, is chosen as the reference unit of length. The algorithm, based on the elastic lattice polymer model analogous to the Rubinstein's repton model,³⁸ works as the following. Two consecutive monomers along the chain can either sit on nearest neighbour lattice sites or occupy the same lattice site, with strictly no more than two monomers on the same site. At the same time, double occupancy of the same site by non-consecutive monomers is not permitted due to excluded volume. As such, the bond length b between nearest neighbor monomers takes the two possible values, = a or = 0, and in the latter case the bond is said to host a unit of stored length. For a generic ring polymer with N_ring bonds, the total contour length L_c ≡ N_ring〈b〉 < N_ringa where 〈b〉 is the average bond length. This numerical “trick” makes the polymer elastic.

Fig. 2 Two-dimensional illustration of the lattice polymer model and kinetic Monte Carlo moves. Monomers (filled dots) occupy the spatial positions of the regular triangular lattice (empty dots) and two nearest neighbor monomers are connected by a black line representing the polymer bond between them. For two nearest neighbor monomers occupying the same lattice site the bond that connects them (represented by an arc) represents a unit of “stored length”. Lattice positions connected by the double arrows are examples of allowed MC moves: (i) a unit of stored length unfolding to a normal bond, (ii) a bond folding into a unit of stored length, (iii) a Rouse-like move. Lattice positions connected by the double arrows with the cross are examples of forbidden MC moves: (iv) a monomer moving to the left that would occupy a lattice position with two monomers already present, (v) two non-nearest neighbor monomers on the same lattice site violating the excluded volume constraint.

To explore the effects of chain composition, we have considered polymer chains of different stiffnesses which have been modelled by introducing the following energy penalty term³¹,

for bending, where κ_bend represents the bending stiffness parameter and t⃑_i ≡ r⃑_i+1 − r⃑_i is the oriented bond vector between monomers i and i + 1 having spatial coordinates r⃑_i and r⃑_i+1.† It should be noticed that, since bond vectors are obviously ill-defined when two monomers form a stored length, the sum in eqn (1) is restricted to the effective bonds of the chains. By increasing κ_bend, the energy term eqn (1) makes polymers stiffer.

Ring polymers are let to evolve in time, and thus reach equilibrium, by proposing moves that obey the occupancy constraints described above. In more detail, a single move consists in randomly picking a monomer of one of the chains in the system and attempting its displacement towards one of the nearest lattice sites (see Fig. 2 for an illustration of possible situations that may occur). As in ref. 31, the move is then accepted based on the standard Metropolis–Hastings criterion³⁹ taking into account the energy term eqn (1) and with the additional constraints of maintaining chain connectivity and that the destination lattice site is either (1) empty or (2) occupied at most by one of the nearest neighbor monomers along the chain. In analogy with classical^37,40,41 polymer dynamics, case (1) is an example of Rouse-like move while case (2) is a reptation-like move (essentially the move produces mass drift along the contour length of the chain, as occurring in reptation dynamics). By these two moves only, the algorithm reproduces known^40,41 features of polymer dynamics and, thanks to the stored length “trick” which integrates local fluctuations of the chain density, remains efficient even when it is applied to the equilibration of very dense systems.³⁶

2.2 Simulation details

We have considered bulk solutions of M closed (ring) polymer chains, each chain made of N_ring monomers or bonds, with values N_ring × M = [40 × 1920, 80 × 960, 160 × 480, 320 × 240, 640 × 120, 1280 × 60], with constant total number of monomers = 76 800. Bulk conditions are implemented through the enforcement of periodic boundary conditions in a simulation box of total surface S, where the linear sizes of the box, equal to has been fixed based on the monomer number density that guarantees melt conditions. All these systems have been studied for bending stiffness parameters κ_bend = 0, 1, 1.5 (eqn (1)). For completeness, we report the measured MC acceptance rate as a function of the chain stiffness: 0.125 (κ_bend = 0), 0.021 (κ_bend = 1), 0.018 (κ_bend = 1.5). Notice that they are in line with the reported³¹ acceptance rates for 3d ring melts.

Each system is initialized by arranging the rings inside the simulation box in a conformation that avoids unphysical overlaps between monomers and nested states. To achieve this goal, the initial condition consisted of a melt of M randomly generated, linear, self-avoiding walks of total number of monomers N_ring/2. Then, we arrange ring monomers 1 and N_ring on the spatial position of one of the walk's ends, and continue by placing ring monomers 2 and N_ring − 1 on the next node of the walk, and so on until we get a perfectly, tightly double-folded ring conformation. Following this, we remove any residual overlap between ring monomers by performing a short kinetic Monte Carlo relaxation run. This procedure ensures that in the beginning configuration each ring is both self-avoiding and non-overlapping with the other rings of the melt while preserving local chain connectivity.

Then, for each MC step we pick a single monomer at random, displace it towards a random lattice neighbor and check if the move satisfies the constraints as described in Section 2.1. By taking as unit of “time” 1 MC sweep ≡ (N_ring×M) MC steps = 76 800 MC steps, we run simulations for ≃10⁷ up to ≃10⁸ MC sweeps. Notice that the length of each individual run is much larger than the characteristic relaxation time of the chain (see definition eqn (30) in Section 5), that guarantees complete equilibration of the chains (see Table S1 in SI). To graphically illustrate this point, we monitor chain relaxation to equilibrium in terms of the monomer mean-square displacement relative to the chain center of mass⁴²,

where . Asymptotically, g₂(t→∞) = 2〈R_g²〉_ring whereis the ring mean-square gyration radius. The appearance of a plateau in the large-time behavior of g₂ (see Fig. S1 in SI) is the signature that our chains have attained complete structural relaxation.

2.3 Reconstructing the primitive shapes of ring polymers: algorithmic definition

We reconstructed the primitive shapes of individual ring polymer configurations in a 2d melt configuration by introducing a numerical algorithm based on the original intuitions of Khokhlov and Nechaev⁷ and Rubinstein,^8,10 according to which each ring polymer behaves like being confined in a lattice of fixed obstacles (representing all the surrounding rings in the melt) and it can be mapped then onto a tree-like structure on the corresponding dual lattice where such obstacles are absent.

In our simulations polymers reside on the 2d triangular lattice of reference unit step = a, whose dual is the honeycomb lattice (panel (a) of Fig. 3) of unit step . For each ring in a typical melt configuration (panel (b) of Fig. 3), we identified first the honeycomb lattice sites residing within its contour which allow us to identify, without ambiguity, the nodes of the associated “backbone”. Then, we connect any two of these nodes of the honeycomb lattice by an edge whenever they are spatial nearest neighbors. Visual inspection on the obtained conformations (black lines in panel (c) of Fig. 3) confirms indeed a pronounced tree-like architecture.

Fig. 3 (a) Illustration of the 2d triangular lattice (orange dots, with 6 unit cells drawn explicitly) of unit step = a where ring polymers are simulated, along with the dual honeycomb lattice (blue dots, with 1 cell drawn explicitly) of unit step where the corresponding trees reside. (b) Snapshots of ring melt configurations for N_ring = 1280 and κ_bend = 0, 1, 1.5 (from left to right, see legend). (c) Examples of single ring polymers (blue/green/red circular contours) isolated from their corresponding melts, along with their primitive tree-like backbones (black lines inside each circular contour). Notice that the trees contain a certain amount of loops. (d) Final, loop-less, primitive shapes obtained by randomly removing one bond in each of the former loops. Shown on the l.h.s. of the figure, zoomed-in selected portions of tree conformations for κ_bend = 0 with loops and without loops after bonds removal.

If a ring was perfectly double-folded, the resulting tree backbone would be, by construction, without closed loops. In our simulated structures, a non-null amount of loops per ring is however always found (panel (c) of Fig. 3 and the zoomed-in portion for κ_bend = 0 on the l.h.s.). On the other hand, after easily verifying that such an amount remains small, essentially one may assume that loops do not compromise in any fundamental manner the statistical description of the trees, only they would make the analysis of tree conformations technically harder. For these reasons, we introduce the rule of removing loops simply by randomly cutting one bond of each of them (panel (d) of Fig. 3 and the zoomed-in portion for κ_bend = 0 on the l.h.s.). As anticipated the fraction of removed bonds is small, at most ≃8% for κ_bend = 1.5 of the original looped conformations (see Fig. S2 in SI).

2.4 Trees analysis

Each reconstructed tree is made of a total number of “honeycomb-lattice” nodes N_tree which, on average, is equal to ≃80% of the original number of “triangular-lattice” nodes N_ring of each ring (l.h.s. panel of Fig. S3 in SI).

As a consequence of the properties of the honeycomb lattice, every tree node may have one, two or, at most, three bonds protruding from it;‡ in the last case, we call that node a “branch-node”. Interestingly, we have found that the mean fraction of branch-nodes depends on the original ring stiffness κ_bend, with flexible rings being sensibly more branched than stiffer rings. Indeed, as shown in the r.h.s. panel of Fig. S3 in SI, the amount of branch-nodes saturates in the asymptotic limit of large rings (i.e., trees) from ≃22% for κ_bend = 0 up to ≃13% for κ_bend = 1.5.

After reconstructing the coordinates and bond connectivities of the ring primitive shapes, we have analyzed tree connectivity and spatial structure by using the “burning” algorithm for trees originally introduced in ref. 32. This algorithm works iteratively: each step consists in removing from the list of all tree nodes those with only one bond protruding from it and updating the number of bonds and the indices of the remaining ones accordingly. The algorithm stops when only one node (the so called “center” of the tree) remains in the list. In this way, by keeping track of the nodes that have been removed, it is possible to obtain information about the mass and the shape of tree branches (Section 3.1.1). The algorithm can be then generalized to detect also the minimal path length [small script l]_ij between any pair of nodes i and j (Section 3.1.1): it is in fact sufficient that both nodes “survive” the burning process. For a more detailed illustration of the burning algorithm and of its applications on trees, see ref. 32–34.

3 Structure of two-dimensional melts of branched polymers: theory and methods

In this Section, we introduce quantities and methods of analysis that will be used (in Section 4) to characterize the tree-like conformational properties of the ring primitive shapes.

3.1 Observables

3.1.1 Definitions. To explore the connection between the primitive shape of our 2d rings and the physics of branched polymers, we follow closely the works^32,33 and discuss the power-law behaviors for the following observables:

1. The mean path length as a function of the mean tree weight 〈N_tree〉§:

〈L〉 ∼ 〈N_tree〉^ρ, (4)

where and [small script l]_ij is the minimal path length connecting the pair of nodes i and j on a tree of weight =N_tree.

2. The mean branch weight as a function of the mean tree weight 〈N_tree〉:

〈N_br〉 ∼ 〈N_tree〉^ε, (5)

where N_br ≡ min(n, N_tree − 1 − n) is the weight of the smallest of the two sub-trees that are obtained by cutting one tree bond.

3. The mean-square end-to-end spatial distance as a function of the mean path length 〈L〉:

〈R_path²〉 ∼ 〈L〉^2ν_path (6)

4. The tree mean-square gyration radius as a function of the mean tree weight 〈N_tree〉:

〈R_g²〉 ∼ 〈N_tree〉^2ν, (7)

where and is the spatial position of the tree center of mass.

Notice the fundamental relations ρ = ε⁴⁴ and ν = ρν_path^32,33 holding between exponents. The complete list of the measured mean values with corresponding error bars for the different ring systems is presented in Table S2 in SI.

3.1.2 Extrapolation of exponents from data for finite-size trees. In order to get reliable values and errors of scaling exponents “ρ, ε, ν_path, ν” in the asymptotic (i.e., 〈N_tree〉 → ∞) limit, we adopt the following strategy, similar to the one employed in previous works,^32,33 which combines together the results obtained from fitting the 〈N_tree〉-dependent data to two functional forms.

To fix the ideas, take the data for the observable 〈R_g²〉 as a function of 〈N_tree〉 and the corresponding exponent ν and use the following expressions:

1. A simple power-law behavior with two (c, ν) fit parameters:

log(〈R_g²〉) = 2ν log(〈N_tree〉) + c, (8)

which we use only on the data with N_ring ≥ 320 (namely for the three largest rings/trees).

2. A power-law behavior with a correction-to-scaling term and four (c, d, Δ, ν) fit parameters:

which we use on the entire N_ring-range.

For the other exponents, analogous expressions to eqn (8) and (9) have been employed. In all cases, best fit parameters were obtained by standard χ²-minimization. The results of the two- (eqn (8)) and four-parameter (eqn (9)) fits are collected in Table S3 in SI. Overall the two procedures lead to similar results for the scaling exponents; we combine then the two results together as:

to obtain a final estimate and the corresponding error bar. See Table 1 for summary, and for comparison with the reference values of the exponents for 2d tree melts reported in ref. 33.

Table 1 Exponents characterizing tree observables (column 1) and relations between them (column 2, see ref. 33 for details). Column 3 reports the values for 2d tree melts estimated in ref. 33. Columns 4 to 6 report the estimated values for the tree-like primitive shapes of the 2d ring melts considered in this work, for different ring bending stiffness κ_bend. For technical details on the derivation of these exponents, see Section 3.1.2

	Relation to other exponents^33	Measured value (ref. 33)	Measured value (this work)
			κ bend = 0	κ bend = 1	κ bend = 1.5
〈L〉 ∼ 〈Ntree〉^ρ	—	0.613 ± 0.007	0.63 ± 0.04	0.61 ± 0.05	0.59 ± 0.07
〈Nbr〉 ∼ 〈Ntree〉^ε	ε = ρ	0.63 ± 0.01	0.64 ± 0.05	0.63 ± 0.04	0.63 ± 0.03
〈Rpath^2〉 ∼ 〈L〉^2νpath	—	0.780 ± 0.005	0.783 ± 0.006	0.79 ± 0.02	0.687 ± 0.006
〈Rg^2〉 ∼ 〈Ntree〉^2ν	ν = ρνpath	0.48 ± 0.02	0.50 ± 0.02	0.499 ± 0.003	0.461 ± 0.002

---

3.2 Distribution functions

3.2.1 Definitions. As originally pointed out in ref. 34, additional insight into the statistics of lattice trees can be gained by moving beyond the mere averages (5)–(7) and considering the related distribution functions that obey insightful scaling relations. Specifically, in this work we will focus on:

1. The distribution p_〈Ntree〉([small script l]) of path lengths [small script l]. Here, data from different 〈N_tree〉 are expected to superimpose, when expressed as functions of the scaled distances, x = [small script l]/〈L〉:

Moreover, the function q = q(x) is expected to be of the so called Redner-des Cloizeaux (RdC) form ref. 45 and 46,

q(x) = C_{[small script l]}x^{θ_{[small script l]}} exp(−(K_{[small script l]}x)^{t_{[small script l]}}), (12)

where the exponents (θ_{[small script l]}, t_{[small script l]}) depend on the trees universality class, while the numerical factors (C_{[small script l]}, K_{[small script l]}) are given by the analytical expressions,which can be found by imposing (i) that q(x) is normalized to 1 (i.e., ) and (ii) that the first moment is the only scaling length (i.e., ).

2. The distribution, p_〈Ntree〉(n), of branch sizes n. For 1 ≪n ≪ 〈N_tree〉, data from different trees are expected to superimpose and to display universal power-law behavior

p_〈Ntree〉(n) ∼ n^−(2−ε). (15)

3. The distribution p_〈Ntree〉(r⃑|〈L〉) of end-to-end spatial distances r⃑ of paths of length = 〈L〉 on trees of mass 〈N_tree〉. The data superimpose, when expressed as functions of the scaled distances, :

and the function q = q(x) is expected to be of the following RdC form:

q(x) = C_pathx^θ_path exp(−(K_pathx)^t_path). (17)

4. The distribution p_〈Ntree〉(r⃑) of vector distances r⃑ between tree nodes. Once again, data from different 〈N_tree〉 are expected to superimpose when expressed as functions of the scaled distances, :

and, again, the function q = q(x) is expected to be of the RdC form,

q(x) = C_treex^θ_tree exp(−(K_treex)^t_tree). (19)

The two sets of exponents (θ_path, t_path) (eqn (17)) and (θ_tree, t_tree) (eqn (19)) depend on the tree universality class, while the two sets of numerical factors (C_path, K_path) and (C_tree, K_tree) are given by the analytical expressions,obtained by imposing (i) that q(x) is normalized to 1 (i.e., ) and (ii) that the second moment is the only scaling length (i.e., ).

Notice that general scaling arguments³⁴ imply the following fundamental relations between the exponents of RdC functions and the exponents characterizing observables: θ_{[small script l]} = 1/ρ − 1, t_{[small script l]} = 1/(1 − ρ), t_path = 1/(1 − ν_path), θ_tree = min(θ_path, 0), t_tree = 1/(1 − ν). The only independent exponent is θ_path, which turns out to be >0 based on numerical estimates. Accordingly, θ_tree = 0. For a detailed characterization of lattice trees in terms of distribution functions and how these are related to RdC functions and their exponents, the reader is referred to ref. 34.

3.2.2 Extrapolation of exponents from data for finite-size trees. For each 〈N_tree〉, estimates of exponents and corresponding errors (θ_{[small script l]} ± Δθ_{[small script l]}, t_{[small script l]} ± Δt_{[small script l]}), (θ_path ± Δθ_path, t_path ± Δt_path) and (θ_tree ± Δθ_tree, t_tree ± Δt_tree) were obtained by best fit of RdC functions eqn (12) (with eqn (13) and (14)), eqn (17) and (19) (with eqn (20) and (21)) to, respectively, data sets for p_〈Ntree〉([small script l]), p_〈Ntree〉(r⃑|〈L〉) and p_〈Ntree〉(r⃑). Results are collected in Table S4 in SI.

Then, to extrapolate the values and errors in the asymptotic regime 〈N_tree〉 → ∞, we have fitted the straight line

y = −Ax + B (22)

to the data (1/〈N_tree〉, θ_{[small script l]}), and so on for the other exponents, for the three largest 〈N_tree〉's and with fit parameters A and B (see Fig. S4 in SI for a graphical illustration of the procedure). The intercept with the y-axis, B, gives the extrapolated value, while error estimate of the extrapolated value is obtained by repeating this procedure to data sets (1/〈N_tree〉, θ_{[small script l]} + Δt_{[small script l]}) and (1/〈N_tree〉, θ_{[small script l]} − Δt_{[small script l]}) and so on for the other exponents. Final results are summarized in Table 2. For comparison, we report also the measured values of the exponents for 2d tree melts published originally in ref. 34.

Table 2 Exponents characterizing the Redner-des Cloizeaux behavior of distribution functions (column 1) and relations to other exponents (column 2, see ref. 34 for details). Column 3 reports the values for 2d tree melts estimated in ref. 34. Columns 4 to 6 report the estimated values for the tree-like primitive shapes of the 2d ring melts considered in this work, for different ring bending stiffness κ_bend. For technical details on the derivation of these exponents, see Section 3.2.2

	Relation to other exponents^34	Measured value (ref. 34)	Measured value (this work)
			κ bend = 0	κ bend = 1	κ bend = 1.5
θ [small script l]		0.593 ± 0.003	0.421 ± 0.001	0.430 ± 0.002	0.413 ± 0.003
t [small script l]		2.35 ± 0.01	2.631 ± 0.006	2.605 ± 0.009	2.57 ± 0.01
θ path	—	0.63 ± 0.04	0.679 ± 0.001	0.711 ± 0.002	0.6602 ± 0.0004
t path		4.2 ± 0.1	6.132 ± 0.006	4.665 ± 0.003	2.94 ± 0.02
θ tree	min(θpath, 0) = 0	−0.14 ± 0.02	−0.09 ± 0.009	−0.07 ± 0.01	−0.06 ± 0.01
t tree		1.857 ± 0.005	1.75 ± 0.02	1.72 ± 0.03	1.65 ± 0.02

---

4 Results

4.1 Observables

We start first with the mean path length 〈L〉 as a function of the mean tree weight 〈N_tree〉 (eqn (4)), and determine the exponent ρ by using the methods of Section 3.1.2. The results for 〈L〉 vs. 〈N_tree〉 for different κ_bend parameters are shown in panel (a) of Fig. 4 (symbols) together with the estimated values and statistical errors for ρ (dashed lines). In particular, the behavior is in very good agreement with the results for 2d tree melts³³ (Table 1). Interestingly, by introducing the asymptotic (i.e., large-〈N_tree〉) fraction of branch-nodes Λ (see r.h.s. panel of Fig. S3 in SI), its inverse Λ⁻¹ corresponds to the mean path length between two branch-nodes⁴⁷ and, so, to the natural length scale beyond which branching becomes effectively relevant. As a consequence, data for individual κ_bend collapse onto a single master curve (Fig. 4(a), inset) by rescaling the x- and y-coordinate as the following:where ρ = 0.61 in eqn (24) is the average¶ of the estimated best values for individual κ_bend (Table 1).

Fig. 4 Conformational properties of trees: observables (symbols) and asymptotic power-law behaviors (dashed lines). (a) 〈L〉 ∼ 〈N_tree〉^ρ, mean path length as a function of the mean tree weight 〈N_tree〉. (b) 〈N_br〉 ∼ 〈N_tree〉^ε, mean branch weight as a function of the mean tree weight 〈N_tree〉. (c) 〈R_path²〉 ∼ 〈L〉^2ν_path, mean-square end-to-end spatial distance of paths of length = 〈L〉. (d) 〈R_g²〉 ∼ 〈N_tree〉^2ν, mean-square gyration radius as a function of the mean tree weight 〈N_tree〉. In each panel the dashed lines express the interval of possible values for the exponent, lying between the minimum lower-bound and the maximum upper-bound for all κ_bend (see Table 1). (Inset) Universal scaling plots for the observables. Here, the value of the exponent used in each plot corresponds to the average of the estimated best values for individual κ_bend (see Table 1).

Then, we consider the mean branch weight 〈N_br〉 as a function of 〈N_tree〉 (eqn (5)), and determine the related exponent ε, see panel (b) of Fig. 4 (symbols and dashed lines). As before, our estimates for ε for individual κ_bend are in agreement with each other and with the corresponding measured values for 2d tree melts³³ (Table 1); moreover, it is important to emphasize that our estimated values satisfy also the relation ε = ρ with great accuracy. Finally, data for individual κ_bend collapse onto a single master curve (Fig. 4(b), inset) by rescaling the x-coordinate according to eqn (23) and, analogously to eqn (24), the y-coordinate as the following:

where ε = 0.63 in eqn (25) is the average of the estimated best values for individual κ_bend (Table 1).

Third, we consider the mean-square end-to-end spatial distance 〈R_path²〉 as a function of the mean path length 〈L〉 (eqn (6)), and determine the related exponent ν_path, see panel (c) of Fig. 4 (symbols and dashed lines). Similarly to previous results, the estimated ν_path for κ_bend = 0 and κ_bend = 1 agree well with each other and with the measured value for 2d tree melt (Table 1); instead, the estimated value for κ_bend = 1.5 now is ≈13% smaller. We argue that this discrepancy may be due to the limited amount of branching observed for κ_bend = 1.5 compared to the structures in the other two cases (r.h.s. panel of Fig. S3 in SI): in fact, for Λ → 0 one recovers the linear chain limit where ν_path = ν = 1/2. At the same time, the reported difference between the value of ν_path for κ_bend = 1.5 and the other values is not large, so we decide to ignore it in the rest of the discussion and proceed with the analysis in the same manner as we did for the other observables. In this regard, we notice that data for individual κ_bend are already on top of each other. Nonetheless, it is useful to rescale the x- and y-coordinates as the following:

〈R_path²〉 → 〈R_path²〉〈L〉^−2ν_path, (27)

where ν_path = 0.75 in eqn (27) is the average of the estimated best values for individual κ_bend (Table 1). The result (see inset of Fig. 4(c)) confirms the universal behavior and suggests that the mean-square end-to-end distance between two branch-nodes is ≃(1/Λ)^2ν_path.

To conclude, we examine the mean-square gyration radius 〈R_g²〉 as a function of 〈N_tree〉 (eqn (7)), and determine the related exponent ν, see panel (d) of Fig. 4 (symbols and dashed lines). Also in this last case, the estimated ν for κ_bend = 0 and κ_bend = 1 agree well with each other and with the measured value for 2d tree melt (Table 1); again the estimated value for κ_bend = 1.5 is slightly (≈8%) smaller, but we disregard this small difference in the analysis. Overall, we confirm the general scaling relation ν = ρν_path relating the different exponents. Finally, based on the results for the previous observables, we expect data for individual κ_bend to collapse on top of each other assuming eqn (23) for the rescaling of the x-coordinate and

for the y-coordinate, and where ν = 0.49 in eqn (28) is the average of the estimated best values for individual κ_bend (Table 1). The plot in the inset of Fig. 4(d) confirms our scaling ansatz.

To summarize this first part, the analysis of the primitive shapes of ring polymers in 2d melt by means of the same fundamental observables used to study melts of randomly branching polymers has shown that these latter systems are fundamentally analogous to the former. In the next Section, we generalize the analysis to the distribution functions.

4.2 Distribution functions

We start with the path length distribution functions, p_〈Ntree〉([small script l]). As illustrated in Fig. 5, the measured p_〈Ntree〉([small script l]) for different tree sizes fall onto universal master curves, when plotted as a function of the rescaled path length x = [small script l]/〈L〉 (eqn (11)), which are well described by the RdC form eqn (12)–(14) (dashed lines). Noticeably, the estimated values of the asymptotic exponents (θ_{[small script l]}, t_{[small script l]}) (see Section 3.2.2 for details) for the different chain flexibilities κ_bend are in good agreement with each other and relatively close to the measured values³⁴ for trees in 2d melts (see Table 2).

Fig. 5 Conformational properties of trees: p_〈Ntree〉([small script l]), distribution functions of linear paths of length [small script l]. The dashed line of each panel corresponds to the predicted RdC functional form for trees, eqn (12)–(21), with the exponents θ_{[small script l]} and t_{[small script l]} as in Table 2. Data of different colors denote different mean tree weight 〈N_tree〉 (see legend).

Next, we consider the distribution functions of branch weight n, p_〈Ntree〉(n). Data for different tree sizes 〈N_tree〉 are in very good agreement (Fig. 6) with the predicted (eqn (15)) power-law behavior p_〈Ntree〉(n) ∼ n^−(2−ε) for trees (dashed lines, ε is as in panel (b) of Fig. 4).

Fig. 6 Conformational properties of trees: p_〈Ntree〉(n), distribution functions of branch weight n. The dashed lines express the interval of possible values for the exponent ε (see Fig. 4(b)), lying between the minimum lower-bound and the maximum upper-bound for all κ_bend (see Table 1). Colorcode/symbols are as in Fig. 5.

Then, we consider the distribution functions p_〈Ntree〉(r⃑|〈L〉) of end-to-end spatial distances r⃑ of linear paths of length =〈L〉. When plotted as a function of the rescaled distance, , the measured curves for different κ_bend fall onto single universal master curves (Fig. 7). Here, however, and contrarily to the results presented so far, the master curves for the different stiffnesses agree less well with each other. This, in particular, can be seen by the fact that the asymptotic RdC ansatz (dashed lines) are characterized by quite distinct values for the exponents θ_path and, especially, t_path (see Table 2).

Fig. 7 Conformational properties of trees: p_〈Ntree〉(r⃑|〈L〉), distribution functions of end-to-end spatial distances r⃑ of paths of length =〈L〉. The dashed line of each panel corresponds to the RdC functional form for trees, eqn (17) with eqn (20) and (21), with the exponents θ_path and t_path as in Table 2. Colorcode/symbols are as in Fig. 5.

Finally, we consider the distribution functions p_〈Ntree〉(r⃑) of vector distances r⃑ between tree nodes. By rescaling the x-axis by the square-root of the corresponding momenta curves from different rings obey universal behavior (Fig. 8). This behavior is well described (dashed lines) by the RdC function eqn (19), with estimated values for the exponents θ_tree and t_tree for the different flexibilities κ_bend in good agreement with each other and close to the estimated values for 2d tree melts (Table 2).

Fig. 8 Conformational properties of trees: p_〈Ntree〉(r⃑), distribution functions of vector distances r⃑ between tree nodes. The dashed line of each panel corresponds to the RdC functional form for trees, eqn (19) with eqn (20) and (21), with the exponents θ_tree and t_tree as in Table 2. Colorcode/symbols are as in Fig. 5.

In conclusion, the results for both observables and distribution functions which have been presented so far demonstrate (besides some discrepancies detected in the distribution functions p_〈Ntree〉(r⃑|〈L〉) for node-to-node spatial distances of linear paths) that the tree-like primitive paths of ring polymers in 2d melts obey the same statistics of 2d melts of randomly branching polymers. We now explore the connection between these two systems in relation to polymer dynamics.

5 A note on ring dynamics

Although it is not the main purpose of this work, we present here a few results on ring dynamics as well. In fact, as originally pointed out by Rubinstein and coworkers^8,10 and reaffirmed at various degrees of sophistications in later works,^16,48–50 fingerprints of the tree-like structure of rings in melt can be found in ring dynamics.

In particular, in ref. 16 it was suggested that the relaxation time of the chains as a function of the monomers number N_ring is given by the power-law behavior

τ_relax = τ_micrN^{2+(1−θ)ρ+θν}_ring, (29)

where τ_micr is a “microscopic” time scale of the order of the time needed for the chain to relax between pairs of neighboring topological obstacles. The ρ and ν exponents in eqn (29) are the same describing the behaviors of observables 〈L〉 (eqn (4)) and 〈R_g²〉 (eqn (7)) as functions of N_ring, while the parameter 0 ≤ θ ≤ 1 (not to be confused with the various exponents θ associated to the RdC functions) quantifies the postulated phenomenon of so called tube dilation: in analogy to linear chains,^40,41 the array of topological obstacles where ring polymers relax can be assimilated to the relaxation inside a tube-like region shaped by the same obstacles. Assuming such obstacles as fixed,⁴⁸ the tube is a static structure and θ = 0; vice versa,¹⁶ assuming that chain relaxation proceeds self-similarly at all scales (as in ordinary Rouse model⁴⁰), the tube “dilate” with time and θ = 1. As shown in ref. 16, eqn (29) applies well to the case of melt of rings in 3d with θ = 1, i.e. in good agreement with the mechanism of complete tube dilation. What about the case of 2d melts, which we should expect to be very different from the 3d case because different chains can not, for instance, thread each other?

To answer this question, we have computed τ_relax for our rings by defining it as the time scale for the entire ring to diffuse a distance equal to its mean-size, i.e.

g₃(t = τ_relax) = 〈R_g²〉_ring, (30)

where

g₃(t) ≡ 〈(r⃑_cm(t) − r⃑_cm(0))²〉 (31)

is the mean-square displacement of the chain centre of mass. Results for τ_relax, normalized by N_ring², as a function of N_ring are shown in the top panel of Fig. 9. Interestingly the data suggest the scaling τ_relax ∼ N_ring², therefore they are not compatible with the power-law (29) whose exponent is, by construction, strictly >2. To validate this result further, we consider the scaling behavior of the chain diffusion coefficient,for which we expect D_ring ∼ 〈R_g²〉/τ_ring ∼ N^−2(1−ν)_ring = N_ring⁻¹. We have estimated diffusion coefficients by best fits of data for g₃(t) to a simple linear function (eqn (32)), with the fit procedure restricted to the time interval t > 5τ_relax (see Fig. S5 in SI). Results for D_ring, each multiplied by the corresponding N_ring, are shown in the bottom panel of Fig. 9 for each stiffness κ_bend, and they confirm the scaling analysis.

Fig. 9 Ring dynamics: (top) Ring relaxation time normalized by the square ring monomer number, τ_relax/N_ring², as a function of N_ring. (bottom) Ring diffusion coefficient normalized by the inverse ring monomer number, D_ring/N_ring⁻¹, as a function of N_ring. Colorcode/symbols are for different chain stiffness κ_bend (see legend, as in Fig. 4). Dotted lines joining the symbols serve as a guide for the eye. Upraising of the data for N_ring = 1280 is due to the limited length of the MC trajectory (see Fig. S5 in SI).

Curiously, we notice that the power-law behaviors τ_relax ∼ N_ring² and D_ring ∼ N_ring⁻¹ reported here are the same that one would have expected for an ordinary Rouse chain in ideal (i.e., no volume interactions) conditions.^40,41 For the polymers considered here this is obviously not the case, and the agreement with the Rouse model must be considered a fortuituous coincidence for which, at the moment, we have no physical intuition. We just conclude by highlighting the fact that the two power-laws are no artefact of the present lattice model, as they were also reported in previous molecular simulations of 2d ring melts.⁵¹

6 Summary and conclusions

In this work we have performed a detailed analysis of the conformational properties of unknotted and non-concatenated ring polymers in 2d melt conditions (Fig. 3), aiming at demonstrating that rings double-fold around randomly branching (tree-like) conformations as originally proposed by Khokhlov, Nechaev, Rubinstein and others.^7,8,10 In particular, after reconstructing the “primitive shapes” of ring conformations from Monte Carlo computer simulations of a generic polymer model on the two-dimensional triangular lattice, we proceeded with an analysis of said branched structures by means of the same observables and distribution functions (Section 3) used for the analysis of lattice trees.

We thus confirm (Section 4, in particular Fig. 4–8) that the conformational properties of the branched structures are, in general, in good agreement with the measured properties^33,34 of two-dimensional branched polymers in melt conditions. At the same time, a closer look reveals also interesting quantitative differences: for instance, some of the exponents which characterize the distribution functions (for instance t_path, see Table 2) show larger discrepancies between themselves and the corresponding predicted values for trees than the exponents for observables (Table 1). We argue that these deviations from tree behavior emerge because our rings are not perfectly double-folded⁵² around their “primitive shape” (see Section 2.3), and then the “mapping” between rings and trees works only approximately. In this regard, it is worth remarking that the analysis of ring dynamics (Section 5) also suggests a rather different behavior than the one expected based on the branched shapes of the rings. In fact, the rings' relaxation time is found to be faster than the theoretical predictions, scaling “only” quadratically with the rings' mass and, therefore, surprisingly close to the result expected for ordinary Rouse motion.

We conclude with an outlook for future work. Besides trying to understand the origin of the reported deviations between rings and the tree description, it would be even more interesting (and challenging) to generalize the analysis of this work to three-dimensional melts of rings. As a matter of fact, how rings fold in 3d is still not fully understood. In fact, while on one side a multi-scale numerical protocol based on the random-tree model²⁹ reveals an accurate similarity between some of the trees' and rings' static properties, on the other the “tree picture” cannot be considered complete since, for example, one of two nearby rings may easily penetrate onto an open loop of the other without violation of topological constraints. Although rare, such so-called threading events^{30,31,53–55} have been observed and are thought to be responsible of some exotic properties in rings' out-of-equilibrium behavior.^56–60 For these reasons, developing an algorithm for the reconstruction of the chain's “primitive shape” along the lines of Section 2.3 would represent an important contribution towards the comprehension of these fundamental systems.

Author contributions

MAU and AR contributed equally to this work.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article are included in the supplementary information (SI). Supplementary information provides: the data supporting this article and some complementing figures. See DOI: https://doi.org/10.1039/d5sm00947b.

Acknowledgements

MAU is supported by an EMBO Postdoctoral Fellowship (ALTF 1305-2024). AR acknowledges S. K. Nechaev for insightful discussions, and financial support from PNRR Grant CN_00000013_CN-HPC, M4C2I1.4, spoke 7, funded by Next Generation EU.

References

1. S. F. Edwards, Proc. Phys. Soc., 1967, 92, 9–16.
2. P. G. de Gennes, J. Chem. Phys., 1971, 55, 572–579.
3. R. Everaers, S. K. Sukumaran, G. S. Grest, C. Svaneborg, A. Sivasubramanian and K. Kremer, Science, 2004, 303, 823.
4. D. J. Read, K. Jagannathan and A. E. Likhtman, Macromolecules, 2008, 41, 6843–6853.
5. N. Uchida, G. S. Grest and R. Everaers, J. Chem. Phys., 2008, 128, 044902.
6. R. Everaers, H. A. Karimi-Varzaneh, F. Fleck, N. Hojdis and C. Svaneborg, Macromolecules, 2020, 53, 1901–1916.
7. A. R. Khokhlov and S. K. Nechaev, Phys. Lett. A, 1985, 112, 156–160.
8. M. Rubinstein, Phys. Rev. Lett., 1986, 57, 3023–3026.
9. S. K. Nechaev, A. N. Semenov and M. K. Koleva, Phys. A, 1987, 140, 506–520.
10. S. P. Obukhov, M. Rubinstein and T. Duke, Phys. Rev. Lett., 1994, 73, 1263–1266.
11. M. G. Brereton and T. A. Vilgis, J. Phys. A: Math. Gen., 1995, 28, 1149–1167.
12. A. L. Kholodenko and T. A. Vilgis, Phys. Rep., 1998, 298, 251–370.
13. T. Sakaue, Phys. Rev. Lett., 2011, 106, 167802.
14. S. Obukhov, A. Johner, J. Baschnagel, H. Meyer and J. P. Wittmer, EPL, 2014, 105, 48005.
15. A. Y. Grosberg, Soft Matter, 2014, 10, 560–565.
16. T. Ge, S. Panyukov and M. Rubinstein, Macromolecules, 2016, 49, 708–722.
17. F. Ferrari, J. Paturej, M. Piatek and Y. Zhao, Nucl. Phys. B, 2019, 945, 114673.
18. M. Kapnistos, M. Lang, D. Vlassopoulos, W. Pyckhout-Hintzen, D. Richter, D. Cho, T. Chang and M. Rubinstein, Nat. Mater., 2008, 7, 997–1002.
19. R. Pasquino, T. C. Vasilakopoulos, Y. C. Jeong, H. Lee, S. Rogers, G. Sakellariou, J. Allgaier, A. Takano, A. R. Brás, T. Chang, S. Gooßen, W. Pyckhout-Hintzen, A. Wischnewski, N. Hadjichristidis, D. Richter, M. Rubinstein and D. Vlassopoulos, ACS Macro Lett., 2013, 2, 874–878.
20. M. Kruteva, J. Allgaier and D. Richter, Macromolecules, 2017, 50, 9482–9493.
21. M. Kruteva, J. Allgaier, M. Monkenbusch, L. Porcar and D. Richter, ACS Macro Lett., 2020, 9, 507–511.
22. M. Kruteva, M. Monkenbusch, J. Allgaier, O. Holderer, S. Pasini, I. Hoffmann and D. Richter, Phys. Rev. Lett., 2020, 125, 238004.
23. M. Kruteva, J. Allgaier and D. Richter, Macromolecules, 2023, 56, 7203–7229.
24. M. E. Cates and J. M. Deutsch, J. Phys. France, 1986, 47, 2121–2128.
25. M. Müller, J. P. Wittmer and M. E. Cates, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1996, 53, 5063–5074.
26. M. Müller, J. P. Wittmer and M. E. Cates, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 61, 4078–4089.
27. J. D. Halverson, W. B. Lee, G. S. Grest, A. Y. Grosberg and K. Kremer, J. Chem. Phys., 2011, 134, 204904.
28. J. D. Halverson, W. B. Lee, G. S. Grest, A. Y. Grosberg and K. Kremer, J. Chem. Phys., 2011, 134, 204905.
29. A. Rosa and R. Everaers, Phys. Rev. Lett., 2014, 112, 118302.
30. R. D. Schram, A. Rosa and R. Everaers, Soft Matter, 2019, 15, 2418–2429.
31. M. A. Ubertini, J. Smrek and A. Rosa, Macromolecules, 2022, 55, 10723–10736.
32. A. Rosa and R. Everaers, J. Phys. A: Math. Theor., 2016, 49, 345001.
33. A. Rosa and R. Everaers, J. Chem. Phys., 2016, 145, 164906.
34. A. Rosa and R. Everaers, Phys. Rev. E, 2017, 95, 012117.
35. V. Hugouvieux, M. A. V. Axelos and M. Kolb, Macromolecules, 2009, 42, 392–400.
36. R. D. Schram and G. T. Barkema, J. Comput. Phys., 2018, 363, 128–139.
37. M. A. Ubertini and A. Rosa, Phys. Rev. E, 2021, 104, 054503.
38. M. Rubinstein, Phys. Rev. Lett., 1987, 59, 1946–1949.
39. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087–1092.
40. M. Doi and S. F. Edwards, The Theory of Polymer Dynamics, Clarendon, Oxford, 1986.
41. M. Rubinstein and R. H. Colby, Polymer Physics, Oxford University Press, New York, 2003.
42. K. Kremer and G. S. Grest, J. Chem. Phys., 1990, 92, 5057–5086.
43. P. H. W. van der Hoek, A. Rosa and R. Everaers, Phys. Rev. E, 2024, 110, 045312.
44. E. J. Janse van Rensburg and N. Madras, J. Phys. A: Math. Gen., 1992, 25, 303–333.
45. S. Redner, J. Phys. A: Math. Gen., 1980, 13, 3525–3541.
46. J. des Cloizeaux and G. Jannink, Polymers in Solution, Oxford University Press, Oxford, 1989.
47. M. Daoud and J. F. Joanny, J. Phys. France, 1981, 42, 1359–1371.
48. J. Smrek and A. Y. Grosberg, J. Phys.: Condens. Matter, 2015, 27, 064117.
49. E. Ghobadpour, M. Kolb, M. R. Ejtehadi and R. Everaers, Phys. Rev. E, 2021, 104, 014501.
50. E. Ghobadpour, M. Kolb, I. Junier and R. Everaers, Macromolecules, 2025, 58, 10632–10646.
51. J. Kim, J. M. Kim and C. Baig, Soft Matter, 2021, 17, 10703–10715.
52. A. Rosa and R. Everaers, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 1–12.
53. D. Michieletto, Soft Matter, 2016, 12, 9485–9500.
54. J. Smrek and A. Y. Grosberg, ACS Macro Lett., 2016, 5, 750–754.
55. J. Smrek, K. Kremer and A. Rosa, ACS Macro Lett., 2019, 8, 155–160.
56. D. Michieletto and M. S. Turner, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 5195–5200.
57. Q. Huang, J. Ahn, D. Parisi, T. Chang, O. Hassager, S. Panyukov, M. Rubinstein and D. Vlassopoulos, Phys. Rev. Lett., 2019, 122, 208001.
58. T. C. O'Connor, T. Ge, M. Rubinstein and G. S. Grest, Phys. Rev. Lett., 2020, 124, 027801.
59. J. Smrek, I. Chubak, C. N. Likos and K. Kremer, Nat. Commun., 2020, 11, 26.
60. C. Micheletti, I. Chubak, E. Orlandini and J. Smrek, ACS Macro Lett., 2024, 13, 124–129.

---

Footnotes

† For ring polymers, it is implicitly assumed the periodic boundary condition along the chain N + 1 ≡ 1.

‡ As noticed in works,^32–34,43 there is no substantial loss of generality by choosing a model for branched polymers where nodes host no more than three bonds compared to models featuring nodes with more than three bonds protruding from them.

§ In the mentioned works^32,33N_tree is not a fluctuating quantity. Here, as shown in the l.h.s. panel of Fig. S3 in SI, the relative fluctuations around the mean value 〈N_tree〉 become rapidly small in the large-N_ring limit, so they are expected to have a negligible effect on the estimation of scaling exponents.

¶ The excellent collapse in the inset of Fig. 4(a) and the agreement between the measured ρ-values for individual κ_bend justifies a posteriori our choice for taking the average. In the rest of this work, we do the same also for the other exponents.

---