# Continuum vs thermodynamical limit in Statistical Mechanics

Giuseppe Florio, Stefano Giordano, Giuseppe Puglisi

► To cite this version:

Giuseppe Florio, Stefano Giordano, Giuseppe Puglisi. Continuum vs thermodynamical limit in Statistical Mechanics. Continuum Mechanics and Thermodynamics, 2025, 37, 10.1007/s00161-025-01406-8. hal-05300814

HAL Id: hal-05300814

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

Submitted on 7 Oct 2025

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
ORIGINAL ARTICLE

Giuseppe Florio · Stefano Giordano · Giuseppe Puglisi

# Continuum vs thermodynamical limit in Statistical Mechanics

Received: 19 February 2025 / Accepted: 17 July 2025 / Published online: 31 July 2025
© The Author(s) 2025

**Abstract** Determining the limiting behavior of discrete systems with a large number of particles in Statistical Mechanics is crucial for developing accurate analytic models, especially when addressing multistability and multiscale effects. Typically, one considers the so called thermodynamical limit or the continuum limit. The guiding principle for selecting the correct limit is to preserve essential properties of the discrete system, including physical attributes such as the interplay between enthalpic and entropic contributions, the influence of boundary conditions, and possible other energetic contributions such as interface effects. In this sense, an important role is played by the fundamental constants. Selecting appropriate rescaling factors for the Planck and Boltzmann constants, according to the specific limit considered, is a key theoretical concern. Despite the importance of this problem, the existing literature often lacks clarity on how different rescalings affect model accuracy. This work aims to clarify these issues by examining classical lattice models – particularly those that exhibit multistable behavior – and by proposing suitable limit rescalings to retain the discrete model’s material response when the number of particles increases.

## 1 Introduction

The important conceptual relation between the atomistic and continuum picture in mechanics and thermodynamics has been extensively discussed from a theoretical point of view. Already at the end of the 19th century, Boltzmann pointed out [1]:

> <<The question whether matter is atomistically constituted or continuous therefore reduces to the question: Which represents the observed properties of matter most accurately, the properties on the assumption of an extremely large finite number of particles, or the limit of the properties if the number grows infinitely large?>>

On the other hand such deep problem has been longly questioned from a mathematical point of view, and is at the root of the sixth Hilbert problem [2]:

G. Florio · G. Puglisi (✉)
Department of Civil Environmental Land Building Engineering and Chemistry, DICATECh, Polytechnic University of Bari, via Orabona 4, 70125 Bari, Italy
E-mail: giuseppe.puglisi@poliba.it

G. Florio
E-mail: giuseppe.florio@poliba.it

G. Florio
INFN, Section of Bari, 70126 Bari, Italy

S. Giordano
University of Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France, UMR 8520 - IEMN - Institut d'Électronique, de Microélectronique et de Nanotechnologie, F-59000 Lille, France
E-mail: stefano.giordano@univ-lille.fr
<<*Thus Boltzmann's work on the principles of mechanics suggests the problem of developing mathematically the limiting processes, there merely indicated, which lead from the atomistic view to the laws of motion of continua.*>>

The question of a correct analysis of this limit has been faced in different fields of science. From a physical point of view the importance of the answer to this question is related both to a correct deduction of the macroscopic thermodynamic laws, starting from the statistical mechanical framework or statistical field theory [3–6], and to the relevant aspect of the modern multiscale approaches to mechanics of materials [7–13]. Also, crucial mathematical aspects are related to this problem, where the theories of homogenization and $\Gamma$-convergences constitute theoretical frameworks aimed at the deduction of continuum, homogenized limits, starting from discrete microstructures [14–20]. In particular, the concept of $\Gamma$-limit, introduced by the Italian mathematician Ennio De Giorgi, helps in solving the important physical problem of determining limit energy functionals having as (global) minima the limit of the energy (global) minima of the discrete problems.

Basically two different approaches can be considered, according to the physical context: the *thermodynamic limit* and the *continuum limit*. The *thermodynamic limit* conceptually considers a system with an increasing number of particles keeping fixed the density of the system (number of particles over volume) and determining the equation of states and the constitutive laws based on the evaluation of the limiting value of the free energy *density*. Historically this large limit behavior has been at the base of the kinetic theory of ideal and non ideal gases and allowed to develop the theoretical framework of equilibrium phase transitions [22, 23]. Roughly speaking we can argue that the main advantage of this approach is the possibility of keeping a real atomistic architecture of the system whereas the limitation can be seen in the loss of boundary conditions effects and in the difficulty of comparing the (growing to infinity) total energy with the energy of local phenomena (interfaces, defects, etc.). This problem is sometimes solved by considering a periodic distribution of interfaces or defects, but this is only an apparent solution since this method typically introduces so called “spurious correlations”.

On the other hand, the *continuum limit* consists in considering a fixed volume with an increasing number of particles. With a correct rescaling of the system parameters the total energy remains finite thus losing the discrete structure of the system. In this case previous limitations are overcome, but the microstructure is neglected as in classical homogenization techniques. Moreover, as will be fully discussed in the following, the correct description of both kinetic energy and entropy leads to a divergent behavior due to the equipartition theorem. This limitation must be solved by rescaling the Planck and Boltzmann constants, as suggested in previous literature [3–5], keeping in the limit finite values for kinetic energy and entropy. Historically, the continuum limit obtained from the molecular theory of elasticity, as proposed by Navier [25] and Cauchy [26], is fully discussed in the classical treatise of Love on the mathematical theory of elasticity [27]. Already in that case the important problem of neglecting boundary conditions leads to the unphysical theory of rariconstant elasticity [28], previously discussed in terms of central (two-body) or multi-body interaction potentials [27]. The possibility of considering boundary effects in the continuum approximation can be crucial in several context such as grain growth and interface effects in phase transitions [21, 29, 30]. We notice that, in absence of a rescaling, the divergence observed in the thermodynamical potentials in the limit of infinite $n$ is consistent with the approach used in the paper based on the framework of classical Statistical Mechanics, in particular with the classical equipartition theorem. On the other hand, in the quantum case the situation would require the extension of the equipartition theorem as observed in Refs. [31, 32].

A related problem, which has been highly debated in the physics literature concerns the statistical ensembles describing the different possible conditions under which a system evolves. As well known, starting from the same Hamiltonian, the boundary conditions define the statistical ensembles and the corresponding exploration of the phase space. Different statistical ensembles lead to different macroscopic behavior for finite number of particles of the system. In this respect, when focussing on the Helmholtz ensemble case (e.g. assigned volume in a gas, extension in a chain, strain in a solid) or on the Gibbs ensemble case (e.g. assigned pressure in a gas, force in a chain, stress is a solid), different results regarding the equivalence or inequivalence of the ensembles in the continuum and thermodynamic limit can be investigated (see [33]). For example, in Ref. [34] the equivalence of the two ensembles in the thermodynamic limit for a chain with a continuous nearest-neighborhood interaction has been demonstrated. This result has been extended to non-nearest neighbor interaction for describing phase transformation is solids with interfacial energy effects [21]. Examples of non equivalences have been observed for more complex lattice geometries such as for adhesion of films [35], adhesion with softening mechanism [36], confined systems [37–40], desorption of polymers [42, 42, 43], and fracture models [51]. Continuous models of adhesion also show the non-equivalence of statistical ensembles, along with the presence of phase transitions [52]. In the perspective of describing the correct limit of a discrete system, whose behavior differs in the two ensembles, we argue that the equivalence in the thermodynamic or continuum limit of the two
<figure><img src="image_2.png" /><figcaption>Fig. 1 Scheme of a system made of *n* non-interacting particles confined on a line</figcaption></figure>

ensembles may come from a non correct consideration of the boundary conditions. In the following, we will discuss in details several examples on how to respect the different limit behavior reflecting the discrete model response. In particular, we show that in the continuum limit we can more realistically account for the boundary conditions, leading to inequivalence for some cases where thermodynamic limits are equivalent. The issue becomes more subtle when systems are also inequivalent in the thermodynamic limits.

To get insight in the described problems, we consider different limits in some explicit paradigmatic examples: a one-dimensional ideal gas, a mass-spring chain, and a multistable mass-spring chain with interfaces in the so called ‘zipper approximation’. While in the first example we have only purely entropic forces, in the second one we also introduce elastic interactions. In the third example, we also introduce bistability in the chain units and consider possible interactions among the units themselves. These models are paradigmatic examples allowing the possibility to obtain analytical formulas and, more important as we will see in the following, to fully understand the role of approximations and rescaling of physical quantities and constants. Moreover, one-dimensional models have been able to capture fundamental properties and relations at the macroscale by considering a proper analysis of costitutive parameters of materials and physical systems in two and three dimensions [21,53]. In all cases, we analyze the Gibbs ensemble characterized by an applied force and the Helmholtz ensemble defined by a prescribed extension of the system. To this aim, we introduce a Laplace-Fourier transform between the Gibbs and Helmholtz partition functions, useful to define a theoretical link between different statistical ensembles. In all cases, we show that the relevant thermodynamic quantities remain finite in both the thermodynamic and continuum limits. In the latter case, we must introduce an appropriate rescaling of the Planck and Boltzmann constants, coherent with the previous literature [3–5]. In the third example, we present the features of a chain made of bistable units with next-to-nearest neighbor interaction allowing the presence of interfaces. In this case, we use the saddle point approximation, thus showing how in the continuum limit and with a proper parameters rescaling, the interface nucleation and coalescence is still present and observable.

## 2 The one-dimensional ideal gas

In this Section we consider the case of a one-dimensional system made of non-interacting particles. We derive the partition function both in the case of fixed external force (Gibbs ensemble) and fixed position of the end-point particle on the right (Helmholtz). In both cases we can obtain the relevant thermodynamical potentials (free energy, internal energy or enthalpy) as well as entropic contribution and derive the relation between force and size of the system, finally considering the different expression in the thermodynamical and continuum limit. In subsection 2.3 we explicit the relation between the partition functions of the two ensembles in terms of Laplace transform.

Consider first the simplest paradigmatic example of a one-dimensional ideal gas composed by *n* particles of mass *m*, moving on the real axis with positions $x_1, x_2, ..., x_n$, and momenta $p_1, p_2, ..., p_n$ (see Figure 1). Let us assume that the particles can collide without passing through each other: $0 \le x_1 \le x_2, x_1 \le x_2 \le x_3, ...,$ and $x_{n-2} \le x_{n-1} \le x_n$. We consider two different boundary conditions (i.e. statistical ensembles).

### 2.1 Gibbs Ensemble

In the first case we apply a (negative) force *f* to the last particle at $x_n$ (isotensional condition, Gibbs ensemble), and in the second case we prescribe a fixed value for $x_n$ (isometric condition, Helmholtz ensemble). Thus, in the Gibbs ensemble we calculate the expected (average) value of the position (volume) $\langle x_n \rangle$ as function of *f*, and in the Helmholtz ensemble we determine $\langle f \rangle$ in terms of $x_n$.

In the case of the Gibbs ensemble the extended Hamiltonian of the system is

$$H_G = H_K - f x_n, \quad (1)$$
where the first sum represents the kinetic energy

$$H_K = \sum_{i=1}^{n} \frac{p_i^2}{2m}, \qquad (2)$$

and the second term describes the potential energy associated to the applied force. The partition function for this system can be calculated as

$$Z_G(f) = \frac{1}{h^n} \int_{\mathbb{R}^n} \int_{\Omega} \exp\left(-\frac{H_G}{k_B T}\right) dx_1 \, dx_2 \dots dx_n \, dp_1 \, dp_2 \dots dp_n, \qquad (3)$$

where $h$ is the Planck constant, $k_B$ is the Boltzmann constant, and $T$ is the temperature.

The presence of Planck's constant makes the partition function dimensionless and can be justified as follows. The exact quantum partition function can be written as $Z_G^Q = \sum_{i=0}^{+\infty} \exp(-\epsilon_i/k_B T)$. Here, $\epsilon_i$ are the eigenvalues of the Hamiltonian operator, defined by the equation $\hat{H}_G \Psi_i = \epsilon_i \Psi_i$, where $\Psi_i$ are the corresponding eigenfunctions. If we develop the expression for $Z_G^Q$ for small values of $h$, the first term is exactly given by Eq.(3). This result and the following terms of the development, representing the first quantum corrections, have been largely studied in the first half of the twentieth century by Wigner [54], Uhlenbeck and Gropper [55], and Kirkwood [56]. More details can be found in the Landau's textbook [57].

In Eq.(3), the set $\Omega \subset \mathbb{R}^n$ is characterized by $0 \le x_1 \le x_2, x_1 \le x_2 \le x_3, \dots, \le x_{n-2} \le x_{n-1} \le x_n$. Therefore, through the change of variables $\xi_1 = x_1, \xi_2 = x_2 - x_1, \dots, \xi_n = x_n - x_{n-1}$, we can consider $\bar{\Omega} = \{\xi_i \ge 0, i = 1, \dots, n\}$ and we get

$$\begin{aligned} Z_G(f) &= \frac{1}{h^n} \int_{\mathbb{R}^n} \exp\left(-\frac{1}{k_B T} \sum_{i=1}^n \frac{p_i^2}{2m}\right) dp_1 \, dp_2 \dots dp_N \int_{\bar{\Omega}} \exp\left(\frac{f \sum_{i=1}^n \xi_i}{k_B T}\right) d\xi_1 \, d\xi_2 \dots d\xi_N \\ &= \left(\frac{\sqrt{2\pi m k_B T}}{h}\right)^n \left(-\frac{k_B T}{f}\right)^n = \left(-\frac{k_B T}{f} \sqrt{\frac{2\pi m k_B T}{h^2}}\right)^n, \end{aligned} \qquad (4)$$

