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

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 onedimensional 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.

The common feature in all these examples is that the observed macroscopic material behavior results from the evolution of the system in a complex multiwells 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. 1(a). 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 multivariables 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 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 phenomenon. For instance, in this class can be inscribed conformational (folded → unfolded) transitions in polymers or protein macromolecules (e.g., β-sheets domain unfolding, α-helix to β-sheet transitions) or martensitic phase changes between different configurations in metallic alloys. The other class corresponds to transitions between unbroken and broken states of breakable units of the system. This process can be reversible, partly reversible or irreversible according to the specific physical phenomenon. Examples of this scheme include unzipping of hairpins, denaturation of macromolecules, fibrillar biological adhesion, cell adhesion, and peeling of films. The (coarse-grained) one-dimensional energy considered in this second case is shown in Fig. 1(b). Here, the unbroken configuration corresponds to a potential well and the broken configuration corresponds to constant energy and zero force. In Ref. [70], the author suggests that this second case can be deduced as a limit of the first one when the natural configuration of the second state is degenerate. Of course the assumption of a noncoercive energy density changes the mathematical structure of the model. In particular, as we show in the following (see the discussion in Sec. 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 studied 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][31][32][33] and macromolecular chains [71][72][73][74][75][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 energy of Fig. 1(a) (continuous line) can be approximated by a biparabolic function (dashed lines) with the switching among the wells described by the spin variable. This strongly facilitates the calculation of the partition function and the corresponding thermodynamic quantities, delivering possible analytical results with a clear physical interpretation. It is important to remark 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 possible 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) negligible 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 adopted to model units with transitions between ground and metastable states [see Fig. 1(a)], the case of units undergoing transition between unbroken and broken states [ Fig. 1(b)] has been investigated, to the best of our knowledge, only in the cases of parallel links [77][78][79]. In the case considered here, the film deposited on a given substrate 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 different 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 examples of biological adhesion and decohesion phenomena the binding enthalpies range from about one k B T (for hydrogen bonds, e.g. controlling the stability of base pairs in DNA) to tens of k B 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 peeling of the film can be induced by prescribing a given extension [Helmholtz ensemble, hard device, Fig. 2 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][86][87][88][89][90][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 → ∞). 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. 1(a)]. Indeed in the case of nondegenerate wells [ Fig. 1(a)] the transition between the two homogeneous (folded↔unfolded) states takes place at a temperature independent force [71][72][73][74][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 = E / x [see Fig. 1(a) for the definition of E and 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 corresponding to the wells of the units. The equilibrium positions can be obtained by ∂U i /∂x = 0 and we get x i = 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 < E / x where x = 2 − 1 . Differently, in the case of breakable links of interest in this paper [with energy function as in Fig. 1(b)], 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. 1(b). Here, we analytically describe this effect and the whole process is explained by means of a phase transition occurring at a given critical temperature. 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 033227-3 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 nonequivalence in the thermodynamic limit, which is an unusual and intriguing behavior.
The paper is organized as it follows. In Sec. II, we introduce the system. In Secs. III and IV, we study the Helmholtz and Gibbs ensembles, respectively. In Sec. V, we study the thermodynamic limit for both ensembles. The conclusions (Sec. 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 onsite 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 leading to a clear physical interpretation. However, 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 ϕ = 1 2 k(y i+1 − y i ) 2 [ Fig. 2(a)], while the vertical ones (elastic constant h) can be broken or unbroken depending on their extension y i [ Fig. 2(b)]. 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 ψ = 1 2 hy 2 i (when |y i | < y M ) and a broken one a contribution equal to ψ = 1 2 hy 2 M (when |y i | > y M ). As anticipated, two different loading conditions are considered. In the first case [ Fig. 2(c)], 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. 2(d)], 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 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. 2(d)] by assuming to have N − ξ broken elements on the right of the chain and ξ 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][96][97][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 ξ belonging to the phase space of the system, and taking its values in the set {0, 1, 2, ..., N}. In this regard, the variable ξ can be considered as a multivalued spin variable.
The aim of this work is to fully analyze the cohesiondecohesion 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. 2(c) (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 ξ of unbroken links. The total potential energy is where y 0 = 0 and y N+1 = y d in the case of imposed displacement on the final element. It is worth noticing that the last term in is not an irrelevant additive constant since it depends implicitly on ξ , which is a variable of the phase space of the system. The assumption of y 0 = 0 reproduces typical boundary conditions in experiments. In addition, 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 can be rearranged by means of the following matrix definition: which is based on the following four submatrices: where is the main nondimensional parameter of the paper, representing the ratio between the elastic constants of vertical and horizontal elastic elements. Moreover, we introduce the vectors v = (0, 0, 0, ..., 0, 1) ∈ R N and y = (y 1 , y 2 , y 3 , ..., y N ) ∈ R N . The energy function can then be rewritten as follows: where y and ξ are the main variables belonging to the phase space of the system. This expression of 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 y, with an additional term independent of y. The partition function of the system analyzed within the Helmholtz ensemble can therefore be written as where β = (k B T ) −1 , k B is the Boltzmann constant and T the absolute temperature. When Eq. (9) is substituted into Eq. (10), we get where we can use the property of the Gaussian integrals holding for any symmetric and positive definite matrix M. Indeed, by considering M = βkA(ξ ) and b = βky d v in Eq. (12), we easily obtain and we have used the definition of 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] which can be rewritten as Another important quantity to describe the decohesion dependence from the temperature is the average value ξ H of unbroken vertical springs. It can be directly evaluated through the expression .
Moreover we can determine the expectation value of the whole position vector y (assigning the deformed configuration of the chain at given displacement): Then, we can use the definition of , thus obtaining the more explicit expression Since, by differentiating Eq. (12) with respect to b we get straightforward calculations give 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.
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 [see Figs. 3(a) and 3(c)] shows the importance of discreteness, exhibiting distinct rupture occurrences in the decreasing steps of the quantity ξ H , and in the peaks of the force-displacement ( f H , y d ) curves along the whole decohesion process. The resulting evolution of the displacements is also reported in Fig. 4(b). 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. 4(a)]. 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 pulloff-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 nonconnected configuration (e.g., completely detached), the partition function corresponds to a nonconvergent 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 Figs. 3(b) and 3(d). 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. 1(a), with respect with the breakable links considered here, is thoroughly studied in Sec. V A where we consider the thermodynamic limit N → ∞.
We remark that the gradual increase of the debonded fraction ξ 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 sawtoothlike 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 (isotensional) ensemble [1], 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.
To conclude, we also observe that the behavior of the forceextension 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 [101][102][103][104][105].

