# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo, Stefano Giordano

► To cite this version:

Andrea Cannizzo, Stefano Giordano. Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems. Symmetry, 2024, 16 (5), pp.632. 10.3390/sym16050632. hal-04593489

HAL Id: hal-04593489

https://hal.science/hal-04593489

Submitted on 29 May 2024

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
# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo, Stefano Giordano

► To cite this version:

Andrea Cannizzo, Stefano Giordano. Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems. Symmetry, 2024, 16 (5), pp.632. 10.3390/sym16050632. hal-04593489

HAL Id: hal-04593489

https://hal.science/hal-04593489v1

Submitted on 29 May 2024

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
<figure><img src="image_2.png" /></figure>

Article

# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo and Stefano Giordano

Special Issue

Feature Papers in Section "Engineering and Materials" 2024

Edited by

Dr. Vasilis K. Oikonomou

<figure><img src="image_3.png" /></figure>
<figure><img src="image_2.png" /></figure>

Article

# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo and Stefano Giordano

Special Issue

Feature Papers in Section "Engineering and Materials" 2024

Edited by
Dr. Vasilis K. Oikonomou

<figure><img src="image_3.png" /></figure>
Article

# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo and Stefano Giordano

CNRS, Centrale Lille Univ. Polytechnique Hauts-de-France UMR 8520, IEMN, Institut d'Électronique de Microélectronique et de Nanotechnologie, University of Lille, F-59000 Lille, France;
andreacannizzo1992@gmail.com

* Correspondence: stefano.giordano@univ-lille.fr

**Abstract:** Systems with a multistable energy landscape are widespread in physics, biophysics, technology, and materials science. They are strongly influenced by thermal fluctuations and external mechanical actions that can be applied at different rates, moving the system from equilibrium to non-equilibrium regimes. In this paper, we focus on a simple system involving a single breaking phenomenon to describe the various theoretical approaches used to study these problems. To begin with, we propose the exact solution at thermodynamic equilibrium based on the calculation of the partition function without approximations. We then introduce the technique of spin variables, which is able to simplify the treatment even for systems with a large number of coordinates. We then analyze the energy balance of the system to better understand its underlying physics. Finally, we introduce a technique based on transition state theory useful for studying the non-equilibrium dynamical regimes of these systems. This method is appropriate for the evaluation of rate effects and hysteresis loops. These approaches are developed for both the Helmholtz ensemble (prescribed extension) and the Gibbs ensemble (applied force) of statistical mechanics. The symmetry and duality of these two ensembles is discussed in depth. While these techniques are used here for a simple system with theoretical purposes, they can be applied to complex systems of interest for several physical, biophysical, and technological applications.

**Keywords:** statistical mechanics; multistable energy; Gibbs ensemble; Helmholtz ensemble; spin variables; rate equations; thermal fluctuations; transition state theory

## 1. Introduction

Systems with multistable energy landscapes or micro-instabilities play a very important role in several artificial and biological investigations. In these systems, we can observe transitions between different well defined states, induced by mechanical actions or thermal fluctuations. These transitions generate a great complexity in the system response, characterized by thermomechanical coupling and strongly nonlinear behaviors. From the point of view of artificial systems, we can mention the following phenomena characterized by multistable behaviors: the peeling of a film from a substrate [1–4], the waves propagation in bistable lattices [5,6], the energy harvesting through multistable chains [7,8], the plasticity and hysteresis in phase transitions and martensitic transformations of solids [9–13], the cracks and dislocations nucleation and propagation in materials and alloys [14–17], and the friction at the nanoscale [18,19]. On the other hand, from the point of view of the micro-mechanical biological phenomena, we can mention the conformational transitions in polymeric and biopolymeric chains [20–30], the attached and detached states of fibrils in cell adhesion [31–34], the unzipping of macromolecular hairpins [35–38], the sarcomeres behavior in skeletal muscles [39–44], the denaturation or degradation of nucleic acids, polypeptide chains or other polymers [45–54], and the macromolecular or soft matter friction [55–58].

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

check for
updates

**Citation:** Cannizzo, A.; Giordano, S. Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems.

*Symmetry* **2024**, *16*, 632. https://doi.org/10.3390/sym16050632

Academic Editor: Vasilis K. Oikonomou

Received: 12 April 2024

Revised: 7 May 2024

Accepted: 14 May 2024

Published: 20 May 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Article

# Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Andrea Cannizzo and Stefano Giordano

CNRS, Centrale Lille Univ. Polytechnique Hauts-de-France UMR 8520, IEMN, Institut d'Électronique de Microélectronique et de Nanotechnologie, University of Lille, F-59000 Lille, France;
andreacannizzo1992@gmail.com

* Correspondence: stefano.giordano@univ-lille.fr

**Abstract:** Systems with a multistable energy landscape are widespread in physics, biophysics, technology, and materials science. They are strongly influenced by thermal fluctuations and external mechanical actions that can be applied at different rates, moving the system from equilibrium to non-equilibrium regimes. In this paper, we focus on a simple system involving a single breaking phenomenon to describe the various theoretical approaches used to study these problems. To begin with, we propose the exact solution at thermodynamic equilibrium based on the calculation of the partition function without approximations. We then introduce the technique of spin variables, which is able to simplify the treatment even for systems with a large number of coordinates. We then analyze the energy balance of the system to better understand its underlying physics. Finally, we introduce a technique based on transition state theory useful for studying the non-equilibrium dynamical regimes of these systems. This method is appropriate for the evaluation of rate effects and hysteresis loops. These approaches are developed for both the Helmholtz ensemble (prescribed extension) and the Gibbs ensemble (applied force) of statistical mechanics. The symmetry and duality of these two ensembles is discussed in depth. While these techniques are used here for a simple system with theoretical purposes, they can be applied to complex systems of interest for several physical, biophysical, and technological applications.

**Keywords:** statistical mechanics; multistable energy; Gibbs ensemble; Helmholtz ensemble; spin variables; rate equations; thermal fluctuations; transition state theory

## 1. Introduction

Systems with multistable energy landscapes or micro-instabilities play a very important role in several artificial and biological investigations. In these systems, we can observe transitions between different well defined states, induced by mechanical actions or thermal fluctuations. These transitions generate a great complexity in the system response, characterized by thermomechanical coupling and strongly nonlinear behaviors. From the point of view of artificial systems, we can mention the following phenomena characterized by multistable behaviors: the peeling of a film from a substrate [1–4], the waves propagation in bistable lattices [5,6], the energy harvesting through multistable chains [7,8], the plasticity and hysteresis in phase transitions and martensitic transformations of solids [9–13], the cracks and dislocations nucleation and propagation in materials and alloys [14–17], and the friction at the nanoscale [18,19]. On the other hand, from the point of view of the micro-mechanical biological phenomena, we can mention the conformational transitions in polymeric and biopolymeric chains [20–30], the attached and detached states of fibrils in cell adhesion [31–34], the unzipping of macromolecular hairpins [35–38], the sarcomeres behavior in skeletal muscles [39–44], the denaturation or degradation of nucleic acids, polypeptide chains or other polymers [45–54], and the macromolecular or soft matter friction [55–58].

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

check for
updates

**Citation:** Cannizzo, A.; Giordano, S. Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems.

*Symmetry* **2024**, *16*, 632. https://doi.org/10.3390/sym16050632

Academic Editor: Vasilis K. Oikonomou

Received: 12 April 2024

Revised: 7 May 2024

Accepted: 14 May 2024

Published: 20 May 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
All these physical situations are characterized by a multibasin energy landscape and the state of the system can shift from one energy well to another due to thermal fluctuations and externally applied forces. The energy wells correspond to stable or meta-stable configurations, depending on the value of the energy minima. Typically, these systems are constituted by a large number of units and each unit can be in one among the admissible states. The transitions between these states for each unit correspond to the exploration of the energy landscape for the overall system. This dynamics governs the macroscopic behavior of the system and, in particular, its thermomechanical response. From the point of view of modeling these systems and their physical behavior, there are different mathematical techniques that can be implemented, as briefly discussed below.

First of all, we have to distinguish between the models belonging to the equilibrium statistical mechanics and those pertaining to the out-of-equilibrium statistical mechanics [59–61]. Concerning the equilibrium statistical mechanics, a first possibility is the exact calculation of the partition function of the system. This can be done only for quite basic and elementary systems, where the total energy assumes a simple mathematical expression, allowing the integral calculation of the partition functions in closed form. We remark that even the exact mathematical analysis of a double-basin energy potential (without approximations) is rather complicated from the integral calculation point of view.

Therefore, the spin variables approach has been introduced to describe the system in a simpler way. In this case, we introduce additional variables (included within the phase space of the system), which are discrete and behave similar to spin variables in quantum mechanics. These quantities (defined for each unit of the system) assume a discrete finite number of values and allow the identification of the energy basin explored by each element. Consequently, each energy well can be described by a simple spring-like (quadratic) potential, strongly facilitating the mathematical analysis of the problem. From a historical point of view, the first analysis based on a spin variable was performed to model skeletal muscles in biomechanics [39,40]. This method has been used recently to study many phenomena with multistability, including macromolecular conformational transitions [62–71], phase transformations in materials science [72,73], friction on solid and soft substrates [74,75], adhesion processes [76,77], sacrificial bonds biomechanics [78], and fracture mechanics [79]. These results generalize many previous theories dealing with polymer physics, and in particular the stretching of different types of linear chains [60,80–83]. As far as macromolecules are concerned, the theoretical results obtained through the spin approach were largely compared with experimental data obtained by the single molecule force spectroscopy method [84–93], obtaining the important experimental confirmation of the statistical mechanics of small systems. A particularly interesting case is the denaturation of nucleic acids for which theoretical models describe the unzipping of DNA and RNA very well [94–97].

However, the spin variables approach is not adapted to study the dynamic regime, where the time behavior of the actions applied to the system is important for the system evolution. As a matter of fact, when we study the out-of-equilibrium dynamics, the representation based on the sequence of basins and spin discrete variables is not sufficient since the relaxation times of the system depend on the energy barriers between the potential wells [61,98]. This is consistent, e.g., with the Kramers rate formula, originally formulated to deal with chemical reaction rates [99–103]. The consideration of the possible rate effects is particularly important when the stretching of the macromolecular chains, or other nanosystems, is performed with a large traction velocity [104,105]. Therefore, in order to take into consideration the out-of-equilibrium dynamics of the system, we have to introduce the probability to be in each admissible state and we must determine all the rate constants describing the number of transitions between two given states per unit time. These constants depend on the temperature and on the energy barrier between the involved states. These ingredients can be inserted into a system of differential equations, the so-called rate equations, governing the dynamics of the state probabilities [70,71]. Typically, this system can be numerically solved and the solution allows the determination of the
All these physical situations are characterized by a multibasin energy landscape and the state of the system can shift from one energy well to another due to thermal fluctuations and externally applied forces. The energy wells correspond to stable or meta-stable configurations, depending on the value of the energy minima. Typically, these systems are constituted by a large number of units and each unit can be in one among the admissible states. The transitions between these states for each unit correspond to the exploration of the energy landscape for the overall system. This dynamics governs the macroscopic behavior of the system and, in particular, its thermomechanical response. From the point of view of modeling these systems and their physical behavior, there are different mathematical techniques that can be implemented, as briefly discussed below.

First of all, we have to distinguish between the models belonging to the equilibrium statistical mechanics and those pertaining to the out-of-equilibrium statistical mechanics [59–61]. Concerning the equilibrium statistical mechanics, a first possibility is the exact calculation of the partition function of the system. This can be done only for quite basic and elementary systems, where the total energy assumes a simple mathematical expression, allowing the integral calculation of the partition functions in closed form. We remark that even the exact mathematical analysis of a double-basin energy potential (without approximations) is rather complicated from the integral calculation point of view.

Therefore, the spin variables approach has been introduced to describe the system in a simpler way. In this case, we introduce additional variables (included within the phase space of the system), which are discrete and behave similar to spin variables in quantum mechanics. These quantities (defined for each unit of the system) assume a discrete finite number of values and allow the identification of the energy basin explored by each element. Consequently, each energy well can be described by a simple spring-like (quadratic) potential, strongly facilitating the mathematical analysis of the problem. From a historical point of view, the first analysis based on a spin variable was performed to model skeletal muscles in biomechanics [39,40]. This method has been used recently to study many phenomena with multistability, including macromolecular conformational transitions [62–71], phase transformations in materials science [72,73], friction on solid and soft substrates [74,75], adhesion processes [76,77], sacrificial bonds biomechanics [78], and fracture mechanics [79]. These results generalize many previous theories dealing with polymer physics, and in particular the stretching of different types of linear chains [60,80–83]. As far as macromolecules are concerned, the theoretical results obtained through the spin approach were largely compared with experimental data obtained by the single molecule force spectroscopy method [84–93], obtaining the important experimental confirmation of the statistical mechanics of small systems. A particularly interesting case is the denaturation of nucleic acids for which theoretical models describe the unzipping of DNA and RNA very well [94–97].

However, the spin variables approach is not adapted to study the dynamic regime, where the time behavior of the actions applied to the system is important for the system evolution. As a matter of fact, when we study the out-of-equilibrium dynamics, the representation based on the sequence of basins and spin discrete variables is not sufficient since the relaxation times of the system depend on the energy barriers between the potential wells [61,98]. This is consistent, e.g., with the Kramers rate formula, originally formulated to deal with chemical reaction rates [99–103]. The consideration of the possible rate effects is particularly important when the stretching of the macromolecular chains, or other nanosystems, is performed with a large traction velocity [104,105]. Therefore, in order to take into consideration the out-of-equilibrium dynamics of the system, we have to introduce the probability to be in each admissible state and we must determine all the rate constants describing the number of transitions between two given states per unit time. These constants depend on the temperature and on the energy barrier between the involved states. These ingredients can be inserted into a system of differential equations, the so-called rate equations, governing the dynamics of the state probabilities [70,71]. Typically, this system can be numerically solved and the solution allows the determination of the
system evolution (force-extension response and other macroscopic quantities). Although the transition state theories are approximate, they allow a fairly easy and correct study of the combination of rate and temperature effects.

These methodologies can be applied within two different ensembles of the statistical mechanics [62]. Experiments performed at constant applied force (realized by means of extremely soft devices) correspond to the Gibbs statistical ensemble, and experiments performed at prescribed extension (realized with extremely hard devices) are realizations of the Helmholtz statistical ensemble [62]. These two statistical sets are perfectly symmetrical and their properties are in strict duality (think for example of the Legendre transform between their free energies) [60]. In reality, depending on the equivalent stiffness of the device, the stretching experiments correspond to a situation placed in between the Gibbs and Helmholtz ensembles of the statistical mechanics [63,64]. In the thermodynamic limit, when the number of the units of the system is very large, the Helmholtz and Gibbs ensembles may or may not be equivalent depending on whether the force-extension relationship is the same or not. This is a problem largely studied in existing literature [106–113].

These issues are clarified in this work by means of a rather simple system consisting of three different elastic springs, one of which can break, generating bistability phenomena in the model (see Figure 1, where the duality between Gibbs and Helmholtz ensembles can be appreciated). Despite the simplicity of the scheme considered, it allows all aspects to be analyzed in great detail and is, therefore, of particular theoretical and pedagogical importance. For both Gibbs and Helmholtz ensembles (corresponding to isotensional and isometric conditions), we perform the exact analysis of the system by calculating the partition functions and the corresponding free energies in closed form without approximations. This allows us to understand the transition of the system in correspondence to the spring breaking. A simpler approach is to introduce a binary spin variable to describe this transition and the calculations are much simpler, although some approximation is introduced, especially for temperatures that are too high. A better understanding of the underlying physics is obtained by performing an energy balance. This procedure makes it possible to perfectly understand the force and extension thresholds at which the rupture, i.e., the transition, occurs. To conclude, we generalize these equilibrium approaches in order to take into account the real dynamics of the system. We identify the energy barrier between the two states of the model and we apply a transition state theory to predict the rate effects induced by the time evolution of the mechanical actions applied to the system. This allows the study of the combined effects of thermal fluctuations and a pulling speed applied to the system.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Scheme of the rupture models in the isometric Helmholtz ensemble (**left** panel) and in the isotensional Gibbs one (**right** panel). In both ensembles, the model consists of three springs linking together a bottom fixed layer located at $y = 0$ to an upper layer located at $y = Y$. In the Helmholtz case, the upper substrate is fixed at $Y$ while in the Gibbs case, it is pulled with a fixed force $f$. While the springs with elastic constants $l$ and $k$ are linear and unbreakable, the spring with elastic constant $h$ can be intact or broken depending on the applied extension.</figcaption></figure>
system evolution (force-extension response and other macroscopic quantities). Although the transition state theories are approximate, they allow a fairly easy and correct study of the combination of rate and temperature effects.

These methodologies can be applied within two different ensembles of the statistical mechanics [62]. Experiments performed at constant applied force (realized by means of extremely soft devices) correspond to the Gibbs statistical ensemble, and experiments performed at prescribed extension (realized with extremely hard devices) are realizations of the Helmholtz statistical ensemble [62]. These two statistical sets are perfectly symmetrical and their properties are in strict duality (think for example of the Legendre transform between their free energies) [60]. In reality, depending on the equivalent stiffness of the device, the stretching experiments correspond to a situation placed in between the Gibbs and Helmholtz ensembles of the statistical mechanics [63,64]. In the thermodynamic limit, when the number of the units of the system is very large, the Helmholtz and Gibbs ensembles may or may not be equivalent depending on whether the force-extension relationship is the same or not. This is a problem largely studied in existing literature [106–113].

These issues are clarified in this work by means of a rather simple system consisting of three different elastic springs, one of which can break, generating bistability phenomena in the model (see Figure 1, where the duality between Gibbs and Helmholtz ensembles can be appreciated). Despite the simplicity of the scheme considered, it allows all aspects to be analyzed in great detail and is, therefore, of particular theoretical and pedagogical importance. For both Gibbs and Helmholtz ensembles (corresponding to isotensional and isometric conditions), we perform the exact analysis of the system by calculating the partition functions and the corresponding free energies in closed form without approximations. This allows us to understand the transition of the system in correspondence to the spring breaking. A simpler approach is to introduce a binary spin variable to describe this transition and the calculations are much simpler, although some approximation is introduced, especially for temperatures that are too high. A better understanding of the underlying physics is obtained by performing an energy balance. This procedure makes it possible to perfectly understand the force and extension thresholds at which the rupture, i.e., the transition, occurs. To conclude, we generalize these equilibrium approaches in order to take into account the real dynamics of the system. We identify the energy barrier between the two states of the model and we apply a transition state theory to predict the rate effects induced by the time evolution of the mechanical actions applied to the system. This allows the study of the combined effects of thermal fluctuations and a pulling speed applied to the system.

<figure><img src="image_6.png" /><figcaption>**Figure 1.** Scheme of the rupture models in the isometric Helmholtz ensemble (**left** panel) and in the isotensional Gibbs one (**right** panel). In both ensembles, the model consists of three springs linking together a bottom fixed layer located at $y = 0$ to an upper layer located at $y = Y$. In the Helmholtz case, the upper substrate is fixed at $Y$ while in the Gibbs case, it is pulled with a fixed force $f$. While the springs with elastic constants $l$ and $k$ are linear and unbreakable, the spring with elastic constant $h$ can be intact or broken depending on the applied extension.</figcaption></figure>
## 2. Problem Statement

We introduce here a single-unit discrete model in order to highlight the fundamental
differences between the Gibbs and the Helmholtz ensembles studied through the principles
of statistical mechanics at equilibrium and out-of-equilibrium. Firstly, we propose the exact
solution of the problem by calculating the pertinent partition function without approxi-
mations. Moreover, the intrinsic simplicity of the model represents the perfect starting
ground in which we can apply the spin approximation and understand its true meaning
and feasibility limits. To generalize the spin methodology to the out-of-equilibrium regime,
we introduce a transition state theory approach leading to specific rate equations. These
developments are separately performed for both Helmholtz and Gibbs ensembles.