where the condition $f < 0$ ensures the convergence of the integral. It follows that the Gibbs free energy is given by

$$\mathcal{G} = -k_B T \ln Z_G = -k_B T n \ln \left( -\frac{k_B T}{hf} \sqrt{2\pi m k_B T} \right). \qquad (5)$$

Moreover, the entropy can be evaluated as

$$\mathcal{S} = -\frac{\partial \mathcal{G}}{\partial T} = n k_B \ln \left( -\frac{k_B T}{hf} \sqrt{2\pi m k_B T} \right) + \frac{3}{2} n k_B, \qquad (6)$$

and the enthalpy of the system is

$$\mathcal{H} = \mathcal{G} + T\mathcal{S} = \frac{3}{2} n k_B T. \qquad (7)$$

The term $\frac{3}{2} n k_B T$, shown in the enthalpy expression, describes two contributions: the quantity $\frac{1}{2} n k_B T$ represents the equipartition theorem applied to the kinetic terms associated to the $n$ particles; the quantity $n k_B T$ represents the thermal fluctuations associated to the potential energy of the applied force (this second term is double the previous one, as it must be).

Consider then the limit of growing number $n$ of particles and, in particular, the thermodynamic and continuum limits to understand how the considered functions converge to finite values in the two cases.

In the thermodynamic limit $n$ tends to infinity by increasing the extension of the system (keeping all physical constants and the particle density fixed). Thus, $h, k_B$, and $m$ are independent of $n$ as well as the densities of Gibbs free energy $\mathcal{G}/n$, entropy $\mathcal{S}/n$, and enthalpy $\mathcal{H}/n$, which converge to a constant for $n \to \infty$. This convergence can be seen immediately from the formulas above.

On the other hand, in the continuum limit $n$ increases while keeping the total size of the system fixed, and thus increasing the particle density. It means that the total mass $M = n m$ must remain constant, and therefore we have a varying particle mass $m = M/n$. In this case, the values of the Gibbs free energy $\mathcal{G}$, the entropy
S, and the enthalpy H have to converge to a constant for $n \to \infty$. This can be obtained only if $nk_B$ and $n^2h$ are constant for $n \to \infty$. These conditions are verified if we apply the rescaling laws $k_B = \tilde{k}_B/n$ and $h = \tilde{h}/n^2$, where $\tilde{k}_B$ and $\tilde{h}$ are constants to be determined by comparison with the real behavior of the discrete system to be approximated by the continuum limit. The result that the Boltzmann constant tends to zero with an increasing particle density is due to the fact that the energy of the thermal fluctuations would tend to infinity. In other words, this indicates that without a correct rescaling we are overestimating this term. Furthermore, the convergence of Planck's constant to zero does not mean that we are applying a classical limit to the study of the system, but that, once again, the particle density is divergent and therefore we must introduce a specific rescaled Planck constant. These arguments are consistent with previous investigations [3–5]. In particular, it is easy to verify that our result $k_B^2/h = \tilde{k}_B^2/\tilde{h}$ valid for one-dimensional geometries can be generalized to $k_B^{d+1}/h^d = \tilde{k}_B^{d+1}/\tilde{h}^d$ for the $d$-dimensional case. We notice that the ratio of the Plank and Boltzmann constants also appears in the context of models for the evaluation of the coefficient of friction at the nanoscale [59]. Also in that context it is possible to observe that the rescaling here introduced can be used to properly define physical quantities in the thermodynamical and continuum limits.

To conclude the analysis of this first paradigmatic example, we can study the force-extension relation for our one-dimensional gas under isotensional condition and compare the behavior obtained in the thermodynamic and continuum limit. It is well known that within the Gibbs ensemble the force-extension relation is given by

$$ \langle x_n \rangle = - \frac{\partial \mathcal{G}}{\partial f} = - \frac{k_B T n}{f}, \qquad (8) $$

where, as before, we consider $f < 0$.

As expected, in the thermodynamic limit the total extension $\langle x_n \rangle$ at assigned force diverges whereas the average strain $\langle x_n \rangle / n$ remains finite. It follows that the force-extension relation assumes the form $\langle x_n \rangle / n = -k_B T / f$, where both left and right hand sides are finite also for $n \to \infty$. On the other hand, in the continuum limit $\langle x_n \rangle$ remains finite and the force-extension relation $\langle x_n \rangle = -k_B T n / f$ is composed by finite left and right hand sides since $k_B n$ is constant for $n \to \infty$.

## 2.2 Helmholtz Ensemble

Let us now consider the Helmholtz ensemble, describing the case of assigned total displacement $x_n$. Since the last particle cannot move, the space of configurations $\Omega \subset \mathbb{R}^{n-1}$ is defined by $0 \le x_1 \le x_2, x_1 \le x_2 \le x_3, ..., \le x_{n-2} \le x_{n-1} \le x_n$, with $x_n$ fixed. The Helmholtz Hamiltonian is simply given by the kinetic energy of the first $n-1$ moving particles

$$ H_H = \sum_{i=1}^{n-1} \frac{p_i^2}{2m}, \qquad (9) $$

and, therefore, the Helmholtz partition function is defined as

$$ Z_H(x_n) = \frac{1}{h^{n-1}} \int_{\mathbb{R}^{n-1}} \int_{\Omega} \exp \left( -\frac{H_H}{k_B T} \right) dx_1 dx_2 \dots dx_{n-1} dp_1 dp_2 \dots dp_{n-1}. \qquad (10) $$

We adopt the change of variables $\xi_1 = x_1, \xi_2 = x_2 - x_1, \dots$, and $\xi_{n-1} = x_{n-1} - x_{n-2}$, and we can write the set $\Omega$ in terms of the new variables $\xi_i$ as $\bar{\Omega} = \{\xi_i \ge 0, i = 1, \dots, n-1, \sum_{i=1}^{n-1} \xi_i \le x_n\}$, and we get

$$ Z_H(x_n) = \left( \frac{\sqrt{2\pi m k_B T}}{h} \right)^{n-1} \int_{\bar{\Omega}} d\xi_1 d\xi_2 \dots d\xi_{n-1}, \qquad (11) $$

where $\bar{\Omega}$ represent a $(n-1)$-dimensional simplex. Since $\int_{\bar{\Omega}} d\xi_1 d\xi_2 \dots d\xi_{n-1} = x_n^{n-1}/(n-1)!$, we have

$$ Z_H(x_n) = \frac{1}{(n-1)!} \left( \sqrt{\frac{2\pi m k_B T}{h^2}} x_n \right)^{n-1}, \qquad (12) $$
<table><caption>**Table 1** Summary of the parameters and constitutive relations for a system of non-interacting particles (one-dimensional gas) in the thermodynamic and continuum limit</caption><thead><tr><th colspan="3">Non-interacting particles (1D ideal gas)</th></tr><tr><th></th><th>Thermodynamic limit</th><th>Continuum limit</th></tr></thead><tbody><tr><td>Mass</td><td>m = M/n fixed</td><td>M fixed</td></tr><tr><td>Planck constant</td><td>h fixed</td><td>ħ = n<sup>2</sup>h fixed</td></tr><tr><td>Boltzmann constant</td><td>k<sub>B</sub> fixed</td><td>k̃<sub>B</sub> = nk<sub>B</sub> fixed</td></tr><tr><td>Parameter (Gibbs)</td><td>⟨x<sub>n</sub>⟩ / n finite</td><td>⟨x<sub>n</sub>⟩ finite</td></tr><tr><td>Relation (Gibbs)</td><td>⟨x<sub>n</sub>⟩<sub>n</sub> = - k<sub>B</sub>T / f</td><td>⟨x<sub>n</sub>⟩ = - k̃<sub>B</sub>T / f</td></tr><tr><td>Parameter (Helmholtz)</td><td>x<sub>n</sub>/n finite</td><td>x<sub>n</sub> finite</td></tr><tr><td>Relation (Helmholtz)</td><td>⟨f⟩ = - k<sub>B</sub>T(n-1) / x<sub>n</sub></td><td>⟨f⟩ ~ - k̃<sub>B</sub>T / x<sub>n</sub></td></tr></tbody></table>

which is defined for $x_n \ge 0$. The Helmholtz free energy is therefore obtained as

$$ \mathcal{F} = -k_B T \ln Z_H = -k_B T(n-1) \ln \left( \frac{\sqrt{2\pi m k_B T}}{h} x_n \right) + k_B T \ln(n-1)! \quad (13) $$

We can use the Stirling approximation $\ln n! \sim n \ln(n/e)$, valid for $n \gg 1$ and, thus, both in the thermodynamic and in continuum limit this approximation leads to

$$ \mathcal{F} \sim -k_B T(n-1) \ln \left( \frac{\sqrt{2\pi m k_B T}}{h} \frac{x_n}{n} e \right). \quad (14) $$

Thus the entropy is given by

$$ \mathcal{S} = -\frac{\partial \mathcal{F}}{\partial T} \sim k_B(n-1) \ln \left( \frac{\sqrt{2\pi m k_B T}}{h} \frac{x_n}{n} e \right) + \frac{1}{2} k_B(n-1), \quad (15) $$

and the internal energy is evaluated as

$$ \mathcal{U} = \mathcal{F} + T\mathcal{S} \sim \frac{1}{2}k_B T(n-1). \quad (16) $$

The value of $\mathcal{U}$ could have been written directly from Eq.(9), based on the energy equipartition theorem. Again, it is clear that in the thermodynamic limit the Helmholtz free energy $\mathcal{F}/n$, the entropy $\mathcal{S}/n$, and the energy $\mathcal{U}/n$, converge to a constant for $n \to \infty$. Note that in this case the extension diverges and thus $x_n/n$ converges to a finite value. It is interesting to note that the $n$ in the denominator of $x_n$ comes from the logarithm of $(n-1)!$, and thus, finally, from the integration over the simplex domain. It is simple to verify that in the continuum limit the Helmholtz free energy $\mathcal{F}$, the entropy $\mathcal{S}$, and the energy $\mathcal{U}$, converge to a constant for $n \to \infty$, provided that $m = M/n$, $k_B = k̃_B/n$, and $h = ħ/n^2$, as before.

The force-extension relation within the Helmholtz ensemble is obtained as

$$ \langle f \rangle = \frac{\partial \mathcal{F}}{\partial x_n} = - \frac{k_B T (n-1)}{x_n}, \quad (17) $$

which has the same form as Eq.(8) except that now we determine the average value of $f$ at assigned $x_n$. As before, this relation is consistent with both thermodynamic and continuum limits.

In this simple example concerning a one-dimensional ideal gas, the Gibbs and Helmholtz force-extension relations assume exactly the same form for $n \to \infty$. It means that the ensembles are equivalent under both thermodynamic and continuum limits. This point can also be underlined by observing that, for $n \to \infty$, $\mathcal{H}$ in the Gibbs ensemble and $\mathcal{U}$ in the Helmholtz one are connected through a Legendre transform $\mathcal{H} = \mathcal{U} - f x_n$, or by using the free energies $\mathcal{G} = \mathcal{F} - f x_n$. Recall that while the Legendre transform between energies is only satisfied, for $n$ large, when the ensembles are equivalent, the Laplace-Fourier transform between partition functions is always valid [23, 24, 34], as discussed in the following Section.

We summarize the obtained results, the correct rescalings and the differences between the thermodynamic and continuum limit in Table 1.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Scheme of a one-dimensional system made of *n* particles interacting by an elastic force. (a) System subject to a fixed force. (b) System subject to a fixed total extension</figcaption></figure>

## 2.3 Relation between Gibbs and Helmholtz partition functions

The general expressions for the partition functions within the Helmholtz and the Gibbs ensembles are given in Eqs.(3) and (10), which are valid for any one-dimensional system. By comparing these equations, it easily follows that

$$Z_G(f) = \sqrt{\frac{2\pi mk_B T}{h^2}} \int_{-\infty}^{+\infty} Z_H(x_n) \exp\left(\frac{f x_n}{k_B T}\right) dx_n, \quad (18)$$

where the square root in front of the integral exactly corresponds to the inverse of the thermal De Broglie wavelength. In general, the domain of integration encompasses the entire real axis but can be restricted according to the cases under study. This relationship, except for the multiplicative constant, corresponds to a Laplace transform between the partition functions. We can verify that this relationship is correct for the partition functions obtained for the ideal gas. One can substitute Eq.(12) into Eq.(18), to get

$$Z_G(f) = \frac{1}{(n-1)!} \left( \frac{\sqrt{2\pi mk_B T}}{h} \right)^{n-1} \int_0^{+\infty} (x_n)^{n-1} \exp\left(\frac{f x_n}{k_B T}\right) dx_n, \quad (19)$$

where we recall that $x_n > 0$ and $f < 0$ (see Eqs.(4) and (12)). By changing variable $f x_n = -t k_B T$, and by using the Euler Gamma function defined as [58,60]

$$\Gamma(z) = \int_{0}^{+\infty} t^{z-1} e^{-t} dt, \quad (20)$$

with $\Gamma(n) = (n - 1)!$, we exactly obtain Eq.(4). Interestingly enough, it is possible to find the Helmholtz partition function from the one in the Gibbs ensemble by evaluating the inverse Laplace transform. Thus, we can verify that Eq.(12) is exactly derived from Eq.(4) by using the same method. This procedure will be generalized in the following sections by using the Fourier transform.

# 3 The one-dimensional mass-spring chain