IV. SOFT DEVICE: GIBBS ENSEMBLE
Consider now the case when the film is loaded by a fixed force f [ Fig. 2(d), isotensional condition]. In this case we introduce the Gibbs ensemble. The total energy of the system can now be written as − f y N+1 , where is the energy function introduced within the Helmholtz ensemble in Eq. (9). Thus, the Gibbs partition function can be written as Of course this corresponds to the Laplace transform of the Helmholtz partition function [85]. By using Eqs. (13) and (14), we get where the integral is Gaussian due to the structure of ξ and, therefore, can be easily evaluated. A straightforward 033227-7 calculation leads to where The expected value of the extension y d G = y N+1 G of the last element at given applied force f can then be evaluated as [85] which gives Similarly, we can also determine the average number of unbroken links Eventually, we can evaluate the average value of the displacement of each elements of the chain. This quantity is defined by the expression The comparison of Eq. (29) with Eqs. (18) and (14) yields where the final integral is Gaussian due to the structure of ξ (β, y N+1 ). The calculation can be performed to give the final result in the form By using again Eqs. (A17)-(A19), 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 ξ 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 extensioncontrolled decohesion, obtained with the hard device. This dissimilarity between the 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 noncooperative 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 y G , plotted in Fig. 6(b), while the sequential or noncooperative ruptures were observed in the isometric behavior of y H , plotted in Fig. 4(b). 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 Sec. V B by considering the thermodynamic limit for N → ∞.
In Fig. 6(a) 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 pulloff-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 → ∞), and we obtain explicit analytic results to describe the system behavior in the thermodynamic limit.

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 033227-8 where we introduced the dimensionless parametersβ and where (see Appendix for details).
To begin, we consider Eq. (32), with D given in Eq. (35). In the limit of large N, we introduce the continuum variable x = ξ/N and substitute the summations with integrals where we have considered a generic function φ.
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 [106]. 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 γ (Nx) and γ (Nx + 1). Under the hypothesis of large values of N, it is easy to see that the function γ in Eq. (37), can be approximated by γ (z) Therefore, after introducing we eventually obtain To simplify this result, we can apply the change of variable s = N (1 − x) + ρ, that allows us to prove that a universal force-extension curve exists in the thermodynamic limit (N → ∞) and its shape is given by the formula

