Temperature controlled decohesion regimes of an elastic chain adhering to a fixed substrate by softening and breakable bonds

Micro-instabilities play a very important role in several artificial and biological mechanical systems. On one hand, concerning artificial systems, we can mention the peeling of a film from a substrate [1–7], the waves propagation in bistabile lattices [8–12], the energy harvesting through multistable chains [13–15], the plasticity and the hysteresis in phase transitions and martensitic transformations of solids [16–26] and the cracks and dislocations nucleation and propagation in materials and alloys [27–35]. On the other hand, micro-mechanical biological phenomena include the conformational transitions in polymeric and biopolymeric chains [36–49], the attached and detached states of fibrils in cell adhesion [50–55], the unzipping of macromolecular hairpins [56–61], the sarcomeres behavior in skeletal muscles [62–68] and the denaturation or degradation of nucleic acids, polypeptidic chains or other macromolecules of biological origin [69–78].

In all these physical situations we eventually find a multi-basin energy landscape and the state of the system can be in stable or meta-stable configurations, identified by the wells of the energy function. As a matter of fact, these systems are constituted by a large number of units characterized by two or more physical states. The transitions between these states or, equivalently, the exploration of the energy landscape, govern the macroscopic behavior of the whole system and, in particular, its static and dynamic features.

In a first class of systems, the intrinsic micro-instabilities may describe bistable units with transitions between one ground state and one metastable state (e.g. for the conformational folded-to-unfolded transitions in macromolecules or martensitic phase changes in metallic alloys). These types of transitions are typically reversible, which means that upon unloading the system follows the same path. Alternatively, in a second class of systems, investigated in this work, the micro-instabilities can model transitions between unbroken (attached) and broken (detached) states of breakable links (e.g. in the unzipping of hairpins, denaturation of macromolecules, fibrillar biological adhesion, peeling of films and cracks propagation). We underline that this kind of transitions can be reversible or irreversible depending on the considered physical system. For example, DNA denaturation is typically reversible whereas crack propagation is typically an irreversible process. In this work, as fully discussed below, since we do not consider cyclic loading, this difference is not crucial. On the other hand, because we refer to equilibrium statistical mechanics approach, we may fix our attention on reversible processes. Of course, the extension to irreversibility could be schematically described by introducing internal variables such as in [79] and assuming that each broken element cannot participate to the subsequent process.

The transitions between the states in all previous systems may be strongly influenced by thermal fluctuations, which can modify the probability of being in a given state or the passage rate between the neighboring energy wells. In particular, entropic energy terms can be relevant in transition and fracture phenomena of soft materials, such as biological and rubberlike materials [80, 81], or in the case of very low dimensional scales systems such as in shape memory nanowires [82]. In this context, the correct framework is the classical statistical mechanics. In particular, the systems exhibiting switching mechanisms between different energy basins can be studied by means of the spin variables approach.

The first theoretical approaches based on this method have been developed to model the biomechanical response of skeletal muscles [62, 63]. Then, this technique has been generalized to study different multi-stable systems [65–68], macromolecular chains [83–89] and phase transformations in nanowires [82]. This approach is based on the introduction of a series of discrete variables (similar to the spins used to deal with magnetic systems), which are able to identify the state associated with a given system unit. In other words, depending on the value of these discrete variables, the system transits between distinct energy wells possibly characterized by different position, shape, stiffness and depth. The introduction of the spin variables frequently simplifies the calculation of the partition function and the analysis of the corresponding averaged thermodynamic quantities. This technique has been largely exploited to investigate bistable systems with transitions between ground and metastable states, with important applications to nanomechanics [82–89]. Similarly, this approach has been considered for breakable materials with the spin variable distinguishing broken and unbroken units to study debonding processes in biological materials [27–29, 80]. In this paper, we introduce an important generalization of the model presented in reference [80], by including a softening mechanism in the breakable units of the system. More specifically, we study a one-dimensional lattice of masses linked by harmonic springs and connected to a substrate by breakable links subjected to a softening mechanism (see figure 1). This geometry has been previously introduced to describe a wide range of phenomena such as peeling of tapes, adhesion of geckos and denaturation of macromolecules [4, 5]. However, the thermal fluctuations and the softening mechanism have been neglected in these works, while in the physical systems recalled above they can play a crucial role that is analyzed in detail in the following development.

Zoom In Zoom Out Reset image size

From the mechanical point of view, we consider two different kinds of loading, which induce the decohesion of the system from the substrate. The decohesion process can be induced by either prescribing a given extension of the last element of the chain (figure 1(a)), or by applying an external force to the last unit of the chain (figure 1(b)). On one hand, the first isometric condition corresponds to the Helmholtz ensemble of statistical mechanics and can be generated by hard devices. On the other hand, the isotensional condition corresponds to the Gibbs ensemble and can be generated by soft devices. Indeed, these two boundary conditions can be deduced as limiting cases of real loading experiments, when the stiffness of the device is large (hard device) or is negligible (soft device) as compared with the loaded system one, respectively [88, 89]. From a theoretical point of view, an intriguing problem concerns the equivalence of the two ensembles in the thermodynamic limit (i.e. for very large systems) [90–96]. Interestingly, we will prove, in our case, the non-equivalence of the ensembles in the thermodynamic limit. Moreover, we are able to study the force necessary to detach the chain from the substrate as function of the temperature and of the external mechanical action applied to the system. In this context, we obtain a critical behavior described by specific phase transitions. Also the evolution of the number of intact, softened and broken elements is thoroughly investigated in both isometric and isotensional conditions. To conclude, the aim of this work is to fully analyze the cohesion–decohesion process with a softening mechanism in both the Helmholtz and Gibbs ensembles, thus providing a complete picture of the effect of the temperature and loading type on this prototypical physical system.

The paper is organized as it follows. In section 2 we define the problem and we simplify the spin variable methods through the zipper assumption. Then, in section 3, we introduce the system behavior under isometric condition, and in section 4, we analyze its thermodynamic limit. Similarly, in section 5, we introduce the system behavior under isotensional condition, and in section 6, we analyze its thermodynamic limit. The conclusions (section 7) and two mathematical appendices close the paper.

Our prototypical system is represented in figure 1. The horizontal springs of the lattice are purely harmonic, characterized by the elastic constant k (figure 1(c)), with a potential energy

$$\begin{equation}W\left({y}_{i+1}-{y}_{i}\right)=\frac{1}{2}k{\left({y}_{i+1}-{y}_{i}\right)}^{2}.\end{equation} \tag{ 1 }$$

On the other hand, the vertical ones can be in three different states, depending on their extension y_i (figure 1(d)). When |y_i | < y_p , they are intact (elastic constant h_e ), when y_p < |y_i | < y_b they are softened (elastic constant h_p < h_e ), and when |y_i | > y_b they are broken (they do not support forces). This scheme represents the softening mechanism of the breakable bonds. We describe such a behavior in a spin formalism by introducing a discrete variable s_i associated to each vertical spring. Hence, we can write the potential energy of the breakable bonds in the form

$$\begin{align}U({y}_{i},{s}_{i})& =\frac{{s}_{i}(1+{s}_{i})}{2}\frac{1}{2}{h}_{e}{y}_{i}^{2}+(1+{s}_{i})(1-{s}_{i})\left(\frac{1}{2}{h}_{p}{y}_{i}^{2}+{\Delta}E\right)\\ & \quad -\frac{{s}_{i}(1-{s}_{i})}{2}\left(\frac{1}{2}{h}_{p}{y}_{b}^{2}+{\Delta}E\right),\end{align} \tag{ 2 }$$

where s_i = +1 corresponds to the intact state, s_i = 0 corresponds to the softened state, and s_i = −1 corresponds to the broken state (i = 1, ..., N). With this assumptions we have a phase space composed on the N continuous variables ${y}_{i}\in \mathbb{R}$ and the N discrete variables ${s}_{i}\in \left\{+1,0,-1\right\}$. Therefore, when we calculate the partition function, we have to integrate over all the continuous variables and to sum over all the discrete ones. So doing, the switching of the variable s_i and their statistics at thermodynamic equilibrium are directly controlled by the statistical ensemble (Helmholtz or Gibbs in our case) imposed to the system.

The assumption of considering independent spins for all the breakable bonds is the most rigorous approach to analyze the system behavior. However, on this assumption, the statistical mechanics of the system under investigation cannot be analytically developed and also its numerical implementation is quite expensive. Indeed, the calculation of the partition function is rather prohibitive because of the sum over all the possible spins combinations. Nevertheless, since we are studying the cohesion–decohesion process under an external mechanical action applied to one end point of the system, suppose the right one, we can simplify the model by assuming to have N − η broken elements on the right region of the chain, η − ξ softened elements on the central region of the chain and ξ intact or unbroken elements on the left region of the chain (see figure 1). In other words we suppose to have two moving interfaces or domain walls between the three regions of the chain. This hypothesis strongly reduces the mathematical complexity of the problem and is reasonable if we work at sufficiently low temperature and sufficiently large applied mechanical load (either extension or force). Indeed, the set of the three-state spin variables s_i is substituted by the two variables η and ξ, assuming values in the set $\left\{0,1,2,\dots ,N\right\}$. In this sense, η and ξ can be viewed as multi-valued spin variables. Moreover, this hypothesis is similar to the one adopted in the so-called zipper model, largely used to describe the helix-coil transitions in proteins, the gel–sol transition of thermo-reversible gels, and the melting or denaturation of DNA [97–100]. By implementing this simplified scheme, we deduce a fully analytical solution of both Helmholtz and Gibbs boundary problems for an arbitrary number N of elements of the chain, and also in the thermodynamic limit (N → ∞).

One of the most important result concerns the analytic expressions of the temperature dependent debonding force (or system strength) in both cases of isotensional and isometric loading. This temperature dependent behavior can be interpreted by observing that the thermal fluctuations may foster the decohesion, allowing the escape from the energy well shown in figure 1(d). More rigorously, we prove that this behavior can be explained in terms of a phase transition occurring at a given critical temperature. In particular, it means that the system can be completely detached from the substrate for supercritical temperatures even without any external mechanical action.

This behavior has been observed in several polymeric systems subjected to a force [101–103]. In particular, several important results have been obtained for the pulling processes of adsorbed polymers on a surface. They were found by means of directed or partially directed walk models of lattice polymers adsorbed at a surface under the influence of an applied force [104, 105]. In this context, a temperature-dependent critical force has been derived for semi-flexible polymers [106], polymers subject to an arbitrarily oriented force [107], heterogeneous adsorption surfaces [108], striped adsorption surfaces [109], and self-avoiding chains [110–112]. From a theoretical point of view, the main difference between these approaches and ours is that they use lattice polymers techniques whereas we adopt continuous geometric variables. Interestingly, the lattice methods can also describe bubbles with different physical states along the chain, without needing the zipper simplifying assumption. In our case neglecting such simplifying assumption would lead to a numerical treatment of the subject and we preferred the analytical clearness of the results. Importantly, in our continuous model we also introduced the softening mechanism, which is the main topic of this work and represents, as we show in the following, an important ingredient to distinguish different observed regimes in the decohesion behavior.

In our system, although the decohesion force for subcritical temperatures is the same for both statistical ensembles, they are not equivalent in the thermodynamic limit since the force–extension curves are different under isometric and isotensional conditions. Importantly, the softening mechanism in the breakable elements generates a strength-temperature curve composed of two branches connected, with continuity, at a given temperature T₀. From the physical point of view, this is the temperature at which all the elements are at least softened. The interesting point is that this transition has been experimentally observed in the strength behavior of some systems including sapphire whiskers [30] and a number of high-entropy and medium-entropy alloys [32, 35]. In these cases, the softening mechanism is originated by the emergence of a population of dislocations beyond a certain threshold of deformation.

