General criterion for harmonicity

TL;DR

This paper introduces a new class of transfer matrices inspired by Markov processes, providing a criterion for perfect harmonicity in systems, exemplified by a polymer model that remains harmonic despite non-Gaussian sub-units.

Contribution

It presents a novel algebraic approach to determine harmonicity in physical systems and constructs a polymer model that remains harmonic with non-Gaussian components.

Findings

Largest eigenvalue determined by explicit algebraic equation

Polymer with non-Gaussian sub-units remains harmonic until fully stretched

Confirmed by Monte Carlo and Langevin simulations

Abstract

Inspired by Kubo-Anderson Markov processes, we introduce a new class of transfer matrices whose largest eigenvalue is determined by a simple explicit algebraic equation. Applications include the free energy calculation for various equilibrium systems and a general criterion for perfect harmonicity, i.e., a free energy that is exactly quadratic in the external field. As an illustration, we construct a "perfect spring", namely a polymer with non-Gaussian, exponentially distributed sub-units which nevertheless remains harmonic until it is fully stretched. This surprising discovery is confirmed by Monte Carlo and Langevin simulations.

Equations

$$X=-d\mathcal{G}/dF.$$

$$Z$$

$${\mathbf{T}}_{{\sigma }^{\prime },\sigma }=e^{-\beta E_{{\sigma }^{\prime },\sigma }}.$$

$$\mathcal{G}=-\frac{\ln Z}{\beta }=-\frac{\mathrm{Tr}{\mathbf{T}}^N}{\beta }\sim -\frac{N\ln \lambda }{\beta },$$

$${\mathbf{A}}_{{\sigma }^{\prime },\sigma }=\left(1-q_{\sigma }\right){\delta }_{\sigma ,{\sigma }^{\prime }}+q_{\sigma }{p}_{{\sigma }^{\prime }}.$$

$${\mathbf{T}}_{{\sigma }^{\prime },\sigma }$$

$$\sum _{\sigma }{\mathbf{T}}_{{\sigma }^{\prime },\sigma }{\varphi }_{\sigma }=\lambda {\varphi }_{{\sigma }^{\prime }},$$

$$e^{-\beta {ϵ}_{{\sigma }^{\prime }}}(1-q_{{\sigma }^{\prime }}){\varphi }_{{\sigma }^{\prime }}+p_{{\sigma }^{\prime }}\sum _{\sigma }{e}^{-\frac{\beta \left({ϵ}_{\sigma }+{ϵ}_{{\sigma }^{\prime }}\right)}{2}}{q}_{\sigma }{\varphi }_{\sigma }=\lambda {\varphi }_{{\sigma }^{\prime }}.$$

$$\sum _{\sigma }{e}^{-\frac{\beta {ϵ}_{\sigma }}{2}}{q}_{\sigma }{\varphi }_{\sigma }=1.$$

$${\varphi }_{\sigma }=\frac{p_{\sigma }}{\lambda e^{\frac{\beta {ϵ}_{\sigma }}{2}}-e^{-\frac{\beta {ϵ}_{\sigma }}{2}}(1-q_{\sigma })}.$$

$$\sum _{\sigma }\frac{p_{\sigma }{q}_{\sigma }}{e^{\beta {ϵ}_{\sigma }}\lambda -(1-q_{\sigma })}=1.$$

$$\lambda =\sum _{\sigma }{p}_{\sigma }{e}^{-\beta {ϵ}_{\sigma }}.$$

$$\lambda =\frac{\sinh (bF\beta )}{bF\beta }.$$

$${\int }_0^{2\pi }d\varphi {\int }_0^{\pi }d\theta \frac{\sin \theta }{4\pi }\frac{q}{e^{bF\beta \cos \theta }\lambda -(1-q)}=1,$$

$$\lambda =\frac{\left(e^{2bF\beta /q}-1\right)(1-q)}{e^{bF\beta \left(2/q-1\right)}-e^{bF\beta }}.$$

$${\mathcal{L}}_q(y)=\frac{\left(2-q\right)e^{2y/q}\left(e^{2y}-1\right)+q\left(e^{2y}-e^{4y/q}\right)}{q\left(e^{2y/q}-1\right)\left(e^{2y}-e^{2y/q}\right)}.$$

$${\kappa }_q=\frac{q}{(2-q)}{\kappa }_1\text{with}{\kappa }_1=\frac{3}{N{b}^2\beta }.$$

$$q_{\sigma }=k_{\sigma }dl,$$