The structure of the model studied (shown in Figure 1) is the same both while consid-
ering the isometric Helmholtz ensemble (left panel in the figure) and while considering
the isotensional Gibbs ensemble (right panel) and it consists of three springs of different
elastic constants *k*, *l*, and *h* arranged in order to link two separate substrates. Springs *k* and
*h* (in the figure in blue and green, respectively) are connected in parallel and link together
the bottom substrate, fixed at *y* = 0, to the lower end of spring *l* (in yellow), located at
*y* = *y*₁, with which they are connected in series. Spring *l*, in turn, connects the latter
parallel springs to an upper layer located at *y* = *Y*. Since *y*₁ will be the only variable of
the system, from now on, we refer to it simply by *y*. The important aspect of this model is
that, while springs *l* and *k* behave like normal linear springs, spring *h* is a breakable one,
meaning that there exists a threshold elongation *Y*<sub>M</sub> after which the *h* spring “breaks” or,
more precisely, its potential energy *E*<sub>b</sub> is constant resulting in a null force for elongation
|*y*| > *Y*<sub>M</sub>. The behavior of this breakable spring *h* is better described through its potential
energy landscape in the left panel of Figure 2. Here, we observe that the breakable spring
behavior is described by a standard parabolic potential when its elongation |*y*| < *Y*<sub>M</sub>
does not exceed the threshold (blue curve), meaning that in this region the spring behaves
like a standard one, and by a constant potential (1/2 *h* *Y*<sub>M</sub><sup>2</sup> when its elongation overcome the
threshold |*y*| > *Y*<sub>M</sub> (red curves), meaning that the resulting force, in this case, is null (see
right panel of Figure 2).

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Potential energy of a breakable spring of elastic constant *h* (**left panel**) and resulting force-extension relation (**right panel**). The quantity *Y*<sub>M</sub> is the elongation after which the spring breaks, resulting in an exerted force equal to zero.</figcaption></figure>

The main difference between the two statistical ensembles is that in the Helmholtz
ensemble, the system is studied with the position of the upper layer is fixed at Y instead, in
the Gibbs ensemble, this upper layer is free to move and a fixed force *f* is applied. In the
following sections, we will analyze this system in both ensembles adopting four different
approaches: firstly, we study the model considering the exact Hamiltonian of the system,
secondly, we approximate the total energy by introducing the spin variable approach,
thirdly, we perform an analysis based on the energy balances of the system, and lastly, we
generalize previous approaches with the rate effects.
## 2. Problem Statement

We introduce here a single-unit discrete model in order to highlight the fundamental
differences between the Gibbs and the Helmholtz ensembles studied through the principles
of statistical mechanics at equilibrium and out-of-equilibrium. Firstly, we propose the exact
solution of the problem by calculating the pertinent partition function without approxi-
mations. Moreover, the intrinsic simplicity of the model represents the perfect starting
ground in which we can apply the spin approximation and understand its true meaning
and feasibility limits. To generalize the spin methodology to the out-of-equilibrium regime,
we introduce a transition state theory approach leading to specific rate equations. These
developments are separately performed for both Helmholtz and Gibbs ensembles.

The structure of the model studied (shown in Figure 1) is the same both while consid-
ering the isometric Helmholtz ensemble (left panel in the figure) and while considering
the isotensional Gibbs ensemble (right panel) and it consists of three springs of different
elastic constants *k*, *l*, and *h* arranged in order to link two separate substrates. Springs *k* and
*h* (in the figure in blue and green, respectively) are connected in parallel and link together
the bottom substrate, fixed at *y* = 0, to the lower end of spring *l* (in yellow), located at
*y* = *y*₁, with which they are connected in series. Spring *l*, in turn, connects the latter
parallel springs to an upper layer located at *y* = *Y*. Since *y*₁ will be the only variable of
the system, from now on, we refer to it simply by *y*. The important aspect of this model is
that, while springs *l* and *k* behave like normal linear springs, spring *h* is a breakable one,
meaning that there exists a threshold elongation *Y*<sub>M</sub> after which the *h* spring “breaks” or,
more precisely, its potential energy *E*<sub>b</sub> is constant resulting in a null force for elongation
|*y*| > *Y*<sub>M</sub>. The behavior of this breakable spring *h* is better described through its potential
energy landscape in the left panel of Figure 2. Here, we observe that the breakable spring
behavior is described by a standard parabolic potential when its elongation |*y*| < *Y*<sub>M</sub>
does not exceed the threshold (blue curve), meaning that in this region the spring behaves
like a standard one, and by a constant potential (1/2 *h* *Y*<sub>M</sub><sup>2</sup> when its elongation overcome the
threshold |*y*| > *Y*<sub>M</sub> (red curves), meaning that the resulting force, in this case, is null (see
right panel of Figure 2).

<figure><img src="image_7.png" /><figcaption>**Figure 2.** Potential energy of a breakable spring of elastic constant *h* (**left panel**) and resulting force-extension relation (**right panel**). The quantity *Y*<sub>M</sub> is the elongation after which the spring breaks, resulting in an exerted force equal to zero.</figcaption></figure>

The main difference between the two statistical ensembles is that in the Helmholtz
ensemble, the system is studied with the position of the upper layer is fixed at Y instead, in
the Gibbs ensemble, this upper layer is free to move and a fixed force f is applied. In the
following sections, we will analyze this system in both ensembles adopting four different
approaches: firstly, we study the model considering the exact Hamiltonian of the system,
secondly, we approximate the total energy by introducing the spin variable approach,
thirdly, we perform an analysis based on the energy balances of the system, and lastly, we
generalize previous approaches with the rate effects.
### 3. Exact Solution within the Helmholtz Ensemble

In this section, we adopt the equilibrium statistical mechanics framework and we study the system embedded in a thermal bath at temperature T. The objective of this investigation is to examine the thermomechanical response at the thermodynamic equilibrium considering the model in the Helmholtz ensemble corresponding to the application of the isometric conditions. Based on the previous premises, we can write the total energy of the system as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + E_b(y), \quad (1) $$

where the first and second addends are the potential energy contributions related to the l and k springs, respectively, while $E_b(y)$ is the exact potential energy contribution of the breakable spring h, defined as:

$$ E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}hY_M^2 & \text{if } |y| \ge Y_M. \end{cases} \quad (2) $$

We notice that, introducing the following characteristic function:

$$ \chi_A(y) = \begin{cases} 1 & \text{if } y \in A, \\ 0 & \text{if } y \notin A, \end{cases} \quad (3) $$

where $A = (-Y_M, Y_M)$, it is possible to write the previous Hamiltonian as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2\chi_A(y) + \frac{1}{2}hY_M^2(1 - \chi_A(y)). \quad (4) $$

We observe that, by considering the particular value that the characteristic function $\chi_A(y)$ assumes, we are able to establish if the system is in the intact configuration, when $\chi_A(y \in A) = 1$ (meaning that the breakable spring h is still intact), or in the broken configuration, when $\chi_A(y \notin A) = 0$ (meaning that the breakable spring h has overcome its maximum threshold $Y_M$). This is an important feature for the characteristic function because it can be used as a helpful tool in the evaluation of some interesting quantities, such as the average number of intact breakable springs, as we will see in the following discussion. Thanks to the total energy, we can evaluate the partition function of the system defined as:

$$ Z_H(Y) = \int_{\mathbb{R}} e^{-\frac{\Phi_H(y)}{k_B T}} dy. \quad (5) $$

Writing explicitly the definition of $\chi_A(y)$ for the energy $\Phi_H(y)$ in the previous integral, we can separate the partition function into three different domains of integration, namely:

$$ \begin{aligned} Z_H(Y) = & \int_{-\infty}^{-Y_M} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dy \\ & + \int_{-Y_M}^{Y_M} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}y^2} dy \\ & + \int_{Y_M}^{+\infty} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dy. \end{aligned} \quad (6) $$

The last three integrals can be easily evaluated thanks to the following generalized Gaussian integral valid for a generic domain of integration $(a, b) \in \mathbb{R}$:

$$ \int_a^b e^{-\alpha x^2 + \beta x} dx = \frac{\sqrt{\pi} e^{\frac{\beta^2}{4\alpha}}}{2\sqrt{\alpha}} \left\{ \phi \left[ \sqrt{\alpha} \left( b - \frac{\beta}{2\alpha} \right) \right] - \phi \left[ \sqrt{\alpha} \left( a - \frac{\beta}{2\alpha} \right) \right] \right\}, \quad (7) $$
### 3. Exact Solution within the Helmholtz Ensemble

In this section, we adopt the equilibrium statistical mechanics framework and we study the system embedded in a thermal bath at temperature T. The objective of this investigation is to examine the thermomechanical response at the thermodynamic equilibrium considering the model in the Helmholtz ensemble corresponding to the application of the isometric conditions. Based on the previous premises, we can write the total energy of the system as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + E_b(y), \quad (1) $$

where the first and second addends are the potential energy contributions related to the l and k springs, respectively, while $E_b(y)$ is the exact potential energy contribution of the breakable spring h, defined as:

$$ E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}hY_M^2 & \text{if } |y| \ge Y_M. \end{cases} \quad (2) $$

We notice that, introducing the following characteristic function:

$$ \chi_A(y) = \begin{cases} 1 & \text{if } y \in A, \\ 0 & \text{if } y \notin A, \end{cases} \quad (3) $$

where $A = (-Y_M, Y_M)$, it is possible to write the previous Hamiltonian as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2\chi_A(y) + \frac{1}{2}hY_M^2(1 - \chi_A(y)). \quad (4) $$

We observe that, by considering the particular value that the characteristic function $\chi_A(y)$ assumes, we are able to establish if the system is in the intact configuration, when $\chi_A(y \in A) = 1$ (meaning that the breakable spring h is still intact), or in the broken configuration, when $\chi_A(y \notin A) = 0$ (meaning that the breakable spring h has overcome its maximum threshold $Y_M$). This is an important feature for the characteristic function because it can be used as a helpful tool in the evaluation of some interesting quantities, such as the average number of intact breakable springs, as we will see in the following discussion. Thanks to the total energy, we can evaluate the partition function of the system defined as:

$$ Z_H(Y) = \int_{\mathbb{R}} e^{-\frac{\Phi_H(y)}{k_B T}} dy. \quad (5) $$

Writing explicitly the definition of $\chi_A(y)$ for the energy $\Phi_H(y)$ in the previous integral, we can separate the partition function into three different domains of integration, namely:

$$ \begin{aligned} Z_H(Y) = & \int_{-\infty}^{-Y_M} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dy \\ & + \int_{-Y_M}^{Y_M} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}y^2} dy \\ & + \int_{Y_M}^{+\infty} e^{-\frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dy. \end{aligned} \quad (6) $$

The last three integrals can be easily evaluated thanks to the following generalized Gaussian integral valid for a generic domain of integration $(a, b) \in \mathbb{R}$:

$$ \int_a^b e^{-\alpha x^2 + \beta x} dx = \frac{\sqrt{\pi} e^{\frac{\beta^2}{4\alpha}}}{2\sqrt{\alpha}} \left\{ \phi \left[ \sqrt{\alpha} \left( b - \frac{\beta}{2\alpha} \right) \right] - \phi \left[ \sqrt{\alpha} \left( a - \frac{\beta}{2\alpha} \right) \right] \right\}, \quad (7) $$
where $\phi(x)$ is the error function defined as:

$$ \phi(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt, \qquad (8) $$

and, in particular, $\phi(\pm\infty) = \pm 1$ and $\phi(-x) = -\phi(x)$. After some long but straightforward calculations, we can write the partition function as:

$$ Z_H(Y) = e^{-\frac{hY_M^2}{2K_B T} - \frac{lkY^2}{2K_B T(l+k)}} \sqrt{\frac{\pi}{\alpha(l+k)}} \left[ \frac{\Gamma(l+k)+2}{2} \right] - e^{-\frac{l(k+h)Y^2}{2K_B T(l+k+h)}} \sqrt{\frac{\pi}{\alpha(l+k+h)}} \left[ \frac{\Gamma(l+k+h)}{2} \right], \qquad (9) $$

where, for the sake of readability, we introduced the two following functions:

$$ \alpha(p) = \frac{p}{2K_B T}, \qquad (10) $$

$$ \Gamma(p) = \phi \left[ \sqrt{\alpha(p)} \left( -Y_M - \frac{lY}{p} \right) \right] - \phi \left[ \sqrt{\alpha(p)} \left( Y_M - \frac{lY}{p} \right) \right]. \qquad (11) $$

Thanks to the partition function, we are now able to evaluate the average value of some interesting quantities. The first quantity we are interested in is the average force applied to the upper substrate, namely $\langle f \rangle$, that is given by:

$$ \langle f \rangle = \frac{1}{Z_H} \int_{\mathbb{R}} f e^{-\frac{\Phi_H(y)}{K_B T}} dy. \qquad (12) $$

Now, if we observe that the force can be derived from the potential energy as $f = d\Phi_H/dY$, it is easy to prove the validity of the following alternative compact version for Equation (12):

$$ \langle f \rangle = -K_B T \frac{\partial \log Z_H}{\partial Y} = \frac{\partial \mathcal{F}}{\partial Y}, \qquad (13) $$

where $\mathcal{F} = -K_B T \log Z_H$ is the Helmholtz free energy. Another interesting quantity is the average number of intact breakable springs that, thanks to the particular relation between the characteristic function $\chi_A(y)$ and the state of the breakable spring $h$, can be calculated as:

$$ \langle \chi_A(y) \rangle = \frac{1}{Z_H} \int_{\mathbb{R}} \chi_A(y) e^{-\frac{\Phi_H(y)}{K_B T}} dy. \qquad (14) $$

By deriving the logarithm of the partition function, $\log Z_H$, with respect to $Y_M$, moreover, we obtain:

$$ -\frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M} = \frac{1}{Z_H} \int_{\mathbb{R}} (1 - \chi_A(y)) e^{-\frac{\Phi_H(y)}{K_B T}} dy = \langle 1 - \chi_A(y) \rangle, \qquad (15) $$

which can be used to rewrite Equation (14) in the following compact version:

$$ \langle \chi_A(y) \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M}. \qquad (16) $$

In Figure 3, we show the average force $\langle f \rangle$ (left panel) and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (right panel) versus the dimensionless elongation $Y/Y_M$ for a system composed of springs of elastic constants identically equal to 1 ($l = k = h = 1$) while also varying the thermal to elastic energy ratio $K_B T/(hY_M^2)$. In the figure, we can observe the main effects of the presence of a breakable spring $h$ in the system. As one can notice from the force-extension graph (left panel of Figure 3), in fact, one of the main effects
where $\phi(x)$ is the error function defined as:

$$ \phi(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt, \qquad (8) $$

and, in particular, $\phi(\pm\infty) = \pm 1$ and $\phi(-x) = -\phi(x)$. After some long but straightforward calculations, we can write the partition function as:

$$ Z_H(Y) = e^{-\frac{hY_M^2}{2K_B T} - \frac{lkY^2}{2K_B T(l+k)}} \sqrt{\frac{\pi}{\alpha(l+k)}} \left[ \frac{\Gamma(l+k)+2}{2} \right] - e^{-\frac{l(k+h)Y^2}{2K_B T(l+k+h)}} \sqrt{\frac{\pi}{\alpha(l+k+h)}} \left[ \frac{\Gamma(l+k+h)}{2} \right], \qquad (9) $$

where, for the sake of readability, we introduced the two following functions:

$$ \alpha(p) = \frac{p}{2K_B T}, \qquad (10) $$

$$ \Gamma(p) = \phi\left[\sqrt{\alpha(p)}\left(-Y_M - \frac{lY}{p}\right)\right] - \phi\left[\sqrt{\alpha(p)}\left(Y_M - \frac{lY}{p}\right)\right]. \qquad (11) $$

Thanks to the partition function, we are now able to evaluate the average value of some interesting quantities. The first quantity we are interested in is the average force applied to the upper substrate, namely $\langle f \rangle$, that is given by:

$$ \langle f \rangle = \frac{1}{Z_H} \int_{\mathbb{R}} f e^{-\frac{\Phi_H(y)}{K_B T}} dy. \qquad (12) $$

Now, if we observe that the force can be derived from the potential energy as $f = d\Phi_H/dY$, it is easy to prove the validity of the following alternative compact version for Equation (12):

$$ \langle f \rangle = -K_B T \frac{\partial \log Z_H}{\partial Y} = \frac{\partial \mathcal{F}}{\partial Y}, \qquad (13) $$

where $\mathcal{F} = -K_B T \log Z_H$ is the Helmholtz free energy. Another interesting quantity is the average number of intact breakable springs that, thanks to the particular relation between the characteristic function $\chi_A(y)$ and the state of the breakable spring $h$, can be calculated as:

$$ \langle \chi_A(y) \rangle = \frac{1}{Z_H} \int_{\mathbb{R}} \chi_A(y) e^{-\frac{\Phi_H(y)}{K_B T}} dy. \qquad (14) $$

By deriving the logarithm of the partition function, $\log Z_H$, with respect to $Y_M$, moreover, we obtain:

$$ -\frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M} = \frac{1}{Z_H} \int_{\mathbb{R}} (1 - \chi_A(y)) e^{-\frac{\Phi_H(y)}{K_B T}} dy = \langle 1 - \chi_A(y) \rangle, \qquad (15) $$

which can be used to rewrite Equation (14) in the following compact version:

$$ \langle \chi_A(y) \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M}. \qquad (16) $$

In Figure 3, we show the average force $\langle f \rangle$ (left panel) and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (right panel) versus the dimensionless elongation $Y/Y_M$ for a system composed of springs of elastic constants identically equal to 1 ($l = k = h = 1$) while also varying the thermal to elastic energy ratio $K_B T/(hY_M^2)$. In the figure, we can observe the main effects of the presence of a breakable spring $h$ in the system. As one can notice from the force-extension graph (left panel of Figure 3), in fact, one of the main effects
of the breakable spring is that of changing the rigidity of the system depending on the specific elongation parameter $Y/Y_M$. In the force-extension graph, we can notice that the slope of the curves changes from a steeper one to a less steep one in correspondence with the rupture of the breakable spring $h$ at a specific total elongation that, for now, we will simply denote as $Y^*/Y_M$. Trying to explain the system behaviors observed in the figures, from a merely mechanical point of view, one could think that the system reaches its rupture point simply when the breakable spring $h$ elongation reaches its breaking threshold $Y_M$. Keeping this in mind, if we consider the system just before its breaking point, we know that the force exerted by the spring $l$ has to be equal to the one exerted by the combined actions of $k$ and $h$ springs, i.e., $l(Y - y) = (k + h)y$ (where $Y$ is the total elongation of the system and $y$ is the elongation of $k$ and $h$ parallel springs). If we then solve the latter identity for $Y$ and we set the breakable spring elongation equal to its breaking threshold, $y = Y_M$, we obtain a result, $Y/Y_M = 3$ (in the case of $l = k = h = 1$), apparently in disagreement with what we observe from the force-extension curves of Figure 3, where the actual breaking point $Y^*/Y_M$ seems to anticipate our merely mechanical result $Y/Y_M = 3$. This apparent incongruity will be fully analyzed and explained in the following discussion thanks to the study of the same problem from an energy-balance point of view rather than relying on the use of the statistical mechanics principles alone. In the right panel of Figure 3, the average number of intact breakable springs $\langle \chi_A(y) \rangle$ versus the adimensional elongation parameter $Y/Y_M$ is shown. In this graph, the average number of intact springs $h$, $\langle \chi_A(y) \rangle$, drops to zero in correspondence to the complete failure of the system, i.e., when the elongation reaches its rupture point $Y^*/Y_M$. We remark again here that in the model used to obtain the results shown in the figure so far, we used three springs of elastic constant equal to one, $l = k = h = 1$, and an elongation threshold for the breakable spring $h$ of $Y_M = 1$.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Behavior of the two-rigidity model with variable thermal to elastic energy ratio $K_B T/(hY_M^2) = \{0, 0.01, 0.02, 0.03\}$ (purple, blue, yellow, and red curves, respectively). Here, the purple curves correspond to the pure mechanical case at a zero temperature. The average force $\langle f \rangle$ and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ are presented versus the dimensionless elongation $Y/Y_M$, where, in both cases, $l = k = h = 1$ and $Y_M = 1$ (in arbitrary units). In both panels, the dimensionless elongation $Y^*/Y_M$ is shown, responsible for the rupture phenomenon.</figcaption></figure>

The two different slopes of the two linear branches in the force-extension curves correspond to the effective elastic constants of the system in the intact and broken configuration and can be evaluated using the following equation:

$$k_{\text{eff}} = \left( \frac{1}{l} + \frac{1}{k+h} \right)^{-1} = \frac{l(k+h)}{l+k+h}. \qquad (17)$$

The effective elastic constant, then, assumes the value $k_{\text{eff}} = \frac{2}{3}$ when $(Y/Y_M) < (Y^*/Y_M)$ and all the three springs are intact ($l = k = h = 1$) and the smaller value $k_{\text{eff}} = \frac{1}{2}$ when $(Y/Y_M) > (Y^*/Y_M)$ in agreement to the fact that, in this second branch, the $h$ spring has overcome its elongation threshold $Y_M$ resulting in its breaking. Furthermore, looking
of the breakable spring is that of changing the rigidity of the system depending on the specific elongation parameter $Y/Y_M$. In the force-extension graph, we can notice that the slope of the curves changes from a steeper one to a less steep one in correspondence with the rupture of the breakable spring $h$ at a specific total elongation that, for now, we will simply denote as $Y^*/Y_M$. Trying to explain the system behaviors observed in the figures, from a merely mechanical point of view, one could think that the system reaches its rupture point simply when the breakable spring $h$ elongation reaches its breaking threshold $Y_M$. Keeping this in mind, if we consider the system just before its breaking point, we know that the force exerted by the spring $l$ has to be equal to the one exerted by the combined actions of $k$ and $h$ springs, i.e., $l(Y - y) = (k + h)y$ (where $Y$ is the total elongation of the system and $y$ is the elongation of $k$ and $h$ parallel springs). If we then solve the latter identity for $Y$ and we set the breakable spring elongation equal to its breaking threshold, $y = Y_M$, we obtain a result, $Y/Y_M = 3$ (in the case of $l = k = h = 1$), apparently in disagreement with what we observe from the force-extension curves of Figure 3, where the actual breaking point $Y^*/Y_M$ seems to anticipate our merely mechanical result $Y/Y_M = 3$. This apparent incongruity will be fully analyzed and explained in the following discussion thanks to the study of the same problem from an energy-balance point of view rather than relying on the use of the statistical mechanics principles alone. In the right panel of Figure 3, the average number of intact breakable springs $\langle\chi_A(y)\rangle$ versus the adimensional elongation parameter $Y/Y_M$ is shown. In this graph, the average number of intact springs $h$, $\langle\chi_A(y)\rangle$, drops to zero in correspondence to the complete failure of the system, i.e., when the elongation reaches its rupture point $Y^*/Y_M$. We remark again here that in the model used to obtain the results shown in the figure so far, we used three springs of elastic constant equal to one, $l = k = h = 1$, and an elongation threshold for the breakable spring $h$ of $Y_M = 1$.

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Behavior of the two-rigidity model with variable thermal to elastic energy ratio $K_B T/(hY_M^2) = \{0, 0.01, 0.02, 0.03\}$ (purple, blue, yellow, and red curves, respectively). Here, the purple curves correspond to the pure mechanical case at a zero temperature. The average force $\langle f \rangle$ and the average number of intact breakable springs $\langle\chi_A(y)\rangle$ are presented versus the dimensionless elongation $Y/Y_M$, where, in both cases, $l = k = h = 1$ and $Y_M = 1$ (in arbitrary units). In both panels, the dimensionless elongation $Y^*/Y_M$ is shown, responsible for the rupture phenomenon.</figcaption></figure>

The two different slopes of the two linear branches in the force-extension curves correspond to the effective elastic constants of the system in the intact and broken configuration and can be evaluated using the following equation:

$$k_{\text{eff}} = \left(\frac{1}{l} + \frac{1}{k+h}\right)^{-1} = \frac{l(k+h)}{l+k+h}. \qquad (17)$$

