Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy

Abstract

Non-local and non-convex energies represent fundamental interacting effects regulating the complex behavior of many systems in biophysics and materials science. We study one-dimensional, prototypical schemes able to represent the behavior of several biomacromolecules and the phase transformation phenomena in solid mechanics. To elucidate the effects of thermal fluctuations on the non-convex non-local behavior of such systems, we consider three models of different complexity relying on thermodynamics and statistical mechanics: (i) an Ising-type scheme with an arbitrary temperature-dependent number of interfaces between different domains, (ii) a zipper model with a single interface between two evolving domains, and (iii) an approximation based on the stationary phase method. In all three cases, we study the system under both isometric condition (prescribed extension, matching with the Helmholtz ensemble of the statistical mechanics) and isotensional condition (applied force, matching with the Gibbs ensemble). Interestingly, in the Helmholtz ensemble the analysis shows the possibility of interpreting the experimentally observed thermal effects with the theoretical force–extension relation characterized by a temperature-dependent force plateau (Maxwell stress) and a force peak (nucleation stress). We obtain explicit relations for the configurational properties of the system as well (expected values of the phase fractions and number of interfaces). Moreover, we are able to prove the equivalence of the two thermodynamic ensembles in the thermodynamic limit. We finally discuss the comparison with data from the literature showing the efficiency of the proposed model in describing known experimental effects.

Appendices

Appendix A: Non-local behavior: relation between the next-to-nearest-neighbor (NNN) interaction strategy and the Ising scheme

The non-local interactions in discrete elastic chains are typically introduced through next-to-nearest-neighbor (NNN) elastic elements [51, 54]. Under specific hypotheses, we prove here that the scheme with the NNN elements can be reconducted to a typical Ising model, widely adopted in previous studies on similar topics [63, 64]. To begin with, we consider a chain of bistable units with additional NNN linear springs to model non-local effects [54]. The Hamiltonian reads

$$\begin{aligned} H =&\sum _{i=1}^{N}\left\{ Q(S_i)+\frac{K(S_i)\ell ^2}{2}\left[ \lambda _i-\lambda _0(S_{i})\right] ^2\right\} \nonumber \\&+\sum _{i=1}^{N-1}\frac{R\ell ^2}{2}[\lambda _{i+1}+\lambda _{i}-\lambda _0(-1)-\lambda _0(+1)]^2, \end{aligned}$$

(A1)

where R is the elastic constant of the NNN springs, and the other quantities are defined in Sect. 2. In particular, we name here $H_0$ the energy corresponding to the bistable nearest neighbor (NN) elements and $H_I$ the energy of the NNN linear springs, such that $H=H_0+H_I$. We assume that the equilibrium length of the NNN elements is fixed at $L(-1)+L(+1) = (1+\chi )\ell$, and we introduce two different behaviors: an antiferromagnetic one when $R>0$, and a ferromagnetic one when $R<0$. In the former case ($R>0$), two adjacent units entail a lower energy if they are in two different states (folded and unfolded) while, in the latter one ($R<0$), two adjacent units result in a lower energy if they are in the same state (either both folded or both unfolded), as shown in Fig. 15.

Fig. 15

Potential energy U of a linear spring representing an arbitrary NNN element. The two cases correspond to the antiferromagnetic ($R>0$) and ferromagnetic ($R<0$) behaviors of the chain with the NNN linear elastic elements

Supposing that the elastic constants k and $\alpha k$ are sufficiently large, the Hamiltonian $H_I$ can be simplified by assuming that the lengths of the units can be approximated with the equilibrium lengths of the explored wells. Then, by using the relations $\lambda _{i+1}\simeq \lambda _0(S_{i+1})$ and $\lambda _{i}\simeq \lambda _0(S_{i})$, we get

$$\begin{aligned} \lambda _{i+1}+\lambda _{i} \simeq \lambda _0(S_{i+1})+\lambda _0(S_{i}) =-\frac{S_{i+1}}{2}(1-\chi )-\frac{S_{i}}{2}(1-\chi )+(1+\chi ). \end{aligned}$$

(A2)

Thus, the interaction Hamiltonian becomes