Consider now the simple system composed of a mass-spring chain represented in Figure 2. As in the case of non-interacting particles, we consider the case of assigned force (Gibbs ensemble) in subsection 3.1 and assigned elongation $x_n$ in subsection 3.2. After evaluating the partition function in the former case, we derive it in the latter by using the Fourier transform. We thus obtain the constitutive relations both in the thermodynamic and continuum limit.

## 3.1 Gibbs ensemble

We start by analyzing the behavior of the system within the Gibbs ensemble. Therefore, we can write the extended Hamiltonian in the form

$$H_G = \sum_{i=1}^{n} \frac{p_i^2}{2m} + \sum_{i=1}^{n} \frac{1}{2} k_0 (x_i - x_{i-1} - \ell)^2 - f x_n, \quad (21)$$
where $x_0 = 0$, $k_0$ is the elastic constant, and $\ell$ is the resting length (or equilibrium length) of each spring. The Gibbs partition function can be evaluated as

$$Z_G(f) = \frac{1}{h^n} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \exp \left[ -\frac{1}{k_B T} \left( \sum_{i=1}^n \frac{p_i^2}{2m} + \sum_{i=1}^n \frac{1}{2} k_0 (x_i - x_{i-1} - \ell)^2 - f x_n \right) \right] d p_1 d p_2 \dots d p_n \, dx_1 \, dx_2 \dots dx_n, \quad (22)$$

where the integrals can be elaborated through the expression

$$\int_{-\infty}^{+\infty} e^{-\alpha x^2} e^{\beta x} dx = \sqrt{\frac{\pi}{\alpha}} e^{\frac{\beta^2}{4\alpha}} \quad (\text{for } \alpha > 0). \quad (23)$$

Straightforward calculations eventually lead to

$$Z_G(f) = \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right)^n \exp \left( \frac{n\ell f}{k_B T} + \frac{nf^2}{2k_B T k_0} \right). \quad (24)$$

Thus, the Gibbs free energy is given by

$$\mathcal{G} = -k_B T \ln Z_G = -k_B T n \ln \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right) - n \left( \ell f + \frac{f^2}{2k_0} \right). \quad (25)$$

We can determine the entropy

$$\mathcal{S} = -\frac{\partial \mathcal{G}}{\partial T} = k_B n + k_B n \ln \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right), \quad (26)$$

and the enthalpy

$$\mathcal{H} = \mathcal{G} + T\mathcal{S} = k_B T n - n \left( \ell f + \frac{f^2}{2k_0} \right). \quad (27)$$

Consider then also for this system the thermodynamic and continuum limits. As already seen in the case of the ideal gas, the *thermodynamic limit* is characterized by an increasing extension of the system, performed by keeping all the physical constants fixed. Thus, in this case $h, k_B, k_0, \ell$ and $m$ are independent from $n$, whereas the *densities* of Gibbs free energy $\mathcal{G}/n$, entropy $\mathcal{S}/n$, and enthalpy $\mathcal{H}/n$ converge to a constant for $n \to \infty$. This convergence can be easily seen from previously determined expressions.

Otherwise, in the *continuum limit*, the physical parameters scale as follows: $\ell = L/n$, where $L$ is the system length, $m = M/n$, where $M$ is the total mass, $k_B = \tilde{k}_B/n$, $h = \tilde{h}/n^2$, and $k_0 = nAE/L$, where $A$ is the system cross-section and $E$ is the Young's modulus. When such assumptions are taken into account, the values of the *total* Gibbs free energy $\mathcal{G}$, the entropy $\mathcal{S}$, and the enthalpy $\mathcal{H}$ converge to a constant for $n \to \infty$, as expected.

The force-extension relation for the mass-spring chain is then obtained by deriving the Gibbs free energy, as follows

$$\langle x_n \rangle = -\frac{\partial \mathcal{G}}{\partial f} = n \left( \ell + \frac{f}{k_0} \right). \quad (28)$$

On the one hand, in the thermodynamic limit, we can write the average strain $\langle x_n \rangle / n = \ell + f/k_0$, with both left and right hand sides that are convergent to finite values. On the other hand, in the continuum limit, the force-extension relation assumes the form $\sigma = E\epsilon$, where $\epsilon = (\langle x_n \rangle - L)/L$ is the strain, and $\sigma = f/A$ is the stress. We then obtain the classical linear response of the one-dimensional solid mechanics.
## 3.2 Helmholtz ensemble

As already seen in the previous Section, it is possible to find a direct connection between the partition functions in the different ensembles, see Eq.(18). In particular, one can use the Fourier transform to derive the Helmholtz partition function from the Gibbs ensemble [23]. This approach is often adopted with respect to a direct calculation because, typically, the calculation of the Gibbs partition function is simpler.

To this aim, we firstly introduce in Eq.(18) the variable g such that $f = -ik_B T g$ (i is the imaginary unit), and we get

$$Z_G(-ik_B T g) = \sqrt{\frac{2\pi mk_B T}{h^2}} \int_{-\infty}^{+\infty} Z_H(x_n) \exp(-igx_n) dx_n. \quad (29)$$

This expression corresponds to a Fourier transform and can be inverted as it follows

$$Z_H(x_n) = \frac{1}{2\pi} \sqrt{\frac{h^2}{2\pi mk_B T}} \int_{-\infty}^{+\infty} Z_G(-ik_B T g) \exp(+igx_n) dg. \quad (30)$$

A further change of variable $\eta = -k_B T g$ leads to the final expression

$$Z_H(x_n) = \frac{1}{2\pi k_B T} \sqrt{\frac{h^2}{2\pi m k_B T}} \int_{-\infty}^{+\infty} Z_G(i\eta) \exp\left(-i\frac{\eta x_n}{k_B T}\right) d\eta, \quad (31)$$

which allows the determination of $Z_H$ if we know $Z_G$. It is important to underline that this expression is based on the analytic continuation of the Gibbs partition function over the imaginary axis.

We can use Eq.(31) by substituting the Gibbs partition function given by Eq.(24) so that we obtain

$$\begin{aligned} Z_H(x_n) &= \frac{1}{2\pi k_B T} \sqrt{\frac{h^2}{2\pi m k_B T}} \int_{-\infty}^{+\infty} \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right)^n \exp \left( i \frac{n\ell\eta}{k_B T} - \frac{n\eta^2}{2k_B T k_0} - i \frac{\eta x_n}{k_B T} \right) d\eta, \\ &= \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right)^{n-1} \frac{1}{\sqrt{n}} \exp \left[ -\frac{k_0}{2nk_B T} (n\ell - x_n)^2 \right], \end{aligned} \quad (32)$$

where we used the Gaussian integral

$$\int_{-\infty}^{+\infty} e^{-\alpha x^2} e^{i\beta x} dx = \sqrt{\frac{\pi}{\alpha}} e^{-\frac{\beta^2}{4\alpha}}, \quad (\text{for } \alpha > 0). \quad (33)$$

The Helmholtz partition function allows to determine the Helmholtz free energy, entropy, and internal energy of the system, as it follows:

$$\mathcal{F} = -k_B T \ln Z_H = -k_B T (n - 1) \ln \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right) + \frac{k_B T}{2} \ln n + \frac{1}{2} n k_0 \left( \ell - \frac{x_n}{n} \right)^2, \quad (34)$$

$$\mathcal{S} = -\frac{\partial \mathcal{F}}{\partial T} = k_B(n-1) \ln \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right) - \frac{k_B}{2} \ln n + k_B(n-1), \quad (35)$$

$$\mathcal{W} = \mathcal{F} + T\mathcal{S} = k_B T (n-1) + \frac{1}{2} n k_0 \left( \ell - \frac{x_n}{n} \right)^2. \quad (36)$$

We can easily recognize in the internal energy the thermal fluctuation term (which corresponds to $n-1$ kinetic terms and $n-1$ potential terms governed by the equipartition theorem) and the mechanical energy term corresponding to the $n$ chain springs. Regarding the convergence in the *thermodynamic limit*, as in previous case, the densities of free energy, entropy and internal energy are finite in the thermodynamic limit where all physical parameters of the system are constant. When instead we consider the *continuum limit*, the total free energy, entropy and internal energy remain finite if we apply the rescaling $\ell = L/n$, $m = M/n$, $k_B = \tilde{k}_B/n$, $h = h/n^2$, and $k_0 = nAE/L$ (note that the term proportional to $\ln n$ in the free energy and entropy does not play any role since $\ln n/n \to 0$ when $n \to \infty$).

To conclude we can evaluate the force-extension relation for the chain as

$$\langle f \rangle = \frac{\partial \mathcal{F}}{\partial x_n} = k_0 \left( \frac{x_n}{n} - \ell \right), \quad (37)$$
<table><caption>**Table 2** Summary of the parameters and constitutive relations for a mass-spring chain in the thermodynamic and continuum limit</caption><thead><tr><th colspan="3">Mass-spring chain</th></tr><tr><th></th><th>Thermodynamic limit</th><th>Continuum limit</th></tr></thead><tbody><tr><td>Mass</td><td><i>m</i> = <i>M</i>/<i>n</i> fixed</td><td><i>M</i> fixed</td></tr><tr><td>Geometric parameter</td><td><i>l</i> = <i>L</i>/<i>n</i> fixed</td><td><i>L</i> fixed</td></tr><tr><td>Elastic parameter</td><td><i>k</i><sub>0</sub> fixed</td><td><i>E</i> = <sup><i>k</i><sub>0</sub><i>L</i></sup>/<sub><i>n</i><sup>2</sup><i>A</i></sub> fixed</td></tr><tr><td>Planck constant</td><td><i>h</i> fixed</td><td><i>ħ</i> = <i>n</i><sup>2</sup><i>h</i> fixed</td></tr><tr><td>Boltzmann constant</td><td><i>k</i><sub>B</sub> fixed</td><td><i>ħ</i><sub>B</sub> = <i>nk</i><sub>B</sub> fixed</td></tr><tr><td>Parameter (Gibbs)</td><td>⟨<i>x</i><sub>n</sub>⟩ / <i>n</i> finite</td><td>⟨ε⟩ = ⟨<i>x</i><sub>n</sub>⟩ - <i>L</i> / <i>L</i> finite</td></tr><tr><td>Relation (Gibbs)</td><td>⟨<i>x</i><sub>n</sub>⟩ / <i>n</i> = <i>l</i> + <i>f</i>/<i>k</i><sub>0</sub></td><td>⟨ε⟩ = <i>σ</i>/<i>E</i> with σ = <i>f</i>/<i>A</i></td></tr><tr><td>Parameter (Helmholtz)</td><td><i>x</i><sub>n</sub>/<i>n</i> finite</td><td><i>ε</i> = <i>x</i><sub>n</sub><i>L</i>/<i>L</i> finite</td></tr><tr><td>Relation (Helmholtz)</td><td>⟨<i>f</i>⟩ = <i>k</i><sub>0</sub> (⟨<i>x</i><sub>n</sub>⟩ - <i>l</i>)</td><td>⟨σ⟩ = <i>Eε</i> with ⟨σ⟩ = ⟨<i>f</i>⟩ / <i>A</i></td></tr></tbody></table>

which corresponds exactly to what was found in the Gibbs ensemble. This proves that also in this simple elastic extension of previous example there is a statistical equivalence of the ensembles in the thermodynamic and continuum limit. Again, this result is consistent with the fact that there is a Legendre transformation between the two ensembles, as can be easily verified [23].

We summarize the obtained results, the corresponding rescalings, and the differences between the thermodynamic and continuum limit in Table 2.

## 4 One-dimensional multistable mass-spring chain with interfacial energy

We introduce in this Section, a model composed of bistable units with an interaction scheme described by an Ising type energy, see Figure 3. To simplify the analysis of the system, we adopt the zipper assumption, which describes a system divided into two parts, each of them being in a different phase (for instance, folded and unfolded conformational phase). These two parts are separated by a single domain wall that can evolve either due to the external mechanical action or thermal fluctuations. This model can be used to describe phase transformations in the mechanics of solids or the configurational transition of macromolecules in the biophysical context. This system highlights the possible important differences between the continuous and thermodynamic limits. In particular, we observe that the energy associated with the domain wall is localized and must compete with the total energy of the system. Since in the thermodynamic limit the total energy diverges, such competition can only be observed in the continuous limit. This difference is then studied in detail in what follows as a paradigmatic example of the role of local energy terms in the limit of growing number of particles. Thus, this represents a first important example where the two different approaches can lead to important differences in the system response.

### 4.1 Model of a mass-spring chain with bistable units

Consider a chain made of *n* bistable units, each subject to a next-to-nearest neighbor interaction [21]. In particular, the configuration of the chain is described by the Hamiltonian *H*, defined as

$$H = H_K + H_{el} + H_J, \qquad (38)$$

where the different terms are defined as follows. The kinetic energy is given by

$$H_K = \frac{1}{2} \sum_{i=1}^{n} \frac{p_i^2}{m}. \qquad (39)$$

The multistable elastic potential energy $H_{el}$ is defined as

$$H_{el} = \sum_{i=1}^{n} \left\{ \frac{1}{2} k(S_i) [x_i - x_{i-1} - l_0(S_i)]^2 + Q(S_i) \right\}, \qquad (40)$$
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Scheme of a one-dimensional system of n particles interacting by multi-stable elastic forces characterized by different elastic constants and natural lengths. (a) System subject to a fixed force (Gibbs ensemble). (b) System subject to a fixed total extension (Helmholtz ensemble). Bottom figures: Energy profile of a multi-stable elastic energy with each phase defined by the value of the spin variable $S_i$. Central panel: energy wells with the same elastic constant ($\alpha = 1$). Left and right panels: energy wells with different elastic constant ($\alpha < 1$ representing a softening transition behavior and $\alpha > 1$, representing a hardening transition behavior, respectively)</figcaption></figure>

