# Role of temperature in the decohesion of an elastic chain tethered to a substrate by onsite breakable links

Giuseppe Florio, Giuseppe Puglisi, Stefano Giordano

► To cite this version:

Giuseppe Florio, Giuseppe Puglisi, Stefano Giordano. Role of temperature in the decohesion of an elastic chain tethered to a substrate by onsite breakable links. Physical Review Research, 2020, 2 (3), 10.1103/PhysRevResearch.2.033227. [https://doi.org/10.1103/PhysRevResearch.2.033227](https://doi.org/10.1103/PhysRevResearch.2.033227)

HAL Id: hal-02947920

https://hal.science/hal-02947920

Submitted on 24 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the role of temperature in the decohesion of an elastic chain tethered to a substrate by on-site breakable links

Giuseppe Florio*

Politecnico di Bari, (DMMM) Dipartimento di Meccanica, Matematica e Management, Via Re David 200, I-70125 Bari, Italy and INFN, Sezione di Bari, I-70126, Italy

Giuseppe Puglisi†

Politecnico di Bari, (DICAR) Dipartimento di Scienza dell'Ingegneria Civile e dell'Architettura, Via Re David 200, I-70125 Bari, Italy

Stefano Giordano‡

Institute of Electronics, Microelectronics and Nanotechnology - UMR 8520, Univ. Lille, CNRS, Centrale Lille, ISEN, Univ. Valenciennes, LIA LICS/LEMAC, F-59000 Lille, France

(Dated: July 23, 2020)

We analyze the role of temperature in the rate-independent cohesion and decohesion behavior of an elastic film, mimicked by a one-dimensional mass-spring chain, grounded to an undeformable substrate via a one-dimensional sequence of breakable links. In the framework of equilibrium statistical mechanics, in both isometric (Helmholtz ensemble) and isotensional (Gibbs ensemble) conditions, we prove that the decohesion process can be described as a transition at a load threshold, sensibly depending on temperature. Under the classical assumption of having a single domain wall between attached and detached links (zipper model), we are able to obtain analytical expressions for the temperature dependent decohesion force, qualitatively reproducing important experimental effects in biological adhesion. Interestingly, although the two ensembles exhibit a similar critical behavior, they are not equivalent in the thermodynamic limit since they display dissimilar force-extension curves and, in particular, significantly different decohesion thresholds.

## I. INTRODUCTION

Mechanical systems with micro-instabilities attracted a great deal of interest in recent years. In particular, the relationship between microstructural instabilities and the macroscopic (homogenized) response represents a fundamental issue in describing the behavior of complex materials and structures. For instance, the relation between the energetic favorable transition strategy and the corresponding behavior of macroscopic observables represents a crucial aspect in the modeling of such multistable systems. As a matter of fact, micro-instabilities are at the origin of several important phenomena observed in many biological materials and artificial structures. Examples of the former case are the conformational transitions in polymeric and biopolymeric chains [1–14], the attachment and detachment of fibrils in cell adhesion [15–20], the unzipping of macromolecular hairpins [21–26], the sarcomeres behavior in skeletal muscles [27–33] and the denaturation or degradation of nucleic acids, polypep-tidic chains or other macromolecules of biological origin [34–43]. On the other hand, with regard to artificial systems, we recall the peeling of a film from a substrate [44–50], the waves propagation in bistable lattices [51–55], the energy harvesting through multistable chains [56–58],

the plasticity and the hysteresis in phase transitions and martensitic transformations of solids [59–69].

The common feature in all these examples is that the observed macroscopic material behavior results from the evolution of the system in a complex multi-wells energy landscape associated to different configurations of the system at the microscale. Indeed, in all previous example the micro (molecular) systems are composed by units characterized by two (or more) distinct physical states with different static and dynamic features. As a result, under external (mechanical, thermal or electromechanical) loading, the system can experience a transition between different configurations with resulting variable macroscopic properties.

Schematically, we can identify two main classes of microinstabilities in multistable materials. On one side, we may observe a bistable (or multistable) behavior between one ground state and one (or more) metastable state. These states represent different, yet mechanically resistant conformations. This case can be represented in a one dimensional setting by introducing an effective two-wells potential energy $U$, as schematized in Fig.1a. We observe that the choice of a simple one dimensional system is aimed at a purely analytical investigation of the phenomenon. Thus, our one dimensional two-wells energy can be thought to be deduced based on a coarse-graining approach, applied to the real multi-variables structure of a unit. In other words, $x$ represents an effective order parameter adopted to describe the transition between the different wells, whereas all other variables

<footnote>* giuseppe.florio@poliba.it</footnote>

<footnote>† giuseppe.puglisi@poliba.it</footnote>

<footnote>‡ stefano.giordano@univ-lille.fr</footnote>
can be considered to be minimized out. Of course, this
is a simplification useful to provide an easier physical
interpretation of the underlying complex full scale phe-
nomenon. For instance, in this class can be inscribed
conformational (folded → unfolded) transitions in poly-
mers or protein macromolecules (e.g., β-sheets domain
unfolding, α-helix to β-sheet transitions) or martensitic
phase changes between different configurations in metal-
lic alloys. The other class corresponds to transitions be-
tween unbroken and broken states of breakable units of
the system. This process can be reversible, partly re-
versible or irreversible according to the specific physical
phenomenon. Examples of this scheme include unzip-
ping of hairpins, denaturation of macromolecules, fibril-
lar biological adhesion, cell adhesion, and peeling of films.
The (coarse-grained) one dimensional energy considered
in this second case is shown in Fig.1b. Here, the un-
broken configuration corresponds to a potential well and
the broken configuration corresponds to constant energy
and zero force. In Ref.[70], the authors suggest that this
second case can be deduced as a limit case of the first
one when the natural configuration of the second state is
degenerate. Of course the assumption of a non-coercive
energy density changes the mathematical structure of the
model. In particular, as we show in the following (see the
discussion in Sect. III), we solve the possible integrability
problems of the partition function by excluding, thanks
to the considered boundary conditions, the fully detached
configuration.

The statistical mechanics of such systems can be stud-
ied by means of the spin variables approach. The first
models based on this technique have been developed for
describing the response of skeletal muscles [27, 28]. More
recently, this approach has been generalized to study
different allosteric systems [30–33], and macromolecular
chains [71–76]. The main idea consists in introducing a
series of discrete (spin) variables to identify the state of
the system units. For example, the bistable potential en-
ergy of Fig.1a (continuous line) can be approximated by
a bi-parabolic function (dashed lines) with the switch-
ing among the wells described by the spin variable. This
strongly facilitates the calculation of the partition func-
tion and the corresponding thermodynamic quantities,
delivering possible analytical results with a correspond-
ing clear physical interpretation. It is important to re-
mark that, as discussed in detail in Refs. [71, 74], the
evaluation of the partition function based on the spin
approach assumes that for both configurations all possi-
ble deformations (values of x in Fig.1) can be attained
by the system. This corresponds to the assumption of a
multivalued energy function (see superposition of dashed
lines in Fig.1). As shown numerically in Refs. [71, 74],
with typical experimental temperatures, the effect of this
approximation can be considered (statistically) negligi-
ble since these artificial configurations (superposition of
dashed curves) have an energy sensibly higher than real
configurations (continuous lines).

While the use of spin variables has been largely

<figure><img src="image_2.png" /><figcaption>FIG. 1. Two different classes of micro-instabilities: bistability between a folded (ground) state and an unfolded (metastable) state in panel (a), and link transition between attached and detached state in panel (b).</figcaption></figure>

adopted to model units with transitions between ground
and metastable states (see Fig.1a), the case of units un-
dergoing transition between unbroken and broken states
(Fig.1b) has been investigated, to the best of our knowl-
edge, only in the cases of *parallel* links [77–79]. In the
case considered here, the film deposited on a given sub-
strate is represented as a lattice of masses connected *in
series* by harmonic springs and linked to a substrate by
breakable links. This assumption may reproduce differ-
ent types of adhesion forces as in the case of biological
and cellular adhesion [47, 80]. Furthermore, such type
of systems has been recently applied to study thermal
effects in the mechanical denaturation of DNA [81].

It is important to observe that in the recalled exam-
ples of biological adhesion and decohesion phenomena the
binding enthalpies range from about one k<sub>B</sub>T (for hydro-
gen bonds, e.g. controlling the stability of base pairs in
DNA) to tens of k<sub>B</sub>T (for covalent bonds, e.g. linking
neighboring bases in a DNA strand or protein structures)
[82, 83]. As a result, the role of temperature cannot be
neglected so that the framework of statistical mechanics
here proposed represents a proper theoretical setting [84].

Another important effect considered in this paper, is
the possibility of two types of loading. Indeed, the peel-
ing of the film can be induced by prescribing a given
extension (Helmholtz ensemble, hard device, Fig.2c), or
by applying an external force to the last unit of the chain
<figure><img src="image_3.png" /><figcaption>FIG. 2. Energies and force-displacement relation of the harmonic (horizontal) shear springs in panel (a), and breakable (vertical) links in panel (b). Different loading conditions: assigned force (Gibbs ensemble, soft device) in panel (c), and assigned end displacement (Helmholtz ensemble, hard device) in panel (d).</figcaption></figure>

(Gibbs ensemble, soft device, Fig.2d). More specifically, as described in Refs.[74, 75], these two boundary conditions can be deduced as limiting conditions of real loading experiments, when the stiffness of the device is large (hard device) or is negligible (soft device), respectively. From a theoretical point of view, one important problem regards the analysis of the equivalence of the two ensembles in the thermodynamic limit (i.e. for very large systems) [85–91].

Here, based on the proposed spin variables approach, we are able to deduce a fully analytical solution of the two boundary problems for an arbitrary number $N$ of elements of the chain and also in the thermodynamic limit ($N \to \infty$). The main result of this paper is the deduction of analytic expressions of the temperature dependent debonding force in both cases of isotensional and isoextensional loading. Interestingly we obtain a new effect as compared with the other case of bistable unbreakable elements (with units as in Fig.1a). Indeed in the case of non-degenerate wells (Fig.1a) the transition between the two homogeneous (folded↔unfolded) states takes place at a temperature independent force [71–75] (see also Ref. [92] for the case when the two wells have different elastic constants). More specifically, the conformational transition corresponds to a force-displacement diagram with a temperature dependent slope, but a temperature independent average force (corresponding to the Maxwell force of the two wells energy). This result

can be schematically interpreted in the framework of the Bell relation $f = \Delta E/\Delta x$ (see Fig.1a for the definition of $\Delta E$ and $\Delta x$), discovered in the context of cell adhesion [15, 16, 89], where this value of the average force represents the threshold necessary to make the unfolding rate equal to the (reverse) folding one. This threshold force can be explained as follows. We consider two potential energies $U_1(x) = \frac{1}{2}k(x-\ell_1)^2 - fx$ and $U_2(x) = \Delta E + \frac{1}{2}k(x-\ell_2)^2 - fx$, corresponding to the wells of the units. The equilibrium positions can be obtained by $\partial U_i/\partial x = 0$ and we get $x_i = \ell_i + f/k$ ($i = 1, 2$). Hence, the unfolded configuration is more favorable of the folded one when $U_2(x_2) < U_1(x_1)$, which corresponds to $f < \Delta E/\Delta x$ where $\Delta x = \ell_2 - \ell_1$. Differently, in the case of breakable links of interest in this paper (with energy function as in Fig.1b), we obtain an unusual temperature dependent force threshold. More specifically, in this case not only the slope of the force-displacement diagram grows with temperature, but importantly the average decohesion force decreases as the temperature increases. Remarkably, this decohesion force becomes zero for a given critical temperature. This result describes the expected effect that thermal fluctuations may anticipate the decohesion, allowing the escape from the energy well and, therefore, the exploration of the whole configurations in Fig.1b. Here, we analytically describe this effect and the whole process is explained by means of a phase transition occurring at a given critical temperature.
ture. In particular, we obtain that the system is able to undergo a complete decohesion even without any external mechanical action. As we show, the decohesion force threshold may be significantly different in the Helmholtz and Gibbs ensembles also in the thermodynamical limit. We also remark that a similar effect of a temperature dependent denaturation force is observed in DNA both experimentally and theoretically [81, 93].

Summarizing, we underline that from the statistical mechanics point of view, the system here investigated is particularly interesting for three reasons: (i) it can be analytically solved within both statistical ensembles; (ii) it shows a phase transition at a critical temperature that can be calculated in closed form; (iii) it exhibits the ensembles non-equivalence in the thermodynamic limit, which is an unusual and intriguing behavior.

The paper is organized as it follows. In Section II, we introduce the system. In Sections III and IV, we study the Helmholtz and Gibbs ensembles, respectively. In Section V, we study the thermodynamic limit for both ensembles. The conclusions (Section VI) and a mathematical Appendix close the paper.

## II. PROBLEM STATEMENT

We schematize the decohesion of a layer from a substrate by considering a one-dimensional chain of elements embedded in an on-site potential reproducing the behavior of detachable links. This type of models and its continuous version describing reversible decohesion has been previously introduced to describe a wide range of phenomena such as peeling of tapes, adhesion of geckos and denaturation of DNA or other chemical structures [47, 48, 94]. However, in those cases thermal effects are typically neglected. On the contrary, as we already stated in the introduction, these effects may play a central role in the decohesion behavior of biological materials. It is important to remark that our paradigmatic system has the advantage of analytical simplicity letting a clear physical interpretation. On the other hand, it can be generalized in various ways even though this would need a numerical treatment. For example, the assumption of biparabolic energy can be extended by adding a persistence length to the elastic chain (changing from flexible to semi-flexible), a more realistic two-dimensional lattice tethered to the substrate can be considered, a deformable substrate can be introduced, and so forth. Of course, each new ingredient can modify the system response and, in particular, can complicate the deduction of its critical behavior.

The horizontal springs of the one-dimensional lattice (elastic constant $k$) are purely harmonic with potential energy $\varphi = \frac{1}{2}k(y_{i+1} - y_i)^2$ (Fig.2a), while the vertical ones (elastic constant $h$) can be broken or unbroken depending on their extension $y_i$ (Fig.2b). When $|y_i| > y_M$ they are broken and when $|y_i| < y_M$ they are unbroken. Therefore, an unbroken spring leads to a contribution

to the potential energy equal to $\psi = \frac{1}{2}hy_i^2$ (when $|y_i| < y_M$) and a broken one a contribution equal to $\psi = \frac{1}{2}hy_M^2$ (when $|y_i| > y_M$). As anticipated, two different loading conditions are considered. In the first case (Fig.2c), the process is controlled by the prescribed position $y_{N+1} = y_d$ of the last element of the chain (isometric condition within the Helmholtz ensemble). In the second case (Fig.2d), the process is controlled by the applied force $f$ (isotensional condition within the Gibbs ensemble).

The most important point, on which is grounded our approach, is that each vertical element is characterized by two different states (broken and unbroken configurations). Therefore, we associate each unit with a spin variable and the energy potential of each vertical spring can be written as

$$ \psi = \frac{1}{4}(1 + s_i)hy_i^2 + \frac{1}{4}(1 - s_i)hy_M^2, \quad (1) $$

where $s_i = +1$ corresponds to the unbroken 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$ and the $N$ discrete variables $s_i$. The switching of the variable $s_i$ and their statistics at thermodynamic equilibrium are directly controlled by the statistical ensemble (Helmholtz and Gibbs in our case) imposed to the system.

The statistical mechanics analysis of this system cannot be analytical and it is computational expensive. Nevertheless, since we are studying the cohesion-decohesion process under an external mechanical action applied to one end point of the system, we can simplify the model (see Fig.2d) by assuming to have $N-\xi$ broken elements on the right of the chain and $\xi$ unbroken elements on the left of the chain. In other words we suppose to have a single moving interface or domain wall between the attached region and detached region. This is a plausible hypothesis, especially if we work at not too large temperature values and not too low values of the force. Indeed, in the zero-temperature case, this configuration is the only energy minimizer of the system [47, 48]. Moreover this hypothesis coincides with 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 [95-98]. The analysis of other regimes, for high values of the temperature is the subject of a forthcoming paper and it is out of the aim of this one [99].

Under these hypotheses, the set of the two-state spin variables $s_i$ is substituted by the single variable $\xi$ belonging to the phase space of the system, and taking its values in the set $\{0, 1, 2, ..., N\}$. In this regard, the variable $\xi$ can be considered as a multivalued spin variable.

The aim of this work is to fully analyze the cohesion-decohesion process 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.
### III. HARD DEVICE: HELMHOLTZ ENSEMBLE

Consider first the case of a prescribed extension $y_{N+1} = y_d$ of the last element of the chain, as represented in Fig.2c (isometric condition). As previously anticipated, the variables belonging to the phase space of this system are the extensions $y_i$ of the vertical springs ($i = 1, ..., N$), and the number $\xi$ of unbroken links. The total potential energy is

$$ \Phi = \sum_{i=0}^{N} \frac{1}{2} k (y_{i+1} - y_i)^2 + \sum_{i=1}^{\xi} \frac{1}{2} h y_i^2 + \sum_{i=1+\xi}^{N} \frac{1}{2} h y_M^2, \quad (2) $$

where $y_0 = 0$ and $y_{N+1} = y_d$ in the case of assigned displacement is imposed. It is worth noticing that the last term in $\Phi$ is not an irrelevant additive constant since it depends implicitly on $\xi$, which is a variable of the phase space of the system. The assumption of $y_0 = 0$ reproduces typical boundary conditions of experimental cases. On the other hand, as can be deduced by the following analysis, such an assumption also solves the possible integrability problems concerning the partition function calculation, that can come from the non coercivity of the potential energy of the breakable links. We also observe that in an analogy with the classical Griffith approach to fracture mechanics [100], here we model the decohesion behavior of the layer by considering the energetic competition between the elastic energy of the chain and the unbinding energetic contribution (fracture energy).

The energy function $\Phi$ can be rearranged by means of the following matrix definition

$$ \mathcal{A}(\xi) = \begin{bmatrix} \mathcal{A}_{11} & \mathcal{A}_{12} \\ \mathcal{A}_{21} & \mathcal{A}_{22} \end{bmatrix} \in \mathcal{M}_{N,N}(\mathbb{R}), \quad (3) $$

which is based on these four sub-matrices

$$ \mathcal{A}_{11} = \begin{bmatrix} 2+\eta & -1 & 0 & \cdots \\ -1 & 2+\eta & -1 & \cdots \\ 0 & -1 & 2+\eta & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \in \mathcal{M}_{\xi,\xi}(\mathbb{R}), \quad (4) $$

$$ \mathcal{A}_{22} = \begin{bmatrix} 2 & -1 & 0 & \cdots \\ -1 & 2 & -1 & \cdots \\ 0 & -1 & 2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \in \mathcal{M}_{N-\xi, N-\xi}(\mathbb{R}), \quad (5) $$

$$ \mathcal{A}_{12} = \begin{bmatrix} 0 & \cdots & 0 & 0 \\ 0 & \cdots & 0 & 0 \\ 0 & \ddots & \vdots & \vdots \\ -1 & 0 & 0 & 0 \end{bmatrix} \in \mathcal{M}_{N, N-\xi}(\mathbb{R}), \quad (6) $$

$$ \mathcal{A}_{21} = \begin{bmatrix} 0 & 0 & 0 & -1 \\ \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \end{bmatrix} \in \mathcal{M}_{N-\xi, N}(\mathbb{R}), \quad (7) $$

where

$$ \eta = \frac{h}{k} \quad (8) $$

is the main non dimensional parameter of the paper, representing the ratio between the elastic constants of vertical and horizontal elastic elements. Moreover we introduce the vectors $\mathbf{v} = (0, 0, 0, ..., 0, 1) \in \mathbb{R}^N$ and $\mathbf{y} = (y_1, y_2, y_3, ..., y_N) \in \mathbb{R}^N$. The energy function can then be rewritten as follows

$$ \Phi(\mathbf{y}, \xi; y_d) = \frac{1}{2} k \mathcal{A}(\xi) \mathbf{y} \cdot \mathbf{y} - k y_d \mathbf{v} \cdot \mathbf{y} + \frac{1}{2} k y_d^2 + \frac{1}{2} k \eta (N - \xi) y_M^2, \quad (9) $$

where $\mathbf{y}$ and $\xi$ are the main variables belonging to the phase space of the system.

This expression of $\Phi$ is useful in the following Gaussian integration for the partition function since it is constituted by the sum of a quadratic form and a linear form in $\mathbf{y}$, with an additional term independent of $\mathbf{y}$. The partition function of the system analyzed within the Helmholtz ensemble can therefore be written as

$$ Z_H(\beta, y_d) = \sum_{\xi=0}^{N} \int_{\mathbb{R}^N} e^{-\beta \Phi(\mathbf{y}, \xi; y_d)} d\mathbf{y}. \quad (10) $$

where $\beta = (k_B T)^{-1}$, $k_B$ is the Boltzmann constant and $T$ the absolute temperature. When Eq.(9) is substituted in Eq.(10), we get

$$ Z_H(\beta, y_d) = e^{-\beta \frac{k}{2} y_d^2} \sum_{\xi=0}^{N} e^{-\beta \frac{k \eta}{2} (N-\xi) y_M^2} \times \int_{\mathbb{R}^N} e^{-\beta \frac{k}{2} \mathcal{A}(\xi) \mathbf{y} \cdot \mathbf{y}} e^{\beta k y_d \mathbf{v} \cdot \mathbf{y}} d\mathbf{y}, \quad (11) $$

where we can use the property of the Gaussian integrals

$$ \int_{\mathbb{R}^N} e^{-\frac{1}{2} \mathcal{M} \mathbf{y} \cdot \mathbf{y} + \mathbf{b} \cdot \mathbf{y}} d\mathbf{y} = \sqrt{\frac{(2\pi)^N}{\det \mathcal{M}}} e^{\frac{1}{2} \mathcal{M}^{-1} \mathbf{b} \cdot \mathbf{b}}, \quad (12) $$

holding for any symmetric and positive definite matrix $\mathcal{M}$. Indeed, by considering $\mathcal{M} = \beta k \mathcal{A}(\xi)$ and $\mathbf{b} = \beta k y_d \mathbf{v}$ in Eq.(12), we easily obtain

$$ \begin{aligned} Z_H(\beta, y_d) &= \left(\frac{2\pi}{\beta k}\right)^{N/2} \sum_{\xi=0}^{N} e^{-\beta \frac{k y_M^2}{2} \eta (N-\xi)} \\ & \times \frac{1}{\sqrt{\det \mathcal{A}(\xi)}} e^{-\beta \frac{k y_d^2}{2} \{1-\mathcal{A}^{-1}(\xi) \mathbf{v} \cdot \mathbf{v}\}} \\ &= \left(\frac{2\pi}{\beta k}\right)^{N/2} e^{-N\beta \frac{k y_M^2}{2} \eta} \sum_{\xi=0}^{N} \Gamma_{\xi}(\beta, y_d), \quad (13) \end{aligned} $$

where

$$ \Gamma_{\xi}(\beta, y_d) = \frac{e^{\beta \frac{k y_M^2}{2} \eta \xi}}{\sqrt{\det \mathcal{A}(\xi)}} e^{-\beta \frac{k y_d^2}{2} \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}} \quad (14) $$