$$\begin{aligned} H_I&\simeq \sum _{i=1}^{N-1}\frac{R\ell ^2}{2}\left[ -\frac{(1-\chi )}{2}(S_{i+1}+S_i)\right] ^2\nonumber \\&= \frac{R}{4}\ell ^2(1-\chi )^2(N-1) - J\sum _{i=1}^{N-1}S_{i+1}S_i, \end{aligned}$$

(A3)

where

$$\begin{aligned} J=-\frac{R}{4}\ell ^2(1-\chi )^2. \end{aligned}$$

(A4)

This proves that we can approximate the behavior of the chain with the NNN elements by means of a classical Ising chain. We observe that when $J>0$ ($R<0$) we are in a ferromagnetic case, and when $J<0$ ($R>0$) in the antiferromagnetic one. For this reason, in this work we adopted the overall Hamiltonian given by Eq. (4), where we neglected the irrelevant constant in $H_I$.

Appendix B: Helmholtz and Gibbs ensembles in the Ising Model

In this appendix, we show the details of the calculation of the Gibbs and Helmholtz partition functions.

1.1 B.1 Gibbs ensemble

For the Gibbs ensemble, using Eq. (5), the partition function can be evaluated as

$$\begin{aligned} \begin{aligned} Z_G({\tilde{f}}) =\,\,&\ell ^N\displaystyle {\sum _{\{S_i\}}}\int _{{\mathbb {R}}^{N}} e^{-{\tilde{\beta }}({\tilde{H}}-{\tilde{f}}\sum _{i=1}^N\lambda _i)}\mathrm {d}\lambda _1 \dots \mathrm {d}\lambda _N\\ =\,\,&\ell ^N\displaystyle {\sum _{\{S_i\}}} e^{-{\tilde{\beta }}\left( \sum _{i=1}^{N} {\tilde{Q}}(S_i)-\sum _{i=1}^{N-1} S_iS_{i+1}\right) } \\&\times \displaystyle {\int _{{\mathbb {R}}^{N}}}e^{ -{\tilde{\beta }}\sum _{i=1}^{N}\left[ \frac{{\tilde{K}}(S_i)}{2}\left( \lambda _i -\lambda _0(S_{i})\right) ^2- {\tilde{f}}\lambda _i\right] } \mathrm {d}\lambda _1\dots \mathrm {d}\lambda _N, \end{aligned} \end{aligned}$$

(B5)

where ${\tilde{f}}$ is the dimensionless force applied to the last unit of the chain. Each sum on $S_i$ ($i=1,\dots ,N$) must be interpreted as a sum over the values $+1$ and $-1$. By a Gaussian integration we obtain

$$\begin{aligned} Z_G({\tilde{f}}) =\ell ^N\sum _{\{S_i\}}e^{{\tilde{\beta }}\sum _{i=1}^{N-1}S_iS_{i+1}}\prod _{i=1}^{N}\sqrt{\frac{2\pi }{{\tilde{\beta }} {\tilde{K}}(S_i)}}e^{{\tilde{\beta }}\left( \frac{{\tilde{f}}^2}{2{\tilde{K}}(S_i)}+\lambda _0(S_{i}){\tilde{f}}-{\tilde{Q}}(S_i)\right) }. \end{aligned}$$

(B6)

We can define

$$\begin{aligned} c(S_i)=\sqrt{\frac{2\pi }{{\tilde{\beta }} {\tilde{K}}(S_i)}}e^{{\tilde{\beta }}\left( \frac{{\tilde{f}}^2}{2{\tilde{K}}(S_i)}+\lambda _0(S_{i}){\tilde{f}}-{\tilde{Q}}(S_i)\right) }, \end{aligned}$$

(B7)

so that we obtain

$$\begin{aligned} \begin{aligned} Z_G({\tilde{f}})= \ell ^N\displaystyle {\sum _{\{S_i\}}} \sqrt{c(S_1)}\left[ \prod _{i=1}^{N-1}e^{{\tilde{\beta }} S_i S_{i+1}} \sqrt{c(S_i)c(S_{i+1})}\right] \sqrt{c(S_N)}. \end{aligned} \end{aligned}$$

(B8)

In order to explicitly evaluate the summation, we can use the transfer matrix method [86]. We obtain