Another fundamental effect induced by the bonds softening is the possibility of describing a temperature dependent process zone that, for enough high temperatures, anticipates the propagation front. As a result we may obtain a transition between different regimes, experimentally observed for example in polymeric materials [113, 114].

Summarizing, the proposed system is particularly important for the following reasons: (i) once previous crucial assumptions are considered, the model can be analytically solved in both isometric and isotensional statistical ensembles; (ii) the solution shows a phase transition at a critical temperature that can be calculated in closed form; (iii) the studied system manifestly shows the ensembles non-equivalence in the thermodynamic limit, which is an unusual feature in statistical mechanics; (iv) finally, the model describe the existence of a two-branch curve for the strength-temperature behavior, observed in real materials, as a result of a softening anticipating breaking effects.

We consider the cohesion–decohesion process in the system represented in figure 1(a), where the detachment is generated by imposing the extension y_N+1 of the last element of the chain. This condition corresponds to the Helmholtz ensemble of statistical mechanics and it is obtained in the case of loading with hard devices (with a very large, ideally infinity, intrinsic elastic constant). We identify the longitudinal springs with the potential energy $W(z)=\frac{1}{2}k{z}^{2}$ for any $z\in \mathbb{R}$, and the transverse springs with the potential energy $U(y)=\frac{1}{2}{h}_{e}{y}^{2}$ if |y| ⩽ y_p (intact), $U(y)=\frac{1}{2}{h}_{p}{y}^{2}+{\Delta}E$ if y_p < |y| ⩽ y_b (softened), and $U(y)=\frac{1}{2}{h}_{p}{y}_{b}^{2}+{\Delta}E$ if |y| > y_b (broken). The values ±y_p of the extension correspond to the softening points A and B of the breakable spring, where the elastic constant switch from h_e to h_p < h_e . The force jump at y = ±y_p is given by $\sqrt{2{\Delta}E({h}_{e}-{h}_{p})}$ (see figure 1(d)). We simply calculate that ${y}_{p}=\sqrt{2{\Delta}E/({h}_{e}-{h}_{p})}$, and we must always impose the inequality y_b > y_p , where ±y_b are the extensions at the breaking points C and D of the transverse elements. As anticipated, we assume that a first group of ξ elements of the chain are intact (i = 1...ξ), a second group of η − ξ elements are softened (i = ξ + 1...η), and a third group of N − η elements are broken (i = η + 1...N). These premises allow us to write the total potential energy of the system in the following form

$$\begin{align}{\Phi}& =\sum\limits _{i=0}^{N}W\left({y}_{i+1}-{y}_{i}\right)+\sum\limits _{i=1}^{N}U\left({y}_{i}\right)\\ & =\sum\limits _{i=0}^{N}\frac{1}{2}k{\left({y}_{i+1}-{y}_{i}\right)}^{2}+\sum\limits _{i=1}^{\xi }\frac{1}{2}{h}_{e}{y}_{i}^{2}\\ & \quad +\sum\limits _{i=1+\xi }^{\eta }\left[\frac{1}{2}{h}_{p}{y}_{i}^{2}+{\Delta}E\right]+\sum\limits _{i=1+\eta }^{N}\left[\frac{1}{2}{h}_{p}{y}_{b}^{2}+{\Delta}E\right],\end{align} \tag{ 3 }$$

where the applied extension y_N+1 is considered as a parameter and we fix y₀ = 0. The variables belonging to the phase space are the extensions y₁, ..., y_N and the two configurational numbers ξ and η, characterizing the state of the breakable elements. For further convenience, the energy function Φ can be further arranged as follows

$$\begin{align}{\Phi}\left(\boldsymbol{y},\xi ,\eta ;{y}_{N+1}\right)& =\frac{1}{2}k\mathcal{A}(\xi ,\eta )\boldsymbol{y}\cdot \boldsymbol{y}-k{y}_{N+1}\boldsymbol{v}\cdot \boldsymbol{y}\\ & \quad +\frac{1}{2}k{y}_{N+1}^{2}+(N-\xi ){\Delta}E+(N-\eta )\frac{1}{2}{h}_{p}{y}_{b}^{2},\end{align} \tag{ 4 }$$

where we defined the vector $\boldsymbol{y}=({y}_{1},\dots ,{y}_{N})\in {\mathbb{R}}^{N}$, with each component representing the displacement of an element of the chain, the constant vector $\boldsymbol{v}=(0,0,0,\dots ,0,1)\in {\mathbb{R}}^{N}$ and the tridiagonal (symmetric and positive definite) matrix

$$\begin{equation}\mathcal{A}(\xi ,\eta )=\left[\begin{matrix}\hfill {a}_{1}\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill \\ \hfill -1\hfill & \hfill {a}_{2}\hfill & \hfill -1\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {a}_{N-1}\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill {a}_{N}\hfill \end{matrix}\right]\in {\mathcal{M}}_{N,N}(\mathbb{R}).\end{equation} \tag{ 5 }$$

We observe that ξ and η must fulfill the condition 0 ⩽ ξ ⩽ η ⩽ N. Moreover, this matrix has all the subdiagonal and superdiagonal elements equal to −1 and the diagonal elements defined as follows

$$\begin{equation}\begin{cases}{a}_{1}=\cdots ={a}_{\xi }=2+\alpha ,\\ {a}_{1+\xi }=\cdots ={a}_{\eta }=2+\beta ,\\ {a}_{1+\eta }=\cdots ={a}_{N}=2,\end{cases}\end{equation} \tag{ 6 }$$

where we introduced the rescaled elastic constants

$$\begin{equation}\alpha ={h}_{e}/k,\hspace{25.0pt}\beta ={h}_{p}/k.\end{equation} \tag{ 7 }$$

The algebraic properties of the matrix $\mathcal{A}(\xi ,\eta )$ are studied in appendix  A, where we evaluate the quantities needed to obtain the temperature dependent relations among the observables during the cohesion–decohesion process.

The equilibrium statistical mechanics of a system in contact with a reservoir at temperature T is described by the canonical distribution. In this context, the expectation values of physical observables can be obtained by evaluating the partition function that, within the Helmholtz ensemble, can be written as

$$\begin{equation}{Z}_{\mathrm{H}}\left({y}_{N+1}\right)=\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }{\int }_{{\mathbb{R}}^{N}}{\mathrm{e}}^{-\frac{{\Phi}\left(\boldsymbol{y},\xi ,\eta ;{y}_{N+1}\right)}{{K}_{B}T}}\enspace \mathrm{d}\boldsymbol{y},\end{equation} \tag{ 8 }$$

where we have to integrate the continuous variables represented by the vector y , and to sum over the discrete variables ξ and η, describing the state (intact, softened or broken) of the elements of the chain (0 ⩽ ξ ⩽ η ⩽ N). Using equation (4), Z_H can be evaluated using the property of Gaussian integrals

$$\begin{equation}{\int }_{{\mathbb{R}}^{N}}{\mathrm{e}}^{-\frac{1}{2}\mathcal{M}\boldsymbol{y}\cdot \boldsymbol{y}}\enspace {\mathrm{e}}^{\boldsymbol{b}\cdot \boldsymbol{y}}\enspace \mathrm{d}\boldsymbol{y}=\sqrt{\frac{{(2\pi )}^{N}}{\mathrm{det}\enspace \mathcal{M}}}\enspace {\mathrm{e}}^{\frac{1}{2}{\mathcal{M}}^{-1}\boldsymbol{b}\cdot \boldsymbol{b}},\end{equation} \tag{ 9 }$$

holding for any symmetric and positive definite matrix $\mathcal{M}$. In particular, we can define $\mathcal{M}=\frac{k}{{K}_{B}T}\mathcal{A}$ and $\boldsymbol{b}=\frac{k{y}_{N+1}}{{K}_{B}T}\boldsymbol{v}$ to obtain

$$\begin{equation}{Z}_{\mathrm{H}}\left({y}_{N+1}\right)={\left(\frac{2\pi {K}_{B}T}{k}\right)}^{N/2}\enspace {\mathrm{e}}^{-N\frac{{\Delta}E}{{K}_{B}T}}\enspace {\mathrm{e}}^{-N\frac{k\beta {y}_{b}^{2}}{2{K}_{B}T}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1}),\end{equation} \tag{ 10 }$$

where

$$\begin{equation}{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})=\frac{{\text{e}}^{\xi \frac{{\Delta}E}{{K}_{B}T}}\enspace {\text{e}}^{\eta \frac{k\beta {y}_{b}^{2}}{2{K}_{B}T}}}{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}}\enspace {\mathrm{e}}^{-\frac{k{y}_{N+1}^{2}}{2{K}_{B}T}\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}},\end{equation} \tag{ 11 }$$

and ${\mathcal{A}}_{NN}^{-1}(\xi ,\eta )={\mathcal{A}}^{-1}(\xi ,\eta )\boldsymbol{v}\cdot \boldsymbol{v}$. We can thus determine the average force associated to the vertical extension y_N+1 of the last element of the chain as [90]

$$\begin{equation}\left\langle f\right\rangle =-{K}_{B}T\frac{\partial \mathrm{log}\enspace {Z}_{\mathrm{H}}}{\partial {y}_{N+1}}=-\frac{{K}_{B}T}{{Z}_{\mathrm{H}}}\frac{\partial {Z}_{\mathrm{H}}}{\partial {y}_{N+1}}.\end{equation} \tag{ 12 }$$

We obtain

$$\begin{equation}\left\langle f\right\rangle =\frac{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})}{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})}k{y}_{N+1}.\end{equation} \tag{ 13 }$$

This expression represents the force–extension relation for the system within the Helmholtz ensemble, or equivalently, under isometric condition. It is also important to calculate the average value of the number of broken elements $\left\langle N-\eta \right\rangle$ and the average number of softened elements $\left\langle \eta -\xi \right\rangle$. These quantities can be directly evaluated through the expressions

$$\begin{align}\left\langle N-\eta \right\rangle & =\frac{1}{{Z}_{\mathrm{H}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }{\int }_{{\mathbb{R}}^{N}}(N-\eta ){\mathrm{e}}^{-\frac{{\Phi}\left(\boldsymbol{y},\xi ,\eta ;{y}_{N+1}\right)}{{K}_{B}T}}\enspace \mathrm{d}\boldsymbol{y}\\ & =\frac{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }(N-\eta ){{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})}{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})},\end{align} \tag{ 14 }$$

$$\begin{align}\left\langle \eta -\xi \right\rangle & =\frac{1}{{Z}_{\mathrm{H}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }{\int }_{{\mathbb{R}}^{N}}(\eta -\xi ){\mathrm{e}}^{-\frac{{\Phi}\left(\boldsymbol{y},\xi ,\eta ;{y}_{N+1}\right)}{{K}_{B}T}}\enspace \mathrm{d}\boldsymbol{y}\\ & =\frac{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }(\eta -\xi ){{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})}{{\sum }_{\eta =0}^{N}{\sum }_{\xi =0}^{\eta }{{\Gamma}}_{\xi ,\eta }(T,{y}_{N+1})}.\end{align} \tag{ 15 }$$

In order to simplify the analysis it is useful to rescale the applied extension y_N+1 with respect to y_b (i.e. the extension corresponding to the breaking of the link) and define the non-dimensional extension