and we have used the definition of $\mathbf{v}$.

The knowledge of the partition function allows us to determine the expectation value of the force conjugated to the assigned displacement $y_d$ [85]

$$ \langle f \rangle_H = -\frac{1}{\beta} \frac{\partial \ln Z_H}{\partial y_d} = -\frac{1}{\beta} \frac{1}{Z_H} \frac{\partial Z_H}{\partial y_d}, \quad (15) $$
which can be rewritten as

$$ \langle f \rangle_H = \frac{\sum_{\xi=0}^{N} \{1 - A_{NN}^{-1}(\xi)\} \Gamma_{\xi}(\beta, y_d)}{\sum_{\xi=0}^{N} \Gamma_{\xi}(\beta, y_d)} k y_d \\ = \left( 1 - \frac{\sum_{\xi=0}^{N} A_{NN}^{-1}(\xi) \Gamma_{\xi}(\beta, y_d)}{\sum_{\xi=0}^{N} \Gamma_{\xi}(\beta, y_d)} \right) k y_d. \quad (16) $$

Another important quantity to describe the decohesion dependence from the temperature is the average value $\langle \xi \rangle_H$ of unbroken vertical springs. It can be directly evaluated through the expression

$$ \langle \xi \rangle_H = \frac{1}{Z_H} \sum_{\xi=0}^{N} \int_{\mathbb{R}^N} \xi e^{-\beta \Phi(y, \xi; y_d)} dy = \frac{\sum_{\xi=0}^{N} \xi \Gamma_{\xi}(\beta, y_d)}{\sum_{\xi=0}^{N} \Gamma_{\xi}(\beta, y_d)}. \quad (17) $$

Moreover we can determine the expectation value of the whole position vector $\mathbf{y}$ (assigning the deformed configuration of the chain at given displacement):

$$ \langle \mathbf{y} \rangle_H = \frac{1}{Z_H} \sum_{\xi=0}^{N} \int_{\mathbb{R}^N} \mathbf{y} e^{-\beta \Phi(\mathbf{y}, \xi; y_d)} d\mathbf{y}. \quad (18) $$

Then, we can use the definition of $\Phi$, thus obtaining the more explicit expression

$$ \langle \mathbf{y} \rangle_H = e^{-\beta \frac{k}{2} y_d^2} \frac{1}{Z_H} \sum_{\xi=0}^{N} e^{-\beta \frac{k\eta}{2} (N-\xi) y_M^2} \quad (19) \\ \times \int_{\mathbb{R}^N} e^{-\beta \frac{k}{2} A(\xi) \mathbf{y} \cdot \mathbf{y}} e^{\beta k y_d \mathbf{v} \cdot \mathbf{y}} \mathbf{y} d\mathbf{y}. $$

Since, by differentiating (12) with respect to $\mathbf{b}$ we get

$$ \int_{\mathbb{R}^N} e^{-\frac{1}{2} \mathcal{A} \mathbf{y} \cdot \mathbf{y} + \mathbf{b} \cdot \mathbf{y}} \mathbf{y} d\mathbf{y} = \sqrt{\frac{(2\pi)^N}{\det \mathcal{A}}} e^{\frac{1}{2} \mathcal{M}^{-1} \mathbf{b} \cdot \mathbf{b}} \mathcal{A}^{-1} \mathbf{b}, \quad (20) $$

straightforward calculations give

$$ \langle \mathbf{y} \rangle_H = \frac{\sum_{\xi=0}^{N} \{\mathcal{A}^{-1}(\xi)\mathbf{v}\} \Gamma_{\xi}(\beta, y_d)}{\sum_{\xi=0}^{N} \Gamma_{\xi}(\beta, y_d)} y_d. \quad (21) $$

In Eqs.(32)-(35) we report the fully analytical expressions of the expectation values of the force, debondend fraction and displacement vector obtained by using the explicit formulas in Eqs.(96)-(98) of the quantities $\det \mathcal{A}(\xi)$, $1 - A_{NN}^{-1}(\xi)$, $\{\mathcal{A}^{-1}(\xi)\mathbf{v}\}_i$ reported in the Appendix.

In Fig.3 we illustrate the obtained behavior for a system under isometric loading for a 'short' chain with only $N = 6$ elements and for a 'longer' chain with $N = 30$ elements. The case with $N = 6$ [left panels (a) and (c) of Fig.3] shows the importance of discreteness, exhibiting distinct rupture occurrences in the decreasing steps of the quantity $\langle \xi \rangle_H$, and in the peaks of the force-displacement $\langle (f)_H, y_d \rangle$ curves along the whole decohesion process. The resulting evolution of the displacements is also reported in Fig.4b. As expected, for large values of the temperature $T$ the debonding process is less "localized",

with the system exploring more configurations. As a result, the curves are smoother and it is more difficult to recognize the single ruptures.

Similar results are obtained also if we change the adhesion energy (see Fig.4a). Indeed, as $h$ is increased the detachment process is more localized and the peak decreases with the system passing from a peeling type decohesion to a pull-off type decohesion, as observed in the case of athermal decohesion [47, 48].