$$\begin{aligned} Z_G({\tilde{f}})={\varvec{w}}^\intercal {\varvec{T}}^{N-1}{\varvec{w}}, \end{aligned}$$

(B9)

where we have defined the transfer matrix ${\varvec{T}}$ and the vector ${\varvec{w}}$ (taking care of the boundary conditions) as follows:

$$\begin{aligned} \begin{aligned} {\varvec{T}}&=\begin{bmatrix} e^{{\tilde{\beta }} }c_- &{} e^{-{\tilde{\beta }} }\sqrt{c_+c_-} \\ e^{-{\tilde{\beta }} }\sqrt{c_+c_-} &{} e^{{\tilde{\beta }} }c_+ \end{bmatrix}, \\ {\varvec{w}}&=\begin{pmatrix}\sqrt{c_-}\\ \sqrt{c_+ } \end{pmatrix}, \end{aligned} \end{aligned}$$

(B10)

with, see Eq. (B7),

$$\begin{aligned} c_+\triangleq c(+1), \quad c_-\triangleq c(-1). \end{aligned}$$

(B11)

By using the standard matrix functions theory [128, 129], we have

$$\begin{aligned} {\varvec{T}}^{N-1}=\frac{{{\hat{\lambda }}_1}^{N-1}-{\hat{\lambda }_2}^{N-1}}{{{\hat{\lambda }}_1}-{\hat{\lambda }_2}}{\varvec{T}}+\frac{{{\hat{\lambda }}_1}{\hat{\lambda }_2}^{N-1}-{{\hat{\lambda }}_1}^{N-1}{{\hat{\lambda }}_2}}{{\hat{\lambda }_1}-{{\hat{\lambda }}_2}}{\varvec{I}}, \end{aligned}$$

(B12)

where ${\varvec{I}}$ is the $2\times 2$ identity matrix and $\hat{\lambda }_{1,2}$ are the eigenvalues of ${\varvec{T}}$, namely

$$\begin{aligned} {\hat{\lambda }}_{1,2}= & {} \frac{e^{{\tilde{\beta }}}}{2} \left[ c_+ + c_- \pm \sqrt{ \left( c_+-c_-\right) ^2+4 c_+c_-e^{-4{\tilde{\beta }}}} \right] \nonumber \\= & {} \frac{e^{{\tilde{\beta }}}}{2} \left( c_+ + c_- \pm \sqrt{ \Delta } \right) . \end{aligned}$$

(B13)

In Eq. (B13), ${{\hat{\lambda }}_1}$ (${{\hat{\lambda }}_2}$) corresponds to the $+$ (−) sign and we have also defined

$$\begin{aligned} \Delta =\left( c_+-c_-\right) ^2+4 c_+c_-e^{-4{\tilde{\beta }}}. \end{aligned}$$

(B14)

By substituting ${\hat{\lambda }}_{1,2}$ into Eqs. (B9) and (B12), we get the partition function given in Eq. (8).

1.2 B.2 Helmholtz ensemble

We consider here the case with fixed $x_N$ (isometric condition), described by the Helmholtz ensemble, and evaluate the canonical partition function $Z_H(x_N)$. The partition functions in the Gibbs and Helmholtz ensembles are linked by a Laplace transform [87]

$$\begin{aligned} Z_G(f)=\int _{-\infty }^{+\infty }Z_H(x_N)\,e^{\beta f x_N}\mathrm {d}x_N. \end{aligned}$$

(B15)

Thus, one can write $Z_H(x_N)$, using the change of variable $f\rightarrow - i \omega /\beta$, as an inverse Fourier transform in the complex plane

$$\begin{aligned} Z_H(x_N)=\frac{1}{2\pi }\int _{-\infty }^{+\infty }Z_G\left( -\frac{i\omega }{\beta }\right) \,e^{ i\omega x_N } \mathrm {d}\omega . \end{aligned}$$

(B16)

To simplify the notation and perform the calculation, we define

$$\begin{aligned} \delta _{1,2}=e^{-{\tilde{\beta }}}\,{\hat{\lambda }}_{1,2}=\frac{1}{2} \left( c_+ + c_- \pm \sqrt{ \Delta } \right) . \end{aligned}$$

(B17)

Thus, the Gibbs partition function can be written as