$$\begin{equation}\mathcal{Y}={y}_{N+1}/{y}_{b}.\end{equation} \tag{ 16 }$$

To simplify the formulas previously derived, we can consider the rescaled energies

$$\begin{equation}\delta =\frac{k\enspace {y}_{b}^{2}}{2{K}_{B}T},\hspace{25.0pt}\varphi =\frac{{\Delta}E}{{K}_{B}T}.\end{equation} \tag{ 17 }$$

Thus, introducing equations (16) and (17) into the expression of $\left\langle f\right\rangle$, $\left\langle N-\eta \right\rangle$ and $\left\langle \eta -\xi \right\rangle$, we obtain

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\frac{\left\langle f\right\rangle }{k{y}_{b}}=\frac{\mathcal{Y}}{{\mathcal{Z}}_{\mathrm{H}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )}{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}}\enspace {\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}},\end{equation} \tag{ 18 }$$

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{1}{{\mathcal{Z}}_{\mathrm{H}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{N-\eta }{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}}\enspace {\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}},\end{equation} \tag{ 19 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{1}{{\mathcal{Z}}_{\mathrm{H}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{\eta -\xi }{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}}\enspace {\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}},\end{equation} \tag{ 20 }$$

where we have used equations (10) and (11) to define the (rescaled) partition function

$$\begin{equation}{\mathcal{Z}}_{\mathrm{H}}(\mathcal{Y})=\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{{\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }}{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}}\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}}.\end{equation} \tag{ 21 }$$

Using the properties of $\mathcal{A}(\xi ,\eta )$ and the relations in equations (A.18) and (A.19), it is possible to find explicit expressions for $\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )$ and ${\mathcal{A}}_{NN}^{-1}(\xi ,\eta )$ so to obtain the behavior of the relevant physical quantities of the system.

Previous results are represented in figure 2, where we study the dependence of the relevant quantities on temperature and size of the system. In particular, we plotted $\left\langle \mathcal{F}\right\rangle$, $\left\langle \xi \right\rangle$, $\left\langle \eta -\xi \right\rangle$ and $\left\langle N-\eta \right\rangle$ versus the applied rescaled extension $\mathcal{Y}$ (changing the temperature) for two different cases: the first one concerns a short chain with N = 5 (panel (a)), and the second one a larger chain with N = 50 (panel (b)). For N = 5 (panel (a)), and for low values of the temperature, we easily identify the partial and complete rupture occurrences (see the peaks and/or steps in the curves of $\left\langle \mathcal{F}\right\rangle$, $\left\langle \xi \right\rangle$, $\left\langle \eta -\xi \right\rangle$ and $\left\langle N-\eta \right\rangle$). However, for larger values of the temperature, the curves are smoother because of the increasing thermal fluctuations. We also observe that the detachment of the elements of the chain occurs in a progressive way in response to the increasing extension $\mathcal{Y}$. This is a specific feature of the Helmholtz ensemble, typically observed in the folding/unfolding process of bistable macromolecules stretched by hard devices [83–86, 88, 89]. For N = 50 (panel (b)), we can identify two important properties of the detachment process in large systems: first, increasing $\mathcal{Y}$ the curve approximately exhibits a temperature dependent plateau; second, during the decohesion process the number of softened elements is temperature dependent and approximately constant for small values of temperature in a certain range of $\mathcal{Y}$. Although these features can be observed in the panel (b) of figure 2 only for low values of the temperature, we will deduce in the thermodynamic limit (i.e. when N → ∞) that this behavior can be explained in terms of a phase transition (see section 4). The existence of a constant finite softened domain shows the possibility of analytically describing the existence of a process zone regulating the advancing of the decohesion front as observed in cohesive fracture phenomena [113]. Interestingly, the size of the process zone increases with temperature, so that in accordance with the experimental behavior the decohesion is 'fragile' at low temperature and 'ductile' at higher values of the temperature [114].

Zoom In Zoom Out Reset image size

We study here the thermodynamic limit within the Helmholtz ensemble, defined by a very large number of elements in the chain (macroscopic limit valid for large systems). To this aim we have to study the behavior of equations (18)–(21) for N → ∞. In order to perform the analysis in this limit, we formally substitute the sums appearing in these expressions by integrals. In particular, we will use the approximation ${\sum }_{\chi =0}^{n}\phi (\chi )\simeq {\int }_{0}^{n}\phi (\chi )\mathrm{d}\chi$ for an generic function ϕ(χ). A better approximation could be adopted by using the Euler–McLaurin formula [115]. However, this approach, discussed in reference [80], is not necessary for the purposes of the present analysis. Moreover, due to the fact that we are considering the case of N large, it is possible to use the simplified expressions of $\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )$ and ${\mathcal{A}}_{NN}^{-1}(\xi ,\eta )$, obtained in equations (A.22) and (A.23).

From equations (18) and (21) we get the following explicit result for the force–extension relation

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\mathcal{Y}\frac{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{-3/2}\enspace \mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\mathrm{d}\xi }{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{-1/2}\enspace \mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\mathrm{d}\xi },\end{equation} \tag{ 22 }$$

where

$$\begin{equation}{\alpha }_{0}=\frac{2+\alpha +\sqrt{{\alpha }^{2}+4\alpha }}{2},\hspace{25.0pt}{\beta }_{0}=\frac{2+\beta +\sqrt{{\beta }^{2}+4\beta }}{2}.\end{equation} \tag{ 23 }$$

To better understand the behavior of the force–extension response for large values of N, we firstly perform the integration on the variable ξ, and then apply the change of variable N − η + β₀/(β₀ − 1) = s, delivering the following result

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\mathcal{Y}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{3/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}}.\end{equation} \tag{ 24 }$$

The behavior of the previous expression depends on the sign of the quantity $\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}$, appearing in the exponential term within the square brackets in both numerator and denominator. From equation (17) we have

$$\begin{equation}{\Theta}=\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\gtrless 0\enspace {\Rightarrow}\enspace T\lessgtr \frac{2{\Delta}E}{{K}_{B}\enspace \mathrm{log}\enspace \frac{{\alpha }_{0}}{{\beta }_{0}}}{:=}{T}_{0},\end{equation} \tag{ 25 }$$

where we have defined the temperature T₀, corresponding to a transition in the system behavior (we will prove that all elements are at least softened for T > T₀). If T < T₀ and we consider large values of N (thermodynamic limit), the exponential term in the square brackets of equation (24) is dominant. On the other hand, if T > T₀ the exponential term in the square brackets of equation (24) is negligible for large values of N. Thus, for N → ∞ we obtain

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\mathcal{Y}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{3/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\quad \text{for}\enspace T< {T}_{0}\end{equation} \tag{ 26 }$$

and

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\mathcal{Y}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{3/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\quad \text{for}\enspace T > {T}_{0}.\end{equation} \tag{ 27 }$$

In order to ensure the convergence of the integrals in equations (26) and (27) it is necessary that the following conditions hold

$$\begin{equation}\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi > 0\quad \text{if}\enspace T< {T}_{0},\end{equation} \tag{ 28 }$$

$$\begin{equation}\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0} > 0\quad \text{if}\enspace T > {T}_{0}.\end{equation} \tag{ 29 }$$

These conditions have a physical interpretation that will be discussed in the following.

The integrals appearing in equations (26) and (27) can be evaluated in terms of error functions as discussed in appendix  B. To simplify the form of the results, we define

$$\begin{equation}{\lambda }_{1}=\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi ,\hspace{25.0pt}{\lambda }_{2}=\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0},\hspace{25.0pt}c=\frac{{\beta }_{0}}{{\beta }_{0}-1}.\end{equation} \tag{ 30 }$$

We thus obtain

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\sqrt{\frac{{\lambda }_{1}}{\delta }}\left[\frac{2}{1+{\mathrm{e}}^{4\mathcal{Y}\sqrt{{\lambda }_{1}\delta }}\frac{1-\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}+\sqrt{{\lambda }_{1}c}\right)}{1+\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}-\sqrt{{\lambda }_{1}c}\right)}}-1\right]\quad \text{for}\enspace T< {T}_{0},\end{equation} \tag{ 31 }$$

$$\begin{equation}\left\langle \mathcal{F}\right\rangle =\sqrt{\frac{{\lambda }_{2}}{\delta }}\left[\frac{2}{1+{\mathrm{e}}^{4\mathcal{Y}\sqrt{{\lambda }_{2}\delta }}\frac{1-\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}+\sqrt{{\lambda }_{2}c}\right)}{1+\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}-\sqrt{{\lambda }_{2}c}\right)}}-1\right]\quad \text{for}\enspace T > {T}_{0}.\end{equation} \tag{ 32 }$$

It is possible to see that the quantities in the brackets in equations (31) and (32) converge to 1 when $\mathcal{Y}\to \infty$. Thus, the expectation value of the force $\left\langle \mathcal{F}\right\rangle$ is characterized by an asymptotic force ${\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}$ for large values of $\mathcal{Y}$ (coherently with figure 2(b)). In particular, we get the following asymptotic values

$$\begin{equation}{\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}=\sqrt{\frac{{\lambda }_{1}}{\delta }}=\sqrt{\frac{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi }{\delta }}=\sqrt{\frac{2{\Delta}E}{k{y}_{b}^{2}}+\frac{{h}_{p}}{k}}\sqrt{1-\frac{T}{{T}_{b}}}\quad \text{for}\enspace T< {T}_{0},\end{equation} \tag{ 33 }$$

where

$$\begin{equation}{T}_{b}=\frac{2{\Delta}E+{h}_{p}{y}_{b}^{2}}{{K}_{B}\enspace \mathrm{log}\enspace {\alpha }_{0}}=\frac{2{\Delta}E+{h}_{p}{y}_{b}^{2}}{{K}_{B}\enspace \mathrm{log}\enspace \frac{2+\alpha +\sqrt{{\alpha }^{2}+4\alpha }}{2}},\end{equation} \tag{ 34 }$$

and

$$\begin{equation}{\langle \mathcal{F}\rangle }_{\mathrm{as}}=\sqrt{\frac{{\lambda }_{2}}{\delta }}=\sqrt{\frac{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}}{\delta }}=\sqrt{\frac{{h}_{p}}{k}}\sqrt{1-\frac{T}{{T}_{c}}}\quad \text{for}\enspace T > {T}_{0},\end{equation} \tag{ 35 }$$

where

$$\begin{equation}{T}_{c}=\frac{{h}_{p}{y}_{b}^{2}}{{K}_{B}\enspace \mathrm{log}\enspace {\beta }_{0}}=\frac{{h}_{p}{y}_{b}^{2}}{{K}_{B}\enspace \mathrm{log}\enspace \frac{2+\beta +\sqrt{{\beta }^{2}+4\beta }}{2}}.\end{equation} \tag{ 36 }$$

The parameters T_b and T_c represent two critical temperatures for the system with ${\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}=0$ in equations (33) and (35), respectively. As we will discuss in detail, it means that all the bonds are certainly broken for supercritical temperatures. These results also show that the conditions in equation (28), introduced to ensure the convergence of the integrals, correspond to requiring that the system is working at subcritical temperatures. It is possible to observe the emergence of two different cases depending on the parameters of the system. From equations (25) and (36), we find

$$\begin{equation}\frac{{h}_{p}{y}_{b}^{2}}{2{\Delta}E}\frac{{T}_{0}}{{T}_{c}}=\frac{\mathrm{log}\enspace {\beta }_{0}}{\mathrm{log}\enspace {\alpha }_{0}-\mathrm{log}\enspace {\beta }_{0}}.\end{equation} \tag{ 37 }$$