The effective elastic constant, then, assumes the value $k_{\text{eff}} = \frac{2}{3}$ when $(Y/Y_M) < (Y^*/Y_M)$ and all the three springs are intact ($l = k = h = 1$) and the smaller value $k_{\text{eff}} = \frac{1}{2}$ when $(Y/Y_M) > (Y^*/Y_M)$ in agreement to the fact that, in this second branch, the $h$ spring has overcome its elongation threshold $Y_M$ resulting in its breaking. Furthermore, looking
at both panels of Figure 3, it is interesting to observe that the breaking elongation point $Y^*/Y_M$ is independent of the temperature $T$ while the only effect of the temperature on the system is that of smoothing out the curves. The slope of the curve in the transition region is indeed strongly influenced by temperature. The independence from the temperature of the breaking elongation point $Y^*/Y_M$, controlling the transition of the system from the intact to the broken configuration, is the direct result of the fact that this system is composed of a single unit. It is indeed well known that the temperature does play a crucial role in the transition of the system when this is composed of many units interacting with each other (collective phenomena with phase transition).

4. Helmholtz Ensemble with Spin Approximation

In this section, we adopt an alternative statistical mechanics approach to study the system at equilibrium while introducing the spin variable approximation to deal with the non-convex potential energy of the breakable spring $h$. This means that, instead of considering the exact profile of the potential energy $E_b(y)$ of the breakable spring as defined in Equation (2), we approximate it by introducing a new discrete variable called spin, $S$, in the phase space. This variable can assume only two values depending on the particular state of the spring $h$, namely $S = 1$ when the $h$ spring is intact ($|y| < Y_M$) or $S = 0$ when it is in the broken state ($|y| > Y_M$). This new discrete variable $S$ is then used to approximate the potential energy of the spring $h$ with that of a standard intact linear spring of elastic constant $h$ and elongation $y$, when $S = 1$ (blue dashed curve in the left panel of Figure 4), or with a constant potential energy $\frac{1}{2}hY_M^2$ that corresponds to a null force, when $S = 0$ (red dashed curve in the left panel of Figure 4).

<figure><img src="image_9.png" /><figcaption>Figure 4. Spin approach on the potential energy of a breakable spring of elastic constant $h$ (left panel) and resulting force (right panel). The spin assumes value $S = 1$ when the spring is in the intact conformation and behaves similar to a standard linear spring (blue dashed curves) and $S = 0$ when it is in the broken configuration (red dashed curves). The continuous curves are for reference and correspond to the exact potential energy profile.</figcaption></figure>

It is important to remark that the evaluation of the partition function based on the spin variable approach assumes that for both configurations, broken and intact, all the possible deformations $y$ can be attained by the system. This, in general, could represent a problem in the calculation of the partition function but, as shown numerically in Ref. [62], with typical experimental temperatures, the effects of this approximation can be considered statistically negligible since these artificial configurations (superposition of dashed curves in the left panel of Figure 4) have energy sensibly higher than real configurations (continuous curves of Figure 2). The new definition for the breakable spring contribution to the total potential energy is then described by:

$$E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } S = 1, \\ \frac{1}{2}hY_M^2 & \text{if } S = 0. \end{cases} \qquad (18)$$
at both panels of Figure 3, it is interesting to observe that the breaking elongation point $Y^*/Y_M$ is independent of the temperature $T$ while the only effect of the temperature on the system is that of smoothing out the curves. The slope of the curve in the transition region is indeed strongly influenced by temperature. The independence from the temperature of the breaking elongation point $Y^*/Y_M$, controlling the transition of the system from the intact to the broken configuration, is the direct result of the fact that this system is composed of a single unit. It is indeed well known that the temperature does play a crucial role in the transition of the system when this is composed of many units interacting with each other (collective phenomena with phase transition).

4. Helmholtz Ensemble with Spin Approximation

In this section, we adopt an alternative statistical mechanics approach to study the system at equilibrium while introducing the spin variable approximation to deal with the non-convex potential energy of the breakable spring $h$. This means that, instead of considering the exact profile of the potential energy $E_b(y)$ of the breakable spring as defined in Equation (2), we approximate it by introducing a new discrete variable called spin, $S$, in the phase space. This variable can assume only two values depending on the particular state of the spring $h$, namely $S = 1$ when the $h$ spring is intact ($|y| < Y_M$) or $S = 0$ when it is in the broken state ($|y| > Y_M$). This new discrete variable $S$ is then used to approximate the potential energy of the spring $h$ with that of a standard intact linear spring of elastic constant $h$ and elongation $y$, when $S = 1$ (blu dashed curve in the left panel of Figure 4), or with a constant potential energy $\frac{1}{2}hY_M^2$ that corresponds to a null force, when $S = 0$ (red dashed curve in the left panel of Figure 4).

<figure><img src="image_9.png" /><figcaption>Figure 4. Spin approach on the potential energy of a breakable spring of elastic constant $h$ (left panel) and resulting force (right panel). The spin assumes value $S = 1$ when the spring is in the intact conformation and behaves similar to a standard linear spring (blue dashed curves) and $S = 0$ when it is in the broken configuration (red dashed curves). The continuous curves are for reference and correspond to the exact potential energy profile.</figcaption></figure>

It is important to remark that the evaluation of the partition function based on the spin variable approach assumes that for both configurations, broken and intact, all the possible deformations $y$ can be attained by the system. This, in general, could represent a problem in the calculation of the partition function but, as shown numerically in Ref. [62], with typical experimental temperatures, the effects of this approximation can be considered statistically negligible since these artificial configurations (superposition of dashed curves in the left panel of Figure 4) have energy sensibly higher than real configurations (continuous curves of Figure 2). The new definition for the breakable spring contribution to the total potential energy is then described by:

$$E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } S = 1, \\ \frac{1}{2}hY_M^2 & \text{if } S = 0. \end{cases} \qquad (18)$$
We can now use the spin variable S in the same way we previously used the characteristic function $\chi_A(y)$ in order to write the total potential energy of the system as follows:

$$ \Phi_H(y, S) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2S + \frac{1}{2}hY_M^2(1-S), \qquad (19) $$

depending on y and S (both belonging to the phase space). Observing Equation (19), it is evident the similarity between the spin variable S and the characteristic function $\chi_A(y)$ defined in Equation (3). The crucial difference between the two is that, now, the spin variable S belongs to the phase space. The parallelism between S and $\chi_A(y)$ will become clear once we evaluate the average number of intact breakable spring h for the spin approximation case. To evaluate the partition function of the system, we have to sum the spin variable values:

$$ Z_H(Y) = \sum_{S \in \{0,1\}} \int_{\mathbb{R}} e^{-\frac{\Phi_H(y,S)}{K_B T}} dy. \qquad (20) $$

Writing explicitly the expression for the total energy and summing over the spin variable values, we obtain:

$$ Z_H(Y) = \int_{\mathbb{R}} e^{-\frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dy \\ + \int_{\mathbb{R}} e^{-\frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dy, \qquad (21) $$

which is easily solved, using the Gaussian integral, eventually obtaining:

$$ Z_H(Y) = \sqrt{\frac{2\pi K_B T}{l+k}} e^{-\frac{l k Y^2}{2 K_B T(l+k)}} - \frac{h Y_M^2}{2 K_B T} + \sqrt{\frac{2\pi K_B T}{l+k+h}} e^{-\frac{l(k+h)Y^2}{2 K_B T(l+k+h)}}. \qquad (22) $$

We note that thanks to the introduction of the spin variable approximation, the integrals that appeared in the calculation of the partition function were now evaluated on the entire real axis $\mathbb{R}$ instead of on limited domains of integration, as seen previously in the exact potential energy calculations. This represents a huge simplification for the calculation of the partition function and, consequently, for the evaluation of interesting quantities, such as the average force $\langle f \rangle$, previously introduced in Equation (13), or the average number of intact breakable springs that, now, corresponds exactly to $\langle S \rangle$. As briefly introduced before, we notice that the spin variable is defined exactly as the characteristic function $\chi_A(y)$ introduced in Equation (3) and, since the spin variable appears in the definition of the potential energy $\Phi_H$, in the same way that $\chi_A(y)$ did, it is easy to see that $\langle S \rangle$ is exactly the average number of intact breakable springs h, namely:

$$ \langle S \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M}. \qquad (23) $$

As shown in Figure 5, both the average force $\langle f \rangle$ and the average number of intact breakable springs $\langle S \rangle$ versus the dimensionless elongation $Y/Y_M$ of the system obtained with the spin approximation (dashed curves) are in good agreement with the results obtained with the exact calculations (continuous curves). We notice that the spin approximation is more precise for low temperatures. Furthermore, we observe that even in this case the rupture elongation does not correspond to the one someone would expect, but, instead, is shifted as in the exact case. We remark again that the spin variable approximation provides a simpler and quicker way to evaluate the partition function of the system, which is even more useful while evaluating partition functions for a system made of several units as discussed in several investigations concerning friction, adhesion, and fracture processes (see Introduction). In the following section, thanks to the energy-balance approach, we will understand the reason behind the apparent discrepancy between the real
We can now use the spin variable $S$ in the same way we previously used the characteristic function $\chi_A(y)$ in order to write the total potential energy of the system as follows:

$$ \Phi_H(y, S) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2S + \frac{1}{2}hY_M^2(1-S), \quad (19) $$

depending on $y$ and $S$ (both belonging to the phase space). Observing Equation (19), it is evident the similarity between the spin variable $S$ and the characteristic function $\chi_A(y)$ defined in Equation (3). The crucial difference between the two is that, now, the spin variable $S$ belongs to the phase space. The parallelism between $S$ and $\chi_A(y)$ will become clear once we evaluate the average number of intact breakable spring $h$ for the spin approximation case. To evaluate the partition function of the system, we have to sum the spin variable values:

$$ Z_H(Y) = \sum_{S \in \{0,1\}} \int_{\mathbb{R}} e^{-\frac{\Phi_H(y,S)}{K_B T}} dy. \quad (20) $$

Writing explicitly the expression for the total energy and summing over the spin variable values, we obtain:

$$ Z_H(Y) = \int_{\mathbb{R}} e^{-\frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dy \\ + \int_{\mathbb{R}} e^{-\frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dy, \quad (21) $$

which is easily solved, using the Gaussian integral, eventually obtaining:

$$ Z_H(Y) = \sqrt{\frac{2\pi K_B T}{l+k}} e^{-\frac{l k Y^2}{2 K_B T(l+k)}} - \frac{h Y_M^2}{2 K_B T} + \sqrt{\frac{2\pi K_B T}{l+k+h}} e^{-\frac{l(k+h)Y^2}{2 K_B T(l+k+h)}}. \quad (22) $$

We note that thanks to the introduction of the spin variable approximation, the integrals that appeared in the calculation of the partition function were now evaluated on the entire real axis $\mathbb{R}$ instead of on limited domains of integration, as seen previously in the exact potential energy calculations. This represents a huge simplification for the calculation of the partition function and, consequently, for the evaluation of interesting quantities, such as the average force $\langle f \rangle$, previously introduced in Equation (13), or the average number of intact breakable springs that, now, corresponds exactly to $\langle S \rangle$. As briefly introduced before, we notice that the spin variable is defined exactly as the characteristic function $\chi_A(y)$ introduced in Equation (3) and, since the spin variable appears in the definition of the potential energy $\Phi_H$, in the same way that $\chi_A(y)$ did, it is easy to see that $\langle S \rangle$ is exactly the average number of intact breakable springs $h$, namely:

$$ \langle S \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_H}{\partial Y_M}. \quad (23) $$

As shown in Figure 5, both the average force $\langle f \rangle$ and the average number of intact breakable springs $\langle S \rangle$ versus the dimensionless elongation $Y/Y_M$ of the system obtained with the spin approximation (dashed curves) are in good agreement with the results obtained with the exact calculations (continuous curves). We notice that the spin approximation is more precise for low temperatures. Furthermore, we observe that even in this case the rupture elongation does not correspond to the one someone would expect, but, instead, is shifted as in the exact case. We remark again that the spin variable approximation provides a simpler and quicker way to evaluate the partition function of the system, which is even more useful while evaluating partition functions for a system made of several units as discussed in several investigations concerning friction, adhesion, and fracture processes (see Introduction). In the following section, thanks to the energy-balance approach, we will understand the reason behind the apparent discrepancy between the real
breaking point of the system $Y^*/Y_M$ and the one wrongfully predicted adopting a merely mechanistic approach.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Comparison between the spin-approximated quantities (dashed curves) and the corresponding exact results (continuous curves) for $\langle f \rangle$ and $\langle S \rangle$ versus the adimensional elongation $Y/Y_M$. In both cases, the elastic constants are $l = k = h = 1$ and $Y_M = 1$, while the thermal to elastic energy ratio assumes the values $K_B T/(hY_M^2) = \{0.01, 0.2\}$ (respectively, in blue and in red).</figcaption></figure>

**5. Helmholtz Ensemble from an Energetic Perspective**

In this section, we study the model in the isometric configuration adopting the principles of energy minimization, which are valid under the zero temperature assumptions. By doing so, we can understand why it looks like the system tries to anticipate its breaking point in the force-extension curves seen before in both the exact and the approximated calculation cases (see Figures 3 and 5).

In order to understand when the change in rigidity occurs, we need to understand in which state the system prefers to be from an energetic point of view given a prescribed elongation $Y$. The total potential energy of the system is as follows:

$$ \Phi_H(y) = \begin{cases} \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2 & \text{if } |y| > Y_M. \end{cases} \qquad (24) $$

Given a prescribed elongation $Y$, we search for the value of $y$, $y_{min}$, which minimizes the potential energy. We start considering the $|y| < Y_M$ case, where the potential energy is defined as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2. \qquad (25) $$

Deriving $\Phi_H(y)$ by $y$ and setting it equal to zero, $\frac{\partial \Phi_H(y)}{\partial y} = 0$, we obtain the following minimizing value for $y$:

$$ y_{min} = \frac{l}{l + k + h} Y, \qquad (26) $$

with the corresponding minimal potential energy:

$$ \Phi_H(y_{min}) = \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^2. \qquad (27) $$

If we go through the same process for the second range, $|y| > Y_M$, where the potential energy is defined as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2, \qquad (28) $$

eventually, we obtain the following minimizing value for $y$:

$$ y_{min} = \frac{l}{l + k} Y, \qquad (29) $$
breaking point of the system $Y^*/Y_M$ and the one wrongfully predicted adopting a merely mechanic approach.

<figure><img src="image_10.png" /><figcaption>**Figure 5.** Comparison between the spin-approximated quantities (dashed curves) and the corresponding exact results (continuous curves) for $\langle f \rangle$ and $\langle S \rangle$ versus the adimensional elongation $Y/Y_M$. In both cases, the elastic constants are $l = k = h = 1$ and $Y_M = 1$, while the thermal to elastic energy ratio assumes the values $K_B T/(hY_M^2) = \{0.01, 0.2\}$ (respectively, in blue and in red).</figcaption></figure>

## 5. Helmholtz Ensemble from an Energetic Perspective

In this section, we study the model in the isometric configuration adopting the principles of energy minimization, which are valid under the zero temperature assumptions. By doing so, we can understand why it looks like the system tries to anticipate its breaking point in the force-extension curves seen before in both the exact and the approximated calculation cases (see Figures 3 and 5).

In order to understand when the change in rigidity occurs, we need to understand in which state the system prefers to be from an energetic point of view given a prescribed elongation $Y$. The total potential energy of the system is as follows:

$$ \Phi_H(y) = \begin{cases} \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2 & \text{if } |y| > Y_M. \end{cases} \qquad (24) $$

Given a prescribed elongation $Y$, we search for the value of $y$, $y_{min}$, which minimizes the potential energy. We start considering the $|y| < Y_M$ case, where the potential energy is defined as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2. \qquad (25) $$

Deriving $\Phi_H(y)$ by $y$ and setting it equal to zero, $\frac{\partial \Phi_H(y)}{\partial y} = 0$, we obtain the following minimizing value for $y$:

$$ y_{min} = \frac{l}{l + k + h} Y, \qquad (26) $$

with the corresponding minimal potential energy:

$$ \Phi_H(y_{min}) = \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^2. \qquad (27) $$

If we go through the same process for the second range, $|y| > Y_M$, where the potential energy is defined as:

$$ \Phi_H(y) = \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2, \qquad (28) $$

eventually, we obtain the following minimizing value for $y$:

$$ y_{min} = \frac{l}{l + k} Y, \qquad (29) $$
with the corresponding minimal potential energy:

$$ \Phi_H(y_{\min}) = \frac{1}{2} \frac{lk}{l+k} Y^2 + \frac{1}{2} h Y_M^2. \qquad (30) $$

Now, by comparing when the first minimal potential energy, valid for the $|y| < Y_M$ case, is greater than the minimal potential energy valid for $|y| > Y_M$, we can find the total elongation $Y^*$ that defines the threshold in which the system switches from the hard-rigidity to the soft-rigidity configuration:

$$ \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^{*2} > \frac{1}{2} \frac{lk}{l+k} Y^{*2} + \frac{1}{2} h Y_M^2, \qquad (31) $$

which, after some simple calculations, eventually gives:

$$ Y^* > \sqrt{\frac{(l+k)(l+k+h)}{l^2}} Y_M. \qquad (32) $$

This value is in agreement with the results for the average force $\langle f \rangle$ and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (and $\langle S \rangle$) shown in Figure 3 (and Figure 5), where we can confirm that, given a system where $l = k = h = 1$ and $Y_M = 1$, the rupture elongation parameter is equal to $Y^*/Y_M = \sqrt{6}$. To conclude, this energetic approach (without thermal effects) is in agreement with both the exact and approximated approaches and offers a useful tool to be able to evaluate the correct threshold point of the entire system $Y^*$. Moreover, this energetic introduction is useful to study the out-of-equilibrium behavior of the system, as discussed in the next section.

## 6. Rate Effect on the Rupture Model under Helmholtz Conditions

In practical applications, systems with multistable energies are often used by applying external time-varying mechanical actions, and thus, often work out of thermodynamic equilibrium. It is, therefore, important to know how to study the combined rate and temperature effects. To understand the out-of-equilibrium dynamics of our rupture model under Helmholtz conditions, we observe the shape of the energy profiles defined in Equation (24). Of course, we can see an energy profile for $|y| < Y_M$ (intact spring) and another one for $|y| > Y_M$ (broken spring). It is interesting to compare these energy landscapes for different values of $Y$ and $y$. This can be found in Figure 6, where we can observe how these shapes change with the variable applied extension $Y$. In the previous section, we determined the threshold value $Y^*$ corresponding to the transition (i.e., the breaking) point. We see in Figure 6 that for $Y \ll Y^*$ and for $Y \gg Y^*$, the minimum values of the two energy profiles (intact and broken system) are significantly different and the two curves do not generate a potential barrier (the state of the system is well identified); this is true for example for $Y = 1$ and $Y = 4$, shown in the left panel of Figure 6. Otherwise, we see in this same panel that for $Y = Y^*$ the two potential minima are equal to each other and a potential barrier is generated between the two energy profiles. This is exactly the transition point that governs the dynamics of the system as we shall see in a moment. When the values of the applied extension $Y$ are close to $Y^*$, it can be seen that there is always an energy barrier between the two profiles and that the energy minima may be slightly different, as shown in the right panel of Figure 6. The values of the energy barrier $B$ is obtained by substituting $y = Y_M$ in one of the two profiles in Equation (24) for $|y| < Y_M$ or $|y| > Y_M$. The result is as follows:

$$ B(Y) = \Phi_H(Y_M) = \frac{1}{2} l (Y - Y_M)^2 + \frac{1}{2} k Y_M^2 + \frac{1}{2} h Y_M^2. \qquad (33) $$
with the corresponding minimal potential energy:

$$ \Phi_H(y_{\min}) = \frac{1}{2} \frac{lk}{l+k} Y^2 + \frac{1}{2} h Y_M^2 \quad (30) $$

Now, by comparing when the first minimal potential energy, valid for the $|y| < Y_M$ case, is greater than the minimal potential energy valid for $|y| > Y_M$, we can find the total elongation $Y^*$ that defines the threshold in which the system switches from the hard-rigidity to the soft-rigidity configuration:

$$ \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^{*2} > \frac{1}{2} \frac{lk}{l+k} Y^{*2} + \frac{1}{2} h Y_M^2 \quad (31) $$

which, after some simple calculations, eventually gives:

$$ Y^* > \sqrt{\frac{(l+k)(l+k+h)}{l^2}} Y_M \quad (32) $$

This value is in agreement with the results for the average force $\langle f \rangle$ and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (and $\langle S \rangle$) shown in Figure 3 (and Figure 5), where we can confirm that, given a system where $l = k = h = 1$ and $Y_M = 1$, the rupture elongation parameter is equal to $Y^*/Y_M = \sqrt{6}$. To conclude, this energetic approach (without thermal effects) is in agreement with both the exact and approximated approaches and offers a useful tool to be able to evaluate the correct threshold point of the entire system $Y^*$. Moreover, this energetic introduction is useful to study the out-of-equilibrium behavior of the system, as discussed in the next section.

## 6. Rate Effect on the Rupture Model under Helmholtz Conditions

In practical applications, systems with multistable energies are often used by applying external time-varying mechanical actions, and thus, often work out of thermodynamic equilibrium. It is, therefore, important to know how to study the combined rate and temperature effects. To understand the out-of-equilibrium dynamics of our rupture model under Helmholtz conditions, we observe the shape of the energy profiles defined in Equation (24). Of course, we can see an energy profile for $|y| < Y_M$ (intact spring) and another one for $|y| > Y_M$ (broken spring). It is interesting to compare these energy landscapes for different values of $Y$ and $y$. This can be found in Figure 6, where we can observe how these shapes change with the variable applied extension $Y$. In the previous section, we determined the threshold value $Y^*$ corresponding to the transition (i.e., the breaking) point. We see in Figure 6 that for $Y \ll Y^*$ and for $Y \gg Y^*$, the minimum values of the two energy profiles (intact and broken system) are significantly different and the two curves do not generate a potential barrier (the state of the system is well identified); this is true for example for $Y = 1$ and $Y = 4$, shown in the left panel of Figure 6. Otherwise, we see in this same panel that for $Y = Y^*$ the two potential minima are equal to each other and a potential barrier is generated between the two energy profiles. This is exactly the transition point that governs the dynamics of the system as we shall see in a moment. When the values of the applied extension $Y$ are close to $Y^*$, it can be seen that there is always an energy barrier between the two profiles and that the energy minima may be slightly different, as shown in the right panel of Figure 6. The values of the energy barrier $B$ is obtained by substituting $y = Y_M$ in one of the two profiles in Equation (24) for $|y| < Y_M$ or $|y| > Y_M$. The result is as follows:

$$ B(Y) = \Phi_H(Y_M) = \frac{1}{2} l (Y - Y_M)^2 + \frac{1}{2} k Y_M^2 + \frac{1}{2} h Y_M^2 \quad (33) $$
This is the value assumed by the two energy profiles at their intersection. Therefore, it represents a true barrier only for values of Y close to Y*. More precisely, we can see that there is a barrier only if:

$$ \frac{l+k}{l} Y_M < Y < \frac{l+k+h}{l} Y_M. \qquad (34) $$

In addition to the value of the energy barrier, it is interesting to consider the energy of the two minima, corresponding to the two profiles for $|y| < Y_M$ and for $|y| > Y_M$. We get $\Phi_1$ from Equation (27):

$$ \Phi_1(Y) = \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^2, \qquad (35) $$

for $|y| < Y_M$, and $\Phi_2$ from Equation (30):

$$ \Phi_2(Y) = \frac{1}{2} \frac{lk}{l+k} Y^2 + \frac{1}{2} h Y_M^2, \qquad (36) $$

for $|y| > Y_M$. Given the three parameters $B(Y)$, $\Phi_1(Y)$, and $\Phi_2(Y)$, we can define the two following kinetic coefficients describing the dynamics of the system:

$$ k_{12} = r_0 e^{-\frac{B(Y)-\Phi_1(Y)}{k_B T}}, \qquad (37) $$

$$ k_{21} = r_0 e^{-\frac{B(Y)-\Phi_2(Y)}{k_B T}}, \qquad (38) $$

where $r_0$ is a prefactor measured in s⁻¹, which can be obtained through the original Kramers formula or one of its generalizations [98,99]. Here, it is considered as an additional parameter of the system. The two parameters $k_{12}$ and $k_{21}$ represent the number of transitions per unit time induced by thermal fluctuations from state 1 to state 2, and from state 2 to state 1, respectively.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** (Left panel) Helmholtz energy profiles with different values of Y and y. Red lines correspond to the energies $\Phi_H$ for $|y| < Y_M$, while blue lines correspond to $\Phi_H$ for $|y| > Y_M$. (Right panel) zoom with three values of Y closer to Y = √6, useful to identify the barriers B. In all cases, the elastic constants are $l = k = h = 1$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

We can now define the probability $p_1$ to be in the intact state 1 and the probability $p_2$ to be in the broken state 2 of the system. These two probabilities evolve in time by means of the following system of differential Equations [61]:

$$ \frac{dp_1(t)}{dt} = -k_{12}p_1(t) + k_{21}p_2(t), \qquad (39) $$

$$ \frac{dp_2(t)}{dt} = -k_{21}p_2(t) + k_{12}p_1(t), \qquad (40) $$
This is the value assumed by the two energy profiles at their intersection. Therefore, it represents a true barrier only for values of Y close to Y*. More precisely, we can see that there is a barrier only if:

$$ \frac{l+k}{l} Y_M < Y < \frac{l+k+h}{l} Y_M. \qquad (34) $$

In addition to the value of the energy barrier, it is interesting to consider the energy of the two minima, corresponding to the two profiles for $|y| < Y_M$ and for $|y| > Y_M$. We get $\Phi_1$ from Equation (27):

$$ \Phi_1(Y) = \frac{l}{2} \left( \frac{k+h}{l+k+h} \right) Y^2, \qquad (35) $$

for $|y| < Y_M$, and $\Phi_2$ from Equation (30):

$$ \Phi_2(Y) = \frac{1}{2} \frac{lk}{l+k} Y^2 + \frac{1}{2} h Y_M^2, \qquad (36) $$

for $|y| > Y_M$. Given the three parameters $B(Y)$, $\Phi_1(Y)$, and $\Phi_2(Y)$, we can define the two following kinetic coefficients describing the dynamics of the system:

$$ k_{12} = r_0 e^{-\frac{B(Y)-\Phi_1(Y)}{k_B T}}, \qquad (37) $$

$$ k_{21} = r_0 e^{-\frac{B(Y)-\Phi_2(Y)}{k_B T}}, \qquad (38) $$

where $r_0$ is a prefactor measured in s⁻¹, which can be obtained through the original Kramers formula or one of its generalizations [98,99]. Here, it is considered as an additional parameter of the system. The two parameters $k_{12}$ and $k_{21}$ represent the number of transitions per unit time induced by thermal fluctuations from state 1 to state 2, and from state 2 to state 1, respectively.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** (Left panel) Helmholtz energy profiles with different values of Y and y. Red lines correspond to the energies $\Phi_H$ for $|y| < Y_M$, while blue lines correspond to $\Phi_H$ for $|y| > Y_M$. (Right panel) zoom with three values of Y closer to Y = √6, useful to identify the barriers B. In all cases, the elastic constants are $l = k = h = 1$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

We can now define the probability $p_1$ to be in the intact state 1 and the probability $p_2$ to be in the broken state 2 of the system. These two probabilities evolve in time by means of the following system of differential Equations [61]:

$$ \frac{dp_1(t)}{dt} = -k_{12}p_1(t) + k_{21}p_2(t), \qquad (39) $$

$$ \frac{dp_2(t)}{dt} = -k_{21}p_2(t) + k_{12}p_1(t), \qquad (40) $$
where the trajectory $Y(t)$, characterizing $B(Y)$, $\Phi_1(Y)$, and $\Phi_2(Y)$, must be imposed to obtain the system dynamics. The Helmholtz partition function obtained trough the spin approach, see Equation (22), can be written as $Z_H = Z_{H,1} + Z_{H,2}$, where:

$$Z_{H,1}(Y) = \sqrt{\frac{2\pi K_B T}{l + k + h}} e^{-\frac{l(k+h)Y^2}{2K_B T(l+k+h)}}, \quad (41)$$

$$Z_{H,2}(Y) = \sqrt{\frac{2\pi K_B T}{l + k}} e^{-\frac{lkY^2}{2K_B T(l+k)} - \frac{hY^2}{2K_B T}}, \quad (42)$$

from which we obtain the force-extension response pertinent to each state:

$$\langle f \rangle_1 = -K_B T \frac{\partial \log Z_{H,1}(Y)}{\partial Y} = \frac{l(k+h)}{l+k+h} Y, \quad (43)$$

$$\langle f \rangle_2 = -K_B T \frac{\partial \log Z_{H,2}(Y)}{\partial Y} = \frac{lk}{l+k} Y. \quad (44)$$

We can finally determine the overall force-extension response of the system by averaging previous relations with the probabilities $p_1$ and $p_2$:

$$\langle f \rangle(t) = p_1(t) \frac{l(k+h)}{l+k+h} Y(t) + p_2(t) \frac{lk}{l+k} Y(t). \quad (45)$$

Finally, by solving the rate equations stated in Equations (39) and (40), we can obtain the force-extension response of the system. An example of solution for $p_1$ and $p_2$ is shown in Figure 7, where we considered a traction with constant velocity $v$, i.e., $Y = vt$ ($v$ is measured in m/s). We remark that different velocities generate different dynamical regimes for the system response. This effect on the force-extension behavior can be found in Figure 8. Importantly, we note that the force for producing a transition (in this case, the breaking of the spring) is always an increasing function of the traction velocity.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** (Left panel) Helmholtz probability $p_1$ versus $Y$. (Right panel) Helmholtz probability $p_2$ versus $Y$. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The dynamical model proposed can also be used to obtain the hysteretic behavior of the Helmholtz system. An example is shown in Figure 9, where one can find the response of the system to a sinusoidal applied extension of the type $Y(t) = A \sin(\omega t)$. We can observe both the evolution of the probability and of the force-extension curve. In both cases, we remark a hysteretic behavior generated by the memory effect described by the rate constants and rate equations. We plotted different curves corresponding to different values of the extension amplitude. It can be seen that for small amplitudes the transition is only just started and not finished, while with larger amplitudes, the transition can be completed. In any case, the forward and reverse paths are always different because of the memory effect. It is important to remark that these two paths become increasingly similar if
where the trajectory $Y(t)$, characterizing $B(Y)$, $\Phi_1(Y)$, and $\Phi_2(Y)$, must be imposed to obtain the system dynamics. The Helmholtz partition function obtained trough the spin approach, see Equation (22), can be written as $Z_H = Z_{H,1} + Z_{H,2}$, where:

$$Z_{H,1}(Y) = \sqrt{\frac{2\pi K_B T}{l+k+h}} e^{-\frac{l(k+h)Y^2}{2K_B T(l+k+h)}}, \quad (41)$$

$$Z_{H,2}(Y) = \sqrt{\frac{2\pi K_B T}{l+k}} e^{-\frac{lkY^2}{2K_B T(l+k)} - \frac{hY^2}{2K_B T}}, \quad (42)$$

from which we obtain the force-extension response pertinent to each state:

$$\langle f \rangle_1 = -K_B T \frac{\partial \log Z_{H,1}(Y)}{\partial Y} = \frac{l(k+h)}{l+k+h} Y, \quad (43)$$

$$\langle f \rangle_2 = -K_B T \frac{\partial \log Z_{H,2}(Y)}{\partial Y} = \frac{lk}{l+k} Y. \quad (44)$$

We can finally determine the overall force-extension response of the system by averaging previous relations with the probabilities $p_1$ and $p_2$:

$$\langle f \rangle(t) = p_1(t) \frac{l(k+h)}{l+k+h} Y(t) + p_2(t) \frac{lk}{l+k} Y(t). \quad (45)$$

Finally, by solving the rate equations stated in Equations (39) and (40), we can obtain the force-extension response of the system. An example of solution for $p_1$ and $p_2$ is shown in Figure 7, where we considered a traction with constant velocity $v$, i.e., $Y = vt$ ($v$ is measured in m/s). We remark that different velocities generate different dynamical regimes for the system response. This effect on the force-extension behavior can be found in Figure 8. Importantly, we note that the force for producing a transition (in this case, the breaking of the spring) is always an increasing function of the traction velocity.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** (Left panel) Helmholtz probability $p_1$ versus $Y$. (Right panel) Helmholtz probability $p_2$ versus $Y$. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The dynamical model proposed can also be used to obtain the hysteretic behavior of the Helmholtz system. An example is shown in Figure 9, where one can find the response of the system to a sinusoidal applied extension of the type $Y(t) = A \sin(\omega t)$. We can observe both the evolution of the probability and of the force-extension curve. In both cases, we remark a hysteretic behavior generated by the memory effect described by the rate constants and rate equations. We plotted different curves corresponding to different values of the extension amplitude. It can be seen that for small amplitudes the transition is only just started and not finished, while with larger amplitudes, the transition can be completed. In any case, the forward and reverse paths are always different because of the memory effect. It is important to remark that these two paths become increasingly similar if
we reduce the value of $\omega$. Indeed, if $\omega$ is quite small, the rate effects are negligible and the forward and reverse trajectories must asymptotically coincide. Recall also that the effect of velocity observed in the uniform traction case is to increase the force required to impose the Helmholtz constraint on the final position of the system. In the case of the periodic driving, we observe a similar phenomenon and the force is an increasing function of the applied frequency.

<figure><img src="image_13.png" /><figcaption>**Figure 8.** (Left panel) Helmholtz force-extension response of the system with different traction velocities. (Right panel) zoom on the transition region. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1$, $r_0 = 120$, $K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Hysteretic behavior of the Helmholtz model. (Left panel) probability $p_2$ versus the prescribed extension $Y$ (we note that $p_1 + p_2 = 1$). (Right panel) average force $\langle f \rangle$ versus extension $Y$. We used an applied extension $Y(t) = A \sin(\omega t)$ with different amplitudes $A = 2.45, 2.5, 2.55, 2.6, 2.7, 3, 3.5$ and $\omega = 2.5$ (arbitrary units). In all cases, we adopted the parameters $l = k = h = 1$, $r_0 = 120$, $K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The approach discussed is based on the Kramers formula (or similar expressions), giving the rates of the transitions between the two states in terms of the parameters of the system. There are, therefore, several limitations, as discussed below. First of all, the Kramers formula has been obtained through a steady-state assumption (slow quasi-stationary escape of particles over the barrier), and therefore, it can be used in our context only for velocities or frequencies that do not exceed a certain threshold [98,99]. There is, therefore, a first important limitation to the dynamics of the system. In addition, for obtaining the Kramers formula, it is important to be able to identify the two states, although embedded into the thermal fluctuations. It means that the difference between the energy of each state and the energy barrier (i.e., the energy jump to be overcome in both directions) must be much larger than the thermal energy $K_B T$. The second limitation, thus, concerns the system temperature, which must be limited. Another limitation concerns the shape of the potential wells and the energy barrier, which must be easily approximated by parabolic curves. We also remark that different approximations exist for the prefactor in front of the Boltzmann exponential
we reduce the value of $\omega$. Indeed, if $\omega$ is quite small, the rate effects are negligible and the forward and reverse trajectories must asymptotically coincide. Recall also that the effect of velocity observed in the uniform traction case is to increase the force required to impose the Helmholtz constraint on the final position of the system. In the case of the periodic driving, we observe a similar phenomenon and the force is an increasing function of the applied frequency.

<figure><img src="image_13.png" /><figcaption>**Figure 8.** (Left panel) Helmholtz force-extension response of the system with different traction velocities. (Right panel) zoom on the transition region. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1$, $r_0 = 120$, $K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Hysteretic behavior of the Helmholtz model. (Left panel) probability $p_2$ versus the prescribed extension $Y$ (we note that $p_1 + p_2 = 1$). (Right panel) average force $\langle f \rangle$ versus extension $Y$. We used an applied extension $Y(t) = A \sin(\omega t)$ with different amplitudes $A = 2.45, 2.5, 2.55, 2.6, 2.7, 3, 3.5$ and $\omega = 2.5$ (arbitrary units). In all cases, we adopted the parameters $l = k = h = 1$, $r_0 = 120$, $K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The approach discussed is based on the Kramers formula (or similar expressions), giving the rates of the transitions between the two states in terms of the parameters of the system. There are, therefore, several limitations, as discussed below. First of all, the Kramers formula has been obtained through a steady-state assumption (slow quasi-stationary escape of particles over the barrier), and therefore, it can be used in our context only for velocities or frequencies that do not exceed a certain threshold [98,99]. There is, therefore, a first important limitation to the dynamics of the system. In addition, for obtaining the Kramers formula, it is important to be able to identify the two states, although embedded into the thermal fluctuations. It means that the difference between the energy of each state and the energy barrier (i.e., the energy jump to be overcome in both directions) must be much larger than the thermal energy $K_B T$. The second limitation, thus, concerns the system temperature, which must be limited. Another limitation concerns the shape of the potential wells and the energy barrier, which must be easily approximated by parabolic curves. We also remark that different approximations exist for the prefactor in front of the Boltzmann exponential
in the Kramers formula [98]. These different prefactors depend on the friction experienced by the system motion and one can adopt the overdamped or the underdamped assumption, or more general intermediate conditions [102]. In our context, it is quite difficult to know or measure reasonable values for the friction coefficients of the complex systems under investigation and, therefore, the prefactor is typically used as a fitting parameter adopted to give an interpretation of experimental data (by doing so, we have no substantial differences between the Arrhenius formula, transition state theory, and the different formulations of the Kramers formula). In this theory, all quantum effects are considered negligible. From the mathematical point of view, the validity of the Kramers law and similar expressions is discussed in Ref. [100]. Moreover, recent developments in Kramers rate theory can be found in Ref. [101].

We would also mention that the mechanism of stochastic resonance is closely related to the Kramers escape rate. This concept is based on the physical synchronization between the Kramers escape timescale and the periodic time of a weak modulation externally applied and acting on the system. This means that sometimes the noise-induced hopping between potential wells may become synchronized with the weak periodic forcing, generating the so-called stochastic resonance [114–116].

The origin of the modern reaction-rate theory goes back to the empirical Arrhenius law introduced to give an interpretation of chemical reaction experimental data [98,102]. It is based on a Boltzmann exponential where the energy term corresponds to the threshold energy for the chemical activation. This result was placed in the framework of nonequilibrium statistical mechanics by Kramers, who obtained a similar expression as an approximate solution of the Smoluchovski Equation [99]. The various approximations of this result are now part of the so-called transition state theory, which has applications to many natural and engineered systems, such as chemical reactions, biological processes, combustion processes, and energy conversion devices [98]. Of course, all models of bi- or multistable systems are based on these theories, and classic examples include the following. One of the first models for biopolymer extensibility was based on a two-level system, able to describe both the polysaccharide dextran and the muscle protein titin behaviors [92]. In these cases, the rate dependence of the response is perfectly described by a transition state theory implementation. A more general model has been developed for arbitrary structures having a multiwell energy landscape [70]. This approach has been applied to study the rate effects in a chain of bistable elements [71]. This paradigmatic case is useful for understanding several physical phenomena including nanoindentation, plasticity, and force-extension experiments in macromolecular materials, such as DNA, spider silk, or artificial elastomers [71]. Interesting generalizations will concern adhesion and fracture processes.

## 7. Exact Solution within the Gibbs Ensemble

We introduce the Gibbs ensemble, characterized by the presence of a prescribed force $f$ (see right panel of Figure 1), by writing the corresponding total energy:

$$ \Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + E_b(y), \quad (46) $$

By observing this function, we can notice that the addition of the potential energy contribution deriving from the applied prescribed force $f$ to the upper substrate is what differs from the Helmholtz case along with the fact that, now, $Y$ belongs to the phase space as $y$ does. The definition of $E_b(y)$ instead is equal to the one introduced in Equation (2), namely:

$$ E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}hY_M^2 & \text{if } |y| \ge Y_M. \end{cases} \quad (47) $$

In the same spirit, adopting the same characteristic function $\chi_A(y)$, see Equation (3), introduced for the Helmholtz case, we can write the total energy of the system as:
in the Kramers formula [98]. These different prefactors depend on the friction experienced by the system motion and one can adopt the overdamped or the underdamped assumption, or more general intermediate conditions [102]. In our context, it is quite difficult to know or measure reasonable values for the friction coefficients of the complex systems under investigation and, therefore, the prefactor is typically used as a fitting parameter adopted to give an interpretation of experimental data (by doing so, we have no substantial differences between the Arrhenius formula, transition state theory, and the different formulations of the Kramers formula). In this theory, all quantum effects are considered negligible. From the mathematical point of view, the validity of the Kramers law and similar expressions is discussed in Ref. [100]. Moreover, recent developments in Kramers rate theory can be found in Ref. [101].

We would also mention that the mechanism of stochastic resonance is closely related to the Kramers escape rate. This concept is based on the physical synchronization between the Kramers escape timescale and the periodic time of a weak modulation externally applied and acting on the system. This means that sometimes the noise-induced hopping between potential wells may become synchronized with the weak periodic forcing, generating the so-called stochastic resonance [114–116].

The origin of the modern reaction-rate theory goes back to the empirical Arrhenius law introduced to give an interpretation of chemical reaction experimental data [98,102]. It is based on a Boltzmann exponential where the energy term corresponds to the threshold energy for the chemical activation. This result was placed in the framework of nonequilibrium statistical mechanics by Kramers, who obtained a similar expression as an approximate solution of the Smoluchovski Equation [99]. The various approximations of this result are now part of the so-called transition state theory, which has applications to many natural and engineered systems, such as chemical reactions, biological processes, combustion processes, and energy conversion devices [98]. Of course, all models of bi- or multistable systems are based on these theories, and classic examples include the following. One of the first models for biopolymer extensibility was based on a two-level system, able to describe both the polysaccharide dextran and the muscle protein titin behaviors [92]. In these cases, the rate dependence of the response is perfectly described by a transition state theory implementation. A more general model has been developed for arbitrary structures having a multiwell energy landscape [70]. This approach has been applied to study the rate effects in a chain of bistable elements [71]. This paradigmatic case is useful for understanding several physical phenomena including nanoindentation, plasticity, and force-extension experiments in macromolecular materials, such as DNA, spider silk, or artificial elastomers [71]. Interesting generalizations will concern adhesion and fracture processes.

## 7. Exact Solution within the Gibbs Ensemble

We introduce the Gibbs ensemble, characterized by the presence of a prescribed force $f$ (see right panel of Figure 1), by writing the corresponding total energy:

$$ \Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + E_b(y), \quad (46) $$

By observing this function, we can notice that the addition of the potential energy contribution deriving from the applied prescribed force $f$ to the upper substrate is what differs from the Helmholtz case along with the fact that, now, $Y$ belongs to the phase space as $y$ does. The definition of $E_b(y)$ instead is equal to the one introduced in Equation (2), namely:

$$ E_b(y) = \begin{cases} \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ \frac{1}{2}hY_M^2 & \text{if } |y| \ge Y_M. \end{cases} \quad (47) $$

In the same spirit, adopting the same characteristic function $\chi_A(y)$, see Equation (3), introduced for the Helmholtz case, we can write the total energy of the system as:
$$ \Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2\chi_A(y) + \frac{1}{2}hY_M^2(1 - \chi_A(y)). \quad (48) $$

This potential energy can then be used to evaluate the partition function of the system:

$$ Z_G(f) = \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{-\frac{\Phi_G(y, Y)}{K_B T}} dY, \quad (49) $$

which, now, is a function of the prescribed force $f$. We briefly observe that it is possible to swap the two integrals appearing in the definition of the partition function thanks to Fubini's theorem. Using the Gaussian integral introduced earlier in Equation (7), we can perform the integrals present in $Z_G(f)$, eventually obtaining:

$$ Z_G(f) = \frac{\pi K_B T}{\sqrt{l k}} e^{-\frac{h Y_M^2}{2 K_B T} + \frac{f^2}{2 K_B T} \frac{l+k}{lk}} [\Lambda(k) + 2] \\ - \frac{\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2 K_B T} \frac{l+k+h}{l(k+h)}} \Lambda(k+h), \quad (50) $$