$${\mathbf{A}}_{{\sigma }^{\prime },\sigma }=1+{\mathbf{K}}_{{\sigma }^{\prime },\sigma }dl,{\mathbf{K}}_{{\sigma }^{\prime },\sigma }=k_{\sigma }\left(p_{{\sigma }^{\prime }}-{\delta }_{\sigma ,{\sigma }^{\prime }}\right).$$

$${ϵ}_{\sigma }=F_{\sigma }dl,\lambda =1+\mu dl,$$

$$\sum _{\sigma }\frac{p_{\sigma }{k}_{\sigma }}{\mu +\beta F_{\sigma }+k_{\sigma }}=1.$$

$$\mathcal{G}\sim -\frac{L\mu }{\beta },$$

$$\mu =l(\beta F{)}^2,$$

$$\sum _{\sigma }{p}_{\sigma }{e}^{-\lambda a_{\sigma }}=\sqrt{\frac{\bar{l}}{l}}\frac{I_1\left(2\sqrt{\frac{l}{\bar{l}}}\lambda \right)}{\lambda }.$$

$$x=X/L=-\frac{1}{L}\frac{d}{dF}\mathcal{G}=2lF\beta ,x\le 1.$$

$$\kappa =\frac{1}{2Ll\beta }\left(1-2\frac{dl}{\bar{l}}\right).$$

$$I_1(\lambda )=\frac{\lambda }{\pi }{\int }_0^{\pi }d\theta e^{-\lambda \cos \theta }{\sin }^2\theta ,$$

$$p_{\mathbf{\Omega }}$$

$$\frac{p_{+}{q}_{+}}{e^{\beta {ϵ}_{+}}\lambda -(1-q_{+})}+\frac{p_{-}{q}_{-}}{e^{\beta {ϵ}_{-}}\lambda -(1-q_{-})}=1.$$

$$\lambda =\frac{\left(e^{-\beta {ϵ}_{-}}\left(1-p_{+}{q}_{-}\right)+e^{-\beta {ϵ}_{+}}\left(1-p_{-}{q}_{+}\right)\right)}{2}\left(1+\sqrt{1-\frac{4e^{-\beta \left({ϵ}_{+}+{ϵ}_{-}\right)}\left(1-p_{+}{q}_{-}-p_{-}{q}_{+}\right)}{{\left(e^{-\beta {ϵ}_{+}}\left(p_{+}{q}_{-}-1\right)+e^{-\beta {ϵ}_{-}}\left(p_{-}{q}_{+}-1\right)\right)}^2}}\right).$$

Full text

The stretching of an (ideal) polymer provides one of the most beautiful illustrations of thermodynamics and equilibrium statistical physics. A force is needed because the number of polymer configurations corresponding to a stretched state is (exponentially) smaller than that of a coiled state. More generally, the extension $X$ versus force $F$ relation is obtained from the derivative, at constant temperature $T$, of the free energy $\mathcal{G}$ 111$\mathcal{G}=U-TS-FX$, with U the internal energy and S the entropy:

[TABLE]

Under isotropic conditions, one expects that the free energy $\mathcal{G}$ is even in $F$, hence quadratic for $F$ small. The extension $X$ is then harmonic, i.e., linear in the force $F$. For larger forces, the polymer is expected to stiffen as it approaches full extension. For example, for a freely-jointed chain in three dimensions, consisting of $N$ units of fixed length $b$, the fractional extension $x=X/Nb$ is given by $x=\mathcal{L}(bF\beta)$ with $\mathcal{L}$ the famous Langevin function $\mathcal{L}(y)=\coth(y)-1/y$ Kuhn and Grun (1942) and $\beta=1/(k_{B}T)$, see also Fig. 1(a).

The issue of elasticity is an important one in materials science. The width of the (harmonic) elastic regime, which can be estimated from the ratio of the tensile strength over the Young’s modulus, varies greatly from very small, for example for nanotubes, to very large for rubber and elastine Howard et al. (2001). It is therefore natural to ask whether variations of a basic microscopic model, such as the freely-jointed chain, can lead to a predominant or even perfect harmonic response. The immediate answer appears to be yes: the Rouse model Rouse Jr (1953) supposes bonds that are perfectly harmonic and hence so is the entire chain. But the assumption of perfectly harmonic bonds is unphysical as it would for example imply that both the bond and chain can be infinitely stretched. And deriving harmonicity from harmonicity is not exactly a great feat. The surprising finding of this letter is the discovery of a simple random walk model for a polymer with non-harmonic bonds which is and remains perfectly harmonic up to full extension (and not beyond).