**Remark.** In passing, we recall that we assume, both in the hard and soft device, $y_0 = 0$, thus fixing the left end point. Of course this choice can be varied according with the specific experimental phenomenon to be reproduced. As a result, in these and later force-displacement curves, after the force plateaux, when the system is fully detached, we have a new elastic branch. In the opposite assumption of free left end point, such an effect is not observed and the force decreases to zero when the full detachment is attained. However, as discussed above, if the system is in a non-connected configuration (e.g., completely detached), the partition function corresponds to a non-convergent integral, generating some technical difficulty. Of course, it is possible to find some ways to introduce the rupture but we preferred to follow a simpler procedure for the sake of simplicity and clarity.

The case with a larger number of elements ($N = 30$) is represented in the right panels (b) and (d) of Fig.3. Observe that the system is initially characterized by load oscillations corresponding to the first debonding effects, but the diagram rapidly converges to a constant force plateau as the debonding front propagates far from the loading end point. The increasing of temperature cancels out also this initial discreteness effect.

The main interesting feature, previously anticipated, is the observation of the temperature dependent unfolding plateau, with a detachment force threshold sensibly decreasing as the temperature grows. This behavior, which marks the strong difference of the bistable elements of the type reported in Fig.1a, with respect with the breakable links considered here, is thoroughly studied in Section VA where we consider the thermodynamic limit $N \to \infty$.

We remark that the gradual increase of the debonded fraction $\langle \xi \rangle_H$ obtained in the hard device is in agreement with the behavior observed in the unfolding of proteins or other macromolecules under isometric conditions [71, 74, 89]. Indeed, the force spectroscopy of protein chains under isometric conditions is characterized by a sawtooth-like force-extension response showing that the domains unfold progressively in reaction to the prescribed increasing extension [4] (similarly to the response seen in Fig.3, at least for small values of $N$). A completely different behavior is observed in the Gibbs (isotential) ensemble [101], studied in the following section, where a much more cooperative behavior is established and the vertical elements break quite simultaneously, at a given critical force.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Average force $\langle f \rangle_H$ in top panels (a) and (b), and average number of unbroken elements $\langle \xi \rangle_H$ in bottom panels (c) and (d), versus $y_d$ for a decohesion process ($y_d$ increasing) under isometric conditions (Helmholtz ensemble). We adopted the parameters $N = 6$ and $y_M = 4$ for panels (a) and (c), and $N = 30$ and $y_M = 2$ for panels (b) and (d). Other parameters are common: $k = 5$, $h = 20$, and six values of $\beta^{-1} = k_B T = 4, 7.2, 10.4, 13.6, 16.8$, and 20 (in arbitrary units).</figcaption></figure>

To conclude, we also observe that the behavior of the force-extension curves, characterized by a force plateau, is in good qualitative agreement with several results obtained in peeling experiments and simulations concerning adhesive films and bioclinical structures [102–106].

## IV. SOFT DEVICE: GIBBS ENSEMBLE

Consider now the case when the film is loaded by a fixed force $f$ (Fig.2b, isotensional condition). In this case we introduce the Gibbs ensemble. The total energy of the system can now be written as $\Phi - f y_{N+1}$, where $\Phi$ is the energy function introduced within the Helmholtz ensemble in Eq.(9). Thus, the Gibbs partition function

can be written as

$$ 
\begin{aligned}
Z_G (f) &= \int_{\mathbb{R}^{N+1}} \sum_{\xi=0}^{N} e^{-\beta \Phi(y, \xi; y_{N+1})} e^{\beta f y_{N+1}} dy dy_{N+1} \\
&= \int_{-\infty}^{+\infty} Z_H (y_{N+1}) e^{\beta f y_{N+1}} dy_{N+1}.
\end{aligned} 
\quad (22) $$

Of course this corresponds to the Laplace transform of the Helmholtz partition function [85]. By using Eqs.(13)-(14), we get

$$ Z_G (f) = \left(\frac{2\pi}{\beta k}\right)^{N/2} e^{-N\beta \frac{k y_M^2}{2}} \eta \\ \times \sum_{\xi=0}^{N} \int_{-\infty}^{+\infty} \Gamma_{\xi}(\beta, y_{N+1}) e^{\beta f y_{N+1}} dy_{N+1}, \quad (23) $$

where the integral is Gaussian due to the structure of $\Gamma_{\xi}$ and, therefore, can be easily evaluated. A straightfor-
<figure><img src="image_5.png" /><figcaption>FIG. 4. Panel (a): Variable decohesion behavior at different value of adhesion energy ($h = 1 - 30, \eta = 0.2 \to 6$). Here we assumed $k = 5, N = 30, k_B T = 3$ and $y_M = 1.2$. Panel (b): Average displacements evolution under isometric loading for the same chain considered in the left panels (a) and (c) of Fig.3.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

ward calculation leads to

$$ 
\begin{aligned}
Z_G(f) &= \left(\frac{2\pi}{\beta k}\right)^{\frac{N+1}{2}} \sum_{\xi=0}^{N} \frac{1}{\sqrt{\det \mathcal{A}(\xi) \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}}} \\
&\quad \times e^{-\frac{\beta k y_M^2}{2}} \eta(N-\xi) e^{\beta \frac{f^2}{2k} \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}^{-1}} \\
&= \left(\frac{2\pi}{\beta k}\right)^{\frac{N+1}{2}} e^{-N \frac{\beta k y_M^2}{2}} \eta \sum_{\xi=0}^{N} \Omega_{\xi}(\beta, f),
\end{aligned} 
\quad (24) $$

where

$$ \Omega_{\xi}(\beta, f) = \frac{e^{\beta \frac{k y_M^2}{2} \eta \xi} e^{\beta \frac{f^2}{2k} \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}^{-1}}}{\sqrt{\det \mathcal{A}(\xi) \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}}}. \quad (25) $$

The expected value of the extension $\langle y_d \rangle = \langle y_{N+1} \rangle$ of the last element of at given applied force $f$ can then be evaluated as [85]

$$ \langle y_d \rangle_G = \frac{1}{\beta} \frac{\partial \ln Z_G}{\partial f} = \frac{1}{\beta} \frac{1}{Z_G} \frac{\partial Z_G}{\partial f}, \quad (26) $$

that gives

$$ \langle y_d \rangle_G = \frac{\sum_{\xi=0}^{N} \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}^{-1} \Omega_{\xi}(\beta, f)}{\sum_{\xi=0}^{N} \Omega_{\xi}(\beta, f)} \frac{f}{k}. \quad (27) $$

Similarly, we can also determine the average number of unbroken links

$$ 
\begin{aligned}
\langle \xi \rangle_G &= \frac{1}{Z_G} \int_{\mathbb{R}^{N+1}} \sum_{\xi=0}^{N} \xi e^{-\beta \Phi(\mathbf{y}, \xi; y_{N+1})} e^{\beta f y_{N+1}} dy dy_{N+1} \\
&= \frac{\sum_{\xi=0}^{N} \xi \Omega_{\xi}(\beta, f)}{\sum_{\xi=0}^{N} \Omega_{\xi}(\beta, f)}.
\end{aligned} 
\quad (28) $$

Eventually, we can evaluate the average value of the displacement of each elements of the chain. This quantity is defined by the expression

$$ \langle \mathbf{y} \rangle_G = \frac{1}{Z_G} \sum_{\xi=0}^{N} \int_{\mathbb{R}^{N+1}} \mathbf{y} e^{-\beta \Phi(\mathbf{y}, \xi; y_{N+1})} e^{\beta f y_{N+1}} dy dy_{N+1}. \quad (29) $$

The comparison of Eq.(29) with Eqs.(18) and (14) yields

$$ 
\begin{aligned}
\langle \mathbf{y} \rangle_G = \left(\frac{2\pi}{\beta k}\right)^{N/2} \frac{e^{-N\beta \frac{k y_M^2}{2} \eta}}{Z_G} \sum_{\xi=0}^{N} \int_{-\infty}^{+\infty} y_{N+1} \Gamma_{\xi}(\beta, y_{N+1}) \\
\times e^{\beta f y_{N+1}} \mathcal{A}^{-1}(\xi) \mathbf{v} dy_{N+1},
\end{aligned} 
\quad (30) $$

where the final integral is Gaussian due to the structure of $\Gamma_{\xi}(\beta, y_{N+1})$. The calculation can be performed to give the final result in the form

$$ \langle \mathbf{y} \rangle_G = \frac{\sum_{\xi=0}^{N} \{1 - \mathcal{A}_{NN}^{-1}(\xi)\}^{-1} \mathcal{A}^{-1}(\xi) \mathbf{v} \Omega_{\xi}(\beta, f)}{\sum_{\xi=0}^{N} \Omega_{\xi}(\beta, f)}. \quad (31) $$

By using again Eqs.(96)-(98), obtained in the Appendix, we obtain the analytic expressions given in Eqs.(58)-(61).

The obtained results for the isotensional loading are illustrated in Fig.5 for the same chains considered in Fig.3, where the case of isometric loading is described. Important differences between the Helmholtz and the Gibbs responses can be recognized. In particular, the analysis of the evolution of $\langle \xi \rangle_G$ within the Gibbs ensemble shows that the detachment process corresponds to a cooperative breaking of the vertical elements. This result can be compared with the sequential unfolding behavior, observed with an extension controlled decohesion, obtained with the hard device. This dissimilarity between the
<figure><img src="image_7.png" /><figcaption>FIG. 5. Average position $\langle y_d \rangle_G$ in top panels (a) and (b), and average number of unbroken elements $\langle \xi \rangle_G$ in bottom panels (c) and (d), versus $f$ for a decohesion process ($f$ increasing) under isotensional conditions (Gibbs ensemble). We adopted the parameters $N=6$ and $y_M = 4$ for panels (a) and (c), and $N=30$ and $y_M = 2$ for panels (b) and (d). Other parameters are common: $k=5$, $h=20$, and six values of $\beta^{-1} = k_B T = 4, 7.2, 10.4, 13.6, 16.8$, and 20 (in arbitrary units).</figcaption></figure>

Helmholtz and the Gibbs ensembles has been already observed in the unfolding of proteins where the isotensional condition produces a cooperative response whereas the isometric condition generates a non-cooperative response [71, 74, 89].

Another important difference concerns the shape of the force-extension curves measured within the two statistical ensembles. While the isometric case leads to a series of peaks corresponding to the rupture occurrences, the isotensional case is characterized by a monotone force-extension curve. Also this feature can be explained by the quite simultaneous rupture of all the elements observed within the Gibbs ensemble. Again, the simultaneous or cooperative ruptures can be identified in the isotensional behavior of $\langle y \rangle_G$, plotted in Fig.6b, while the sequential or non-cooperative ruptures were observed in the isometric behavior of $\langle y \rangle_H$, plotted in Fig.4b. Also in this case we may observe the significant effect of temperature on the decohesion threshold (particularly evident for large values of $N$). This behavior is thoroughly studied in

Section V B by considering the thermodynamic limit for $N \to \infty$.

In Fig.6a we report also in the case of assigned force the effect of variable adhesion energy. Again by increasing the adhesion energy the system decohesion changes from a peeling type to a pull-off type decohesion [47, 48].

V. THERMODYNAMIC LIMIT

In this Section, we study the behavior of the system, under both isometric and isotensional conditions, for a large chain length (ideally, $N \to \infty$), and we obtain explicit analytic results to describe the system behavior in the thermodynamic limit.
<figure><img src="image_8.png" /><figcaption>FIG. 6. Panel (a): Variable decohesion behavior in the Gibbs ensemble at different value of adhesion energy ($h = 1 \to 30, \eta = 0.2 \to 6$). Here we assumed $k = 5, N = 30, k_B T = 3$ and $y_M = 1.2$. Panel (b): Average displacements evolution in the isotensional loading for the same chain considered in the left panels (a) and (c) of Fig.5.</figcaption></figure>

A. The thermodynamic limit in the Helmholtz ensemble

To analyze the behavior in the thermodynamic limit within the Helmholtz ensemble, we first report the final

expressions for the expected values of the force, for the average number of attached elements, and for displacement vector, as given in Eqs.(16),(17), and (21). After nondimensionalization and the use of Eqs.(96)-(98), we obtain