$$\begin{aligned} Z_G({\tilde{f}}) =\frac{(\ell e^{{\tilde{\beta }}})^N}{2\cosh {\tilde{\beta }}} \left[ \delta _1^N \left( 1+e^{-2{\tilde{\beta }}}\frac{\delta _1+\delta _2}{\delta _1-\delta _2} \right) +\delta _2^N\left( 1-e^{-2{\tilde{\beta }}}\frac{\delta _1+\delta _2}{\delta _1-\delta _2}\right) \right] . \end{aligned}$$

(B18)

By using the Newton binomial rule to expand the powers $\delta _1^N$ and $\delta _2^N$, we find

$$\begin{aligned} \begin{aligned} Z_G({\tilde{f}}) =&\left( \frac{\ell e^{{\tilde{\beta }}}}{2}\right) ^N\frac{1}{2\cosh {\tilde{\beta }}} \left\{ \sum _{k=0}^N \left( {\begin{array}{c}N\\ k\end{array}}\right) (c_++c_-)^{N-k}\Delta ^{\frac{k}{2}} \left[ 1+(-1)^k \right] \right. \\&\left. +\sum _{k=0}^N\left( {\begin{array}{c}N\\ k\end{array}}\right) (c_++c_-)^{N-k+1}e^{-2{\tilde{\beta }}}\Delta ^{\frac{k-1}{2}} \left[ 1-(-1)^k\right] \right\} . \end{aligned} \end{aligned}$$

(B19)

Here, we can separate the even and odds terms as follows

$$\begin{aligned} \sum _{k=0}^Na_k=\sum _{k=0}^{\left[ \frac{N}{2}\right] }a_{2k}+\sum _{k=0}^{\left[ \frac{N-1}{2}\right] }a_{2k+1}, \end{aligned}$$

(B20)

where the square brackets in the sums stand for the floor function defined as $\left[ x\right] =\max \left\{ n\in {\mathbb {Z}}\vert n\le x \right\}$. Then, we get

$$\begin{aligned} \begin{array}{ll} Z_{G}({\tilde{f}})= &{} \frac{(\ell e^{{\tilde{\beta }}})^N}{2^{N}\cosh {\tilde{\beta }}}\left\{ \sum _{k=0}^{\left[ \frac{N}{2}\right] }\left( {\begin{array}{c}N\\ 2k\end{array}}\right) (c_++c_-)^{N-2k}\Delta ^k\right. \\ &{}+\left. e^{-2{\tilde{\beta }}}\sum _{k=0}^{\left[ \frac{N-1}{2}\right] }\left( {\begin{array}{c}N\\ 2k+1\end{array}}\right) (c_++c_-)^{N-2k}\Delta ^{k}\right\} . \end{array} \end{aligned}$$

(B21)

We can further develop the powers $\Delta ^k$ and $(c_++c_-)^{N-2j}$ through the Newton binomial rule, obtaining

$$\begin{aligned} \begin{array}{ll} Z_G({\tilde{f}}) =&{}\frac{(\ell e^{{\tilde{\beta }}})^N}{2^{N}\cosh {\tilde{\beta }}}\Biggl \{\sum _{k=0}^{\left[ \frac{N}{2}\right] }\sum _{j=0}^k\sum _{s=0}^{N-2j}\left( {\begin{array}{c}N\\ 2k\end{array}}\right) \left( {\begin{array}{c}k\\ j\end{array}}\right) \left( {\begin{array}{c}N-2j\\ s\end{array}}\right) \Biggr .\\ &{}\times c_-^{N-j-s}c_+^{j+s}(-1)^j4^{j}\left( 1-e^{-4{\tilde{\beta }}}\right) ^j \\ &{}+e^{-2{\tilde{\beta }}}\sum _{k=0}^{\left[ \frac{N-1}{2}\right] }\sum _{j=0}^k\sum _{s=0}^{N-2j}\left( {\begin{array}{c}N\\ 2k+1\end{array}}\right) \left( {\begin{array}{c}k\\ j\end{array}}\right) \left( {\begin{array}{c}N-2j\\ s\end{array}}\right) \\ &{}\Biggl .\times c_-^{N-j-s}c_+^{j+s}(-1)^j4^{j}\left( 1-e^{-4{\tilde{\beta }}}\right) ^j\Biggr \}. \end{array} \end{aligned}$$