Before proceeding to the more technical derivation, we comment on the route that led to this discovery and the additional results that were obtained. The statistical physics literature on polymers is huge, but exact results can only be derived for some very simple models such as the freely-jointed chain Flory et al. (1969); Birshtein (1966). One of the main tools to arrive at these results is the evaluation of the partition function via a transfer matrix method, essentially by identifying the largest eigenvalue. Such transfer matrices have positive entries and are therefore reminiscent of Markov matrices, which describe the dynamics of Markov chains. We introduce a special class of transfer matrices whose structure is inspired by a specific type of Markov process. We coin the name of “Kubo-Anderson” transfer matrices in referral to two early papers (on line-width problems) in which such Markov processes have been introduced KA . The bonus is that the largest eigenvalue of such a transfer matrix is determined by a simple explicit algebraic equation. Applications include the free energy calculation for various equilibrium systems, including a simple model for polymer chains with persistence. Furthermore we derive, after an additional simplification, a general criterion for perfect harmonicity, i.e., a free energy that is exactly quadratic in the external field. We will discuss it here in the context of a simple random walk model for a polymer, but the results are equally valid for other systems, such as magnetic systems with exactly linear susceptibility. The application to simple polymer models leads to the discovery of the “perfect spring”, i.e., a random walk model for a polymer with non-Gaussian sub-units which nevertheless remains harmonic until it is fully stretched.

As is well known Baxter (2007), many problems in equilibrium statistical mechanics, including the celebrated Onsager solution of the two-dimensional Ising model, can be formulated in terms of a transfer matrix. One supposes that the energy of the system can be written as a sum $E=\sum_{i=1}^{N}E_{\mathbf{\sigma}_{i+1},\mathbf{\sigma}_{i}}$, where $\mathbf{\sigma}_{i}$ represents the state of the “i-th layer” and $E_{\mathbf{\sigma}_{i+1},\mathbf{\sigma}_{i}}$ is the interaction energy between layers $i+1$ and $i$. For notational simplicity, we consider periodic boundary conditions with layer $N+1$ identified with layer $1$. The central quantity is the partition function:

[TABLE]

where the sum over $\{{\mathbf{\sigma}}\}$ runs over all configurational states of the system. The transfer matrix $\mathbf{T}$ is defined by its elements:

[TABLE]

In the thermodynamic limit $N\rightarrow\infty$, the free energy $\mathcal{G}$ is obtained as:

[TABLE]

where $\sim$ refers to an equality to dominant order in $N$, and $\lambda$ is the largest eigenvalue of the transfer matrix $\mathbf{T}$.

A transfer matrix has positive elements. The elements of a Markov matrix represent probabilities, which are obviously also positive but obey in addition a normalization condition. The relation between transfer matrices and Markov matrices has been noticed a long time ago, see, e.g., Miller (1961); Fisher (1964), and has been revisited more recently in the context of large deviations for conditioned Markov processes Chetrite and Touchette (2015). In this letter, we introduce a new class of transfer matrices inspired by “Kubo-Anderson” Markov processes. The latter have been used for a detailed analytic description of a large variety of physical processes KA . In the context of transfer matrices, the bonus is a simple explicit relation for its largest eigenvalue. Our starting point is the following Markov matrix:

[TABLE]

$\mathbf{A}_{\mathbf{\sigma}^{\prime},\mathbf{\sigma}}$ is the transition probability to go from state $\mathbf{\sigma}$ to $\mathbf{\sigma}^{\prime}$. Its explicit form can be explained as follows. With a probability $1-q_{\mathbf{\sigma}}$ the system remains in its present state $\mathbf{\sigma}$. With probability $q_{\mathbf{\sigma}}$ a novel state is selected. This novel state is $\mathbf{\sigma}^{\prime}$ with probability $p_{\mathbf{\sigma}^{\prime}}$. As is explained in the introduction with the stretching of a polymer, the addition of an external field allows to explore the regions of “higher free energy”. We do something similar here by adding an extra Boltzmann factor to the transfer matrix, with $\epsilon_{\mathbf{\sigma}}$ representing the energy contribution due to such a field on layer $\mathbf{\sigma}$, leading finally to the following “Kubo-Anderson” transfer matrix:

[TABLE]

The connection with the corresponding interaction energy $E_{{\mathbf{\sigma}}^{\prime},{\mathbf{\sigma}}}$ is obtained by comparison with Eq. (3). Conversely, this relation establishes the dependence of the quantities $q_{\mathbf{\sigma}}$ and $p_{\mathbf{\sigma}}$ on the energies $E_{{\mathbf{\sigma}}^{\prime},{\mathbf{\sigma}}}$ and the temperature.