$$ \langle \mathcal{F} \rangle_H = \beta y_M \langle f \rangle_H = 2 \frac{\tilde{\beta}}{\eta} \frac{\mathcal{Y}}{D} \sum_{\xi=0}^{N} e^{\tilde{\beta}\xi} \frac{\gamma(\xi+1) - \gamma(\xi)}{[(N-\xi+1)\gamma(\xi+1) - (N-\xi)\gamma(\xi)]^{3/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{\gamma(\xi+1)-\gamma(\xi)}{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}} \mathcal{Y}^2, \quad (32) $$

$$ \langle \xi \rangle_H = \frac{1}{D} \sum_{\xi=0}^{N} \xi e^{\tilde{\beta}\xi} \frac{1}{[(N-\xi+1)\gamma(\xi+1) - (N-\xi)\gamma(\xi)]^{1/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{\gamma(\xi+1)-\gamma(\xi)}{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}} \mathcal{Y}^2, \quad (33) $$

$$ \langle \mathcal{Y}_i \rangle_H = \frac{\langle y_i \rangle_H}{y_M} = \frac{\mathcal{Y}}{D} \sum_{\xi=0}^{i-1} e^{\tilde{\beta}\xi} \frac{(i-\xi)\gamma(\xi+1) - (i-\xi-1)\gamma(\xi)}{[(N-\xi+1)\gamma(\xi+1) - (N-\xi)\gamma(\xi)]^{3/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{\gamma(\xi+1)-\gamma(\xi)}{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}} \mathcal{Y}^2 \\ + \frac{\mathcal{Y}}{D} \sum_{\xi=i}^{N} e^{\tilde{\beta}\xi} \frac{\gamma(i)}{[(N-\xi+1)\gamma(\xi+1) - (N-\xi)\gamma(\xi)]^{3/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{\gamma(\xi+1)-\gamma(\xi)}{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}} \mathcal{Y}^2, \quad i=1, \dots, N \quad (34) $$

$$ \mathcal{D} = \sum_{\xi=0}^{N} e^{\tilde{\beta}\xi} \frac{1}{[(N-\xi+1)\gamma(\xi+1) - (N-\xi)\gamma(\xi)]^{1/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{\gamma(\xi+1)-\gamma(\xi)}{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)} \mathcal{Y}^2}, \quad (35) $$

where we introduced the dimensionless parameters

and where

$$ \tilde{\beta} = \frac{h y_M^2}{2k_B T}, \quad \mathcal{Y} = \frac{y_d}{y_M}, \quad (36) $$

$$ \gamma(z) = \frac{1}{\sqrt{\eta^2 + 4\eta}} \left( \frac{2 + \eta + \sqrt{\eta^2 + 4\eta}}{2} \right)^z \\ - \frac{1}{\sqrt{\eta^2 + 4\eta}} \left( \frac{2 + \eta - \sqrt{\eta^2 + 4\eta}}{2} \right)^z, \quad (37) $$

(see Appendix for details).
To begin, we consider Eq.(32), with $\mathcal{D}$ given in Eq.(35). In the limit of large $N$, we introduce the continuum variable $x = \xi/N$ and substitute the summations with integrals

$$ \sum_{\xi=0}^{N} \phi(\xi) \approx \int_{0}^{N} \phi(\xi) d\xi + \frac{\phi(0) + \phi(N)}{2} \\ = N \int_{0}^{1} \phi(Nx) dx + \frac{\phi(0) + \phi(N)}{2}, \quad (38) $$

where we have considered a generic function $\phi$.

*Remark* Observe that this expression corresponds to the trapezoidal rule for approximating a definite integral or, equivalently, to the Euler-Maclaurin formula with only one remainder term [107]. We point out that this more

refined approximation is essential, as proved in the following, to show the important different behavior in the hard and soft device, also in the thermodynamic limit.

After this substitution, we observe that in the new version of Eq.(32) we have several terms of the form $\gamma(Nx)$ and $\gamma(Nx + 1)$. Under the hypothesis of large values of $N$, it is easy to see that the function $\gamma$ in Eq.(37), can be approximated by $\gamma(z) \approx \frac{1}{\sqrt{\eta^2+4\eta}} b^z$, with

$$ b = \frac{2 + \eta + \sqrt{\eta^2 + 4\eta}}{2} > 1. \quad (39) $$

Therefore, after introducing

$$ \Lambda = \tilde{\beta} - \frac{1}{2} \ln b, \quad \rho = \frac{b}{b-1}, \quad (40) $$

we eventually obtain

$$ \langle \mathcal{F} \rangle_H = 2 \frac{\tilde{\beta} N \int_0^1 \frac{e^{\frac{\tilde{\beta}}{N+\rho} x} e^{-\frac{\tilde{\beta}}{N(1-x)+\rho} \frac{y^2}{\eta}}}{[N(1-x)+\rho]^{3/2}} dx + \frac{e^{-\frac{\tilde{\beta}}{N+\rho} \frac{y^2}{\eta}}}{2(N+\rho)^{3/2}} + \frac{e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}}}{2\rho^{3/2}}}{\eta N \int_0^1 \frac{e^{\Lambda N x} e^{-\frac{\tilde{\beta}}{N(1-x)+\rho} \frac{y^2}{\eta}}}{[N(1-x)+\rho]^{1/2}} dx + \frac{e^{-\frac{\tilde{\beta}}{N+\rho} \frac{y^2}{\eta}}}{2(N+\rho)^{1/2}} + \frac{e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}}}{2\rho^{1/2}}} \mathcal{Y}, \quad (41) $$

To simplify this result, we can apply the change of variable $s = N(1-x) + \rho$, that allows us to prove that a universal force-extension curve exists in the thermodynamic limit ($N \to \infty$) and its shape is given by the formula

$$ \langle \mathcal{F} \rangle_H = 2 \frac{\tilde{\beta} \int_{\rho}^{+\infty} \frac{1}{s^{3/2}} e^{-\Lambda s} e^{-\frac{\tilde{\beta}}{s} \frac{y^2}{\eta}} ds + \frac{1}{2} \rho^{-3/2} e^{-\Lambda \rho} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}}}{\eta \int_{\rho}^{+\infty} \frac{1}{s^{1/2}} e^{-\Lambda s} e^{-\frac{\tilde{\beta}}{s} \frac{y^2}{\eta}} ds + \frac{1}{2} \rho^{-1/2} e^{-\Lambda \rho} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}}} \mathcal{Y} \quad (42) $$

where it is not difficult to verify that both integrals are well defined provided that

$$ \Lambda > 0, \quad (43) $$

a condition thoroughly discussed below. Now, a direct calculation of the two integrals gives

$$ \langle \mathcal{F} \rangle_H = 2 \sqrt{\frac{\tilde{\beta}}{\eta} \Lambda} \frac{\sqrt{\pi} \left\{ [1 - g^{-}(\mathcal{Y})] e^{-2\sqrt{\frac{\tilde{\beta}}{\eta}} \Lambda \mathcal{Y}} - [1 - g^{+}(\mathcal{Y})] e^{2\sqrt{\frac{\tilde{\beta}}{\eta}} \Lambda \mathcal{Y}} \right\} + \sqrt{\frac{\tilde{\beta}}{\eta}} \rho^{-3/2} e^{-\Lambda \rho} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}} \mathcal{Y}}{\sqrt{\pi} \left\{ [1 - g^{-}(\mathcal{Y})] e^{-2\sqrt{\frac{\tilde{\beta}}{\eta}} \Lambda \mathcal{Y}} + [1 - g^{+}(\mathcal{Y})] e^{2\sqrt{\frac{\tilde{\beta}}{\eta}} \Lambda \mathcal{Y}} \right\} + \sqrt{\Lambda} \rho^{-1/2} e^{-\Lambda \rho} e^{-\frac{\tilde{\beta}}{\rho} \frac{y^2}{\eta}}}, \quad (44) $$

where the functions $g_{\pm}$ are defined by means of the error function, as

$$ g^{\pm}(\mathcal{Y}) = \operatorname{erf}\left(\sqrt{\Lambda \rho \pm \sqrt{\frac{\tilde{\beta}}{\rho \eta} \mathcal{Y}}\right). \quad (45) $$

Recalling that $\langle \mathcal{F} \rangle_H = \beta y_M \langle f \rangle_H$, the force is given by

$$ \langle f \rangle_H = y_M \sqrt{kh} \sqrt{1 - \frac{T}{T_c} \frac{\sqrt{\pi} \left[ (1 - h^{-}(y_d)) e^{-\Omega(y_d)} - (1 - h^{+}(y_d)) e^{\Omega(y_d)} \right] + \frac{1}{\rho} \sqrt{\Theta_1(y_d)} e^{-\Theta_1(y_d)} e^{-\Theta_2}}{\sqrt{\pi} \left[ (1 - h^{-}(y_d)) e^{-\Omega(y_d)} + (1 - h^{+}(y_d)) e^{\Omega(y_d)} \right] + \frac{1}{\rho} \sqrt{\Theta_2} e^{-\Theta_2(y_d)} e^{-\Theta_2}}}, \quad (46) $$

where

$$ \Omega(y_d) = \frac{y_d y_M \sqrt{kh}}{k_B T} \sqrt{1 - \frac{T}{T_c}}, \quad \Theta_1(y_d) = \frac{ky_d^2}{2k_B T} \frac{1}{\rho}, \quad \Theta_2 = \frac{hy_M^2}{2k_B T} \rho, \quad (47) $$
$$h^{\pm}(y_d) = \operatorname{erf}\left(\sqrt{\frac{ky_M^2}{2k_B T}}\left[\sqrt{\eta\rho}\sqrt{1-\frac{T}{T_c}} \pm \sqrt{\frac{1}{\rho}}\frac{y_d}{y_M}\right]\right), \quad (48)$$

and, for later convenience, we introduced the critical temperature

$$T_c = \frac{hy_M^2}{k_B \ln b} = \frac{hy_M^2}{k_B \ln \frac{2+\eta+\sqrt{\eta^2+4\eta}}{2}}. \quad (49)$$

Importantly, in Eq.(44) or Eq.(46), we obtained the closed form expression for the force-extension behavior during the detachment process in the thermodynamic limit under isometric condition. It is possible to see that the main fraction in Eqs.(44) and (46) converges to 1 when $\mathcal{Y} \to +\infty$ and $y_d \to +\infty$, respectively.

Thus, we can obtain the asymptotic value

$$\langle \mathcal{F} \rangle_{as} = \lim_{\mathcal{Y} \to \infty} \langle \mathcal{F} \rangle_H = 2\sqrt{\frac{\tilde{\beta}}{\eta}\left(\tilde{\beta} - \frac{1}{2}\ln b\right)}. \quad (50)$$

Equivalently, we have

$$\langle f \rangle_{as} = \lim_{y_d \to \infty} \langle f \rangle_H = \sqrt{kh y_M} \sqrt{1 - \frac{T}{T_c}}. \quad (51)$$

This formula describes the asymptotic value of the force plateau in terms of the temperature and the material parameters of the system. Based on this equation, we may observe that the previously required condition, stated in Eq.(43) and assuring the convergence of the above integrals, corresponds to the requirement of subcritical temperatures $T < T_c$.

The obtained results are illustrated in Fig.7a, where the force displacement relation is plotted for different values of the temperature. In particular, we compare the force-displacement response of a discrete system with $N = 100$ elements given by Eq.(16) (colored dashed curves) with the result of Eq.(46) obtained in the thermodynamic limit (black continuous curves). Moreover, in Fig.7a, we also reported the asymptotic value of the decohesion force provided by Eq.(51) (horizontal straight lines). We can observe that the obtained result in the thermodynamic limit is able to represent the first force peak, which is the peculiar property of the system loaded under isometric condition. For low values of the temperature, the solution given in Eq.(46) is only an approximation

of the system response for $N \to \infty$. However, this representation can be arbitrarily improved by adding more Euler-Mclaurin remainder terms in Eq.(38) that could be substituted with

$$\sum_{\xi=0}^{N} \phi(\xi) \approx \int_{0}^{N} \phi(\xi) d\xi + \frac{\phi(0) + \phi(N)}{2} + \sum_{k=1}^{p} \frac{B_{2k}}{(2k)!} [\phi^{(2k-1)}(N) - \phi^{(2k-1)}(0)], \quad (52)$$

where $B_2 = 1/6$, $B_4 = -1/30$, $B_6 = 1/42$, $B_8 = -1/30$, and so on, are the Bernoulli numbers, and $p$ is an integer representing the number of additional remainder terms [107]. Although we tested the better approximations obtained through Eq.(52) for increasing values of $p$, the final results are much more cumbersome than Eqs.(44) and (46), and therefore in this work, for simplicity, we only use Eq.(38). In any case, the asymptotic result given in Eq.(51) is in perfect agreement with the numerical result obtained with a discrete system analyzed through Eq.(16). This behavior explains the variation of the plateau force with the temperature in terms of a phase transition occurring at the critical temperature $T_c$. This point is further discussed in the following development.

In Fig.7b, we plot the obtained temperature dependent value of the asymptotic unfolding force for different values of the non dimensional parameter $\eta$. We note that the force approaches zero at the critical temperature $T_c$, and that $T_c$ is an increasing function of the ratio $\eta$. It is important anyway to remark that our single domain wall assumption becomes questionable as we approach the critical temperature and slightly lower values of $T_c$ should be expected with an exact multi-walls solution. This aspect is out of the aim of this paper and is the subject of future work [99]. It is also worth to point out that the value of the critical force for $T = 0$ is $\langle f \rangle_{as} = \sqrt{kh y_M}$, perfectly coherent with the result obtained in the case of pure mechanical peeling [47].

To better elucidate the decohesion behavior, the meaning of the critical temperature and of the associated phase transition we further analyze the behavior of the variable $\langle \xi \rangle_H$ given in Eq.(33), measuring the average number of unbroken vertical springs. As before, by assuming a large value of $N$ and using (38), we can substitute the sums with the corresponding integrals and we find

$$\langle \xi \rangle_H = \frac{N \int_0^1 \frac{N x e^{\Lambda N x}}{[N(1-x)+\rho]^{1/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{y^2}{\eta}} \frac{1}{N(1-x)+\rho} dx + \frac{1}{2} \rho^{-1/2} N e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\eta} \frac{y^2}{\eta}}}{N \int_0^1 \frac{e^{\Lambda N x}}{[N(1-x)+\rho]^{1/2}} e^{-\frac{\tilde{\beta}}{\eta} \frac{y^2}{\eta}} \frac{1}{N(1-x)+\rho} dx + \frac{1}{2} \rho^{-1/2} e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\eta} \frac{y^2}{\eta}}}. \quad (53)$$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Panel (a): comparison between the force-extension curves given in Eq.(16) (colored dashed curves) with the thermodynamic limit obtained in Eq.(46) (black continuous curves) and the asymptotic force value given in Eq.(51) (horizontal dashed straight lines). We adopted the parameters $N = 100$, $k = 5$, $h = 20$, $y_M = 1.2$, and $\beta^{-1} = k_B T = 4.0, 5.57, 7.14, 8.71, 10.2, 11.8, 13.4$, and 15 (in arbitrary units). Panel (b): critical behavior of the asymptotic force within the Helmholtz ensemble. We plotted $\langle f \rangle_{as}$ versus $k_B T = \beta^{-1}$ (see Eq.(51)) for different values of the ratio $\eta = h/k$ between the elastic constants of the vertical and horizontal elements. Here $k = 1$, $y_M = 1$ (in arbitrary units) and we plotted ten curves with $1/10 \le \eta \le 10$.</figcaption></figure>

We can now apply the change of variable $s = N(1-x) + \rho$ and get