As a consequence, we have

$$\begin{equation}\frac{{h}_{p}{y}_{b}^{2}}{2{\Delta}E}\lessgtr \frac{\mathrm{log}\enspace {\beta }_{0}}{\mathrm{log}\enspace {\alpha }_{0}-\mathrm{log}\enspace {\beta }_{0}}{\Rightarrow}{T}_{c}\lessgtr {T}_{0}.\end{equation} \tag{ 38 }$$

We remark that this analysis is valid if y_b > y_p , i.e. if the softening region exists. It is possible to verify the continuity of ${\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}$ for T = T₀ from equations (33) and (35). As a matter of fact, from the definitions of T₀, T_b and T_c in equations (25), (34) and (36) one can verify the equality

$$\begin{equation}\sqrt{\frac{2{\Delta}E}{k{y}_{b}^{2}}+\frac{{h}_{p}}{k}}\sqrt{1-\frac{{T}_{0}}{{T}_{b}}}=\sqrt{\frac{{h}_{p}}{k}}\sqrt{1-\frac{{T}_{0}}{{T}_{c}}}.\end{equation} \tag{ 39 }$$

Moreover, from this equality we obtain that

$$\begin{equation}\frac{k{y}_{b}^{2}}{2{\Delta}E}=\frac{{T}_{c}}{{T}_{0}}\frac{{T}_{b}-{T}_{0}}{{T}_{c}-{T}_{b}}.\end{equation} \tag{ 40 }$$

Thus, for the case with T_c > T₀ (i.e. when T₀ represent the crossover temperature for the force–extension curves) in equation (38), the positivity of $k{y}_{b}^{2}/2{\Delta}E$ implies that T_c > T_b > T₀.

This scenario is illustrated in figure 3, where equation (18) is compared, for a large value of N, with equations (33) and (35), corresponding to the thermodynamic limit. In order to show the effect of the condition in equation (38), we have fixed the values of the constitutive parameters but different values of the breaking extension y_b in the top and bottom panels (smaller value in top panels). This choice is reflected by the different size of the softening region of the potential energy U in panels (b) and (e). In panel (a) and (d) we plotted the force–extension curves from equation (18) for increasing temperatures and compared it with the asymptotic values obtained from equations (33) and (35). The good agreement observed proves that the theoretical procedure adopted to analyze the thermodynamic limit is correct and that the decreasing of the asymptotic force when the temperature increases is explained by the presence of a phase transition. The critical point depends on the constitutive parameters of the system. When $\frac{{h}_{p}{y}_{b}^{2}}{2{\Delta}E}< \frac{\mathrm{log}\enspace {\beta }_{0}}{\mathrm{log}\enspace {\alpha }_{0}-\mathrm{log}\enspace {\beta }_{0}}$, we have that T₀ > T_c and, thus, the force–extension curve follows equation (33). As a consequence, the phase transition appears at the critical temperature T_b . On the other hand, when $\frac{{h}_{p}{y}_{b}^{2}}{2{\Delta}E} > \frac{\mathrm{log}\enspace {\beta }_{0}}{\mathrm{log}\enspace {\alpha }_{0}-\mathrm{log}\enspace {\beta }_{0}}$, we have T₀ < T_c . In this case, the force–extension follows equation (33) until the temperature reaches the value T₀; when temperature is increased, the behavior of the force is described by equation (35) and the phase transition occurs at the critical temperature T_c . From panels (c) and (f) of figure 3, we deduce that the force needed to obtain the complete detachment of the system monotonically decreases to zero for increasing values of the temperature. This process terminates at T = T_b in panel (c), and at T = T_c in panel (f). It means that the thermal fluctuations are able to promote the fracture of the breakable elements, and that the critical temperatures T_b and T_c are sufficient to induce complete denaturation of the system even if the resulting force experienced by the system is zero. This phenomenon can be better clarified by studying the behavior of $\left\langle N-\eta \right\rangle$ (average number of broken elements) and $\left\langle \eta -\xi \right\rangle$ (average number of softened elements).

Zoom In Zoom Out Reset image size

We will evaluate $\left\langle N-\eta \right\rangle$ and $\left\langle \eta -\xi \right\rangle$ in the limiting case of N → ∞. By performing the analysis used for the force–extension relation, starting from equations (19)–(21), we obtain

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{{\int }_{0}^{N}(N-\eta ){\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}\frac{\mathrm{d}\eta }{{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{1/2}}{\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi }{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}\frac{\mathrm{d}\eta }{{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{1/2}}{\int }_{0}^{\eta }\enspace {\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi },\end{equation} \tag{ 41 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{{\int }_{0}^{N}{\int }_{0}^{\eta }(\eta -\xi ){\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace {\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}\frac{\mathrm{d}\xi \enspace \mathrm{d}\eta }{{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{1/2}}}{{\int }_{0}^{N}{\int }_{0}^{\eta }\enspace {\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace {\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{Y}}^{2}\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}}\frac{\mathrm{d}\xi \enspace \mathrm{d}\eta }{{\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}^{1/2}}}.\end{equation} \tag{ 42 }$$

After the evaluation of the integral on ξ, we apply the change of variable N − η + β₀/(β₀ − 1) = s. We get

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left(s-\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\end{equation} \tag{ 43 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[s-N-\frac{{\beta }_{0}}{{\beta }_{0}-1}+\frac{1}{{\Theta}}\enspace {\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-\frac{1}{{\Theta}}\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\end{equation} \tag{ 44 }$$

where we used ${\Theta}=\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}$ defined in equation (25). As for the force–extension relation, in the limit with N → ∞ we can consider separately the cases with Θ > 0 (T < T₀, exponential term in the brackets is dominant) and Θ < 0 (T > T₀, exponential term in the brackets is negligible). As a result, we get

$$\begin{equation}T< {T}_{0}:\begin{cases}\left\langle N-\eta \right\rangle & =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }\left(s-\frac{{\beta }_{0}}{{\beta }_{0}-1}\right){\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\\ \left\langle \eta -\xi \right\rangle & =\frac{1}{{\Theta}}=\frac{1}{\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}}.\end{cases}\end{equation} \tag{ 45 }$$

$$\begin{equation}T > {T}_{0}:\begin{cases}\left\langle N-\eta \right\rangle & =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }\left(s-\frac{{\beta }_{0}}{{\beta }_{0}-1}\right){\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\mathrm{e}}^{-\frac{1}{s}\delta {\mathcal{Y}}^{2}}\frac{\mathrm{d}s}{{s}^{1/2}}},\\ \left\langle \eta -\xi \right\rangle & \to +\infty .\end{cases}\end{equation} \tag{ 46 }$$

We observe that the average number of softened elements $\left\langle \eta -\xi \right\rangle$ is independent of $\mathcal{Y}$ and is varying only with the temperature T for T < T₀. It means that the two domain walls between intact and softened elements, and between softened and broken elements, move simultaneously with increasing $\mathcal{Y}$ conserving a constant distance between them (which means a constant number of softened elements). On the other hand, the value of $\left\langle \eta -\xi \right\rangle$ is divergent to infinity when T > T₀. This divergence means that we have no intact elements and the vertical springs can be either softened or broken. Therefore, the transition at the temperature T₀ can be characterized by the partial breaking of all the elements of the chain induced by the thermal fluctuations.

The two integrals appearing in equations (45) and (46) can be calculated by using the formulas given in appendix  B. Using the definitions of λ₁, λ₂ and c in equation (30), we obtain

$$\begin{equation}\left\langle N-\eta \right\rangle =\begin{cases}\frac{{\Upsilon}({\lambda }_{1})}{{\Psi}({\lambda }_{1})}-c\quad \text{for}\enspace T< {T}_{0},\quad \\ \frac{{\Upsilon}({\lambda }_{2})}{{\Psi}({\lambda }_{2})}-c\quad \text{for}\enspace T > {T}_{0},\quad \end{cases}\end{equation} \tag{ 47 }$$

where the functions ϒ(λ) and Ψ(λ) are defined as follows

$$\begin{align}{\Upsilon}(\lambda )& =\frac{1}{2}\sqrt{\frac{\pi }{\lambda }}\left\{{\mathrm{e}}^{2\mathcal{Y}\sqrt{\lambda \delta }}\left[1-\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}+\sqrt{\lambda c}\right)\right]\right.\\ & \quad \left.+{\mathrm{e}}^{-2\mathcal{Y}\sqrt{\lambda \delta }}\left[1+\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}-\sqrt{\lambda c}\right)\right]\right\},\end{align} \tag{ 48 }$$

$$\begin{align}{\Psi}(\lambda )& =\frac{1}{4\lambda }\sqrt{\frac{\pi }{\lambda }}\left\{{\mathrm{e}}^{2\mathcal{Y}\sqrt{\lambda \delta }}\left(1-2\mathcal{Y}\sqrt{\lambda \delta }\right)\left[1-\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}+\sqrt{\lambda c}\right)\right]\right.\\ & \quad \left.+{\mathrm{e}}^{-2\mathcal{Y}\sqrt{\lambda \delta }}\left(1+2\mathcal{Y}\sqrt{\lambda \delta }\right)\left[1+\text{erf}\left(\sqrt{\frac{\delta }{c}}\mathcal{Y}-\sqrt{\lambda c}\right)\right]\right\}+\frac{\sqrt{c}}{\lambda }\enspace {\mathrm{e}}^{-\lambda c}\enspace {\mathrm{e}}^{-\frac{\delta }{c}{\mathcal{Y}}^{2}}.\end{align} \tag{ 49 }$$

We can evaluate the asymptotic behavior of $\left\langle N-\eta \right\rangle$ for large values of $\mathcal{Y}$ for both T < T₀ and T > T₀. By means of previous results, and using the asymptotic expression $1-\text{erf}(x)\sim {\mathrm{e}}^{-{x}^{2}}/(x\sqrt{\pi })$ for x → ∞, we easily get the asymptotic formulas for $\mathcal{Y}\to \infty$

$$\begin{equation}\left\langle N-\eta \right\rangle \sim \begin{cases}\mathcal{Y}\sqrt{\frac{\delta }{{\lambda }_{1}}}\quad \text{for}\enspace T< {T}_{0},\quad \\ \mathcal{Y}\sqrt{\frac{\delta }{{\lambda }_{2}}}\quad \text{for}\enspace T > {T}_{0}.\quad \end{cases}\end{equation} \tag{ 50 }$$

The final asymptotic result can be therefore written as

$$\begin{equation}T< {T}_{0}:\begin{cases}\left\langle N-\eta \right\rangle & \sim \mathcal{Y}\sqrt{\frac{\delta }{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi }}=\frac{\mathcal{Y}}{\sqrt{\frac{2{\Delta}E}{k{y}_{b}^{2}}+\frac{{h}_{p}}{k}}\sqrt{1-\frac{T}{{T}_{b}}}}\quad \text{for}\enspace \mathcal{Y}\to \infty ,\\ \left\langle \eta -\xi \right\rangle & =\frac{1}{\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}}=\frac{{K}_{B}T}{{\Delta}E}\frac{1}{1-\frac{T}{{T}_{0}}},\end{cases}\end{equation} \tag{ 51 }$$