033227-10
where it is not difficult to verify that both integrals are well defined provided that a condition thoroughly discussed below. Now, a direct calculation of the two integrals gives where the functions g ± are defined by means of the error function, as Recalling that F H = βy M f H , the force is given by where and, for later convenience, we introduced the critical temperature 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 Y → +∞ and y d → +∞, respectively. Thus, we can obtain the asymptotic value Equivalently, we have 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. 7(a), where the force displacement relation is plotted for different values of the temperature. In particular, we compare the forcedisplacement 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. 7(a), 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 → ∞. However, this representation can be arbitrarily improved by adding more Euler-Mclaurin remainder terms in Eq. (38) that could be substituted with 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 [106]. 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. 7(b), we plot the obtained temperature dependent value of the asymptotic unfolding force for different values of the nondimensional parameter η. We note that the force approaches zero at the critical temperature T c , and that T c is an increasing function of the ratio η. 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 multiwalls 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 f as = √ khy 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 ξ H given in Eq. (33), measuring the average number of unbroken vertical springs. As before, by assuming a large value of N and using Eq. (38), we can substitute the sums with the corresponding integrals and we find We can now apply the change of variable s = N (1 − x) + ρ and get Then, we can determine the expectation number ζ H = N − ξ H of broken links in the thermodynamic limit (N → ∞), eventually obtaining This result proves that for a long chain the number of broken elements must depends 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 → ∞, describing the full 033227-12 It means that we have a linear relation between ζ H and Y for Y → ∞. This expression can be written with the physical parameters of the system as This behavior is confirmed in Fig. 8, where we plotted ζ H calculated through Eq. (17) or Eq. (33) for a discrete system with N = 100, and the quantity ζ H | y d →∞ /y d , obtained from Eq. (57). We observe that the attainment of the vertical asymptotes 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 → 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][35][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. (A17)-(A19), we obtain whereβ is defined in Eq. (36), the adimensional force is given by F = β f y M , and the function γ 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 γ defined in Eq. (37) as γ (z)

033227-13
First, we get By considering the change of variable s = N (1 − x) + ρ, we obtain In the thermodynamic limit (N → ∞), we have Both integrals in Eq. (64) converge ifβ This condition will be thoroughly discussed in the following. By integrating, we eventually obtain the force displacement relation in the thermodynamic limit It is easy to verify that Y G is an increasing function of F and it diverges when the force 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: As before, y d 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 thus interpreting the condition in Eq. (65) or, equivalently, 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 → ∞, obtaining the asymptotic form As before, we can use the change of variable s = N (1 − x) + ρ, which delivers From now on, we consider the average number of broken vertical elements ζ G = N − ξ G given by This result proves that the limit of ζ for N → ∞ exists and can be written as follows: where the critical temperature T c is given in Eq. (49). Coherently with previous results we observe that the number of detached elements monotonically grows with 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 forceextension relation in the corresponding Helmholtz and Gibbs ensembles reported in Fig. 10. In particular, in Fig. 10(a), 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 . Moreover, in Fig. 10(b), we plotted Eqs. (46) and (67), representing the force-extension behavior in the thermodynamic limit. Observe the 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 f 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 only the 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 Fig. 11(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 f as given in Eq. (51). In addition, the values of these Helmholtz force peaks, given by max { f H }, have been plotted in Fig. 11(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 } = f as . This quantity is represented by the blue curve in Fig. 11(b). 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 nonequivalence 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) 033227-15 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 calculations give for the Helmholtz ensemble and the Gibbs ensemble, respectively. Here, for the sake of simplicity, we defined the quantity In terms of physical parameters the effective stiffnesses are The plot of the effective stiffness defined in Eqs. (77) and (78) versus the temperature can be found in Fig. 12(a), from which we can deduce again the nonequivalence 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 Helmholtz with ρ 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. 12(b)]. 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 cohesion/decohesion process related to a film deposited on a substrate, by focusing on the effects of thermal fluctuations. The paradigmatic system adopted is composed of a onedimensional 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 [107,108], 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 analysis of the thermodynamic limit. Indeed, when N → ∞, 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. 11(b).
The existence of finite-temperature phase transitions in low-dimensional many-body models is a subject of large interest in theoretical physics [109,110]. 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 [111]. 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 tensioninduced binding of two parallel semi-flexible polymers [112]. 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. In contrast, the observation of a genuine phase transition in our model is due to important differences as compared with Ref. [112]: the force is transversal in our case whereas the binding tension is longitudinal in Ref. [112]; moreover, we adopted a flexible chain whereas a two semi-flexible chains are introduced in Ref. [112]; finally, our model is discrete, while the model of Ref. [112] 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][87][88][89][90][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 nonequivalence between the defined statistical ensembles [113][114][115][116][117][118]. In particular, in Refs. [113,114] an end-tethered polymer chain compressed between two pistons is considered and shows nonequivalence of the ensembles and a phase transition corresponding to the escape from the gap between the pistons. Reference [115] 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. [116] a Gaussian chain is tethered on a rigid planar surface at one end and the ensemble nonequivalence 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. [117]. Finally, in Ref. [118], the desorption of a self-avoiding polymer chain from a surface has been studied yielding the nonequivalence 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 nonequivalence 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 [116,117,119]. In our system the strong interactions, which are nonlocal and long-range, are generated by the ladder network geometry composed of linear and breakable springs. As discussed in Refs. [116,117,119], the resulting fluctuations of the macroscopic observables cause a nonuniform convergence to the thermodynamic limit eventually producing the ensembles nonequivalence.
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 terms 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 versus 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 nonequivalence 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 [120][121][122]. 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. ACKNOWLEDGMENTS S.G. acknowledges the region "Hauts de France" for the financial support under Project MEPOFIB. G.F. and G.P. 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 Project RAEE SUD-PVP. G.F. is also supported by Istituto Italiano di Fisica Nucleare through Project "QUANTUM" and by a FFABR research grant (MIUR).

APPENDIX: PROPERTIES OF A ξ N (η)
We prove here some properties concerning the matrix 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 [107,108] that the elements of the inverse matrix T −1 can be represented as where the sequences ϑ i and ϕ i are given by the recursive relations and While Eq. (A3) is an increasing recursive relation going from i = 1 to i = N, Eq. (A4) is a decreasing recursive relation going from i = N to i = 1. We also remember that det T = ϑ N [107,108].
In the case of the matrix A ξ N (η) in Eq. (3), we have that b i = c i = −1∀i. Under this hypothesis, the general result can be simplified as follows: where the sequences ϑ i and ϕ i are given by the simplified recursive relations ϑ i = a i ϑ i−1 − ϑ i−2 , i = 2, ..., N, and 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) Therefore, we consider Eq. (A6) with a 1 = · · · = a ξ = 2 + η and a ξ +1 = · · · = a N = 2. Then, for i ξ we have the difference equation ϑ i = (2 + η)ϑ i−1 − ϑ i−2 , whose general solution can be written as with = η 2 + 4η and where the coefficients p and q must be fixed through the condition ϑ 0 = 1 and ϑ 1 = a 1 . A straightforward calculation leads to the explicit solution for i ξ , For i > ξ, we have the simpler difference equation ϑ i = 2ϑ i−1 − ϑ i−2 , with the general solution ϑ i = r + si. In this case, the coefficients r and s must be obtained by imposing ϑ ξ −1 and ϑ ξ by means of Eq. (A11). Hence, the result for i > ξ can be eventually found as where ϑ ξ −1 and ϑ ξ are given by Eq. (A11). The obtained results can be summarized through the final expression holding for 1 i N where the function γ (z) is defined as follows: A more compact form of the solution can be also written as ϑ i = [γ (i + 1) − (i − ξ + 1)γ (ξ + 1) where ϑ i is given in Eq. (A16) and the function γ (z) in Eq. (A15).