(B22)

Finally, Eq. (B16) can be evaluated, yielding the canonical partition function within the Helmholtz ensemble given in Eqs. (14) and (15).

Appendix C: Asymptotic force–extension relation for the zipper model under isometric condition

Let us then introduce the average chain stretch ${\bar{\lambda }} = \frac{x_N}{\ell N}$, prescribed to the chain under isometric condition. By using Eqs. (16) and (51), we obtain the force–extension relation in the form

$$\begin{aligned} \langle {\tilde{f}} \rangle = \frac{ (F_0G_0+F_NG_N)\left( 1-e^{-2{\tilde{\beta }}}\right) +e^{-2{\tilde{\beta }}}\displaystyle \sum \nolimits _{\xi =0}^NF(\xi )G(\xi ) }{ (F_0+F_N)\left( 1-e^{-2{\tilde{\beta }}}\right) +e^{-2{\tilde{\beta }}}\displaystyle \sum \nolimits _{\xi =0}^NF(\xi ) }, \end{aligned}$$

(C23)

where we introduced the following functions:

$$\begin{aligned} \begin{aligned} G(\xi )=&{\hat{k}}\left( \frac{\xi }{N}\right) \frac{(N-\xi )({\bar{\lambda }}-1)+ \xi ({\bar{\lambda }}-\chi )}{N},\\ F(\xi )=&\frac{\exp \left\{ -{\tilde{\beta }}\left[ {\hat{k}}\left( \frac{\xi }{N}\right) \frac{[(N-\xi )({\bar{\lambda }}-1)+ \xi (\bar{\lambda }-\chi )]^2}{2\,N}+\Delta {\tilde{E}}\xi \right] \right\} }{\sqrt{N\alpha ^{\xi }/{\hat{k}}\left( \frac{\xi }{N}\right) }}, \end{aligned} \end{aligned}$$

(C24)

we defined the rescaled global stiffness of the system (with $0\le t\le 1$)

$$\begin{aligned} {\hat{k}}(t)=\left( \frac{1-t}{{\tilde{k}}}+\frac{t}{\alpha {\tilde{k}}}\right) ^{-1}, \end{aligned}$$

(C25)

and we used the compact notations $F_0=F(0),G_0=G(0),F_N=F(N),G_N=G(N)$. Let us now consider the behavior in the thermodynamical limit $N\rightarrow \infty$. By following the approach suggested in Refs. [65, 66] to obtain explicit analytical results, not only for the stress plateau but also for the stress peak, we may consider the Euler–MacLaurin (EM) approximation for a given function $\phi$

$$\begin{aligned} \sum _{\xi =0}^N\phi (\xi ) \simeq \int _0^N\phi (\xi )\mathrm {d}\xi +\frac{\phi (0)+\phi (N)}{2}, \end{aligned}$$

(C26)

where higher-order terms of the EM approximation would lead to more detailed, but analytically cumbersome results. Thus, from Eq. (C23) we get

$$\begin{aligned} \langle {\tilde{f}} \rangle = \frac{ (F_0G_0+F_NG_N)\left( 1-\frac{e^{-2{\tilde{\beta }}}}{2}\right) +e^{-2{\tilde{\beta }}}\int _{0}^N F(\xi )G(\xi )\mathrm {d}\xi }{ (F_0+F_N)\left( 1-\frac{e^{-2{\tilde{\beta }}}}{2}\right) +e^{-2{\tilde{\beta }}}\int _{0}^N F(\xi )\mathrm {d}\xi }. \end{aligned}$$

(C27)

To simplify the calculation, we may rewrite the integrals as

$$\begin{aligned} \int _0^NF(\xi )\mathrm {d}\xi =&\,\sqrt{N}\int _0^1e^{Ng(\eta )}\sqrt{{\hat{k}}(\eta )}\mathrm {d}\eta , \end{aligned}$$

(C28)

$$\begin{aligned} \int _0^NF(\xi )G(\xi )\mathrm {d}\xi =&\,\sqrt{N}\int _0^1e^{Ng(\eta )}\sqrt{{\hat{k}}(\eta )}f(\eta )\mathrm {d}\eta , \end{aligned}$$