$$\begin{equation}T > {T}_{0}:\enspace \begin{cases}\left\langle N-\eta \right\rangle & \sim \mathcal{Y}\sqrt{\frac{\delta }{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}}}=\frac{\mathcal{Y}}{\sqrt{\frac{{h}_{p}}{k}}\sqrt{1-\frac{T}{{T}_{c}}}}\quad \text{for}\enspace \mathcal{Y}\to \infty ,\\ \left\langle \eta -\xi \right\rangle & \to +\infty .\end{cases}\end{equation} \tag{ 52 }$$

This result shows that the increasing of the broken elements is linear with $\mathcal{Y}$ (for $\mathcal{Y}$ large), confirming the progressive detachment process induced by the Helmholtz (isometric) condition. We can combine equations (33) and (35) with equations (51) and (52) to obtain the simple relation $\left\langle N-\eta \right\rangle {\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}=\mathcal{Y}$, which is valid for any temperature T, and for large values of the extension $\mathcal{Y}$. This relation, rewritten in terms of the real physical quantities for large y_N+1, reads

$$\begin{equation}{\left\langle f\right\rangle }_{\mathrm{a}\mathrm{s}}=\frac{k}{\left\langle N-\eta \right\rangle }{y}_{N+1}.\end{equation} \tag{ 53 }$$

The behavior of the system within the Helmholtz ensemble in the thermodynamic limit is summarized in figure 4, where we plot the most important quantities versus the temperature by considering the two cases identified in equation (38). In particular, we fix the values of the constitutive parameters of the system but consider two different values of y_b . For small values of y_b (left panels) the temperature T₀ is larger than T_b and T_c . Therefore, only the equations that are valid for T < T₀ must be considered (panel (a)). On the other hand, for larger values of y_b (right panels), T₀ is smaller that T_b and T_c : thus, the force–extension exhibits different regimes for T < T₀ and T > T₀ (panel (b)). In panel (a) (T_c < T_b < T₀), the asymptotic force ${\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}$ shows a monotone decreasing behavior terminating at the critical temperature T_b corresponding to a phase transition. The number of the softened elements $\left\langle \eta -\xi \right\rangle$ is constant with respect to $\mathcal{Y}$ but is an increasing function of the temperature until the point A, corresponding to the phase transition. It means that the region of softened elements exhibits a fixed length and is gradually moved leftward by the increasing applied extension $\mathcal{Y}$. The extension of this region represents the temperature dependent analytic measure of the process zone. Coherently, the number of broken elements $\left\langle N-\eta \right\rangle$ is linearly increasing with $\mathcal{Y}$ and nonlinearly increasing with T (for this reason we plotted the ratio $\left\langle N-\eta \right\rangle /\mathcal{Y}$ which is only temperature dependent). We observe that $\left\langle N-\eta \right\rangle /\mathcal{Y}$ diverges at the critical temperature T_b . Therefore, $\left\langle N-\eta \right\rangle \to \infty$ for T → T_b for any value of the extension $\mathcal{Y}$ (for large values of $\mathcal{Y}$). It means that the phase transition can be explained through the total detachment of the system induced by the strong thermal fluctuations. In panel (b) (T₀ < T_b < T_c ), the asymptotic force ${\left\langle \mathcal{F}\right\rangle }_{\mathrm{a}\mathrm{s}}$ is a decreasing function of the temperature with a first transition at the point B (temperature T₀) and a second final transition at temperature T_c , where we have the complete detachment. The transition at T₀ can be explained by observing the behavior of the number of softened elements $\left\langle \eta -\xi \right\rangle$ versus the temperature. We see that $\left\langle \eta -\xi \right\rangle \to \infty$ for T → T₀, proving that all elements are (at least) softened for temperatures larger than T₀. Therefore, we have the first transition at T₀ where all the elements are softened, and a second transition at T_c where all the elements are broken. This is coherent with the plot of $\left\langle N-\eta \right\rangle /\mathcal{Y}$ versus the temperature, from which we deduce that $\left\langle N-\eta \right\rangle \to \infty$ for T → T_c for any value of the extension $\mathcal{Y}$. The point C shows the slope transition in the evolution of the broken elements at the temperature T₀. To complete the description of figure 4, we underline that the vertical dashed lines in all panels are useful to easily identify the characteristic temperatures T_c , T_b and T₀. Moreover, the dashed curved lines correspond to solutions that are not relevant to the behavior of the system. More specifically, they correspond to equations (35) and (52) in panel (a) and to equations (33) and (51) in panel (b).

Zoom In Zoom Out Reset image size

To conclude, the behavior of the system under isometric condition can be usefully described by analyzing the phase diagrams represented in figure 5. Two cases are described: T_c < T_b < T₀ (panel (a)); T₀ < T_b < T_c (panel (b)). Observe that, since the system is symmetric (bilateral), in the figure we considered both positive and negative values of $\mathcal{Y}$. We can identify three possible phases for the system: $\mathcal{L}$ corresponds to all the elements broken; $\mathcal{M}$ corresponds to a combination of softened and broken elements; finally, $\mathcal{N}$ corresponds to a mix of intact, softened and broken elements. These phases are separated by vertical lines corresponding to the phase transitions. Therefore, the phase diagrams explain clearly that, for Helmholtz isometric condition, we have no finite values of the prescribed extension $\mathcal{Y}$ able to induce the whole system detachment for subcritical temperatures. We will see that the Gibbs isotensional condition leads to a different behavior.

Zoom In Zoom Out Reset image size

In this section we determine the properties of the detachment process realized within the Gibbs ensemble, which corresponds to a soft device able to apply a force to the last element of the chain. If we consider a force f applied to the last element with vertical extension y_N+1, the total potential energy of the system is given by Φ − fy_N+1, where Φ is the Helmholtz potential energy defined in equation (4). It follows that the Gibbs partition function is given by the Laplace transform of the Helmholtz partition function [92]. Moreover, the knowledge of the Gibbs partition function allows the calculation of the average extension of the last element of the chain by using classical thermodynamic relations [90]. We can introduce the non-dimensional quantities (analogous the those used in section 3)

$$\begin{equation}\mathcal{F}=\frac{f}{k{y}_{b}},\hspace{25.0pt}\left\langle \mathcal{Y}\right\rangle =\frac{\left\langle {y}_{N+1}\right\rangle }{{y}_{b}},\end{equation} \tag{ 54 }$$

representing the rescaled applied force and average extension, respectively. Thus, we can write the relations

$$\begin{equation}{\mathcal{Z}}_{\mathrm{G}}(\mathcal{F})={\int }_{-\infty }^{+\infty }{\mathcal{Z}}_{\mathrm{H}}(\mathcal{Y}){\mathrm{e}}^{2\delta \mathcal{F}\mathcal{Y}}\enspace \mathrm{d}\mathcal{Y},\end{equation} \tag{ 55 }$$

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\frac{1}{2\delta }\frac{1}{{\mathcal{Z}}_{\mathrm{G}}(\mathcal{F})}\frac{\partial {\mathcal{Z}}_{\mathrm{G}}(\mathcal{F})}{\partial \mathcal{F}}.\end{equation} \tag{ 56 }$$

We can easily determine the Gibbs partition function through equations (21) and (55) as a Gaussian integral. We obtain (up to a multiplicative constant)

$$\begin{equation}{\mathcal{Z}}_{\mathrm{G}}(\mathcal{F})=\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{{\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }}{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}}}\enspace {\mathrm{e}}^{\frac{\delta {\mathcal{F}}^{2}}{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )}},\end{equation} \tag{ 57 }$$

where δ and φ have been defined in equation (17). By applying equation (56), we directly obtain the force–extension relation as

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\frac{\mathcal{F}}{{\mathcal{Z}}_{\mathrm{G}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{{\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }}{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}{\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}}^{3/2}}\enspace {\mathrm{e}}^{\frac{\delta {\mathcal{F}}^{2}}{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )}}.\end{equation} \tag{ 58 }$$

We can finally calculate the average value of the number of broken elements $\left\langle N-\eta \right\rangle$ and the average number of softened elements $\left\langle \eta -\xi \right\rangle$, eventually obtaining

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{1}{{\mathcal{Z}}_{\mathrm{G}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{N-\eta }{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}}}\enspace {\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }\enspace {\mathrm{e}}^{\frac{\delta {\mathcal{F}}^{2}}{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )}},\end{equation} \tag{ 59 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{1}{{\mathcal{Z}}_{\mathrm{G}}}\sum\limits _{\eta =0}^{N}\sum\limits _{\xi =0}^{\eta }\frac{\eta -\xi }{\sqrt{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )\left\{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )\right\}}}{\mathrm{e}}^{\varphi \xi }\enspace {\text{e}}^{\beta \delta \eta }\enspace {\mathrm{e}}^{\frac{\delta {\mathcal{F}}^{2}}{1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )}}.\end{equation} \tag{ 60 }$$

An example of application of these results can be found in figure 6, where we plotted $\left\langle \mathcal{Y}\right\rangle$, $\left\langle \xi \right\rangle$, $\left\langle \eta -\xi \right\rangle$ and $\left\langle N-\eta \right\rangle$ versus the applied force $\mathcal{F}$ (using the temperature as a parameter) for two different cases: the first one concerns a short chain with N = 5 (panel (a)), and the second one a larger chain with N = 50 (panel (b)). Concerning the case with N = 5, it is interesting to remark the differences between the Gibbs and the Helmoltz ensemble. We therefore compare the panel (a) of figure 2 with the panel (a) of figure 6. While the detachment of the elements is quite simultaneous within the Gibbs ensemble (cooperative process occurring at a given threshold force), the same phenomenon is sequential within the Helmholtz ensemble, where we observe a gradual and progressive bonds breaking (non-cooperative response). This different behavior is well known in the context of the stretching of macromolecules by means of force-spectroscopy techniques [83–86, 88, 89]. Another difference concerns the fact that in the Gibbs ensemble we are not able to identify peaks or steps in the curves of figure 6, corresponding to the single breaking occurrences. Also this point is clearly related to the cooperative response characterizing the Gibbs ensemble. If we look at the case with N = 50 (panel (b) of figure 6), we can see a well identified force plateau corresponding to the cooperative detachment of the elements. This force plateau is temperature dependent, as already observed within the Helmholtz ensemble. Moreover, we also deduce that the number of softened elements remains approximately constant and only temperature dependent before the complete rupture of the elements. These properties of the detachment process will be thoroughly analyzed in the next section, by introducing the thermodynamic limit of the Gibbs ensemble.

Zoom In Zoom Out Reset image size

In this section, we analyze the thermodynamic limit within the Gibbs ensemble and consider the behavior of equations (57)–(60) for N → ∞. Following the procedure shown in section 4 and using the formulas in appendix A, equations (57) and (58) give the force–extension response of the system within the Gibbs ensemble as

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\frac{{\int }_{0}^{N}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right){\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\mathrm{d}\xi }{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\mathrm{d}\xi },\end{equation} \tag{ 61 }$$

where α₀ and β₀ have been defined in equation (23). After the integration over ξ we can apply the change of variable N − η + β₀/(β₀ − 1) = s thus obtaining

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}s\enspace \mathrm{d}s}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}\enspace \mathrm{d}s}.\end{equation} \tag{ 62 }$$

As before, this expression depends on the sign of the quantity ${\Theta}=\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}$. In the limit with N → ∞ we consider separately the cases with Θ > 0 (T < T₀, exponential term in the brackets is dominant) and Θ < 0 (T > T₀, exponential term in the brackets is negligible). We eventually get (for N → ∞)

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}s\enspace \mathrm{d}s}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi \right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}\enspace \mathrm{d}s}\quad \text{for}\enspace T< {T}_{0},\end{equation} \tag{ 63 }$$

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}s\enspace \mathrm{d}s}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{+\infty }{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)s}\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}s}\enspace \mathrm{d}s}\quad \text{for}\enspace T > {T}_{0}.\end{equation} \tag{ 64 }$$

