Monte Carlo simulations of discrete Rouse dynamics on a 2D lattice: emergence of global behavior in a polymer chain from local constraints

Monte Carlo Simulations of Discrete Rouse Dynamics on a 2D Lattice Emergence of Global Behavior in a Polymer Chain from Local Constraints Arpan Dey M1 Physics (IDIL) Université de Montpellier Acknowledgments I would like to express my sincere thanks to Dr. Jean-Charles Walter, Prof. Andrea Parmeggiani and Linda Delimi for their guidance and support in the preparation and refinement of this report. This report was prepared during my bibliographic studies in advance of my internship at the Laboratoire Charles Coulomb (L2C), Université de Montpellier. Table of Contents Abstract 1 The Rouse Model 2 Monte Carlo Implementation of Discrete Rouse Dynamics 3 Polymer Dynamics under Heterogeneous Geometric Rules Discussion and Conclusion References Appendix: End-Monomer Interactions as a Minimal Model for Polymer Collapse Abstract A polymer, in simple words, is a connected chain of individual units called monomers. Unlike gas particles, polymers are connected chains, which makes the modeling of polymer interactions more difficult – each monomer is connected to the monomers immediately before and after it along the chain, which puts constraints on the ways the monomers (and hence the polymer) can move, even without considering external effects from the surrounding. The Rouse model provides a minimal theoretical framework for studying polymer dynamics. It describes monomer motion as harmonically coupled Brownian oscillators (restoring force plus random thermal noise!). In this report, we explore a discrete, lattice-based analogue of the Rouse model by enforcing a single local constraint – preservation of the nearest-neighbor distances along the polymer chain – on a two-dimensional lattice, along with self-avoidance. Using Monte Carlo simulations, we study the evolution of the chain under this constraint. We quantify this behavior by tracking the squared end-to-end distance 𝑅𝑒 2, the squared radius of gyration 𝑅𝑔 2 and the mean squared displacement ⟨|𝑟𝑖(𝑡)−𝑟𝑖(0)|2〉, the latter exhibiting a characteristic subdiffusive (~ 𝑡1/2) to diffusive (~ 𝑡) crossover at a timescale of the order of 𝑁2 if we drop self-avoidance (𝑁 is the number of monomers). This observation is consistent with Rouse scaling. Despite its simplicity and lack of explicit energy terms, the discrete distancepreserving rule encapsulates essential aspects of polymer dynamics, illustrating how complex global behavior can emerge from minimal microscopic constraints. To further explore the behavior of polymers on a 2D lattice beyond Rouse dynamics, we introduce a set of alternating-rule toy models (alternating copolymers) and a block copolymer toy model to illustrate how simple, spatially heterogeneous geometric constraints can drive the system away from detailed balance. By modifying only the kinds of moves different monomers are permitted to attempt – without introducing forces, potentials or energetic biases – we create a family of simple, nonequilibrium models that generate diverse, global behavior. In some models, deviations from standard Rouse dynamics are very weak, because geometric constraints suppress the heterogeneous update rules; others produce strong nonequilibrium expansion of polymers, driven by fluctuations in the bond lengths. Together, these results show that nonuniform geometric rules can encode a wide spectrum of dynamical responses while retaining conceptual simplicity. Overall, this report demonstrates that complex global phenomena in polymer dynamics can emerge from minimal local constraints, and that lattice models provide an instructive platform for dissecting how equilibrium and nonequilibrium features arise from purely geometric principles. Computational note: The Monte Carlo simulations performed involve a large number of update attempts (up to the order of 107) per run and large ensembles (up to 100 independent realizations per chain length). Executing this serially, either locally or on browser-based platforms such as Google Colab, becomes extremely slow and impractical due to single-core execution and process overhead. Thus, the simulations were run locally using Python installed via the Miniforge Conda distribution, and parallelized using Python’s multiprocessing module, with each run executed on a separate CPU core. This approach preserves identical Monte Carlo dynamics while greatly reducing total runtime. The Rouse Model The Rouse model provides one of the simplest descriptions of polymer dynamics. It treats a polymer chain as a sequence of 𝑁 monomers connected by harmonic springs. Each monomer is immersed in a viscous medium and experiences two opposing influences: friction from the surrounding fluid and random thermal kicks due to molecular collisions. The combination of these forces gives rise to stochastic motion of the entire chain (Brownian oscillators!). In its continuum form, the time evolution of the position 𝑟𝑖(𝑡) of the 𝑖-th monomer is given by the equation: 𝜁𝑑𝑟𝑖 𝑑𝑡 =𝑘(𝑟𝑖+1+𝑟𝑖−1−2𝑟𝑖)+𝜉𝑖(𝑡) where 𝜁 is the friction coefficient, 𝑘 is the spring constant, and 𝜉𝑖(𝑡) represents the random thermal forces with zero mean. The term involving the neighboring monomers (𝑟𝑖−1 and 𝑟𝑖+1) expresses a harmonic restoring force – each monomer behaves like a coupled oscillator, connected to its neighbors, while still subject to random thermal motion. The Rouse model originally assumes the polymer to be an ideal Gaussian chain. In simple words, this means we allow more than one monomer to “overlap” – they do not repel each other and do not collide (no volume exclusion). Thus, even two distant parts of the chain are allowed – mathematically speaking – to pass arbitrarily close or overlap in space. The Rouse equations describe a long phantom chain that can freely crumple and fold without getting in its own way – something that is, of course, not true for real polymers but still accurately captures many large-scale dynamical properties. Under this assumption, each step of the polymer (basically the vector from monomer 𝑖 to monomer 𝑖+1) – can treated as an independent random step. There is no preferred direction and no correlation between different steps. Thus, the ideal polymer just follows the statistics of a random walk. If we take 𝑁 such steps, each of length 𝑎, the end-to-end vector of this 𝑁segment behaves like a sum of independent random vectors. We define the steps using the step vector Δ𝑟𝑖, 𝑖=1,2…𝑁, each vector of fixed size |Δ𝑟𝑖|=𝑎. The average of each step is zero and steps are uncorrelated, hence we estimate the square of the size of the segment. As a reminder, we do not assume self-avoidance, but we assume the steps are independent of direction and identical. Then the end-to-end vector would be (please mind the notation!): Δ𝑟=∑Δ𝑟𝑖 𝑁 𝑖=1 Since we have assumed isotropy and independence of steps, we expect the average displacement to be zero: ⟨Δ𝑟𝑖〉=0. We also have: ⟨Δ𝑟𝑖.Δ𝑟𝑗〉=𝑎2𝛿𝑖𝑗 𝛿𝑖𝑗 is the Kronecker delta, which means the mean of the dot product of 𝑟𝑖 and 𝑟𝑗 is always zero, except when 𝑖=𝑗 (when we are talking about the same step, or the same monomer on the lattice), in which case their dot product yields 𝑎2, as expected. This does not mean that when 𝑖≠𝑗, the 𝑖-th and 𝑗-th steps are orthogonal – it simply means there is no correlation between them over a large number of steps (since we are considering the mean of the dot product). Given that we have assumed isotropy and independence, this makes sense. Now, keeping all that in mind, let us calculate the mean squared end-to-end distance of the chain: ⟨Δ𝑟2〉=⟨Σ𝑖Δ𝑟𝑖.Σ𝑗Δ𝑟𝑗〉=Σ𝑖,𝑗⟨Δ𝑟𝑖.Δ𝑟𝑗〉=Σ𝑖𝑎2=𝑎2𝑁 This is characteristic of Gaussian polymers: the variance of the end-to-end vector ⟨Δ𝑟2⟩ grows linearly with the number of steps 𝑁. This means an ideal polymer does not stretch out stiffly; it meanders. If you zoom in on any chunk of monomers, that chunk performs a miniature random walk of its own. Doubling the number of steps 𝑁 makes the mean-squared size twice as large, and so on; and so the typical linear size grows only as √𝑁. This square-root scaling essentially sets the natural length scale for any internal motion of the chain. If some portion of the polymer has had time to internally relax (lose correlations with its neighbors), the spatial fluctuations it explores are of the order Δ𝑟(𝑁)∼𝑁1/2, not 𝑁 itself. The Gaussian viewpoint makes the Rouse model solvable and underpins all its scaling laws. It is important to keep in mind that in this picture, the chain’s internal “wiggles” have the same statistical structure as a random walk, and the dynamics simply govern how fast different parts of that walk can rearrange. In this work, we reduce the continuous Rouse model to a discrete, lattice-based version that still preserves its essential physical idea. Instead of explicit harmonic forces or stochastic differential equations, we enforce a single geometric rule – each monomer can move only in ways that preserve the distance to its nearest neighbors along the chain. This discrete constraint effectively replaces the harmonic potential with a distance-preserving rule. Further, we also drop the assumption that monomers can overlap, and incorporate volume exclusion in our model. Using Monte Carlo simulations on a 2D lattice, we study how such purely local moves give rise to emergent large-scale behaviors consistent with Rouse dynamics, including relaxation from stretched and random initial configurations, fluctuations in polymer size (through both the end-to-end distance and the radius of gyration) and the subdiffusive-todiffusive crossover in monomer motion. We can move the polymer chain (we are talking about a Gaussian polymer now!) in many different ways – stretch it, wiggle it in the middle, wiggle only one end – and the small, local disturbances fade away quickly, whereas the larger, slower disturbances spread throughout the polymer and take time to die out. Thus, there are many relaxation times for the different kinds of motion, and the slowest of these relaxation times grow as 𝑁2, where 𝑁 is the length of the polymer (the number of monomers): 𝜏∝𝑁2 This is simply the order of time it takes for the polymer to completely “forget” its initial configuration. This is visualized by plotting ⟨|𝑟𝑖(𝑡)−𝑟𝑖(0)|2〉 – the mean squared displacement of an arbitrary monomer 𝑖 (no monomer is special in the chain, so this monomer is a representative of the average monomer) – over time. Even though we are tracking a single monomer, this gives a measure of how much the polymer has moved from its initial configuration (time 𝑡=0) at some later time 𝑡=𝑡. At short times, the monomers have a strong memory of the initial configurations, and the motion of a single monomer is subdiffusive – the polymer cannot diffuse freely since the motion of the monomers are heavily constrained by the connections to the neighboring monomers. The mean square displacement, in this regime, scales as: ⟨Δ𝑟2(𝑡)⟩∼𝑡1/2 Over a sufficiently long period of time, the polymer as a whole diffuses freely: ⟨Δ𝑟2(𝑡)⟩∼𝑡 The crossover from the subdiffusive to diffusive regime takes place at a timescale of about 𝑁2. First, we understand why ⟨Δ𝑟2(𝑡)⟩ scales as 𝑡1/2 for small times. We focus on an arbitrary monomer in the bulk of the polymer (not at the edges). When this monomer is perturbed by random thermal fluctuations, it slightly “drags” its neighbors (springs pull!). Those neighbors, in turn, drag their neighbors and so on. Over time, the initial monomer effectively gets coupled to a growing chunk of the chain around it. Let the size of this chain, which is a function of time, be 𝑁(𝑡). In the Rouse model, the dynamics of each monomer are overdamped because the surrounding solvent exerts a strong viscous drag that dominates over inertial effects (the term 𝜁 in the original equation). At the microscopic scale, a monomer is tiny and does not accumulate momentum – any significant displacement is immediately opposed by a frictional force proportional to its velocity, leading to slow relaxational motion rather than linear wave-like propagation. As a consequence, a local change in the position of one monomer is transmitted to its neighbors only gradually, in a way that is mathematically analogous to a diffusive process along the monomers (just like the diffusion equation, the Rouse differential equation contains a first order time derivative and a second order space derivative – the spring forces term! – along the monomer index). The influence of a disturbance spreads outward with the characteristic diffusion law 𝑁(𝑡)∝√𝑡. Thus, after time 𝑡, only a segment of roughly 𝑁(𝑡) monomers around the tagged bead has responded to the perturbation. Now, as we have seen using random walk statistics: Δ𝑟2 ~ 𝑎2𝑁. After time 𝑡, only the 𝑁(𝑡)- monomer chunk has had time to respond to the fluctuations of the initial monomer. So the mean-square displacement of the initial monomer after time 𝑡 is of the same order as the equilibrium variance of that chunk: ⟨Δ𝑟2⟩=𝑎2𝑁 ~ 𝑎2√𝑡∝𝑡1/2 So the typical displacement scales as: Δ𝑟(𝑡) ~ 𝑡1/4 This is the distance a monomer moves due to internal chain motions before the chain starts drifting collectively. We have already seen that, for a random walk polymer, ⟨Δ𝑟2〉=𝑎2𝑁. To estimate how the size of this random walk polymer scales, we look at the square-root of ⟨Δ𝑟2〉: √⟨Δ𝑟2〉 ~ √𝑁 ~ 𝑁1/2 Thus, the size of an ideal polymer scales as Δ𝑟 ~ 𝑁1/2 under random walk statistics. This means the monomer loses “memory” of its initial position only when it has diffused roughly the full chain size: Δ𝑟(𝑡) ~ 𝑁1/2 Substitute this result in Δ𝑟(𝑡) ~ 𝑡1/4, we get: Δ𝑟(𝑡) ~ 𝑡1/4 ~ 𝑁1/2 It is important to note that Δ𝑟(𝑡) ~ 𝑡1/4 holds in the subdiffusive regime, and since 𝑁1/2 is a measure of the full size of the chain, Δ𝑟(𝑡) ~ 𝑁1/2 indicates that the monomer has already diffused through the entire chain, and hence fully lost the “memory” of its initial position. Hence, the order of time we get by comparing 𝑡1/4 and 𝑁1/2 would give us a measure of the timescale of the crossover of the polymer from subdiffusive to diffusive regime: 𝑡1/4 ~ 𝑁1/2 → 𝑡 ≡𝜏 ~ 𝑁2 At timescales 𝜏 ~ 𝑁2, the monomer’s local subdiffusive motion has carried it a distance comparable to the overall chain size, and the polymer diffuses freely since all “memory” of the initial configuration is lost. Essentially with the passage of time, the more the polymer “relaxes” (the farther the polymer moves from its initial configuration), the more freely it can diffuse – this is one of the defining features of the Rouse model. At first glance, it may seem counterintuitive that shorter polymers (small 𝑁) diffuse more rapidly (since 𝜏 ~ 𝑁2). Monomers in a shorter chain are more strongly correlated, so we would expect their motions tightly constrained, and the polymer should move more slowly and take more time to diffuse. However in the Rouse picture, these correlations actually make the entire chain move together more coherently and diffuse faster. The total friction experienced by the polymer scales with the number of monomers 𝑁, so the diffusion coefficient of the center of mass of the polymer is of the order: 𝐷𝑐𝑚 ∼1 𝑁 A shorter polymer, having fewer monomers, experiences less total drag and thus diffuses more easily as a whole. While longer chains exhibit slower, less coordinated motion due to internal fluctuations, shorter chains behave almost like a single rigid body, leading to faster overall diffusion despite their stronger internal correlations. We can intuitively try to understand why 𝐷𝑐𝑚 scales as 1/𝑁. For a polymer consisting of 𝑁 monomers, the center of mass 𝑟𝑐𝑚 is given by: 𝑟𝑐𝑚 =1 𝑁∑𝑟𝑖 𝑁 𝑖=1 Each monomer feels some friction (friction coefficient 𝜁) from the surrounding medium and random thermal noise 𝜉𝑖(𝑡). The Rouse equation for one bead is: 𝜁𝑑𝑟𝑖 𝑑𝑡 =(spring forces) +𝜉𝑖(𝑡) Now, we consider all 𝑁 monomers and sum over all the equations. All the internal spring forces cancel – every spring acts equal and opposite on its neighbors. This leaves us with: 𝜁∑𝑑𝑟𝑖 𝑑𝑡 𝑖=∑𝜉𝑖(𝑡) 𝑖 Using 𝑟𝑐𝑚 =1 𝑁∑𝑟𝑖 𝑁 𝑖=1 , we may rewrite the left hand side of the above expression: 𝑁𝜁𝑑𝑟𝑐𝑚 𝑑𝑡 =∑𝜉𝑖(𝑡) 𝑖 Now, each random thermal force 𝜉𝑖(𝑡) is independent, with variance 𝜎2 (say). When we sum 𝑁 of them, the variances add up and the variance of the total noise is the individual variance multiplied by 𝑁: ⟨(Σi𝜉𝑖−⟨Σi𝜉𝑖〉)2⟩∝𝑁𝜎2 Since for thermal noise ⟨Σi𝜉𝑖〉=0, we get: ⟨(Σi𝜉𝑖)2⟩∝𝑁𝜎2 This means: Σi𝜉𝑖 ~ √𝑁𝜎 Using this result in 𝑁𝜁𝑑𝑟𝑐𝑚 𝑑𝑡 =∑ 𝜉𝑖(𝑡) 𝑖, we get: 𝑁𝜁𝑑𝑟𝑐𝑚 𝑑𝑡 ~ √𝑁𝜎 From this, we see that the typical velocity of the center of mass is: 𝑣𝑐𝑚 =𝑑𝑟𝑐𝑚 𝑑𝑡 ~√𝑁𝜎 𝑁𝜁 In the above expression, √𝑁𝜎 is a measure of the typical amplitude of noise that the polymer is subjected to – the noise pushes the particle away from its initial position and contributes to its diffusion. The denominator 𝑁𝜁 represents the friction, which is obviously inversely related to the velocity. We next define a noise correlation time 𝜏𝑐, which is simply the maximum time for which the noise remains correlated (for times greater than 𝜏𝑐, the system can be treated to be independent of earlier random kicks). In other words, the polymer takes an uncorrelated step after time intervals of 𝜏𝑐. Over the correlation time 𝜏𝑐, the displacement of the center of mass would scale as: Δ𝑟𝑐𝑚 ∼𝑣𝑐𝑚𝜏𝑐 And thus, the mean squared displacement per kick scales as: ⟨(Δ𝑟𝑐𝑚)2〉 ~ 𝑣𝑐𝑚 2𝜏𝑐 2 ➢ If it is an end monomer, allow it to “orbit” its neighbor by one lattice step, while respecting self-avoidance (end flips) ➢ If it is not an end monomer, create a list of the sites that keep both neighbor bonds at length 1 ➢ Return those candidate positions for the Monte Carlo step to choose from (equiprobably) 5. Once the move has been decided, update the configuration and repeat (and generate visuals of the polymer chain on the 2D lattice at regular intervals, as specified by the plot interval). First, we start with a straight chain (fully stretched initial configuration) – as a test case. We have chosen lattice size 𝐿=70 (so we have a 70X70 lattice, and 4900 lattice sites in total), 𝑁=50 (the number of monomers, or length of the polymer), number of steps = 5000, plot interval = 100. We look at some snapshots at different Monte Carlo steps to visualize the behavior of the polymer. Figure 2: Monte Carlo simulation of a polymer chain of length 50 under discrete Rouse dynamics on 70X70 2D lattice for straight initial configuration Starting from a fully stretched (straight) initial configuration, we see that the polymer undergoes visible relaxation over successive Monte Carlo steps, gradually adopting more compact and irregular conformations. Even though we do not have any explicit attractive potential here, the chain appears to “shrink” slightly as it explores configurations that satisfy the nearest-neighbor distance constraint while avoiding self-intersections. This reflects the entropic tendency of the polymer to sample a wide range of conformations rather than remain extended. It should also be recalled that the simulation is inherently stochastic – at each step, each monomer move is selected randomly from the set of allowed options – and hence we would get slightly different conformations if we rerun the same simulation without changing any of the parameter values. However, we would still observe a similar overall relaxation behavior. We now slightly change the parameters: 𝐿=50, 𝑁=30, and start with a fully stretched configuration like before: Figure 3: Monte Carlo simulation of a polymer chain of length 30 under discrete Rouse dynamics on 50X50 2D lattice for straight initial configuration When a shorter polymer chain (𝑁=30) is simulated under the same Monte Carlo conditions, we see that the overall relaxation occurs more rapidly and the chain displays more pronounced bending within the same number of steps. This is not surprising, because shorter chains have fewer internal degrees of freedom and hence, stronger correlations with the neighboring monomers – any local move affects a significant portion of the entire polymer, effectively bending it quickly and appreciably even if the initial configuration was fully straight. The lattice size (𝐿=50) plays no significant role here, since it remains much larger than the polymer itself (we have chosen 𝐿=50 here instead of 𝐿=70 for visual clarity). We now visualize the same polymer as above (𝑁=30, 𝐿=50), but after removing the constraint of volume exclusion. We clearly see that the Gaussian polymer overlaps with and intersects itself, and for the same number of Monte Carlo steps the Gaussian polymer collapses to a much greater extent than the self-avoiding polymer. Figure 4: Monte Carlo simulation of a Gaussian polymer chain of length 30 under discrete Rouse dynamics on 50X50 2D lattice for straight initial configuration Let us now run the simulation starting from a random initial configuration. This time, we choose an intermediate polymer length 𝑁=40, and since we are starting from a random (bent) configuration, we choose a slightly smaller lattice for better visual clarity, 𝐿=30. We keep the number of steps = 5000 and plot interval = 100. When generating a random initial configuration, we must ensure to avoid situations in which the generating process becomes trapped in a dead end. In the procedure used here, the chain is initiated near the center of the lattice and extended step-by-step in both directions by choosing an available nearest-neighbor site (according to the rule of preserving nearest-neighbor distances). Although this method usually succeeds for moderate chain lengths, it is still possible for the growth to terminate prematurely (for example, we may specify 𝑁=100, but the process terminates at 𝑁=93). This is especially important in two dimensions – due to the limited availability of free neighboring sites in 2D, self-avoidance strongly restricts the space of accessible configurations. Growing the chain outward in both directions reduces the likelihood of getting trapped, but it does not eliminate the risk entirely, and it is always a good idea to check that the algorithm successfully generates a full-length configuration. Figure 5: Monte Carlo simulation of a polymer chain of length 40 under discrete Rouse dynamics on 30X30 2D lattice for random initial configuration We see that for a random initial configuration, the polymer undergoes a more equilibrated evolution, showing no strong tendency to either swell or collapse. The polymer keeps exploring a wide range of conformations, fluctuating around an average size determined by the balance between connectivity constraints and self-avoidance. The absence of any energetic bias means the dynamics are purely entropic – the chain rearranges locally while preserving bond lengths and avoiding overlaps. We now repeat the above simulation for a Gaussian polymer with 𝑁= 40, 𝐿=30 as above. This time, we only show three snapshots instead of six. Figure 6: Monte Carlo simulation of a Gaussian polymer chain (N=40, L=30) under discrete Rouse dynamics for random initial configuration Next, we plot the square of the end-to-end distance of the polymer over time (or number of Monte Carlo sweeps). The end-to-end distance is simply the distance between the first and last monomer on the chain: 𝑅𝑒=|𝑟𝑁−𝑟1| We also plot the cumulative mean of 𝑅𝑒 2, which means at each time step we also plot the mean of all prior values of 𝑅𝑒 2 up to that step. It is more natural to analyze ⟨𝑅𝑒 2⟩ rather than ⟨𝑅𝑒⟩ because the squared distance averages in a statistically clean manner under thermal fluctuations and adds linearly across independent segments, making it directly comparable to analytical scaling laws. For the following plot, we have taken a very large number of Monte Carlo steps (1000000); this ensures the polymer has sufficient time to explore its configuration space, and that in turn ensures smoother plots with more statistical accuracy. Figure 7: Time series of squared end-to-end distance and its mean, for a polymer chain of length 40 under discrete Rouse dynamics, 1000000 Monte Carlo steps, 30X30 2D lattice, random initial configuration The curve keeps fluctuating. This is the behavior we observed when we were visually examining the different polymer configurations on the lattice over time for a random initial configuration – the polymer did not definitively lean toward either collapsing or swelling, but kept fluctuating and exploring different conformations. In the above plot, we have plotted the square of the distance between the first and last monomers as the polymer evolves through 1000000 Monte Carlo steps, and 𝑅𝑒 2 keeps fluctuating as the polymer explores different conformations. In short, the polymer keeps wiggling forever. We would expect it to spend equal time stretched and compressed – something which may not be visible directly on the plot due to a number of reasons. First, even though 1000000 steps might seem like a large enough number, it may not always be sufficient to exclude effects like the size of the box (lattice), biases in the initial configuration and slow mixing. To illustrate this, we run the same simulation as above once more, and see that the cumulative mean converges to a different value this time. Figure 8: Rerun of the above simulation, same parameters, same number of steps, random initial configuration In the above plot, we visually see that ⟨𝑅𝑒 2⟩ converges to a value slightly above 250. For the above simulation, we have considered a self-avoiding polymer with 𝑁=40 in two dimensions, (𝜈=3/4), hence we would expect ⟨𝑅𝑒 2⟩ ~ 𝑁2𝜈 ~ 𝑁3/2 ~ 253, which is indeed close to the converging value that we visually see on the plot. In a different run of the same simulation (figure 7), however, ⟨𝑅𝑒 2⟩ converges to a very different (and much smaller) value. This variation arises because each simulation begins from a single, randomly generated initial configuration and thus explores only one stochastic trajectory through the polymer’s configuration space. Although long Monte Carlo runs allow the system to sample many conformations, the sampling remains limited to the subset accessible from that particular initial configuration point within a finite simulation time. In principle, true ensemble convergence would require averaging over many independently initialized polymers or running extremely long trajectories. Obviously, in order to obtain reliable estimates of ⟨𝑅𝑒 2⟩ and its cumulative mean, it is essential that the observation time of the simulation significantly exceeds the polymer’s characteristic correlation time, 𝜏𝑐∼𝑁2𝜈+1. Only beyond this timescale do successive configurations become effectively uncorrelated, allowing the running mean to converge toward the true ensemble average. This temporal decorrelation requirement is particularly strong for long polymers, since 𝜏𝑐 grows rapidly with 𝑁. However, even after ensuring long simulation times and averaging over multiple independent realizations, we see that the mean value of ⟨𝑅𝑒 2⟩ does not necessarily equal the simple scaling 𝑁2𝜈 (figures 9, 12). This is because the correct relation is ⟨𝑅𝑒 2⟩=𝐾𝑎2𝑁2𝜈, where 𝐾 is a proportionality constant that depends on microscopic details of the model. The factor 𝐾 is not necessarily 1; in lattice models it can be substantially smaller, leading to numerically lower ⟨𝑅𝑒 2⟩ values than the ideal scaling estimate. When ⟨𝑅𝑒 2⟩ is plotted against 𝑁 on a log-log scale, we should generally obtain a straight line, and the slope provides an empirical measure of 2𝜈. In figure 11, we obtain a slope of about 1.2 instead of the expected 1.5 (corresponding to 𝜈=3/4 in 2D). This deviation can possibly be attributed to incomplete relaxation, limited ensemble sampling, etc. In figure 14, we create the same plot for six polymers (with 𝑁=20,30,40,50,60,70), and we get a slightly better fitted slope of about 1.3. Figure 9: Rerun of the above simulation for longer time, after averaging over 10 independent runs, same parameters, same number of steps, random initial configuration Figure 10: Time series of the cumulative means of 𝑅𝑒 2 for 10 independent runs, and their ensembleaveraged curve Figure 11: log ⟨𝑅𝑒 2⟩-vs-log𝑁 plot of the above simulation Figure 12: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ averaged over 100 independent runs Figure 13: 100 independent runs in groups of 10, with error bands around ensemble-averaged curve Figure 14: log ⟨𝑅𝑒 2⟩-vs-log𝑁 plot for 6 polymers It would also be a good idea to look at the log ⟨𝑅𝑒 2⟩-vs-log𝑁 plot for the above six polymers without self-avoidance; in other words, we assume all the above polymers to be Gaussian this time. We could expect a slope of exactly one, since for Gaussian polymers ⟨𝑅𝑒 2⟩ ~ 𝑁, and this time we get an almost perfect match (figure 15). Figure 15: log ⟨𝑅𝑒 2⟩-vs-log𝑁 plot for 6 Gaussian polymers We now look at the probability distribution 𝑃(𝑅𝑒) of the end-to-end distance for the original (non-Gaussian) polymers (figure 16). For moderate chain lengths, the distribution is slightly skewed, with a non-Gaussian tail rather than a perfect symmetric bell shape. This asymmetry arises because the polymer spends most of its time in compact, coiled configurations, while occasionally exploring rare stretched conformations that extend the end-to-end distance far beyond the mean. Also, here we have averaged over only eight independent runs, hence significant noise remains. Figure 16: Probability distribution of 𝑅𝑒 for three polymers of lengths 20, 40 and 80 We now look at plots of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ as functions of time for a smaller polymer (𝑁=20) in a larger lattice (𝐿=50); this would eliminate effects of the lattice boundaries on the motion of the polymer, and the polymer would decorrelate faster. This time, we plot for 500000 steps for proper comparison (1000000 steps for 𝑁=40 would be equivalent to 500000 steps for 𝑁= 20). Figure 17: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ for a polymer with N=20, L=50, 500000 Monte Carlo steps, random initial configuration Figure 18: Rerun of the above simulation, same parameters, same number of steps, random initial configuration We clearly see that now the differences between the converging value of ⟨𝑅𝑒 2⟩ over reruns of the code are smaller. And ⟨𝑅𝑒 2⟩ does not converge to a value near 𝑁, but always greater than 𝑁, since we have included the effect of volume exclusion in the above simulations. We now plot 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ as functions of time for a Gaussian polymer. average over multiple independent samples. In figure 26, we have plotted the ensembleaveraged MSD curve of the above polymer over 10 independent runs. We see a cleaner crossover from subdiffusive to diffusive regimes at 𝑡 ~ 102. For a larger polymer (𝑁=40), we see a slower crossover from subdiffusive to diffusive regime (figure 27), as expected from 𝐷𝑐𝑚 ~ 1/𝑁. In figure 27, the crossover from subdiffusive to diffusive regime takes place at 𝑡>103. This makes sense, because for this plot 𝑁=40, which means 𝑁2.5 =10119.29, which is of the order of 104 (this is a self-avoiding polymer). We see a much denser saturation here because this time the polymer is significantly long as compared to the lattice size (𝑁= 40, 𝐿=30). Figure 27: MSD-vs-time curve for a longer polymer (N=40, L=30), 1000000 steps, random initial configuration Figure 28: MSD-vs-time curve for a longer polymer (N=40, L=80) chain, 10000000 steps, random initial configuration In figure 28, we simulate a self-avoiding polymer of the same length, but in a larger lattice and for a much higher number of steps (𝑁=40, 𝐿=80, 10000000 steps). In this case too, we see similar behavior (the order of 𝜏𝑐 seems to be around 103), although there is significant noise and saturation. The saturation persists even in a bigger lattice and over longer times because our Monte Carlo scheme is local, and our moves relax the polymer’s internal modes without efficiently displacing its center of mass. Now, to compare the behavior of polymers with different lengths, we rescale the time axis by the characteristic relaxation time of the chain, and plot MSD as a function of 𝑡/𝑁2 (figure 29). In the subdiffusive regime – where the monomer motion follows MSD ∼𝑡1/2 – overlaps well for chains of different lengths (considering this curve has not been averaged over different independent runs). Beyond the crossover however, the long-time diffusive regime (MSD ∼𝑡) does not align perfectly under this scaling. This mismatch arises because after decorrelation, the overall center-of-mass diffusion scales as 𝐷𝑐𝑚 ∼1/𝑁, meaning the scaled time should instead be 𝑡/𝑁, not 𝑡/𝑁2. In other words, there are two distinct dynamic scalings – 𝑡/𝑁2 for subdiffusive internal relaxation and 𝑡/𝑁 for center-of-mass diffusion – and a single plot cannot simultaneously rescale both regimes into perfect overlap. Figure 29: MSD-vs-scaled time (t/N2) curve for Gaussian polymers of different lengths (N=10, 20, 30, L=60), 500000 steps, random initial configuration The important point is that plotting the MSD of different chain lengths against this scaled time 𝑡/𝑁2 removes the trivial shift in the crossover position that would otherwise occur when using unscaled time 𝑡, where longer chains appear to relax much more slowly simply because their intrinsic timescale is larger. Overall, this demonstrates that the underlying relaxation mechanism is the same across different chain lengths, and that the apparent differences in unscaled time arise solely from the 𝑁2 scaling of the global Rouse relaxation time. In the plot, we also see that the crossover seems to be not exactly at 𝑡/𝑁2∼1 (which we would ideally expect), but at a slightly larger value. This shift is not surprising, because we have implemented a local Monte Carlo scheme – only local, distance-preserving moves are allowed – and thus many attempted moves result in no displacement at all. These microscopic details effectively multiply the physical relaxation time by a constant 𝐾 (say), and the slowest mode really relaxes at 𝜏𝑐 ~ 𝐾𝑁2, leading to a crossover at t/N2 ~ 𝐾 rather than exactly 1. But clearly, after rescaling by 𝑁2, the three curves in the above plot collapse onto the same trajectory and exhibit the same subdiffusive 𝑡1/2 behavior and eventual diffusive drift. Thus, the underlying Rouse scaling is captured correctly. The noticeable noise in the high-𝑡 (diffusive) region arises from large fluctuations and limited statistics at long times (independent configurations are fewer and MSD accumulates stochastic variance). To compensate this, we next plot the ensembleaveraged MSD (averaged over multiple independent runs with different random initial configurations), for polymers of different lengths (figures 30, 31, 32). Figure 30: Ensemble-averaged MSD-vs-scaled time (t/N2) curve over 10 independent runs, for Gaussian polymers of lengths N=10, 20, 30, L=60 Figure 31: Ensemble-averaged MSD-vs-scaled time (t/N2) curve over 10 independent runs, for Gaussian polymers of lengths N=20, 40, 80, L=100 Figure 32: Ensemble-averaged MSD-vs-scaled time (t/N2) curve over 100 independent runs, for Gaussian polymers of lengths N=20, 30, 40, 50, 60, 70, L=100 In the above plots, the random fluctuations in the diffusive region are significantly suppressed, and we get a smoother and more statistically representative depiction of the MSD across both regimes. Finally, we briefly examine the dynamics of the center of mass of a Gaussian polymer. By tracking the instantaneous position of the center of mass, 𝑟𝑐𝑚 =1 𝑁∑𝑟𝑖 𝑁 𝑖=1 , we can visualize how the entire polymer drifts as a single entity over time. The plotted trajectory of the center of mass (figure 33) for a polymer (𝑁=20) reveals a smooth, wandering path, reminiscent of a two-dimensional random walk, confirming that although individual monomers undergo complex correlated motions, their collective motion behaves diffusively on long timescales. Figure 33: Trajectory of the center of mass of a Gaussian polymer with N=20, under discrete Rouse dynamics We now plot the time series of the MSD of the center of mass, ⟨|𝑟𝑐𝑚(𝑡)−𝑟𝑐𝑚(0)|2⟩, for polymers of three different lengths (𝑁=20,40,80), on the log-log scale (figure 34). All three curves show an almost perfect linear dependence on time, reflecting the expected diffusive scaling MSDcm ∼𝑡. The longer polymers lie systematically below the shorter ones, indicating a slower diffusion rate (recall 𝐷cm ∼1/𝑁). The center of mass of a longer polymer diffuses more sluggishly because its overall friction increases proportionally to the number of monomers, while the thermal driving force remains constant per monomer. Figure 34: Ensemble-averaged MSD-vs-time curves over 10 independent runs, for Gaussian polymers of lengths N=20, 40, 80, L=120 In the MSD plot, very short times were deliberately excluded. This is because at early times the polymer’s internal relaxation dominates, and the center of mass barely moves. Removing the very early time points helps isolate the regime of steady diffusive motion, where the center of mass dynamics have decorrelated from the initial configuration and the expected 𝑡1 scaling becomes clear and smooth. Figure 35: Ensemble-averaged MSD-vs-time curves over 10 independent runs, for self-avoiding polymers of lengths N=20, 40, 80, L=120 For a self-avoiding polymer, the center of mass motion exhibits slower diffusion (figure 35). The shortest polymer (𝑁=20) seems to grossly follow the expected 𝑀𝑆𝐷∼𝑡 scaling; however there is noticeable flattening for longer polymers, owing mainly to excluded volume effects – monomers can now no longer freely pass through one another, and the longer chains are significantly impacted because of this; they feel an increased effective friction which causes them to diffuse more slowly, consistent with the scaling 𝐷cm ∼1/𝑁. In conclusion, the transition from a 𝑡1/2 scaling to a 𝑡 scaling in the mean squared displacement marks the crossover from subdiffusive to diffusive behavior in the polymer’s dynamics. In the short-time regime (subdiffusive), each monomer is trapped by its neighbors along the chain, and its motion is hindered. At longer times, once these internal modes have relaxed sufficiently, the entire polymer moves collectively and the motion becomes diffusive, dominated by the free Brownian diffusion of the polymer’s center of mass. Polymer Dynamics under Heterogeneous Geometric Rules In this section, we explore the effect of simple, but spatially heterogeneous modifications to the geometric rules on the behavior of the polymer. Unlike the standard Rouse dynamics we have so far explored – where every monomer obeys identical distance-preserving moves and translational symmetry is preserved along the chain – we now assign distinct local rules to alternating monomers. Apart from alternating copolymers, we also build a toy model of a block copolymer, where we divide the polymer into two continuous halves and apply different geometric rules to each half. The objective is not to replicate any specific microscopic mechanism, but to explore how nonuniform mobility, imposed purely through geometry and bond-preserving lattice moves, can qualitatively model some aspects of the nonequilibrium polymer dynamics. Based on what we have explored so far, we aim to build an extremely lightweight analogue of situations in which local polymer environments differ in temperature, effective viscosity, crowding or other mechanical constraints. In biological systems, polymers such as DNA, chromatin fibers and RNA regularly encounter heterogeneous and spatially varying surroundings: some regions encounter strong steric confinement, others transient anchoring and still others enhanced agitation. While such systems are far more complex than anything we model here, our constructions provide a simple way to introduce local variability without invoking explicit forces, potentials or energetic parameters. All deviations from equilibrium arise solely from asymmetries in the allowed geometric moves for different monomers. It is important to note that breaking the uniformity of the geometric rules breaks the effective translational symmetry along the chain – even if the polymer starts in a fully symmetric conformation, the differing local mobilities disrupt detailed balance at the algorithmic level, often pushing the system into a driven steady state. In the original Rouse model, every monomer follows the same geometric rules, has the same mobility and attempts the same kinds of moves equiprobably. Thus, any microscopic change in configuration has an equally probable reverse change. This one-to-one symmetry between the probability of forward and backward transitions ensures detailed balance. Even though the chain undergoes diffusion, there is no built-in preference for moving in any particular direction through the configuration space. However in our alternating-rule models, we will deliberately assign different kinds of moves to different monomers, and because of this, certain configuration changes become more probable than their exact reverses. This breaks detailed balance, even though all bonds remain intact and no external forces are applied. The system no longer relaxes to the same equilibrium statistics as the uniform Rouse model, but instead settles into a nonequilibrium steady state driven solely by differences in the permitted geometric motions along the chain. To introduce a simple form of spatial heterogeneity into the lattice polymer dynamics, we first try constructing a toy nonequilibrium model in which different monomers obey different geometric update rules. The polymer is initialized as a self-avoiding walk on a 40X40 lattice (figure 36), with all bond lengths restricted to one lattice unit. We permanently fix a small fraction of monomers (~ 3%) to their initial lattice site to mimic local environmental pinning. All remaining monomers must preserve the bond length to their neighbors at every step. The alternating rules are as follows: 1. Monomers with even-number indices follow the familiar rule – they can only attempt moves that preserve their distances from their nearest neighbors. 2. Odd-indexed monomers are allowed to attempt moves to any of the eight surrounding lattice sites (including diagonal moves) and are additionally allowed to attempt two-step moves chosen at random, within the same 8-neighbor geometry. The inclusion of two-step moves in the rule set was intended as a purely geometric way to encode heterogeneous local mobility, loosely analogous to segments of a chain experiencing regions of differing local viscosity, crowding or agitation. However, despite this deliberately heterogeneous mechanism, the resulting polymer dynamics is very constrained – the polymer hardly changes its configuration significantly even after a large number of steps. The key reason is that the larger move set available to the odd monomers is almost entirely suppressed by the bond-preservation condition – which is a very tight constraint. Any candidate move must maintain a Manhattan distance of exactly one to both neighboring monomers, but diagonal moves typically increase this distance to two or more, and two-step moves almost always violate the bond-length constraint. Once volume exclusion and lattice boundaries are also imposed, only a very small subset of the nominally available odd-monomer moves remain viable. In practice, even though odd monomers attempt more varied moves, almost all of these proposals are rejected. Figure 36: Monte Carlo simulation of a self-avoiding polymer chain of length 100 under the rules of our alternating-rule toy model 1 on 40X40 2D lattice for random initial configuration As a result, the alternating-rule design in our first toy model does not generate the intended nonequilibrium effects, even though we start from a random initial configuration. The chain is essentially “caged” by its self-avoiding structure and the strict bond-length requirement, so the heterogeneity in the proposal rules does not translate into appreciable heterogeneity in accepted moves on the lattice. This highlights an important lesson: under strong geometric constraints, blindly introducing heterogeneous rules and modifying only the proposal distribution is insufficient to induce substantial dynamical differences. If we remove the condition of self-avoidance from the above model, we see a very slight increase in the variations in the configuration of the polymer through the same number of steps (figure 37), but this effect is extremely weak and the polymer still remains significantly rigid. Figure 37: Monte Carlo simulation of a Gaussian polymer of length 100 under the rules of our alternating-rule toy model 1 on 40X40 2D lattice for random initial configuration Now, we try a different alternating-rule toy model. First, we start with a Gaussian polymer (no self-avoidance). Like before, even-indexed monomers still follow the familiar rule: after every move the nearest-neighbor distances must be preserved. But this time, for odd monomers we do not check that the nearest-neighbor distances remain constant post-move. This means, bonds to neighbors can stretch or compress arbitrarily for odd-indexed monomers! Further, we assign three possible rules to the odd-indexed monomers, which can be thought to represent regions of the polymer experiencing heterogeneous or “active’’ local environments, and at each step, one out of these three rules is equiprobably and randomly chosen and implemented on the oddindexed monomers. The three possible behaviors of the odd-indexed monomers are: 1. They remain completely fixed for that update attempt. 2. They move to one of the four nearest lattice sites, overlap allowed. 3. They move to any neighboring lattice site that lies one or two steps away (including diagonal steps). These modes intentionally differ in their local dynamical freedom, and the random switching breaks full translational and dynamical symmetry along the chain. As in all previous simulations, the polymer is confined by hard boundaries at the edges of the lattice. Under this second alternating-rule model, the most striking feature is the appearance of large, irregular and often physically unrealistic fluctuations in the shape of the polymer (figure 38). This behavior arises primarily because odd-indexed monomers are no longer required to maintain a fixed bond length with their neighbors, allowing the polymer to momentarily “stretch’’ or “collapse’’ unpredictably. Even though the even-indexed monomers still obey strict nearest-neighbor constraints – providing a partial scaffold – the odd-indexed monomers periodically break local structure, producing abrupt jumps, distorted segments and large deviations from standard Rouse dynamics. It is clear that under this model, the polymer keeps expanding overall. As a consequence, the squared gyration radius keeps steadily increasing (figure 39). As we have previously seen, in a normal Rouse or self-avoiding chain, 𝑅𝑔 2(𝑡) fluctuates around an equilibrium value, because the dynamics satisfy detailed balance. Here however, the fluctuating and unbounded bond lengths of odd monomers explicitly break detailed balance. The polymer is definitively “pushed’’ outward by stretches that have no corresponding restoring moves. As a result, 𝑅𝑔 2(𝑡) increases steadily rather than stabilizing, reflecting a genuine nonequilibrium expansion driven by these asymmetric geometric rules. Figure 38: Monte Carlo simulation of a Gaussian polymer of length 40 under the rules of our alternating-rule toy model 2 on 80X80 2D lattice for random initial configuration Figure 39: Time series of gyration radius squared for a Gaussian polymer of length 40 under the rules of our alternating-rule toy model 2 on 80X80 2D lattice for random initial configuration One may ask how fluctuating bond lengths are even possible in this model, given that every odd monomer is directly connected to two even monomers whose bond lengths are strictly In figure 46, we examine the MSD of both the driven and constrained blocks of the polymer separately. We plot the running mean of all the MSDs of all the monomers within each half (orange: mean over first 25; blue: mean over last 25), taken relative to the run’s initial positions. Averaging across the whole half reduces single-site noise and highlights the systematic mobility difference between the blocks – the orange mean MSD rises more steeply because the driven monomers have larger accessible move sets and higher attempt frequency as compared to the constrained (blue) monomers. Plotting the MSD of a single “representative” monomer from each block can give ambiguous or misleading signals (for instance, if the chosen monomer is very close to the ends). Figure 46: MSD-vs-time plots for both driven and constrained blocks of the block copolymer, averaged over all monomers in each half The shaded translucent regions around the curves represent the standard error of the mean across those six independent runs. We see that the driven half genuinely explores more configuration space (so its averaged MSD is larger), but global geometry and exclusion prevent unbounded swelling (𝑅𝑔 2 therefore fluctuates around a finite range). This model is a minimal, geometric toy that captures the intuition of a localized mobility/temperature gradient: stronger local activity (the driven block) increases local fluctuations and diffusion, yet connectivity and excluded-volume couple those fluctuations to the rest of the chain and keep the whole polymer in a constrained, quasi-steady state. It is clear that the MSD curves for the orange and blue monomers clearly diverge: the driven block exhibits a much steeper growth in MSD (as compared to the constrained block), reflecting its higher local mobility and larger accessible configurational volume. It is interesting to focus on what happens at the interface – the bond between the last orange and first blue monomers. The last orange monomer is free to attempt all the long-range “driven” moves like the other orange monomers, but each move must still maintain a valid bond length with its immediate blue neighbor (valid according to the rules followed by the blue monomers), so the blue block enforces stricter geometric constraints on the last orange monomer that limit which of those moves actually get accepted. In our Monte Carlo scheme, bond validity is checked only locally; each time a move is attempted on a particular monomer, the constraints on the monomer are checked only against its two immediate neighbors. Thus restrictions at the interface propagate only one bond deep, and do not globally suppress mobility across the entire orange block. The rest of the orange segment remains highly mobile because once a move satisfies local geometry, it incurs no penalty from distant constraints. In figure 47, the black curve tracks the MSD of the last orange monomer (the interface orange monomer directly bonded to the first blue monomer). As expected, its diffusion lies in between the fully driven orange block and the constrained blue block. Figure 47: MSD-vs-time plots for the driven block, the constrained blocks (averaged over all monomers in each half) and the last driven monomer at the interface of the block copolymer, averaged over 100 independent runs Figure 48: MSD of the center of mass of the block copolymer One might expect that over long times the entire polymer diffuses as a single object, eventually causing the MSD of all monomers (orange as well as blue) to converge. However in this model, the differing dynamical rules are applied continuously, not just at initialization, and the driven and constrained blocks remain kinetically distinct at all times. This is why the blue MSD curve never merges with the orange curve, even after long times, and our copolymer cannot be trated as a single Brownian particle (it is a composite object with a sustained mobility gradient). While the center of mass of the entire polymer undergoes diffusion governed by all monomers collectively (see figure 48), the internal relative motion continues to reflect the imposed heterogeneous dynamics. Thus, the blue block remains systematically slower, and its MSD curve does not asymptotically merge with the orange MSD. This dynamical asymmetry is not something that “averages out”, because it is actively enforced at every update step. In figure 49 (the higher noise is because here we have not averaged over 100 independent runs), we remove the asymmetry in the attempt-rate; orange monomers are no longer chosen to update twice as often as the blue ones, now all monomers are selected for updates with equal probability. From the position of the black curve, we see that now the last orange monomer at the interface almost does not exhibit the intermediate diffusivity observed earlier. Although it still has the larger geometric move set characteristic of the driven block, most of those moves are rejected because it is now not chosen more frequently than the blue monomers, and as before, the bond to its blue neighbor must remain within the constrained set of allowed distances. Without being updated more frequently, the interface bead becomes limited by the mobility of the slower block it is tethered to, and its extra freedom is never fully expressed in the dynamics. Consequently, its MSD lies much closer to the blue curve rather than interpolating between orange and blue as in figure 47. Figure 49: MSD-vs-time plots for the driven block, the constrained blocks (averaged over all monomers in each half) and the last driven monomer at the interface of the block copolymer, both orange and blue monomers now equally likely to be chosen at each step In a physical or biological context, such a heterogeneous block copolymer could correspond to a chain embedded in an inhomogeneous environment, for example one experiencing a gradual temperature gradient or a spatial variation in local crowding or solvent density from its one end to the other. The driven block would then represent a “hotter” and less crowded region where thermal agitation is stronger and monomer motion is enhanced, while the constrained block would correspond to a “colder,” more viscous or sterically hindered region. Similar situations occur in living cells, where biopolymers can span regions with significantly different microenvironments. Thus, this toy model – despite being purely geometric – captures an essential idea: how local heterogeneity in effective temperature or crowding can shape global polymer conformation and dynamics, producing nonequilibrium steady states that still respect strong geometric constraints. Figure 50: Snapshot from a different run, with slightly different parameters, of the block copolymer model discussed above In the end, it is important to note that all the models explored in this section are not intended to mimic any specific biological mechanism, but they provide good illustrations of how heterogeneous local mobilities can generate nonequilibrium behavior even in an otherwise structureless polymer. It would be interesting in future work to test additional combinations and variations of heterogenous geometric rules and constraints. It is also crucial to note that the resulting dynamics in our alternate-rule model do not guarantee increased fluctuations in general; depending on the specific choice of the rules, some portions of the chain may exhibit larger configuration-space excursions while other portions may become more constrained (as compared to standard Rouse-like dynamics). This flexibility could be exploited to probe a relatively wide spectrum of behavior – from unexpectedly stiffened motion to extremely quick relaxation – within a single unifying geometric framework. Overall, these toy models provide a controlled environment for examining how spatial heterogeneity alone, implemented through simple lattice-based rules, can produce departures from the equilibrium-like behavior. But of course, these are toy models and are not intended to provide quantitative predictions for real, nonequilibrium polymers. Discussion and Conclusion In this report, we have explored Rouse model analytically, and then simulated a simplified, discrete version of the Rouse model with self-avoidance in order to understand how global polymer behavior can emerge from simple local rules. By constraining each monomer to move in ways that preserve the nearest-neighbor distances and avoid overlap, we have observed how a polymer chain evolves on a 2D lattice through purely entropic dynamics. Despite the absence of explicit energetic terms or solvent interactions, the system displayed key qualitative features of real polymer motion – relaxation from stretched or random initial configurations, fluctuations in end-to-end distance and gyration radius, and a clear crossover from subdiffusive (𝑡1/2) to diffusive (𝑡) regimes over time. These results highlight how the interplay of connectivity and self-avoidance alone can encode much of the essential physics of polymer motion. Then we introduced heterogeneity into the update rules to probe the onset of nonequilibrium behavior in a controlled manner. The alternating-rule models and the block copolymer model demonstrated that modifying only the proposal mechanisms – while keeping all moves purely geometric – can break detailed balance and produce a wide range of interesting behavior. Some heterogeneous models produced minimal change in global behavior, because the fixed bondlength and self-avoidance constraints reject most asymmetric proposals. However, in models in which bond length preservation was relaxed for a subset of monomers, we saw strong nonequilibrium effects, including persistent expansion and fluctuations in the gyration radius. Intermediate models combining enhanced mobility, intermittent immobilization and geometric constraints of fixed bond length and self-avoidance produced nonequilibrium steady states that remained close to equilibrium, with the polymer “breathing’’ around a bounded size rather than drifting outward indefinitely. In real systems, polymer conformations arise from a competition between energy and entropy: attractive interactions between monomers in poor solvents (solvents that repel monomers) cause the chain to collapse, while thermal motion in good solvents (solvents that attract monomers) promotes swelling. Although the models we have explored in this report do not include explicit attraction or repulsion, the geometric constraint of self-avoidance encode a behavior similar to entropic resistance to collapse in the dynamics. Ultimately, these results illustrate that even a minimal rule set – purely local, geometric and sometimes heterogeneous – can give rise to collective, emergent dynamics at the macroscopic level. This reflects one of the central insights in the study of complex systems: global behavior often stems not from intrinsically complex rules or intricate microscopic laws, but from the repeated logic of simple, local constraints. References • Schiessel, H. (2014). Biophysics for Beginners: A Journey through the Cell Nucleus. Pan Stanford Publishing. • Padding, J. T. (2005). Theory of Polymer Dynamics. University of Cambridge. • Li, B., Madras, N., & Sokal, A. (1994). Critical Exponents, Hyperscaling and Universal Amplitude Ratios for Twoand Three-Dimensional Self-Avoiding Walks. arXiv:heplat/9409003. • Van Leeuwen, J., & Drzewiński, A. (2009). Stochastic Lattice Models for the Dynamics of Linear Polymers. Physics Reports, 475(5–6), 53–90. • Newman, M. E. J., & Barkema, G. T. (1999). Monte Carlo Methods in Statistical Physics. Oxford University Press. Appendix: End-Monomer Interaction as a Minimal Model for Polymer Collapse In this appendix, we explore a very simplified polymer model, where we explicitly introduce energetic interactions, but only between the two end monomers (all internal monomers feel no energetic bias at all). This minimal setup lets us probe how global conformational transitions get modified by even just boundary interactions. So every internal monomer does not experience any energetic bias, and in context of our previous discussions, still undergoes purely geometric Rouse-like moves. There is an (attractive or repulsive) interaction only at the end points (between the two ends of the chain). As the interaction strength 𝐽 varies, the polymer can neatly transition from an extended, swollen state into a collapsed configuration where the ends prefer to meet. In a sense, the polymer’s long-term behavior is largely determined by the end monomers, and the rest of the chain negotiates with whatever the ends demand. To formalize this, we write a Hamiltonian involving only the end monomers (we use a Kronecker delta contact term between the two end monomers indexed 1 and 𝑁): 𝐻=−𝐽𝛿𝑟1,𝑟𝑁 where 𝑟1 and 𝑟𝑁 denote the lattice positions of the first and last monomer. According to the definition of the Kronecker delta, 𝛿𝑟1,𝑟𝑁=1 if 𝑟1=𝑟𝑁 (the two end monomers occupy the same lattice site), 𝛿𝑟1,𝑟𝑁=0 otherwise. This means when the end monomers interact, the Hamiltonian is 𝐻=−𝐽, and if 𝐽>0, we see that there is end-to-end interactions would be energetically favored, since they lower the overall energy of the chain by −𝐽. (Throughout this appendix, we assume 𝐽>0. If we want the end monomers to repel, we can penalize end-toend interaction by setting 𝐽<0.) Thus in a way, the entire energy of the polymer comes from a single binary event: either the ends meet (and the system gains an energy −𝐽), or they do not (and the Hamiltonian is zero). No internal monomer has any pairwise or self-interaction term. Despite its simplicity, this Hamiltonian represents a clear physical competition between entropy and energy. At 𝐽=0, the polymer behaves just like an ideal 2D random walk, with no reason for the ends to attract each other. In that case, the probability that the chain spontaneously forms a loop scales as: 𝑝0 ~ 1 𝜋𝑁 Intuitively, we can easily see why the above scaling makes sense for an ideal, Gaussian polymer in 2D. We have already looked at the analytics of a 2D random walk. Assume the walk starts at a fixed point (the fixed end). Then the position of the free end after 𝑁 steps can be thought of as the endpoint of the 2D random walk. As we have seen, after 𝑁 steps the typical end-toend distance scales as Δ𝑟 ~ √⟨Δ𝑟2〉∝√𝑁 (because ⟨Δ𝑟2〉=𝑎2𝑁, assume 𝑎=1 for simplicity). Now we can think of the endpoint as being distributed over a roughly circular region with a radius of the order of √𝑁. Of course, this circular-cloud picture is not completely accurate, and hence this is not a rigorous derivation. We are simply interested in an intuitive justification of why the probability of the end monomers interacting and forming a loop spontaneously (𝐽=0) scales as 1/𝜋𝑁 (for a polymer with 𝑁 monomers). So if the typical distance scale is √𝑁, then the characteristic region explored by the endpoint has an area on the order of 𝜋(√𝑁)2=𝜋𝑁 (assuming a circular region, of course). This area represents the “target zone’’ in which the endpoint is most likely to be found, and although of course the exact density is not necessarily uniform, for the purpose of our estimate, the endpoint may be regarded as being spread roughly evenly over this region. On the 2D lattice, a single site occupies an area of one unit cell, and hence the probability that the free end of the polymer lands precisely on the specific lattice site occupied by the first monomer (the other end) is approximately the ratio of the area of one lattice cell to the total area of the endpoint cloud: 𝑝0 ~ 1 𝜋𝑁 Clearly, the return probability decreases inversely with 𝑁: for longer polymers, the endpoint wanders over a correspondingly larger region, and the chance that it returns exactly to the origin becomes proportionally smaller. Even for 𝑁=40, 𝑝0 ~ 1/𝜋𝑁≈0.008. Hence, without explicit energetic interactions, it is very unlikely that the end monomers interact and the polymer forms a loop spontaneously. We now explore what happens when 𝐽≠0. Initially, we tried to directly model the collapse of the polymer using our local Monte Carlo scheme on the 2D lattice, in order to visualize how the normalized end-to-end squared distance, ⟨𝑅𝑒 2⟩/𝑁, varies with the interaction strength 𝐽. When 𝐽 is very small, we would not expect the polymer to collapse, and ⟨𝑅𝑒 2⟩/𝑁 would not vary significantly (since it is the normalized meansquared end-to-end distance). We would thus expect a plateau on the ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 plot for small 𝐽. As we keep increasing 𝐽, we would expect a relatively sharp drop toward zero (once the polymer decides to form a loop, the end-to-end distance becomes zero). Also, for longer polymers (larger 𝑁), this transition should occur at larger 𝐽 (because the entropic penalty of bringing two distant ends together increases with chain length). However, the Monte Carlo simulations did not reproduce this expected behavior (figure 51). Despite extensive sampling and averaging over independent runs, the plots of ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 for three polymers (of different lengths) were nearly flat on average. Figure 51: ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 plots using Monte Carlo simulations (averaged over 2 independent runs, error bars plotted) for three polymers of lengths 20, 40 and 80, with end-only attraction The reason is because the Hamiltonian 𝐻=−𝐽𝛿𝑟1,𝑟𝑁 only affects the two end monomers. For most polymer configurations, this contact event is extremely rare (𝑝0∼1/𝜋𝑁). Even for moderate chain lengths, 𝑝0 is already so small that the effective temperature range explored by even 𝐽≤500 is far too weak to overcome the entropic cost. Hence, the above simulation did not show a collapse transition. To see that directly, we would need extremely large 𝐽 values plus extremely long simulation times, which is computationally expensive and timeconsuming. To improve exploration of conformations, we next incorporated pivot moves, where a randomly chosen monomer acts as a hinge and the remaining chain segment is rotated or reflected by a lattice symmetry operation. Pivot moves do not alter bond lengths or violate chain connectivity, but they can dramatically change global structure in a single step, and increase the probability of the end monomers meeting. Despite this, and despite longer runs over multiple independent histories, and relatively short polymers (𝑁=20), the end-to-end contact remained effectively unobserved under our Monte Carlo scheme. Again, this is because the Hamiltonian contains only a Kronecker-delta interaction localized at the two ends, and the energetic gain occurs only when those ends precisely occupy the same lattice site. Since this represents a single microstate out of an exponentially large configuration space, unbiased Monte Carlo rarely encounters it. In order to capture the correct qualitative physics without such computational expense, we analytically simplified the system into a two-state model. The polymer is treated as either “noncontact,” where we assume the ideal result ⟨𝑅𝑒 2⟩≈𝑁; or “contact,” where the two ends overlap and 𝑅𝑒 2≈0. The relative statistical weight of these two states is governed by the Boltzmann factor associated with the contact energy −𝐽. If 𝑝0 is the contact probability at 𝐽=0, then the contact probability at arbitrary 𝐽 becomes: 𝑝𝐽=𝑝0𝑒𝛽𝐽 (1−𝑝0)+𝑝0𝑒𝛽𝐽 The ensemble-averaged end-to-end distance is then approximated as: ⟨𝑅𝑒 2⟩=𝑁(1−𝑝𝐽) This formula captures the correct physical trend. For low 𝐽 (when 𝐽→0), 𝑒𝛽𝐽 →1 and 𝑝𝐽→ 𝑝0. Since 𝑝0 is usually negligibly small (even for moderate length polymers), for low 𝐽 we effectively have 𝑝𝐽→0. Since the contact probability is vanishingly small at low 𝐽, from ⟨𝑅𝑒 2⟩=𝑁(1−𝑝𝐽) we have ⟨𝑅𝑒 2⟩→𝑁, and hence ⟨𝑅𝑒 2⟩/𝑁→1 (extended state, entropy dominates!). For large 𝐽 (when effectively 𝐽→∞), 𝑝0𝑒𝛽𝐽 →∞, and hence (1−𝑝0)+𝑝0𝑒𝛽𝐽 ≈ 𝑝0𝑒𝛽𝐽, which means 𝑝𝐽→1. Combining this with ⟨𝑅𝑒 2⟩=𝑁(1−𝑝𝐽) gives us ⟨𝑅𝑒 2⟩/𝑁→0 (collapsed state, energy wins!). We now plot ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 using the reduced two-state model, for the same three polymers (figure 52). We see that the plot of ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 cleanly displays the expected crossover from a swollen, coil-like configuration characterized by ⟨𝑅𝑒 2⟩/𝑁≈1, to a compact, looped state characterized by ⟨𝑅𝑒 2⟩/𝑁≈0. For relatively small values of 𝐽, the curve remains flat for a wide range, reflecting purely entropic behavior. Then the curve drops abruptly over a relatively narrow interval, indicating a strongly cooperative transition (which resembles a first-order phase transition). The above plot is actually a semilog plot, because we have plotted 𝐽 on a logarithmic axis for better visualization (the collapsed state extends over several orders of magnitude in 𝐽, and using a linear scale would compress the entire transition into an invisible sliver). Figure 52: ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 plots for three polymers (N=20,40,80) with end-only attraction, using the simplified two-state analytic model It is important to note that there is no true first-order phase transition in our model, even if the transition appears first-order-like. This is simply because for a true thermodynamic first-order transition, the observable should become non-analytic in some extreme limit (for example, 𝑁→∞). Here, however, the polymer is a finite chain with only a single energetic interaction between its two end monomers; no matter how sharply the end-to-end distance changes, the free energy remains smooth and analytic for all finite 𝑁. This steep sigmoidal profile also resembles the behavior of molecular switches often exploited in biological systems, for instance, a small change in binding energy can flip an entire polymer between distinct structural states. The parameter 𝐽 here acts as an effective control knob, tuning the competition between entropy and end-only attraction and thereby triggering a switch-like collapse. Now, we look at plots of the contact probability 𝑝𝐽 as functions of the interaction strength 𝐽, for the three polymers as before (figure 53). The 𝑝𝐽-vs-𝐽 plot shows the corresponding growth of the end-to-end contact probability as the interaction strength increases. For small 𝐽, 𝑝𝐽 remains close to zero because the entropic cost of bringing the two ends together dominates; once 𝐽 approaches the crossover point, the probability rises sharply in a sigmoidal, switch-like manner. This behavior is fully consistent with the ⟨𝑅𝑒 2⟩/𝑁 curves: the collapse of the polymer coincides with the rapid increase in the likelihood of end-to-end binding.