We now show how the largest eigenvalue (and corresponding eigenfunction) of a “Kubo-Anderson” transfer matrix $T$ can be obtained. From the eigenvalue equation

[TABLE]

one finds that

[TABLE]

The largest eigenvalue $\lambda$ and its corresponding eigenfunction $\mathbf{\phi}$ can be identified by invoking the Perron-Frobenius theorem: the eigenvalue $\lambda$ is unique and positive and all components of $\mathbf{\phi}$ have the same sign, and can thus be chosen to be positive. As a consequence $\sum_{{\mathbf{\sigma}}}\exp\left(-\beta\epsilon_{\mathbf{\sigma}}/2\right)q_{\mathbf{\sigma}}\mathbf{\phi}_{\mathbf{\sigma}}$ is also positive. Furthermore, since $\mathbf{\phi}$ is only determined up to a constant factor, we can assume the following normalisation:

[TABLE]

With this constraint, we find from Eq. (8) the following explicit expression for the eigenvector $\mathbf{\phi}$:

[TABLE]

By substitution of this expression in Eq. (9), one concludes that $\lambda$ is determined by:

[TABLE]

This simple explicit algebraic equation for the dominant eigenvalue is the first important result of this letter.

We mention a few classes of systems which can be solved exactly. As a first example, we consider systems without “persistence”, i.e., $q_{\mathbf{\sigma}}\equiv 1$, leading to:

[TABLE]

The 3-d freely-jointed chain is obtained with the following identifications: $\mathbf{\sigma}$ represents the space angle $\mathbf{\Omega}$, specifying the orientation of each subunit (of fixed length $b$). Since this orientation is random, one has $p_{\mathbf{\Omega}}=1/\left(4\pi\right)$, while the sum becomes an integral over the space angle, $\sum_{\mathbf{\sigma}}\rightarrow\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta\sin\theta$. Furthermore $\epsilon_{\mathbf{\Omega}}=Fb\cos\theta$ represents the effect of an external field, with $\theta$ the angle between the monomer and this field. Eq. (12) gives the result:

[TABLE]

By combination with Eq. (1) and Eq. (4), one recovers the aforementioned result $x=X/(Nb)=\mathcal{L}(bF\beta)$.

The above analysis can be reproduced for state-independent persistence, i.e., $q\leq 1$ but independent of $\mathbf{\sigma}=\mathbf{\Omega}$, by starting from Eq. (11) rather than from Eq. (12). The integral determining the eigenvalue $\lambda$,

[TABLE]

can still be solved:

[TABLE]

The fractional extension $x=\mathcal{L}_{q}(bF\beta)$ is described by a generalized Langevin function, see also Fig. 1(a):

[TABLE]

The freely-jointed chain is retrieved in the limit $q\rightarrow 1$, with $\mathcal{L}_{q}$ converging to $\mathcal{L}$. In the limit of strong persistence, $q\rightarrow 0$, $\mathcal{L}_{q}$ converges to the sign function $\mathcal{L}_{0}(y)=\textrm{sgn}(y)$. For small forces, the polymer behaves as a harmonic spring, with spring constant $\kappa_{q}=F/X$ given by:

[TABLE]

$\kappa_{1}$ is the spring constant of the 3-d freely-jointed chain. Note the weakening of the spring for increasing persistence, corresponding to a decreasing value of $q$. These predictions have been verified using Langevin simulations, cf. Fig. 1(a). Other applications, e.g., 1-d polymers and systems with only two $\mathbf{\sigma}$-states, are presented in the supplemental material.

Depending on the interpretation of the model, the above transfer matrix describes discrete steps taking place in space (for example in an Ising spin or polymer chain problem), in time (for a Markov chain) or in another possibly more abstract coordinate (for example an angle coordinate or a higher dimensional vectorial coordinate). To simplify further the eigenvalue equation, we focus on a “hydrodynamic” limit. For simplicity, we present it as taking place in the context of a single scalar spatial variable. We associate an elementary spatial displacement of length $dl$ (corresponding to the bond length $b$ in the above polymer problem) to each discrete step. This length will be small compared to the average length of “straight” segments, belonging to a given state $\mathbf{\sigma}$, $\bar{l}_{\mathbf{\sigma}}$, i.e., $dl/\bar{l}_{\mathbf{\sigma}}\rightarrow 0$. Meanwhile, the total length $L=Ndl$ will become large compared to the typical lengths of the problem, $\bar{l}_{\mathbf{\sigma}}/L\rightarrow 0$. In this limit, the “jump probabilities” $q_{\mathbf{\sigma}}$ are replaced by transition probabilities per unit length $k_{\mathbf{\sigma}}=1/\bar{l}_{\mathbf{\sigma}}$:

[TABLE]

resulting in “straight” segments with lengths that are exponentially distributed. The matrix $\mathbf{A}$ is replaced by a transition matrix $\mathbf{K}$:

[TABLE]

Consistent with this limit, the energy $\epsilon_{\mathbf{\sigma}}$ and the eigenvalue $\lambda$ converge as follows to [math] and $1$, respectively:

[TABLE]

with $F_{\mathbf{\sigma}}$ representing the external “force” or “energy density” when the system is in the state $\mathbf{\sigma}$. The eigenvalue equation Eq. (11) thus reduces to the following algebraic relation for $\mu$:

[TABLE]

We note that the dependence on the external field contribution $F_{\mathbf{\sigma}}$ is no longer via an exponential function. The corresponding free energy, cf. Eq. (4), becomes:

[TABLE]

with a simple proportionality to the eigenvalue $\mu$. The extensivity in $N$ is replaced by an extensivity in $L$. Eq. (22) is an equality to dominant order in $L$.

We now turn to the search for a “perfect spring”, defined as a system for which the free energy Eq. (22) is exactly quadratic in the external field amplitude. More precisely, we introduce the overall force $F$ via the specification $F_{\mathbf{\sigma}}=a_{\mathbf{\sigma}}F$, with $a_{\mathbf{\sigma}}$ an $F$-independent amplitude, and require that

[TABLE]

with $l$ an $F$-independent reference length scale. The question thus reduces to finding probability distributions $p_{\mathbf{\sigma}}$ such that Eq. (21) holds under this constraint. By Taylor expansion in $F$, and under the assumption that $k_{\mathbf{\sigma}}={1}/{\bar{l}}$ is independent of the state $\mathbf{\sigma}$, one finds the following explicit expression for the moment generating function associated with $p_{\mathbf{\sigma}}$ (see supplemental materials):

[TABLE]

$I_{1}$ is the modified Bessel function. This general criterion for harmonicity is our second major result. The probability distribution $p_{\mathbf{\sigma}}$ can be found from it, depending on the topology of the phase space and the form of $a_{\mathbf{\sigma}}$, by an inverse (integral) transform. Eq. (24) obviously has no solution for a finite state space of $\mathbf{\sigma}$. In particular, two-state models (corresponding for example to a polymer model in $d=1$) can not be turned into fully harmonic springs.

Combining Eq. (1), Eq. (22) and Eq. (23), one finds that the corresponding stretching fraction $x$ is given by:

[TABLE]

The extension is thus exactly linear in $F$. The above result is only valid up to $x=1$, i.e., until the polymer is fully stretched. The reason for this limitation is that the Taylor expansion of Eq. (21) has a radius of convergence given by $F_{c}=1/(2l\beta)$. For larger values of $F$, $x$ stays put at its maximal value $x=1$, hence the system undergoes a second order phase transition at $F=F_{c}$ (discontinuous second derivative of $\mathcal{G}$). The spring constant $\kappa$, corresponding to the harmonic law Eq. (25), is given by $\kappa={1}/(2Ll\beta)$. The above result has been derived in the limit $dl/\bar{l}\rightarrow 0$. In the supplemental material, we evaluate the first order correction in $dl$ and conclude that the harmonic behavior prevails, but with a modified spring constant:

[TABLE]

As a concrete application of the above harmonicity criterion, we return to the polymer problem in d-dimensional Euclidean space, with the identification of $\mathbf{\sigma}$ as a d-dimensional spatial angle $\mathbf{\Omega}$. Identifying $b$ with $dl$, one finds from $\epsilon_{\mathbf{\Omega}}=Fb\cos\theta$ that $a_{\mathbf{\Omega}}=\cos\theta$. If we furthermore assume that $p_{\mathbf{\Omega}}$ only depends on $\theta$, the integral Eq. (24) can be solved by inverting the integral transform. Anticipating that $l=\bar{l}/4$ and comparing the harmonicity criterion, Eq. (24), with the following integral representation of the Bessel function Weisstein (2002):

[TABLE]

we conclude (remembering that the Jacobian of the $d$-sphere features the factor $\sin^{d-2}\theta$):

[TABLE]