(C29)

where we introduced the phase fraction $\eta =\xi /N$ and the auxiliary functions

$$\begin{aligned} g(\eta )=&\frac{\eta }{2}\log \frac{1}{\alpha }-{\tilde{\beta }}\left( \eta \Delta {\tilde{E}}+{\hat{k}}(\eta )\frac{[(1-\eta )({\bar{\lambda }}-1)+ \eta ({\bar{\lambda }}-\chi )]^2}{2}\right) , \end{aligned}$$

(C30)

$$\begin{aligned} f(\eta )=&{\hat{k}}(\eta )\left[ (1-\eta )({\bar{\lambda }}-1)+ \eta ({\bar{\lambda }}-\chi )\right] . \end{aligned}$$

(C31)

Now, we use the following general results [130, 131]. Let

$$\begin{aligned} I(x)=\int _a^b e^{\,x\,\,{\mathcal {G}}(t)}{\mathcal {F}}(t)\mathrm {d}t, \end{aligned}$$

(C32)

then, we have

1. Suppose ${\mathcal {F}}$ is bounded and continuous on $(a,\,b)$, and ${\mathcal {F}}(a){\mathcal {F}}(b)\ne 0$. Suppose also that ${\mathcal {G}}$ is strictly monotone and differentiable and that $\frac{{\mathcal {F}}(a)}{{\mathcal {G}}'(a)}$ and $\frac{{\mathcal {F}}(b)}{{\mathcal {G}}'(b)}$ both exist as finite reals, defined as limits if either endpoint is infinite. Assume also that the integral in Eq. (C32) exists for all $x>0$. Then, we have that

$$\begin{aligned} I(x)\underset{x\rightarrow \infty }{\sim }\frac{1}{x}\frac{{\mathcal {F}}(b)}{{\mathcal {G}}'(b)}e^{\,x\,\,{\mathcal {G}}(b)}-\frac{1}{x}\frac{{\mathcal {F}}(a)}{{\mathcal {G}}'(a)}e^{\,x\,\,{\mathcal {G}}(a)}. \end{aligned}$$

(C33)

2. Suppose ${\mathcal {F}}$ is bounded and continuous on $(a,\,b)$, that ${\mathcal {G}}$ has unique maximum at some c in the open interval $(a,\,b)$, ${\mathcal {G}}$ is differentiable in some neighborhood of c, ${\mathcal {G}}''(c)$ exists and is ${\mathcal {G}}''(c)<0$, and that ${\mathcal {F}}(c)\ne 0$. Then, we have that

$$\begin{aligned} I(x)\underset{x\rightarrow \infty }{\sim }\frac{\sqrt{2\pi }{\mathcal {F}}(c)e^{\,x\,\,{\mathcal {G}}(c)}}{\sqrt{-x{\mathcal {G}}''(c)}}. \end{aligned}$$

(C34)

We use the results in Eqs. (C33) and (C34) to analyze the behavior of the integrals defined in Eqs. (C28) and (C29). First, by using Eq. (C30), we obtain

$$\begin{aligned} \begin{aligned} \frac{\partial g}{\partial \eta }=&{\tilde{\beta }}{\hat{k}}(\eta )\left[ (1-\eta )({\bar{\lambda }}-1)+ \eta ({\bar{\lambda }}-\chi )\right] (\chi -1)-{\tilde{\beta }}\Delta {\tilde{E}}+\frac{1}{2}\log \frac{1}{\alpha }\\&+\frac{{\tilde{\beta }}[{\hat{k}}(\eta )]^2}{2}[(1-\eta )({\bar{\lambda }}-1)+ \eta (\bar{\lambda }-\chi )]^2\left( \frac{1}{\alpha {\tilde{k}}}-\frac{1}{{\tilde{k}}}\right) . \end{aligned} \end{aligned}$$

(C35)

Solving $\frac{\partial g}{\partial \eta }=0$, we obtain the solution $\eta _0$ that represents the stationary point to be computed. Thus, following the definition in Eq. (C31), we also introduce the value

$$\begin{aligned} {\tilde{f}}_0={\tilde{f}}(\eta _0)={\hat{k}}(\eta _0)\left[ (1-\eta _0)(\bar{\lambda }-1)+ \eta _0({\bar{\lambda }}-\chi )\right] . \end{aligned}$$