Finally, we obtain the force–extension formulas for N → ∞ under the same conditions for the Helmholtz ensemble (see equations (28) and (29)) as

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\left(\frac{{\beta }_{0}}{{\beta }_{0}-1}+\frac{1}{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi -\delta {\mathcal{F}}^{2}}\right)\quad \text{for}\enspace T< {T}_{0},\end{equation} \tag{ 65 }$$

$$\begin{equation}\left\langle \mathcal{Y}\right\rangle =\mathcal{F}\left(\frac{{\beta }_{0}}{{\beta }_{0}-1}+\frac{1}{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}}\right)\quad \text{for}\enspace T > {T}_{0}.\end{equation} \tag{ 66 }$$

We remark that these expressions reveal an asymptotic behavior corresponding to a force plateau, as already observed within the Helmholtz ensemble. In particular, the values of the asymptotic force are the same for the two statistical ensembles. Indeed, for T < T₀, we deduce from equation (65) that the asymptotic force is solution of the quadratic equation $\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi -\delta {\mathcal{F}}^{2}=0$, which is exactly solved by the Helmholtz result stated in equation (33). Similarly, for T > T₀, we deduce from equation (66) that the asymptotic force is solution of the quadratic equation $\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}=0$, which is exactly solved by the Helmholtz result stated in equation (35). While the asymptotic forces are the same for the Helmholtz and Gibbs ensembles, it is not difficult to verify that the shape of the force–extension response is sensibly different for the two ensembles. This can be done by comparing the analytical solutions for the force–extension response obtained at the thermodynamic limit in equations (31) and (32) (Helmholtz) and equations (65) and (66) (Gibbs). It means that the two statistical ensembles are nonequivalent in the thermodynamic limit. For instance, this point can be shown by calculating the slopes of the force–extensions curves at the origin (for small forces and extensions) for both ensembles and by observing that these quantities are different. More details, concerning the case of the adhesion-decohesion process without the softening mechanism, can be found in [80].

To complete the analysis, we determine the average value of the number of broken elements $\left\langle N-\eta \right\rangle$ and the average number of softened elements $\left\langle \eta -\xi \right\rangle$ in the limiting case of N → ∞. To this aim, we apply the same procedure already used for the force–extension relation. From equations (59) and (60), combined with the partition function given in equation (57), we obtain

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{{\int }_{0}^{N}(N-\eta ){\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi }{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi },\end{equation} \tag{ 67 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }(\eta -\xi ){\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi }{{\int }_{0}^{N}{\mathrm{e}}^{\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}\right)\eta }\enspace {\text{e}}^{-\delta {\mathcal{F}}^{2}\left(N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}\mathrm{d}\eta {\int }_{0}^{\eta }{\mathrm{e}}^{\left(\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}\right)\xi }\enspace \mathrm{d}\xi }.\end{equation} \tag{ 68 }$$

After the integration over ξ and the change of variable N − η + β₀/(β₀ − 1) = s we obtain

$$\begin{equation}\left\langle N-\eta \right\rangle =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left(s-\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}\right)s}\enspace \mathrm{d}s}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}\right)s}\enspace \mathrm{d}s},\end{equation} \tag{ 69 }$$

$$\begin{equation}\left\langle \eta -\xi \right\rangle =\frac{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[-N+s-\frac{{\beta }_{0}}{{\beta }_{0}-1}+\frac{1}{{\Theta}}\enspace {\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-\frac{1}{{\Theta}}\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}\right)s}\enspace \mathrm{d}s}{{\int }_{\frac{{\beta }_{0}}{{\beta }_{0}-1}}^{\frac{{\beta }_{0}}{{\beta }_{0}-1}+N}\left[{\mathrm{e}}^{{\Theta}\left(N-s+\frac{{\beta }_{0}}{{\beta }_{0}-1}\right)}-1\right]{\mathrm{e}}^{-\left(\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}\right)s}\enspace \mathrm{d}s},\end{equation} \tag{ 70 }$$

where, as in section 4, we considered ${\Theta}=\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}$. To complete the analysis, we consider the cases with Θ > 0 (T < T₀, exponential term in the brackets is dominant) and Θ < 0 (T > T₀, exponential term in the brackets is negligible). We get

$$\begin{equation}T< {T}_{0}:\begin{cases}\left\langle N-\eta \right\rangle & =\frac{1}{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi -\delta {\mathcal{F}}^{2}},\\ \left\langle \eta -\xi \right\rangle & =\frac{1}{{\Theta}}=\frac{1}{\varphi +\frac{1}{2}\enspace \mathrm{log}\enspace \frac{{\beta }_{0}}{{\alpha }_{0}}},\end{cases}\end{equation} \tag{ 71 }$$

$$\begin{equation}T > {T}_{0}:\begin{cases}\left\langle N-\eta \right\rangle & =\frac{1}{\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}},\\ \left\langle \eta -\xi \right\rangle & \to +\infty .\end{cases}\end{equation} \tag{ 72 }$$

These results can be used to understand the behavior of the detachment process within the Gibbs ensemble in the thermodynamic limit as shown in figures 7 and 8. In both figures, we have fixed the constitutive parameters of the model but the values of y_b (small value in the left panel, larger value in the right panel). We recall that in the Helmholtz ensemble we obtained two different responses depending on the parameters, as stated in equation (38). In panel (a) of figures 7 and 8 we show the results for the case with T_c < T_b < T₀ (y_b = 1.5y_p ), and in panel (b) those of the case with T₀ < T_b < T_c (y_b = 3.5y_p ). While in figure 7 we show the behavior of the relevant quantities as function of the applied force (using the temperature as a parameter), in figure 8 we show the same quantities as function of the temperature (now using the applied force as a parameter).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 7 we draw a comparison between equations (57)–(60) with a large value of N and equations (65), (66), (71) and (72) obtained for N → ∞. This comparison prove that the mathematical treatment of the thermodynamic limit here developed is consistent with the model of detachment proposed. Indeed, we find a very good agreement between the general results, applied to system with a large value of N, and the theoretical ones obtained in the thermodynamic limit. Concerning panel (a), the force extension-curves exhibit the force plateau discussed above, which is very similar to the behavior observed within the Helmholtz ensemble. The average number of softened elements remains constant with $\mathcal{F}$ and depends on the temperature T. Finally, the number of broken elements diverges to infinity when the force approaches its critical value. The situation is more complicated in panel (b), where we have an intermediate transition describing the partial breaking of all the element of the system. In this case, when T < T₀ (black solid curves), the behavior is very similar to the one described in panel (a). On the other hand, when T > T₀ (red solid curves), all the elements are softened due to the thermal fluctuations and therefore $\left\langle \eta -\xi \right\rangle \to \infty$ (for this reason no red curves are shown in the plot of $\left\langle \eta -\xi \right\rangle$ versus $\mathcal{F}$). It means that the detachment process consists in the complete break of all the elements already softened.

Also in figure 8, we separately show the behavior of the system with T_c < T_b < T₀ (panel (a) with y_b = 1.5y_p ), and the system with T₀ < T_b < T_c (panel (b) with y_b = 3.5y_p ). We used a logarithmic scale to better appreciate the shape of the plots. The representation of the extension $\left\langle \mathcal{Y}\right\rangle$, the average number of softened elements $\left\langle \eta -\xi \right\rangle$, and average number of broken elements $\left\langle N-\eta \right\rangle$ versus the temperature is particularly useful to identify a crucial difference between the Helmholtz and the Gibbs ensembles. For the isometric condition (Helmholtz ensemble), we have no finite values of the prescribed extension $\mathcal{Y}$ able to completely detach the system for subcritical temperatures. Indeed, since the detachment progress is gradual or progressive in response to the prescribed extension, only an infinite value of extension can cause a complete detachment of the chain (see figure 4). Differently, for the isotensional condition (Gibbs ensemble), for any subcritical temperature, we can identify a value of force inducing the complete detachment of the chain. This can be seen by looking at the force–extension curves in panel (a) and (b) of figure 8. Here, we observe that for any value of the temperature, we can identify a force producing an asymptotic behavior of $\left\langle \mathcal{Y}\right\rangle$, which corresponds to the total detachment of the system. The relationships between force and temperature are given by $\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\alpha }_{0}+\varphi -\delta {\mathcal{F}}^{2}=0$ for T < T₀, and by $\beta \delta -\frac{1}{2}\enspace \mathrm{log}\enspace {\beta }_{0}-\delta {\mathcal{F}}^{2}=0$ for T > T₀, as stated in equations (65) and (66) (we remember that δ and φ, defined in equation (17), depend on T). This is coherent with the fact that within the Gibbs ensemble all the element break simultaneously or cooperatively, thus allowing the force to completely detach the chain independently of the (subcritical) temperature of the system. The average number of softened elements $\left\langle \eta -\xi \right\rangle$ is independent of the force and therefore we find a single curve ($\forall \enspace \mathcal{F}$) in the corresponding plot. However, for any applied force we have a different critical temperature of the system and then the curve of $\left\langle \eta -\xi \right\rangle$ must be considered valid up to the red circle symbols shown in figure 8. Indeed, these symbols correspond to the critical temperature, which is, in turn, determined by the applied force. Finally, the curves of the average number of broken elements $\left\langle N-\eta \right\rangle$ show a divergence at the critical temperature of the system, representing the complete detachment of the chain from the substrate. The important difference between panel (a) and panel (b) consists in the presence of the additional transition at the temperature T₀ in the case with T₀ < T_b < T_c (large breaking extension y_b = 3.5y_p ). In this case, we observe that for T → T₀ we have $\left\langle \eta -\xi \right\rangle \to \infty$ even if $\left\langle N-\eta \right\rangle$ is finite. It means that at the temperature T ⩾ T₀ all the vertical elements are spontaneously softened (due to the thermal fluctuations), and the detachment process is realized by completing their rupture thanks to the applied force.

To conclude the discussion concerning the thermodynamic limit in the Gibbs ensemble, in figure 9 we show the phase diagram of the system in the force–temperature plane. We separately considered the two cases with T_c < T_b < T₀ (panel (a)) and with T₀ < T_b < T_c (panel (b)). Again, due to symmetry, we considered both positive and negative values of $\mathcal{F}$ (see the scheme in figure 1 (panel (b))). Moreover, we observe the same three possible cases obtained in the hard device: the region $\mathcal{L}$ corresponds to a phase with all broken elements; the region $\mathcal{M}$ corresponds to a combination of softened and broken elements; finally, the region $\mathcal{N}$ corresponds to a mixture of intact, softened and broken elements. In figure 9 we can identify the critical points (C_b in panel (a), and C_c in panel (b)), whose meaning is now quite clear: for temperatures larger than those of the critical points, the chain is always completely detached from the substrate. In panel (b) of figure 9 we can see the triple points P⁺ and P⁻, where the three phases of the system coexist. In any case, we observe that for a value of the applied force $\mathcal{F}$ (in a suitable range), we have a value of the temperature able to detach the chain. When $\mathcal{F}=0$, such a temperature corresponds to the critical temperature of the system. In other words, for any subcritical temperature, we can identify a value of force inducing the complete detachment of the chain. As previously discussed, the critical behavior of the Helmholtz ensemble is different. In fact, for isometric condition, we have no finite values of the prescribed extension able to detach the system for subcritical temperatures, as shown in figure 5. This difference between the Helmholtz and the Gibbs ensembles is consistent with their non-cooperative and cooperative interpretation, respectively.