($\mathcal{N}_{d}$ is a normalisation constant.) This is our third major result. We conclude that a polymer with persistence, consisting of exponentially distributed straight segments, is perfectly harmonic until full stretching in $d=4$. In $d=3$ one needs an additional field, orthogonal to the stretching direction, which induces a biased distribution $p_{\mathbf{\Omega}}\sim{\sin\theta}$. Such a field can be realised by application of an electromagnetic force, which is often used in the experimental stretching of polymers. We have verified the latter prediction via Monte Carlo and Langevin simulations. The numerical results are in perfect agreement with the theory, cf. Fig. 1(b) and supplemental material for more details. These conclusions are of course not restricted to polymer models.

In conclusion, we have introduced “Kubo-Anderson” transfer matrices, for which the largest eigenvalue can be obtained from the simple, explicit algebraic equation, Eq. (11). The latter simplifies upon taking a continuous limit to Eq. (21). By assuming a uniform transition rate, a simple, explicit criterion for harmonicity results, cf. Eq. (24). As a concrete example, we show that a polymer chain with persistence can behave, until fully stretched, as a perfect harmonic spring. Although our model is rather theoretical, and perhaps not experimentally feasible, its discovery can serve as the starting point for the construction of complex polymer systems with enhanced harmonicity. These results can be easily mapped on other systems. For example, one could construct an ising-like magnet with perfectly linear susceptibility.

I Two state model

For a two-state system, $\sigma=\pm$, Eq. (11) for the eigenvalue $\lambda$ reduces to:

[TABLE]

The relevant solution of this quadratic equation in $\lambda$ is given by:

[TABLE]

The free energy follows by combination with Eq. (4), and the force versus extension relation is then obtained from Eq. (1).

As an example, we consider an unbiased random walker with persistence. Setting $p_{+}=p_{-}=1/2$, $q_{+}=q_{-}=q$ and $\epsilon_{+}=-\epsilon_{-}=Fb$ in the above expression, one finds:

[TABLE]

The resulting fractional extension reads:

[TABLE]

For small forces, one can associate a spring constant to this system equal to

[TABLE]

In the absence of persistence, i.e., in the limit $q\rightarrow 1$, one reproduces the $d=1$ random walk result $x=\tanh\left(bF\beta\right)$, while $x=\textrm{sgn}(bF\beta)$ in the ”complementary” limit $q\rightarrow 0$.

II Temperature and Energy dependence

To illustrate the dependence of the parameters $q_{\mathbf{\sigma}}$ and $p_{\mathbf{\sigma}}$ on energy and temperature, we consider the case of constant $q_{\mathbf{\sigma}}=q$ and $p_{\mathbf{\sigma}}=p$, independent of ${\mathbf{\sigma}}$. Introducing the interaction energies $E^{+}$, upon staying in the same state, and $E^{-}$, upon switching states, both again assumed to be independent of the state (and $E^{+}\leq E^{-}$), one readily finds:

[TABLE]

III Freely-jointed 4d chain

Applying the eigenvalue Eq. (11) to the four dimensional freely jointed chain, one finds

[TABLE]

For the case without persistence, $q=1$, the integral simplifies to

[TABLE]

and the fractional extension becomes

[TABLE]

By considering the continuous limit, using Eq. (21), as departure point, one finds for the eigenvalue $\mu$:

[TABLE]

Using the following definite integral:

[TABLE]

we recover the result of the main text:

[TABLE]

IV Finite $dl$ correction

To show the robustness of the perfectly harmonic regime, we evaluate the effect of a finite value of $dl/\bar{l}$, by calculating a first order correction. We expand Eq. (11) for $\lambda$, as in the main text. $\lambda$ is expanded one order further in $dl$:

[TABLE]

Making a Taylor expansion up to second order in $dl$, one finds:

[TABLE]

The first equation determines the value of $\mu$, For the perfect spring model, discussed in the main text, this gives:

[TABLE]

under the assumption that $l$ is inversely proportional to $k$, $l=\alpha/k$. If we assume that $k_{\mathbf{\sigma}}=k$, and that $k,p_{\mathbf{\sigma}}$ and $F_{\mathbf{\sigma}}$ are independent of $\beta$, one can show, by taking the derivatives of the first equation to $k$ and $\beta$,

[TABLE]

[TABLE]

This allows to rewrite Eq. (S15), determining the value of $\nu$ as

[TABLE]

giving

[TABLE]

From Eq. (S16), one deduces

[TABLE]

and therefore,

[TABLE]

or

[TABLE]

V Moments of perfect spring

The set of equations for the moments of the probability distribution for the harmonicity criterion are found by expanding:

[TABLE]

in a Taylor series of $F$ followed by a binomial expansion:

[TABLE]