(C36)

We observe that the equation $\frac{\partial g}{\partial \eta }=0$ can be written in terms of ${\tilde{f}}_0$ as

$$\begin{aligned} \frac{1}{2}\log \left( \frac{1}{\alpha }\right) -{\tilde{\beta }}\left[ \Delta {\tilde{E}}-(\chi -1){\tilde{f}}_0-\left( \frac{1}{\alpha }-1\right) \frac{{\tilde{f}}_0^2}{2{\tilde{k}}}\right] =0, \end{aligned}$$

(C37)

that coincides with Eq. (29) (i.e., with $\delta =0$). Thus, we identify ${\tilde{f}}_0$ with the Maxwell force ${\tilde{f}}_M$. Now, we can use Eq. (C34) (stationary phase theorem) to simplify Eqs. (C28) and (C29) only if $0<\eta _0<1$ (the stationary point must be within the integration interval). For $\eta _0=0$, we have ${\tilde{f}}_0=(\bar{\lambda }-1){\tilde{k}}$, and for $\eta _0=1$, we have ${\tilde{f}}_0=(\bar{\lambda }-\chi )\alpha {\tilde{k}}$, as given by Eq. (C36). Hence, the interval $0<\eta _0<1$ is equivalent to $1+\frac{{\tilde{f}}_0}{{\tilde{k}}}<\bar{\lambda }<\chi +\frac{{\tilde{f}}_0}{\alpha {\tilde{k}}}$, which corresponds to the force plateau region between the two elastic branches. Therefore, only in the plateau interval we have a stationary point and we can approximate Eqs. (C28) and (C29) with Eq. (C34). The application of the stationary phase method is further justified by the relation $g''(\eta _0)<0$, simply proved by a direct evaluation

$$\begin{aligned} g''(\eta _0)= -{\tilde{\beta }}{\hat{k}}(\eta _0)\left[ \left( \chi +\frac{{\tilde{f}}_0}{\alpha {\tilde{k}}}\right) -\left( 1+\frac{{\tilde{f}}_0}{{\tilde{k}}}\right) \right] ^2<0. \end{aligned}$$

(C38)

Note that from now on, we use the notation ${\tilde{f}}_M$ for the quantity ${\tilde{f}}_0$, to be consistent with previous Sections. We obtain from Eq. (C27) the expression

$$\begin{aligned} \langle {\tilde{f}} \rangle = \frac{{\mathcal {C}}_N\left[ e^{\,\Delta g_0 N}({\bar{\lambda }}-1){\tilde{k}}^{\frac{1}{2}}+e^{\,\Delta g_1 N}(\bar{\lambda }-\chi )(\alpha {\tilde{k}})^{\frac{3}{2}}\right] +{\tilde{f}}_M}{{\mathcal {C}}_N\left[ \sqrt{{\tilde{k}}}\,\,e^{\,\Delta g_0 N}+\sqrt{\alpha {\tilde{k}}}\,\,e^{\,\Delta g_1 N}\right] +1}, \end{aligned}$$

(C39)

where

$$\begin{aligned} \Delta g_0=&g(0)-g(\eta _0)=-\frac{{\tilde{\beta }}{\tilde{k}}}{2}\left[ {\bar{\lambda }}-\left( 1+\frac{{\tilde{f}}_M}{{\tilde{k}}}\right) \right] ^2,\end{aligned}$$

(C40)

$$\begin{aligned} \Delta g_1=&g(1)-g(\eta _0)=-\frac{{\tilde{\beta }}\alpha {\tilde{k}}}{2}\left[ {\bar{\lambda }}-\left( \chi +\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}\right) \right] ^2, \end{aligned}$$

(C41)

$$\begin{aligned} {\mathcal {C}}_N=&\frac{2 \,e^{\,2{\tilde{\beta }}}-1}{2 }\sqrt{\frac{{\tilde{\beta }}}{2\pi \,N}}\left[ \chi +\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}-\left( 1+\frac{{\tilde{f}}_M}{{\tilde{k}}}\right) \right] . \end{aligned}$$

(C42)