$$ \langle \xi \rangle_H = \frac{\int_{\rho}^{N+\rho} (N-s+\rho) s^{-1/2} e^{\Lambda(N-s+\rho)} e^{-\frac{1}{s} \frac{\tilde{\beta}}{\eta} Y^2} ds + \frac{1}{2} \rho^{-1/2} N e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\rho} \frac{Y^2}{\eta}}}{\int_{\rho}^{N+\rho} s^{-1/2} e^{\Lambda(N-s+\rho)} e^{-\frac{1}{s} \frac{\tilde{\beta}}{\eta} Y^2} ds + \frac{e^{-\frac{\tilde{\beta}}{\eta} Y^2} e^{\frac{1}{N+\rho}}}{2(N+\rho)^{1/2}} + \frac{1}{2} \rho^{-1/2} e^{\Lambda N} e^{-\frac{\tilde{\beta}}{\rho} \frac{Y^2}{\eta}}}. \quad (54) $$

Then, we can determine the expectation number $\langle \zeta \rangle_H = N - \langle \xi \rangle_H$ of broken links in the thermodynamic limit ($N \to \infty$), eventually obtaining

$$ \begin{aligned} \langle \zeta \rangle_H &= \frac{\int_{\rho}^{+\infty} s^{1/2} e^{-\Lambda s} e^{-\frac{1}{s} \alpha Y^2} ds - \rho \int_{\rho}^{+\infty} s^{-1/2} e^{-\Lambda s} e^{-\frac{1}{s} \alpha Y^2} ds}{\int_{\rho}^{+\infty} s^{-1/2} e^{-\Lambda s} e^{-\frac{1}{s} \alpha Y^2} ds + \frac{1}{2} \left(\frac{b-1}{b}\right)^{1/2} e^{-\Lambda \rho} e^{-\alpha \frac{b-1}{b} Y^2}} \\ &= \frac{1}{\Lambda} \frac{\sqrt{\pi} \left\{ e^{-2\sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y} \left( \Lambda\rho + \frac{1}{2} + \sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y \right) [1-g^{-}(Y)] + e^{2\sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y} \left( \Lambda\rho + \frac{1}{2} - \sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y \right) [1-g^{+}(Y)] \right\} + 2\sqrt{\Lambda\rho}e^{-\Lambda\rho}e^{-\frac{\tilde{\beta}}{\rho}\frac{Y^2}{\eta}}}{\sqrt{\pi} \left\{ e^{-2\sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y} [1-g^{-}(Y)] + e^{2\sqrt{\frac{\tilde{\beta}\Lambda}{\eta}}Y} [1-g^{+}(Y)] \right\} + \sqrt{\Lambda\rho}^{-1/2}e^{-\Lambda\rho}e^{-\frac{\tilde{\beta}}{\rho}\frac{Y^2}{\eta}}}. \end{aligned} \quad (55) $$

This result proves that for a long chain the number of broken elements must depend only on the extension $Y$ and on the temperature $T$.

To get a further physical interpretation of the unbroken to broken transition corresponding to the critical temperature $T_c$, we consider the limit when $Y \to \infty$, describing the full detachment state in the thermodynamic limit. From Eq.(55) we can eventually obtain, for large values of $Y$, the asymptotic relation

$$ \langle \zeta \rangle_H |_{Y \to \infty} \sim Y \left( \frac{\tilde{\beta}}{\eta \Lambda} \right)^{\frac{1}{2}}. \quad (56) $$

It means that we have a linear relation between $\langle \zeta \rangle_H$ and $Y$ for $Y \to \infty$. This expression can be written with the physical parameters of the system as

$$ \frac{\langle \zeta \rangle_H |_{y_d \to \infty}}{y_d} \sim \frac{1}{y_M} \sqrt{\frac{k}{h}} \frac{1}{\sqrt{1 - \frac{T}{T_c}}}. \quad (57) $$

This behavior is confirmed in Fig.8, where we plotted $\langle \zeta \rangle_H$ calculated through Eq.(17) or Eq.(33) for a discrete system with $N = 100$, and the quantity $\langle \zeta \rangle_H |_{y_d \to \infty}/y_d$, obtained from Eq.(57).

We observe that the attainment of the vertical asymp-
<figure><img src="image_10.png" /><figcaption>FIG. 8. Asymptotic behavior of the number of broken elements $\langle \zeta \rangle_H$ within the Helmholtz ensemble and with $y_d \to \infty$. The number of broken elements $\langle \zeta \rangle_H$ (colored continuous lines), calculated through Eq.(17) or Eq.(33), and $\langle \zeta \rangle_H|_{y_d \to \infty}$ (colored dashed straight lines), obtained from Eq.(57), are represented versus the (deformation type) quantity $y_d/\langle \zeta \rangle_H$. They show the same asymptotic behavior. We adopted the parameters $N = 100$, $k = 5$, $h = 20$, $y_M = 1.2$, and $\beta^{-1} = k_B T = 6.0, 7.8, 9.6, 11.4, 13.2$, and $15$ (in arbitrary units).</figcaption></figure>

totes indicate that as we increase $y_d$ the total displacement and the number of broken links grows to infinity with a fixed limit ratio (measuring the deformation of the detached portion) depending on the temperature. We therefore conclude that for $T \to T_c$ we have a phase transition corresponding to the rupture of all the vertical elements of the chain, i.e. to the complete detachment of the chain from the substrate. For this reason, the critical temperature can also be referred to as the denaturation temperature of the system (terminology frequently used in the biological context [34–36]). Once again we remark that the obtained value of the critical temperature can be slightly overestimated due to the neglect of solutions with more domain walls [99].

## B. Thermodynamic limit in the Gibbs ensemble

Consider now the thermodynamic limit in the case of isotensional loading. As in the previous case, we first report the expressions for the expected values of the average extension, average number of attached elements, and displacement vector, given in Eqs.(27), (28), and (31). After nondimensionalization and the use of Eqs.(96)-(98), we obtain

$$ \langle Y \rangle_G = \frac{\langle y_d \rangle}{y_M} = \frac{\eta \mathcal{F}}{2 \tilde{\beta} \mathcal{E}} \sum_{\xi=0}^{N} e^{\tilde{\beta} \xi} \frac{(N - \xi + 1)\gamma(\xi + 1) - (N - \xi)\gamma(\xi)}{[\gamma(\xi + 1) - \gamma(\xi)]^{3/2}} e^{\frac{\eta}{4\tilde{\beta}} \frac{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}{\gamma(\xi+1)-\gamma(\xi)} \mathcal{F}^2}, \quad (58) $$

$$ \langle \xi \rangle_G = \frac{1}{\mathcal{E}} \sum_{\xi=0}^{N} \xi e^{\tilde{\beta} \xi} \frac{1}{[\gamma(\xi + 1) - \gamma(\xi)]^{1/2}} e^{\frac{\eta}{4\tilde{\beta}} \frac{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}{\gamma(\xi+1)-\gamma(\xi)} \mathcal{F}^2}, \quad (59) $$

$$ \langle Y_i \rangle_G = \frac{\langle y_i \rangle_G}{y_M} = \frac{\eta \mathcal{F}}{2 \tilde{\beta} \mathcal{E}} \sum_{\xi=0}^{i-1} e^{\tilde{\beta} \xi} \frac{(i - \xi)\gamma(\xi + 1) - (i - \xi - 1)\gamma(\xi)}{[\gamma(\xi + 1) - \gamma(\xi)]^{3/2}} e^{\frac{\eta}{4\tilde{\beta}} \frac{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}{\gamma(\xi+1)-\gamma(\xi)} \mathcal{F}^2} \\ + \frac{\eta \mathcal{F}}{2 \tilde{\beta} \mathcal{E}} \sum_{\xi=i}^{N} e^{\tilde{\beta} \xi} \frac{\gamma(i)}{[\gamma(\xi + 1) - \gamma(\xi)]^{3/2}} e^{\frac{\eta}{4\tilde{\beta}} \frac{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}{\gamma(\xi+1)-\gamma(\xi)} \mathcal{F}^2}, \quad i = 1, \dots, N \quad (60) $$

$$ \mathcal{E} = \sum_{\xi=0}^{N} e^{\tilde{\beta} \xi} \frac{1}{[\gamma(\xi + 1) - \gamma(\xi)]^{1/2}} e^{\frac{\eta}{4\tilde{\beta}} \frac{(N-\xi+1)\gamma(\xi+1)-(N-\xi)\gamma(\xi)}{\gamma(\xi+1)-\gamma(\xi)} \mathcal{F}^2} \quad (61) $$

where $\tilde{\beta}$ is defined in Eq.(36), the adimensional force is given by $\mathcal{F} = \beta f y_M$, and the function $\gamma(z)$ is defined in Eq.(37).

Consider the force-extension relation in Eq.(58). In

the limit of large $N$, we substitute the summations with the corresponding integrals, as described in Eq.(38). As before, we can approximate the function $\gamma(z)$ defined in Eq.(37) as $\gamma(z) \approx \frac{1}{\sqrt{\eta^2+4\eta}} b^z$. First, we get

$$ \langle Y \rangle_G = \frac{N \int_0^1 e^{\Lambda N x} [N(1-x) + \rho] e^{\frac{\eta}{4\beta} \mathcal{F}^2 [N(1-x)+\rho]} dx + \frac{1}{2} (N+\rho) e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} \rho e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}{N \int_0^1 e^{\Lambda N x} e^{\frac{\eta}{4\beta} \mathcal{F}^2 [N(1-x)+\rho]} dx + \frac{1}{2} e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}} \frac{\eta \mathcal{F}}{2 \tilde{\beta}}. \quad (62) $$
By considering the change of variable $s = N(1 - x) + \rho$, we obtain

$$ \langle \mathcal{Y} \rangle_G = \frac{\int_{\rho}^{N+\rho} s e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} (N+\rho) e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} \rho e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}{\int_{\rho}^{N+\rho} e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}} \frac{\eta \mathcal{F}}{2\tilde{\beta}}. \quad (63) $$

In the thermodynamic limit ($N \to \infty$), we have

$$ \langle \mathcal{Y} \rangle_G = \frac{\int_{\rho}^{+\infty} s e^{-(\tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\beta} \mathcal{F}^2)s} ds + \frac{1}{2} \rho e^{-(\tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\beta} \mathcal{F}^2)\rho}}{\int_{\rho}^{+\infty} e^{-(\tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\beta} \mathcal{F}^2)s} ds + \frac{1}{2} e^{-(\tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\beta} \mathcal{F}^2)\rho}} \frac{\eta \mathcal{F}}{2\tilde{\beta}}. \quad (64) $$

Both integrals in Eq.(64) converge if

$$ \tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\tilde{\beta}} \mathcal{F}^2 > 0. \quad (65) $$

This condition will be thoroughly discussed in the following. By integrating, we eventually obtain the force displacement relation in the thermodynamic limit

$$ \langle \mathcal{Y} \rangle_G = \left( \rho + \frac{1}{\tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\tilde{\beta}} \mathcal{F}^2} \times \frac{2}{2 + \tilde{\beta} - \frac{1}{2} \ln b - \frac{\eta}{4\tilde{\beta}} \mathcal{F}^2} \right) \frac{\eta \mathcal{F}}{2\tilde{\beta}}. \quad (66) $$

It is easy to verify that $\langle \mathcal{Y} \rangle_G$ is an increasing function of $\mathcal{F}$ and it diverges when the forces attains the same asymptotic value given in Eq.(50), obtained in the case of assigned displacement (Helmholtz ensemble). The force-extension relation given in Eq.(66) can be written in terms of the temperature and the other physical parameters of the system, as follows

$$ \langle y_d \rangle_G = \left[ \frac{1}{1 - e^{-\frac{hy_M^2}{k_B T_c}}} + \frac{1}{\frac{hy_M^2 \beta}{2} \left(1 - \frac{T}{T_c} - \frac{f^2}{khy_M^2}\right)} \times \frac{2}{2 + \frac{hy_M^2 \beta}{2} \left(1 - \frac{T}{T_c} - \frac{f^2}{khy_M^2}\right)} \right] \frac{f}{k}. \quad (67) $$

As before, $\langle y_d \rangle_G$ is an increasing function of $f$ and it diverges when the force attains the same asymptotic value given in Eq.(51), obtained for the Helmholtz ensemble with

$$ \mathcal{F} \le \mathcal{F}_{as} = 2 \sqrt{\frac{\tilde{\beta}}{\eta} \left(\tilde{\beta} - \frac{1}{2} \ln b\right)}, \quad (68) $$

thus interpreting the condition in Eq.(65) or, equivalently

$$ f \le f_{as} = \sqrt{kh y_M} \sqrt{1 - \frac{T}{T_c}}. \quad (69) $$

The other important parameter of the system is the number of unbroken vertical elements, which can be calculated through Eq.(59). This expression can be elaborated in the limit $N \to \infty$, obtaining the asymptotic form

$$ \langle \xi \rangle_G = \frac{N \int_0^1 N x e^{\Lambda N x} e^{\frac{\eta}{4\beta} \mathcal{F}^2 [N(1-x)+\rho]} dx + \frac{1}{2} N e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}{N \int_0^1 e^{\Lambda N x} e^{\frac{\eta}{4\beta} \mathcal{F}^2 [N(1-x)+\rho]} dx + \frac{1}{2} e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}. \quad (70) $$

As before, we can use the change of variable $s = N(1 - x) + \rho$, which delivers

$$ \langle \xi \rangle_G = \frac{\int_{\rho}^{N+\rho} (N+\rho-s) e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} N e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}{\int_{\rho}^{N+\rho} e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}. \quad (71) $$

From now on, we consider the average number of broken vertical elements $\langle \zeta \rangle_G = N - \langle \xi \rangle_G$ given by