Zoom In Zoom Out Reset image size

We considered the cohesion/decohesion process under thermal fluctuations by means of a paradigmatic system represented by a linear mass-spring chain, connected to a substrate through breakable elements. These elements are described by a softening mechanism, i.e. by three operating regions as follows: the first linear region corresponds to intact elements (small deformation), the second region corresponds to softened elements with a reduced stiffness (intermediate deformation), and finally, the third region corresponds to broken elements (large deformation). This paradigmatic system has been studied under the assumption to be embedded in a thermal bath at a given temperature and to be loaded by an external mechanical action. This latter can be either a prescribed extension (hard device) or an applied force (soft device) to the last element of the chain. These two configurations correspond to the Helmholtz ensemble (isometric condition) and the Gibbs ensemble (isotensional condition) of the equilibrium statistical mechanics, respectively. In order to cope with the mathematical analysis of the problem, we adopted the spin variable approach for defining the state of any breakable element [83–86, 88, 89], we introduced the zipper model hypothesis with two domain walls for describing the three operating regions of the bonds chain [97–100], and we finally exploited some analytical properties of the tridiagonal matrices [116, 117] to obtain closed form results under isometric and isotensional conditions. In both cases we obtained the general solution of the problem for a chain composed of N elements. These solutions, in both cases, lead to a constant decohesion force, which represent the strength of the system. Importantly, this strength is temperature dependent, and its value is decreasing with an increasing temperature. The physical interpretation of this behavior can be done by observing that the thermal fluctuations foster the damaging of the breakable elements, thus reducing the strength of the system. Also the number of broken and softened elements have been determined as function of the temperature and of the external mechanical action for both statistical ensembles. Interestingly, we proved that the number of softened elements depends only on the temperature and it is independent of the external mechanical action. The increasing size of the softened region with the temperature gives a thorough interpretation of the fragile to ductile transition of several materials. This analysis provides a complete picture of the system behavior with a finite number N of elements. To better understand the meaning of these results, we analyzed the system behavior under the hypothesis of thermodynamic limit, i.e. for N → ∞. In this case, the sums appearing in the partition functions of both ensembles and in the macroscopic averages of the relevant thermodynamic variables can be substituted by integrals, eventually leading to useful closed form expressions. The analysis of these results delivers a clear explanation of the temperature dependent strength in terms of phase transitions. Indeed, a first critical temperature exists for the system and corresponds to the partial breaking of all the elements of the chain. Then, a second critical temperature describes the complete breaking of all the bonds. This complex critical scenario is able to give a rigorous explanation of the relationship between strength and temperature in system where a transition in the breaking mechanism occurs (here represented by the bonds softening). In several experimental measurements, a two-branch curve for the strength-temperature relation has been indeed observed [30, 32, 35]. Another interesting result concerns the non-equivalence of the ensembles in the thermodynamic limit. It means that the force–extension curves are different for the Helmholtz and the Gibbs ensembles, also for N → ∞. It is well known that the equivalence of the Gibbs and Helmholtz statistical ensemble has been proved for a large class of physical systems including, e.g. single flexible polymer chains without confinement effects and with a continuous pairing interaction potential between neighboring monomers [91–96]. However, it is difficult to prove the ensembles equivalence for more general systems or structures. For instance, other problems, concerning for example the escape of a polymer confined between two planes or the desorption of a polymer tethered on a surface, lead to the non-equivalence between the defined statistical ensembles [118–123]. Interestingly, these investigations are coherent with our results and definitively prove the possibility to have non-equivalence between different statistical ensembles in relevant physical systems. Generalizations of the theory may include the dynamic behavior of cohesion/decohesion processes (where the barriers between the wells of the energetic landscape play an important role [124–126]), or the multi-state softening i.e. composed of several steps with different stiffness (described by several domain walls in the chain).

GF and GP have been supported by the Italian Ministry MIUR-PRIN project 'Mathematics of active materials: from mechanobiology to smart devices' (2017KL4EF3). GF and GP are supported by GNFM (INdAM) and by the Italian Ministry MISE through the project RAEE SUD-PVP. GF is also supported by INFN through the project QUANTUM, by the FFABR research grant (MIUR) and the PON S.I.ADD. SG has been supported by 'Central Lille' and 'Région Hauts-de-France' under projects MePoFib and MiBaMs.

No new data were created or analysed in this study.

In this appendix we will recall some properties of tridiagonal matrices used in the text. In particular, we consider the matrix $\mathcal{A}(\xi ,\eta )$ defined in equation (5) and obtain the formulas for $\mathrm{det}\left[\mathcal{A}(\xi ,\eta )\right]$ and $1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )$.

For sake's of generality, let us consider $\mathcal{T}$ defined as

$$\begin{equation}\mathcal{T}=\left[\begin{matrix}\hfill {a}_{1}\hfill & \hfill {b}_{1}\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill {c}_{1}\hfill & \hfill {a}_{2}\hfill & \hfill {b}_{2}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill 0\hfill & \hfill {c}_{2}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill 0\hfill \\ \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {a}_{N-1}\hfill & \hfill {b}_{N-1}\hfill \\ \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill {c}_{N-1}\hfill & \hfill {a}_{N}\hfill \end{matrix}\right]\in {\mathcal{M}}_{N,N}(\mathbb{R}),\end{equation} \tag{ A.1 }$$

where the diagonal is composed by the elements (a₁, ⋯, a_N ), the superdiagonal by (b₁, ⋯, b_N−1) and the subdiagonal by (c₁, ⋯, c_N−1).

The elements of the inverse matrix ${\mathcal{T}}^{-1}$ can be written as [116, 117]

$$\begin{equation}{\mathcal{T}}_{ij}^{-1}=\begin{cases}\frac{1}{{\vartheta }_{N}}{(-1)}^{i+j}{b}_{i}\times \cdots \times {b}_{j-1}{\vartheta }_{i-1}{\varphi }_{j+1},& i< j\\ \frac{1}{{\vartheta }_{N}}{\vartheta }_{i-1}{\varphi }_{i+1},& i=j\\ \frac{1}{{\vartheta }_{N}}{(-1)}^{i+j}{c}_{j}\times \cdots \times {c}_{i-1}{\vartheta }_{j-1}{\varphi }_{i+1},& i > j\end{cases},\end{equation} \tag{ A.2 }$$

where the sequences ϑ_i and φ_i are obtained using the recursive rules

$$\begin{equation}\begin{cases}{\vartheta }_{i}={a}_{i}{\vartheta }_{i-1}-{b}_{i-1}{c}_{i-1}{\vartheta }_{i-2},\quad i=1,\cdots ,N\quad \\ {\vartheta }_{-1}=0,\hspace{25.0pt}{\vartheta }_{0}=1,\quad \end{cases}\end{equation} \tag{ A.3 }$$

and

$$\begin{equation}\begin{cases}{\varphi }_{i}={a}_{i}{\varphi }_{i+1}-{b}_{i}{c}_{i}{\varphi }_{i+2},\quad i=N,\cdots ,1\quad \\ {\varphi }_{N+1}=1,\hspace{25.0pt}{\varphi }_{N+2}=0.\quad \end{cases}\end{equation} \tag{ A.4 }$$

Notice that while equation (A.3) is an increasing recursive rule going from i = 1 to i = N, equation (A.4) is a decreasing recursive rule for i = N to i = 1. It is possible to show [116, 117] that

$$\begin{equation}\mathrm{det}\enspace \mathcal{T}={\vartheta }_{N}.\end{equation} \tag{ A.5 }$$

These formulas can be specialized to the case where b_i = c_i = −1 ∀ i (such as for matrix $\mathcal{A}(\xi ,\eta )$ in equation (5)). The general result can thus be simplified as

$$\begin{equation}{\mathcal{T}}_{ij}^{-1}=\begin{cases}\frac{1}{{\vartheta }_{N}}{\vartheta }_{i-1}{\varphi }_{j+1},& i< j\\ \frac{1}{{\vartheta }_{N}}{\vartheta }_{i-1}{\varphi }_{i+1},& i=j\\ \frac{1}{{\vartheta }_{N}}{\vartheta }_{j-1}{\varphi }_{i+1},& i > j\end{cases},\end{equation} \tag{ A.6 }$$

where the sequences ϑ_i and φ_i are given by the reduced recursive laws

$$\begin{equation}\begin{cases}{\vartheta }_{i}={a}_{i}{\vartheta }_{i-1}-{\vartheta }_{i-2},\quad i=1,\cdots ,N,\quad \\ {\vartheta }_{-1}=0,\hspace{25.0pt}{\vartheta }_{0}=1,\quad \end{cases}\end{equation} \tag{ A.7 }$$

and

$$\begin{equation}\begin{cases}{\varphi }_{i}={a}_{i}{\varphi }_{i+1}-{\varphi }_{i+2},\quad i=N,\cdots ,1,\quad \\ {\varphi }_{N+2}=0,\hspace{25.0pt}{\varphi }_{N+1}=1.\quad \end{cases}\end{equation} \tag{ A.8 }$$

In order to evaluate the physical quantities in sections 3 and 5 we need to compute $\mathrm{det}\left[\mathcal{A}(\xi ,\eta )\right]$ and ${\mathcal{A}}_{NN}^{-1}(\xi ,\eta )$. We have that

$$\begin{equation}\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )={\vartheta }_{N},\hspace{25.0pt}{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )={\vartheta }_{N-1}{\varphi }_{N+1}/{\vartheta }_{N}={\vartheta }_{N-1}/{\vartheta }_{N}.\end{equation} \tag{ A.9 }$$

Let us consider equation (A.7) with a₁ = ⋯ = a_ξ = 2 + α, a_ξ+1 = ⋯ = a_η = 2 + β, and a_η+1 = ⋯ = a_N = 2. Then, for i ⩽ ξ we have the difference equation ϑ_i = (2 + α)ϑ_i−1 − ϑ_i−2, whose general solution can be written as

$$\begin{equation}{\vartheta }_{i}=p{\left(\frac{2+\alpha +\sqrt{{\Delta}}}{2}\right)}^{i}+q{\left(\frac{2+\alpha -\sqrt{{\Delta}}}{2}\right)}^{i},\end{equation} \tag{ A.10 }$$

with Δ = α² + 4α and where the coefficients p and q must be fixed using the conditions ϑ₋₁ = 0 and ϑ₀ = 1. A straightforward calculation leads to the explicit solution for i ⩽ ξ

$$\begin{equation}{\vartheta }_{i}=\frac{1}{\sqrt{{\Delta}}}{\left(\frac{2+\eta +\sqrt{{\Delta}}}{2}\right)}^{i+1}-\frac{1}{\sqrt{{\Delta}}}{\left(\frac{2+\eta -\sqrt{{\Delta}}}{2}\right)}^{i+1}.\end{equation} \tag{ A.11 }$$

Let us define the function

$$\begin{equation}\mathcal{G}(\gamma ,z)=\frac{1}{\sqrt{{\gamma }^{2}+4\gamma }}\left[{\left(\frac{2+\gamma +\sqrt{{\gamma }^{2}+4\gamma }}{2}\right)}^{z}-{\left(\frac{2+\gamma -\sqrt{{\gamma }^{2}+4\gamma }}{2}\right)}^{z}\right],\end{equation} \tag{ A.12 }$$