This is the main result given in Eqs. (52)–(55) of the main text. Considering that $\Delta g(0)<0$ and $\Delta g(1)<0$ (since $g(\eta _0)$ is the maximum value of $g(\eta )$ in the interval $0<\eta <1$), we may prove that in the thermodynamic limit

$$\begin{aligned} \lim _{N\rightarrow \infty }\langle {\tilde{f}}\rangle = {\tilde{f}}_M, \end{aligned}$$

(C43)

which is valid for $1+\frac{{\tilde{f}}_M}{{\tilde{k}}}<\bar{\lambda }<\chi +\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}$, and meaning that the Maxwell force is the same for both the Helmholtz and Gibbs ensembles, proving their equivalence for the zipper model.

In order to conclude the analysis, we have to simplify Eq. (C27) also for the external regions, i.e., for ${\bar{\lambda }}<1+\frac{{\tilde{f}}_M}{{\tilde{k}}}$ and $\bar{\lambda }>\chi +\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}$. Hence, we suppose that the critical point $\eta _0$ is external to the interval (0, 1). We can have either $\eta _0<0$ or $\eta _0>1$. In these cases, the asymptotic behavior of the integrals in Eqs. (C28) and (C29) is described by Eq. (C33). Then, Eq. (C27) assumes the form

$$\begin{aligned} \langle {\tilde{f}}\rangle = \frac{e^{Ng(0)}{\tilde{k}}^{\frac{3}{2}}(\bar{\lambda }-1)+e^{Ng(1)}(\alpha {\tilde{k}})^{\frac{3}{2}}(\bar{\lambda }-\chi )+\frac{2}{2e^{2{\tilde{\beta }}}-1}{\mathcal {A}}}{e^{Ng(0)}\sqrt{{\tilde{k}}}+e^{Ng(1)}\sqrt{\alpha {\tilde{k}}}+\frac{2}{2e^{2{\tilde{\beta }}}-1}{\mathcal {B}}}, \end{aligned}$$

(C44)

where

$$\begin{aligned} {\mathcal {A}}=&\frac{(\alpha {\tilde{k}})^{\frac{3}{2}}}{g'(1)}e^{Ng(1)}({\bar{\lambda }}-\chi )-\frac{{\tilde{k}}^{\frac{3}{2}}}{g'(0)}e^{Ng(0)}({\bar{\lambda }}-1),\end{aligned}$$

(C45)

$$\begin{aligned} {\mathcal {B}}=&\frac{\sqrt{\alpha {\tilde{k}}}}{g'(1)}e^{Ng(1)}-\frac{\sqrt{{\tilde{k}}}}{g'(0)}e^{Ng(0)}. \end{aligned}$$

(C46)

Now, by using Eq. (C37), we determine the quantity $g(1)-g(0)$ that eventually reads

$$\begin{aligned} g(1)-g(0)=\frac{{\tilde{\beta }}{\tilde{k}}}{2}\left[ \left( \bar{\lambda }-1-\frac{{\tilde{f}}_M}{{\tilde{k}}}\right) ^2 -\alpha \left( \bar{\lambda }-\chi -\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}\right) ^2\right] . \end{aligned}$$

(C47)

Accordingly, in the limit of $N\rightarrow \infty$ we obtain

$$\begin{aligned} g(1)-g(0) \left\{ \begin{aligned}&>0 \,\,\,\,\,\, \text {if} \,\,\,\,\,\, {\bar{\lambda }}>\chi +\frac{{\tilde{f}}_M}{\alpha {\tilde{k}}}\Rightarrow \langle {\tilde{f}}\rangle =\alpha {\tilde{k}} ({\bar{\lambda }}-\chi ),\\&<0 \,\,\,\,\,\, \text {if} \,\,\,\,\,\, {\bar{\lambda }}<1+\frac{{\tilde{f}}_M}{{\tilde{k}}}\Rightarrow \langle {\tilde{f}}\rangle ={\tilde{k}} ({\bar{\lambda }}-1),\\ \end{aligned} \right. \end{aligned}$$

(C48)

a result representing the elastic branches in the external regions, corresponding to the fully folded (left) and fully unfolded (right) phases. This completes the proof of the equivalence of the ensembles in the thermodynamic limit for the zipper model.