where we consider the boundary condition $x_0 = 0$. Each term in $H_{el}$ can exhibit two phases, each characterized by using the spin variable $S_i$ that can assume the values $\pm 1$, where $S_i = -1$ corresponds to the first well (folded element), whereas $S_i = +1$ corresponds to the second well (unfolded element). In general, the two phases are characterized by different values of the elastic constants $k(S_i)$, energy gaps $\Delta E$ between the minima of the potential wells measuring the transition energy, and different natural lengths $l_0(S_i)$. In particular, we fix

$$l_0(S_i) = \begin{cases} l, & \text{if } S_i = -1 \\ \chi l, & \text{if } S_i = +1 \end{cases}, \quad k(S_i) = k_0\alpha(S_i) = \begin{cases} k_0, & \text{if } S_i = -1 \\ \alpha k_0, & \text{if } S_i = +1 \end{cases}, \quad Q(S_i) = \begin{cases} 0, & \text{if } S_i = -1 \\ \Delta E, & \text{if } S_i = +1 \end{cases}, \quad (41)$$

where $\chi > 1$, $l > 0$, $k_0 > 0$, $\alpha > 0$, $\Delta E > 0$. In Figure 3 we show the energetic profile of a single bistable unit.

Finally, we model the non-local interaction by considering the Ising energy term

$$H_J = -J \sum_{i=1}^{n-1} S_i S_{i+1}, \quad (42)$$

where the coupling constant $J > 0$ characterizes the interaction strength. The choice of a positive value for $J$ corresponds to favor phase coalescence in the configuration of the chain (interface energy penalization).

In order to simplify the analysis and obtain analytically simple expressions, we consider here a modified version of the Hamiltonian. In particular, as previously anticipated, we focus on the so-called *zipper model* [61–63]. Within this approximation, a single interface between the folded and unfolded parts of the chain is nucleated and propagated along the system (single domain wall approximation). Thus, we introduce the discrete variable $\xi \in \{0, \dots, n\}$ representing the number of elements in the unfolded state (correspondingly, we have $n - \xi$ units in the folded state).

In the following we distinguish the cases when the two wells have the same elastic constant ($\alpha = 1$) or different elastic constant ($\alpha \neq 1$) because they lead to a different transition behavior.

## 4.2 Identical wells

Consider first the case with $\alpha = 1$, when the elastic constants of the different potential wells coincide and only the reference lengths of the two material phases are different.
### 4.2.1 Zipper model in the Gibbs ensemble

When the external force $f$ is assigned, the Hamiltonian describing the system reads

$$H_G(\xi) = \sum_{i=1}^{n} \frac{p_i^2}{2m} + \sum_{i=1}^{n-\xi} \left\{ \frac{k_0}{2} (x_i - x_{i-1} - l)^2 \right\} + \sum_{i=n-\xi+1}^{n} \left\{ \Delta E + \frac{\alpha k_0}{2} (x_i - x_{i-1} - \chi l)^2 \right\} - J[n-1-2I(\xi)] - f x_n \quad (43)$$

where $I(\xi)$ represents the number of interfaces corresponding to the assigned value of $\xi$. In particular, within the zipper model, we have $I(0) = I(n) = 0$ (no interface) and $I(\xi) = 1$ if $1 \le \xi \le N-1$ (only one interface). We also stress that the quantity $n-\xi$ represents the position of the domain wall in the chain. Finally, we notice that $H_G$ depends explicitly on the number of unfolded elements $\xi$.

The partition function for the Gibbs ensemble can be evaluated as in the previous Sections. We have

$$Z_G = \frac{1}{h^n} \sum_{\xi=0}^{n} \int_{\mathbb{R}^n} \int_{\mathbb{R}^n} \exp \left[ -\frac{H_G(\xi)}{k_B T} \right] d p_1 d p_2 \dots d p_n dx_1 dx_2 \dots dx_n, \quad (44)$$

where we included the summation over $\xi$ in order to take into account the different configurations of the interface between folded and unfolded sections of the chain. An explicit calculations gives

$$Z_G = \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right)^n \exp \left[ \frac{J(n-1)}{k_B T} \right] \sum_{\xi=0}^{n} \exp \left[ -\frac{2J}{k_B T} I(\xi) \right] \Gamma_G(\xi), \quad (45)$$

where

$$\Gamma_G(\xi) = \exp \left[ -\frac{\Delta E}{k_B T} \xi \right] \exp \left\{ \frac{1}{k_B T} \left[ \frac{nf^2}{2k_0} + fl[n + (\chi - 1)\xi] \right] \right\}. \quad (46)$$

As in the previous case, the Gibbs free energy is

$$\mathcal{G} = -k_B T \ln Z_G. \quad (47)$$

The relation between the applied force and the expectation value of the total elongation of the chain is

$$\langle x_n \rangle = - \frac{\partial \mathcal{G}}{\partial f}, \quad (48)$$

finally giving

$$\frac{\langle x_n \rangle}{n} = \frac{f}{k_0} + l \left(1 - \frac{\langle \xi \rangle}{n}\right) + \chi l \frac{\langle \xi \rangle}{n}, \quad (49)$$

where the expectation value of the unfolded fraction $\langle \rho \rangle = \langle \xi \rangle/n$ is defined as

$$\langle \rho \rangle = \frac{\langle \xi \rangle}{n} = \frac{1}{n} \frac{\sum_{\xi=0}^{n} \xi \exp \left[ -\frac{2J}{k_B T} I(\xi) \right] \Gamma_G(\xi)}{\sum_{\xi=0}^{n} \exp \left[ -\frac{2J}{k_B T} I(\xi) \right] \Gamma_G(\xi)}. \quad (50)$$

We notice that the relation between force and extension remains valid both in the thermodynamic limit and in the continuum limit where it assumes the form

$$\frac{\sigma}{E} = \frac{\langle x_n \rangle}{L} - (1 - \langle \rho \rangle) - \chi \langle \rho \rangle. \quad (51)$$
### 4.2.2 Zipper model in the Helmholtz ensemble

Consider now the case of the Helmholtz ensemble when the value of the extension $x_n$ is assigned and we want to determine the expectation value of the force $\langle f \rangle$ experienced by the chain. We proceed by using the Fourier transform method outlined in Section 3. We obtain

$$ 
\begin{aligned}
Z_H(x_n) &= \frac{1}{2\pi k_B T} \sqrt{\frac{h^2}{2\pi m k_B T}} \int_{-\infty}^{+\infty} Z_G(i\eta) \exp\left(-i \frac{\eta x_n}{k_B T}\right) d\eta \\
&= \frac{1}{\sqrt{n}} \left(\frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}}\right)^{n-1} \exp\left[\frac{J(n-1)}{k_B T}\right] \sum_{\xi=0}^{n} \exp\left[-\frac{2J}{k_B T} I(\xi)\right] \Gamma_H(\xi),
\end{aligned} 
\quad (52) $$

where

$$ \Gamma_H(\xi) = \exp\left[-\frac{\Delta E}{k_B T} \xi\right] \exp\left\{-\frac{k_0}{2nk_B T} \left[x_n - l(n + (\chi - 1)\xi)\right]^2\right\}. \quad (53) $$

Evaluating the Helmholtz free energy

$$ \mathcal{F} = -k_B T \ln Z_H, \quad (54) $$

gives us the possibility to obtain the force-elongation relation in this ensemble as

$$ \langle f \rangle = \frac{\partial \mathcal{F}}{\partial x_n}. \quad (55) $$

We get

$$ \frac{\langle f \rangle}{k_0} = \frac{x_n}{n} - l \left(1 - \frac{\langle \xi \rangle}{n}\right) - \chi l \frac{\langle \xi \rangle}{n} = \frac{x_n}{n} - l(1 - \langle \rho \rangle) - \chi l \langle \rho \rangle, \quad (56) $$

where the expectation value of the unfolded fraction $\langle \rho \rangle = \langle \xi \rangle / n$ in the Helmholtz ensemble is defined as

$$ \langle \rho \rangle = \frac{\langle \xi \rangle}{n} = \frac{1}{n} \frac{\sum_{\xi=0}^{n} \xi \exp \left[ -\frac{2J}{k_B T} I(\xi) \right] \Gamma_H(\xi)}{\sum_{\xi=0}^{n} \exp \left[ -\frac{2J}{k_B T} I(\xi) \right] \Gamma_H(\xi)}. \quad (57) $$

After some algebraic manipulation we obtain

$$ \langle \rho \rangle = \frac{\langle \xi \rangle}{n} = \frac{1}{n} \frac{n \Gamma_H(N) \left[ \exp\left(\frac{2J}{k_B T}\right) - 1 \right] + \sum_{\xi=0}^{n} \xi \Gamma_H(\xi)}{[\Gamma_H(0) + \Gamma_H(N)] \left[ \exp\left(\frac{2J}{k_B T}\right) - 1 \right] + \sum_{\xi=0}^{n} \Gamma_H(\xi)}. \quad (58) $$

We finally have

$$ \frac{\langle f \rangle}{E A} = \lambda - (1 - \langle \rho \rangle) - \chi \langle \rho \rangle = \lambda - 1 - (\chi - 1) \langle \rho \rangle, \quad (59) $$

where

$$ \lambda = \frac{x_n}{L} \quad (60) $$

and the function $\Gamma_H$ reads

$$ \Gamma_H(\xi) = \exp\left[-\frac{q}{k_B T} \frac{\xi}{n}\right] \exp\left\{-\frac{A E L}{2k_B T} \left[1 + (\chi - 1)\frac{\xi}{n} - \lambda\right]^2\right\}, \quad (61) $$

with $q = n \Delta E$. This rescaling is consistent with the fact that the energy corresponding to the transition of the whole system $n \Delta E$ remains fixed as $n$ grows. Thus, this represents the correct scaling in the continuum limit.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Panels (a)-(b): expectation value of $\sigma/E$ versus $\lambda$ (with $\sigma = f/A$). In all cases we fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$. Panel (a): we used $n = 100, 200, 300, 400$, $\beta = 100$, $J = 0.01$. Panel (b): we used $n = 100, 200, 300, 400$, $\beta/n = 1$, $J = 0.01$. Panel (c): expectation value of $f/k_0$ versus $x_n/n$, see Eq.(56); we fixed $k_0 = l = \beta = \Delta E = J = 1$, $\chi = 5$, and $n = 20, 50, 80$</figcaption></figure>

In Figure 4, we show the stress-stretch behavior of the system for increasing number of elements. In order to show the behavior in the $n \to \infty$ limit, we use the rescaling adopted in Section 3 for the material parameters in Figure 4 (a) (for the thermodynamic limit) and in Figure 4 (b) (for the continuum limit). We also remark that since in the thermodynamic limit the total stiffness decreases to zero, we rescale with $L$ the total stress. Observe that the nucleation stress peak decreases in the thermodynamic limit rescaling with increasing $n$, see Figure 4 (a). This is an artificial effect of the rescaling making the contribution of the the ‘localized’ interfacial energy negligible with respect to the total energy of the system, which diverges to infinity when $n \to \infty$. As a result, the nucleation peak disappears in the thermodynamic limit. On the other hand, the behavior in the continuum limit rescaling, represented in Figure 4 (b) shows that the nucleation peak is slightly increasing with $n$. It is possible to show that this peak converges to a given value for $n \to \infty$, as discussed in the next Section. For comparison, in Figure 4 (c) we plot the expectation value of $f/k_0$ versus $x_n/n$, see Eq.(56). Having in mind the thermodynamic limit with $n, L \to \infty$, we fix $l$ and consider increasing values of $n$. Also in this case we observe that the nucleation peak decreases for larger values of $n$.

It is worth noticing that in phase transition phenomena the presence of stress peaks, distinguishing nucleation vs interfaces propagation stresses represents a crucial aspect to be described [21,29]. Therefore, models that want to describe such phenomena must be able to represent this behavior. For this reason, the continuum limit becomes essential in this field of applications.

### 4.2.3 Continuum limit in the Helmholtz ensemble

Consider now the continuum limit of the multistable system analyzed in terms of the saddle point approximation. Specifically, consider the classical Euler-MacLaurin approximation for a given function $\phi$ and substitute
the summation over $\xi$ with integrals as

$$ \sum_{\xi=0}^{n} \phi(\xi) \approx \int_{0}^{n} \phi(\xi) d\xi, \qquad (62) $$

where we neglected higher order terms in order to simplify the calculations and to make the final formulas more transparent. To exhibit the important role of the continuum limit as compared with the thermodynamic limit when interfacial effects – and, similarly, when boundary conditions take an important role in defining phase nucleation and evolution for systems with non-local interactions [29] –, in the following we evaluate explicitly the limit with $n \to +\infty$ with L fixed. We can thus rewrite the expectation value of the fraction $\langle\rho\rangle$ in Eq.(58) as

$$ \langle\rho\rangle \approx \frac{\Gamma_H(n) (e^{2Jn\tilde{\beta}} - 1) + n \int_0^1 \rho e^{-n\tilde{\beta}g(\rho)} d\rho}{(\Gamma_H(0) + \Gamma_H(n)) (e^{2Jn\tilde{\beta}} - 1) + n \int_0^1 e^{-n\tilde{\beta}g(\rho)} d\rho}, \qquad (63) $$

where we defined

$$ g(\rho) = q\rho + \frac{AEL}{2} [1 + (\chi - 1)\rho - \lambda]^2, \qquad (64) $$

and we introduced the rescaled inverse temperature