$$ \langle \zeta \rangle_G = \frac{\int_{\rho}^{N+\rho} (s-\rho) e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} N e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)}}{\int_{\rho}^{N+\rho} e^{\Lambda(N+\rho-s)} e^{\frac{\eta}{4\beta} \mathcal{F}^2 s} ds + \frac{1}{2} e^{\frac{\eta}{4\beta} \mathcal{F}^2 (N+\rho)} + \frac{1}{2} e^{\Lambda N} e^{\frac{\eta}{4\beta} \mathcal{F}^2 \rho}}. \quad (72) $$
<figure><img src="image_11.png" /><figcaption>FIG. 9. Behavior of the detachment process within the Gibbs ensemble in terms of force-extension relation, panel (a), and average number of broken elements, panel (b). We compare the results in the thermodynamic limit given in Eqs.(67) and (73) (colored dashed curves) with their discrete counterparts given in Eqs.(58) and (59) (black curves), applied to a system with N = 100. The horizontal, in panel (a), and vertical, in panel (b), dashed straight lines represent the values of the force that induces complete decohesion of the system. We adopted the parameters k = 5, h = 20, y<sub>M</sub> = 1.2, and β<sup>-1</sup> = k<sub>B</sub>T = 4.0, 5.57, 7.14, 8.71, 10.2, 11.8, 13.4, and 15 (in arbitrary units).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. Comparison between the Helmholtz force-extension curves (⟨f⟩<sub>H</sub> versus y<sub>d</sub>, solid lines) and the Gibbs force-extension curves (f versus ⟨y<sub>d</sub>⟩<sub>G</sub>, dashed lines) proving the non-equivalence of the two statistical ensembles in the thermodynamic limit. Panel (a): comparison between the results for a discrete system with N = 100 based on Eqs.(16) and (27). Panel (b): comparison between the results in the thermodynamic limit based on Eqs.(46) and (67). In both panels, the horizontal dashed straight lines correspond to ⟨f⟩<sub>as</sub> given in Eq.(51). We adopted the parameters k = 5, h = 20, y<sub>M</sub> = 1.2, and β<sup>-1</sup> = k<sub>B</sub>T = 4.0, 5.57, 7.14, 8.71, 10.2, 11.8, 13.4, and 15 (in arbitrary units).</figcaption></figure>

This result proves that the limit of ζ for N → ∞ exists and can be written as follows

$$
\begin{align}
\langle \zeta \rangle_G &= \frac{\int_{\rho}^{+\infty} (s-\rho) e^{-(\tilde{\beta}-\frac{1}{2}\ln b-\frac{\eta}{4\beta}\mathcal{F}^2)} ds}{\int_{\rho}^{+\infty} e^{-(\tilde{\beta}-\frac{1}{2}\ln b-\frac{\eta}{4\beta}\mathcal{F}^2)} ds + \frac{1}{2} e^{-(\tilde{\beta}-\frac{1}{2}\ln b-\frac{\eta}{4\beta}\mathcal{F}^2)} \rho} \nonumber \\
&= \frac{1}{\tilde{\beta}-\frac{1}{2}\ln b-\frac{\eta}{4\beta}\mathcal{F}^2} \times \frac{2}{2+\tilde{\beta}-\frac{1}{2}\ln b-\frac{\eta}{4\beta}\mathcal{F}^2} \nonumber \\
&= \frac{1}{\frac{hy_M^2\beta}{2}\left(1-\frac{T}{T_c}-\frac{f^2}{khy_M^2}\right)} \times \frac{2}{2+\frac{hy_M^2\beta}{2}\left(1-\frac{T}{T_c}-\frac{f^2}{khy_M^2}\right)}, \tag{73}
\end{align}
$$
where the critical temperature $T_c$ is given in Eq.(49). Coherently with previous results we observe that the number of detached element monotonically grows with $\mathcal{F}$ and diverges for the same force threshold in Eq.(50) or Eq.(51), corresponding to a transition to the fully detached configuration. This point confirms the presence of the phase transition for the isotensional condition as well.

In Fig.9 we compare the results obtained in the thermodynamic limit, Eqs.(67) and (73), with their discrete counterparts given in Eqs.(58) and (59), applied to a system with $N = 100$. Again, this figure shows the practical utility of the fully explicit expressions in the thermodynamic limit since the two behaviors are significantly superimposed also in the case of isotensional loading.

### C. Comparison of the Helmholtz and Gibbs ensembles

In this section we focus on the important differences of the decohesion behavior that we obtained under the different loading conditions. The comparison between the hard and soft device loading can be drawn by observing the force-extension relation in the corresponding Helmholtz and Gibbs ensembles reported in Fig.10. In particular, in Fig.10a, we plotted Eqs.(16) and (27), representing the force-extension response for a discrete system, where we used a large value of $N$ and different values of the temperature $T$. On the other hand, in Fig.10b, we plotted Eqs.(46) and (67), representing the force-extension behavior in the thermodynamic limit. Observe the force-extension behavior is markedly different in the two ensembles. In the case of hard device the system starts to unfold at a threshold force larger than $\langle f \rangle_{as}$ and (at low temperatures) oscillates around this value approaching the asymptotic limit only after the initial discrete debonding process. Differently, in the case of isotensional loading the force monotonically increases attaining the limit value given in Eq.(51) only when the total displacement diverges and the whole detachment is observed. This behavior can be seen in both panels of Fig.10, proving that the analytic expressions found in the case of thermodynamic limit are well adapted to represent the system behavior for large values of $N$ also in the case when the only first correction term of the Euler-McLaurin formula is considered.

The difference between the Helmholtz and Gibbs ensemble can be clearly appreciated also in Fig.11. In panel (a), we show a zoom on the first peak of the force extension relation within the Helmholtz ensemble for several temperatures in the range between zero and the critical temperature. These curves are compared with the asymptotic value $\langle f \rangle_{as}$ given in Eq.(51). In addition, the values of these Helmholtz force peaks, given by $\max\{\langle f \rangle_H\}$, have been plotted in panel (b) versus the system temperature (red curve). For the Gibbs case, the largest force corresponds to the asymptotic one and we

can therefore write: $\sup\{f\} = \langle f \rangle_{as}$. This quantity is represented by the blue curve in panel (b) of Fig.11. To conclude, we underline that the decohesion force threshold for the Helmholtz case is sensibly larger than the same quantity for the Gibbs ensemble, especially for low values of the temperature.

The non-equivalence of the ensembles in the decohesion process can be further appreciated by comparing the initial slope of the force extension curves, represented by Eqs.(46) and (66) in the thermodynamical limit. Indeed, we can calculate the effective stiffness observed by the device at the beginning of the peeling process, i.e., for small values of force and extension. It means that we have to determine the derivatives of the curves given in Eqs.(46) and (66) for small force or extension. Their calculation give

$$ \left. \frac{\partial \langle \mathcal{F} \rangle_H}{\partial \mathcal{Y}} \right|_{\mathcal{Y}=0} = 4 \frac{\tilde{\beta}}{\eta} \Lambda \frac{2 \left( 1 + \frac{1}{4\rho} \right) e^{-\Lambda\rho} - \Theta}{\Lambda e^{-\Lambda\rho} + \Theta}, \quad (74) $$

$$ \left( \left. \frac{\partial \langle \mathcal{Y} \rangle_G}{\partial \mathcal{F}} \right|_{\mathcal{F}=0} \right)^{-1} = 2 \frac{\tilde{\beta}}{\eta} \Lambda \frac{1}{\Lambda\rho + \frac{2}{2+\Lambda}}, \quad (75) $$

for the Helmholtz ensemble and the Gibbs ensemble, respectively. Here, for the sake of simplicity, we defined the quantity

$$ \Theta = 2\sqrt{\pi\Lambda\rho[1 - \erf(\Lambda\rho)]}. \quad (76) $$

In terms of physical parameters the effective stiffnesses are

$$ k_{eff}^{H} = \left. \frac{\partial \langle f \rangle_H}{\partial y_d} \right|_{y_d=0} = \frac{k_B T}{y_M^2} \left. \frac{\partial \langle \mathcal{F} \rangle_H}{\partial \mathcal{Y}} \right|_{\mathcal{Y}=0}, \quad (77) $$

$$ k_{eff}^{G} = \left( \left. \frac{\partial \langle y_d \rangle_G}{\partial f} \right|_{f=0} \right)^{-1} = \frac{k_B T}{y_M^2} \left( \left. \frac{\partial \langle \mathcal{Y} \rangle_G}{\partial \mathcal{F}} \right|_{\mathcal{F}=0} \right)^{-1} . \quad (78) $$

The plot of the effective stiffness defined in Eqs.(77) and (78) versus the temperature can be found in Fig.12a, from which we can deduce again the non-equivalence of the ensembles. The behavior of the stiffness is similar for the two ensembles for low values of the temperature. In this case, we have the common limiting value

$$ k_{eff}^{H}(T=0) = k_{eff}^{G}(T=0) = \frac{k}{\rho}, \quad (79) $$

with $\rho$ defined in Eq.(40), in agreement with the results obtained in Ref.[47]. On the contrary, the ratio between the two stiffnesses diverges as the temperature converges to its critical value (see Fig.12b). The behavior of the two stiffness coefficients plotted in Fig.12 is responsible for the markedly different initial slope of the force extension curves shown in Fig.10.

## VI. DISCUSSION AND CONCLUSION

In this paper, we elaborated a model to describe the decohesion/decohesion process related to a film deposited
<figure><img src="image_13.png" /><figcaption>FIG. 11. Comparison between the Helmholtz and the Gibbs response. Panel (a): zoom of the first peak in the force extension curves for the Helmholtz ensemble, based on Eq.(16) (colored curves), and asymptotic value $\langle f \rangle_{as}$ given in Eq.(51) (black horizontal lines). Panel (b): critical behavior of both ensembles, represented by the quantity $\max\{\langle f \rangle_H\}$ for the isometric condition (red curve) and by $\sup\{f\} = \langle f \rangle_{as}$ for the isotensional one (blue curve). We adopted the parameters $N = 100$, $k = 5$, $h = 20$, $y_M = 1.2$, and $\beta^{-1} = k_B T$ assumes 48 values from zero to $k_B T_c$ (in arbitrary units).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 12. Effective stiffness within the Helmholtz and Gibbs ensembles, in panel (a), and their ratio in panel (b). We adopted the parameters $k = 5$, $h = 20$, $y_M = 1.2$ (in arbitrary units).</figcaption></figure>

on a substrate, by focusing on the effects of thermal fluctuations. The paradigmatic system adopted is composed of a one-dimensional elastic chain grounded to a substrate through a series of breakable links. We analyzed this process under an additional mechanical action, represented by either an external force or a prescribed extension applied to the end-terminal of the chain. These two conditions correspond to the Gibbs and the Helmholtz ensembles in the framework of equilibrium statistical mechanics, respectively. Based on a spin variables approach [71, 72] and using some known properties of tridiagonal matrices [108, 109], we are able to obtain an analytical description of the decohesion behavior of the system with a resulting clear physical interpretation of the results. We firstly developed the theory for the Helmholtz ensemble

and then we obtained the results for the Gibbs one with the Laplace transform describing the relationship between the partition functions of the two ensembles [85]. Eventually, for both statistical ensembles, we obtained explicit force-extension relations, the average number of broken units, and the average extension of all the elements of the chain as function of the temperature and of the external mechanical action (a force for the Gibbs ensemble and an extension for the Helmholtz one). These achievements, summed up in Eqs.(32)-(35) (Helmholtz) and Eqs.(58)-(61) (Gibbs), and obtained for an arbitrary number $N$ of elements of the chain, are useful to fully understand the behavior of the system, to compare with existing results concerning the folding/unfolding of macromolecular bistable chains, and finally to perform the anal-
ysis of the thermodynamic limit. Indeed, when $N \to \infty$, the sums in Eqs.(32)-(35) (Helmholtz) and Eqs.(58)-(61) (Gibbs) can be substituted by suitable integrals. As we showed, to obtain a more detailed description of the results and in particular to describe the existence of a force peak anticipating the decohesion force in the hard device we adopted the Euler-McLaurin approximation formula. This approach will be in our opinion useful in many other limit analysis in Statistical Mechanics and indeed we have shown its importance in describing the fundamental differences of the behavior in the thermodynamic limit for the hard and soft device. In particular, such an analysis allows us to prove that in the thermodynamic limit the decohesion of the film from the substrate takes place at a given critical force, which is temperature dependent. This is an important difference between the cohesion/decohesion process and the folding/unfolding of bistable chains, where the transition force is temperature independent. More explicitly, in the case of bistable chains only the slope of the transition path changes with the temperature, while keeping fixed the average value of the transition plateau. In the cohesion/decohesion process the origin of the temperature dependent peeling force is explained through the observation of a phase transition taking place at a given critical temperature able to fully detach the film from the substrate. In the subcritical regime, the thermal fluctuations promote the detachment of the film and, therefore, a lower peeling force is needed for higher temperatures. The decreasing trend of the peeling force with the temperature is the same for both statistical ensembles. However, these ensembles are nonequivalent in the thermodynamic limit since they show a different force-extension curve. In particular, the force extension curve for the Helmholtz case is characterized by a force peaks followed by some oscillations before reaching the asymptotic force value. On the contrary, the force extension curve for the Gibbs case is always monotonically increasing from zero to the asymptotic force value. This observation leads to two different critical behaviors, as shown in Fig.11b.

The existence of finite-temperature phase transitions in low-dimensional many-body models is a subject of large interest in theoretical physics [110, 111]. It is useful to remark that the existence of a genuine phase transition in our system (at thermodynamic limit and for both isometric and isotensional conditions) is coherent with the observation (both theoretical and numerical) of different kinds of phase transitions in the Peyrard model for the DNA thermal denaturation [112]. In both cases, the system is one-dimensional along the longitudinal direction with a one-dimensional series of interactions in the transverse direction. Another interesting example of phase transition in low-dimensional systems concerns the tension-induced binding of two parallel semi-flexible polymers [113]. In this case, the mean-field theory predicts a phase transition describing the discontinuous increasing bonding of the chains with an increasing applied force. However, the authors explain that the transition

turns into a crossover if the mean-field theory is substituted with the exact solution. On the other hand, the observation of a genuine phase transition in our model is due to important differences as compared with Ref.[113]: the force is transversal in our case whereas the binding tension is longitudinal in Ref.[113]; moreover, we adopted a flexible chain whereas a two semi-flexible chains are introduced in Ref.[113]; finally, our model is discrete, while the model of Ref.[113] has continuous worm like chains with discrete links. As we show in this paper, these differences can lead to different critical behaviors and different universal classes.