where, for the sake of readability, we use the function $\alpha(p)$, previously defined in Equation (10), and we introduce a new function $\Lambda(p)$, similar to the function $\Gamma(p)$ defined in Equation (11), as:

$$ \alpha(p) = \frac{p}{2K_B T}, \quad (51) $$

$$ \Lambda(p) = \phi \left[ \sqrt{\alpha(p)} \left( -Y_M - \frac{f}{p} \right) \right] - \phi \left[ \sqrt{\alpha(p)} \left( Y_M - \frac{f}{p} \right) \right]. \quad (52) $$

Thanks to the partition function, we can now evaluate the average value of the elongation of the system $\langle Y \rangle$ and the average number of intact springs $\langle \chi_A(y) \rangle$, given a fixed applied force $f$:

$$ \langle Y \rangle = K_B T \frac{\partial \log Z_G}{\partial f} = - \frac{\partial \mathcal{G}}{\partial f}, \quad (53) $$

$$ \langle \chi_A(y) \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_G}{\partial Y_M}, \quad (54) $$

where $\mathcal{G} = -K_B T \log Z_G$ is the Gibbs free energy.

In Figure 10, we show the average adimensional elongation parameter $\langle Y \rangle / Y_M$ (left panel) and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (right panel) versus the prescribed force $f$ for a system composed of springs of elastic constants $l = k = h = 1$ and for different thermal to elastic energy ration $K_B T / (hY_M^2)$. As we can see from the force-extension panel of the figure, the mechanical response to the variation of the prescribed force $f$, especially for low temperatures, presents a force plateau that corresponds to a transition for the system's rigidity. This change in rigidity occurs when the force $f$ reaches a specific value $f = f^*$ (Maxwell force, see Ref. [73]) that does not correspond to the result one could obtain basing the analysis only on mechanical principles (as seen before for the Helmholtz case). If we think about the system from a purely mechanical point of view, in fact, the force necessary to break the breakable spring $h$ should be $f = (k+h)Y_M$ that, in the case shown in Figure 10 where the elastic constants are $l = k = h = 1$ and the elongation threshold is $Y_M = 1$, is equal to $f = 2$ that is higher than the correct one shown in the figure ($f = 2 > f^*$). This apparent discrepancy will be fully understood and explained in the following discussion while adopting an energy-balance approach as previously done for the Helmholtz case.
$$ \Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2\chi_A(y) + \frac{1}{2}hY_M^2(1 - \chi_A(y)). \quad (48) $$

This potential energy can then be used to evaluate the partition function of the system:

$$ Z_G(f) = \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{-\frac{\Phi_G(y, Y)}{K_B T}} dY, \quad (49) $$

which, now, is a function of the prescribed force $f$. We briefly observe that it is possible to swap the two integrals appearing in the definition of the partition function thanks to Fubini's theorem. Using the Gaussian integral introduced earlier in Equation (7), we can perform the integrals present in $Z_G(f)$, eventually obtaining:

$$ Z_G(f) = \frac{\pi K_B T}{\sqrt{l k}} e^{-\frac{h Y_M^2}{2 K_B T} + \frac{f^2}{2 K_B T} \frac{l+k}{lk}} [\Lambda(k) + 2] \\ - \frac{\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2 K_B T} \frac{l+k+h}{l(k+h)}} \Lambda(k+h), \quad (50) $$

where, for the sake of readability, we use the function $\alpha(p)$, previously defined in Equation (10), and we introduce a new function $\Lambda(p)$, similar to the function $\Gamma(p)$ defined in Equation (11), as:

$$ \alpha(p) = \frac{p}{2K_B T}, \quad (51) $$

$$ \Lambda(p) = \phi \left[ \sqrt{\alpha(p)} \left( -Y_M - \frac{f}{p} \right) \right] - \phi \left[ \sqrt{\alpha(p)} \left( Y_M - \frac{f}{p} \right) \right]. \quad (52) $$

Thanks to the partition function, we can now evaluate the average value of the elongation of the system $\langle Y \rangle$ and the average number of intact springs $\langle \chi_A(y) \rangle$, given a fixed applied force $f$:

$$ \langle Y \rangle = K_B T \frac{\partial \log Z_G}{\partial f} = - \frac{\partial \mathcal{G}}{\partial f}, \quad (53) $$

$$ \langle \chi_A(y) \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_G}{\partial Y_M}, \quad (54) $$

where $\mathcal{G} = -K_B T \log Z_G$ is the Gibbs free energy.

In Figure 10, we show the average adimensional elongation parameter $\langle Y \rangle / Y_M$ (left panel) and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ (right panel) versus the prescribed force $f$ for a system composed of springs of elastic constants $l = k = h = 1$ and for different thermal to elastic energy ration $K_B T / (hY_M^2)$. As we can see from the force-extension panel of the figure, the mechanical response to the variation of the prescribed force $f$, especially for low temperatures, presents a force plateau that corresponds to a transition for the system's rigidity. This change in rigidity occurs when the force $f$ reaches a specific value $f = f^*$ (Maxwell force, see Ref. [73]) that does not correspond to the result one could obtain basing the analysis only on mechanical principles (as seen before for the Helmholtz case). If we think about the system from a purely mechanical point of view, in fact, the force necessary to break the breakable spring $h$ should be $f = (k+h)Y_M$ that, in the case shown in Figure 10 where the elastic constants are $l = k = h = 1$ and the elongation threshold is $Y_M = 1$, is equal to $f = 2$ that is higher than the correct one shown in the figure ($f = 2 > f^*$). This apparent discrepancy will be fully understood and explained in the following discussion while adopting an energy-balance approach as previously done for the Helmholtz case.
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Behavior of the two-rigidity model with variable thermal to elastic energy ratio $K_B T/(hY_M^2) = \{0, 0.01, 0.02, 0.03\}$ (purple, blue, yellow, and red curves, respectively). The purple curves correspond to the pure mechanical case at zero temperature. The average elongation $\langle Y \rangle$ and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ are presented versus the fixed force $f$, where, in both cases, $l = k = h = 1$ and $Y_M = 1$. In both panels, the force responsible for the change in rigidity, $f^*$, is shown.</figcaption></figure>

The difference between the two rigidity regimes emerges from the difference in slope of the two curve branches before and after the breaking force point and the same considerations done in the Helmholtz case hold in this case as well. The threshold force $f^*$, responsible for the transition between the two rigidity, is even more clear when observing the right panel of Figure 10 where the number of intact breakable spring $h$ changes from 1 to 0 exactly at $f = f^*$. Furthermore, we stress the fact that, even in the Gibbs ensemble, $f^*$ does not depend on the temperature that, instead, is responsible only for smoothing the curves while increasing. This temperature independence of the threshold force $f^*$ is expected since the model is made of a single unit.

## 8. Gibbs Ensemble with Spin Approximation

As introduced previously in the Helmholtz case, in this section, we adopt the spin variable approach to study the system in the Gibbs ensemble. As in previous studies, we introduce a spin variable $S$ to help us approximate the total potential energy of the system. The spin variable will then assume the value $S = 0$ when the breakable spring is in the broken state, and $S = 1$ when it is in the intact conformation. By doing so, we are able to express the total energy as:

$$ \Phi_G(y, Y, S) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2S + \frac{1}{2}hY_M^2(1 - S), \quad (55) $$

which, now, is a function of $y$, $Y$ and $S$. We proceed by evaluating the partition function of the system that, in the spin variable approximation, assumes the following definition:

$$ Z_G(f) = \sum_{S=0,1} \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{-\frac{\Phi_G(y, Y, S)}{k_B T}} dY. \quad (56) $$

By introducing here the expression for the total energy and explicitly expressing the sum over the spin values, we obtain:

$$ Z_G(f) = \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{\frac{fY}{k_B T} - \frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dY \\ + \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{\frac{fY}{k_B T} - \frac{l}{2k_B T}(Y-y)^2 - \frac{k}{2k_B T}y^2 - \frac{h}{2k_B T}Y_M^2} dY, \quad (57) $$

which is easily solved using the Gaussian integral, eventually obtaining:
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Behavior of the two-rigidity model with variable thermal to elastic energy ratio $K_B T/(hY_M^2) = \{0, 0.01, 0.02, 0.03\}$ (purple, blue, yellow, and red curves, respectively). The purple curves correspond to the pure mechanical case at zero temperature. The average elongation $\langle Y \rangle$ and the average number of intact breakable springs $\langle \chi_A(y) \rangle$ are presented versus the fixed force $f$, where, in both cases, $l = k = h = 1$ and $Y_M = 1$. In both panels, the force responsible for the change in rigidity, $f^*$, is shown.</figcaption></figure>

The difference between the two rigidity regimes emerges from the difference in slope of the two curve branches before and after the breaking force point and the same considerations done in the Helmholtz case hold in this case as well. The threshold force $f^*$, responsible for the transition between the two rigidity, is even more clear when observing the right panel of Figure 10 where the number of intact breakable spring $h$ changes from 1 to 0 exactly at $f = f^*$. Furthermore, we stress the fact that, even in the Gibbs ensemble, $f^*$ does not depend on the temperature that, instead, is responsible only for smoothing the curves while increasing. This temperature independence of the threshold force $f^*$ is expected since the model is made of a single unit.

## 8. Gibbs Ensemble with Spin Approximation

As introduced previously in the Helmholtz case, in this section, we adopt the spin variable approach to study the system in the Gibbs ensemble. As in previous studies, we introduce a spin variable $S$ to help us approximate the total potential energy of the system. The spin variable will then assume the value $S = 0$ when the breakable spring is in the broken state, and $S = 1$ when it is in the intact conformation. By doing so, we are able to express the total energy as:

$$ \Phi_G(y, Y, S) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2S + \frac{1}{2}hY_M^2(1 - S), \quad (55) $$

which, now, is a function of $y$, $Y$ and $S$. We proceed by evaluating the partition function of the system that, in the spin variable approximation, assumes the following definition:

$$ Z_G(f) = \sum_{S=0,1} \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{-\frac{\Phi_G(y, Y, S)}{K_B T}} dY. \quad (56) $$

By introducing here the expression for the total energy and explicitly expressing the sum over the spin values, we obtain:

$$ Z_G(f) = \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{\frac{fY}{K_B T} - \frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dY \\ + \int_{\mathbb{R}} dy \int_{\mathbb{R}} e^{\frac{fY}{K_B T} - \frac{l}{2K_B T}(Y-y)^2 - \frac{k}{2K_B T}y^2 - \frac{h}{2K_B T}Y_M^2} dY, \quad (57) $$

which is easily solved using the Gaussian integral, eventually obtaining:
$$Z_G(f) = \frac{2\pi K_B T}{\sqrt{lk}} e^{\frac{f^2}{2K_B T} \left(\frac{l+k}{lk}\right)} - \frac{hY_M^2}{2K_B T} + \frac{2\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2K_B T} \left(\frac{l+k+h}{l(k+h)}\right)}. \quad (58)$$

Thanks to this approximated partition function, we can evaluate the following quantities:

$$\langle Y \rangle = K_B T \frac{\partial \log Z_G}{\partial f}, \quad (59)$$

$$\langle S \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_G}{\partial Y_M}. \quad (60)$$

In Figure 11, we can see that the approximations obtained both for $\langle Y \rangle / Y_M$ and $\langle S \rangle$ versus the fixed force $f$ are in good agreement with the results obtained with the exact calculations. As in the Helmholtz spin approximation, once again the spin approximation is more accurate for low temperatures as we should expect. Furthermore, we observe that even in this case, the rupture force does not correspond to the one we would expect but, instead, is shifted as in the exact case. In the following section, we will solve this apparent dilemma by adopting an energy-balance approach.

<figure><img src="image_16.png" /><figcaption>**Figure 11.** Comparison between the spin-approximated quantities (dashed curves) and the corresponding exact results (continuous curves) for $\langle Y \rangle / Y_M$ and $\langle S \rangle$ versus the applied force $f$. In both cases, the elastic constants are $l = k = h = 1$ and $Y_M = 1$, while the thermal to elastic energy ratio assumes the values $K_B T / (hY_M^2) = \{0.01, 0.2\}$ (respectively, in blue and in red).</figcaption></figure>

## 9. Gibbs Ensemble from an Energetic Perspective

In this section, we study the model in the isotensional configuration adopting the principles of energy minimization in order to obtain the exact rupture force $f^*$ responsible for the rigidity change of the system. The total potential energy of the system is given by:

$$\Phi_G(y, Y) = \begin{cases} -fY + \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ -fY + \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2 & \text{if } |y| > Y_M. \end{cases} \quad (61)$$

Given a fixed force $f$, we search for the quantities $y$ and $Y$ that minimize the potential energy. We start by considering our system in the $|y| < Y_M$ case, where the potential energy is defined as:

$$\Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2. \quad (62)$$

We derive $\Phi_G(y, Y)$ by $y$ and $Y$ separately, obtaining:

$$\frac{\partial \Phi_G(y, Y)}{\partial y} = -l(Y - y) + hy + ky = 0, \quad (63)$$

$$\frac{\partial \Phi_G(y, Y)}{\partial Y} = -f + l(Y - y) = 0. \quad (64)$$
$$Z_G(f) = \frac{2\pi K_B T}{\sqrt{l k}} e^{\frac{f^2}{2K_B T} \left(\frac{l+k}{l k}\right)} e^{-\frac{h Y_M^2}{2K_B T}} + \frac{2\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2K_B T} \left(\frac{l+k+h}{l(k+h)}\right)}. \quad (58)$$

Thanks to this approximated partition function, we can evaluate the following quantities:

$$\langle Y \rangle = K_B T \frac{\partial \log Z_G}{\partial f}, \quad (59)$$

$$\langle S \rangle = 1 + \frac{K_B T}{h Y_M} \frac{\partial \log Z_G}{\partial Y_M}. \quad (60)$$

In Figure 11, we can see that the approximations obtained both for $\langle Y \rangle / Y_M$ and $\langle S \rangle$ versus the fixed force $f$ are in good agreement with the results obtained with the exact calculations. As in the Helmholtz spin approximation, once again the spin approximation is more accurate for low temperatures as we should expect. Furthermore, we observe that even in this case, the rupture force does not correspond to the one we would expect but, instead, is shifted as in the exact case. In the following section, we will solve this apparent dilemma by adopting an energy-balance approach.

<figure><img src="image_16.png" /><figcaption>Figure 11. Comparison between the spin-approximated quantities (dashed curves) and the corresponding exact results (continuous curves) for $\langle Y \rangle / Y_M$ and $\langle S \rangle$ versus the applied force $f$. In both cases, the elastic constants are $l = k = h = 1$ and $Y_M = 1$, while the thermal to elastic energy ratio assumes the values $K_B T / (h Y_M^2) = \{0.01, 0.2\}$ (respectively, in blue and in red).</figcaption></figure>

## 9. Gibbs Ensemble from an Energetic Perspective

In this section, we study the model in the isotensional configuration adopting the principles of energy minimization in order to obtain the exact rupture force $f^*$ responsible for the rigidity change of the system. The total potential energy of the system is given by:

$$\Phi_G(y, Y) = \begin{cases} -fY + \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2 & \text{if } |y| < Y_M, \\ -fY + \frac{1}{2}l(Y-y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hY_M^2 & \text{if } |y| > Y_M. \end{cases} \quad (61)$$

Given a fixed force $f$, we search for the quantities $y$ and $Y$ that minimize the potential energy. We start by considering our system in the $|y| < Y_M$ case, where the potential energy is defined as:

$$\Phi_G(y, Y) = -fY + \frac{1}{2}l(Y - y)^2 + \frac{1}{2}ky^2 + \frac{1}{2}hy^2. \quad (62)$$

We derive $\Phi_G(y, Y)$ by $y$ and $Y$ separately, obtaining:

$$\frac{\partial \Phi_G(y, Y)}{\partial y} = -l(Y - y) + hy + ky = 0, \quad (63)$$

$$\frac{\partial \Phi_G(y, Y)}{\partial Y} = -f + l(Y - y) = 0. \quad (64)$$
After some straightforward calculations, we obtain the two minimizing values:

$$y_{\min} = \frac{f}{h+k'} \qquad (65)$$

$$Y_{\min} = \frac{f}{h+k} + \frac{f}{l} \qquad (66)$$

that we use to evaluate the resulting minimal potential energy, obtaining:

$$\Phi_G(y_{\min}, Y_{\min}) = -\frac{1}{2} \frac{f^2}{l} - \frac{1}{2} \frac{f^2}{h+k}. \qquad (67)$$

Repeating the same process for the $|y| > Y_M$ case, we obtain the following minimized potential energy:

$$\Phi_G(y_{\min}, Y_{\min}) = -\frac{1}{2} \frac{f^2}{k} - \frac{1}{2} \frac{f^2}{l} + \frac{1}{2} h Y_M^2, \qquad (68)$$

and, comparing the two minimized potential energy, we obtain the exact rupture force at which the system changes its rigidity, namely:

$$f^* = \sqrt{k(h+k)}Y_M \qquad (69)$$

We observe that using the same parameters of Figures 10 and 11, we have $f^* = \sqrt{2}$ in agreement with the results shown in the figures. Therefore, as in the Helmholtz ensemble, also in the present Gibbs case, the energetic approach (without thermal effects) is in perfect agreement with the exact calculation based on statistical mechanics and with the spin approximation. This confirms the value of the rupture threshold given in Equation (69) and underlines the fact that the correct approach to determining the rupture thresholds is the one based on energetic assumptions. These consideration are also useful to introduce the out-of-equilibrium behavior of the system, as discussed in the next section.

## 10. Rate Effect on the Rupture Model under Gibbs Conditions

We introduce here the theoretical development necessary to describe the rate effects in the rupture Gibbs model. The energy profiles for this model have been obtained in Equation (61) for both $|y| < Y_M$ (intact spring) and $|y| > Y_M$ (broken spring). Unlike the Helmholtz case, in this Gibbs ensemble, we have that the energy profiles depend on two variables, namely $y$ and $Y$. To apply transition state theory, i.e., Kramers formula, we must minimize these energies with respect to the variable $Y$ that does not participate directly in the transition. We see from Equation (64) that the relation $\frac{\partial\Phi_G(y,Y)}{\partial Y} = -f + l(Y - y) = 0$ is valid for both the regions $|y| < Y_M$ (intact spring) and $|y| > Y_M$ (broken spring). Therefore, we can eliminate the variable $Y$ by means of the relation $Y = y + f/l$ in both energy profiles. This allows us to graphically represent the two profiles as shown in Figure 12. In particular, we can see how these profiles change with a varying applied force $f$. In the previous section, we obtained the value of the transition force $f^*$, which corresponds to $\sqrt{2}$ for the parameters adopted in the figures (arbitrary units). It can be seen that for values of the force far from this threshold, the two potential minima are noticeably different from each other, and thus, the state of the system is well identified; see the left panel of Figure 12. When, on the other hand, the force values are close to the transition threshold, an energy barrier $B$ is identified which must be crossed to effect the transition; see the right panel of Figure 12. The value of the barrier $B$ is calculated as:

$$B(f) = \Phi_G(Y_M, Y_M + \frac{f}{l}) = -f(Y_M + \frac{f}{l}) + \frac{1}{2}l\left(\frac{f}{l}\right)^2 + \frac{1}{2}kY_M^2 + \frac{1}{2}hY_M^2, \qquad (70)$$
After some straightforward calculations, we obtain the two minimizing values:

$$y_{\min} = \frac{f}{h+k'} \qquad (65)$$

$$Y_{\min} = \frac{f}{h+k} + \frac{f}{l} \qquad (66)$$

that we use to evaluate the resulting minimal potential energy, obtaining:

$$\Phi_G(y_{\min}, Y_{\min}) = -\frac{1}{2} \frac{f^2}{l} - \frac{1}{2} \frac{f^2}{h+k}. \qquad (67)$$

Repeating the same process for the $|y| > Y_M$ case, we obtain the following minimized potential energy:

$$\Phi_G(y_{\min}, Y_{\min}) = -\frac{1}{2} \frac{f^2}{k} - \frac{1}{2} \frac{f^2}{l} + \frac{1}{2} h Y_M^2, \qquad (68)$$

and, comparing the two minimized potential energy, we obtain the exact rupture force at which the system changes its rigidity, namely:

$$f^* = \sqrt{k(h+k)}Y_M \qquad (69)$$

We observe that using the same parameters of Figures 10 and 11, we have $f^* = \sqrt{2}$ in agreement with the results shown in the figures. Therefore, as in the Helmholtz ensemble, also in the present Gibbs case, the energetic approach (without thermal effects) is in perfect agreement with the exact calculation based on statistical mechanics and with the spin approximation. This confirms the value of the rupture threshold given in Equation (69) and underlines the fact that the correct approach to determining the rupture thresholds is the one based on energetic assumptions. These consideration are also useful to introduce the out-of-equilibrium behavior of the system, as discussed in the next section.

## 10. Rate Effect on the Rupture Model under Gibbs Conditions

We introduce here the theoretical development necessary to describe the rate effects in the rupture Gibbs model. The energy profiles for this model have been obtained in Equation (61) for both $|y| < Y_M$ (intact spring) and $|y| > Y_M$ (broken spring). Unlike the Helmholtz case, in this Gibbs ensemble, we have that the energy profiles depend on two variables, namely $y$ and $Y$. To apply transition state theory, i.e., Kramers formula, we must minimize these energies with respect to the variable $Y$ that does not participate directly in the transition. We see from Equation (64) that the relation $\frac{\partial\Phi_G(y,Y)}{\partial Y} = -f + l(Y - y) = 0$ is valid for both the regions $|y| < Y_M$ (intact spring) and $|y| > Y_M$ (broken spring). Therefore, we can eliminate the variable $Y$ by means of the relation $Y = y + f/l$ in both energy profiles. This allows us to graphically represent the two profiles as shown in Figure 12. In particular, we can see how these profiles change with a varying applied force $f$. In the previous section, we obtained the value of the transition force $f^*$, which corresponds to $\sqrt{2}$ for the parameters adopted in the figures (arbitrary units). It can be seen that for values of the force far from this threshold, the two potential minima are noticeably different from each other, and thus, the state of the system is well identified; see the left panel of Figure 12. When, on the other hand, the force values are close to the transition threshold, an energy barrier $B$ is identified which must be crossed to effect the transition; see the right panel of Figure 12. The value of the barrier $B$ is calculated as:

$$B(f) = \Phi_G(Y_M, Y_M + \frac{f}{l}) = -f(Y_M + \frac{f}{l}) + \frac{1}{2}l\left(\frac{f}{l}\right)^2 + \frac{1}{2}kY_M^2 + \frac{1}{2}hY_M^2, \qquad (70)$$
and we see that it is dependent on the applied force $f$. We also take into consideration the values of the two energy minima, which have been calculated in Equations (67) and (68). For $|y| < Y_M$ (intact spring) we have:

$$ \Phi_1(f) = -\frac{1}{2} \frac{f^2}{l} - \frac{1}{2} \frac{f^2}{h+k}. \qquad (71) $$

and for $|y| > Y_M$ (broken spring), we get:

$$ \Phi_2(2) = -\frac{1}{2} \frac{f^2}{k} - \frac{1}{2} \frac{f^2}{l} + \frac{1}{2} h Y_M^2. \qquad (72) $$

<figure><img src="image_17.png" /><figcaption>Figure 12. (Left panel) Gibbs energy profiles with different values of $f$ and $y$. Red lines correspond to the energies $\Phi_G$ for $|y| < Y_M$, while blue lines correspond to $\Phi_G$ for $|y| > Y_M$. (Right panel) zoom with three values of $Y$ closer to $Y = \sqrt{2}$, useful to identify the barriers $B$. In all cases, the elastic constants are $l = k = h = 1$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

As before, given the three parameters $B(f)$, $\Phi_1(f)$, and $\Phi_2(f)$, we can define the two following kinetic coefficients describing the dynamics of the Gibbs system:

$$ k_{12} = r_0 e^{-\frac{B(Y)-\Phi_1(Y)}{K_B T}}, \qquad (73) $$

$$ k_{21} = r_0 e^{-\frac{B(Y)-\Phi_2(Y)}{K_B T}}, \qquad (74) $$

where $r_0$ is a prefactor measured in $s^{-1}$, which can be obtained through the original Kramers formula or one of its generalizations [98,99]. As already discussed, it is considered as an additional parameter of the system.

We can now define the probability $p_1$ to be in the intact state and the probability $p_2$ to be in the broken state of the system. These two probabilities evolve in time by means of the following system of differential Equations [61]:

$$ \frac{d p_1(t)}{dt} = -k_{12} p_1(t) + k_{21} p_2(t), \qquad (75) $$

$$ \frac{d p_2(t)}{dt} = -k_{21} p_2(t) + k_{12} p_1(t), \qquad (76) $$

where the trajectory $f(t)$, characterizing $B(f)$, $\Phi_1(f)$, and $\Phi_2(f)$, must be imposed to obtain the system dynamics. The Gibbs partition function obtained through the spin approach, see Equation (58), can be written as $Z_G = Z_{G,1} + Z_{G,2}$, where:

$$ Z_{G,1}(f) = \frac{2\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2K_B T} \left( \frac{l+k+h}{l(k+h)} \right)}, \qquad (77) $$

$$ Z_{G,2}(f) = \frac{2\pi K_B T}{\sqrt{l k}} e^{\frac{f^2}{2K_B T} \left( \frac{l+k}{l} \right) - \frac{h Y_M^2}{2K_B T}}. \qquad (78) $$
and we see that it is dependent on the applied force $f$. We also take into consideration the values of the two energy minima, which have been calculated in Equations (67) and (68). For $|y| < Y_M$ (intact spring) we have:

$$ \Phi_1(f) = -\frac{1}{2} \frac{f^2}{l} - \frac{1}{2} \frac{f^2}{h+k}. \qquad (71) $$

and for $|y| > Y_M$ (broken spring), we get:

$$ \Phi_2(2) = -\frac{1}{2} \frac{f^2}{k} - \frac{1}{2} \frac{f^2}{l} + \frac{1}{2} h Y_M^2. \qquad (72) $$

<figure><img src="image_17.png" /><figcaption>Figure 12. (Left panel) Gibbs energy profiles with different values of $f$ and $y$. Red lines correspond to the energies $\Phi_G$ for $|y| < Y_M$, while blue lines correspond to $\Phi_G$ for $|y| > Y_M$. (Right panel) zoom with three values of $Y$ closer to $Y = \sqrt{2}$, useful to identify the barriers $B$. In all cases, the elastic constants are $l = k = h = 1$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

As before, given the three parameters $B(f)$, $\Phi_1(f)$, and $\Phi_2(f)$, we can define the two following kinetic coefficients describing the dynamics of the Gibbs system:

$$ k_{12} = r_0 e^{-\frac{B(Y)-\Phi_1(Y)}{K_B T}}, \qquad (73) $$

$$ k_{21} = r_0 e^{-\frac{B(Y)-\Phi_2(Y)}{K_B T}}, \qquad (74) $$

where $r_0$ is a prefactor measured in $s^{-1}$, which can be obtained through the original Kramers formula or one of its generalizations [98,99]. As already discussed, it is considered as an additional parameter of the system.

We can now define the probability $p_1$ to be in the intact state and the probability $p_2$ to be in the broken state of the system. These two probabilities evolve in time by means of the following system of differential Equations [61]:

$$ \frac{d p_1(t)}{dt} = -k_{12} p_1(t) + k_{21} p_2(t), \qquad (75) $$

$$ \frac{d p_2(t)}{dt} = -k_{21} p_2(t) + k_{12} p_1(t), \qquad (76) $$

where the trajectory $f(t)$, characterizing $B(f)$, $\Phi_1(f)$, and $\Phi_2(f)$, must be imposed to obtain the system dynamics. The Gibbs partition function obtained through the spin approach, see Equation (58), can be written as $Z_G = Z_{G,1} + Z_{G,2}$, where:

$$ Z_{G,1}(f) = \frac{2\pi K_B T}{\sqrt{l(k+h)}} e^{\frac{f^2}{2K_B T} \left( \frac{l+k+h}{l(k+h)} \right)}, \qquad (77) $$

$$ Z_{G,2}(f) = \frac{2\pi K_B T}{\sqrt{l k}} e^{\frac{f^2}{2K_B T} \left( \frac{l+k}{l} \right) - \frac{h Y_M^2}{2 K_B T}}. \qquad (78) $$
from which we obtain the force-extension response pertinent to each state:

$$ \langle Y \rangle_1 = K_B T \frac{\partial \log Z_{G,1}(f)}{\partial f} = \frac{l+k+h}{l(k+h)} f, \qquad (79) $$

$$ \langle Y \rangle_2 = K_B T \frac{\partial \log Z_{G,2}(f)}{\partial f} = \frac{l+k}{lk} f. \qquad (80) $$

We can, therefore, determine the average value of the extension $\langle Y \rangle$ as function of the applied force $f$. The resulting force-extension curve is given by:

$$ \langle Y \rangle(t) = p_1(t) \frac{l+k+h}{l(k+h)} f(t) + p_2(t) \frac{l+k}{lk} f(t), \qquad (81) $$

where the probabilities $p_1(t)$ and $p_2(t)$ are obtained by solving the system of differential equation stated in Equations (75) and (76), for a given applied force $v(t)$. An example of solution for $p_1$ and $p_2$ is shown in Figure 13, where we considered a time evolution for the applied force given by $f = vt$, where $v$ is the force rate (with units N/s). As before, we remark that different force rates generate different dynamical regimes for the system response. This effect on the force-extension behavior can be found in Figure 14. Importantly, we note that the force for producing a transition (in this case, the breaking of the spring) is always an increasing function of the traction velocity.

<figure><img src="image_18.png" /><figcaption>Figure 13. (Left panel) Gibbs probability $p_1$ versus $Y$. (Right panel) Gibbs probability $p_2$ versus $Y$. We used different force rates $v = 0.03, 0.2, 0.5, 1.2, 2.5, 5$ (arbitrary units), defining the trajectory $f = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14. (Left panel) force-extension response of the system with different traction velocities. (Right panel) zoom on the transition region. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

As before, the dynamical model can also be used to obtain the hysteretic behavior of the Gibbs system. An example is shown in Figure 15, where we can find the response of
from which we obtain the force-extension response pertinent to each state:

$$ \langle Y \rangle_1 = K_B T \frac{\partial \log Z_{G,1}(f)}{\partial f} = \frac{l+k+h}{l(k+h)} f, \qquad (79) $$

$$ \langle Y \rangle_2 = K_B T \frac{\partial \log Z_{G,2}(f)}{\partial f} = \frac{l+k}{lk} f. \qquad (80) $$

We can, therefore, determine the average value of the extension $\langle Y \rangle$ as function of the applied force $f$. The resulting force-extension curve is given by:

$$ \langle Y \rangle(t) = p_1(t) \frac{l+k+h}{l(k+h)} f(t) + p_2(t) \frac{l+k}{lk} f(t), \qquad (81) $$

where the probabilities $p_1(t)$ and $p_2(t)$ are obtained by solving the system of differential equation stated in Equations (75) and (76), for a given applied force $v(t)$. An example of solution for $p_1$ and $p_2$ is shown in Figure 13, where we considered a time evolution for the applied force given by $f = vt$, where $v$ is the force rate (with units N/s). As before, we remark that different force rates generate different dynamical regimes for the system response. This effect on the force-extension behavior can be found in Figure 14. Importantly, we note that the force for producing a transition (in this case, the breaking of the spring) is always an increasing function of the traction velocity.

<figure><img src="image_18.png" /><figcaption>Figure 13. (Left panel) Gibbs probability $p_1$ versus $Y$. (Right panel) Gibbs probability $p_2$ versus $Y$. We used different force rates $v = 0.03, 0.2, 0.5, 1.2, 2.5, 5$ (arbitrary units), defining the trajectory $f = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 14. (Left panel) force-extension response of the system with different traction velocities. (Right panel) zoom on the transition region. We used different velocities $v = 0.5, 1.2, 2.5, 5, 10$ (arbitrary units), defining the trajectory $Y = vt$. In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

As before, the dynamical model can also be used to obtain the hysteretic behavior of the Gibbs system. An example is shown in Figure 15, where we can find the response of
the system to a sinusoidal applied force of the type $f(t) = A \sin(\omega t)$. We can observe both the evolution of the probability and of the force-extension curve. In both cases, we remark a hysteretic behavior generated by the memory effect described by the rate constants and rate equations. We plotted different curves corresponding to different values of the force amplitude. As in the Helmholtz case, we see that for small amplitudes the transition is only just started and not finished, while with larger amplitudes, the transition can be completed. In any case, the forward and reverse paths are always different because of the memory effect. The difference, as before, depend on the velocity of variation of the force, which is controlled by $\omega$.

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Hysteretic behavior of the Gibbs model. (Left panel) probability $p_1$ versus the applied force $f$ (we note that $p_1 + p_2 = 1$). (Right panel) force $f$ versus average extension $\langle Y \rangle$. We used an applied force $f(t) = A \sin(\omega t)$ with different amplitudes $A = 1.55, 1.6, 1.65, 1.7, 1.8, 2, 2.4$ and $\omega = 2.5$ (arbitrary units). In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The limitations in velocity, frequency, and temperature already discussed in Section 6 apply, of course, to both the Gibbs ensemble and the Helmholtz ensemble.

**11. Discussion and Conclusions**

In this paper, we reviewed the various approaches that can be taken into account to study the evolution of systems characterized by multistable energy. In particular, these approaches have been applied to a simple model exhibiting a breaking phenomenon and are developed for the two symmetric and dual Gibbs and Helmholtz statistical ensembles. The result of all these approaches consists of the force-extension relation characterizing the system response, which is typically depending on the combination between temperature and pulling speed (or driving frequency). The symmetry between the two statistical ensembles can be deduced through the physical observables we use to describe this force-extension response. While in the Helmholtz ensemble we have a deterministic prescribed extension $Y$ and a stochastic force described by its average value $\langle f \rangle$, in the Gibbs ensemble, we have a deterministic applied force and a stochastic extension characterized by the average value $\langle Y \rangle$ (e.g., see Figures 3 and 10, left panels). In both cases, at thermodynamic equilibrium, the determination of the force-extension relation is based on the calculation of the partition function and the corresponding free energy (see Equations (9) and (50)). The two partition functions corresponding to the Helmholtz and Gibbs ensembles are related through a Laplace transform, which is always an exact relationship valid for any physical system. For a large system (i.e., large number of elementary units), this property can be used to prove that the free energies are related through a Legendre transform [106,107]. However, in the example discussed in this paper, only the Laplace transform can be used, since the system consists of only one unit and the thermodynamic limit cannot be applied. For arbitrary systems, when the thermodynamic limit can be defined, the concept of equivalence or nonequivalence of statistical sets is introduced. Two sets are said to be equivalent when, in the thermodynamic limit, the two force-extension relations for Helmholtz and Gibbs are asymptotically coincident. In general, for a given system, it is difficult to say a priori
the system to a sinusoidal applied force of the type $f(t) = A \sin(\omega t)$. We can observe both the evolution of the probability and of the force-extension curve. In both cases, we remark a hysteretic behavior generated by the memory effect described by the rate constants and rate equations. We plotted different curves corresponding to different values of the force amplitude. As in the Helmholtz case, we see that for small amplitudes the transition is only just started and not finished, while with larger amplitudes, the transition can be completed. In any case, the forward and reverse paths are always different because of the memory effect. The difference, as before, depend on the velocity of variation of the force, which is controlled by $\omega$.

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Hysteretic behavior of the Gibbs model. (Left panel) probability $p_1$ versus the applied force $f$ (we note that $p_1 + p_2 = 1$). (Right panel) force $f$ versus average extension $\langle Y \rangle$. We used an applied force $f(t) = A \sin(\omega t)$ with different amplitudes $A = 1.55, 1.6, 1.65, 1.7, 1.8, 2, 2.4$ and $\omega = 2.5$ (arbitrary units). In all cases, we adopted the parameters $l = k = h = 1, r_0 = 120, K_B T = 0.01$, and $Y_M = 1$ (arbitrary units).</figcaption></figure>

The limitations in velocity, frequency, and temperature already discussed in Section 6 apply, of course, to both the Gibbs ensemble and the Helmholtz ensemble.

## 11. Discussion and Conclusions

In this paper, we reviewed the various approaches that can be taken into account to study the evolution of systems characterized by multistable energy. In particular, these approaches have been applied to a simple model exhibiting a breaking phenomenon and are developed for the two symmetric and dual Gibbs and Helmholtz statistical ensembles. The result of all these approaches consists of the force-extension relation characterizing the system response, which is typically depending on the combination between temperature and pulling speed (or driving frequency). The symmetry between the two statistical ensembles can be deduced through the physical observables we use to describe this force-extension response. While in the Helmholtz ensemble we have a deterministic prescribed extension $Y$ and a stochastic force described by its average value $\langle f \rangle$, in the Gibbs ensemble, we have a deterministic applied force and a stochastic extension characterized by the average value $\langle Y \rangle$ (e.g., see Figures 3 and 10, left panels). In both cases, at thermodynamic equilibrium, the determination of the force-extension relation is based on the calculation of the partition function and the corresponding free energy (see Equations (9) and (50)). The two partition functions corresponding to the Helmholtz and Gibbs ensembles are related through a Laplace transform, which is always an exact relationship valid for any physical system. For a large system (i.e., large number of elementary units), this property can be used to prove that the free energies are related through a Legendre transform [106,107]. However, in the example discussed in this paper, only the Laplace transform can be used, since the system consists of only one unit and the thermodynamic limit cannot be applied. For arbitrary systems, when the thermodynamic limit can be defined, the concept of equivalence or nonequivalence of statistical sets is introduced. Two sets are said to be equivalent when, in the thermodynamic limit, the two force-extension relations for Helmholtz and Gibbs are asymptotically coincident. In general, for a given system, it is difficult to say a priori
whether equivalence is met or not. Nevertheless, there are rather important classes of
systems for which the character of equivalence is known [110]. For example, it is known
that polymer chains with continuous two-body interactions satisfy the equivalence between
the two statistical ensembles considered (without confinement effects) [106–108]. This can
be directly applied to the freely jointed chain model and the worm-like chain model as
well. Other problems exhibit an interesting nonequivalence between the defined statis-
tical ensembles. In particular, a polymer chain compressed between two pistons shows
nonequivalence of the ensembles and a phase transition corresponding to the escape from
the gap between the pistons [109]. The desorption of a single chain from a substrate without
excluded volume interactions and other confined systems also exhibit nonequivalence of
the ensembles and the emergence of a phase transition [111–113]. Adhesion and fracture
processes show nonequivalence as well [76,77,79]. When discussing duality between statis-
tical ensembles, therefore, one must always remember the possibility of having equivalence
or nonequivalence between them, in the thermodynamic limit.

Concerning the equilibrium statistical mechanics, the most accurate method is based
on the evaluation of the partition function in closed form, without approximations. This
can be performed only for a simple system, such as the one considered in this work.
When multistable energies are considered, this technique can be applied analytically only
in some particular cases, but it can be always implemented numerically with different
computational approaches.

If we want to adopt analytical techniques even in the presence of multistable energies, we must use the method of spin variables, which adds a set of discrete variables to the phase space but simplifies the calculation of integrals in partition functions. Indeed, in this case, the partition function is defined by summing over all the spin variable and integrating over all the continuous coordinates (see Equations (22) and (58)). The integral are simpler since they are calculated over the whole domain for quadratic energy well (Gaussian integration). The results are quite accurate if excessive temperatures are not considered (see Figures 5 and 11 for the force-extension response and the spin variable behavior in both ensembles).

Often, to better understand the transition strategies of multistable systems, it is advan-
tageous to perform a zero-temperature energy study. This makes it easier to identify the
thresholds of applied forces or extensions responsible for the system’s internal transitions.
This procedure was also applied in the examples studied in this work and enabled us to
better understand the transition thresholds and especially to show that they often do not
take on the values that one might imagine with simple mechanical observations of the
system. This point in the analysis performed is perfectly dual between the two Helmholtz
and Gibbs ensembles.

A last approach concerns the out-of-equilibrium regime for our system. To deal with
this problem, we can mention two alternative approach. The first one, the most rigorous,
consists of determining the Hamiltonian function for the system and the corresponding
equation of motion. To introduce the out-of-equilibrium statistical mechanics, we can add
to this equation suitable noise and viscous friction terms, which are able to mimics the
dynamics of the system by considering temperature and rate effect. The resulting stochastic
differential equation is the so-called Langevin Equation [102]. The probability density
describing the evolution of the system state is associated with the Langevin equation
and it is referred to as the Fokker-Planck Equation [102]. For complicated systems, with
multistable energies, both the Langevin and Fokker-Planck equations can be solved only
through quite costly computational techniques. The second approach, used to simplify the
numerical solution of the problem, consists of implementing the transition state theory [98].
This method is based on the definition of the probability to be in a given state and on
the rate equations describing these probabilities. The rate constants adopted in these
equations are obtained through the Kramers formula based on the energy barrier between
the states [99]. A central point of this analysis is therefore to identify the energy barriers
in the system (see Figures 6 and 12, corresponding two the two statistical ensembles).
whether equivalence is met or not. Nevertheless, there are rather important classes of systems for which the character of equivalence is known [110]. For example, it is known that polymer chains with continuous two-body interactions satisfy the equivalence between the two statistical ensembles considered (without confinement effects) [106–108]. This can be directly applied to the freely jointed chain model and the worm-like chain model as well. Other problems exhibit an interesting nonequivalence between the defined statistical ensembles. In particular, a polymer chain compressed between two pistons shows nonequivalence of the ensembles and a phase transition corresponding to the escape from the gap between the pistons [109]. The desorption of a single chain from a substrate without excluded volume interactions and other confined systems also exhibit nonequivalence of the ensembles and the emergence of a phase transition [111–113]. Adhesion and fracture processes show nonequivalence as well [76,77,79]. When discussing duality between statistical ensembles, therefore, one must always remember the possibility of having equivalence or nonequivalence between them, in the thermodynamic limit.

Concerning the equilibrium statistical mechanics, the most accurate method is based on the evaluation of the partition function in closed form, without approximations. This can be performed only for a simple system, such as the one considered in this work. When multistable energies are considered, this technique can be applied analytically only in some particular cases, but it can be always implemented numerically with different computational approaches.

If we want to adopt analytical techniques even in the presence of multistable energies, we must use the method of spin variables, which adds a set of discrete variables to the phase space but simplifies the calculation of integrals in partition functions. Indeed, in this case, the partition function is defined by summing over all the spin variable and integrating over all the continuous coordinates (see Equations (22) and (58)). The integral are simpler since they are calculated over the whole domain for quadratic energy well (Gaussian integration). The results are quite accurate if excessive temperatures are not considered (see Figures 5 and 11 for the force-extension response and the spin variable behavior in both ensembles).

Often, to better understand the transition strategies of multistable systems, it is advantageous to perform a zero-temperature energy study. This makes it easier to identify the thresholds of applied forces or extensions responsible for the system's internal transitions. This procedure was also applied in the examples studied in this work and enabled us to better understand the transition thresholds and especially to show that they often do not take on the values that one might imagine with simple mechanical observations of the system. This point in the analysis performed is perfectly dual between the two Helmholtz and Gibbs ensembles.

A last approach concerns the out-of-equilibrium regime for our system. To deal with this problem, we can mention two alternative approach. The first one, the most rigorous, consists of determining the Hamiltonian function for the system and the corresponding equation of motion. To introduce the out-of-equilibrium statistical mechanics, we can add to this equation suitable noise and viscous friction terms, which are able to mimics the dynamics of the system by considering temperature and rate effect. The resulting stochastic differential equation is the so-called Langevin Equation [102]. The probability density describing the evolution of the system state is associated with the Langevin equation and it is referred to as the Fokker-Planck Equation [102]. For complicated systems, with multistable energies, both the Langevin and Fokker-Planck equations can be solved only through quite costly computational techniques. The second approach, used to simplify the numerical solution of the problem, consists of implementing the transition state theory [98]. This method is based on the definition of the probability to be in a given state and on the rate equations describing these probabilities. The rate constants adopted in these equations are obtained through the Kramers formula based on the energy barrier between the states [99]. A central point of this analysis is therefore to identify the energy barriers in the system (see Figures 6 and 12, corresponding two the two statistical ensembles).
The solution of the rate equations is much simpler than the solution of the Langevin or Fokker-Planck equation. Indeed, in this second approach, we have to solve an ordinary system of differential equations which is more advantageous from the computational point of view. With this method, the combined effects of thermal fluctuations and field velocity applied to the system can be easily observed (see Figures 8 and 14). Furthermore, the emergence of hysteresis loops can be investigated (see Figures 9 and 15).

We can also compare the results obtained within the Helmholtz and the Gibbs ensembles. Concerning the solution obtained through the equilibrium statistical mechanics, we note the following important features. On the one hand, we remark that within the isometric Helmholtz ensemble in the force-extension response, we always see in correspondence to the transition a force peak (see Figures 3 and 5). On the other hand, within the Gibbs isotensional ensemble, we observe a force plateau in correspondence to the breaking transition (see Figures 10 and 11). This is a signature characteristic that is also observed in numerous experiments, where the systems are composed of several units. The Helmholtz response (sawtooth-like, with many peaks) shows that the units experience the transitions progressively in reaction to the increasing extension (non-cooperative process), as observed, e.g., in protein unfolding [62–64]. The Gibbs response (plateau-like) is interpreted by supposing that the transitions occur simultaneously for all the units at a given threshold force (cooperative process), as observed, e.g., in DNA overstretching [62–64]. The duality between peak force and plateau force is explained by the fact that the isometric condition is more constraining, and therefore, the system reacts with a larger force while the isotensional condition permits more freedom to the system, which can, hence, accommodate the force more easily. This trend is also confirmed by the results of out-of-equilibrium statistical mechanics, which also shows important rate effects: the force at the transition (peak or plateau depending on the ensemble considered) is significantly increasing with the speed of the applied mechanical action. This is also a result that corresponds to various force spectroscopy experiments, conducted with via high-speed atomic force microscope [104,105]. The models based on rate equations are also able to describe the hysteretic effect of several multistable system, as discussed for our paradigmatic example. We emphasize that the simple model proposed is useful in showing all aspects of symmetry between the Helmholtz and Gibbs dual statistical ensembles and allows us to explore in detail all possible approaches to study the statistical mechanics of multistable energy systems.

These methods, of course, are important when applied to more complex systems than those analyzed in this paper. In general, we deal with a system composed of an arbitrary number of units or elements, and each of them can be bistable or multistable. This is true for problems of adhesion, fracture, and friction, but also for the unzipping of DNA and RNA or the folding/unfolding of macromolecules, such as proteins [62,76,77,79]. In all these cases, the first step is to write an explicit form of the Hamiltonian function of the system in terms of the generalized coordinates and of the spin variables associated with the units of the system. For complex systems, it is almost always impossible to make an exact study of the partition function in any statistical ensemble considered. For this reason, when we assume the thermodynamic equilibrium, we use the spin-variable approach. In this scheme, all the generalized coordinates and the spin variables belong to the phase space and are integrated and summed in the partition function. Interestingly, the evolution of the spin variables is automatically obtained in terms of the externally applied fields to the structure. When interested in the out-of-equilibrium behavior of the system, the approach discussed for one bistable element can be easily generalized to an arbitrary number of spin variables. Firstly, we have to define the states of the system and we must associate a probability to each of these states. Then, we have to identify, thanks to the Hamiltonian function of the system, the energy barriers among each couple of states. These two points allow the writing of the system of differential equations for all the state probabilities controlled by rate coefficients obtained by the Kramers formula. The (numerical) solution of this system allows the study of the dynamics of the system for an arbitrary loading condition.
The solution of the rate equations is much simpler than the solution of the Langevin or Fokker-Planck equation. Indeed, in this second approach, we have to solve an ordinary system of differential equations which is more advantageous from the computational point of view. With this method, the combined effects of thermal fluctuations and field velocity applied to the system can be easily observed (see Figures 8 and 14). Furthermore, the emergence of hysteresis loops can be investigated (see Figures 9 and 15).

We can also compare the results obtained within the Helmholtz and the Gibbs ensembles. Concerning the solution obtained through the equilibrium statistical mechanics, we note the following important features. On the one hand, we remark that within the isometric Helmholtz ensemble in the force-extension response, we always see in correspondence to the transition a force peak (see Figures 3 and 5). On the other hand, within the Gibbs isotensional ensemble, we observe a force plateau in correspondence to the breaking transition (see Figures 10 and 11). This is a signature characteristic that is also observed in numerous experiments, where the systems are composed of several units. The Helmholtz response (sawtooth-like, with many peaks) shows that the units experience the transitions progressively in reaction to the increasing extension (non-cooperative process), as observed, e.g., in protein unfolding [62–64]. The Gibbs response (plateau-like) is interpreted by supposing that the transitions occur simultaneously for all the units at a given threshold force (cooperative process), as observed, e.g., in DNA overstretching [62–64]. The duality between peak force and plateau force is explained by the fact that the isometric condition is more constraining, and therefore, the system reacts with a larger force while the isotensional condition permits more freedom to the system, which can, hence, accommodate the force more easily. This trend is also confirmed by the results of out-of-equilibrium statistical mechanics, which also shows important rate effects: the force at the transition (peak or plateau depending on the ensemble considered) is significantly increasing with the speed of the applied mechanical action. This is also a result that corresponds to various force spectroscopy experiments, conducted with via high-speed atomic force microscope [104,105]. The models based on rate equations are also able to describe the hysteretic effect of several multistable system, as discussed for our paradigmatic example. We emphasize that the simple model proposed is useful in showing all aspects of symmetry between the Helmholtz and Gibbs dual statistical ensembles and allows us to explore in detail all possible approaches to study the statistical mechanics of multistable energy systems.

These methods, of course, are important when applied to more complex systems than those analyzed in this paper. In general, we deal with a system composed of an arbitrary number of units or elements, and each of them can be bistable or multistable. This is true for problems of adhesion, fracture, and friction, but also for the unzipping of DNA and RNA or the folding/unfolding of macromolecules, such as proteins [62,76,77,79]. In all these cases, the first step is to write an explicit form of the Hamiltonian function of the system in terms of the generalized coordinates and of the spin variables associated with the units of the system. For complex systems, it is almost always impossible to make an exact study of the partition function in any statistical ensemble considered. For this reason, when we assume the thermodynamic equilibrium, we use the spin-variable approach. In this scheme, all the generalized coordinates and the spin variables belong to the phase space and are integrated and summed in the partition function. Interestingly, the evolution of the spin variables is automatically obtained in terms of the externally applied fields to the structure. When interested in the out-of-equilibrium behavior of the system, the approach discussed for one bistable element can be easily generalized to an arbitrary number of spin variables. Firstly, we have to define the states of the system and we must associate a probability to each of these states. Then, we have to identify, thanks to the Hamiltonian function of the system, the energy barriers among each couple of states. These two points allow the writing of the system of differential equations for all the state probabilities controlled by rate coefficients obtained by the Kramers formula. The (numerical) solution of this system allows the study of the dynamics of the system for an arbitrary loading condition.
For the examples presented in this work, the main result is represented by the force-extension relationship obtained both at thermodynamic equilibrium and in the out-of-equilibrium regime. However, when we apply these methodologies to real-world systems, for instance to study adhesion and fracture phenomena, a very important result concerns the emergence of phase transitions characterized by a typical critical temperature [76,77,79]. In fact, for these phenomena, a particular value of the temperature exists for which the system is completely detached (or broken) without the application of external loads (as in the thermal denaturation of nucleic acids). In the simple systems studied in this work, the independence from the temperature on the behavior of the transition from the intact to the broken configuration is the direct consequence of the fact that this system is composed of a single unit. Differently, the temperature does play a crucial role in the transition of the system when this is composed of many units interacting with each other through collective phenomena, at the origin of the critical behavior. Two important generalizations will be studied in the near future. Firstly, it could be interesting to perform the continuum limit of the discrete systems describing adhesion and fracture, in order to study the emergence of critical behaviors in continuum systems, which is an important topic in statistical mechanics. Secondly, it is important to generalize the out-of-equilibrium statistical mechanics to adhesion and fracture phenomena in order to see the possible dynamic phase transitions in these situations. Although the rate equations are approximated, this approach can yield new insights into this crucial issue.

The authors hope that the detailed description of all methods applicable to systems with multibasin energy is useful to foster a better understanding and wider dissemination within the scientific community. Utilizing these methods in relatively simple systems is intended to have a pedagogical nature, thereby facilitating their application to more complex systems.

**Author Contributions:** Conceptualization, A.C. and S.G.; methodology, A.C. and S.G.; software, A.C.; writing—original draft preparation, S.G.; writing—review and editing, S.G.; project administration and funding acquisition, S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Central Lille and Region Hauts-de-France under project MIBAMS (Modeling micro Instabilities in Biophysics And Materials Science), and project STAMENA (STAtistical MEchanics for macromolecular structures of NAnotechnology).

**Data Availability Statement:** The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Cortet, P.-P.; Ciccotti, M.; Vanel, L. Imaging the stick-slip peeling of an adhesive tape under a constant load. *J. Stat. Mech.* **2007**, P03005. [CrossRef]

2. Dalbe, M.-J.; Cortet, P.-P.; Ciccotti, M.; Vanel, L.; Santucci, S. Multiscale Stick-Slip Dynamics of Adhesive Tape Peeling. *Phys. Rev. Lett.* **2015**, *115*, 128301. [CrossRef]

3. Maddalena, F.; Percivale, D.; Puglisi, G.; Truskinovsky, L. Mechanics of reversible unzipping. *Contin. Mech. Thermodyn.* **2009**, *21*, 251. [CrossRef]

4. Puglisi, G.; Truskinovsky, L. Cohesion-decohesion asymmetry in geckos. *Phys. Rev. E* **2013**, *87*, 032714. [CrossRef]

5. Katz, S.; Givli, S. Solitary waves in a bistable lattice. *Extrem. Mech. Lett.* **2018**, *22*, 106–111. [CrossRef]

6. Katz, S.; Givli, S. Solitary waves in a nonintegrable chain with double-well potentials. *Phys. Rev. E* **2019**, *100*, 032209. [CrossRef]

7. Harne, R.L.; Schoemaker, M.E.; Dussault, B.E.; Wang, K.W. Wave heave energy conversion using modular multistability. *Appl. Energy* **2014**, *130*, 148–156. [CrossRef]

8. Harne, R.L.; Schoemaker, M.E.; Wang, K.W. Multistable chain for ocean wave vibration energy harvesting. *Proc. of SPIE Act. Passiv. Smart Struct. Integr. Syst.* **2014**, *9057*, 1–12. [CrossRef]

9. Benichou, I.; Givli, S. Structures undergoing discrete phase transformation. *J. Mech. Phys. Solids* **2013**, *61*, 94–113. [CrossRef]

10. Caruel, M.; Allain, J.-M.; Truskinovsky, L. Mechanics of collective unfolding. *J. Mech. Phys. Solids* **2015**, *76*, 237–259. [CrossRef]

11. Caruel, M.; Truskinovsky, L. Bi-stability resistant to fluctuations. *J. Mech. Phys. Solids* **2017**, *109*, 117–141. [CrossRef]

12. Puglisi, G. Hysteresis in multi-stable lattices with non-local interactions. *J. Mech. Phys. Solids* **2006**, *54*, 2060–2088. [CrossRef]
For the examples presented in this work, the main result is represented by the force-extension relationship obtained both at thermodynamic equilibrium and in the out-of-equilibrium regime. However, when we apply these methodologies to real-world systems, for instance to study adhesion and fracture phenomena, a very important result concerns the emergence of phase transitions characterized by a typical critical temperature [76,77,79]. In fact, for these phenomena, a particular value of the temperature exists for which the system is completely detached (or broken) without the application of external loads (as in the thermal denaturation of nucleic acids). In the simple systems studied in this work, the independence from the temperature on the behavior of the transition from the intact to the broken configuration is the direct consequence of the fact that this system is composed of a single unit. Differently, the temperature does play a crucial role in the transition of the system when this is composed of many units interacting with each other through collective phenomena, at the origin of the critical behavior. Two important generalizations will be studied in the near future. Firstly, it could be interesting to perform the continuum limit of the discrete systems describing adhesion and fracture, in order to study the emergence of critical behaviors in continuum systems, which is an important topic in statistical mechanics. Secondly, it is important to generalize the out-of-equilibrium statistical mechanics to adhesion and fracture phenomena in order to see the possible dynamic phase transitions in these situations. Although the rate equations are approximated, this approach can yield new insights into this crucial issue.

The authors hope that the detailed description of all methods applicable to systems with multibasin energy is useful to foster a better understanding and wider dissemination within the scientific community. Utilizing these methods in relatively simple systems is intended to have a pedagogical nature, thereby facilitating their application to more complex systems.

**Author Contributions:** Conceptualization, A.C. and S.G.; methodology, A.C. and S.G.; software, A.C.; writing—original draft preparation, S.G.; writing—review and editing, S.G.; project administration and funding acquisition, S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Central Lille and Region Hauts-de-France under project MIBAMS (Modeling micro Instabilities in Biophysics And Materials Science), and project STAMENA (STAtistical MEchanics for macromolecular structures of NAnotechnology).

**Data Availability Statement:** The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## References

1. Cortet, P.-P.; Ciccotti, M.; Vanel, L. Imaging the stick-slip peeling of an adhesive tape under a constant load. *J. Stat. Mech.* **2007**, P03005. [CrossRef]

2. Dalbe, M.-J.; Cortet, P.-P.; Ciccotti, M.; Vanel, L.; Santucci, S. Multiscale Stick-Slip Dynamics of Adhesive Tape Peeling. *Phys. Rev. Lett.* **2015**, *115*, 128301. [CrossRef]

3. Maddalena, F.; Percivale, D.; Puglisi, G.; Truskinovsky, L. Mechanics of reversible unzipping. *Contin. Mech. Thermodyn.* **2009**, *21*, 251. [CrossRef]

4. Puglisi, G.; Truskinovsky, L. Cohesion-decohesion asymmetry in geckos. *Phys. Rev. E* **2013**, *87*, 032714. [CrossRef]

5. Katz, S.; Givli, S. Solitary waves in a bistable lattice. *Extrem. Mech. Lett.* **2018**, *22*, 106–111. [CrossRef]

6. Katz, S.; Givli, S. Solitary waves in a nonintegrable chain with double-well potentials. *Phys. Rev. E* **2019**, *100*, 032209. [CrossRef]

7. Harne, R.L.; Schoemaker, M.E.; Dussault, B.E.; Wang, K.W. Wave heave energy conversion using modular multistability. *Appl. Energy* **2014**, *130*, 148–156. [CrossRef]

8. Harne, R.L.; Schoemaker, M.E.; Wang, K.W. Multistable chain for ocean wave vibration energy harvesting. *Proc. of SPIE Act. Passiv. Smart Struct. Integr. Syst.* **2014**, 9057, 1–12. [CrossRef]

9. Benichou, I.; Givli, S. Structures undergoing discrete phase transformation. *J. Mech. Phys. Solids* **2013**, *61*, 94–113. [CrossRef]

10. Caruel, M.; Allain, J.-M.; Truskinovsky, L. Mechanics of collective unfolding. *J. Mech. Phys. Solids* **2015**, *76*, 237–259. [CrossRef]

11. Caruel, M.; Truskinovsky, L. Bi-stability resistant to fluctuations. *J. Mech. Phys. Solids* **2017**, *109*, 117–141. [CrossRef]

12. Puglisi, G. Hysteresis in multi-stable lattices with non-local interactions. *J. Mech. Phys. Solids* **2006**, *54*, 2060–2088. [CrossRef]
13. Puglisi, G.; Truskinovsky, L. Thermodynamics of rate-independent plasticity. *J. Mech. Phys. Solids* **2005**, *53*, 655–679. [CrossRef]

14. Blumberg Selinger, R.L.; Wang, Z.-G.; Gelbart, W.M.; Ben-Shaul, A. Statistical-thermodynamic approach to fracture. *Phys. Rev. A* **1991**, *43*, 4396–4400. [CrossRef] [PubMed]

15. Borja da Rocha, H.; Truskinovsky, L. Equilibrium unzipping at finite temperature. *Arch. Appl. Mech.* **2019**, *89*, 535–544. [CrossRef]

16. Borja da Rocha, H.; Truskinovsky, L. Rigidity-Controlled Crossover: From Spinodal to Critical Failure. *Phys. Rev. Lett.* **2020**, *124*, 015501. [CrossRef] [PubMed]

17. Buche, M.R.; Grutzik, S.J. Statistical mechanical model for crack growth. *Phys. Rev. E* **2024**, *109*, 015001. [CrossRef] [PubMed]

18. Manini, N.; Mistura, G.; Paolicelli, G.; Tosatti, E.; Vanossi, A. Current trends in the physics of nanoscale friction. *Adv. Phys. X* **2017**, *2*, 569–590. [CrossRef]

19. Vanossi, A.; Manini, N.; Urbakh, M.; Zapperi, S.; Tosatti, E. Colloquium: Modeling friction: From nanoscale to mesoscale. *Rev. Mod. Phys.* **2013**, *85*, 529–552. [CrossRef]

20. Bleha, T.; Cifra, P. Stretching and compression of DNA by external forces under nanochannel confinement. *Soft Matter* **2018**, *14*, 1247–1259. [CrossRef]

21. Bleha, T.; Cifra, P. Energy/entropy partition of force at DNA stretching. *Biopolymers* **2022**, *113*, e23487. [CrossRef] [PubMed]

22. Prados, A.; Carpio, A.; Bonilla, L.L. Sawtooth patterns in force-extension curves of biomolecules: An equilibrium-statistical-mechanics theory. *Phys. Rev. E* **2013**, *88*, 012704. [CrossRef] [PubMed]

23. Bonilla, L.L.; Carpio, A.; Prados, A. Theory of force-extension curves for modular proteins and DNA hairpins. *Phys. Rev. E* **2015**, *91*, 052712. [CrossRef] [PubMed]

24. De Tommasi, D.; Millardi, N.; Puglisi, G.; Saccomandi, G. An energetic model for macromolecules unfolding in stretching experiments. *J. R. Soc. Interface* **2013**, *10*, 20130651. [CrossRef] [PubMed]

25. Dudko, O.K.; Mathé, J.; Szabo, A.; Meller, A.; Hummer, G. Extracting Kinetics from Single-Molecule Force Spectroscopy: Nanopore Unzipping of DNA Hairpins. *Biophys. J.* **2007**, *92*, 4188–4195. [CrossRef]

26. Dudko, O.K. Decoding the mechanical fingerprints of biomolecules. *Q. Rev. Biophys.* **2016**, *49*, e3. [CrossRef] [PubMed]

27. Hoffmann, T.; Dougan, L. Single molecule force spectroscopy using polyproteins. *Chem. Soc. Rev.* **2012**, *41*, 4781–4796. [CrossRef] [PubMed]

28. Hughes, M.L.; Dougan, L. The physics of pulling polyproteins: A review of single molecule force spectroscopy using the AFM to study protein unfolding. *Rep. Prog. Phys.* **2016**, *79*, 076601. [CrossRef] [PubMed]

29. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F.; Colombo, L. Two-state theory of single-molecule stretching experiments. *Phys. Rev. E* **2013**, *87*, 032705. [CrossRef]

30. Benetatos, P.; Razbin, M. Orientational Fluctuations and Bimodality in Semiflexible Nunchucks. *Polymers* **2021**, *13*, 2031. [CrossRef]

31. Bell, G.I. Models for the Specific Adhesion of Cells to Cells. *Science* **1978**, *200*, 618–627. [CrossRef] [PubMed]

32. Bell, G.I.; Dembo, M.; Bongr, P. Cell adhesion. Competition between nonspecific repulsion and specific bonding. *Biophys. J.* **1984**, *45*, 1051–1064. [CrossRef] [PubMed]

33. Gumbiner, B.M. Cell Adhesion: The Molecular Basis of Tissue Architecture and Morphogenesis. *Cell* **1996**, *84*, 345–357. [CrossRef] [PubMed]

34. Schwarz, U.S.; Safran, S.A. Physics of adherent cells. *Rev. Mod. Phys.* **2013**, *85*, 1327–1381. [CrossRef]

35. Hyeon, C.; Thirumalai, D. Mechanical unfolding of RNA hairpins. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 6789–6794. [CrossRef] [PubMed]

36. Mishra, G.; Giri, D.; Li, M.S.; Kumar, S. Role of loop entropy in the force induced melting of DNA hairpin. *J. Chem. Phys.* **2011**, *135*, 035102. [CrossRef] [PubMed]

37. Stephenson, W.; Keller, S.; Santiago, R.; Albrecht, J.E.; Asare-Okai, P.N.; Tenenbaum, S.A.; Zuker, M.; Li, P.T.X. Combining temperature and force to study folding of an RNA hairpin. *Phys. Chem. Chem. Phys.* **2014**, *16*, 906–917. [CrossRef] [PubMed]

38. De Lorenzo, S.; Ribezzi-Crivellari, M.; Arias-Gonzalez, J.R.; Smith, S.B.; Ritort, F. A Temperature-Jump Optical Trap for Single-Molecule Manipulation. *Biophys. J.* **2015**, *108*, 2854–2864. [CrossRef] [PubMed]

39. Hill, T.L. Theory of Muscular Contraction Extended to Groups of Actin Sites. *Proc. Natl. Acad. Sci. USA* **1973**, *70*, 2732–2736. [CrossRef]

40. Huxley, A.F.; Simmons, R.M. Proposed Mechanism of Force Generation in Striated Muscle. *Nature* **1971**, *233*, 533–538. [CrossRef]

41. Caruel, M.; Allain, J.-M.; Truskinovsky, L. Muscle as a Metamaterial Operating Near a Critical Point. *Phys. Rev. Lett.* **2013**, *110*, 248103. [CrossRef] [PubMed]

42. Caruel, M.; Moireau, P.; Chapelle, D. Stochastic modeling of chemical-mechanical coupling in striated muscles. *Biomech. Model. Mechanobiol.* **2019**, *18*, 563–587. [CrossRef] [PubMed]

43. Caruel, M.; Truskinovsky, L. Statistical mechanics of the Huxley-Simmons model. *Phys. Rev. E* **2016**, *93*, 062407. [CrossRef] [PubMed]

44. Caruel, M.; Truskinovsky, L. Physics of muscle contraction. *Rep. Prog. Phys.* **2018**, *81*, 036602. [CrossRef] [PubMed]

45. Kittel, C. Phase Transition of a Molecular Zipper. *Am. J. Phys.* **1969**, *37*, 917–920. [CrossRef]

46. Peyrard, M.; Bishop, A.R. Statistical mechanics of a nonlinear model for DNA denaturation. *Phys. Rev. Lett.* **1989**, *62*, 2755–2758. [CrossRef] [PubMed]

47. Peyrard, M. Nonlinear dynamics and statistical physics of DNA. *Nonlinearity* **2004**, *17*, R1. [CrossRef]
13. Puglisi, G.; Truskinovsky, L. Thermodynamics of rate-independent plasticity. *J. Mech. Phys. Solids* **2005**, *53*, 655–679. [CrossRef]

14. Blumberg Selinger, R.L.; Wang, Z.-G.; Gelbart, W.M.; Ben-Shaul, A. Statistical-thermodynamic approach to fracture. *Phys. Rev. A* **1991**, *43*, 4396–4400. [CrossRef] [PubMed]

15. Borja da Rocha, H.; Truskinovsky, L. Equilibrium unzipping at finite temperature. *Arch. Appl. Mech.* **2019**, *89*, 535–544. [CrossRef]

16. Borja da Rocha, H.; Truskinovsky, L. Rigidity-Controlled Crossover: From Spinodal to Critical Failure. *Phys. Rev. Lett.* **2020**, *124*, 015501. [CrossRef] [PubMed]

17. Buche, M.R.; Grutzik, S.J. Statistical mechanical model for crack growth. *Phys. Rev. E* **2024**, *109*, 015001. [CrossRef] [PubMed]

18. Manini, N.; Mistura, G.; Paolicelli, G.; Tosatti, E.; Vanossi, A. Current trends in the physics of nanoscale friction. *Adv. Phys. X* **2017**, *2*, 569–590. [CrossRef]

19. Vanossi, A.; Manini, N.; Urbakh, M.; Zapperi, S.; Tosatti, E. Colloquium: Modeling friction: From nanoscale to mesoscale. *Rev. Mod. Phys.* **2013**, *85*, 529–552. [CrossRef]

20. Bleha, T.; Cifra, P. Stretching and compression of DNA by external forces under nanochannel confinement. *Soft Matter* **2018**, *14*, 1247–1259. [CrossRef]

21. Bleha, T.; Cifra, P. Energy/entropy partition of force at DNA stretching. *Biopolymers* **2022**, *113*, e23487. [CrossRef] [PubMed]

22. Prados, A.; Carpio, A.; Bonilla, L.L. Sawtooth patterns in force-extension curves of biomolecules: An equilibrium-statistical-mechanics theory. *Phys. Rev. E* **2013**, *88*, 012704. [CrossRef] [PubMed]

23. Bonilla, L.L.; Carpio, A.; Prados, A. Theory of force-extension curves for modular proteins and DNA hairpins. *Phys. Rev. E* **2015**, *91*, 052712. [CrossRef] [PubMed]

24. De Tommasi, D.; Millardi, N.; Puglisi, G.; Saccomandi, G. An energetic model for macromolecules unfolding in stretching experiments. *J. R. Soc. Interface* **2013**, *10*, 20130651. [CrossRef] [PubMed]

25. Dudko, O.K.; Mathé, J.; Szabo, A.; Meller, A.; Hummer, G. Extracting Kinetics from Single-Molecule Force Spectroscopy: Nanopore Unzipping of DNA Hairpins. *Biophys. J.* **2007**, *92*, 4188–4195. [CrossRef]

26. Dudko, O.K. Decoding the mechanical fingerprints of biomolecules. *Q. Rev. Biophys.* **2016**, *49*, e3. [CrossRef] [PubMed]

27. Hoffmann, T.; Dougan, L. Single molecule force spectroscopy using polyproteins. *Chem. Soc. Rev.* **2012**, *41*, 4781–4796. [CrossRef] [PubMed]

28. Hughes, M.L.; Dougan, L. The physics of pulling polyproteins: A review of single molecule force spectroscopy using the AFM to study protein unfolding. *Rep. Prog. Phys.* **2016**, *79*, 076601. [CrossRef] [PubMed]

29. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F.; Colombo, L. Two-state theory of single-molecule stretching experiments. *Phys. Rev. E* **2013**, *87*, 032705. [CrossRef]

30. Benetatos, P.; Razbin, M. Orientational Fluctuations and Bimodality in Semiflexible Nunchucks. *Polymers* **2021**, *13*, 2031. [CrossRef]

31. Bell, G.I. Models for the Specific Adhesion of Cells to Cells. *Science* **1978**, *200*, 618–627. [CrossRef] [PubMed]

32. Bell, G.I.; Dembo, M.; Bongr, P. Cell adhesion. Competition between nonspecific repulsion and specific bonding. *Biophys. J.* **1984**, *45*, 1051–1064. [CrossRef] [PubMed]

33. Gumbiner, B.M. Cell Adhesion: The Molecular Basis of Tissue Architecture and Morphogenesis. *Cell* **1996**, *84*, 345–357. [CrossRef] [PubMed]

34. Schwarz, U.S.; Safran, S.A. Physics of adherent cells. *Rev. Mod. Phys.* **2013**, *85*, 1327–1381. [CrossRef]

35. Hyeon, C.; Thirumalai, D. Mechanical unfolding of RNA hairpins. *Proc. Natl. Acad. Sci. USA* **2005**, *102*, 6789–6794. [CrossRef] [PubMed]

36. Mishra, G.; Giri, D.; Li, M.S.; Kumar, S. Role of loop entropy in the force induced melting of DNA hairpin. *J. Chem. Phys.* **2011**, *135*, 035102. [CrossRef] [PubMed]

37. Stephenson, W.; Keller, S.; Santiago, R.; Albrecht, J.E.; Asare-Okai, P.N.; Tenenbaum, S.A.; Zuker, M.; Li, P.T.X. Combining temperature and force to study folding of an RNA hairpin. *Phys. Chem. Chem. Phys.* **2014**, *16*, 906–917. [CrossRef] [PubMed]

38. De Lorenzo, S.; Ribezzi-Crivellari, M.; Arias-Gonzalez, J.R.; Smith, S.B.; Ritort, F. A Temperature-Jump Optical Trap for Single-Molecule Manipulation. *Biophys. J.* **2015**, *108*, 2854–2864. [CrossRef] [PubMed]

39. Hill, T.L. Theory of Muscular Contraction Extended to Groups of Actin Sites. *Proc. Natl. Acad. Sci. USA* **1973**, *70*, 2732–2736. [CrossRef]

40. Huxley, A.F.; Simmons, R.M. Proposed Mechanism of Force Generation in Striated Muscle. *Nature* **1971**, *233*, 533–538. [CrossRef]

41. Caruel, M.; Allain, J.-M.; Truskinovsky, L. Muscle as a Metamaterial Operating Near a Critical Point. *Phys. Rev. Lett.* **2013**, *110*, 248103. [CrossRef] [PubMed]

42. Caruel, M.; Moireau, P.; Chapelle, D. Stochastic modeling of chemical-mechanical coupling in striated muscles. *Biomech. Model. Mechanobiol.* **2019**, *18*, 563–587. [CrossRef] [PubMed]

43. Caruel, M.; Truskinovsky, L. Statistical mechanics of the Huxley-Simmons model. *Phys. Rev. E* **2016**, *93*, 062407. [CrossRef] [PubMed]

44. Caruel, M.; Truskinovsky, L. Physics of muscle contraction. *Rep. Prog. Phys.* **2018**, *81*, 036602. [CrossRef] [PubMed]

45. Kittel, C. Phase Transition of a Molecular Zipper. *Am. J. Phys.* **1969**, *37*, 917–920. [CrossRef]

46. Peyrard, M.; Bishop, A.R. Statistical mechanics of a nonlinear model for DNA denaturation. *Phys. Rev. Lett.* **1989**, *62*, 2755–2758. [CrossRef] [PubMed]

47. Peyrard, M. Nonlinear dynamics and statistical physics of DNA. *Nonlinearity* **2004**, *17*, R1. [CrossRef]
48. Grinza, P.; Mossa, A. Topological Origin of the Phase Transition in a Model of DNA Denaturation. *Phys. Rev. Lett.* **2004**, *92*, 158102.
[CrossRef] [PubMed]

49. Palmeri, J.; Manghi, M.; Destainville, N. Thermal denaturation of fluctuating finite DNA chains: The role of bending rigidity in bubble nucleation. *Phys. Rev. E* **2008**, *77*, 011913. [CrossRef]

50. Perret, G.; Lacornerie, T.; Manca, F.; Giordano, S.; Kumemura, M.; Lafitte, N.; Jalabert, L.; Tarhan, M.C.; Lartigau, E.F.; Cleri, F.; et al. Real-time mechanical characterization of DNA degradation under therapeutic X-rays and its theoretical modeling. *Microsyst. Nanoeng.* **2016**, *2*, 16062. [CrossRef]

51. Perret, G.; Lacornerie, T.; Manca, F.; Giordano, S.; Kumemura, M.; Lafitte, N.; Jalabert, L.; Tarhan, M.C.; Lartigau, E.F.; Cleri, F.; et al. Mesure de la dégradation bio-mécanique d'une fibre d'ADN sous l'effet des rayons X thérapeutiques. *Med. Sci.* **2017**, *33*, 1026–1029. [CrossRef] [PubMed]

52. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F. Scaling Shift in Multicracked Fiber Bundles. *Phys. Rev. Lett.* **2014**, *113*, 255501.
[CrossRef]

53. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F. Stochastic mechanical degradation of multi-cracked fiber bundles with elastic and viscous interactions. *Eur. Phys. J. E* **2015**, *38*, 44. [CrossRef] [PubMed]

54. Buche, M.R.; Silberstein, M.N. Chain breaking in the statistical mechanical constitutive theory of polymer networks. *J. Mech. Phys. Solids* **2021**, *156*, 104593. [CrossRef]

55. Bormuth, V.; Varga, V.; Howard, J.; Schäffer, E. Protein Friction Limits Diffusive and Directed Movements of Kinesin Motors on Microtubules. *Science* **2009**, *325*, 870–873. [CrossRef] [PubMed]

56. Sahli, R.; Pallares, G.; Ducottet, C.; Ben Ali, I.E.; Al Akhrass, S.; Guibert, M.; Scheibert, J. Evolution of real contact area under shear and the value of static friction of soft materials. *Proc. Natl. Acad. Sci. USA* **2018**, *115*, 471–476. [CrossRef] [PubMed]

57. Liamas, E.; Connell, S.D.; Ramakrishna, S.N.; Sarkar, A. Probing the frictional properties of soft materials at the nanoscale. *Nanoscale* **2020**, *12*, 2292–2308. [CrossRef] [PubMed]

58. Daub, E.G.; Carlson, J.M. Friction, Fracture, and Earthquakes. *Annu. Rev. Condens. Matter Phys.* **2010**, *1*, 397–418. [CrossRef]

59. Gibbs, J.W. *Elementary Principles in Statistical Mechanics*; Charles Scribner's Sons: New York, NY, USA, 1902.

60. Weiner, J.H. *Statistical Mechanics of Elasticity*; Dover Publications Inc.: New York, NY, USA, 2002.

61. Zwanzig, R. *Nonequilibrium Statistical Mechanics*; Oxford University Press: New York, NY, USA, 2001.

62. Giordano, S. Spin variable approach for the statistical mechanics of folding and unfolding chains. *Soft Matter* **2017**, *13*, 6877–6893.
[CrossRef]

63. Bellino, L.; Florio, G.; Puglisi, G. The influence of device handles in single-molecule experiments. *Soft Matter* **2019**, *15*, 8680–8690.
[CrossRef]

64. Florio, G.; Puglisi, G. Unveiling the influence of device stiffness in single macromolecule unfolding. *Sci. Rep.* **2019**, *9*, 4997.
[CrossRef] [PubMed]

65. Benedito, M.; Giordano, S. Thermodynamics of small systems with conformational transitions: The case of two-state freely jointed chains with extensible units. *J. Chem. Phys.* **2018**, *149*, 054901. [CrossRef] [PubMed]

66. Benedito, M.; Giordano, S. Isotensional and isometric force-extension response of chains with bistable units and Ising interactions. *Phys. Rev. E* **2018**, *98*, 052146. [CrossRef]

67. Benedito, M.; Manca, F.; Giordano, S. Full Statistics of Conjugated Thermodynamic Ensembles in Chains of Bistable Units. *Inventions* **2019**, *4*, 19. [CrossRef]

68. Benedito, M.; Giordano, S. Unfolding pathway and its identifiability in heterogeneous chains of bistable units. *Phys. Lett. A* **2020**, *384*, 126124. [CrossRef]

69. Benedito, M.; Manca, F.; Palla, P.L.; Giordano, S. Rate-dependent force–extension models for single-molecule force spectroscopy experiments. *Phys. Biol.* **2020**, *17*, 056002. [CrossRef] [PubMed]

70. Benichou, I.; Givli, S. Rate Dependent Response of Nanoscale Structures Having a Multiwell Energy Landscape. *Phys. Rev. Lett.* **2015**, *114*, 095504. [CrossRef] [PubMed]

71. Benichou, I.; Zhang, Y.; Dudko, O.K.; Givli, S. The rate dependent response of a bistable chain at finite temperature. *J. Mech. Phys. Solids* **2016**, *95*, 44–63. [CrossRef]

72. Bellino, L.; Florio, G.; Giordano, S.; Puglisi, G. On the competition between interface energy and temperature in phase transition phenomena. *Appl. Eng. Sci.* **2020**, *2*, 100009. [CrossRef]

73. Cannizzo, A.; Bellino, L.; Florio, G.; Puglisi, G.; Giordano, S. Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy. *Eur. Phys. J. Plus* **2022**, *137*, 569. [CrossRef]

74. Giordano, S. Statistical mechanics of rate-independent stick-slip on a corrugated surface composed of parabolic wells. *Contin. Mech. Thermodyn.* **2022**, *34*, 1343–1372. [CrossRef]

75. Giordano, S. Temperature dependent model for the quasistatic stick-slip process on a soft substrate. *Soft Matter* **2023**, *19*, 1813.
[CrossRef]

76. Florio, G.; Puglisi, G.; Giordano, S. Role of temperature in the decohesion of an elastic chain tethered to a substrate by onsite breakable links. *Phys. Rev. Res.* **2020**, *2*, 033227. [CrossRef]

77. Cannizzo, A.; Florio, G.; Puglisi, G.; Giordano, S. Temperature controlled decohesion regimes of an elastic chain adhering to a fixed substrate by softening and breakable bonds. *J. Phys. A Math. Theor.* **2021**, *54*, 445001. [CrossRef]
48. Grinza, P.; Mossa, A. Topological Origin of the Phase Transition in a Model of DNA Denaturation. *Phys. Rev. Lett.* **2004**, *92*, 158102.
[CrossRef] [PubMed]

49. Palmeri, J.; Manghi, M.; Destainville, N. Thermal denaturation of fluctuating finite DNA chains: The role of bending rigidity in bubble nucleation. *Phys. Rev. E* **2008**, *77*, 011913. [CrossRef]

50. Perret, G.; Lacornerie, T.; Manca, F.; Giordano, S.; Kumemura, M.; Lafitte, N.; Jalabert, L.; Tarhan, M.C.; Lartigau, E.F.; Cleri, F.; et al. Real-time mechanical characterization of DNA degradation under therapeutic X-rays and its theoretical modeling. *Microsyst. Nanoeng.* **2016**, *2*, 16062. [CrossRef]

51. Perret, G.; Lacornerie, T.; Manca, F.; Giordano, S.; Kumemura, M.; Lafitte, N.; Jalabert, L.; Tarhan, M.C.; Lartigau, E.F.; Cleri, F.; et al. Mesure de la dégradation bio-mécanique d'une fibre d'ADN sous l'effet des rayons X thérapeutiques. *Med. Sci.* **2017**, *33*, 1026–1029. [CrossRef] [PubMed]

52. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F. Scaling Shift in Multicracked Fiber Bundles. *Phys. Rev. Lett.* **2014**, *113*, 255501.
[CrossRef]

53. Manca, F.; Giordano, S.; Palla, P.L.; Cleri, F. Stochastic mechanical degradation of multi-cracked fiber bundles with elastic and viscous interactions. *Eur. Phys. J. E* **2015**, *38*, 44. [CrossRef] [PubMed]

54. Buche, M.R.; Silberstein, M.N. Chain breaking in the statistical mechanical constitutive theory of polymer networks. *J. Mech. Phys. Solids* **2021**, *156*, 104593. [CrossRef]

55. Bormuth, V.; Varga, V.; Howard, J.; Schäffer, E. Protein Friction Limits Diffusive and Directed Movements of Kinesin Motors on Microtubules. *Science* **2009**, *325*, 870–873. [CrossRef] [PubMed]

56. Sahli, R.; Pallares, G.; Ducottet, C.; Ben Ali, I.E.; Al Akhrass, S.; Guibert, M.; Scheibert, J. Evolution of real contact area under shear and the value of static friction of soft materials. *Proc. Natl. Acad. Sci. USA* **2018**, *115*, 471–476. [CrossRef] [PubMed]

57. Liamas, E.; Connell, S.D.; Ramakrishna, S.N.; Sarkar, A. Probing the frictional properties of soft materials at the nanoscale. *Nanoscale* **2020**, *12*, 2292–2308. [CrossRef] [PubMed]

58. Daub, E.G.; Carlson, J.M. Friction, Fracture, and Earthquakes. *Annu. Rev. Condens. Matter Phys.* **2010**, *1*, 397–418. [CrossRef]

59. Gibbs, J.W. *Elementary Principles in Statistical Mechanics*; Charles Scribner's Sons: New York, NY, USA, 1902.

60. Weiner, J.H. *Statistical Mechanics of Elasticity*; Dover Publications Inc.: New York, NY, USA, 2002.

61. Zwanzig, R. *Nonequilibrium Statistical Mechanics*; Oxford University Press: New York, NY, USA, 2001.

62. Giordano, S. Spin variable approach for the statistical mechanics of folding and unfolding chains. *Soft Matter* **2017**, *13*, 6877–6893.
[CrossRef]

63. Bellino, L.; Florio, G.; Puglisi, G. The influence of device handles in single-molecule experiments. *Soft Matter* **2019**, *15*, 8680–8690.
[CrossRef]

64. Florio, G.; Puglisi, G. Unveiling the influence of device stiffness in single macromolecule unfolding. *Sci. Rep.* **2019**, *9*, 4997.
[CrossRef] [PubMed]

65. Benedito, M.; Giordano, S. Thermodynamics of small systems with conformational transitions: The case of two-state freely jointed chains with extensible units. *J. Chem. Phys.* **2018**, *149*, 054901. [CrossRef] [PubMed]

66. Benedito, M.; Giordano, S. Isotensional and isometric force-extension response of chains with bistable units and Ising interactions. *Phys. Rev. E* **2018**, *98*, 052146. [CrossRef]

67. Benedito, M.; Manca, F.; Giordano, S. Full Statistics of Conjugated Thermodynamic Ensembles in Chains of Bistable Units. *Inventions* **2019**, *4*, 19. [CrossRef]

68. Benedito, M.; Giordano, S. Unfolding pathway and its identifiability in heterogeneous chains of bistable units. *Phys. Lett. A* **2020**, *384*, 126124. [CrossRef]

69. Benedito, M.; Manca, F.; Palla, P.L.; Giordano, S. Rate-dependent force–extension models for single-molecule force spectroscopy experiments. *Phys. Biol.* **2020**, *17*, 056002. [CrossRef] [PubMed]

70. Benichou, I.; Givli, S. Rate Dependent Response of Nanoscale Structures Having a Multiwell Energy Landscape. *Phys. Rev. Lett.* **2015**, *114*, 095504. [CrossRef] [PubMed]

71. Benichou, I.; Zhang, Y.; Dudko, O.K.; Givli, S. The rate dependent response of a bistable chain at finite temperature. *J. Mech. Phys. Solids* **2016**, *95*, 44–63. [CrossRef]

72. Bellino, L.; Florio, G.; Giordano, S.; Puglisi, G. On the competition between interface energy and temperature in phase transition phenomena. *Appl. Eng. Sci.* **2020**, *2*, 100009. [CrossRef]

73. Cannizzo, A.; Bellino, L.; Florio, G.; Puglisi, G.; Giordano, S. Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy. *Eur. Phys. J. Plus* **2022**, *137*, 569. [CrossRef]

74. Giordano, S. Statistical mechanics of rate-independent stick-slip on a corrugated surface composed of parabolic wells. *Contin. Mech. Thermodyn.* **2022**, *34*, 1343–1372. [CrossRef]

75. Giordano, S. Temperature dependent model for the quasistatic stick-slip process on a soft substrate. *Soft Matter* **2023**, *19*, 1813.
[CrossRef]

76. Florio, G.; Puglisi, G.; Giordano, S. Role of temperature in the decohesion of an elastic chain tethered to a substrate by onsite breakable links. *Phys. Rev. Res.* **2020**, *2*, 033227. [CrossRef]

77. Cannizzo, A.; Florio, G.; Puglisi, G.; Giordano, S. Temperature controlled decohesion regimes of an elastic chain adhering to a fixed substrate by softening and breakable bonds. *J. Phys. A Math. Theor.* **2021**, *54*, 445001. [CrossRef]