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 Interaction 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. 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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 (Generated using Python on Google Colab) 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𝜈 (figure 9). 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. 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 (Generated using Python on Google Colab) Figure 10: Time series of the cumulative means of 𝑅𝑒 2 for 10 independent runs, and their ensembleaveraged curve (Generated using Python on Google Colab) Figure 11: log ⟨𝑅𝑒 2⟩-vs-log𝑁 plot of the above simulation (Generated using Python on Google Colab) We now look at the probability distribution 𝑃(𝑅𝑒) of the end-to-end distance (figure 12). 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. Figure 12: Probability distribution of 𝑅𝑒 for three polymers of lengths 20, 40 and 80 (Generated using Python on Google Colab) We now look at plots of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ as functions of time for a smaller polymer (𝑁=20) in a larger lattice (𝐿=50); we would expect this to 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 13: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ for a polymer with N=20, L=50, 500000 Monte Carlo steps, random initial configuration (Generated using Python on Google Colab) Figure 14: Rerun of the above simulation, same parameters, same number of steps, random initial configuration (Generated using Python on Google Colab) We clearly see that now the differences between the converging value of ⟨𝑅𝑒 2⟩ over reruns of the code are minimal. 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. Figure 15: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ for a Gaussian polymer, N=20, L=50, 500000 Monte Carlo steps, random initial configuration (Generated using Python on Google Colab) Multiple reruns of the above code produced the same result – ⟨𝑅𝑒 2⟩ converges to 20 (or very close to 20), which is the value of 𝑁 we considered. We notice that smaller values of ⟨𝑅𝑒 2⟩ seem to be more heavily favored in the above plot, which makes sense for a Gaussian polymer. Since monomers are permitted to overlap, the squared end-to-end distance 𝑅𝑒 2 fluctuates asymmetrically around its mean – configurations with small end-to-end separation occur more frequently since the polymer readily folds back on itself. In particular, 𝑅𝑒 2 can occasionally reach zero, when the two ends occupy the same lattice site. At the same time, we may have less frequent, but highly stretched configurations in which a large portion of the chain becomes extended in a single direction (the spikes in the plot). These two features – frequent compact states and occasional stretched states – produce the asymmetric time series observed in the simulation. Despite these fluctuations, the cumulative mean converges to the theoretical value ⟨𝑅𝑒 2⟩∝𝑁. We see that this holds accurately for longer Gaussian polymers too (figures 16 and 17). Figure 16: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ for a Gaussian polymer, N=40, L=60, 1000000 Monte Carlo steps, random initial configuration (Generated using Python on Google Colab) Figure 17: Time series of 𝑅𝑒 2 and ⟨𝑅𝑒 2⟩ for a Gaussian polymer, N=80, L=120, 2000000 Monte Carlo steps, random initial configuration (Generated using Python on Google Colab) Notice that for Gaussian polymers, the time series of 𝑅𝑒 2 is heavily crowded near the bottom (lower values of 𝑅𝑒2), because the chain frequently folds back on itself, bringing the two ends close together. Due to the fact that overlapping is allowed in Gaussian polymers, most configurations are compact, with the ends separated by a relatively shorter distance, while occasional random fluctuations stretch the chain into elongated states that appear as sharp spikes in the 𝑅𝑒 2 time series. The end-to-end distance depends only on the two end monomers and large-scale extensions are rare; in contrast smaller values of 𝑅𝑒 2 are far more common. 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. The important point is that 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 ensemble-averaged MSD (averaged over ten independent runs with different random initial configurations), for polymers of different lengths (figures 26, 27). Figure 26: Ensemble-averaged MSD-vs-scaled time (t/N2) curve over 10 independent runs, for Gaussian polymers of lengths N=10, 20, 30, L=60 (Generated using Python on Google Colab) Figure 27: Ensemble-averaged MSD-vs-scaled time (t/N2) curve over 10 independent runs, for Gaussian polymers of lengths N=20, 40, 80, L=100 (Generated using Python on Google Colab) In the above plots, the random fluctuations in the diffusive region is 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 28) 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 28: Trajectory of the center of mass of a Gaussian polymer with N=20, under discrete Rouse dynamics (Generated using Python on Google Colab) 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 29). 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. 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 29: Ensemble-averaged MSD-vs-time curve over 10 independent runs, for Gaussian polymers of lengths N=20, 40, 80, L=120 (Generated using Python on Google Colab) The above simulations assume a Gaussian polymer. For a self-avoiding polymer, the center of mass motion exhibits slower diffusion (figure 30). Over long times, the polymer grossly follows the expected expected 𝑀𝑆𝐷∼𝑡 scaling; however there is noticeable flattening for longer polymers, owing to excluded-volume interactions and the increased effective friction on the chain, causing polymers to diffuse more slowly, consistent with the scaling 𝐷cm ∼1/𝑁. Figure 30: Ensemble-averaged MSD-vs-time curve over 10 independent runs, for self-avoiding polymers of lengths N=20, 40, 80, L=120 (Generated using Python on Google Colab) 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 31), 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 31: 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 (Generated using Python on Google Colab) 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 32), but this effect is extremely weak and the polymer still remains significantly rigid. Figure 32: 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 (Generated using Python on Google Colab) 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 33). 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 34). 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 33: 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 (Generated using Python on Google Colab) Figure 34: 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 (Generated using Python on Google Colab) 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 preserved. The answer lies in the asymmetry of the update rules themselves: the bond-length constraint is enforced only when an even monomer attempts a move, not when its odd neighbor moves. An even monomer may update its position only if it can remain exactly one lattice step away from its nearest neighbors, but an odd monomer is free to jump to any allowed lattice site without checking this condition. Consequently, whenever an odd monomer moves, it may temporarily place itself at a distance much larger than one step from its neighboring even monomers. The even monomer can restore the bond only on its own update attempt, and only if a legal position exists that re-establishes the unit distance. In many configurations such a position does not exist, and so the stretched bond persists. At each Monte Carlo step, we pick a monomer at random and attempt to move it; the type of move depends on whether the monomer that got picked has an even index or an odd index. This means in effect, the rigidity of the backbone applies in one direction (even to odd) but not in the reverse (odd to even), creating a dynamic asymmetry in which some bonds behave as fixed springs while others intermittently behave as completely flexible links. This algorithmic structure explains the emergence of fluctuating and often unphysical bond lengths, despite the apparent adjacency of odd and even monomers, and underlies the large-scale nonequilibrium expansion seen in the simulations. For good measure, it is also important to clarify another subtle aspect of the rules used in this model. We said that the odd-indexed monomers either remain in their position, or move to one of the four nearest lattice sites, or we allow them to move up to two steps. Even-indexed monomers, on the other hand, are constrained to remain at a fixed distance of one from their nearest neighbors. This means whenever an even bead moves, both of its bonds with its nearest neighbors on both sides (assuming it is not an end monomer) are restored to unit length. Oddindexed monomers, however, update under a different rule: they may not move at all, or they may hop to any nearest-neighbor site without needing to maintain a fixed bond length of one with its nearest neighbors, or it may move to any site one or two lattice steps away, again without any requirement to preserve the distances to adjacent monomers. As a result, a single odd-indexed monomer move can instantly stretch one or both nearest-neighbor bonds even if the odd monomer in question moves only one step on the lattice. As we have mentioned, the even monomers attempt to compensate for this on their own updates, but that is possible only when a legal unit-distance position actually exists, and this becomes unlikely after large oddbead jumps. Consequently, bond-length violations accumulate and persist. In order to obtain behavior that is less explosive and more physically grounded, we implement a modification of the rules of the previous model by selectively taming the dynamics of the “active” (odd-indexed) monomers. The key change is that the chain is now strictly selfavoiding, so no two monomers are allowed to occupy the same lattice site. This single constraint already makes the model more realistic: the polymer can no longer collapse into unphysical configurations or teleport through itself! The even-indexed monomers retain their original Rouse-like nearest-neighbor-distance preserving rule. However, the odd-indexed monomers are now restricted to either remaining in place or attempting a single-step nearestneighbor jump (all two-step moves removed entirely). Further and importantly, these two behaviors are not chosen equiprobably; the odd monomer remains immobile with probability 11/12≈91.7%. This simulates a weak, spatially heterogeneous environment that is more realistic compared to the strongly nonequilibrium forcing of the previous model. The resulting configurations under these updated rules evolve much more gently (figure 35). The visual snapshots show that even though bond lengths fluctuate, the polymer does not undergo the wild, unbounded distortions characteristic of the previous model. For one, selfavoidance blocks most geometrically incompatible rearrangements, and also because odd monomers now attempt moves much less often. Thus, the chain explores configuration space more slowly and locally. Yet the broken symmetry between evenand odd-indexed monomers still injects a bias into the dynamics, and the polymer never settles into a purely equilibriumlike state (detailed balance is still broken). Quantitatively, the squared gyration radius squared 𝑅𝑔 2(𝑡) continues to increase over time (figure 36), indicating that the chain is still gradually expanding under persistent nonequilibrium drive. However, the growth is now visibly milder (compare with figure 34). The polymer is still being pushed outward by the fluctuations in the bond length of the odd monomers, but the reduced mobility and self-avoidance limit how quickly this expansion can occur. Figure 35: Monte Carlo simulation of a self-avoiding polymer of length 40 under the rules of our updated alternating-rule toy model 3 on 80X80 2D lattice for random initial configuration (Generated using Python on Google Colab) Figure 36: Time series of gyration radius squared for a self-avoiding polymer of length 40 under the updated rules of our alternating-rule toy model 3 on 80X80 2D lattice for random initial configuration (Generated using Python on Google Colab) In the final alternating-rule toy model, we combine two forms of heterogeneity – variable mobility for evenand odd-indexed monomers and irregularly-spaced local immobilization of odd-indexed monomers – while restoring the realistic constraints of fixed bond length and selfavoidance. All monomers, even and odd, must maintain strict bond lengths of one lattice unit to their nearest neighbors at all times, and self-avoidance applies to both even and odd monomers. Even-indexed monomers behave exactly as in the Rouse-type dynamics considered earlier: whenever an even monomer is selected, it may move to any lattice site that simultaneously maintains unit distance to its neighbors and respects self-avoidance. Odd-indexed monomers, however, have a higher intrinsic mobility: they are chosen to attempt an update twice as frequently as even ones (this could very crudely mimic local differences in effective viscosity or friction along the chain). Moreover, when an odd monomer is selected, it has three possible behaviors: it may remain fixed with probability 0.2 (this could mimic transient local crowding or temporary binding); it may execute a standard one-step Rouse move (like the even monomers) with probability 0.1; or if allowed by the constraints, it may perform a two-step move in the same Monte Carlo attempt with probability 0.7, executing up to two consecutive legal unit moves (if this is not allowed by the constraints, the sequence is immediately aborted). Together, these rules inject asymmetry in both how often odd monomers attempt to move and how far they may move in a single attempt. Despite these heterogeneities, the global behavior of the chain remains remarkably close to an equilibrium-like state (figure 37). This is because the dynamics still enforce strict bond length preservation and self-avoidance, both of which are very strong constraints. Even though odd beads are more active and can occasionally take two consecutive steps, every accepted move must satisfy the same geometric constraints as the standard Rouse monomer. As a result, the chain’s configurational changes are predominantly local: loops expand and contract, but the global size of the polymer remains bounded. 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 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 42). 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 42: ⟨𝑅𝑒 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 (Generated using Python on Google Colab) 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. 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 43). Figure 43: ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 plots for three polymers (L=20,40,80) with end-only attraction, using the simplified two-state analytic model (Generated using Python on Google Colab) 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). 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 44). Figure 44: 𝑝𝐽-vs-𝐽 plots for three polymers (L=20,40,80) with end-only attraction, using the simplified two-state analytic model (Generated using Python on Google Colab) The above 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. The crossover point would be the point where the polymer is equally likely to be in the extended or collapsed states. Let’s denote the critical value of 𝐽 at the crossover by 𝐽𝑐. So if we increase 𝐽 beyond 𝐽𝑐, the polymer definitely collapses, and for values of 𝐽<𝐽𝑐, the polymer never collapses. At 𝐽=𝐽𝑐, we have: 𝑝𝐽=𝐽𝑐=1/2 1/2+1/2=1 2 We now solve for 𝐽𝑐 by replacing 𝐽→𝐽𝑐 and 𝑝𝐽→1/2 in the general expression for 𝑝𝐽: 𝑝0𝑒𝛽𝐽𝑐 (1−𝑝0)+𝑝0𝑒𝛽𝐽𝑐=1 2 This gives: 2𝑝0𝑒𝛽𝐽𝑐=(1−𝑝0) + 𝑝0𝑒𝛽𝐽𝑐 Simplifying, we get: 𝑝0𝑒𝛽𝐽𝑐=1−𝑝0 For simplicity, we now set 𝛽=1, and some trivial algebra gives us: 𝑒𝐽𝑐=1−𝑝0 𝑝0 Hence: 𝐽𝑐=ln(1−𝑝0 𝑝0) Substituting 𝑝0∼1/𝜋𝑁, we obtain: 𝐽𝑐=ln(𝜋𝑁−1 1) We approximate 𝜋𝑁−1≈𝜋𝑁 (which is a reasonable approximation for large and even moderate 𝑁), and using ln(1)=0, we arrive at an analytic scaling law: 𝐽𝑐≈ln(𝜋𝑁)=ln(𝜋)+ln(𝑁) This predicts that 𝐽𝑐 increases linearly with ln(𝑁), a result confirmed by 𝐽𝑐-vs-ln(𝑁) plots of the three polymers (figure 45) – the data fall on a perfect straight line with slope close to one. In figure 45, we have plotted the analytically-derived values of 𝐽𝑐 against ln(𝑁) and performed a linear regression of the form 𝐽𝑐=𝑎+ 𝑏 ln(𝑁). From the analytic scaling law, we would expect 𝑎=ln(𝜋)≈1.145, and 𝑏=1. The fitting gives 𝑎=1.119961, 𝑏=1.003994, in excellent agreement with the theoretical values. This validates the expression 𝐽𝑐≈ln(𝜋𝑁). Physically, this means that to compensate for the entropic cost of end-to-end contact, the attractive energy must grow logarithmically with polymer length. Also, since the analytic expression for 𝐽𝑐 is now known, we incorporate these threshold values directly into the visualization and replot ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽, where the corresponding 𝐽𝑐 values for all three polymers are marked as vertical dashed lines, clearly indicating where the collapse transition occurs for each polymer (figure 46). Figure 45: 𝐽𝑐-vs-𝑙𝑛(𝑁) plots for three polymers (L=20,40,80) with end-only attraction, using the simplified two-state analytic model (Generated using Python on Google Colab) Figure 46: ⟨𝑅𝑒 2⟩/𝑁-vs-𝐽 plots for the three polymers, 𝐽𝑐 visualized as vertical dashed lines, using the simplified two-state analytic model (Generated using Python on Google Colab) This simplified analytic model gives us smooth, noise-free curves for ⟨𝑅𝑒 2⟩/𝑁 as a function of 𝐽, and it beautifully captures the plateau-to-collapse behavior we expected, including the shift of the transition with chain length (longer changes collapse at higher values of 𝐽). By replacing the full Monte Carlo sampling with an explicit expression for the contact probability, the simplified two-state model sidesteps the problem of end-only contact probability being mostly negligible (even for moderately long polymer chains). The analytic expression for 𝐽𝑐 also makes the underlying physics clearer: a single attractive bond just at the ends is enough to produce a sharp, cooperative change in conformation, and the position of that crossover grows in a predictable way with 𝑁. The tiniest energetic bias at the ends propagates through the entire chain, reshaping its large-scale statistics. It is worth reemphasizing that all of this emerges while the internal monomers (all the monomers except the two end monomers!) remain unaware of any energetic law – they only obey geometry. The full transition is thus encoded in a single energetic link. Although minimal, the model captures the essential balance between entropy and end-binding, and provides a useful baseline against which more detailed interaction schemes or more realistic polymer models can be compared. Throughout this report, using the example of a polymer, we have seen how global behavior can emerge from simple local rules, and this appendix reinforces that same, fundamental idea – but in a slightly different way – by showing that even a single energetic constraint applied only at the ends of the polymer suffices to generate a clear, system-wide conformational response.