Although the equivalence of the Gibbs and Helmholtz statistical ensemble has been proved for a large class of systems (namely, single flexible polymer chains without confinement effects and with a continuous pairing interaction potential between neighboring monomers [86–91]), interestingly, it cannot be assumed for other structures. Indeed, other classes of problems, e.g. concerning the escape of a polymer confined between two surfaces and the desorption of a polymer initially tethered onto a surface, exhibit an unusual non-equivalence between the defined statistical ensembles [114–119]. In particular, in Refs.[114, 115] an end-tethered polymer chain compressed between two pistons is considered and shows non-equivalence of the ensembles and a phase transition corresponding to the escape from the gap between the pistons. Ref.[116] deals with the desorption of a single chain from a substrate without excluded volume interactions: also in this case the equivalence of the ensembles and the emergence of a phase transition are thoroughly discussed. In Ref.[117] a Gaussian chain is tethered on a rigid planar surface at one end and the ensemble non-equivalence comes from a pure confinement effect and does not involve any potential interaction, unlike in our case. Different rectangular, spherical and cylindrical geometries of confinement have been considered in Ref.[118]. Finally, in Ref.[119], the desorption of a self-avoiding polymer chain from a surface has been studied yielding the non-equivalence of the ensembles, the existence of a first-order phase transition without phase coexistence and a quantitative relation between adsorption exponent and adsorption energy. These results have been theoretically proved by means of the grand canonical ensemble method and confirmed by Monte Carlo simulations. Such investigations are coherent with our achievements and prove the possibility to have the non-equivalence between different (canonical) ensembles in statistical mechanics. In systems with strong interactions and low dimensionality, we always observe strong fluctuations, which are persisting also in the thermodynamic limit, thus inducing a different behavior for different statistical ensembles [117, 118, 120]. In our system the strong interactions, which are non-local and long-range, are generated by the ladder network geometry composed of linear and breakable springs. As discussed in Refs.[117, 118, 120], the resulting fluctuations of the macroscopic observables cause a non-uniform convergence to the thermodynamic limit eventually producing
the ensembles non-equivalence.

Another important result concerns the temperature dependent effective stiffness of the system, as reported in Fig.12. In particular, the model describes the decreasing of the stiffness with the temperature. Moreover, also in this case we observe a loading type dependent stiffness with the ratio between Helmholtz and Gibbs stiffness diverging as the temperature approaches its critical value. We remark that the considered dual types of loading processes represent limiting values of real applied conditions as the stiffness of the loading device changes [74, 75]. Thus, we observe that the real response in term of stiffness and decohesion force cannot be considered as an intrinsic property of the system, since it depends on the loading condition. For example, we can think to the case of the variable stiffness of the atomic force microscopy device or to the interacting molecule (RNA vs DNA) and so on. These real cases are placed in-between the ideal Helmholtz and Gibbs ensembles.

It is important to remark that while the presented model is interesting because it leads to a fully analytic approach able to explain the emergence of the phase transition and the non-equivalence of the ensembles, it could be generalized to take into consideration more complex effects. For instance, the spin variable methodology, here employed to calculate the partition functions, is limited to the study of the equilibrium thermodynamics. It could be interesting to generalize it to the dynamic regime, where the whole energy landscape and in particular the energy barriers may play a crucial role [121–123]. Moreover, our approach could be applied to different system configurations including fiber bundles with breakable strands characterized by variable detaching thresholds, buckling of films deposited over substrates, cracks propagation in brittle or plastic solids and rupture phenomena in polymer networks.

## APPENDIX: PROPERTIES OF $\mathcal{A}_N^\xi(\eta)$

We prove here some properties concerning the matrix $\mathcal{A}_N^\xi(\eta)$, defined in Eq.(3). To begin, we consider the following arbitrary tridiagonal matrix $\mathcal{T}$

$$ \mathcal{T} = \begin{bmatrix} a_1 & b_1 & 0 & \cdots & 0 \\ c_1 & a_2 & b_2 & \cdots & \vdots \\ 0 & c_2 & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & a_{N-1} & b_{N-1} \\ 0 & \cdots & 0 & c_{N-1} & a_N \end{bmatrix} \in \mathcal{M}_{N,N}(\mathbb{R}), \quad (80) $$

where the diagonal is composed by the elements $(a_1, ..., a_N)$, the superdiagonal by $(b_1, ..., b_{N-1})$ and the subdiagonal by $(c_1, ..., c_{N-1})$. It has been proved [108, 109] that the elements of the inverse matrix $\mathcal{T}^{-1}$ can be

represented as

$$ [\mathcal{T}^{-1}]_{ij} = \begin{cases} \frac{1}{\partial_N} (-1)^{i+j} b_i \times \dots \times b_{j-1} \partial_{i-1} \varphi_{j+1}, & i < j \\ \frac{1}{\partial_N} \partial_{i-1} \varphi_{i+1}, & i = j \\ \frac{1}{\partial_N} (-1)^{i+j} c_j \times \dots \times c_{i-1} \partial_{j-1} \varphi_{i+1}, & i > j \end{cases} \quad (81) $$

where the sequences $\vartheta_i$ and $\varphi_i$ are given by the recursive relations