which satisfies the properties $\mathcal{G}(\gamma ,0)=0$, $\mathcal{G}(\gamma ,1)=1$, and $\mathcal{G}(\gamma ,-1)=-1$. We can write the final result for i ⩽ ξ as

$$\begin{equation}{\vartheta }_{i}=\mathcal{G}(\alpha ,i+1).\end{equation} \tag{ A.13 }$$

For ξ + 1 ⩽ i ⩽ η, we have the difference equation ϑ_i = (2 + β)ϑ_i−1 − ϑ_i−2, with the general solution

$$\begin{equation}{\vartheta }_{i}=r{\left(\frac{2+\beta +\sqrt{{\Sigma}}}{2}\right)}^{i}+s{\left(\frac{2+\beta -\sqrt{{\Sigma}}}{2}\right)}^{i},\end{equation} \tag{ A.14 }$$

where Σ = β² + 4β. In this case, the coefficients r and s must be obtained by imposing ϑ_ξ−1 and ϑ_ξ by means of equation (A.13). A long but straightforward calculation yields the final solution for ξ + 1 ⩽ i ⩽ η in the form

$$\begin{equation}{\vartheta }_{i}=\mathcal{G}(\beta ,i-\xi +1)\mathcal{G}(\alpha ,\xi +1)-\mathcal{G}(\beta ,i-\xi )\mathcal{G}(\alpha ,\xi ).\end{equation} \tag{ A.15 }$$

For i ⩾ η + 1, we have the simpler difference equation ϑ_i = 2ϑ_i−1 − ϑ_i−2, with the general solution ϑ_i = v + wi. In this case, the coefficients v and w must be obtained by imposing ϑ_η−1 and ϑ_η by means of equation (A.15). Hence, the result for i ⩾ η + 1 can be eventually found as

$$\begin{equation}{\vartheta }_{i}=(i-\eta +1){\vartheta }_{\eta }-(i-\eta ){\vartheta }_{\eta -1},\end{equation} \tag{ A.16 }$$

where ϑ_η−1 and ϑ_η are given by equation (A.15). The result for i ⩾ η + 1 can be therefore written in the following final form

$$\begin{align}{\vartheta }_{i}& =(i-\eta +1)\left[\mathcal{G}(\beta ,\eta -\xi +1)\mathcal{G}(\alpha ,\xi +1)-\mathcal{G}(\beta ,\eta -\xi )\mathcal{G}(\alpha ,\xi )\right]\\ & \quad -(i-\eta )\left[\mathcal{G}(\beta ,\eta -\xi )\mathcal{G}(\alpha ,\xi +1)-\mathcal{G}(\beta ,\eta -\xi -1)\mathcal{G}(\alpha ,\xi )\right].\end{align} \tag{ A.17 }$$

From equation (A.9) we finally find

$$\begin{align}\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )& =(N-\eta +1)\left[\mathcal{G}(\beta ,\eta -\xi +1)\mathcal{G}(\alpha ,\xi +1)-\mathcal{G}(\beta ,\eta -\xi )\mathcal{G}(\alpha ,\xi )\right]\\ & \quad -(N-\eta )\left[\mathcal{G}(\beta ,\eta -\xi )\mathcal{G}(\alpha ,\xi +1)-\mathcal{G}(\beta ,\eta -\xi -1)\mathcal{G}(\alpha ,\xi )\right],\end{align} \tag{ A.18 }$$

and

$$\begin{align}1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )& =\frac{\mathcal{G}(\alpha ,\xi +1)}{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}\left[\mathcal{G}(\beta ,\eta -\xi +1)-\mathcal{G}(\beta ,\eta -\xi )\right]\\ & \quad -\frac{\mathcal{G}(\alpha ,\xi )}{\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )}\left[\mathcal{G}(\beta ,\eta -\xi )-\mathcal{G}(\beta ,\eta -\xi -1)\right].\end{align} \tag{ A.19 }$$

These formulas are used in sections 3 and 5.

The function $\mathcal{G}(\gamma ,z)$ defined in equation (A.12) can be simplified by observing that the second term is defined by a power with a base always strictly less than one. Thus, for large values of z, the second term can be neglected with respect to the first one. More specifically, we have two different cases:

$$\begin{equation}\mathcal{G}(\alpha ,z)\sim \frac{1}{\sqrt{{\Delta}}}{\alpha }_{0}^{z}\quad \text{with}\quad {\alpha }_{0}=(2+\alpha +\sqrt{{\Delta}})/2,\enspace {\Delta}={\alpha }^{2}+4\alpha ,\end{equation} \tag{ A.20 }$$

$$\begin{equation}\mathcal{G}(\beta ,z)\sim \frac{1}{\sqrt{{\Sigma}}}{\beta }_{0}^{z}\quad \text{with}\quad {\beta }_{0}=(2+\beta +\sqrt{{\Sigma}})/2,\enspace {\Sigma}={\beta }^{2}+4\beta .\end{equation} \tag{ A.21 }$$

These results can be used to simplify equations (A.18) and (A.19) leading to

$$\begin{equation}\mathrm{det}\enspace \mathcal{A}(\xi ,\eta )=\frac{{\alpha }_{0}{\beta }_{0}-1}{\sqrt{{\Delta}{\Sigma}}}{\beta }_{0}^{\eta -\xi -1}{\alpha }_{0}^{\xi }\left[(N-\eta +1){\beta }_{0}-(N-\eta )\right],\end{equation} \tag{ A.22 }$$

$$\begin{equation}1-{\mathcal{A}}_{NN}^{-1}(\xi ,\eta )=\frac{1}{N-\eta +\frac{{\beta }_{0}}{{\beta }_{0}-1}}.\end{equation} \tag{ A.23 }$$

These expressions are used to evaluate the formulas in sections 4 and 6.

In order to develop equations (26) and (27) in closed form we have to consider the integrals

$$\begin{equation}{I}_{1}=\int {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\frac{\mathrm{d}x}{{x}^{3/2}},\hspace{25.0pt}{I}_{2}=\int {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\frac{\mathrm{d}x}{{x}^{1/2}},\end{equation} \tag{ B.1 }$$

where a and b are two real and positive constants. For obtaining an explicit expression for these integrals, we define the quantities

$$\begin{equation}{I}^{+}=\sqrt{b}{I}_{1}+\sqrt{a}{I}_{2},\hspace{25.0pt}{I}^{-}=\sqrt{b}{I}_{1}-\sqrt{a}{I}_{2}.\end{equation} \tag{ B.2 }$$

We can therefore write

$$\begin{equation}{I}^{\mp }=\int \left(\sqrt{\frac{b}{{x}^{3}}}\mp \sqrt{\frac{a}{x}}\right){\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\mathrm{d}x.\end{equation} \tag{ B.3 }$$

Now, we can observe that

$$\begin{equation}{\left(\sqrt{\frac{b}{x}}\pm \sqrt{ax}\right)}^{2}=\frac{b}{x}+ax\pm 2\sqrt{ab},\end{equation} \tag{ B.4 }$$

and hence, we can write equation (B.3) in the form

$$\begin{equation}{I}^{\mp }=\int \left(\sqrt{\frac{b}{{x}^{3}}}\mp \sqrt{\frac{a}{x}}\right){\mathrm{e}}^{-{\left(\sqrt{\frac{b}{x}}\pm \sqrt{ax}\right)}^{2}\pm 2\sqrt{ab}}\enspace \mathrm{d}x.\end{equation} \tag{ B.5 }$$

To conclude, by remembering the definition of the error function [115]

$$\begin{equation}\text{erf}(x)=\frac{2}{\sqrt{\pi }}{\int }_{0}^{x}{\mathrm{e}}^{-{t}^{2}}\enspace \mathrm{d}t,\end{equation} \tag{ B.6 }$$

we can observe that

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}x}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}\pm \sqrt{ax}\right)=-\frac{1}{\sqrt{\pi }}\enspace {\mathrm{e}}^{-{\left(\sqrt{\frac{b}{x}}\pm \sqrt{ax}\right)}^{2}}\left(\sqrt{\frac{b}{{x}^{3}}}\mp \sqrt{\frac{a}{x}}\right),\end{equation} \tag{ B.7 }$$

and we directly obtain the first important result

$$\begin{equation}{I}^{\mp }=-\sqrt{\pi }\enspace {\mathrm{e}}^{\pm 2\sqrt{ab}}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}\pm \sqrt{ax}\right).\end{equation} \tag{ B.8 }$$

Since from equation (B.2) we have that ${I}^{+}+{I}^{-}=2\sqrt{b}{I}_{1}$ and ${I}^{+}-{I}^{-}=2\sqrt{a}{I}_{2}$, we get the final results

$$\begin{equation}{I}_{1}=-\frac{1}{2}\sqrt{\frac{\pi }{b}}\left[{\mathrm{e}}^{+2\sqrt{ab}}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}+\sqrt{ax}\right)+{\mathrm{e}}^{-2\sqrt{ab}}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}-\sqrt{ax}\right)\right],\end{equation} \tag{ B.9 }$$

$$\begin{equation}{I}_{2}=+\frac{1}{2}\sqrt{\frac{\pi }{a}}\left[{\mathrm{e}}^{+2\sqrt{ab}}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}+\sqrt{ax}\right)-{\mathrm{e}}^{-2\sqrt{ab}}\enspace \text{erf}\left(\sqrt{\frac{b}{x}}-\sqrt{ax}\right)\right].\end{equation} \tag{ B.10 }$$

The two integrals appearing in equations (45) and (46) can be developed by evaluating

$$\begin{equation}I=\int \sqrt{x}{\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\mathrm{d}x.\end{equation} \tag{ B.11 }$$

Since

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}x}\left(\sqrt{x}\enspace {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\right)={\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\left(-a{x}^{1/2}+\frac{b}{{x}^{3/2}}+\frac{1}{2{x}^{1/2}}\right),\end{equation} \tag{ B.12 }$$

we get the relation

$$\begin{equation}\sqrt{x}\enspace {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}=\int {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\left(-a{x}^{1/2}+\frac{b}{{x}^{3/2}}+\frac{1}{2{x}^{1/2}}\right)\mathrm{d}x.\end{equation} \tag{ B.13 }$$

Therefore, we easily obtain

$$\begin{equation}I+\frac{\sqrt{x}}{a}\enspace {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}=I+\int {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}\left(-{x}^{1/2}+\frac{b}{a{x}^{3/2}}+\frac{1}{2a{x}^{1/2}}\right)\mathrm{d}x=\frac{b}{a}{I}_{1}+\frac{1}{2a}{I}_{2}.\end{equation} \tag{ B.14 }$$

This expression allows to write the final result in the form

$$\begin{align}I& =\frac{\sqrt{\pi }}{2a}\left[{\mathrm{e}}^{+2\sqrt{ab}}\left(\frac{1}{2\sqrt{a}}-\sqrt{b}\right)\text{erf}\left(\sqrt{\frac{b}{x}}+\sqrt{ax}\right)\right.\\ & \quad \left.-{\mathrm{e}}^{-2\sqrt{ab}}\left(\frac{1}{2\sqrt{a}}+\sqrt{b}\right)\text{erf}\left(\sqrt{\frac{b}{x}}-\sqrt{ax}\right)\right]-\frac{\sqrt{x}}{a}\enspace {\mathrm{e}}^{-ax}\enspace {\mathrm{e}}^{-\frac{b}{x}}.\end{align} \tag{ B.15 }$$