By regrouping terms with the same exponent of $\beta F$ we derive the following set of conditions:

[TABLE]

with the brackets referring to an average with respect to $p_{\mathbf{\sigma}}$. A situation of particular interest is the case where $k_{\mathbf{\sigma}}={1}/{\bar{l}}$ is independent of the state $\mathbf{\sigma}$. $\bar{l}$ is the (now $\mathbf{\sigma}$-independent) average persistence length. This set of equations can be solved using the following mathematical identity:

[TABLE]

which immediately implies

[TABLE]

while all odd moments $\langle a^{2n+1}_{\mathbf{\sigma}}\rangle$ are all zero. The associated moment generating function is given by

[TABLE]

Identification with the following expansion of the modified Bessel function:

[TABLE]

leads to the result Eq. (24) of the main text.

VI Monte Carlo simulations

We first generate an array of length $N$ with the directions of the individual monomers chosen at random. At each subsequent step of the Monte Carlo simulation, we randomly choose an element $i$ from the array, and update its state according to the ”familiar” Monte Carlo rules

• •

If $i=1$ or $N$, the site will be in the same state as the neighbouring site with probability

[TABLE]

and in a state drawn from the distribution

[TABLE]

with probability $1-p_{\theta}$.

• •

If the site on the left of $i$ is in the same state, $\theta$, as the site on the right of $i$, $i$ will also be set to this state with probability

[TABLE]

• •

If the state on the left and on the right of $i$ differ, the probability to be set in the state $\theta_{i-1}$ or $\theta_{i+1}$ is equal to

[TABLE]

and with probability $1-p_{\theta_{i-1}}-p_{\theta_{i+1}}$ from the probability distribution Eq. (S34).

This procedure is repeated until convergence, typically for $10^{7}L/(kdl)^{2}$ steps.

VII Langevin simulations

We apply the Langevin formalism to construct a numerical scheme that achieves the most experimentally viable realisation of a perfect spring. We model the spring as a chain of $N+1$ beads with mass $m$ and position vector $\vec{x}_{i}$, with $i=0,1,\dots,N$. We take $\vec{x}_{i}\in{\rm I\!R}^{4}$, this reduces the complexity of the simulation and does not influence the resulting steady state behaviour. The beads interact through a potential $U(\vec{x}_{0},\vec{x}_{1},\dots,\vec{x}_{N})=U(\{\vec{x}_{i}\})$. The Langevin equation describes the molecular dynamics of the beads while they are in contact with a viscous heat-bath with friction coefficient $\gamma$ and at temperature $T$. We assume overdamped dynamics where $m/\gamma\rightarrow 0$, the Langevin equation is then given by

[TABLE]

The second term on the right-hand side, $\vec{\xi}_{i}(t)$, represents the random collisions of the bath particles. This is a stochastic process (uncorrelated between beads) that is Gaussian distributed with zero mean and a correlation that is determined by the fluctuation-dissipation relation: $\langle\xi_{i,\alpha}(t)\xi_{j,\beta}(t^{\prime})\rangle=2\gamma k_{B}T\delta(t-t^{\prime})\delta_{i,j}\delta_{\alpha,\beta}$, with $k_{B}$ the Boltzmann constant and the second (Greek) indices indicates the specific components of the vector.

Notice that this description introduces, in comparison with the theory, a new fundamental parameter $\gamma$, i.e. the friction coefficient. This parameter will, however, only affect the transient dynamics of the spring and not the steady state behaviour which is the relevant quantity for our discussion.

The potential $U(\{\vec{x}_{i}\})$ should describe a perfect spring, it consists of two components: $U=U_{h}+U_{p}$. First, the distance between neighbouring beads should maintain a constant value $b$. We achieve this by linking them with a stiff harmonic spring with rest-length $b$ and a very large spring constant $k$. We have

[TABLE]

with $\vec{r}_{i}=(\vec{x}_{i}-\vec{x}_{i-1})$ a bond-vector and $k\gg 1$. When $U_{h}$ is the only contribution to the total potential $U(\{\vec{x}_{i}\})$, the spring is called a freely-jointed chain.

If, however, we want to describe a perfect spring, the following potential should be added

[TABLE]

with $\beta=(k_{B}T)^{-1}$, $\theta=\cos^{-1}(\hat{r}_{i}\cdot\hat{r}_{i+1})$ and $q$ the chance of a bond to change direction. The Dirac delta function $\delta(\cdot)$ in this equation can not be implemented effectively in a numerical scheme, we therefore consider the following Gaussian which defines the Dirac delta function