$$ \left\{ \begin{array}{l} \vartheta_i = a_i \vartheta_{i-1} - b_{i-1} c_{i-1} \vartheta_{i-2}, \quad i = 2, \dots, N, \\ \vartheta_0 = 1, \vartheta_1 = a_1, \end{array} \right. \quad (82) $$

and

$$ \left\{ \begin{array}{l} \varphi_i = a_i \varphi_{i+1} - b_i c_i \varphi_{i+2}, \quad i = N, \dots, 1, \\ \varphi_{N+2} = 0, \varphi_{N+1} = 1, \varphi_N = a_N. \end{array} \right. \quad (83) $$

While Eq.(82) is an increasing recursive relation going from $i=1$ to $i=N$, Eq.(83) is a decreasing recursive relation going from $i=N$ to $i=1$. We also remember that $\det \mathcal{T} = \vartheta_N$ [108, 109].

In the case of the matrix $\mathcal{A}_N^\xi(\eta)$ in Eq.(3), we have that $b_i = c_i = -1 \forall i$. Under this hypothesis, the general result can be simplified as follows

$$ [\mathcal{T}^{-1}]_{ij} = \begin{cases} \frac{1}{\partial_N} \vartheta_{i-1} \varphi_{j+1}, & i < j \\ \frac{1}{\partial_N} \vartheta_{i-1} \varphi_{i+1}, & i = j \\ \frac{1}{\partial_N} \vartheta_{j-1} \varphi_{i+1}, & i > j \end{cases} \quad (84) $$

where the sequences $\vartheta_i$ and $\varphi_i$ are given by the simplified recursive relations

$$ \left\{ \begin{array}{l} \vartheta_i = a_i \vartheta_{i-1} - \vartheta_{i-2}, \quad i = 2, \dots, N, \\ \vartheta_0 = 1, \vartheta_1 = a_1, \end{array} \right. \quad (85) $$

and

$$ \left\{ \begin{array}{l} \varphi_i = a_i \varphi_{i+1} - \varphi_{i+2}, \quad i = N, \dots, 1, \\ \varphi_{N+2} = 0, \varphi_{N+1} = 1, \varphi_N = a_N. \end{array} \right. \quad (86) $$

These properties of the tridiagonal matrices are used here to determine the determinant and the inverse of our main tridiagonal matrix. Such properties have been exploited to perform analytically these algebraic operations, which can be rather expensive from the numerical point of view. It is important to remark that the tridiagonal structure of our matrices comes from the geometry of the lattice structure used in our investigation (with horizontal and vertical springs). In addition, the physical hypothesis of considering a single domain wall induces a block structure of the tridiagonal matrix, as described in Eq.(3) and following equations. More specifically, the introduction of the single domain wall assumption reduces previous properties to difference equations with piecewise constant coefficients, which can be easily solved as discussed in this Appendix.
Since we need to determine det $\mathcal{A}_N^\xi(\eta)$, $[\mathcal{A}_N^\xi(\eta)]^{-1} \vec{v}$.
$\vec{v} = [\mathcal{A}_N^\xi(\eta)]_{NN}^{-1}$, and $[\mathcal{A}_N^\xi(\eta)]^{-1} \vec{v}$, it is sufficient to calculate only the sequence $\vartheta_i$ associated to the matrix $\mathcal{A}_N^\xi(\eta)$. Indeed, we have that

$$ \det \mathcal{A}_N^\xi(\eta) = \vartheta_N, \quad (87) $$

$$ \begin{aligned} \vec{v}^T [\mathcal{A}_N^\xi(\eta)]^{-1} \vec{v} &= [\mathcal{A}_N^\xi(\eta)]_{NN}^{-1} \\ &= \vartheta_{N-1} \varphi_{N+1} / \vartheta_N = \vartheta_{N-1} / \vartheta_N, \quad (88) \end{aligned} $$

$$ \left\{ \left[ \mathcal{A}_N^\xi(\eta) \right]^{-1} \vec{v} \right\}_i = \vartheta_{i-1} \varphi_{N+1} / \vartheta_N = \vartheta_{i-1} / \vartheta_N. \quad (89) $$

Therefore, we consider Eq.(85) with $a_1 = \dots = a_\xi = 2 + \eta$ and $a_{\xi+1} = \dots = a_N = 2$. Then, for $i \le \xi$ we have the difference equation $\vartheta_i = (2 + \eta)\vartheta_{i-1} - \vartheta_{i-2}$, whose general solution can be written as

$$ \vartheta_i = p \left( \frac{2 + \eta + \sqrt{\Delta}}{2} \right)^i + q \left( \frac{2 + \eta - \sqrt{\Delta}}{2} \right)^i, \quad (90) $$

with $\Delta = \eta^2 + 4\eta$ and where the coefficients $p$ and $q$ must be fixed through the condition $\vartheta_0 = 1$ and $\vartheta_1 = a_1$. A straightforward calculation leads to the explicit solution for $i \le \xi$

$$ \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}. \quad (91) $$

For $i > \xi$, we have the simpler difference equation $\vartheta_i = 2\vartheta_{i-1} - \vartheta_{i-2}$, with the general solution $\vartheta_i = r + si$. In this case, the coefficients $r$ and $s$ must be obtained by imposing $\vartheta_{\xi-1}$ and $\vartheta_\xi$ by means of Eq.(90). Hence, the result for $i > \xi$ can be eventually found as

$$ \vartheta_i = (i - \xi + 1)\vartheta_\xi - (i - \xi)\vartheta_{\xi-1}, \quad (92) $$

where $\vartheta_{\xi-1}$ and $\vartheta_\xi$ are given by Eq.(90). The obtained results can be summarized through the final expression holding for $1 \le i \le N$

$$ \vartheta_i = \begin{cases} \gamma(i+1), & i \le \xi \\ (i - \xi + 1)\gamma(\xi + 1) - (i - \xi)\gamma(\xi), & i > \xi \end{cases} \quad (93) $$

where the function $\gamma(z)$ is defined as follows

$$ \gamma(z) = \frac{1}{\sqrt{\Delta}} \left( \frac{2 + \eta + \sqrt{\Delta}}{2} \right)^z - \frac{1}{\sqrt{\Delta}} \left( \frac{2 + \eta - \sqrt{\Delta}}{2} \right)^z. \quad (94) $$

A more compact form of the solution can be also written as

$$ \vartheta_i = \left[ \gamma(i+1) - (i-\xi+1)\gamma(\xi+1) + (i-\xi)\gamma(\xi) \right] \mathbf{1}(\xi-i) + \left[ (i-\xi+1)\gamma(\xi+1) - (i-\xi)\gamma(\xi) \right], \quad (95) $$

in terms of the Heaviside step function $\mathbf{1}(x)$, defined as $\mathbf{1}(x) = 1$ if $x \ge 0$, and $\mathbf{1}(x) = 0$ if $x < 0$. To conclude, Eqs.(87)-(89) assume the final form

$$ \det \mathcal{A}_N^\xi(\eta) = (N - \xi + 1)\gamma(\xi + 1) - (N - \xi)\gamma(\xi), \quad (96) \\ 1 - \left[ \mathcal{A}_N^\xi(\eta) \right]_{NN}^{-1} = \frac{\gamma(\xi + 1) - \gamma(\xi)}{(N - \xi + 1)\gamma(\xi + 1) - (N - \xi)\gamma(\xi)}, \quad (97) $$

$$ \left\{ \left[ \mathcal{A}_N^\xi(\eta) \right]^{-1} \vec{v} \right\}_i = \frac{\vartheta_{i-1}}{(N - \xi + 1)\gamma(\xi + 1) - (N - \xi)\gamma(\xi)}, \quad (98) $$

where $\vartheta_i$ is given in Eq.(95) and the function $\gamma(z)$ in Eq.(94).

## ACKNOWLEDGMENTS

SG acknowledges the region “Hauts de France” for the financial support under project MEPOFIB. GF and GP have been supported by the Italian Ministry MIUR-PRIN project “Mathematics of active materials: From mechanobiology to smart devices”, 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” and by the FFABR research grant (MIUR).

[1] S. M. Smith, Y. Cui, and C. Bustamante, Science 271, 795 (1996).

[2] M. Rief, F. Oesterhelt, B. Heymann, and H. E. Gaub, Science 275, 28 (1997).

[3] M. Rief, M. Gautel, F. Oesterhelt, J. M. Fernandez, and H. E. Gaub, Science 276, 1109 (1997).

[4] M. Rief, J. M. Fernandez, H. E. Gaub, Phys. Rev. Lett. 81, 4764 (1998).

[5] C. Storm and P. C. Nelson, Phys. Rev. E 67, 051906 (2003).

[6] F. Ritort, J. Phys.: Condens. Matter 18, R531 (2006).

[7] K. C. Neuman and A. Nagy, Nat. Meth. 5, 491 (2008).

[8] S. Kumar and M. S. Li, Phys. Rep. 486, 1 (2010).

[9] T. Hoffmann and L. Dougan, Chem. Soc. Rev. 41, 4781 (2012).

[10] O. K. Dudko, Q. Rev. Biophys. 49, e3 (2016).

[11] M. L. Hughes and L. Dougan, Rep. Prog. Phys. 79, 076601 (2016).

[12] A. Prados, A. Carpio, and L. L. Bonilla, Phys. Rev. E 88, 012704 (2013).

[13] L. L. Bonilla, A. Carpio, and A. Prados, Phys. Rev. E 91, 052712 (2015).

[14] D. De Tommasi, N. Millardi, G. Puglisi, and G. Saccomandi, J. R. Soc. Interface 10, 20130651 (2013).

[15] G. I. Bell, Science 200, 618 (1978).

[16] G. I. Bell, M. Dembo, and P. Bongrand, Biophys J. 45, 1051 (1984).

[17] B. M. Gumbiner, Cell. 84, 345 (1996).
[18] T. Erdmann and U. S. Schwarz, Eur. Phys. J. E **22**, 123 (2007).

[19] H. Gao, J. Qian, and B. Chen, J. R. Soc. Interface **8**, 1217 (2011).

[20] U. S. Schwarz and S. A. Safran, Rev. Mod. Phys. **85**, 1327 (2013).

[21] J. Liphardt, B. Onoa, S.B. Smith, I. Tinoco, and C. Bustamante, Science **292**, 733 (2001).

[22] O. K. Dudko, J. Mathé, A. Szabo, A. Meller, and G. Hummer, Biophys. J. **92**, 4188 (2007).

[23] J. Mathé, H. Visram, V. Viasnoff, Y. Rabin, and A. Meller, Biophys. J. **87**, 3205 (2004).

[24] M. T. Woodside, C. Garcia-Garcia, and S. M. Block, Curr. Opin. Chem. Biol. **12**, 640 (2008).

[25] S. de Lorenzo, M. Ribezzi-Crivellari, J. R. Arias-Gonzalez, S. B. Smith, and F. Ritort, Biophysical Journal **108**, 2854 (2015).

[26] M. Manosas, J. Camunas-Soler, V. Croquette, and F. Ritort, Nature Comm. **8**, 304 (2017).

[27] A. F. Huxley and R. M. Simmons, Nature **233**, 533 (1971).

[28] T. L. Hill, Proc. Nat. Acad. Sci. **70**, 2732 (1973).

[29] M. Epstein and W. Herzog, Theoretical models of skeletal muscle: biological and mathematical considerations (New York, Wiley, 1998).

[30] M. Caruel, J.-M. Allain, and L. Truskinovsky, Phys. Rev. Lett. **110**, 248103 (2013).

[31] M. Caruel, and L. Truskinovsky, Phys. Rev. E **93**, 062407 (2016).

[32] M. Caruel and L. Truskinovsky, Rep. Prog. Phys. **81**, 036602 (2018).

[33] M. Caruel, P. Moireau, D. Chapelle, Biomech. Model. Mechanobiol. **18**, 563 (2019).

[34] M. Peyrard, A. R. Bishop, Phys. Rev. Lett. **62**, 2755 (1989).

[35] N. Theodorakopoulos, M. Peyrard and R. S. MacKay, Phys. Rev. Lett. **93**, 258101 (2004).

[36] P. Grinza, A. Mossa, Phys. Rev. Lett. **92**, 158102 (2004).

[37] M. Peyrard, Nonlinearity **17**, R1 (2004).

[38] V. Ivanov, Y. Zeng, and G. Zocchi, Phys. Rev. E **70**, 051907 (2004).

[39] S. Cocco, J. Yan, J.-F. Léger, D. Chatenay, and J. F. Marko, Phys. Rev. E **70**, 011910 (2004).

[40] J. Palmeri, M. Manghi, and N. Destainville, Phys. Rev. E **77**, 011913 (2008).

[41] A. E. B. Pupo, F. Falo, and A. Fiasconaro, J. Chem. Phys. **139**, 095101 (2013).

[42] F. Manca, S. Giordano, P. L. Palla, and F. Cleri, Phys. Rev. Lett. **113**, 255501 (2014).

[43] G. Perret, T. Lacornerie, F. Manca, S. Giordano, M. Kumemura, N. Lafitte, L. Jalabert, M. C. Tarhan, E. F. Lartigau, F. Cleri, H. Fujita, and D. Collard, Microsystems & Nanoengineering **2**, 16062 (2016).

[44] P.-P. Cortet, M. Ciccotti, and L. Vanel, J. Stat. Mech. **2007**, P03005 (2007).

[45] S. Xia, L. Ponson, G. Ravichandran and K. Bhattacharya, Phys. Rev. Lett. **108**, 196101 (2012).

[46] M.J. Dalbe, P.P. Cortet, M. Ciccotti, L. Vanel, and S. Santucci, Phys. Rev. Lett. **115**, 128301 (2015).

[47] F. Maddalena, D. Percivale, G. Puglisi, L. Truskinovsky, Continuum Mech. Thermodyn. **21**, 251 (2009).

[48] G. Puglisi and L. Truskinovsky, Phys. Rev. E **87**, 032714 (2013).

[49] S. Das, S. Chary, J. Yu, J. Tamelier, K. L. Turner, and J. N. Israelachvili, Langmuir **29**, 15006 (2013).

[50] J. Qian, J. Lin, G.-K. Xu, Y. Lin, H. Gao, J. Mech. Phys. Sol. **101**, 197 (2017).

[51] A. Rafsanjani, A. Akbarzadeh, and D. Pasini, Adv. Mater. **27**, 5931 (2015).

[52] N. Nadkarni, A. F. Arrieta, C. Chong, D. M. Kochmann, and C. Daraio, Phys. Rev. Lett. **116**, 244501 (2016).

[53] M. J. Frazier, D. M. Kochmann, J. Sound. and Vib. **388**, 315 (2017).

[54] S. Katz, S. Givli, Extr. Mech. Lett. **22**, 106 (2018).

[55] S. Katz, S. Givli, Phys. Rev. E **100**, 032209 (2019).

[56] M. Hwang, and A. F. Arrieta, Scientific Reports **8**, 3630 (2018).

[57] R. L. Harne, M. E. Schoemaker, K. W. Wang, Proc. SPIE **9057**, 90570B (2014).

[58] R. L. Harne, M. E. Schoemaker, B. E. Dussault, K. W. Wang, Applied Energy **130**, **148** (2014).

[59] J.L. Ericksen, J. Elasticity **5**, 191 (1975).

[60] I. Müller and P. Villaggio, Arch. Rational Mech. Anal. **65**, 25 (1977).

[61] B. Fedelich and G. Zanzotto, J. Nonlinear Sci. **2**, 319 (1992).

[62] L. Truskinovsky, A. Vainchtein, J. Mech. Phys. Sol. **52**, 1421 (2004).

[63] G. Puglisi and L. Truskinovsky, J. Mech. Phys. Sol. **53**, 655 (2005).

[64] G. Puglisi and L. Truskinovsky, J. Mech. Phys. Sol. **54**, 2060 (2006).

[65] Y. R. Efendiev and L. Truskinovsky, Continuum Mech. Thermodyn. **22**, 679 (2010).

[66] A. Mielke and L. Truskinovsky, Arch. Rational Mech. Anal. **203**, 577 (2012).

[67] I. Benichou, S. Givli, J. Mech. Phys. Sol. **61**, 94 (2013).

[68] M. Caruel, J.-M. Allain, and L. Truskinovsky, J. Mech. Phys. Sol. **76**, 237 (2015).

[69] M. Caruel, L. Truskinovsky, J. Mech. Phys. Sol. **109**, 117 (2017).

[70] L. Truskinovsky, Fracture as a phase transition, in Contemporary research in mechanics and mathematics of materials, Ericksen's Symposium, ed. R. Batra & M. Beatty, CIMNE, Barcelona, pp. 322-332 (1996).

[71] S. Giordano, Soft Matter **13**, 6877 (2017).

[72] M. Benedito and S. Giordano, J. Chem. Phys. **149**, 054901 (2018).

[73] M. Benedito and S. Giordano, Phys. Rev. E **98**, 052146 (2018).

[74] G. Florio, and G. Puglisi, Scientific Reports **9**, 4997 (2019).

[75] L. Bellino, G. Florio and G. Puglisi, Soft Matter **15**, 8680 (2019).

[76] M. Benedito, F. Manca, P. L. Palla and S. Giordano, Phys. Biol. in press, https://doi.org/10.1088/1478-3975/ab97a8 (2020).

[77] R. L. B. Selinger, Z.-G. Wang, W. M. Gelbart, and A.Ben-Shaul, Phys. Rev. A **43**, 4396-4400 (1991).

[78] H. Borja da Rocha, and L. Truskinovsky, Arch. Appl. Mech. **89**, 535-544 (2019).

[79] H. Borja da Rocha, and L. Truskinovsky, Phys. Rev. Lett. **124**, 015501 (2020).

[80] N. Sharon, and H. Lis, Scientific American **268**, 82-98 (1993).

[81] G. Florio, G. Puglisi, in preparation (2020)

[82] E.L. Florin, V.T. Moy, and H.E. Gaub, Science **264**,
415-417 (1994).

[83] J.N. Israelachvili, Intermolecular and Surface Forces (Elsevier, Amsterdam, 2007).

[84] D. Zuckerman and R. Bruinsma, Phys. Rev. Lett. **74**, 3900 (1995).

[85] J.H. Weiner, Statistical Mechanics of Elasticity (Dover Publication Inc., New York, 2002).

[86] R.G. Winkler, Soft Matter **6**, 6183 (2010).

[87] F. Manca, S. Giordano, P. L. Palla, R. Zucca, F. Cleri, and L. Colombo, J. Chem. Phys. **136**, 154906 (2012).

[88] F. Manca, S. Giordano, P. L. Palla, F. Cleri, L. Colombo, J. Chem. Phys. **138**, 157102 (2013).

[89] F. Manca, S. Giordano, P. L. Palla, F. Cleri, L. Colombo, Phys. Rev. E **87**, 032705 (2013).

[90] F. Manca, S. Giordano, P.L. Palla, F. Cleri, F., Phys. A Stat. Mech. its Appl. **395**, 154 (2014).

[91] S. Giordano, Continuum Mech. Thermodyn. **30**, 459 (2018).

[92] L. Bellino L., G. Florio G., G. Puglisi G. In: Proceedings of XXIV AIMETA Conference 2019. Lecture Notes in Mechanical Engineering. Springer, pp. 733-743 (2020).

[93] N. K. Voulgarakis, A. Redondo, A. R. Bishop and K. O. Rasmussen, Phys. Rev. Lett. **96**, 248101 (2006).

[94] K. Kendall, Molecular Adhesion and Its Applications: the Sticky Universe, Kluwer Academic Publishers, New York 2004

[95] J. H. Gibbs, and E. A. DiMarzio, J. Chem. Phys. **30**, 271 (1959).

[96] D. M. Crothers, N. R. Kallenbach, B. H. Zimm, J. Mol. Biol. **11**, 802 (1965).

[97] C. Kittel, Am. J. Phys. **37**, 917 (1969).

[98] K. Nishinari, S. Koide, P.A. Williams, G.O. Phillips, Journal de Physique **51**, 1759 (1990).

[99] G. Florio, G. Puglisi, S. Giordano. On the validity of the zipper hypothesis in the high temperature regimes. Forthcoming (2020).

[100] A. A. Griffith, Phil. Trans. R. Soc. Lond. A **221**, 163-198 (1921).

[101] S.B. Smith, Y. Cui, C. Bustamante. Overstretching B-DNA: the elastic response of individual double-stranded and single-stranded DNA molecules. Science **271** 795-9 (1996).

[102] P. M. McGuiggan, A. Chiche, J. J. Filliben, D. J.

Yarusso, International Journal of Adhesion & Adhesives **28**, 185 (2008).

[103] Z. Peng, C. Wang, L. Chen, S. Chen, Int. J. Sol. Struct. **51**, 4596 (2014).

[104] S. Murray, C. Hillman, M. Pecht, IEEE Trans. Compon. Packag. Technol. **26**, 524 (2003).

[105] F. Bosia, S. Colella, V. Mattoli, B. Mazzolai and N. M. Pugno, RSC Adv. **4**, 25447 (2014).

[106] S. Iliafar, D. Vezenov, and A. Jagota, Langmuir **29**, 1435 (2013).

[107] M. Abramowitz and I. A. Stegun, *Handbook of Mathematical Functions* (Dover Publication Inc., New York, 1970).

[108] R. A. Usmani, Linear Algebra Its Appl. **212/213**, 413 (1994).

[109] R. A. Usmani, Computers Math Applic. **27**, 59 (1994).

[110] N. D. Mermin and H. Wagner, Phys. Rev. Lett. **17**, 1133 (1966).

[111] A. Gelfert and W. Nolting, J. Phys.: Condens. Matter **13**, R505 (2001).

[112] N. Theodorakopoulos, T. Dauxois, and M. Peyrard, Phys. Rev. Lett. **85**, 6 (2000).

[113] P. Benetatos, A. von der Heydt, and A. Zippelius, New J. Phys. **16**, 113037 (2014).

[114] A.M. Skvortsov, L.I. Klushin, F.A.M. Leermakers, J. Chem. Phys. **126**, 024905 (2007).

[115] D.I. Dimitrov, L.I. Klushin, A.M. Skvortsov, A. Milchev, K. Binder, Eur. Phys. J. E **29**, 9 (2009).

[116] A.M. Skvortsov, L.I. Klushin, A.A. Polotsy, K. Binder, Phys. Rev. E **85**, 031803 (2012).

[117] S. Dutta and P. Benetatos, Soft Matter **14**, 6857 (2018).

[118] S. Dutta and P. Benetatos, Soft Matter **16**, 2114 (2020).

[119] A. Milchev, V. G. Rostiashvili, S. Bhattacharya and T. A. Vilgis, Physics Procedia **3**, 1459 (2010).

[120] A. Ivanov, L. I. Klushin and A. M. Skvortsov, Polym. Sci. Ser. A **54**, 602 (2012).

[121] H. A. Kramers, Physica (The Hague) **7**, 284 (1940).

[122] I. Benichou, S. Givli, Phys. Rev. Lett. **114**, 095504 (2015).

[123] I. Benichou, Y. Zhang, O. K. Dudko, S. Givli, J. Mech. Phys. Sol. **95**, 44 (2016).