$$ \tilde{\beta} = \frac{1}{nk_B T} = \frac{1}{\tilde{k}_B T}. \qquad (65) $$

The definition of $\tilde{\beta}$ is consistent with the observations in previous Sections where we considered the continuum limit taking into account the rescaled Boltzmann constant as $k_B = \tilde{k}_B/n$.

In order to obtain the correct limit for $n \to +\infty$, we evaluate the asymptotic expression for $\langle\rho\rangle$ and, thus, the relation obtained from Eq.(59). Let us first determine the expectation value of the phase fraction for two-phases solution with $0 < \rho < 1$ corresponding to the conditions

$$ 1 + \frac{q}{AEL(\chi - 1)} < \lambda < \chi + \frac{q}{AEL(\chi - 1)}. \qquad (66) $$

In this case we can evaluate the integrals in Eq.(63) by the stationary phase method. By considering

$$ \frac{\partial g}{\partial \rho} = 0, \qquad (67) $$

we can find the solution

$$ \rho^* = \frac{\lambda - 1}{\chi - 1} - \frac{q}{AEL(\chi - 1)^2}. \qquad (68) $$

Thus, we obtain the asymptotic expression for the integrals in Eq.(63) as

$$ \int_0^1 \rho e^{-n\tilde{\beta}g(\rho)} d\rho \sim \rho^* \int_{-\infty}^{+\infty} e^{-n\tilde{\beta}[g(\rho^*) + \frac{1}{2}g''(\rho^*)(\rho - \rho^*)^2]} d\rho = \rho^* e^{-n\tilde{\beta}g(\rho^*)} \sqrt{\frac{2\pi}{n\tilde{\beta}g''(\rho^*)}}, \qquad (69) $$

$$ \int_0^1 e^{-n\tilde{\beta}g(\rho)} d\rho \sim \int_{-\infty}^{+\infty} e^{-n\tilde{\beta}[g(\rho^*) + \frac{1}{2}g''(\rho^*)(\rho - \rho^*)^2]} d\rho = e^{-n\tilde{\beta}g(\rho^*)} \sqrt{\frac{2\pi}{n\tilde{\beta}g''(\rho^*)}}, \qquad (70) $$

where

$$ g(\rho^*) = q \frac{\lambda - 1}{\chi - 1}, \quad g''(\rho^*) = AEL(\chi - 1)^2. \qquad (71) $$
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Panel (a): expectation (asymptotic) value $\langle\sigma\rangle_{as}/E$ versus $\lambda$ for increasing values of $n$. We fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$, $J = 0.1$ and $\tilde{\beta} = 0.1$. The curves correspond to $n = 200, 10^3, 10^4, 10^5$. Panel (b): comparison of the expectation value $\langle\sigma\rangle/E$ (solid line) and the asymptotic limit (dashed line) versus $\lambda$. We fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$, $J = 0.1$, $n = 100$ and $\tilde{\beta} = 1$</figcaption></figure>

We finally find

$$ \langle\rho\rangle_{as} \sim \frac{\rho^* + e^{-n\tilde{\beta}f_{\chi}(\lambda)} \sqrt{\frac{\tilde{\beta}g''(\rho^*)}{2\pi n}}}{1 + e^{-n\tilde{\beta}f_{\chi}(\lambda)} (1 + e^{-n\tilde{\beta}(f_1(\lambda)-f_{\chi}(\lambda))}) \sqrt{\frac{\tilde{\beta}g''(\rho^*)}{2\pi n}}} \quad (72) $$

$$ = \frac{\rho^* + (\chi - 1)e^{-n\tilde{\beta}f_{\chi}(\lambda)} \sqrt{\frac{\tilde{\beta}AEL}{2\pi n}}}{1 + (\chi - 1)e^{-n\tilde{\beta}f_{\chi}(\lambda)} (1 + e^{-n\tilde{\beta}(f_1(\lambda)-f_{\chi}(\lambda))}) \sqrt{\frac{\tilde{\beta}AEL}{2\pi n}}} \quad (73) $$

where

$$ f_1(\lambda) = \frac{AEL}{2}(\lambda - 1)^2 - 2J - g(\rho^*), \quad f_{\chi}(\lambda) = q + \frac{AEL}{2}(\lambda - \chi)^2 - 2J - g(\rho^*). \quad (74) $$

On the other hand, when

$$ \lambda \le 1 + \frac{q}{AEL(\chi - 1)} \quad \text{or} \quad \lambda \ge \chi + \frac{q}{AEL(\chi - 1)}, \quad (75) $$

we find, respectively,

$$ \langle\rho\rangle_{as} \sim \frac{(e^{2\tilde{\beta}Jn} - 1) e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2 - n\tilde{\beta}q}}{(e^{2\tilde{\beta}Jn} - 1) \left(e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2} + e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2 - n\tilde{\beta}q}\right) + ne^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - 1)^2}}, \quad (76) $$

$$ \langle\rho\rangle_{as} \sim \frac{(e^{2\tilde{\beta}Jn} - 1) e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2 - n\tilde{\beta}q} + ne^{-n\tilde{\beta}[q + \frac{AEL}{2}(\lambda - \chi)^2]}}{(e^{2\tilde{\beta}Jn} - 1) \left(e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2} + e^{-n\tilde{\beta} \frac{AEL}{2} (\lambda - \chi)^2 - n\tilde{\beta}q}\right) + ne^{-n\tilde{\beta}[q + \frac{AEL}{2}(\lambda - \chi)^2]}} \quad (77) $$

The force-stretch relation in the asymptotic (continuum) limit is thus given by

$$ \frac{\langle f \rangle_{as}}{EA} = \lambda - 1 - (\chi - 1)\langle\rho\rangle_{as}. \quad (78) $$

In Figure 5 (a), we plot the asymptotic expectation value $\langle\sigma\rangle_{as}/E$ versus $\lambda$ given by Eq. (78) for increasing values of $n$. We notice the the use of $\tilde{\beta}$ with the rescaling of the Boltzmann constant as $nk_B$ allows to observe the peaks corresponding to the nucleation and reabsorption (coalescence) of the domain wall generated by the coupling energy $J$.

In order to compare the asymptotic result with the original one, in Figure 5 (b) we plot the expectation value $\langle\sigma\rangle/E$ (solid line) and the asymptotic limit (dashed line) versus $\lambda$ (values of the parameters in the caption
of the figure). In agreement with the rescaling argument of the Boltzmann constant, we have to evaluate the exact result for $\beta = n\tilde{\beta}$. Having this in mind, we notice the excellent agreement of the results.

In order to get insight in the continuum limit, we focus our attention on the limit of the force peak when $n \to +\infty$. The analytic expression of this quantity can be obtained by considering that, for $\lambda \simeq 1$, Eq. (78) can be approximated as

$$ \frac{\langle f \rangle_{as}}{EA} = \frac{q}{EAL(\chi - 1)} + \left[ \lambda - 1 - \frac{q}{EAL(\chi - 1)} \right] \frac{ae^{-n\tilde{\beta}f_1}}{1 + ae^{-n\tilde{\beta}f_1}}, \quad (79) $$

where

$$ a = (\chi - 1)\sqrt{\frac{\tilde{\beta}EAL}{2\pi n}} (e^{2\tilde{\beta}Jn} - 1). \quad (80) $$

We can impose the condition to find the maximum value of the force as

$$ \frac{\partial}{\partial \lambda} \langle f \rangle_{as} = 0. \quad (81) $$

An explicit computation shows that this condition leads to the Lambert equation

$$ we^w = s, \quad (82) $$

where

$$ w = n\tilde{\beta}\frac{AEL}{2}\left(\lambda - 1 - \frac{q}{\chi - 1}\right)^2 - \frac{1}{2}, \quad (83) $$

$$ s = \frac{\chi - 1}{2} \sqrt{\frac{\tilde{\beta}EAL}{2\pi n e}} (e^{2\tilde{\beta}Jn} - 1). \quad (84) $$

The solution can be written as [64,65]

$$ w = W_0(s), \quad (85) $$

where $W_0$ is the so-called Lambert function. Using this result with Eq. (79), we can evaluate the force peak as

$$ \frac{\langle f \rangle_P}{EA} = \frac{q}{EAL(\chi - 1)} + \sqrt{\frac{2}{n\tilde{\beta}EAL}} \frac{W_0(s)}{\sqrt{\frac{1}{2} + W_0(s)}}, \quad (86) $$

When $s \gg 1$ (verified when $\tilde{\beta}$ is fixed for $n \to +\infty$) we have that $W_0(s) \sim \ln s$ and, thus, the force peak in the continuum limit can be written as

$$ \frac{\langle f \rangle_P}{EA} \sim \frac{q}{EAL(\chi - 1)} + 2\sqrt{\frac{J}{AEL}}. \quad (87) $$

