# Isotensional and isometric force-extension response of chains with bistable units and Ising interactions

Manon Benedito, Stefano Giordano

► To cite this version:

Manon Benedito, Stefano Giordano. Isotensional and isometric force-extension response of chains with bistable units and Ising interactions. Physical Review E, 2018, 98 (5), 10.1103/PhysRevE.98.052146. [hal-02337965](https://hal.science/hal-02337965)

HAL Id: hal-02337965

https://hal.science/hal-02337965

Submitted on 24 Sep 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Isotensional and isometric force-extension response of chains with bistable units and Ising interactions

Manon Benedito and Stefano Giordano*

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

(Dated: October 8, 2018)

The combination of bistability and cooperativity plays a crucial role in several biological and artificial micro- and nano-systems. In particular, the exhaustive understanding of the mechanical response of such systems under the effect of thermal fluctuations is essential to elucidate a rich variety of phenomena. Here, a linear chain composed of elastic units, which are bistable (folded or unfolded) and coupled through an Ising-like interaction, is selected as a case study. We assess the macroscopic thermoelastic response of this chain in terms of its microscopic description. For small systems, far from the thermodynamic limit, this response depends on the applied isometric or isotensional boundary conditions, which correspond to the Helmholtz or Gibbs ensembles of the statistical mechanics, respectively. The theoretical analysis is conducted through the spin variables approach, based on a set of discrete quantities able to identify the folded or unfolded state of the chain units. Eventually, this technique yields closed form expressions for the force-extension curves and the average number of unfolded units, as function of the applied fields. In addition, it allows to unveil a critical behavior of such systems, characterizing the operating regions with negative differential stiffness (spinoidal phase).

## I. INTRODUCTION

The thorough understanding and the tuning of the physical properties of bistable and cooperative systems are the object of extensive research activity [1, 2]. While the bistability represents the tendency for a system to be in only one of two distinct states, the cooperativity concerns systems composed of several units with mutual interactions. It means that, in bistable cooperative systems, each unit can undergo a transition between two states and the transition of one unit affects the transition of the others (because of the interactions among units).

This apparently simple scheme is able to induce very complex behaviors in a large number of physical systems. Without the pretension of being exhaustive, we can cite the snapping and unidirectional waves in elastic metamaterials [3, 4], the mechanics of muscle contraction [5, 6], the magnetic, optical, and structural bistability in spin-crossover nanocrystals [7, 8], the information processing in biochemical reactions [9, 10], the protein folding-unfolding processes [11–14], the DNA over-stretching and denaturation [15–18], and the physics of force-spectroscopy experiments on macromolecules [19–21]. This last example is particularly important since force-spectroscopy experiments, conceived to measure the force-extension relation of a single macromolecule, were able for the first time to directly test the thermodynamics and the statistical mechanics of small systems [22, 23]. In particular, devices like atomic-force microscopes, laser optical tweezers, magnetic tweezers and micro-electro-mechanical systems [24–29] have been employed to investigate proteins [30–32], RNA [33, 34], and DNA [35–40].

While the force-extension response of molecular chains without bistability is considered to be well understood [41–45], the real complexity of chains with bistable units has only been revealed through the introduction of the above force-spectroscopy techniques. In this context, the mechanically induced folding/unfolding of the units of a chain, governed by the conformational transition between two states, has been detected in polypeptides, nucleic acids and other molecules. Notably, for relatively short bistable molecular chains (small systems thermodynamics), the applied boundary conditions play an important role in defining their overall response [46–49].

On the one side, isotensional experiments (conducted at constant applied force by soft devices) correspond to the Gibbs statistical ensemble, and lead to a plateau-like force-extension curve with a threshold force characterizing the synchronized unfolding of all chain units [37, 50–55]. On the other side, isometric experiments (conducted at prescribed displacement by hard devices) represent a realization of the Helmholtz statistical ensemble, and the corresponding force-extension curve shows a sawtooth-like shape, proving that the units unfold sequentially in reaction to the increasing extension [30, 32, 54–61]. In any case, the differences between isotensional and isometric force-extension curves disappear whenever the number of units is very large since, in the thermodynamic limit, the Gibbs and Helmholtz ensembles become statistically equivalent [62, 63]. A different point of view about two-state systems driven by hard or soft devices has been introduced to model plasticity, hysteretic behaviors and martensitic transformations in solids [64–71].

The models proposed in the literature in order to explain the behavior of mechanical bistable chains typically disregarded the interactions among the units, thus neglecting the actual cooperativity of the system. Conversely, we propose here a model where a bistable chain

<footnote>* stefano.giordano@iemn.univ-lille1.fr</footnote>
is explicitly coupled to an Ising spin system, allowing the understanding of the combination of bistability and cooperativity. As a consequence, our model paradigmatically represents most features of the above discussed real systems. Previous attempts to integrate Ising chains in mechanical systems only concerned the dynamics of a single harmonic oscillator coupled to a linear chain of spins [72, 73] and the analysis of ripples in strings [74]. To complement this picture, our analysis fully describes the (entropic and/or enthalpic) elastic response and the units transitions in small systems, within both the Helmholtz and the Gibbs ensembles of the statistical mechanics. For instance, this approach can help the explanation of the heterogeneous series of force-peaks observed in the Helmholtz response during a protein unfolding process. A first reported case concerns the tandem repeats in red cell spectrin [75], where two adjacent units unfold at the same time (cooperatively) because of their strong interaction. Moreover, a similar observation has been made on Filamin A, where domain-domain interactions lead to a hierarchy of unfolding forces that may be properly studied by an Ising scheme [76]. From the continuum mechanics point of view, our model may represent plastic phenomena with nonlocal behaviours, i.e. with interactions among the different regions of the plastic body.

The underlying idea of our method consists in associating to each unit a discrete variable (or spin), able to define the folded or unfolded state of the unit itself. So, the bistable potential function of each unit can be approximated by two quadratic potentials and the switching between them is controlled by the corresponding spin variable (see Fig.1 for details). From the historical point of view, the first model based on a discrete quantity, similar to a spin variable, has been performed to predict the response of skeletal muscles [77, 78]. This method has been recently applied to different two-state systems and molecular chains as well [79–82]. Both Gibbs and Helmholtz ensembles can be considered by the spin variables approach, allowing to draw direct comparisons between isotensional and isometric conditions at thermodynamic equilibrium. While the Gibbs ensemble will be studied by means of the classical transfer matrix method [83], typically adopted for one-dimensional interacting models, the Helmholtz ensemble presents major difficulties and will be approached by exploiting the Laplace transform relationship between the Gibbs and Helmholtz partition functions [84]. It is important to remark that we are studying small systems (with the inequivalence of the ensembles), and we need therefore to determine the exact value of the partition functions and not their approximations holding for a large number of units, as usually done for systems attaining the thermodynamic limit. We provide evidence that the cooperativity, measured by the Ising interaction coefficient, strongly modifies the force-extension response of the chain and its configurational properties. In particular, under isometric conditions, we thoroughly analyze the hierarchy of force peaks as function of the interaction coefficient. To complement the

equilibrium picture, we further characterize the criticality of the spinoidal phase, describing the regions with negative differential stiffness.

The structure of the paper is the following. In Section II, we define the system under investigation and its Hamiltonian function. In Section III and IV, we analyze the behavior of the chain with Ising interactions under isotensional (Gibbs) and isometric (Helmholtz) conditions, respectively. Since the problem of the Helmholtz ensemble is solved here through a semi-analytic procedure, we propose in Sections V, VI and VII additional explicit asymptotic results describing the behavior of the system under weak and strong Ising interactions (ferromagnetic-like and antiferromagnetic-like). Finally, in Section VIII, we generalize our results in order to take account of a finite extensibility of the chain units, and we illustrate its effect on the critical behavior of the system.

## II. THE SYSTEM

We take into consideration a chain of m two-state elements (see Fig.1a), each described by a bistable potential energy with a stable folded state and a metastable unfolded state (see Fig.1b). The two potential wells in Fig.1b can be characterized by the elastic constant $k(S_i)$, the equilibrium length $\ell_0(S_i)$ and the basal energy $v(S_i)$, where $S_i$ is a discrete variable (or spin variable) assuming values in $\{-1, +1\}$, used to distinguish one well from the other. We state that $S_i = +1$ corresponds to unfolded elements, whereas $S_i = -1$ corresponds to folded ones. This description suggests that the bistable energy potential can be represented by two quadratic potentials approximating the real wells of the units (see again Fig.1b) [81]. In this case, the discrete variables belong to the phase space of the system and allow to specify the explored well for each unit. The introduction of the discrete or spin variables also allows the direct implementation of an interaction between adjacent elements of the chain, e.g. described by a classical Ising Hamiltonian. The overall Hamiltonian of this system can be therefore written as

$$ H = -\lambda \sum_{i=1}^{m-1} S_i S_{i+1} - \mu \sum_{i=1}^{m} S_i + \sum_{i=1}^{m} \left[ v(S_i) + \frac{1}{2} k(S_i) \left[ \|\vec{r}_i - \vec{r}_{i-1}\| - \ell_0(S_i) \right]^2 \right]. \quad (1) $$

While the first line in Eq.(1) represents the Ising interaction among spin variables, the second line describes the spring-like behavior of each unit placed between positions $\vec{r}_{i-1}$ and $\vec{r}_i$, $\forall i = 1, ..., m$. We remark that $\lambda > 0$ tries to force all elements to be folded or unfolded (ferromagnetic-like interaction), whereas $\lambda < 0$ tries to force all elements to be alternatively folded and unfolded (antiferromagnetic-like interaction). The parameter $\mu$ is a sort of external field or chemical potential ($\mu > 0$ tries to unfold the domains and $\mu < 0$ tries to fold them).
<figure><img src="image_2.png" /><figcaption>FIG. 1. (a) Chain of m two-state units with Ising interactions. While the first end-terminal α is able to tether the first unit to a given substrate, the second one β is able either to apply a force (Gibbs condition) or a position (Helmholtz condition) to the last unit. (b) Potential energy of a single unit of the chain (dashed black curve). The potential wells are approximated through two parabolic (i.e., quadratic) profiles (solid blue curves), identified by $S_i = -1$ (folded state) and $S_i = +1$ (unfolded state).</figcaption></figure>

A first validation of the spin approach to describe multistable potential energies has been performed in Ref.[81], where we directly compared a real two-state system (dashed black curve in Fig.1b) without Ising interactions with the approximation given by two parabolic profiles governed by the spin variables (solid blue curves in Fig.1b), by obtaining a good agreement provided that $\mathcal{E}_b \ge v(1) - v(-1) \gg K_B T$ (where $\mathcal{E}_b$ is the energy barrier). However, the spin variables approach can be adopted only when we work at thermodynamic equilibrium. As a matter of fact, the quadratic potentials and the associated spin variables are not sufficient to describe the dynamic regime since the relaxation times of the system strongly depend on the energy barriers between the potential wells, which are neglected within our approach. This is a well-known result, encoded within the Kramers rate formula, originally formulated to study chemical reactions [85], and recently generalized for arbitrary systems with nonconvex energy landscapes [86, 87].

In the following, we suppose to embed the system in a thermal bath at the temperature $T$, we consider the system at thermodynamic equilibrium and we study the effects of the Ising interactions on the mechanical and configurational behavior within the Gibbs (applied external force) or the Helmholtz (prescribed end-to-end distance) ensembles (see Fig.1a).

### III. TWO-STATE CHAIN WITH ISING INTERACTIONS: THE GIBBS ENSEMBLE

We consider now the extended Hamiltonian

$$ H_G = -\lambda \sum_{i=1}^{m-1} S_i S_{i+1} - \mu \sum_{i=1}^{m} S_i - \vec{f} \cdot \vec{r}_m \quad (2) \\ + \sum_{i=1}^{m} \left[ v(S_i) + \frac{1}{2} k(S_i) \left[ \|\vec{r}_i - \vec{r}_{i-1}\| - \ell_0(S_i) \right]^2 \right], $$

introduced to deal with the isotational conditions. Here, $\vec{f}$ is the force applied to the last unit, identified by its position $\vec{r}_m$. We suppose that quantities $\vec{r}_i \in \mathbb{R}^3$ and $S_i \in \{-1, +1\} \quad \forall i \in \{1...m\}$ belong to the phase space of the system. Moreover, to fix the ideas, we always consider $\vec{r}_0 = 0$. The statistical mechanics of the system can be introduced by calculating the partition function, as follows

$$ \begin{aligned} Z_G &= \sum_{S_1} \dots \sum_{S_m} \int_P \exp \left[ -\frac{H_G(\{S_i\}, \{\vec{r}_i\})}{K_B T} \right] d\vec{r}_1 \dots d\vec{r}_m \\ &= \sum_{S_1} \dots \sum_{S_m} \exp \left( \frac{\lambda}{K_B T} \sum_{i=1}^{m-1} S_i S_{i+1} \right) \\ & \quad \times \exp \left( \frac{\mu}{K_B T} \sum_{i=1}^{m} S_i \right) \exp \left( -\frac{1}{K_B T} \sum_{i=1}^{m} v(S_i) \right) \\ & \quad \times \int_P \exp \left( -\frac{1}{2} \sum_{i=1}^{m} \frac{k(S_i)}{K_B T} \left[ \|\vec{r}_i - \vec{r}_{i-1}\| - \ell_0(S_i) \right]^2 \right) \\ & \quad \times \exp \left( \frac{\vec{f} \cdot \vec{r}_m}{K_B T} \right) d\vec{r}_1 \dots d\vec{r}_m, \end{aligned} \quad (3) $$

where $\mathcal{P} = \mathbb{R}^{3m}$. The integral $I = \int_P \dots d\vec{r}_1 \dots d\vec{r}_m$, shown in the last two lines of Eq.(3), can be developed by means of the change of variables $\vec{\xi}_1 = \vec{r}_1 - \vec{r}_0$, $\vec{\xi}_2 = \vec{r}_2 - \vec{r}_1$, ..., $\vec{\xi}_m = \vec{r}_m - \vec{r}_{m-1}$, giving

$$ I = \int_P \exp \left( -\frac{1}{2} \sum_{i=1}^{m} \frac{k(S_i)}{K_B T} \left[ \|\vec{\xi}_i\| - \ell_0(S_i) \right]^2 \right) \times \exp \left( \frac{\vec{f}}{K_B T} \cdot \sum_{i=1}^{m} \vec{\xi}_i \right) d\vec{\xi}_1 \dots d\vec{\xi}_m. \quad (4) $$

To further simplify this integral, by exploiting the isotropy of the system, we suppose that $\vec{f} = (0, 0, f)$ and we introduce the spherical coordinates for the vectors
$\vec{\xi}_i$, namely $\vec{\xi}_i = (\xi_i \cos \varphi_i \sin \theta_i, \xi_i \sin \varphi_i \sin \theta_i, \xi_i \cos \theta_i)$. Therefore, we easily obtain $\|\vec{\xi}_i\| = \xi_i$, $\vec{f} \cdot \vec{\xi}_i = f\xi_i \cos \theta_i$ and $d\vec{\xi}_i = \xi_i^2 \sin \theta_i d\xi_i d\varphi_i d\theta_i$, and we get for I the expression

$$ 
\begin{aligned}
I &= \int_D \exp \left\{ -\frac{1}{2} \sum_{i=1}^{m} \frac{k(S_i)}{K_B T} [\xi_i - \ell_0(S_i)]^2 \right\} \\
&\quad \times \exp \left( \sum_{i=1}^{m} \frac{f \xi_i \cos \theta_i}{K_B T} \right) \prod_{i=1}^{m} (\xi_i^2 \sin \theta_i d\xi_i d\varphi_i d\theta_i) \\
&= (2\pi)^m \int_{A^m} \exp \left\{ -\frac{1}{2} \sum_{i=1}^{m} \frac{k(S_i)}{K_B T} [\xi_i - \ell_0(S_i)]^2 \right\} \\
&\quad \times \left[ \prod_{i=1}^{m} \int_C \exp \left( \frac{f \xi_i \cos \theta_i}{K_B T} \right) \sin \theta_i d\theta_i \right] \prod_{i=1}^{m} \xi_i^2 d\xi_i \\
&= (4\pi)^m \prod_{i=1}^{m} \int_A \exp \left\{ -\frac{1}{2} \frac{k(S_i)}{K_B T} [\xi_i - \ell_0(S_i)]^2 \right\} \\
&\quad \times \frac{\sinh \left( \frac{f \xi_i}{K_B T} \right)}{\frac{f \xi_i}{K_B T}} \xi_i^2 d\xi_i,
\end{aligned} 
\quad (5) $$

where $\mathcal{D} = \mathcal{A}^m \times \mathcal{B}^m \times \mathcal{C}^m$ with $\mathcal{A} = (0, +\infty)$, $\mathcal{B} = (0, 2\pi)$, and $\mathcal{C} = (0, \pi)$. The integral over $\mathcal{A}$ can be interpreted as a convolution between the enthalpic response (represented by the exponential term) and the entropic one (represented by the function $\sinh(z)/z$). Anyway, the integral $I$ can be strongly simplified if we make the assumption to deal with a freely jointed chain [84], with elements of fixed lengths, i.e. without elasticity. It is equivalent to say that $k(+1) = k(-1) \to +\infty$. This hypothesis will be removed in a successive section of the paper, where we will study an extensible chain with Ising interactions. If we use the property $\sqrt{\frac{\alpha}{\pi}} e^{-\alpha x^2} \to \delta(x)$ for $\alpha \to \infty$, then we simplify the result for $I$ as follows

$$ I = \prod_{i=1}^{m} \int_{\mathcal{A}} \delta(\xi_i - \ell_0(S_i)) \frac{\sinh\left(\frac{f\ell_0(S_i)}{K_B T}\right) \xi_i^2}{\ell^2} d\xi_i, \quad (6) $$

where we omitted a non-influential multiplicative constant and $\ell = \ell_0(-1)$ corresponds to the length of the folded units. We finally obtain

$$ I = \prod_{i=1}^{m} \frac{\sinh\left(\frac{f\ell_0(S_i)}{K_B T}\right)}{\frac{f\ell_0(S_i)}{K_B T}} \frac{\ell_0^2(S_i)}{\ell^2}. \quad (7) $$

Eventually, the partition function assumes the simpler form

$$ Z_G = \sum_{S_1} \dots \sum_{S_m} \left[ \exp \left( \frac{\lambda}{K_B T} \sum_{i=1}^{m-1} S_i S_{i+1} \right) \right. \\
\times \exp \left( \frac{\mu}{K_B T} \sum_{i=1}^{m} S_i \right) \exp \left( -\frac{1}{K_B T} \sum_{i=1}^{m-1} v(S_i) \right) \\
\left. \times \prod_{i=1}^{m} \frac{\sinh \left( \frac{f \ell_0(S_i)}{K_B T} \right)}{\frac{f \ell_0(S_i)}{K_B T}} \frac{\ell_0^2(S_i)}{\ell^2} \right]. \quad (8) $$

We have now to approach the problem of calculating the sums over the spin variables. To this aim, a more symmetric form of Eq.(8) can be obtained by observing that

$$ \prod_{i=1}^{m} c_i = \sqrt{c_1} \left[ \prod_{i=1}^{m-1} \sqrt{c_i c_{i+1}} \right] \sqrt{c_m}, \quad (9) $$

which is a property valid for real numbers $c_i > 0 \forall i$. Accordingly, we have

$$ 
\begin{aligned}
Z_G = & \sum_{S_1, \dots, S_m} \exp \left( \frac{\mu}{2K_B T} S_1 \right) \exp \left( -\frac{v(S_1)}{2K_B T} \right) \sqrt{c_1} \\
& \times \prod_{i=1}^{m-1} \left\{ \exp \left( \frac{\lambda}{K_B T} S_i S_{i+1} \right) \exp \left[ \frac{\mu(S_i + S_{i+1})}{2K_B T} \right] \right. \\
& \left. \times \exp \left[ -\frac{v(S_i)}{2K_B T} - \frac{v(S_{i+1})}{2K_B T} \right] \sqrt{c_i c_{i+1}} \right\} \\
& \times \exp \left( \frac{\mu}{2K_B T} S_m \right) \exp \left( -\frac{v(S_m)}{2K_B T} \right) \sqrt{c_m}, \quad (10)
\end{aligned}
 $$

where we defined

$$ c_i = \frac{\sinh\left(\frac{f\ell_0(S_i)}{K_B T}\right)}{\frac{f\ell_0(S_i)}{K_B T}} \frac{\ell_0^2(S_i)}{\ell^2}. \quad (11) $$

To further elaborate the partition function, we also consider $\ell_0(-1) = \ell$, $v(-1) = 0$, $\ell_0(+1) = \chi\ell$, and $v(+1) = \Delta E$, where $\chi$ is the ratio between unfolded and folded lengths, and $\Delta E$ is the energy jump between the wells (see Fig.1). We can adopt the technique of the transfer matrix [83] and then we can directly write

$$ Z_G = \tilde{w}^T \mathcal{T}^{m-1} \tilde{w}, \quad (12) $$

where we have

$$ \tilde{w} = \sqrt{2} \left( \frac{\sqrt{p}}{\sqrt{q}} \right), \quad (13) $$

$$ \mathcal{T} = 2 \left( \frac{pe^{\frac{\lambda}{K_B T}}}{\sqrt{pqe^{-\frac{\lambda}{K_B T}}}} \quad \sqrt{pqe^{-\frac{\lambda}{K_B T}}} \right), \quad (14) $$

with the parameters

$$ p = \frac{1}{2} e^{-\frac{\mu}{K_B T}} \frac{\sinh(\xi)}{\xi}, \quad (15) $$

$$ q = \frac{1}{2} e^{\frac{\mu}{K_B T}} \chi \phi \frac{\sinh(\chi\xi)}{\xi}, \quad (16) $$

and the coefficients

$$ \phi = e^{-\frac{\Delta E}{K_B T}}, \quad (17) $$

$$ \xi = \frac{f\ell}{K_B T}, \quad (18) $$

representing the Boltzmann factor calculated with $\Delta E$ and the normalized force. Since we are studying the thermodynamics of small systems (small values of $m$), we need to calculate the exact value of the partition function
<figure><img src="image_3.png" /><figcaption>FIG. 2. Average normalized extension (a) and average number of unfolded units (b) versus the applied normalized force for a chain with Ising interactions under isotensional (Gibbs) conditions. The curves have been obtained with $\Delta E = 7K_B T$, $m=4$, $\chi = 2$, $\mu = 0$, and $\lambda = 0, \pm 1, \pm 2K_B T$. The black dashed curves correspond to the chain without interactions ($\lambda = 0$), the dark red (or dark gray) curves to ferromagnetic-like Ising interactions ($\lambda > 0$) and the orange (or light gray) curves to antiferromagnetic-like interactions ($\lambda < 0$).</figcaption></figure>

given in Eq.(12) and not its approximation evaluated for a large value of $m$, corresponding to the thermodynamic limit. To develop the calculation of $Z_G$, we can simply calculate the eigenvalues of $\mathcal{T}$

$$ \lambda_{1,2} = \frac{1}{x} \left[ (p+q) \pm \sqrt{(p+q)^2 - 4pq(1-x^4)} \right], (19) $$

where we also introduced $x = e^{-\frac{\lambda}{K_B T}}$. We underline that if $\lambda_1$ corresponds to the sign '+' and $\lambda_2$ to the sign '-' then we get $\lambda_1 > \lambda_2 > 0$. Now, we need to explicitly determine the matrix power $\mathcal{T}^{m-1}$. Hence, we use the matrix function theory [88], and we get after straightforward calculations

$$ \mathcal{T}^{m-1} = \frac{\lambda_1^{m-1} - \lambda_2^{m-1}}{\lambda_1 - \lambda_2} \mathcal{T} + \frac{\lambda_1 \lambda_2^{m-1} - \lambda_2 \lambda_1^{m-1}}{\lambda_1 - \lambda_2} \mathcal{I}, (20) $$

where $\mathcal{I}$ is the $2 \times 2$ identity matrix. The partition function then assumes the form

$$ Z_G = \frac{\lambda_1^{m-1} - \lambda_2^{m-1}}{\lambda_1 - \lambda_2} \alpha + \frac{\lambda_1 \lambda_2^{m-1} - \lambda_2 \lambda_1^{m-1}}{\lambda_1 - \lambda_2} \beta, (21) $$

where $\alpha$ and $\beta$ can be obtained through long but straightforward calculations as

$$ \alpha = \bar{w}^T \mathcal{T} \bar{w} = \frac{x(1+x^2)(\lambda_1 + \lambda_2)^2 - 2x\lambda_1\lambda_2}{1+x^2}, (22) $$

$$ \beta = \bar{w}^T \bar{w} = x(\lambda_1 + \lambda_2). \quad (23) $$

Finally, the explicit exact form of the Gibbs partition function is

$$ Z_G(f) = \frac{x}{1+x^2} \left[ \lambda_1^m \left( 1 + x^2 \frac{\lambda_1 + \lambda_2}{\lambda_1 - \lambda_2} \right) + \lambda_2^m \left( 1 - x^2 \frac{\lambda_1 + \lambda_2}{\lambda_1 - \lambda_2} \right) \right], (24) $$

written as function of $x$, $\lambda_1$ and $\lambda_2$. This is the most important result of this section and allows us to determine the mechanical and configurational macroscopic behavior of the whole chain under isotensional conditions. As usual, we obtain the force-extension response as

$$ \langle r \rangle = K_B T \frac{\partial \log Z_G}{\partial f} = K_B T \frac{1}{Z_G} \frac{\partial Z_G}{\partial f}, \quad (25) $$

where $\langle r \rangle$ represents the average value of the extension, measured in the direction of the applied force. We also note that the quantity $\frac{S_i+1}{2}$ gives 0 for folded elements and 1 for unfolded elements. Therefore, we have that $\langle \sum_{i=1}^m \frac{S_i+1}{2} \rangle$ is the average number of unfolded elements. On the other hand, the term $\sum_{i=1}^m v(S_i)$ of the Hamiltonian in Eq.(2) can be also written as $\sum_{i=1}^m \frac{S_i+1}{2} \Delta E$ and, therefore, we have that $\sum_{i=1}^m \frac{S_i+1}{2} = \frac{\partial H_G}{\partial \Delta E}$. It follows that $\langle \sum_{i=1}^m \frac{S_i+1}{2} \rangle$ can be evaluated through the expression

$$ \langle \sum_{i=1}^{m} \frac{S_i + 1}{2} \rangle = \frac{\sum_{\vec{S}} \int \frac{\partial H_G}{\partial \Delta E} e^{-\frac{H_G}{K_B T}} d\vec{r}}{\sum_{\vec{S}} \int e^{-\frac{H_G}{K_B T}} d\vec{r}}, \quad (26) $$

where $\vec{S} = (S_1, ..., S_m)$ and $\vec{r} = (\vec{r}_1, ..., \vec{r}_m)$. Then, Eq.(26) can be simplified to give

$$ \langle \sum_{i=1}^{m} \frac{S_i + 1}{2} \rangle = -K_B T \frac{\partial \log Z_G}{\partial \Delta E} = -K_B T \frac{1}{Z_G} \frac{\partial Z_G}{\partial \Delta E}, \quad (27) $$

which is the final expression for the average value of unfolded domains. It is useful to introduce here the Gibbs free energy of the system $G = -K_B T \log Z_G$. The above expected values can be reformulated in terms of this ther-
modynamic function as follows

$$ \langle r \rangle = - \frac{\partial \mathcal{G}}{\partial f}, \qquad (28) $$

$$ \langle \sum_{i=1}^{m} \frac{S_i + 1}{2} \rangle = \frac{\partial \mathcal{G}}{\partial \Delta E}. \qquad (29) $$

The knowledge of $Z_G$ or $\mathcal{G}$ allows therefore the determination of both the average extension of the chain and the average number of unfolded units as function of the applied force and temperature.

An application of Eqs.(28) and (29) can be found in Fig.2. First of all, we discuss the curves for $\lambda = 0$, i.e. without Ising interactions: in the force-extension curve (Fig.2a), we note a force plateau corresponding to the synchronized unfolding of the $m$ units. This behavior is confirmed by the number of unfolded units (Fig.2b), which shows a transition from 0 to $m$, at the same threshold force as the previously mentioned plateau. This force plateau is the classical result of force-spectroscopy experiments conducted with soft devices [50–55]. The Ising interactions modify this scenario as follows: if $\lambda > 0$, the units are favored to be in the same state and the transition is therefore sharper than the one observed for $\lambda = 0$; on the other hand, if $\lambda < 0$, it is easier to unfold the first units and the response is therefore smoother. It is interesting to note that, for strong antiferromagnetic-like interactions, even with $f = 0$, we can have a positive number of unfolded units because of the energetic compromise between $\lambda$ and $\Delta E$. Indeed, when we observe unfolded units with $f = 0$, the cooperativity effect is stronger than the energy jump between folded and unfolded states. We finally remark that the curves in Fig.2 are valid for any value of $\ell$ and $T$. If, as an example, we consider $\ell = 0.4$nm and $T = 300$K, we get a plateau force at $f = 70$pN, which is coherent, e.g., with the DNA overstretching transition [48].

## IV. TWO-STATE CHAIN WITH ISING INTERACTIONS: THE HELMHOLTZ ENSEMBLE

We consider now the Helmholtz ensemble where the last element of the chain is fixed at a given arbitrary position (isometric condition). The Hamiltonian reads

$$ H_H = \lambda \sum_{i=1}^{m-1} S_i S_{i+1} - \mu \sum_{i=1}^{m} S_i \\ + \sum_{i=1}^{m} \left[ v(S_i) + \frac{1}{2} k(S_i) \left[ \|\vec{r}_i - \vec{r}_{i-1}\| - \ell_0(S_i) \right]^2 \right], \qquad (30) $$

where $\vec{r}_m = \vec{r}$ is fixed. The phase space is therefore composed of $\vec{r}_i \forall i = 1, ..., m-1$ and $S_i \forall i = 1, ..., m$.

Hence, the partition function can be written as

$$ \begin{aligned} Z_H &= \sum_{S_1} \dots \sum_{S_m} \int_Q \exp \left[ -\frac{H_H(\{S_i\}, \{\vec{r}_i\})}{K_B T} \right] d\vec{r}_1 \dots d\vec{r}_{m-1} \\ &= \sum_{S_1} \dots \sum_{S_m} \exp \left( \frac{\lambda}{K_B T} \sum_{i=1}^{m-1} S_i S_{i+1} \right) \\ & \quad \times \exp \left( \frac{\mu}{K_B T} \sum_{i=1}^{m} S_i \right) \exp \left( -\frac{1}{K_B T} \sum_{i=1}^{m} v(S_i) \right) \\ & \quad \times \int_Q \exp \left( -\frac{1}{2} \sum_{i=1}^{m} \frac{k(S_i)}{K_B T} \left[ \|\vec{r}_i - \vec{r}_{i-1}\| - \ell_0(S_i) \right]^2 \right) \\ & \quad \times d\vec{r}_1 \dots d\vec{r}_{m-1}, \end{aligned} \qquad (31) $$

where $Q = \mathbb{R}^{3(m-1)}$. It is not difficult to realize that the calculation of $Z_H$ is much more complicated than the one performed for $Z_G$. Indeed, in this case, we can not apply a simple change of variables in order to factorize the multi-dimensional integral. From the physical point of view, this difficulty depends on the fact that the isometric condition induces an effective interaction among the units, fixing the sum of all vectors $\vec{r}_i - \vec{r}_{i-1}$ (for $i$ from 1 to $m$). Then, we have now the combination of two forms of interaction among the units, the first being implicitly encoded in the isometric condition and the second explicitly implemented through the Ising scheme. An useful technique to cope with this difficulty is the following. By comparing Eqs.(3) and (31), we deduce that the two partition functions $Z_G$ and $Z_H$ are related through a three-dimensional bilateral Laplace transform, as follows

$$ Z_G(\vec{f}) = \int_{\mathbb{R}^3} Z_H(\vec{r}) \exp \left( \frac{\vec{r} \cdot \vec{f}'}{K_B T} \right) d\vec{r}, \qquad (32) $$

where, as usual, we neglect the non-influential multiplicative constants in the partition functions. Moreover, by considering the spherical symmetry of the problem, we easily obtain the inverse relationship

$$ Z_H(r) = \int_{-\infty}^{+\infty} Z_G(i\eta) \frac{\eta}{r} \sin \frac{\eta r}{K_B T} d\eta, \qquad (33) $$

where $Z_G(i\eta)$ is the analytic continuation of the partition function $Z_G(f)$ for the Gibbs ensemble, given in Eq.(24). The integral in Eq.(33) can be simplified by the change of variable $y = \frac{\eta\ell}{K_B T}$, leading to

$$ Z_H(r) = \int_{-\infty}^{+\infty} Z_G \left( i y \frac{K_B T}{\ell} \right) \frac{y}{r} \sin \frac{y r}{\ell} dy, \qquad (34) $$

where, as before, we neglected the non-influential multiplicative constant. Coherently with our assumptions, the variables $p$ and $q$ assume the form

$$ p = \frac{1}{2} e^{-\frac{\mu}{K_B T}} \frac{\sin y}{y} = \frac{\tilde{p}}{y}, \qquad (35) $$

$$ q = \frac{1}{2} e^{\frac{\mu}{K_B T}} \chi \phi \frac{\sin \chi y}{y} = \frac{\tilde{q}}{y}. \qquad (36) $$
Accordingly, the eigenvalues of the transfer matrix become

$$ \lambda_{1,2} = \frac{1}{xy} \left[ \tilde{p} + \tilde{q} \pm \sqrt{(\tilde{p} + \tilde{q})^2 - 4\tilde{p}\tilde{q}(1 - x^4)} \right] = \frac{\tilde{\lambda}_{1,2}}{y}, \quad (37) $$

where $\tilde{p} = py$, $\tilde{q} = qy$, and $\tilde{\lambda}_{1,2} = \lambda_{1,2}y$. Hence, the analytic continuation of the Gibbs partition function becomes

$$ Z_G \left( iy \frac{K_B T}{\ell} \right) = \frac{x}{1+x^2} \frac{1}{y^m} \times \left[ \tilde{\lambda}_1^m \left( 1 + x^2 \frac{\tilde{\lambda}_1 + \tilde{\lambda}_2}{\tilde{\lambda}_1 - \tilde{\lambda}_2} \right) + \tilde{\lambda}_2^m \left( 1 - x^2 \frac{\tilde{\lambda}_1 + \tilde{\lambda}_2}{\tilde{\lambda}_1 - \tilde{\lambda}_2} \right) \right], \quad (38) $$

where, importantly, $\tilde{\lambda}_1$ and $\tilde{\lambda}_2$ depend on $y$ only through $\sin y$ and $\sin \chi y$. In particular, when $\chi$ is an integer (or also a rational number), $Z_G$ is composed of a periodic function of $y$ divided by $y^m$. So we have

$$ Z_G \left( iy \frac{K_B T}{\ell} \right) = \frac{1}{y^m} P(y), \quad (39) $$

where $P(y) = P(y + L_y)$ for a given $L_y$ and we have

$$ P(y) = \frac{x}{1+x^2} \times \left[ \tilde{\lambda}_1^m \left( 1 + x^2 \frac{\tilde{\lambda}_1 + \tilde{\lambda}_2}{\tilde{\lambda}_1 - \tilde{\lambda}_2} \right) + \tilde{\lambda}_2^m \left( 1 - x^2 \frac{\tilde{\lambda}_1 + \tilde{\lambda}_2}{\tilde{\lambda}_1 - \tilde{\lambda}_2} \right) \right]. \quad (40) $$

If we consider integer values of $\chi$, $P(y)$ is periodic with a period of $L_y = 2\pi$, and it can be developed in Fourier series, as follows

$$ P(y) = \sum_{k=-\infty}^{+\infty} C_k e^{iky}, \quad (41) $$

where

$$ C_k = \frac{1}{2\pi} \int_0^{2\pi} P(y)e^{-iky} dy. \quad (42) $$

The values of $C_k$ can be obtained numerically by calculating the integrals through classical numerical techniques (we verified that the simple Simpson's rule is sufficient to obtain accurate results). Once determined the $C_k$ coefficients, the Helmholtz partition function can be obtained analytically as follows. To begin, we have from Eq.(33)

$$ Z_H(r) = -i \int_{\Gamma} Z_G \left( iy \frac{K_B T}{\ell} \right) \frac{y}{r} e^{\frac{ry}{\ell}} dy. \quad (43) $$

Here, we used the Euler formula $e^{i\frac{ry}{\ell}} = \cos\frac{ry}{\ell} + i\sin\frac{ry}{\ell}$, and we observed that the integral with $\cos\frac{ry}{\ell}$ is zero since $Z_G(iy\frac{K_B T}{\ell})$ is an even function of $y$. Moreover, since the function to integrate is regular on the real axis and holomorphic on a strip $|Imy| < M$ for an arbitrary $M \in \mathbb{R}$, we can use the path $\Gamma$ shown in Fig.3. This will be useful to elaborate the partition function integral and to

<figure><img src="image_4.png" /><figcaption>FIG. 3. Definition of the contour $\Gamma$ on the complex plane with an arbitrary radius $\epsilon$.</figcaption></figure>

write it in a form with singularities at the origin. Indeed, we have

$$ Z_H(r) = -i \int_{\Gamma} \frac{P(y)}{y^m} \frac{y}{r} e^{\frac{iry}{\ell}} dy \quad (44) $$

$$ \begin{aligned} &= -i \int_{\Gamma} \frac{1}{y^m} \sum_{k=-\infty}^{+\infty} C_k e^{iky} \frac{y}{r} e^{\frac{iry}{\ell}} dy \\ &= -i \sum_{k=-\infty}^{+\infty} C_k \frac{1}{r} \int_{\Gamma} \frac{1}{y^{m-1}} e^{i(k+\frac{r}{\ell})y} dy, \end{aligned} \quad (45) $$

where the last integral is well defined since the path $\Gamma$ excludes the singularity at the origin from the integration. We know that an application of the residue theorem delivers [81]

$$ \int_{\Gamma} \frac{e^{iay}}{y^m} dy = \begin{cases} 0 & \text{if } a > 0, \\ -2\pi i^m \frac{a^{m-1}}{(m-1)!} & \text{if } a \le 0. \end{cases} \quad (46) $$

Therefore,

$$ \begin{aligned} Z_H(r) &= i \sum_{k=-\infty}^{+\infty} \frac{C_k}{r} 2\pi i^{m-1} \frac{(k+\frac{r}{\ell})^{m-2}}{(m-2)!} \mathbf{1} \left(-k-\frac{r}{\ell}\right) \\ &= \frac{2\pi i^m}{r} \sum_{h=-\infty}^{+\infty} C_{-h} \frac{(\frac{r}{\ell}-h)^{m-2}}{(m-2)!} \mathbf{1} \left(h-\frac{r}{\ell}\right), \end{aligned} \quad (47) $$

where $\mathbf{1}(x)$ represents the Heaviside step function, defined as $\mathbf{1}(x) = 1$ if $x \ge 0$, and $\mathbf{1}(x) = 0$ if $x < 0$. Finally,

$$ Z_H(r) = \frac{2\pi i^m}{r(m-2)!} \sum_{h=-\infty}^{+\infty} (C_h)^* \left(\frac{r}{\ell} - h\right)^{m-2} \mathbf{1}\left(h - \frac{r}{\ell}\right), \quad (48) $$

where we used the property stating that $C_{-h} = (C_h)^*$, which is valid for the Fourier coefficients of a real periodic function. The result obtained in Eq.(48) is exact for $\chi \in \mathbb{N}$, but it is based on the numerical computation of the coefficients $C_k$ (semi-analytic procedure). The limitation introduced by considering integer values for $\chi$ does not restrict the physical interpretation of the results. Moreover, this procedure can be easily generalized in order to consider arbitrary rational values for $\chi$ (of course, the function $P(y)$ remains periodic with $\chi \in \mathbb{Z}$). Furthermore, in next sections, we also discuss additional asymptotic results, which are not based on restrictions over the values of the parameter $\chi$.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Force-extension response (a), average number of unfolded units (b) and Helmholtz free energy (c) for a chain with Ising interactions under isometric (Helmholtz) conditions. The curves have been obtained with $\Delta E = 5K_BT$, $m=5$, $\chi = 2$, $\mu = 0$, and $\lambda = 0, \pm 0.5, \pm 1K_BT$. The black dashed curves correspond to the chain without interactions ($\lambda = 0$), the dark red (or dark gray) curves to ferromagnetic-like Ising interactions ($\lambda > 0$) and the orange (or light gray) curves to antiferromagnetic-like interactions ($\lambda < 0$).</figcaption></figure>

It is important to remark that our semi-analytic procedure, leading to Eq.(48) and based on the numerical implementation of Eq.(42), is very efficient for the determination of the Helmholtz partition function. Indeed, the direct numerical calculation of the original integral in Eq.(34), grounded on the knowledge of the Gibbs partition function given in Eq. (38), is a really hard – if not impossible – task since the integrand function is decreasing (as $1/y^{m-1}$) and oscillating for any $r$ in the whole interval between 0 and $m\chi\ell$. Since we need the quantity $\log Z_H(r)$ to analyze the system behavior, all the oscillations of the integrand function (also for large values of $|y|$) play an important role in defining the result. For this reason, our procedure leads to very accurate results, being based on the analytic determination of the integral over $\Gamma$ and on the numerical evaluation of the integrals over $(0, 2\pi)$ defined in Eq.(42), which are much more stable than the one defined in Eq.(34).

By mean of $Z_H$ given in Eq.(48), we can find the force-extension response through the expression

$$ \langle f \rangle = -K_B T \frac{\partial}{\partial r} \log Z_H(r) = \frac{\partial \mathcal{F}}{\partial r}, \quad (49) $$

and the average value of unfolded domains with the relation

$$ \left\langle \sum_{i=1}^{m} \frac{S_i + 1}{2} \right\rangle = -K_B T \frac{\partial}{\partial \Delta E} \log Z_H(r) = \frac{\partial \mathcal{F}}{\partial \Delta E}, \quad (50) $$

where we introduced the Helmholtz free energy of the system $\mathcal{F} = -K_B T \log Z_H$. An example of application can be found in Fig.4, where we show the force extension response, the average number of unfolded units and the Helmholtz free energy for a chain stretched under isometric conditions and with a variable Ising interaction coefficient. First of all, concerning the case with $\lambda = 0$, we observe that the force extension curve is composed of a number of peaks corresponding to the non-synchronized (sequential) unfolding of the units. This is confirmed by the step-wise curve representing the average number of unfolded units versus the chain extension. Each step corresponds to the unfolding of a unit induced by the increasing extension of the chain. This behavior agrees with previous theoretical and experimental results [54–61]. This scenario is modified by the introduction of the Ising interactions. If $\lambda > 0$, the unfolding of the first units requires a larger force peak since the units are favored to remain in the initial folded state. On the contrary, the unfolding of the last units requires a smaller force since most of the units are already unfolded and they prefer to be in the same state as the majority. This interpretation equally holds for the plot of the average value of unfolded units, where we can also note that the antiferromagnetic-like behavior may induce the unfolding of some units also without the applied force. The origin of the non-synchronized transitions can be highlighted in the Helmholtz free energy curves, characterized by a series of cusps able to induce the force peaks in the force
<figure><img src="image_6.png" /><figcaption>FIG. 5. Spinoidal phase diagram showing the end-points of the spinoidal intervals (with $\partial f/\partial r < 0$) versus the temperature of the system ($T > T_0$). Different Ising interaction coefficients $\lambda$ have been considered. Gray arrows indicate the increasing values of $\lambda$. The curves have been obtained with $\Delta E = 5K_B T_0$, $m = 5$, $\chi = 2$ and $\mu = 0$. The parameters $\ell$ and $T_0$ are arbitrary (typical values are $\ell = 1$nm and $T_0 = 300$K).</figcaption></figure>

extension curve. To conclude, the Ising interactions induce a specific cooperativity, which can be detected in the modification of the hierarchy of forces in the sawtooth-like response, as recently observed in force spectroscopy experiments of proteins (e.g., in Filamin A) [76].

A form of criticality can be noticed for the Helmholtz response of the bistable Ising chain. To do this, in the force-extension curves shown in Fig.4a, we can identify the spinoidal regions, characterized by a negative slope or, equivalently, by a negative differential stiffness. It means that, for each force peak observed in Fig.4a, we have a spinoidal interval with $\partial f/\partial r < 0$. It is interesting to study the evolution of these spinoidal regions in terms of the temperature. In general, we can say that the system is or not in a spinoidal phase depending on values of $r$ and $T$. We can therefore determine a sort of phase diagram, as shown in Fig.5, where the end-points of each spinoidal interval (on the extension axis) are shown versus the temperature. While the left end-point corresponds to the maximum of the force peak, the right end-point corresponds to the following minimum. These curves have been represented for different values of the interaction coefficient $\lambda$ to explore the effects of the Ising scheme on this critical behavior. Importantly, we can observe that each spinoidal interval disappears for a given temperature, which is a critical temperature for the system. Hence, for a given chain composed of $m$ units, there are $m$ different critical temperatures, one for each unfolding process. We remark that, for a system with-

out Ising interactions, the critical temperature is larger for the last unfolded units. This contrast among critical temperatures is further amplified for antiferromagnetic-like systems. On the other hand, a given intensity of ferromagnetic-like interactions is able to equilibrate the critical temperatures among the unfolding processes (see, e.g., the curves in Fig.5 corresponding to $\lambda = +1K_B T_0$). This point can be explained by observing that $\lambda > 0$ favors the simultaneous unfolding of the units, thus inducing a similar behavior of these ones. We can say that the ferromagnetic-like interactions induce a resistance to fluctuations within the system. It means that all unfolding processes do not lose their snap-through response for a given range of temperature. The bistability resistant to fluctuations is an important concept for micro-and nano-mechanical systems with non-convex elastic energy, where one attempts to sustain the bistability at possibly large temperatures [2]. We remark that the observation of a negative differential stiffness for subcritical temperatures and of a positive differential stiffness for supercritical temperatures can be interpreted by stating that the system behaves as a metamaterial [79, 80]. In a following section, we will also explore the effect of the intrinsic stiffness of the units on this critical behavior.

## V. HELMHOLTZ RESPONSE UNDER WEAK ISING INTERACTIONS

We investigate in more detail the particular case with weak Ising interaction, i.e. $|\lambda| \ll K_B T$, by considering both ferromagnetic-like and antiferromagnetic-like interactions. Under this condition, we will introduce an asymptotic development yielding a closed form expression of $Z_H$, not needing the numerical calculation of the coefficients $C_h$. To this aim, we observe that $Z_H$ can be written as in Eq.(34), and therefore we can develop $Z_G(iy\frac{K_B T}{\ell})$ in Taylor series with respect to the parameter $\lambda$. We easily find that

$$Z_G\left(iy\frac{K_B T}{\ell}\right) = \left[1 + \frac{\lambda}{K_B T}(m-1)\right](a+b)^m - \frac{4\lambda}{K_B T}(m-1)ab(a+b)^{m-2} + \mathcal{O}(\lambda^2). \quad (51)$$

The first order approximation in Eq.(51) (which is valid for $|\lambda| \ll K_B T$) can be easily integrated to obtain the Helmholtz partition function. Here, we defined

$$a = \frac{\sin y}{y}, \quad (52)$$

$$b = \chi\phi\frac{\sin\chi y}{y}. \quad (53)$$

In this section, we consider $\mu = 0$ to simplify the following calculations. We have to determine

$$Z_H(r) = -i \int_{\Gamma} Z_G \left( iy \frac{K_B T}{\ell} \right) \frac{y}{r} e^{i \frac{ry}{\ell}} dy \\ = \left[ 1 + \frac{\lambda}{K_B T} (m-1) \right] I_1 - \frac{4\lambda}{K_B T} (m-1) I_2,$$
<figure><img src="image_7.png" /><figcaption>FIG. 6. Force-extension response (a-c) and average number of unfolded units (b-d) for weak Ising interactions with ferromagnetic-like behavior (a-b) and antiferromagnetic-like behavior (c-d). In each panel, the response without interactions ($\lambda = 0$, black dashed curves) is shown together with the exact result (orange or light gray curves) obtained through Eq.(48) and the first order approximation (dark red or dark gray curves) given in Eq.(61). The curves have been obtained with $\Delta E = 5K_B T$, $m = 5$, $\chi = 2$, $\mu = 0$ and $\lambda = \pm 0.3K_B T$.</figcaption></figure>

where we used the approximation of $Z_G$ given in Eq.(51). We firstly calculate the quantity $I_1$, as follows

$$ 
\begin{aligned}
I_1 &= -i \int_{\Gamma} (a+b)^m \frac{y}{r} e^{i \frac{ry}{\ell}} dy \\
&= -i \int_{\Gamma} \sum_{k=0}^{m} \binom{m}{k} \left(\frac{\sin y}{y}\right)^{m-k} \left(\chi\phi \frac{\sin \chi y}{y}\right)^k \frac{y}{r} e^{i \frac{ry}{\ell}} dy.
\end{aligned} 
\quad (55) $$

Since

$$ 
\begin{aligned}
\sin^{m-k} y &= \frac{1}{(2i)^{m-k}} \sum_{p=0}^{m-k} \binom{m-k}{p} (-1)^p e^{iy(m-k-2p)}, \\
\sin^k(\chi y) &= \frac{1}{(2i)^k} \sum_{q=0}^{k} \binom{k}{q} (-1)^q e^{i\chi y(k-2q)},
\end{aligned} 
\quad (56) $$

we find

$$ 
\begin{aligned}
I_1 &= -i \int_{\Gamma} \frac{1}{y^m} \sum_{k=0}^{m} \sum_{p=0}^{m-k} \sum_{q=0}^{k} \binom{m}{k} \binom{m-k}{p} \binom{k}{q} \frac{1}{(2i)^m} \\
&\quad \times (-1)^{p+q} (\chi\phi)^k e^{iy(m-k-2p+\chi k-2q\chi)} \frac{y}{r} e^{i\frac{ry}{\ell}} \\
&= \frac{1}{2^{m+i(m+1)}} \frac{1}{r} \sum_{k=0}^{m} \sum_{p=0}^{m-k} \sum_{q=0}^{k} \binom{m}{k} \binom{m-k}{p} \binom{k}{q} \\
&\quad \times (-1)^{p+q} (\chi\phi)^k \int_{\Gamma} \frac{e^{iy(m-k-2p+\chi k-2q\chi)}}{y^{m-1}} dy.
\end{aligned} 
\quad (57) $$

Now, we use the result in Eq.(46) and we obtain

$$ 
\begin{aligned}
I_1 &= \frac{\pi}{2^{m-1}} \frac{1}{r(m-2)!} \sum_{k=0}^{m} \sum_{p=0}^{m-k} \sum_{q=0}^{k} \binom{m}{k} \binom{m-k}{p} \binom{k}{q} \\
&\quad \times (-1)^{p+q} (\chi\phi)^k (-\Lambda)^{m-2} \mathbf{1}(\Lambda),
\end{aligned} 
\quad (58) $$
where $-\Lambda = m - k - 2p + \chi k - 2\chi q + \frac{r}{\ell}$. Similarly, we calculate the integral $I_2$ given by

$$ 
\begin{aligned}
I_2 &= -i \int_{\Gamma} ab(a+b)^{m-2} \frac{y}{r} e^{i \frac{ry}{\ell}} dy \\
&= -i \int_{\Gamma} \sum_{k=0}^{m-2} \binom{m-2}{k} b^{k+1} a^{m-1-k} \frac{y}{r} e^{i \frac{ry}{\ell}} dy \\
&= -i \sum_{k=0}^{m-2} \sum_{q=0}^{k+1} \sum_{p=0}^{m-1-k} \binom{m-2}{k} \binom{k+1}{q} \binom{m-1-k}{p} \\
&\quad \times \int_{\Gamma} \frac{y}{r} e^{i \frac{ry}{\ell}} (-1)^{p+q} \frac{1}{(2iy)^m} (\chi\phi)^{k+1} e^{i\chi y(k+1-q)} \\
&\quad \times e^{-i\chi yq} e^{iy(m-1-k-p)} e^{-iyp} dy,
\end{aligned} 
\quad (59) $$

and, by using again the integral in Eq.(46), we eventually obtain

$$ 
\begin{aligned}
I_2 = & \frac{\pi}{2^{m-1}} \frac{1}{r(m-2)!} \sum_{k=0}^{m-2} \sum_{q=0}^{k+1} \sum_{p=0}^{m-1-k} \binom{m-2}{k} \binom{k+1}{q} \\
& \times \binom{m-1-k}{p} (-1)^{p+q} (\chi\phi)^{k+1} (-\Lambda_0)^{m-2} \mathbf{1}(\Lambda_0),
\end{aligned} 
\quad (60) $$

where $-\Lambda_0 = m - 1 - k - 2p - 2\chi q + \chi k + \chi + \frac{r}{\ell}$. Finally, the partition function reads

$$ 
\begin{aligned}
Z_H(r) = & \frac{\pi}{2^{m-1}r(m-2)!} \left[ 1 + \frac{\lambda}{K_B T} (m-1) \right] \\
& \times \sum_{k=0}^{m} \sum_{p=0}^{m-k} \sum_{q=0}^{k} \binom{m}{k} \binom{m-k}{p} \binom{k}{q} \\
& \times (-1)^{p+q} (\chi\phi)^k (-\Lambda)^{m-2} \mathbf{1}(\Lambda) \\
& - \frac{\pi}{2^{m-1}r(m-2)!} \frac{4\lambda}{K_B T} (m-1) \\
& \times \sum_{k=0}^{m-2} \sum_{q=0}^{k+1} \sum_{p=0}^{m-1-k} \binom{m-2}{k} \binom{k+1}{q} \binom{m-1-k}{p} \\
& \times (-1)^{p+q} (\chi\phi)^{k+1} (-\Lambda_0)^{m-2} \mathbf{1}(\Lambda_0).
\end{aligned} 
\quad (61) $$

This is the final form of the Helmholtz partition function, calculated under the hypothesis of weak Ising interaction. In this expression, there are no limitations concerning the parameter $\chi$. We note that, for $\lambda = 0$, we obtain the partition function of the system without Ising interactions, discussed in recent literature [81]. In Fig.6, one can find some results for $\lambda = 0, \lambda = +0.3K_B T$ and $\lambda = -0.3K_B T$. We plotted the force-extension curves and the average number of unfolded units for both ferromagnetic-like and antiferromagnetic-like interactions. In each case, we compared the approximated result stated in Eq.(61) (dark red or dark gray curves), the exact result given in Eq.(48) (orange or light gray curves), and the response without Ising interactions (black dashed curves). We note a good agreement between approximated and exact results both for $\lambda < 0$ and $\lambda > 0$. Concerning the interpretation of the curves, the discussion

reported at the end of Section IV remains valid for all results of Fig.6.

## VI. HELMHOLTZ RESPONSE UNDER STRONG ISING FERROMAGNETIC-LIKE INTERACTIONS

We introduce here an asymptotic development concerning the case of a strong Ising ferromagnetic-like interaction. If $\lambda \to +\infty$, then $x = e^{-\kappa_B T} \to 0$, and it is not difficult to obtain the asymptotic expression for the analytic continuation of $Z_G$. The result is

$$ 
\begin{aligned}
Z_G \left( iy \frac{K_B T}{\ell} \right) = & \exp \left[ \frac{\lambda(m-1)}{K_B T} \right] \\
& \times \left\{ \left[ \frac{\sin(y)}{y} \right]^m + \left[ \chi\phi \frac{\sin(\chi y)}{y} \right]^m \right\}.
\end{aligned} 
\quad (62) $$

Therefore, we easily determine the Helmholtz partition function, as follows

$$ 
\begin{aligned}
Z_H(r) = & -i \int_{\Gamma} Z_G \left( iy \frac{K_B T}{\ell} \right) \frac{y}{r} e^{i \frac{ry}{\ell}} dy \\
& = -i \int_{\Gamma} \exp \left[ \frac{\lambda(m-1)}{K_B T} \right] \frac{1}{(2i)^m} \sum_{k=0}^{m} \binom{m}{k} (-1)^k \frac{1}{r} \\
& \times \left[ e^{iy(m-2k)} + (\chi\phi)^m e^{i\chi y(m-2k)} \right] \frac{1}{y^{m-1}} e^{i \frac{ry}{\ell}} dy.
\end{aligned} 
\quad (63) $$

To conclude the calculation, we can use the integral in Eq.(46) and we get the final result

$$ 
\begin{aligned}
Z_H(r) = & \frac{\pi}{2^{m-1}(m-2)!r} \exp \left[ \frac{\lambda(m-1)}{K_B T} \right] \sum_{k=0}^{m} \binom{m}{k} (-1)^k \\
& \times \left[ \left( m - 2k + \frac{r}{\ell} \right)^{m-2} \mathbf{1} \left( 2k - m - \frac{r}{\ell} \right) + (\chi\phi)^m \right. \\
& \left. \times \left( \chi m - 2\chi k + \frac{r}{\ell} \right)^{m-2} \mathbf{1} \left( 2\chi k - \chi m - \frac{r}{\ell} \right) \right],
\end{aligned} 
\quad (64) $$

which is valid for strong ferromagnetic-like Ising interactions. An application of this expression is shown in Fig.7. In particular, we compare the approximated result in Eq.(64) (yellow or light gray curve) with the exact response obtained from Eq.(48) for $\lambda = 0.5, 1, 1.5, \dots, 7K_B T$ (red or gray curves) and with the response without Ising interactions (black dashed curve). In Fig.7, one can find the force-extention curves, the average number of unfolded units and the Helmholtz free energy. It is interesting to discuss the evolution of the overall behavior of the system with an increasing interaction coefficient. Indeed, as $\lambda$ is increased, the units are progressively favored to be in the same state, and therefore there is an increasing average number of units which unfolds at $r = m\ell$. It means that the number of unfolding processes at $r = m\ell$ is a growing function of the Ising coefficient $\lambda$, going from 1 with $\lambda = 0$ to $m$ with $\lambda$ approaching infinity. This can be seen in Fig.7b, where this process is represented by
<figure><img src="image_8.png" /><figcaption>FIG. 7. Force-extension response (a), average number of unfolded units (b) and Helmholtz free energy (c) for a chain with strong ferromagnetic-like interactions. In each panel, the response without interactions ($\lambda = 0$, black dashed curves) is shown together with the exact results for an increasing ferromagnetic interaction ($\lambda = 0.5jK_BT \forall j = 1, \dots, 14$, red or gray curves) and the asymptotic expansion ($\lambda \to \infty$, yellow or light gray curves) given in Eq.(64). The curves have been obtained with $\Delta E = 5K_BT$, $m = 5$, $\chi = 2$ and $\mu = 0$.</figcaption></figure>

the series of red curves (or gray) with increasing $\lambda$, and it ends with the yellow (or light gray) curve obtained through Eq.(64). The latter means that all units unfold at the same time at $r = m\ell$ when $\lambda \to \infty$, and this behavior is perfectly caught by the asymptotic development. Accordingly, the peaks in the force-extension curve are strongly modified by increasing $\lambda$: while the first peak becomes more and more pronounced, the others are progressively reduced, as shown in Fig.7a. As a matter of fact, the first peak corresponds to the simultaneous unfolding of the units when $\lambda$ is very large. Hence, in the limiting case of $\lambda \to \infty$, the force-extension curve is composed of only one peak (yellow or light gray curve), as one can see in Fig.7a. Of course, the origin of peaks modification in the force-extension curve and of the steps structure in the average number of unfolded units can be observed in the plot of the Helmholtz free energy, shown in Fig.7c. Here, we can see the evolution of the typical cusps with the increasing Ising coefficient. As an example, the collapse of all the force peaks into a single unfolding event explains the tandem repeats behavior in red cell spectrin, where two units unfold simultaneously because of a strong cooperativity [75].

## VII. HELMHOLTZ RESPONSE UNDER STRONG ISING ANTIFERROMAGNETIC-LIKE INTERACTIONS

We discuss here the development of the theory under strong Ising antiferromagnetic-like interactions. As before, we can develop the analytic continuation of the Gibbs partition function in a power series of $\lambda \to -\infty$, i.e. for $x = e^{-K_BT}/\lambda \to +\infty$. The result can be eventually obtained as

$$ 
\begin{aligned}
Z_G \left( iy \frac{K_B T}{\ell} \right) &= x^{m-1} \left\{ P^{\frac{m}{2}} [1 + (-1)^m] \right. \\
&\qquad \left. + \frac{1}{2} P^{\frac{m-1}{2}} S [1 - (-1)^m] \right\} \quad (65) \\
&= \begin{cases} x^{m-1} P^{\frac{m}{2}} & \text{if } m \text{ is even,} \\ x^{m-1} S P^{\frac{m-1}{2}} & \text{if } m \text{ is odd,} \end{cases}
\end{aligned}
$$

where

$$ S = \frac{\sin(y)}{y} + \chi\phi \frac{\sin(\chi y)}{y}, \quad (66) $$

$$ P = \frac{\sin(y)}{y} \chi\phi \frac{\sin(\chi y)}{y}. \quad (67) $$

We will develop the asymptotic theory for both the cases with $m$ odd and even. We first elaborate the Helmholtz partition function for $m$ odd

$$ 
\begin{aligned}
Z_H(r) &= -i \int_{\Gamma} Z_G \left( iy \frac{K_B T}{\ell} \right) \frac{y}{r} e^{i \frac{r y}{\ell}} dy \quad (68) \\
&= -i \frac{x^{m-1}}{r} \left( \chi\phi \right)^{\frac{m-1}{2}} \int_{\Gamma} \left[ \sin^{\frac{m+1}{2}}(y) \sin^{\frac{m-1}{2}}(\chi y) \right. \\
&\qquad \left. + \chi\phi \sin^{\frac{m-1}{2}}(y) \sin^{\frac{m+1}{2}}(\chi y) \right] \frac{1}{y^{m-1}} e^{i \frac{r y}{\ell}} dy.
\end{aligned}
$$
Here, we use the relation

$$ \sin^k y = \frac{1}{(2i)^k} \sum_{p=0}^{k} \binom{k}{p} (-1)^p e^{iy(k-2p)}, \quad (69) $$

and we get

$$ \begin{aligned} Z_H(r) = & -i \frac{x^{m-1}}{r} (\chi\phi)^{\frac{m-1}{2}} \int_{\Gamma} \frac{1}{y^{m-1}} e^{i \frac{r y}{\ell}} \\ & \times \left[ \frac{1}{(2i)^{\frac{m+1}{2}}} \sum_{p=0}^{\frac{m+1}{2}} \left( \frac{m+1}{p} \right) (-1)^p e^{iy \left( \frac{m-1}{2} - 2p \right)} \right. \\ & \left. \times \frac{1}{(2i)^{\frac{m-1}{2}}} \sum_{q=0}^{\frac{m-1}{2}} \left( \frac{m-1}{q} \right) (-1)^q e^{iy \chi \left( \frac{m-1}{2} - 2q \right)} \right. \\ & \left. + \chi\phi \frac{1}{(2i)^{\frac{m-1}{2}}} \sum_{p=0}^{\frac{m-1}{2}} \left( \frac{m-1}{q} \right) (-1)^p e^{iy \chi \left( \frac{m-1}{2} - 2p \right)} \right] dy. \end{aligned} \quad (70) $$

Then, straightforward calculations deliver

$$ \begin{aligned} Z_H(r) = & -\frac{ix^{m-1}}{r(2i)^m} (\chi\phi)^{\frac{m-1}{2}} \sum_{p=0}^{\frac{m+1}{2}} \sum_{q=0}^{\frac{m-1}{2}} \left(\frac{m+1}{p}\right) \left(\frac{m-1}{q}\right) \\ & \times (-1)^{p+q} \left[ \int_{\Gamma} \frac{e^{-iy\Lambda_a}}{y^{m-1}} dy + \chi\phi \int_{\Gamma} \frac{e^{-iy\Lambda_b}}{y^{m-1}} dy \right], \end{aligned} \quad (71) $$

where

$$ \Lambda_a = 2p - \frac{m+1}{2} - \chi \frac{m-1}{2} + 2\chi q - \frac{r}{\ell}, \quad (72) $$

$$ \Lambda_b = 2q - \frac{m-1}{2} - \chi \frac{m+1}{2} + 2\chi p - \frac{r}{\ell}. \quad (73) $$

Hence, by using Eq.(46), we get the final result

$$ \begin{aligned} Z_H(r) = & \frac{2\pi}{2^m} \frac{\exp\left[-\frac{\lambda}{K_B T}(m-1)\right]}{(m-2)!} \frac{1}{r} (\chi\phi)^{\frac{m-1}{2}} \\ & \times \sum_{p=0}^{\frac{m+1}{2}} \sum_{q=0}^{\frac{m-1}{2}} \left(\frac{m+1}{p}\right) \left(\frac{m-1}{q}\right) (-1)^{p+q} \\ & \times \left[(-\Lambda_a)^{m-2} \mathbf{1}(\Lambda_a) + \chi\phi (-\Lambda_b)^{m-2} \mathbf{1}(\Lambda_b)\right], \end{aligned} \quad (74) $$

which is valid for *m* odd.

We calculate now the same quantity for *m* even

$$ \begin{aligned} Z_H(r) = & -i \int_{\Gamma} Z_G \left( i y \frac{K_B T}{\ell} \right) \frac{y}{r} e^{i \frac{r y}{\ell}} dy \\ = & -i \int_{\Gamma} x^{m-1} \left[ \frac{\sin(y)}{y} \chi\phi \frac{\sin(\chi y)}{y} \right]^{\frac{m}{2}} \frac{y}{r} e^{i \frac{r y}{\ell}} dy \\ = & -\frac{ix^{m-1}}{r(2i)^m} (\chi\phi)^{\frac{m}{2}} \sum_{p=0}^{\frac{m}{2}} \sum_{q=0}^{\frac{m}{2}} \left(\frac{m}{p}\right) \left(\frac{m}{q}\right) (-1)^{p+q} \\ & \times \int_{\Gamma} \frac{1}{y^{m-1}} \exp[-iy\Lambda_c] dy \end{aligned} \quad (75) $$

where $\Lambda_c = 2p - \frac{m}{2} + 2\chi q - \chi \frac{m}{2} - \frac{r}{\ell}$. Finally, by using again Eq.(46), we get

$$ \begin{aligned} Z_H(r) = & \frac{2\pi}{2^m} \frac{\exp\left[-\frac{\lambda}{K_B T}(m-1)\right]}{(m-2)!} \frac{1}{r} (\chi\phi)^{\frac{m}{2}} \\ & \times \sum_{p=0}^{\frac{m}{2}} \sum_{q=0}^{\frac{m}{2}} \left(\frac{m}{p}\right) \left(\frac{m}{q}\right) (-1)^{p+q} (-\Lambda_c)^{m-2} \mathbf{1}(\Lambda_c), \end{aligned} \quad (76) $$

which is valid for *m* even. The solutions given in Eqs.(74) and (76) represent the most important result of this section and must be discussed as follows.

An example of application of Eq.(74) for *m* odd can be found in Fig.8, where force-extension curves, average number of unfolded units and Helmholtz free energy are represented for several values of $\lambda$. When we consider an increasing value of $|\lambda|$ ($\lambda < 0$), we observe that the first force peaks tend to disappear while the last ones become more and more pronounced (Fig.8a). This is coherent with the assumption that, in an antiferromagnetic-like system, the favored states are alternatively folded and unfolded. Accordingly, with an increasing value of $|\lambda|$ ($\lambda < 0$), we have an increasing number of unfolded units in the initial configuration with $r = 0$. Clearly, the maximum value of this number of unfolded units with $r = 0$ is $(m-1)/2$ and not $(m+1)/2$ since the unfolded units are costly from the energetic point of view and the system chooses the configuration with the smallest number of unfolded units between $(m-1)/2$ and $(m+1)/2$. This overall interpretation of Fig.8 must be improved to better understand the yellow (or light gray) curves corresponding to $\lambda \to \infty$. So, if we look at Eq.(74), we note that $Z_H(r)$ is different from 0 when $\Lambda_a > 0$ or $\Lambda_b > 0$. It means

$$ 2p - \frac{m+1}{2} - \chi \frac{m-1}{2} + 2\chi q - \frac{r}{\ell} > 0, \quad (77) $$

or

$$ 2q - \frac{m-1}{2} - \chi \frac{m+1}{2} + 2\chi p - \frac{r}{\ell} > 0. \quad (78) $$

Now, the maximum values of *p* and *q* are $(m+1)/2$ and $(m-1)/2$, respectively. Hence, we have

$$ \frac{r}{\ell} < 2 \frac{m+1}{2} - \frac{m+1}{2} - \chi \frac{m-1}{2} + 2\chi \frac{m-1}{2} \quad (79) $$

or

$$ \frac{r}{\ell} < 2 \frac{m-1}{2} - \frac{m-1}{2} - \chi \frac{m+1}{2} + 2\chi \frac{m+1}{2}. \quad (80) $$

Equivalently,

$$ \frac{r}{\ell} < \frac{m+1}{2} + \chi \frac{m-1}{2} < m\chi \quad (81) $$

or

$$ \frac{r}{\ell} < \frac{m-1}{2} + \chi \frac{m+1}{2} < m\chi, \quad (82) $$

where $m\chi$ is the upper limit of $\frac{r}{\ell}$, attained when all elements are unfolded. It means that the total length
<figure><img src="image_9.png" /></figure>

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

<figure><img src="image_11.png" /><figcaption>FIG. 8. Force-extension response (a), average number of unfolded units (b) and Helmholtz free energy (c) for a chain with strong antiferromagnetic-like interactions and an odd number of units. In each panel, the response without interactions ($\lambda = 0$, black dashed curve) is shown together with the exact results for an increasing antiferromagnetic interaction ($\lambda = -0.25jK_BT \quad \forall j = 1, ..., 6$, red or gray curves) and the asymptotic expansion ($\lambda \to -\infty$, yellow or light gray curve) given in Eq.(74). The curves have been obtained with $\Delta E = 5K_BT$, $m = 5$, $\chi = 2$ and $\mu = 0$.</figcaption></figure>

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

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

<figure><img src="image_14.png" /><figcaption>FIG. 9. Force-extension response (a), average number of unfolded units (b) and Helmholtz free energy (c) for a chain with strong antiferromagnetic-like interactions and an even number of units. In each panel, the response without interactions ($\lambda = 0$, black dashed curve) is shown together with the exact results for an increasing antiferromagnetic interaction ($\lambda = -0.25jK_BT \quad \forall j = 1, ..., 6$, red or gray curves) and the asymptotic expansion ($\lambda \to -\infty$, yellow or light gray curve) given in Eq.(76). The curves have been obtained with $\Delta E = 5K_BT$, $m = 4$, $\chi = 2$ and $\mu = 0$.</figcaption></figure>
of the unit can not exceed $\frac{m-1}{2} + \chi \frac{m+1}{2}$ (which is the highest value between Eqs.(81) and (82)). The value $\frac{m+1}{2} + \chi \frac{m-1}{2}$ corresponds to $\frac{m+1}{2}$ folded domains and $\frac{m-1}{2}$ unfolded domains: this is the starting configuration for the strong antiferromagnetic system, which can be represented as $\vec{S} = \{\uparrow\downarrow\uparrow\downarrow\uparrow\}$ for $m=5$, where $\{\uparrow\}$ represents a folded unit and $\{\downarrow\}$ an unfolded one. This configuration is stable with $r=0$. When we apply a sufficient extension, the unfolded domains become more stable, and we have an inversion in the alternating disposition leading to $\vec{S} = \{\downarrow\uparrow\downarrow\uparrow\downarrow\}$. This transition can be observed through the single peak in the force-extension yellow (or light gray) curve (Fig.8a) and in the shift from 2 to 3 of the average number of unfolded units (yellow or light gray curve in Fig.8b). As a matter of fact, in this case, we can not unfold all elements because of the hypothesis of large (ideally infinite) antiferromagnetism ($\lambda \to -\infty$). For this reason, $Z_H(r)$ is defined for $r < \frac{m-1}{2}\ell + \frac{m+1}{2}\chi\ell$ and not for $r < m\chi\ell$, as in previous cases.

A similar discussion holds for the case with $m$ even. An example of application of Eq.(76) is shown in Fig.9. While the red (or gray) curves with an increasing value of $|\lambda|$ ($\lambda < 0$) are similar to those shown in Fig.8, the yellow (or light gray) curves representing the asymptotic behavior for $\lambda \to -\infty$ are different and require some comments. In this case, $Z_H(r)$ is different from 0 if $\Lambda_c > 0$, or if $2p - \frac{m}{2} + 2\chi q - \chi \frac{m}{2} - \frac{r}{\ell} > 0$. Now, $p$ and $q$ assume the same maximum value $\frac{m}{2}$ and therefore we have $\frac{r}{\ell} < 2\frac{m}{2} - \frac{m}{2} + 2\chi \frac{m}{2} - \chi \frac{m}{2}$ or equivalently, $\frac{r}{\ell} < \frac{m}{2} + \chi \frac{m}{2}$. It is evident that $(\frac{m}{2} + \chi \frac{m}{2})\ell$ is the length of the chain with $\frac{m}{2}$ folded domains and $\frac{m}{2}$ unfolded domains. So, if $\lambda \to -\infty$, the number of unfolded domains is a constant equal to $\frac{m}{2}$ and $Z_H(r)$ is defined only for $r < (\frac{m}{2} + \chi \frac{m}{2})\ell$. In this case, we have no inversion of the configuration, which corresponds to $\vec{S} = \{\uparrow\downarrow\uparrow\downarrow\uparrow\}$ or to $\vec{S} = \{\uparrow\downarrow\uparrow\downarrow\downarrow\}$ for $\lambda \to -\infty$ and $m=4$. These two configuration are indeed completely equivalent from the energetic point of view and they are therefore indistinguishable. As a conclusion, if $\lambda \to -\infty$, the force extension curve is a monotonically increasing function of $r$ (without peaks), the average number of unfolded units is a constant equal to $m/2$ and the Helmholtz free energy is without cusps (yellow or light gray curves in Fig.9).

VIII. CHAIN WITH ISING INTERACTIONS AND EXTENSIBLE UNITS

We consider now a chain of bistable units characterized by a finite elastic constant. We start the analysis by considering the Gibbs ensemble defined through the extended Hamiltonian given in Eq.(2). Here, for the sake of simplicity, we suppose that the folded and unfolded basins of the potential energy shown in Fig.1 exhibit the same elastic constant $k(-1) = k(+1) = \kappa$. In this case, it is not difficult to prove that Eq.(7) of Section III can

be substituted by the approximated result

$$I = \exp \left( \frac{f^2}{2\kappa K_B T} \right) \prod_{i=1}^{m} \frac{\sinh \left( \frac{f \ell_0(S_i)}{K_B T} \right)}{\frac{f \ell_0(S_i)}{K_B T} \ell^2} \frac{\ell_0^2(S_i)}{\ell^2}, \quad (83)$$

which can be easily proved by calculating the exact solution of Eq.(5) and by considering its approximation for finite relatively high values of $\kappa$. The details of this calculation can be found in recent literature [82]. Typical values of $\kappa$ for real macromolecules are compatible with such approximations [89]. Coherently with Eq.(83), the Gibbs partition function can be obtained as

$$Z_G(f) = \frac{x}{1+x^2} \exp\left(\frac{m\alpha\xi^2}{2}\right) \quad (84) \\ \times \left[ \lambda_1^m \left(1 + x^2 \frac{\lambda_1 + \lambda_2}{\lambda_1 - \lambda_2}\right) + \lambda_2^m \left(1 - x^2 \frac{\lambda_1 + \lambda_2}{\lambda_1 - \lambda_2}\right) \right],$$

where, with respect to Eq.(24) of Section III, we added an exponential term, which is quadratic in the normalized force $\xi = \frac{f\ell}{K_B T}$ and defined through a coefficient $\alpha = \frac{K_B T}{\kappa \ell^2}$. This latter represents the ratio between the thermal energy and the elastic one, thus measuring the compromise between the enthalpic contribution and fluctuations. The parameters $x$, $\lambda_1$ and $\lambda_2$ remain defined as in Section III.

Concerning the Helmholtz ensemble, the system is described by the Hamiltonian in Eq.(31), where, as before, we assume that $k(-1) = k(+1) = \kappa$. The calculation of the corresponding partition function is based on the Laplace transform relation between the Gibbs and the Helmholtz ensembles, summed up in Eq.(34). If we consider integer values for the parameter $\chi$ (ratio between unfolded and folded length), the analytic continuation of the Gibbs partition function assumes the form

$$Z_G \left( iy \frac{K_B T}{\ell} \right) = \frac{1}{y^m} \exp \left( - \frac{m \alpha y^2}{2} \right) P(y), \quad (85)$$

where $P(y)$ is the periodic function defined in Eq.(40) and described by the Fourier coefficients given in Eq.(42). The Helmholtz partition function can be finally determined as follows

$$\begin{align*}
Z_H(r) &= -i \int_{\Gamma} Z_G \left( iy \frac{K_B T}{\ell} \right) \frac{y}{r} e^{\frac{yr}{\ell}} dy \\
&= -i \int_{\Gamma} \frac{P(y)}{y^m} \frac{y}{r} e^{-\frac{m\alpha y^2}{2}} e^{\frac{iry}{\ell}} dy \\
&= -i \sum_{k=-\infty}^{+\infty} C_k \frac{1}{r} \int_{\Gamma} \frac{1}{y^{m-1}} e^{-\frac{m\alpha y^2}{2}} e^{i(k+r/\ell)y} dy \\
&= -i \sum_{k=-\infty}^{+\infty} C_k \frac{1}{r} J_{m-1} \left( k + \frac{r}{\ell}, \frac{m\alpha}{2} \right) \quad (86)
\end{align*}$$

where we used the Fourier development of the function $P(y)$ in order to perform the calculation. To complete
<figure><img src="image_15.png" /><figcaption>FIG. 10. Force-extension response for a chain of interacting units with finite intrinsic stiffness $κ$. We used the values $κℓ^2 = 100K_B T$ (a), and $κℓ^2 = 200K_B T$ (b). In each panel, the response without interactions ($λ = 0$, black curves) is shown together with the results with $λ = +1K_B T$ (dark red or dark gray curves) and $λ = -1K_B T$ (orange or light gray curves). The Gibbs and Helmholtz responses correspond to dashed and solid lines, respectively. The curves have been obtained with $ΔE = 5K_B T$, $m = 5$, $χ = 2$ and $μ = 0$.</figcaption></figure>

the task, we have to calculate a sequence of integral of the form

$$J_N(a, b) = \int_{\Gamma} e^{iay} e^{-by^2} \frac{dy}{y^N}, \quad (87)$$

where $a \in \mathbb{R}$, $b > 0$, $N \in \mathbb{N}$ and the path $\Gamma$ is given in Fig.3. An application of the complex variable method allows us to obtain the closed form expression for this integral, as follows [82]

$$J_N(a, b) = \sqrt{\pi} \left\{ b^{\frac{N-1}{2}} e^{-\frac{a^2}{4b}} (2i)^N H_{-N} \left( -\frac{a}{2\sqrt{b}} \right) + \mathbf{1}(-a) \sum_{h=0}^{N-1} \frac{(ia)^{N-1-h} (-4b)^{h/2}}{ih! (N-1-h)!} \Gamma \left( \frac{h+1}{2} \right) [1+(-1)^h] \right\}, \quad (88)$$

where $\mathbf{1}(x)$ represents the Heaviside step function, and $H_{-N}(z)$ are the generalization of the Hermite polynomials, obtained by considering negative indices. They can be defined as [82, 90]

$$H_{-N}(z) = \frac{1}{(2i)^N \sqrt{\pi}} \int_{-\infty}^{+\infty} \frac{e^{-\eta^2} d\eta}{(\eta - iz)^N}, \quad (89)$$

and they can be obtained recursively through the following formula [82, 90]

$$H_{-(N+1)}(z) = \frac{1}{2N} H_{-(N-1)}(z) - \frac{z}{N} H_{-N}(z), \quad (90)$$

initialized with $H_0(z) = 1 \forall z \in \mathbb{R}$ and $H_{-1}(z) = \frac{\sqrt{\pi}}{2} e^{z^2} [1 - \operatorname{erf}(z)]$ for $z > 0$ and $H_{-1}(z)$ odd, i.e. $\dot{H}_{-1}(-z) = -H_{-1}(z)$ [82, 90]. This recursive law allows us to affirm that $H_{-N}(z)$ is odd if $N$ is odd and

that $H_{-N}(z)$ is even if $N$ is even. We underline that the Hermite elements with $N < 0$ are not polynomials, contrarily to the classical Hermite polynomials, defined with $N > 0$. Many other interesting properties of this sequence of functions can be found in Ref. [90].

An application of the Gibbs and Helmholtz partition functions, stated in Eqs. (84) and (86), respectively, is presented in Fig.10, where we show the force-extension curves for two values of the constant $κ$, and for three values of the coefficient $λ$. Firstly, we note that the constant slope of the final part of the force-extension curves represents the finite effective stiffness of the chain, after the unfolding processes. Moreover, it is interesting to remark that the softer systems exhibit a sensibly reduced force peaks in the Helmholtz response. This point can be also noticed by drawing a comparison between Fig.4a, obtained for $κ \to \infty$, and Fig.10, corresponding to finite values of $κ$. A similar phenomenon can be also observed in the phase diagram showing the critical behavior of the spinoidal response of the system. Indeed, we plotted in Fig.11 four phase diagrams corresponding to four different values of the elastic constant. We observe that the critical temperature of the unfolding processes is an increasing function of $κ$, similarly to the previously discussed force peaks of the Helmholtz response. Besides, as already seen in Fig.5 concerning the case with $κ \to \infty$, antiferromagnetic-like interactions amplify the dissimilarity among the critical temperatures, while ferromagnetic-like interactions reduce this contrast, eventually producing a more uniform response of the unfolding processes.
<figure><img src="image_16.png" /><figcaption>FIG. 11. Spinoidal phase diagrams obtained for four different values of the intrinsic stiffness $\kappa$ of the units (defined by $\kappa\ell^2 = 100, 133, 166, 200K_BT_0$). They show the end-points of the spinoidal intervals versus the temperature of the system ($T > T_0$). The curves have been obtained with $\Delta E = 5K_BT_0$, $m = 5$, $\chi = 2$, $\mu = 0$ and five different values of the Ising interaction coefficient $\lambda = 0, \pm 0.5, \pm 1K_BT_0$. Gray arrows indicate increasing values of $\lambda$. The parameters $\ell$ and $T_0$ are arbitrary (typical values are $\ell = 1$nm and $T_0 = 300$K).</figcaption></figure>

## IX. CONCLUSIONS

We investigated the properties of a chain of two-state units coupled through an Ising interaction scheme, providing a paradigmatic description of the effects of bistability and cooperativity in biological and artificial micro- and nano-systems. Accordingly, we studied our model by means of the statistical mechanics of small systems, i.e. far from the thermodynamic limit. It means that, for a limited number *m* of units of the chain, the Gibbs and Helmholtz statistical ensembles are not equivalent and we ultimately obtain two different isotensional and isometric responses, well recognized, e.g., in force-spectroscopy experiments.

Some of the most interesting findings of this paper concern the influence of the cooperativity, measured by the Ising coefficient $\lambda$, on the mechanical behavior and on the configurational features of the system. In particular, we analyzed the force-extension curve under isotensional conditions, obtaining a sharper or smoother transition depending on $\lambda$, and under isometric conditions, getting a variable hierarchy of force peaks as function of the cooperativity. Also, the unfolding processes of the units have been characterized by plotting the number of unfolded units versus the mechanical quantity ($f$ or $r$) inducing the chain stretching. This point allows the interpretation of the unfolding processes as synchronized or simultaneous under isotensional conditions, and as non-synchronized or sequential under isometric conditions. This result underlines the convenience of the spin variables to investigate the configurational properties of the system.

From the methodological point of view, we underline that the spin variables approach is useful to elaborate semi-analytic or closed-form expressions for the relevant observables. More specifically, to solve the problem within the Gibbs ensemble, we coupled this spin variables approach with the classic transfer matrix technique to take account of the interactions. On the other hand, to overcome the complexity induced by the isometric conditions, we made use of the Laplace transform between Gibbs and Helmholtz partition functions. This expedient consents to study the Helmholtz ensemble response by solving a given integral, based on the analytic continuation of the Gibbs partition function on the imaginary axis. In addition, we proposed explicit asymptotic results describing the behavior of the system under weak and strong Ising interactions (for both the ferromagnetic-like and the antiferromagnetic-like schemes).

To give a complete picture of the equilibrium behavior of the system, we also investigated a form of criticality exhibited by the system. In particular, our analysis highlights the critical behavior of the spinoidal regions, characterizing the part of the isometric response showing a negative differential elastic stiffness. We prove that each unfolding process exhibits a critical temperature defined by stating that we measure a negative differential stiffness for subcritical temperatures and a positive differential stiffness for supercritical temperatures. This behavior is influenced by the cooperativity, which has the capability to make the critical temperatures of the unfolding processes more uniform. We can therefore state that a positive cooperativity increases the resistance to fluctuations, making the spinoidal intervals equally stable to temperature variations.

While being a paradigmatic model for the understanding of several phenomena, our chain with Ising interactions should be improved to better represent more realistic situations. One drawback concerns the uniformity of all parameters defining the properties of the units. Indeed, in order to correctly model the actual mechanical behavior of heterogeneous structures, such as proteins,
we would have the possibility to freely choose these pa-
rameters for each unit. Nevertheless, this heterogeneity
consists in a form of quenched disorder, which is much
more complicated to be taken into account by classical
statistical mechanics methods. However, it should be im-
portant to introduce this point since it could allow to de-
termine the full unfolding pathway, which depends on the
system microstructure. As an example, this is directly
related to the biological function of a protein. Another
improvement concerns the dynamics of the unfolding pro-
cesses, which should be studied in the context of the out-
of-equilibrium statistical mechanics. It is worth noting
here that the spin variable approach can be used for de-
coupling two kinds of characteristic times: (i) the purely

mechanical times induced by the stiffness of each basin
of the potential energy, and (ii) the times induced by the
transition rates between the basins, which depend on the
energy barrier as classically described by the Kramers
theory. This approach should permit to consider out-
of-equilibrium unfolding processes, typically induced in
isometric force-extension experiments conducted at fixed
pulling velocity of the tethered chain.

ACKNOWLEDGMENTS

We acknowledge the region “Hauts de France” for the
financial support under project MEPOFIB.

[1] C. A. Hunter and H. L. Anderson, Angew. Chem. Int. Ed. **48**, 7488 (2009).

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

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

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

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

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

[7] I. Boldog, A. B. Gaspar, V. Martínez, P. Pardo-Ibañez, V. Ksenofontov, A. Bhattacharjee, P. Gütlich, and J. A. Real, Angew. Chem. Int. Ed. **47**, 6433 (2008).

[8] H. Banerjee, S. Chakraborty, and T. Saha-Dasgupta, In-organics **5**, 47 (2017).

[9] H. Qian, Annu. Rev. Biophys. **41**, 179 (2012).

[10] E. Agliari, A. Barra, L. Dello Schiavo, and A. Moro, Sci. Rep. **6**, 36314 (2016).

[11] A. Bakk and J. Hye, Phys. A: Stat. Mech. App. **323**, 504 (2003).

[12] D. A. Smith, D. J. Brockwell, R. C. Zinober, A. W. Blake, G. S. Beddard, P. D. Olmsted, S. E. Radford, Phil. Trans. R. Soc. Lond. A **361**, 713 (2003).

[13] B. Chakrabarti and A. J. Levine, Phys. Rev. E **71**, 031905 (2005).

[14] P. Malhotra, and J. B. Udgaonkar, Protein Sci. **25**, 1924 (2016).

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

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

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

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

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

[20] O. K. Dudko, Quarterly Reviews of Biophysics **49**, 1 (2016).

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

[22] C. Bustamante, J. Liphardt and F. Ritort, Physics Today **58**, 43 (2005).

[23] E. Dieterich, J. Camunas-Soler, M. Ribezzi-Crivellari, U. Seifert, and F. Ritort, Phys. Rev. E **94**, 012107 (2016).

[24] T. R. Strick, M.-N. Dessinges, G. Charvin, N. H. Dekker, J.-F. Allemand, D. Bensimon and V. Croquette, Rep. Progr. Phys. **66**, 1 (2003).

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

[26] K. C. Neuman, A. Nagy, Nature Meth. **5**, 491 (2008).

[27] S. Kumar, M. S. Li, Physics Reports **486**, 1 (2010).

[28] H. Miller, Z. Zhou, J. Shepherd, A. J. M. Wollman, and M. C. Leake, Rep. Prog. Phys. **81**, 024601 (2018).

[29] C. Yamahata, D. Collard, B. Legrand, T. Takekawa, M. Kumemura, G. Hashiguchi, H. Fujita, J. Microelectromech. Syst. **17**, 623 (2008).

[30] T. E. Fisher, A. F. Oberhauser, M. Carrion-Vazquez, P. E. Marszalek, and J. M. Fernandez, Trends Biochem. Sci. **24**, 379 (1999).

[31] H. Li, A. F. Oberhauser, S. B. Fowler, J. Clarke, and J. M. Fernandez, Proc. Nat. Acad. Sci. **97**, 6527 (2000).

[32] A. Imparato, F. Sbrana, and M. Vassalli, Europhys Lett. **82**, 58006 (2008).

[33] M. Bonin, R. Zhu, Y. Klaue, J. Oberstrass, E. Oesterschulze, and W. Nellen, Nucleic Acids Res. **30**, e81 (2002).

[34] J. Lipfert, G. M. Skinner, J. M. Keegstra, T. Hensgens, T. Jager, D. Dulin, M. Kber, Z. Yu, S. P. Donkers, F.-C. Chou, R. Das, and N. H. Dekker, Proc. Nat. Acad. Sci. **111**, 15408 (2014).

[35] S. B. Smith, L. Finzi, and C. Bustamante, Science **258**, 1122 (1992).

[36] J. F. Marko, E. D. Siggia, Macromolecules **28**, 8759 (1995).

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

[38] K. R. Chaurasiya, T. Paramanathan, M. J. McCauley, and M. C. Williams, Phys. Life Rev. **7**, 299 (2010).

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

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

[41] R.G. Winkler, J. Chem. Phys. **118**, 2919 (2003).
[42] F. Manca, S. Giordano, P. L. Palla, R. Zucca, F. Cleri, and L. Colombo, J. Chem. Phys. **136**, 154906 (2012).

[43] F. Manca, S. Giordano, P. L. Palla, F. Cleri, L. Colombo, J. Chem. Phys. **137**, 244907 (2012).

[44] T. Su and P. K. Purohit, J. Mech. Phys. Solids **58**, 164 (2010).

[45] J. Kierfeld, O. Niamploy, V. Sa-Yakanit, and R. Lipowsky, Eur. Phys. J. E **14**, 17 (2004).

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

[47] H. J. Kreuzer, S. H. Payne, Phys. Rev. E **63**, 021906 (2001).

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

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

[50] N. Bosaeus, A. H. El-Sagheer, T. Brown, S. B. Smith, B. Akerman, C. Bustamante, and B. Nordén, Proc. Natl. Acad. Sci. USA **109**, 15179 (2012).

[51] J. Wang, T. B. Kouznetsova, R. Boulatov, and S. L. Craig, Nature Comm. **7**, 13433 (2016).

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

[53] S. Cocco, J.F. Marko, R. Monasson, A. Sarkar, and J. Yan, Eur. Phys. J. E **10**, 249 (2003).

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

[55] F. Hanke and H. J. Kreuzer, Eur. Phys. J. E **22**, 163 (2007).

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

[57] D. B. Staple, S. H. Payne, A. L. C. Reddin, and H. J. Kreuzer, Phys. Biol. **6**, 025005 (2009).

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

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

[60] D. De Tommasi, N. Millardi, G. Puglisi, and G. Sacco-mandi, J. R. Soc. Interface **10**, 20130651 (2013).

[61] D. E. Makarov, Biophys. J. **96**, 2160 (2009).

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

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

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

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

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

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

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

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

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

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

[72] A. Prados, L. L. Bonilla and A. Carpio, J. Stat. Mech. **2010**, P06016 (2010).

[73] L. L. Bonilla, A. Prados and A. Carpio, J. Stat. Mech. **2010**, P09019 (2010).

[74] L. L. Bonilla, A. Carpio, A. Prados, and R. R. Rosales, Phys. Rev. E **85**, 031125 (2012).

[75] R. Law, P. Carl, S. Harper, P. Dalhaimer, D. W. Speicher, and D. E. Discher, Biophys. J. **84**, 533 (2003).

[76] T. Xu, H. Lannon, S. Wolf, F. Nakamura, and J. Brujic, Biophys. J. **104**, 2022 (2013).

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

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

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

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

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

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

[83] R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic Press, London, 1982).

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

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

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

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

[88] F. R. Gantmacher, The theory of matrices (Chelsea Publishing Company, New York, 1960).

[89] M. Radiom and M. Borkovec, Phys. Rev. E **96**, 062501 (2017).

[90] Y. Hu, T. Lindstrøm, B. Øksendal, J. Ubøe, and T. Zhang, Inverse Powers of White Noise, Proceedings of Symposia in Pure Mathematics, Vol. 57: Stochastic Analysis (Michael C. Cranston, Mark A. Pinsky, Eds.), American Mathematical Society (1995).