[TABLE]

This introduces a parameter $\epsilon$ which should be as small as possible. The “perfect” potential can thus be approximated by

[TABLE]

with $\epsilon\ll 1$. We also used $(\hat{r}_{i}-\hat{r}_{i+1})^{2}=2-2\cos(\theta_{i})$.

The (first order) numerical integrator we used to find the time-evolution of the beads, i.e. $\vec{x}_{i}(t)$, resulting from equation (S38) is given by Brańka and Heyes (1999)

[TABLE]

The stochastic term $\eta(t)$ is a Gaussian random variable with zero mean and a correlation given by $\langle\eta(t)\eta(t^{\prime})\rangle=2k_{B}T\delta(t-t^{\prime})\Delta t/\gamma$. Due to the large forces arising from small $\epsilon$ and large $k$, one should take the time-step $\Delta t$ appropriately small to avoid large numeric errors.

To examine the performance of our numerical scheme, we devised two tests. Firstly, we should check that the bond-lengths do not deviate greatly from the $b$-value. In the left graph of figure S1, we plot the evolution of the relative bond-length, $|\vec{r}_{i}|/b$, of one bond as a function of time. When taking $k=5\cdot 10^{4}$ and $\Delta t=10^{-6}$, we find that the relative bond-length never exceeds a deviation of more than $10\%$ from its desired value of one. Secondly, one can investigate the distribution $\rho(\theta)$ of the angles between two consecutive bonds. From the “perfect” potential $U_{p}$ this distribution can be found using the Boltzmann factor: $\rho(\theta)=\mathcal{N}\exp(-\beta U_{p})\sin^{2}(\theta)$, with $\mathcal{N}$ a normalisation factor. Note that here we take $U_{p}$ to be between only two bond-vectors, i.e. equation (S42) without the summation. It is easy to see that $\mathcal{N}=1$, we therefore find

[TABLE]

This distribution is plotted (full line) in the right graph of figure S1, for $\epsilon=b=q=0.1$. From this figure it is clear that the potential is sharply peaked around small angles, while it has a smaller - yet free - distribution for larger angles. Simulated distributions are also plotted in this figure for $\Delta t=10^{-6}$ and $\Delta t=10^{-7}$. One can see that a smaller time-step yields a better agreement with the analytical curve, while for larger $\Delta t$ the distribution of the bulk angles increases. This is due to the inaccuracy of the finite integration step, which is only completely resolved for $\Delta t\rightarrow 0$. To reach steady state in a reasonable computation time the time-step can not be too small, we therefore choose to work with $\Delta t=10^{-6}$. The deviation from the desired probability distribution can be fixed by assuming an effective $q$. From figure S1, one can see that if we plot the distribution (dashed line) for $q=0.23$, it corresponds exactly with the less-accurate simulated distribution of $q=0.1$. We therefore accept $q=0.23$ as the effective $q$ of the simulated perfect spring.

In order to extract the force-extension relation of the perfect spring, we need to introduce two new forces in equation (S38). We add a constant force $F$ in the $x$-direction to one of the end-beads (for example to $\vec{x}_{N}$), its integration step becomes

[TABLE]

The other end-bead, $\vec{x}_{0}$, should remain fixed in the origin. We therefore assume that any force it experiences from the chain or the heat-bath will be counteracted by the origin. Or in other words, there is no net force acting upon this bead, so

[TABLE]

The initial configuration of the spring is chosen random, only enforcing the constant $b$-length of the bonds. When the spring is allowed to evolve under the influence of the constant force, it will first go through a transient regime where it increases its extension. Thereafter, it arrives in a steady state where its average extension does not chance. The average extension $X(t)$ through time is computed as follows

[TABLE]

where $L_{x}(t)=x_{N,x}(t)-x_{0,x}(t)$ and $L=|\vec{x}_{N}(t)-\vec{x}_{0}(t)|$. The averaging $\langle\cdot\rangle$ is done over different realisations of the spring. When the spring has reached steady state we write the average extension as $X=X(t\gg 1)$.

From the time-plot (left graph in figure S2) of the averaged extension, one can see that for all shown forces, $X(t)$ is able to reach steady state. When in steady state, we can obtain $X$ from the simulation. To acquire even more averaging, we also take different (uncorrelated) samples from the steady state regime. The right graph in figure S2 (see also main text) shows the force-extension for different values of $N$. For small forces they correspond well with the theoretical prediction (and even better with the first order correction). Deviations for larger forces are due to the smoothing of the delta function (i.e. $\epsilon\neq 0$) and the fact that $N$ is finite.