The second term in the right-hand side corresponds to the force gap between the plateau (corresponding to the first term in the summation of the right-hand side) and the peak. Interestingly, we have obtained that the gap does not depend on temperature if we use the rescaled inverse temperature $\tilde{\beta}$ in the continuum limit. This result is drastically different with respect to the analysis of the continuum limit without considering $\tilde{k}_b$ constant when the force peaks goes to zero. We also point out the important difference in the behavior between the thermodynamic and the continuum limit. In the former case, the nucleation peak disappears when $n \to \infty$. In the second case the peak tends to a constant that depends on the Ising interaction factor, which represents the energy associated with the domain wall.
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Comparison of the zipper model (blue curve) with the results obtained from the total energy including the term in Eq.(42) with all possible configurations (red curve), see Appendix A. Panel (a), (b) and (c) corresponds to decreasing values of temperature (increasing values of $\beta$ as reported in each figure. We considered $n = 30$, $J = 0.1$, $E = 1$, $A = 1$, $L = 1$, $q = 0.5$, $\chi = 5$</figcaption></figure>

<table><caption>**Table 3** Summary of the parameters and costitutive relations for a multi-stable chain in the thermodynamic and continuum limit</caption><caption>Mass-spring bistable chain with Ising interactions</caption><thead><tr><th></th><th>Thermodynamic limit</th><th>Continuum limit</th></tr></thead><tbody><tr><td>Mass</td><td>$m = M/n$ fixed</td><td>$M$ fixed</td></tr><tr><td>Geometric parameter</td><td>$l = L/n$ fixed</td><td>$L$ fixed</td></tr><tr><td>Elastic parameter</td><td>$k_0$ fixed</td><td>$E = \frac{k_0 L}{n A}$ fixed</td></tr><tr><td>Planck constant</td><td>$h$ fixed</td><td>$\tilde{h} = n^2 h$ fixed</td></tr><tr><td>Boltzmann constant</td><td>$k_B$ fixed</td><td>$\tilde{k}_B = n k_B$ fixed</td></tr><tr><td>Parameter (Gibbs)</td><td>$\langle x_n \rangle / n$ finite</td><td>$\langle \lambda \rangle = \frac{\langle x_n \rangle}{L}$ finite</td></tr><tr><td>Relation (Gibbs)</td><td>$\frac{\langle x_n \rangle}{n} = \frac{f}{k_0} + l \left(1 - \frac{\langle \xi \rangle}{n}\right) + \chi l \frac{\langle \xi \rangle}{n}$</td><td>$\langle \lambda \rangle = \frac{\sigma}{E} + (1 - \langle \rho \rangle) + \chi \langle \rho \rangle$ with $\sigma = f/A$</td></tr><tr><td>Parameter (Helmholtz)</td><td>$x_n / n$ finite</td><td>$\lambda = \frac{x_n}{L}$ finite</td></tr><tr><td>Relation (Helmholtz)</td><td>$\frac{\langle f \rangle}{k_0} = \frac{x_n}{n} - l \left(1 - \frac{\langle \xi \rangle}{n}\right) - \chi l \frac{\langle \xi \rangle}{n}$</td><td>$\frac{\langle \sigma \rangle}{E} = \lambda - (1 - \langle \rho \rangle) - \chi \langle \rho \rangle$ with $\sigma = f/A$</td></tr><tr><td>Nucleation peak (Helmholtz)</td><td>$\frac{\langle f \rangle_p}{E A} - \frac{q}{E A L(\chi-1)} = 0$</td><td>$\frac{\langle f \rangle_p}{E A} - \frac{q}{E A L(\chi-1)} = 2\sqrt{\frac{J}{A E L}}$</td></tr></tbody></table>

### 4.2.4 Remark and summary

As a remark, we stress that the zipper model is an approximation of the real behavior of the system that includes the elastic energy and the total non-local interaction energy in Eq.(42). One should actually include the possibility of more barriers (i.e., more domain walls) and consider their effect in the evaluation of the macroscopic quantities. On the other hand, this would make the analytical treatment more difficult and the results less transparent in their physical interpretation. The analysis of the complete model is reported in Appendix A. In order to show the validity of the zipper approximation, we plot in Figure 6 the comparison of the two models. We can appreciate that the zipper model is in very good agreement with the complete model for the typical values of the parameters considered in our discussion.

In Figure 7 we summarize the behavior of the $\sigma - \lambda$ relation depending on the limit performed. The results obtained and the differences between the thermodynamic and continuum limit are reported in Table 3.

### 4.3 Different energy wells: the case with $\alpha \neq 1$

In the previous example the elastic constants of the different phases were equal, i.e., the parameter $\alpha$ in Eq.(40) was equal to one. In the case of $\alpha \neq 1$, it is possible to observe a different phenomenology. In particular, one observes that the Maxwell force phase transition plateau is temperature dependent [21]. We are going to consider the effect of the difference of the elastic constants also for the continuum limit using the rescaled inverse temperature $\tilde{\beta} = \beta/n$.

Starting from Eq.(43), we can obtain the partition function in the Gibbs ensemble as in the case with $\alpha = 1$. In particular, the results found for $\alpha = 1$ can be generalized as

$$Z_G^{(\alpha)} = \left( \frac{2\pi k_B T}{h} \sqrt{\frac{m}{k_0}} \right)^n \exp \left[ \frac{J(n-1)}{k_B T} \right] \sum_{\xi=0}^{n} \Omega_G^{(\alpha)}(\xi), \quad (88)$$
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Summary of the behavior of the $\sigma - \lambda$ relation. As explained in the main text, depending on the limit performed we observe different results in the diagrams. Starting from a finite value of $n$ (left branch of the diagram), performing the limit $n \to \infty$ with $k_B$ constant (right upper branch of the diagram) corresponds to the thermodynamic limit whereas the limit $n \to \infty$ with $\tilde{k}_B = nk_B$ constant (right lower branch of the diagram) corresponds to the continuum limit</figcaption></figure>

where

$$ \Omega_G^{(\alpha)}(\xi) = \frac{1}{\alpha^{\alpha/2}} \exp\left[-\frac{2JI(\xi)}{k_B T}\right] \exp\left[-\frac{\Delta E \xi}{k_B T}\right] \exp\left\{\frac{1}{k_B T} \left[\frac{nf^2}{2k_0}\left(1-\frac{\xi}{n}+\frac{\xi}{\alpha n}\right)+nfl\left(1+(\chi-1)\frac{\xi}{n}\right)\right]\right\}. \quad (89) $$

The Gibbs free energy and the force-elongation relation can be evaluated as

$$ \mathcal{G}^{(\alpha)} = -k_B T \ln Z_G^{(\alpha)}, \quad \langle x_n \rangle_{\alpha} = -\frac{\partial \mathcal{G}^{(\alpha)}}{\partial f}, \qquad (90) $$

where we defined the ($\alpha$ dependent) expectation value

$$ \langle \bullet \rangle_{\alpha} = \frac{\sum_{\xi=0}^{n} \bullet \Omega_{G}^{(\alpha)}(\xi)}{\sum_{\xi=0}^{n} \Omega_{G}^{(\alpha)}(\xi)}. \qquad (91) $$

We thus obtain the generalized expression for the force-elongation relation in the Gibbs ensemble as

$$ \frac{\langle x_n \rangle_{\alpha}}{n} = \frac{f}{k_0} \left( 1 + \frac{1-\alpha}{\alpha} \frac{\langle \xi \rangle_{\alpha}}{n} \right) + l \left( 1 - \frac{\langle \xi \rangle_{\alpha}}{n} \right) + \chi l \frac{\langle \xi \rangle_{\alpha}}{n}, \qquad (92) $$
which, having in mind the idea of analyzing the continuum limit, can be rewritten as

$$ \frac{\langle x_n \rangle_{\alpha}}{L} = \frac{\sigma}{E} \left( 1 + \frac{1-\alpha}{\alpha} \langle \rho \rangle_{\alpha} \right) + (1 - \langle \rho \rangle_{\alpha}) + \chi \langle \rho \rangle_{\alpha}, \quad (93) $$

where, as before, we defined $\rho = \xi/n$, $\sigma = f/A$, and we recall that $k_0 = nAE/L$.

In the Helmholtz ensemble (where the value of the extension $x_n$ is assigned) we can generalize the expression of the partition function obtained for $\alpha = 1$ as

$$ Z_{H}^{(\alpha)}(x_n) = \frac{1}{2\pi k_B T} \sqrt{\frac{h^2}{2\pi m k_B T}} \int_{-\infty}^{+\infty} Z_{G}^{(\alpha)}(i\eta) \exp\left(-i\frac{\eta x_n}{k_B T}\right) d\eta = \left(\frac{2\pi k_B T}{h}\sqrt{m}\right)^{n-1} e^{\frac{J(n-1)}{k_B T}} \sum_{\xi=0}^{n} \Omega_{H}^{(\alpha)}(\xi), \quad (94) $$

where

$$ \Omega_{H}^{(\alpha)}(\xi) = \frac{\exp\left[-\frac{2J}{k_B T} I(\xi)\right] \exp\left[-\frac{\Delta E}{k_B T} \xi\right]}{\sqrt{\left(\frac{n-\xi}{k_0} + \frac{\xi}{\alpha k_0}\right) k_0^{n-\xi} (\alpha k_0)^{\xi}}} \exp\left\{-\frac{1}{2k_B T} \left(\frac{n-\xi}{k_0} + \frac{\xi}{\alpha k_0}\right)^{-1} [x_n - l(n + (\chi - 1)\xi)]^2\right\}. \quad (95) $$

The Helmholtz free energy and the force-elongation relation in this ensemble read

$$ \mathcal{F}^{(\alpha)} = -k_B T \ln Z_H^{(\alpha)}, \quad \langle f \rangle_{\alpha} = \frac{\partial \mathcal{F}^{(\alpha)}}{\partial x_n}. \quad (96) $$

We obtain

$$ \frac{\langle f \rangle_{\alpha}}{k_0} = \left\langle \frac{1}{1 - \frac{\xi}{n} + \frac{\xi}{\alpha n}} \right\rangle_{\alpha} \frac{x_n}{n} - l \left\langle \frac{1 - \frac{\xi}{n}}{1 - \frac{\xi}{n} + \frac{\xi}{\alpha n}} \right\rangle_{\alpha} - \chi l \left\langle \frac{\frac{\xi}{n}}{1 - \frac{\xi}{n} + \frac{\xi}{\alpha n}} \right\rangle_{\alpha}, \quad (97) $$

where, from now on, we consider

$$ \langle \bullet \rangle_{\alpha} = \frac{\sum_{\xi=0}^{n} \bullet \Omega_{H}^{(\alpha)}(\xi)}{\sum_{\xi=0}^{n} \Omega_{H}^{(\alpha)}(\xi)}. \quad (98) $$

Having in mind the continuum limit, we can rewrite Eq. (97) as

$$ \frac{\langle \sigma \rangle_{\alpha}}{E} = \left\langle \frac{1}{1 - \rho + \frac{\rho}{\alpha}} \right\rangle_{\alpha} \lambda - l \left\langle \frac{1 - \rho}{1 - \rho + \frac{\rho}{\alpha}} \right\rangle_{\alpha} - \chi l \left\langle \frac{\rho}{1 - \rho + \frac{\rho}{\alpha}} \right\rangle_{\alpha}. \quad (99) $$

As in the case with $\alpha = 1$, we can evaluate the continuum limit by using the saddle point approximation in the limit $n \to +\infty$ with $L$ fixed. We find the asymptotic expectation value of $\sigma$ as

$$ \frac{\langle \sigma \rangle_{\alpha \text{ as}}}{E} \sim \frac{C_n \left[ (\lambda - 1) e^{n g_1} + (\lambda - \chi) \alpha^{3/2} e^{n g_\chi} \right] + \frac{f^*}{EAL}}{C_n \left[ e^{n g_1} + \alpha^{1/2} e^{n g_\chi} \right] + 1}, \quad (100) $$

where

$$ \frac{f^*}{AEL} = \frac{1}{1 - \rho^* + \rho^*/\alpha} \left[ (1 - \rho^*)(\lambda - 1) + \rho^*(\lambda - \chi) \right], \quad (101) $$

and $\rho^*$ is the solution of the saddle point equation $\frac{\partial g}{\partial \rho} = 0$, with

$$ g(\rho) = \tilde{\beta} q \rho + \frac{\tilde{\beta}}{2} \frac{AEL}{1 - \rho + \rho/\alpha} [\lambda - 1 - (\chi - 1)\lambda]^2 + \frac{\rho}{2} \ln \alpha, \quad \tilde{\beta} = \frac{1}{nk_B T}. \quad (102) $$
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Panel (a): expectation value $\langle\sigma\rangle_{as}/E$ versus $\lambda$ for increasing values of $\tilde{\beta}$. We fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$, $J = 0.01$, $n = 320$ and $\tilde{\beta} = 10, 7, 4$. Panel (b): asymptotic expectation value $\langle\sigma\rangle/E$ (solid line) versus $\lambda$. We fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$, $J = 0.01$, $n = 3000$ and $\tilde{\beta} = 10, 7, 4$. Inset: comparison of the expectation value $\langle\sigma\rangle/E$ (solid line) and the asymptotic limit (dashed line) versus $\lambda$. We fixed $E = A = L = 1$, $q = 0.5$, $\chi = 5$, $J = 0.01$, $n = 320$ and $\tilde{\beta} = 10$</figcaption></figure>

In these equations, we defined

$$C_n = \sqrt{\frac{\tilde{\beta} A E L}{2 \pi n}} \left(e^{2 \tilde{\beta} J n} - 1\right) \sqrt{\left[\left(1 + \frac{f^*}{A E L}\right) - \left(\chi + \frac{f^*}{\alpha A E L}\right)\right]^2}, \quad (103)$$

$$g_1 = -\frac{\tilde{\beta}}{2} A E L \left[ \lambda - \left(1 + \frac{f^*}{A E L}\right) \right]^2, \quad (104)$$

$$g_{\chi} = -\frac{\tilde{\beta}}{2} \alpha A E L \left[ \lambda - \left( \chi + \frac{f^*}{\alpha A E L} \right) \right]^2. \quad (105)$$

Focussing on the nucleation peak as in the case $\alpha = 1$, we can use again the Lambert function so to obtain in the continuum limit

$$\frac{\langle\sigma\rangle_{\alpha P}}{E} \sim \frac{f^*}{E A L} + 2\sqrt{\frac{J}{A E L}}. \quad (106)$$

This result is analogous to the case with $\alpha = 1$. In particular, the force gap with respect to the plateau is the same and, again, temperature independent. In Figure 8 we plot the relation between $\sigma$ and $\lambda$ in the case with $\alpha \neq 1$. In particular, in panel (a) we consider the case with $\alpha < 1$ where we observe that the stress plateau decreases by increasing the temperature (determined by the rescaled inverse temperature $\tilde{\beta}$), the case with $\alpha > 1$ being analogous but with an increasing plateau for increasing temperature. In panel (b) we plot the asymptotic value of the expectation value of $\sigma$ for a larger value of $n$. In this case we observe that, even considering larger values of temperature, the stress peak does not disappear. In the inset we show that the asymptotic formula of the stress perfectly matches the non approximated expression.

The important point is that the difference between the elastic constants of the two energy wells of the bistable unit generates a temperature-dependent plateau, a fact that is experimentally observed in systems that exhibit phase transformations. This is the novelty with respect to the case with $\alpha = 1$. In several experimentally observed cases, however, the nucleation peak has also been seen to be temperature dependent. This phenomenon cannot be reproduced from the model under consideration studied with the zipper approximation. It would be interesting to study whether the Ising-type interaction model, studied without the zipper approximation, i.e., with an arbitrary number of domain walls induced by thermal fluctuations, is able to provide a temperature-dependent nucleation peak (always in the continuous limit). Although this seems to us to be true from the analysis of some numerical results, we do not currently have a rigorous demonstration of this fact. We therefore leave as a future perspective the study of the continuum limit of the exact Ising model i.e., without zipper approximation.
## 5 Conclusions

In this paper, we dealt with the problem of moving from discrete to continuous models, including the effects of temperature and entropy. By considering systems made of an increasing number of elements *n*, we studied two different limits, i.e. the thermodynamical and the continuum limits. We analyzed two types of boundary conditions, representing two different statistical ensembles. The isotensional conditions correspond to the Gibbs ensemble whereas the isometric conditions correspond to the Helmholtz ensemble. As a matter of fact, the analysis of the behavior of discrete systems for increasing number of elements represents a fundamental topic in several theoretical and applied fields. From one side system with a large number of elements represent the standard case in several artificial and biological systems. Second the analysis of these limits are crucial to obtain analytical expressions especially in the field of multiscale phenomena (see e.g. [66]). Two different approaches have been considered in the literature. From one side the so-called thermodynamic limit corresponds to increasing the number of elements, with the costitutive parameters of the system independent on the rescaling. As a result all the extensive quantities (volume, entropy, energy, etc) grow to infinity keeping finite the corresponding densities. The main drawback of this approach is the difficulty of considering the competition between the local energy terms (such as the interface energy as considered in this work) and the role of boundary conditions. It is exactly this last aspect that leads to the classical results of statistical equivalence in the large *n* limit of the Gibbs and Helmholtz statistical ensemble in the thermodynamical limit. The second approach is to consider the so-called continuum limit, when the rescaling corresponds to increasing *n* while keeping finite both the densities and the total amount of the extensive energetic components. As we show in this work a non careful rescaling of this type can induce an overestimation of entropic and kinetic terms. Indeed, when we consider the continuum limit we have to introduce a rescaling of the two most important physical constants, namely the Boltzmann and Planck constants, to obtain the convergence of all the energetic functions. This result is consistent with the fact that as *n* increases the internal energy of the system grows indefinitely because of the equipartition theorem that associates a constant amount of energy to each degree of freedom. This is the ultimate reason that requires us to rescale the fundamental constants to regain convergence of the energy functions. As already stated this result is fully consistent with the framework of classical equilibrium Statistical Mechanics.

Specifically, by explicitly scrutinizing three distinct models of growing complexity, we revealed the interplay between these limits and the ensembles. Starting with a one-dimensional gas of non-interacting particles, we found that energy, enthalpy, free energy, and entropy densities converge in the thermodynamic limit. On the other hand, the convergence in the continuum limit requires the rescaling of the key constants due to the infinite energy growth predicted by the equipartition theorem. Thus, we observed that the densities of energy, enthalpy, free energy, and entropy converge in the thermodynamic and (constant rescaled) continuum limit. In this first example, the forces are purely entropic (there are no elastic forces in the system) and the two statistical ensembles are equivalent.

Adding elastic forces in a second model, we observed similar convergence criteria, with both Gibbs and Helmholtz ensembles yielding equivalent results.

To introduce greater complexity, in the third case, we introduced a chain of bistable units and considered an Ising-type interaction between them. By adopting the terminology used to describe macromolecules, particularly proteins, the two states of each unit of the system units mimic folded and unfolded configurations. We assume a positive constant *J* in the coupling constant of the units to represent ferromagnetic behavior, meaning that adjacent elements prefer to be in the same state (folded or unfolded). In this model, the general solution is not straightforward in terms of transparency of the results and, therefore, we used the so-called zipper assumption which introduces a single domain wall between the folded and unfolded regions. This approach largely simplifies the analysis of the problem. The results remain valid in the case of a strongly ferromagnetic system, i.e., with a large value of *J*. This system immediately exhibits more complex departing behavior with the Helmholtz ensemble showing a response with nucleation and coalescence force peaks. Notably, the thermodynamic limit eliminates these peaks while they remain in the continuum limit with rescaled Planck and Boltzmann constants due to the finite total energy of the system. Specifically, in the Helmholtz ensemble, the force peak is due to the nucleation of a domain wall, while in the Gibbs ensemble, this peak is always absent due to force control. We therefore focused on the isometric Helmholtz ensemble, which effectively describes this nucleation peak (as well as the opposed final coalescence peak). Other important differences emerge when observing the thermodynamic and continuum limits. Let us begin by discussing the case describing the two energy wells of the bistable units with two equal elastic constants. In this case, the nucleation (and coalescence) peak disappears in the thermodynamic limit as *n* tends to infinity. This can be understood by
observing that the domain wall energy is finite and localized, while the total energy of the system diverges in this limit. Therefore, the competition between the domain wall energy and other extensive energy components cannot manifest. However, in the continuum limit, the total energy of the system remains finite, and, thus, the nucleation peak also remains finite as *n* tends to infinity. As a result, we are able to keep the important difference between the nucleation and propagation observed in the Helmholtz ensemble when the continuum limit is considered. For this system, the height of the nucleation peak does not depend on temperature in the continuum limit. Furthermore, in this configuration, we observe that the plateau force – that is, the force that propagates the domain wall along the chain (also known as the Maxwell force) – also remains constant with respect to the system temperature.

We then considered the case of two energy wells with different elastic constants. In this case the situation changes partially. Indeed, we find a plateau that depends on temperature, while the nucleation peak remains constant and equal to the previously found value. It is interesting to note that in numerous experiments concerning phase transformations, it has been observed that both the Maxwell force and the nucleation peak depend on temperature. This fact leads us to the following reasoning. In studying this system, we introduced the zipper assumption to simplify the analytical treatment. This may have slightly altered the evaluation of the system's entropy because we reduced the number of domain walls to one. It is possible that by considering the exact complete model with an arbitrary number of domain walls induced by thermal fluctuations, a temperature-dependent nucleation peak would be found. Although some numerical evidence presented, we currently lack a theoretical proof on this point, which we leave as a future perspective for further exploration.

Our results show that the choice of the ensemble and the model assumptions, such as the zipper approximation, are essential in capturing realistic thermo-mechanical behaviors. Ultimately, this study highlights that both continuum and thermodynamic limits, alongside ensemble selection, shape the emergent physical behaviors of discrete models. These factors are fundamental when modeling systems that bridge microscopic and macroscopic scales, offering insights for theoretical approaches for modelling complex physical and mechanical phenomena.

We believe that the important theoretical differences of the thermodynamic and continuum limits (with non rescaled constants or with rescaled constants) represent different approximations that should be analyzed carefully. Indeed, the recalled rescaling properties of the different limits can induce several important conceptual drawbacks, leading to overestimation or underestimation of important physical quantities characterizing the phenomena with a resulting non effectiveness of the deduced limits. Since the analysis of these important aspects is in many case not well exploited in the scientific literature, we think that the results discussed in the analyzed examples can shed light onto this field.

**Acknowledgements** G.F. and G.P. have been supported by “Gruppo Nazionale per la Fisica Matematica” (GNFM) under “Istituto Nazionale di Alta Matematica” (INdAM). G.F. and G.P.’s research is funded by the European Union–Next Generation EU. G.F. and G.P. are supported by PNRR, National Center for HPC, Big Data and Quantum Computing–M4C2–I 1.4 (grant number CN0000013, CUP D93C22000430001)–Spoke 5 (Environment and Natural Disasters). G.P. is supported by the Project of National Relevance (PRIN), financed by European Union–Next-GenerationEU–NRRP–M4C2–I 1.1, CALL PRIN 2022 (Project 2022XLBLRX, CUP D53D23006020006) and CALL PRIN 2022 PNRR (Project P2022KHFBNB, CUP D53D23018910001) granted by the Italian MUR. G.F. is supported by the PRIN, financed by European Union–Next-Generation EU–NRRP–M4C2–I 1.1, CALL PRIN 2022 PNRR (Project P2022MXCJ2, CUP D53D23018940001) and CALL PRIN 2022 (Project 2022MKB7MM, CUP D53D23005900006) granted by the Italian MUR. G.F. is also supported by “Istituto Nazionale di Fisica Nucleare” (INFN) through the project QUANTUM. S.G. has been supported by “Central Lille” and “Région Hauts de France” under project Sta-MeNa (Statistical mechanics for macromolecular structures of nanotechnology).

**Author contributions** All the authors contributed equally

**Funding** Open access funding provided by Politecnico di Bari within the CRUI-CARE Agreement.

**Data Availability** No datasets were generated or analysed during the current study.

**Declarations**

**Competing interests** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## Appendix A Complete model without the zipper approximation

In this Appendix we consider the complete model beyond the zipper approximation used in Section 4. Having in mind the continuum limit of the discrete model, we define, in accordance with the choice made in Section 4.1

$$l_0(S_i) = \begin{cases} l, & \text{if } S_i = -1 \\ \chi l, & \text{if } S_i = +1 \end{cases}, \quad k(S_i) = k_0\alpha(S_i) = \begin{cases} k_0, & \text{if } S_i = -1 \\ \alpha k_0, & \text{if } S_i = +1 \end{cases}, \quad Q(S_i) = \begin{cases} 0, & \text{if } S_i = -1 \\ \Delta E, & \text{if } S_i = +1 \end{cases} \quad (A1)$$

where $\chi > 1$, $l > 0$, $k_0 > 0$, $\alpha > 0$, $\Delta E > 0$ and the relation between the Young modulus $E$ and the elastic constant $k_0$ is obtained by introducing the total length of the chain $L = nl$ and the section of the sample $A$ as

$$k_0 = \frac{EAn}{L}. \quad (A2)$$

Thus, we consider the rescaled energy obtained from the Hamiltonian in Eq. (38)

$$\bar{\Phi} = \frac{nH}{kL^2} = \frac{H}{EAL} = \frac{n}{2MEAL} \sum_{i=1}^{n} p_i^2 + \sum_{i=1}^{n} \left\{ \frac{1}{2} \frac{\alpha(S_i)}{n} [\lambda_i - \lambda_0(S_i)]^2 + \frac{1}{n} \frac{nQ(S_i)}{EAL} \right\} - \bar{J} \sum_{i=1}^{n-1} S_i S_{i+1}, \quad (A3)$$

where $\lambda_i = (x_i - x_{i-1})/l$, $\lambda_0(S_i) = l_0(S_i)/l$ and we defined $\bar{J} = \frac{J}{EAL}$.

We notice that if we apply a fixed force $f$ to the right end of the chain, we have to include the work performed by the force. Thus, it would be necessary to include an energy term of the form

$$H_f = f \frac{L}{n} \sum_{i=1}^{n} \lambda_i. \quad (A4)$$

By a rescaling analogous to the one used for $\bar{\Phi}$ we obtain

$$\bar{\Phi}_f = \frac{H_f}{EAL} = \frac{f}{nEA} \sum_{i=1}^{n} \lambda_i = \frac{f}{EA} \lambda. \quad (A5)$$

### Appendix A.1 Fixed force: Gibbs ensemble

In order to introduce the temperature effects, we have to compute the canonical partition in the case of assigned force $f$, within the Gibbs ensemble function:

$$\bar{Z}_G(f) = \frac{1}{h^n} \left(\frac{L}{n}\right)^n \sum_{\{S_i\}} \int_{\mathbb{R}^{2n}} e^{-\bar{\beta}(\bar{\Phi}-\bar{\Phi}_f)} d\lambda_1 \dots d\lambda_n d p_1 \dots d p_n. \quad (A6)$$

where we have defined the non dimensional rescaled inverse temperature

$$\bar{\beta} = \frac{EAL}{k_B T}, \quad (A7)$$

with $h$ the Planck constant, $k_B$ the Boltzmann constant and $T$ the absolute temperature. Following the approach described in [21] we find that

$$\bar{Z}_G(f) = \frac{1}{2 \cosh \bar{\beta} \bar{J}} \left( \frac{2\pi L \sqrt{MEAL}}{nh\bar{\beta}} \right)^n \left[ \lambda_+^n + \lambda_-^n + e^{-2\bar{\beta} \bar{J}} (\lambda_+^n - \lambda_-^n) \frac{\lambda_+ + \lambda_-}{\lambda_+ - \lambda_-} \right], \quad (A8)$$
where $\lambda_{\pm}$ are the eigenvalues of the transfer matrix $T$ defined as

$$T = \begin{bmatrix} e^{\bar{\beta}J} c_{-} & e^{-\bar{\beta}J} \sqrt{c_{+}c_{-}} \\ e^{-\bar{\beta}J} \sqrt{c_{+}c_{-}} & e^{\bar{\beta}J} c_{+} \end{bmatrix}, \qquad (A9)$$

and coefficients

$$c_{\pm} = c(\pm 1), \qquad (A10)$$

obtained from

$$c(S_i) = \sqrt{\frac{1}{\alpha(S_i)}} e^{\frac{\bar{\beta}}{n} \left[ \frac{1}{2\alpha(S_i)} \left(\frac{f}{EA}\right)^2 + \lambda_0(S_i) \frac{f}{EA} - \frac{nq(S_i)}{EA} \right]} \qquad (A11)$$

The rescaled Gibbs free energy is

$$\bar{G} = \frac{G}{EAL} = -\frac{1}{\bar{\beta}} \ln \bar{Z}_G. \qquad (A12)$$

It is thus possible to obtain the relation between the expectation value of the average stretch and the rescaled force as

$$\langle \lambda \rangle = -EA \frac{\partial \bar{G}}{\partial f}. \qquad (A13)$$

## Appendix A.2 A special case

Let us consider the case where $J = 0$ and $\alpha = 1$. Thus, we are neglecting the interface energy and considering equal stiffness for the two phases. Adapting the general expression of the partition function and Gibbs free energy, form Eq. (A13) we find that

$$\langle \lambda \rangle = \frac{f}{EA} + \frac{1 + \chi e^{\frac{\bar{\beta}}{n} \frac{f}{EA} (\chi - 1)} e^{-\frac{\bar{\beta}}{n} \frac{g}{EAL}}}{1 + e^{\frac{\bar{\beta}}{n} \frac{f}{EA} (\chi - 1)} e^{-\frac{\bar{\beta}}{n} \frac{g}{EAL}}}, \qquad (A14)$$

where we defined $q = n\Delta E$. We can rephrase this result in terms of the strain $\langle \epsilon \rangle = (\langle x_n \rangle - L)/L$ and stress $\sigma = f/A$ as

$$\langle \epsilon \rangle = \frac{\sigma}{E} + \frac{\chi - 1}{1 + e^{-\frac{\bar{\beta}}{nE} \left[ \sigma(\chi - 1) - \frac{g}{AL} \right]}}. \qquad (A15)$$

We notice that the trivial limit $n \to \infty$ corresponds to a shift of the usual relation with a correction given by $(\chi - 1)/2$. On the other hand, the same limit, keeping $\bar{\beta}/n$ as a constant when $n \to \infty$, allows to show the temperature effects. This behavior will be even more evident in the general case of applied fixed stretch. It is interesting to rephrase Eq. (A14) to show that it is possible to capture the behaviour of the thermodynamical limit. By a proper rescaling we see that

$$\frac{\langle x_n \rangle}{n} = \frac{f}{k_0} + l + l \frac{\chi - 1}{1 + e^{-\frac{1}{k_B T} [f l (\chi - 1) - \Delta E]}}. \qquad (A16)$$

## Appendix A.3 Fixed extension: Helmholtz ensemble

We consider the case of fixed total elongation of the chain obtained by assigning $\delta = x_n/l = \sum_i \lambda_i = n\lambda$ thus considering the Helmholtz ensemble. As already presented in the main text, the connection between the ensembles is given by the Fourier transform allowing to obtain the canonical partition function in the Helmholtz ensemble, starting from $\bar{Z}_G(f)$. We have

$$\bar{Z}_H(\delta) = \frac{B_n}{2\pi} \int_{-\infty}^{+\infty} \bar{Z}_g(i\omega) e^{-i \frac{\bar{\beta}}{n} \frac{\delta}{EA} \omega} d\omega, \qquad (A17)$$
<figure><img src="image_11.png" /><figcaption>Fig. 9 (a): applied average stretch versus expectation value of rescaled force for $\bar{\beta} = 30$ and $n = \{10, 30, 40, 50, 60\}$. (b): applied average stretch versus expectation value of rescaled force for $n = 30$ and $\bar{\beta} = \{30, 50, 80, 100\}$. In both figures the parameters are $M = 1$, $E = 1$, $A = 1$, $L = 1$, $\alpha = 1$, $\chi = 5$, $J = 0.2$, $q = 3$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10 Applied average stretch versus expectation value of rescaled force for a constant value of $\bar{\beta}/n = 3$ and $\bar{\beta} = \{30, 90, 120, 150\}$. The parameters are $M = 1$, $E = 1$, $A = 1$, $L = 1$, $\alpha = 1$, $\chi = 5$, $J = 0.2$, $q = 3$</figcaption></figure>

where $B_n$ is a constant and we considered the change of variable $f \to i\omega$. A direct calculation, following [21], gives

$$ \bar{Z}_H(\delta) = \frac{A_n}{\cosh \bar{\beta} \bar{J}} \left\{ \sum_{k=0}^{\lfloor \frac{n}{2} \rfloor} \binom{n}{2k} \mathcal{W}_k + e^{-2\bar{\beta} \bar{J}} \sum_{k=0}^{\lfloor \frac{n-1}{2} \rfloor} \binom{n}{2k+1} \mathcal{W}_k \right\}, \quad (A18) $$

where $A_n$ is a constant, non relevant for the subsequent results, and

$$ \mathcal{W}_k = \sum_{j=0}^{k} \sum_{s=0}^{n-2j} \binom{k}{j} \binom{n-2j}{s} e^{-(s+j) \frac{\bar{\beta}}{n} \frac{q}{EAL}} \sqrt{\frac{1}{\alpha^{s+j}}} \sqrt{\frac{1}{n-s-j+\frac{s+j}{\alpha}}} (-1)^j 4^j (1-e^{-4\bar{\beta}\bar{J}})^j \\ \times \exp \left\{ -\frac{\bar{\beta}}{2n} \frac{[\delta - (n-s-j) - \chi(s+j)]^2}{(n-s-j+\frac{s+j}{\alpha})} \right\}. \quad (A19) $$

The rescaled Helmholtz free energy is

$$ \bar{F} = \frac{F}{EAL} = -\frac{1}{\bar{\beta}} \ln \bar{Z}_H. \quad (A20) $$

Thus, the relation between the expectation value of the rescaled force and the average stretch is

$$ \frac{\langle f \rangle}{E A} = \frac{\partial \bar{F}}{\partial \lambda} = n \frac{\partial \bar{F}}{\partial \delta}. \quad (A21) $$
In Figure 9, we show two examples of the force-stretch relation for different values of *n* for constant *β* (a) and the case of *n* constant for different values of $\bar{\beta}$ (b). In Figure 10, we show the force-stretch relation for increasing values of *n* with $\bar{\beta}/n$ kept constant. We observe that the oscillations are suppressed but the initial peak is still present.

References

1. Boltzmann, L.: Theoretical Physics and Philosophical Problems, translated by P. Foulkes (Reidel, Dordrecht/Boston, 1974), p. 169

2. Hilbert, D.: Mathematical Problems, Lecture delivered before the International Congress of Mathematicians at Paris in 1900. Bull. Am. Math. Soc. **8**, 437–479 (1902)

3. Compagner, A.: Thermodynamics as the continuum limit of statistical mechanics. Am. J. Phys. **57**, 106 (1989)

4. Bruin, C., Compagner, A.: A linear lattice gas with variable mesh size subject to gravity. Physica **68**, 171–179 (1973)

5. Uhlenbeck, G. E.: in Fundamental Problems in Statistical Mechanics II. Edited by E. G. D. Cohen (North-Holland, Amsterdam, 1968)

6. Parisi, G.: Statistical Field Theory. Addison-Wesley Publishing Company Inc, New York (1988)

7. Todhunter, I.: A History of the Theory of Elasticity and of the Strength of Materials from Galilei to the Present Time. Cambridge University Press, Cambridge, UK (1893)

8. Torquato, S.: Random heterogeneous materials: Microstructure and macroscopic properties. Springer-Verlag, New York (2002)

9. Sahimi, M.: Heterogeneous materials I, linear transport and optical properties. Springer-Verlag, New York (2003)

10. Sahimi, M.: Heterogeneous materials II, nonlinear and breakdown properties and atomistic modeling. Springer-Verlag, New York (2003)

11. Kanaun, S., Levin, V.: Self-consistent methods for composites. Static problems: Vol. 1. Dordrecht: Springer, (2008)

12. Kanaun, S., Levin, V.: Self-consistent methods for composites. Wave propagation in heterogeneous materials: Vol. 2. Dordrecht: Springer, (2008)

13. Colombo, L., Giordano, S.: Nonlinear elasticity in nanostructured materials. Report on Progress in Physics **74**, 116501 (2011)

14. Braides, A.: Γ-convergence for beginners. Oxford University Press, (2002)

15. Dal Maso, G.: An introduction to Γ-convergence. Birkhäuser, Basel (1993)

16. Francfort, G.A., Murat, F.: Homogenization and optimal bounds in linear elasticity. Arch. Ration. Mech. Anal. **94**, 307–334 (1986)

17. Milton, G.W.: The theory of composites. Cambridge University Press, Cambridge (2004)

18. Tartar, L.: The general theory of homogenization: A personalized introduction. Springer-Verlag, Berlin (2009)

19. Giordano, S.: Order and disorder in heterogeneous material microstructure: electric and elastic characterization of dispersions of pseudo oriented spheroids. Int. J. Eng. Sci. **43**, 1033–1058 (2005)

20. Giordano, S.: Nonlinear effective properties of heterogeneous materials with ellipsoidal microstructure. Mech. Mater. **105**, 16–28 (2017)

21. Cannizzo, A., Bellino, L., Florio, G., Puglisi, G.: Giordano, Thermal control of nucleation and propagation transition stresses in discrete lattices with non-local interactions and non-convex energy. European Physical Journal Plus **137**, 569 (2022)

22. Gibbs, J.W.: Elementary Principles in Statistical Mechanics. Charles Scribner's Sons, New York (1902)

23. Weiner, J.H.: Statistical Mechanics of Elasticity. Dover Publication Inc., New York (2002)

24. Manca, F., Giordano, S., Palla, P.L., Zucca, R., Cleri, F., Colombo, L.: Elasticity of flexible and semiflexible polymers with extensible bonds in the Gibbs and Helmholtz ensembles. J. Chem. Phys. **136**, 154906 (2012)

25. Navier, L.: Sur les lois de l'équilibre et du mouvement des corps solides élastiques. Bulletin des sciences par la Société Philomatique de Paris, 177-181, (in French) (1823)

26. Cauchy, A.: Sur l'équilibre et le mouvement d'un système de points matériels sollicités par des forces d'attraction ou de répulsion mutuelle. Exercices de Mathématiques **3**, 188-212 (1828). (in French)

27. Love, A. E. H.: A Treatise on the Mathematical Theory of Elasticity (Dover Books on Engineering)

28. Lancia, M.R., Vergara Caffarelli, G., Podio-Guidugli, P.: Null lagrangians in linear elasticity. Math. Models Methods Appl. Sci. **05**, 415–427 (1995)

29. Puglisi, G.: Hysteresis in multi-stable lattices with non-local interactions. J. Mech. Phys. Solids **54**, 2060–2088 (2006)

30. Gadomski, A., Łuczka, J., Rudnicki, R.: Finite volume effects in a model grain growth. XXPhys. A **325**, 284–291 (2003)

31. Bialas, P., Spiechowicz, J., Luczka, J.: Quantum analogue of energy equipartition theorem. Journal of Physics A: Mathematical and Theoretical **52**, 15LT01 (2019)

32. Luczka, J.: Quantum Counterpart of Classical Equipartition of Energy. J. Stat. Phys. **179**, 839–845 (2020)

33. Bellino, L., Florio, G., Puglisi, G.: The influence of device handles in single-molecule experiments. Soft Matter **43**, 8680–8690 (2019)

34. Manca, F., Giordano, S., Palla, P.L., Cleri, F.: On the equivalence of thermodynamics ensembles for flexible polymer chains. XXPhys. A **395**, 154–170 (2014)

35. 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. Research **2**, 033227 (2020)

36. 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. and Theor. **54**, 445001 (2021)

37. Skvortsov, A.M., Klushin, L.I., Leermakers, F.A.M.: Negative compressibility and nonequivalence of two statistical ensembles in the escape transition of a polymer chain. J. Chem. Phys. **126**, 024905 (2007)

38. Dimitrov, D.I., Klushin, L.I., Skvortsov, A.M., Milchev, A., Binder, K.: The escape transition of a polymer: A unique case of non-equivalence between statistical ensembles. Eur. Phys. J. E **29**, 9 (2009)
39. Dutta, S., Benetatos, P.: Inequivalence of fixed-force and fixed-extension statistical ensembles for a flexible polymer tethered to a planar substrate. Soft Matter **14**, 6857 (2018)

40. Dutta, S., Benetatos, P.: Statistical ensemble inequivalence for flexible polymers under confinement in various geometries. Soft Matter **16**, 2114 (2020)

41. Skvortsov, A.M., Klushin, L.I., Polotsy, A.A., Binder, K.: Mechanical desorption of a single chain: Unusual aspects of phase coexistence at a first-order transition. Phys. Rev. E **85**, 031803 (2012)

42. Milchev, A., Rostiashvili, V.G., Bhattacharya, S., Vilgis, T.A.: Polymer desorption under pulling a 1st-order phase transition without phase coexistence. Phys. Procedia **3**, 1459 (2010)

43. Ivanov, A., Klushin, L.I., Skvortsov, A.M.: How to understand the ensemble equivalence during stretching of a single macromolecule. Polym. Sci., Ser. A **54**, 602 (2012)

44. Giordano, S.: Spin variable approach for the statistical mechanics of folding and unfolding chains. Soft Matter **13**, 6877–6893 (2017)

45. Benedito, M., Giordano, S.: Isotensional and isometric force-extension response of chains with bistable units and Ising interactions. Phys. Rev. E **98**, 052146 (2018)

46. 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. **149**, 054901 (2018)

47. Florio, G., Puglisi, G.: Unveiling the influence of device stiffness in single macromolecule unfolding. Sci. Rep. **9**, 4997 (2019)

48. Florio, G., Puglisi, G.: A predictive model for the thermomechanical melting transition of double stranded DNA. Acta Biomater. **157**, 225–235 (2023)

49. Giordano, S.: Statistical mechanics of rate-independent stick-slip on a corrugated surface composed of parabolic wells. Continuum Mech. Thermodyn. **34**, 1343–1372 (2022)

50. Giordano, S.: Temperature dependent model for the quasistatic stick-slip process on a soft substrate. Soft Matter **19**, 1813 (2023)

51. Cannizzo, A., Giordano, S.: Thermal effects on fracture and the brittle-to-ductile transition. Phys. Rev. E **107**, 035001 (2023)

52. Binetti, C., Cannizzo, A., Florio, G., Pugno, N.M., Puglisi, G., Giordano, S.: Exploring the impact of thermal fluctuations on continuous models of adhesion. Int. J. Eng. Sci. **208**, 104194 (2025)

53. Gadomski, A., Kruszewska, N.: Matter-Aggregating Low-Dimensional Nanostructures at the Edge of the Classical vs. Quantum Realms, Entropy **25**, 1 (2025)

54. Wigner, E.: On the Quantum Correction For Thermodynamic Equilibrium. Phys. Rev. **40**, 749–759 (1932)

55. Uhlenbeck, G.E., Gropper, L.: The Equation of State of a Non-ideal Einstein-Bose or Fermi-Dirac Gas. Phys. Rev. **41**, 79–90 (1932)

56. Kirkwood, J.G.: Quantum Statistics of Almost Classical Assemblies. Phys. Rev. **44**, 31–37 (1933)

57. Landau, L.D., Lifshitz, E.M.: Statistical Physics. Pergamon Press, Oxford (1980)

58. Abramowitz, M., Stegun, I.A.: Handbook of Mathematical Functions. Dover Publication, New York (1970)

59. Gadomski, A., Karpinski, K.: Erasure by friction: an over-scales-manifesting estimation of the nanoscale (quantum) coefficient of friction. J. Phys. D Appl. Phys. **58**, 135310 (2025)

60. Olver, F.W.J., Lozier, D.W., Boisvert, R.F., Clark, C.W.: NIST Handbook of Mathematical Functions. National Institute of Standards and Technology and Cambridge University Press, New York (2010)

61. Gibbs, J.H., DiMarzio, E.A.: Statistical Mechanics of Helix-Coil Transitions in Biological Macromolecules. J. Chem. Phys. **30**, 271 (1959)

62. Crothers, D.M., Kallenbach, N.R., Zimm, B.H.: The melting transition of low-molecular-weight DNA: Theory and experiment. J. Mol. Biol. **11**, 802 (1965)

63. Kittel, C.: Phase Transition of a Molecular Zipper. Am. J. Phys. **37**, 917 (1969)

64. Lambert, J.H.: Observationes variae in mathesin puram. Acta Helvetica **3**, 128–168 (1958)

65. Corless, R.M., Gonnet, G.H., Hare, D.E.G., Jeffrey, D.J., Knuth, D.E.: On the Lambert W function. Adv. Comput. Math. **5**, 329–359 (1996)

66. De Tommasi, D., Puglisi, G., Saccomandi, G.: Multiscale mechanics of macromolecular materials with unfolding domains. J. Mech. Phys. Sol. **78**, 154–172 (2015)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.