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

HAL
open science

# Thermalization of a driven bi-stable FPU chain

Yalchin Efendiev, Lev Truskinovsky

► To cite this version:

Yalchin Efendiev, Lev Truskinovsky. Thermalization of a driven bi-stable FPU chain. Continuum Mechanics and Thermodynamics, 2010, 22, pp.679-698. 10.1007/s00161-010-0166-5. hal-00547961

HAL Id: hal-00547961

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

Submitted on 25 Feb 2021

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Yalchin R. Efendiev · Lev Truskinovsky

# Thermalization of a driven bi-stable FPU chain

**Abstract** We study Hamiltonian dynamics of a Fermi–Pasta–Ulam (FPU) chain with bi-stable elements. We show, that a quasi-static driving of a ‘cold’ chain beyond the spinodal threshold leads to complex dynamical behavior involving equipartition which suggests thermalization. The subsequent quasi-static cycling between the two energy wells produces reversible temperature oscillations which we link to the release (or absorbtion) of the latent heat. By adopting canonical distribution we obtain a thermodynamical description of the chain which agrees well with numerically computed time-averaged behavior of the corresponding dynamical system.

**Keywords** Fermi–Pasta–Ulam system · Bi-stable elements · Thermolization · Ergodicity · Statistical mechanics · Latent heat

## 1 Introduction

About 60 years ago Fermi et al. [9] asked whether a nonlinear chain with fixed ends and nonzero initial velocities can reach a state of thermalization and over the years this problem attracted a lot of attention (see reviews [2, 4, 10, 11, 13, 25]). It was soon realized that for springs with convex energies the thermalization can be reached only if the initial kinetic energy is above a certain threshold. Below this threshold the nonlinearity is weak and the behavior of the system remains close to the quasi-periodic behavior of the integrable linear system.

The important aspect of the classical FPU nonlinearity was its linearizability. In this paper we consider non-convex springs which implies that there exists at least one state (spinodal) where the linearized model degenerates. In this state one can expect non-integrability even at infinitesimal energy of initial excitations. To support this intuition we provide numerical and analytical evidence that for the springs with double-well energy the state of thermalization or, at least of equipartition, can be reached with zero kinetic energy input. To reach the spinodal state we subjected the chain to infinitely slow driving in a hard loading device. In this case the ‘cold’ system will first accumulate a certain amount of potential energy while remaining in a metastable state, and then release it through a dynamic spinodal instability. Reaching thermalization in this way would then mean that the system was brought exactly to the boundary of the KAM island of regular behavior.

To be more specific, consider a mass-spring chain governed by a dimensionless system of Newton equations. Denote by $x_i(t)$, $i = 0, \dots, N$ the coordinates of the particles and by $e(\epsilon)$ the elastic energy of the springs. Then

Y. R. Efendiev (✉)

Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA

E-mail: efendiev@math.tamu.edu
$$ \delta \frac{d^2 x_i}{dt^2} = \sigma \left( \frac{x_{i+1} - x_i}{\delta} \right) - \sigma \left( \frac{x_i - x_{i-1}}{\delta} \right) \qquad (1) $$

where $\sigma(\epsilon) = e'(\epsilon)$. We assume zero initial velocities and apply the following time-dependent boundary conditions:

$$ x_0 = 0, \quad x_N(t) = x_N^0 + vt. \qquad (2) $$

The system (1) depends on the two dimensionless parameters. The first parameter, measuring the degree of discreteness of the system, is $\delta = h/L = N^{-1}$, where h is the interatomic distance and L is the macroscopic length of the system. The second parameter, measuring the slowness of the driving, is $v = V/(\sqrt{E/\rho})$, where $\rho$ is mass density, E is the characteristic elastic modulus and V is the velocity of driving. We are interested in the behavior of the system (1) at $\delta \to 0$ (thermodynamic limit) and $v \to 0$ (quasi-static driving).

Non-convex energies at a micro-scale appear often in the modeling of materials and distributed structures with internal instabilities. A typical example is a shape memory alloy where the presence of a phase transition evokes description in terms multi-well snap-springs [16]. The energy landscape for such chains was found to be extremely rugged (e.g. [18]). First numerical experiments with a fully inertial bi-stable chain [1] showed that the system becomes highly oscillatory after reaching the spinodal threshold. The use of a micro-canonical ensemble for a thermodynamical description of a minimal chain with two snap-springs was rationalized in [5]. A symmetric bi-stable chain in a soft device was studied numerically in [26] and a canonical distribution was found to be adequate for the prediction of the macroscopic behavior. Different numerical methods of finding equilibrium properties for generic 1D lattices were recently reviewed in [3]. Other contributions to the thermalization problem in discrete systems with non-convexity can be found in [7, 12, 14, 15, 24].

In the present work, most of which was completed in 2001 (see [26]), we study the emergence of ergodicity in an asymmetric bi-stable FPU chain driven quasi-statically in a hard device. In our numerical experiment we first bring the ‘cold’ system to the point of instability and study the equipartition of the energy. We then perform a cyclic quasi-static loading of the ‘self-heated’ chain through the phase transition. This leads to reversible and well reproducible temperature oscillations which we associate with the release (or absorbtion) of the latent heat. To justify the ergodicity hypothesis we compare the time-averaged behavior of the system with the statistical computations based on canonical distribution. For the case of a bi-quadratic springs we obtain explicit expression for the free energy and show that the time-averaged behavior of the chain agrees almost perfectly with the computed adiabats. We then study the emergence of hysteretic behavior in the damped chain; the overdamped limit was analyzed in [19].

The paper is organized as follows: In Sect. 2 we report the results of our numerical experiments, demonstrate the equipartition of the energy, and clarify distinction between reversible and irreversible heating. In Sect. 3 we use the classical methods of statistical mechanics to derive the thermodynamic response of the system. We show that in case of bi-quadratic energy the statistical problem can be reduced to a 1D Ising problem for a special diamagnetic lattice. In Sect. 4 we compare the results of the thermodynamic analysis with direct numerical simulations. The final Sect. 5 contains a brief discussion of the under-damped dynamics.

## 2 Numerical studies

Assume that $\delta = 1$ and denote by $\epsilon_j = x_{j+1} - x_j$ the strain in the jth spring, $j = 0, \dots, N-1$. To make analytical computations possible, we assume that the function $e(\epsilon)$ is bi-quadratic

$$ e(\epsilon) = \begin{cases} \frac{\mu_1}{2} (\epsilon - a_1)^2 & \text{if } \epsilon \le a_0 \\ \frac{\mu_2}{2} (\epsilon - a_2)^2 & \text{if } \epsilon > a_0, \end{cases} \qquad (3) $$

where the spinodal strain is $a_0 = (1 - \kappa)/(1 + \kappa)$, $\kappa = \sqrt{\mu_2/\mu_1}$ (see Fig. 1). While one of the parameters $\mu_1$ or $\mu_2$ can be chosen to be equal to 1, we preserve both of them for the symmetry of the analytical expressions. We say that the interval $\epsilon < a_0$ belongs to Phase I and the interval $\epsilon > a_0$—to Phase II. The spinodal region is represented in this model by a single-point $\epsilon = a_0$.

We first prepare a “cold” chain with all springs initially at rest (in phase I) and then drive it quasi-statically with a constant strain rate $v$. In the actual numerical experiments we diminished the value of $v$ until the averaged behavior stabilized. For numerical integration of the system (1) we used a symplectic numerical method widely known as Stormer–Verlet algorithm (see [20] for more details).
<figure><img src="image_3.png" /><figcaption>Fig. 1 Stress–strain and energy–strain relations for a bi-linear snap-spring. Available energy is released during spinodal instability</figcaption></figure>

To quantify the observed phenomena, we need to introduce a measure of intensity of the micro-oscillations. Consider, for instance, the stage of the stretching process where all springs oscillate in a single quadratic energy well. Due to the linearity of the system we can decompose the motion into the average and the fluctuations. The average deformation is both homogeneous and quasi-static

$$ \overline{x_i(t)} = (i/N) (x_N^0 + vt). $$

The general motion can then be represented in the form

$$ x_i(t) = \overline{x_i(t)} + \omega_i(t), $$

where $\omega_i(t)$ is the oscillatory contribution with zero average over a time interval which is much larger than the period of the microscopic oscillations ~1 and much smaller than the characteristic time of the loading ~ $v^{-1}$. Next, we can decompose the oscillatory motion into normal modes in the quadratic regions. Recall that $\omega_k(t)$ is the solution of the linear system

$$ \frac{d^2}{dt^2} \omega(t) = A\omega, $$

where A is a tri-diagonal matrix. To define the modes we need to diagonalize the matrix A, which can be done explicitly (e.g. [2]). If we assume $\mu = 1$ (otherwise we can re-scale the time) we obtain that the eigenvalues of A are $\lambda_l = -2 + 2 \cos(\frac{l\pi}{N+1})$ and the corresponding eigenvectors are $(\sin(\frac{l\pi}{N+1}), \sin(\frac{2l\pi}{N+1}), \dots, \sin(\frac{ln\pi}{N+1}))^t$. Once the eigenvectors and eigenvalues are known, the matrix A can be represented as

$$ D = U^T A U, \qquad (4) $$

where D is a diagonal matrix consisting of eigenvalues and the columns of U are corresponding eigenvectors. The system can be represented as a set of noninteracting harmonic oscillators

$$ \frac{d}{dt^2} \phi = D\phi, $$

where $\phi = U\omega$. The energy of each mode (each oscillator) is then

$$ E_l = \frac{1}{2} \left( \left( \frac{d}{dt} \phi_l \right)^2 + \lambda_l^2 \phi_l^2 \right). $$

Consider now the symmetric bi-quadratic energy with $\mu_1 = \mu_2 = \mu$. In this case the definition of the modes can be formally extended into the regime where the springs change phase and where oscillators are
<figure><img src="image_4.png" /><figcaption>Fig. 2 Total energy of the modes for the system of five springs during quasistatic stretching of an initially ‘cold’ chain</figcaption></figure>

no longer independent. In Fig. 2 we show the evolution of the total energy of the modes for the system of five springs which is plotted as a function of the total elongation. As we have already mentioned, the chain is initially “cold” and the energy of the modes is equal to zero. Before the system reaches the spinodal point the springs remain in Phase I and the energy of oscillations remains equal to zero which is expected in view of the quasi-static nature of the driving. At large elongations all springs are again oscillating in a single well (Phase II) and the energy of oscillations again remains constant because the average motion decouples. At this stage the system regains its linearity and can be represented by a set of noninteracting harmonic oscillators with a randomness in initial conditions inherited from the mixed-phase dynamics. As the chain enters the nonlinear, interphase regime the springs start to oscillate between phases I and II. In this regime the micro and macro motion become coupled and, in particular, during the spinodal instability a finite amount of energy is converted from the average motion into the oscillatory motion. This energy is marked as ‘available energy’ in Fig. 1. One half of the available energy goes into kinetic energy of oscillations and another half goes into potential energy.

Observe that due to simplicity of the model where only the nearest neighbors interact and the elastic energy is permutationally invariant we do not see phase boundaries propagating in an organized fashion. Adding linear interaction of next to nearest neighbors of the ferromagnetic type (see [22]) fixes this problem and localizes transformation in several phase boundaries, whose dynamics is expected to remain complex.

The next experiment is intended to show the details of the self-heating process. We prepare a ‘cold’ chain with $\epsilon_i = a_0$ (for all $i$) and then we slightly perturb one element. Note that we do not load the chain in this case. In Fig. 3 we show the ensuing motion. One can see that after the regular nucleation stage described in [1,17,21] the motion becomes increasingly complex due to multiple interaction of the waves reflected from the boundaries. We demonstrate the equipartition by plotting in Fig. 4 the time dependence of the averaged kinetic energies of the particles

$$ \frac{1}{2}kT_i = \frac{1}{\tau} \int_{t_0}^{t_0+\tau} \left( \frac{dx_i(t)}{dt} \right)^2 dt, $$

where $k$ is a dimensional (Boltzmann) constant and $\tau$ is the averaging window. One can see that after a finite time the system achieves the state of equipartition with

$$ T_1 = \dots = T_{N-1} = T. $$

The same conclusion can be reached when one considers the energy of the modes. Therefore, our numerical experiment indicates an important sign of ergodicity of the infinitesimally perturbed system with regular behavior.
<figure><img src="image_5.png" /><figcaption>Fig. 3 The motion of the chain with 80 particles after spinodal instability</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 Kinetic energy of the mass points in the chain with eight springs. The ‘cold’ system was perturbed near the threshold of spinodal instability</figcaption></figure>

We can now study the time evolution of the temperature during several cycles of loading and unloading through the phase transformation (several subsequent elongations and contractions). For the case of equal moduli of the phases the total energy of the linear modes is shown in Fig. 5. Note that in both phases, the coordinate transformation *U* (see (4)) is the same and this transformation is used in the mixed region to compute the total energy of the linear modes. Notice again that after the initial irreversible “self heating” the system does not heat or cool any more, and the total temperature of the modes remains practically constant except in the mixed-phase regions. This is consistent with an assumption that the quasi-static driving in the hard device can be represented at the macroscale as the adiabatic evolution of the quasi-static thermodynamical system. The continuity of the temperature in the pure phases can then be viewed as the consequence of the fact that the entropies of the phases are equal. In mechanical terms one can argue that in the “hot” chain instead of massive spinodal instability we see fluctuations-induced tunneling through the barrier. This process does not require
<figure><img src="image_7.png" /><figcaption>Fig. 5 Total energy of linear modes versus time for the system with $\mu_1 = \mu_2 = 1$</figcaption></figure>

overstepping and is not accompanied by the irreversible macroscopic energy conversion into microscopic vibrations.

To confirm our thermodynamical interpretations, we studied the case when the elastic moduli in the two energy wells are different. Then the phase with more shallow energy minimum will have a higher entropy than the phase with the more narrow minimum, and transition from former to latter in adiabatic conditions would require a decrease of temperature. Indeed, the experiments for the chains with $\mu_1 \neq \mu_2$ show that the value of the temperature (represented by the total energy of the modes defined in each of the phase) oscillates in the predicted way as one switch between the phases (see Fig. 6). More precisely, our numerical experiments indicate that the ratio $E^{\mu_1}/E^{\mu_2} \approx \sqrt{\mu_1/\mu_2}$, where $E^{\mu_i}$ is the total energy of the modes in the single phase stage of deformation with stiffness $\mu_i$, ($i = 1, 2$). In terms of the temperature of the particles this relation translates into

$$T^{\mu_1}/T^{\mu_2} \approx \sqrt{\mu_1/\mu_2}.$$

In the next section we show that this numerical result is fully supported by the computations of the time-averaged behavior based on ergodicity hypothesis.

To check the ergodicity of the system one would like to compare statistical mechanical computations with the macroscopic stress–strain and energy–strain relations obtained from the direct numerical simulations for the self-heated chain driven quasi-statically through the phase transition. To perform time averaging and filter oscillations, we used the window size equal to about 30 oscillations. The results for the system with symmetric energy wells are presented in Fig. 7, where we show the macroscopic stress and potential energy for the system of 15 springs. According to our assumption those curves describe the adiabats. For comparison, we plot them against the corresponding curves obtained by the global minimization of the energy (e.g. [8]). Notice that in the two-phase region, where the purely mechanical description predicts zero elastic modulus, our numerical data suggest nonzero “entropic” elasticity. The numerically computed macroscopic potential energy is close to being convex (as can be expected in a 1D system) and is always above the energy of the mechanical ground state due to the potential energy of oscillations (half of the ‘available energy’). The goal of the statistical theory presented in the next section is to qualitatively reproduce these numerical results.

## 3 Statistical modeling

Under the assumption of ergodicity at the microscale one can use the techniques of statistical mechanics to compute the macroscopic constitutive response of the chain. For our finite ‘self-heated’ system the exact nature of the invariant measure remains to be established; however, the presence of equipartition suggests the applicability of the canonical distribution at least when the number of elements is sufficiently large (e.g [2]).
<figure><img src="image_8.png" /><figcaption>Fig. 6 Time dependence of the total energy of the linear modes for a ‘self-heated’ chain driven cyclically (and quasi-statically) through the phase transition. Left figure $\mu_1 = 9$, $\mu_2 = 1$. Right figure $\mu_1 = 0.25$, $\mu_2 = 1$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 Macroscopic behavior of a quasi-statically driven chain with 15 springs. Left figure solid line corresponds to the global energy minimum response, dotted line time averaging of the dynamical behavior. Right figure solid line corresponds to elastic energy in the global minimum, dotted line shows elastic energy obtained from direct numerical simulations</figcaption></figure>

To define the ensemble we fix the total strain $\Delta = \sum_{i=1}^{N} \epsilon_i$ and denote the mean strain by $\bar{\epsilon} = \Delta/N$. We introduce a canonical measure by assuming that a microstate with the energy $E = \sum_{i=1}^{N} \frac{1}{2} p_i^2 + e(\epsilon_i)$ has a probability $f(p_i, \epsilon_i) = \frac{1}{Z} \exp(-\beta E)$, where $\beta = (kT)^{-1}$ is the new macroscopic parameter ($k$ is the Boltzmann constant). The statistical sum $Z$ is defined as a normalization of the probability distribution

$$Z(\bar{\epsilon}, \beta) = \int \exp(-\beta E) d p_1, \dots, d p_{N-1} d\epsilon_1, \dots, d\epsilon_{N-1}$$

and the integration is carried over all micro-configurations characterized by $N-1$ strains $\epsilon_i$ and $N-1$ momenta $p_i$ (we assume that $p_0 = 0$ and $p_N = 0$).

The partition function can be immediately factorized into kinetic and configurational parts $Z = Z_k Z_p$. Due to the independence of momenta, the computation of the kinetic energy contribution is straightforward
$Z_k(\beta) = (2\pi/\beta)^{N-1/2}$. The configurational contribution $Z_p$ is given by

$$Z_p(\bar{\epsilon}, \beta) = \int_{\epsilon_1, \dots, \epsilon_N, \sum_i \epsilon_i = N\bar{\epsilon}} \exp\left(-\frac{1}{2}\beta E_p\right) d\epsilon_1, \dots, \epsilon_N. \qquad (5)$$

where $E_p = \sum_{i=1}^N e(\epsilon_i)$ is the total potential energy of the system.

After the statistical sum is known, the free energy can be computed as follows:

$$F(\bar{\epsilon}, \beta) = -\frac{1}{\beta} \log(Z) = F_k + F_p, \qquad (6)$$

where

$$F_k = -\frac{1}{\beta} \frac{(N-1)}{2} \log\left(\frac{(2\pi)}{\beta}\right), \quad F_p = -\frac{1}{\beta} \log(Z_p).$$

The expression for the entropy takes the form

$$S(\bar{\epsilon}, \beta) = k \frac{N-1}{2} \left(1 + \log \frac{2\pi}{\beta}\right) + k \left(\log(Z_p) - \beta \frac{\partial \log(Z_p(\beta, N))}{\partial \beta}\right). \qquad (7)$$

## 3.1 Single phase system

To illustrate the above formulas we first consider quadratic potential with $e(\epsilon) = \frac{1}{2}\mu\epsilon^2$. In this case the problem reduces to the computation of a Gaussian integral

$$Z_p = \int_{\epsilon_1, \dots, \epsilon_N, \sum_i \epsilon_i = N\bar{\epsilon}} \exp\left(-\frac{1}{2}\beta\mu\sum_i \epsilon_i^2\right) d\epsilon_1, \dots, \epsilon_N. \qquad (8)$$

As is common in such calculations we introduce the new variables $r_0 = 0$, $r_k = \sum_{i=1}^k \epsilon_i$, and then use the known formulas

$$\int_{r_i} \exp\left(-\beta \frac{\mu}{2} \sum_{i=1}^{N} (r_i - r_{i-1})^2\right) dr_1, \dots, dr_N = N^{1/2} \left(\frac{2\pi}{\beta\mu}\right)^{\frac{N-1}{2}} \int_{r_i} \exp\left(-\beta \frac{\mu}{2N} (r_N - r_0)^2\right) dr_1, \dots, dr_N. \qquad (9)$$

to obtain

$$Z_p = N^{1/2} \left(\frac{2\pi}{\beta\mu}\right)^{\frac{N-1}{2}} \exp\left(-\frac{1}{2}\beta\mu N\bar{\epsilon}^2\right).$$

The average potential energy takes the form $\bar{E}_p = (N-1)\bar{\epsilon} + \frac{1}{2}N\mu\bar{\epsilon}^2$. At large N the energy density $\bar{E}_p/N$ splits into the macroscopic elastic energy $\frac{1}{2}\mu\epsilon^2$ and the microscopic energy $1/(2\beta) = kT/2$. The microscopic energy is due to fluctuations which behave as an ideal gas. The specific heat at constant strain is equal to $c_{\bar{\epsilon}} = (\partial\bar{E}/\partial T)_{\bar{\epsilon}} = kN = R$, with half of it due to kinetic energy of oscillations and half due to their potential energy. For the free energy function we obtain

$$F_p(\bar{\epsilon}, \beta) = \frac{1}{(2\beta)} \log(N) - \frac{N-1}{(2\beta)} \log\left(\frac{2\pi}{(\beta\mu)}\right) + \frac{1}{2} N \mu \bar{\epsilon}^2.$$

One can see that the stress–strain relation at constant temperature

$$\sigma = \partial(F_p/N)/\partial\bar{\epsilon} = \mu\epsilon$$
remains purely mechanical. The expression for the finite size entropy takes the form

$$S = k \left( \frac{(N-1)}{2} \log \frac{(4\pi^2)}{(\beta^2\mu)} - \frac{1}{2} \log(N) + N - 1 \right). \quad (10)$$

We can use the last formula to explain the equilibrium temperature jump as the “hot” asymmetric chain is driven quasi-statically through the phase transition. Indeed, due to adiabatic nature of the loading the entropies of the two limiting linear regimes must be identical. If the corresponding elastic moduli are different the equilibrium temperatures will also be different due to the release or absorption of the latent heat. From (10) we obtain

$$\beta_1^2 \mu_1 = \beta_2^2 \mu_2.$$

which agrees with the relation between the limiting temperatures found in our numerical experiments.

## 3.2 Two-phase system

For the two-phase system with a bi-quadratic energy we can write

$$E_p = \sum_{i=1}^{N} \frac{\mu_i}{2} (\epsilon_i + \alpha_i)^2.$$

Here the discrete parameter $\alpha$ which locates the bottoms of the energy wells takes only two values $-a_1$ and $-a_2$. Similarly, $\mu_i$ takes values $\mu_1$ or $\mu_2$ depending on the value of $\epsilon_i$. Note that both $\mu_i$ and $\alpha_i$ depend on $\epsilon_i$.

For the sake of analytical transparency we can start with a simplified model where we expand the quadratic energy wells beyond the spinodal point into the whole real line and neglect the nonlinear switching condition at $\epsilon_i = \epsilon_0$. Such approximation implying the multivalued energy can be justified in the limit of small temperatures where the contribution from the non single-valued response is negligible. Our numerical results confirm that the simplified model is accurate for the values of temperature generated during the ‘self-heating’ of the ‘cold’ chain loaded quasi-statically through the phase transition.

More specifically, the simplified model adopts the following approximate expression for the configurational partition function:

$$Z_p(\bar{\epsilon}, \beta) = \sum_{\alpha_i} \int_{\epsilon_1, \dots, \epsilon_N, \sum_i \epsilon_i = N\bar{\epsilon}} \exp\left(-\frac{1}{2}\beta E_p\right) d\epsilon_1, \dots, \epsilon_N. \quad (11)$$

To compute the integrals, we can again introduce the variables $r_k$

$$r_0 = 0, \quad r_k = \sum_{i=1}^{k} \epsilon_i + \left( \sum_{i=1}^{k} \alpha_i \right). \quad (12)$$

In terms of these new variables, the function $Z_p$ can be written in the form

$$Z_p = \sum_{\alpha_i} \int_{r_1} \dots \int_{r_{N-1}} dr_1, \dots, dr_{N-1} \exp \left( -\frac{\beta}{2} \sum_{i=1}^{N} \mu_i (r_i - r_{i-1})^2 \right). \quad (13)$$

To proceed further, we denote

$$P(\mu) = \begin{bmatrix} \mu_1 + \mu_2 & \mu_2 & 0 & \dots & 0 \\ \mu_2 & \mu_2 + \mu_3 & \mu_3 & \dots & 0 \\ 0 & \mu_3 & \mu_3 + \mu_4 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & \mu_{N-1} & \mu_{N-1} + \mu_N \end{bmatrix}, \\ \phi(\mu) = (0, \dots, 0, \mu_N r_N). \quad (14)$$
Then (13) can be rewritten as

$$Z_p = \sum_{\alpha_i} \int_{r_1} \dots \int_{r_{N-1}} \exp \left( -\frac{\beta}{2} \left( (P(\mu)r, r) + \beta(\phi(\mu), r) \right) \right) dr_1, \dots, dr_{N-1} \exp \left( -\frac{\beta}{2} \mu_N r_N^2 \right), \quad (15)$$

where $(\cdot, \cdot)$ denotes scalar product and $r = (r_1, \dots, r_{N-1})$. Next, we use again the well known property of Gaussian integrals

$$\int \exp \left( -\frac{1}{2} \sum_{i,j} P_{i,j} x_i x_j + \sum_k \phi_k x_k \right) dx_1, \dots, dx_N = \frac{(2\pi)^{N/2}}{\sqrt{\det(P)}} \exp \left( \frac{1}{2} \sum_{i,j} P_{i,j}^{-1} \phi_i \phi_j \right), \quad (16)$$

to obtain

$$Z_p = \sum_{\alpha_i} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}} \sqrt{\det(P)}} \exp \left( \frac{1}{2} \beta (\mu_N P_{N-1,N-1}^{-1} - 1) \mu_N r_N^2 \right), \quad (17)$$

where $P_{N-1,N-1}^{-1}$ is the last entry (last column and last row) of the inverse of $P(\mu)$. Now, noting that $r_N = \Delta + \sum_i \alpha_i$, the expression for $Z_p$ can be rewritten as

$$Z_p = \sum_{\alpha_i} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}} \sqrt{\det(P)}} \exp \left( \frac{1}{2} \beta (\mu_N P_{N-1,N-1}^{-1} - 1) \mu_N (\Delta + \sum_i \alpha_i)^2 \right). \quad (18)$$

Every term that involves $P$ depends on $\mu_i$ and consequently on $\alpha_i$. Thus $\det(P)$ can be evaluated explicitly $\det(P) = \sum_{i=1}^N D/\mu_i$, where $D = \mu_1 \dots \mu_N$. The term $P_{N-1,N-1}^{-1}$ can be written in the closed form (see [6])

$$P_{N-1,N-1}^{-1} = \mu_{N-1}\theta_{N-1}/\theta_N,$$

where $\theta_i = (\mu_i + \mu_{i+1})\theta_{i-1} - \mu_{i-1}^2\theta_{i-2}$, $\theta_0 = 1$, $\theta_{-1} = 0$. Denoting $\eta_i = \theta_i/\theta_{i-1}$, we obtain the following recursive relation

$$\eta_i = \frac{1}{\mu_i + \mu_{i+1} - \mu_{i-1}^2 \eta_{i-1}}, \quad \eta_1 = \frac{1}{\mu_1 + \mu_2}$$

We then note that

$$P_{N-1,N-1} = \mu_{N-1}\eta_N,$$

and by induction argument show that

$$P_{N-1,N-1} = \frac{\prod_{i=1}^{N-1} \mu_i \sum_{i=1}^{N-1} \frac{1}{\mu_i}}{\prod_{i=1}^{N} \mu_i \sum_{i=1}^{N} \frac{1}{\mu_i}}.$$

Using these expressions we finally obtain

$$Z_p = \sum_{\alpha_i} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}} \left( \prod_{i=1}^{N} \mu_i \sum_{i=1}^{N} \frac{1}{\mu_i} \right)^{1/2}} \exp \left( -\frac{1}{2} \beta \left( \sum_{i=1}^{N} \frac{1}{\mu_i} \right)^{-1} (\Delta + \sum_i \alpha_i)^2 \right). \quad (19)$$

Notice that we were able to eliminate all linear variables (strains) and reduce the problem to solving a particular (diamagnetic) one-dimensional Ising model. The mean field nature of this model is due to absence of interactions beyond nearest neighbors, when the discrete statistical sum (19) is invariant with respect to the permutation of $\mu_i$ and can be rewritten in terms of the number of springs in the first phase $p$. First we observe that

$$\prod_{i=1}^{N} \mu_i = \mu_1^p \mu_2^{N-p}, \quad \sum_{i=1}^{N} \frac{1}{\mu_i} = \frac{p}{\mu_1} + \frac{N-p}{\mu_2}, \quad \sum \alpha_i = pa_1 + (N-p)a_2.$$
Then, we can write the following explicit expression capturing the size effect:

$$Z_p(\bar{\epsilon}, \beta) = N^{-1/2} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}}} \sum_{p=0}^{N} \frac{N!}{p!(N-p)!} (\mu_1^{p/N} \mu_2^{1-p/N})^{-N/2} \left( \frac{p/N}{\mu_1} + \frac{1-p/N}{\mu_2} \right)^{-1/2} \\ \times \exp \left( -\frac{\beta}{2} N \left( \frac{p/N}{\mu_1} + \frac{1-p/N}{\mu_2} \right)^{-1} \left( \bar{\epsilon} + \frac{p}{N} a_1 + \left( 1 - \frac{p}{N} \right) a_2 \right)^2 \right). \quad (20)$$

### 3.3 Mechanical limit

To make connection with the purely mechanical theory based on the global energy minimization of the elastic energy of the chain (ground state response) we need to consider the limit $\beta \to \infty$ when the simplified model becomes asymptotically exact. In this limit the sum (20) can be computed analytically using Laplace method. More precisely, we can use the fact that in the limit of large $\beta$ one can approximate the partition function $Z_p$ by

$$\int \exp(-\beta P(\mathbf{x})) d\mathbf{x} = \exp(-\beta P(\mathbf{x}_0)) \sqrt{\frac{(2\pi)^N}{\det(A)}}, \quad (21)$$

where the matrix $A$ is defined by

$$A_{ij} = -\beta \frac{\partial^2}{\partial x_i \partial x_j} P(\mathbf{x})|_{\mathbf{x}=\mathbf{x}_0}.$$

and $\mathbf{x}_0$ is the point of the minimum of $P(\mathbf{x})$. Now, suppose for simplicity that we are dealing with equal moduli case and assume that $E(\epsilon_1, \dots, \epsilon_N)$ has minimum $E_{\min}$. In the neighborhood of the minimum $E$ has quadratic behavior with

$$\frac{\partial^2}{\partial \epsilon_i \partial \epsilon_j} E = \mu \delta^{ij},$$

therefore, by using the Laplace formula (21) we obtain

$$Z_p = \left(\frac{2\pi}{\beta\mu}\right)^{\frac{N}{2}} \sum_{\alpha_i} \exp(-\beta E_{\min}). \quad (22)$$

By using the explicit expression for the ground state configurations [8, 18], we can rewrite (22) in the form

$$Z_p = \left(\frac{2\pi}{\beta\mu}\right)^{\frac{N}{2}} \sum_{\alpha_i} \exp \left( -\beta \frac{\mu}{2N} \left( \Delta + \sum_{i=1}^{N} \alpha_i \right)^2 \right).$$

In terms of the average strain $\bar{\epsilon} = \Delta/N$ we obtain

$$Z_p(\bar{\epsilon}, \beta) = \left(\frac{2\pi}{\beta\mu}\right)^{\frac{N}{2}} \sum_{\alpha_i} \exp\left(-\frac{N\beta\mu}{2}(\epsilon + (2p/N - 1))^2\right). \quad (23)$$

The corresponding stress-strain relation reads

$$\sigma = \mu\epsilon - \mu(1 - 2p/N),$$

where $p$ is a piece wise constant function of $\Delta$ (see [8, 18] for more detail).
### 3.4 Thermodynamic limit

To study the limiting case $N \to \infty$, we introduce the new variable $x = p/N$ and write

$$Z_p = N^{-1/2} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}}} \sum_{x=0}^{1} \left(\frac{x}{\mu_1} + \frac{1-x}{\mu_2}\right)^{-1/2} \times \exp\left(-N\left(\frac{\beta}{2}\left(\frac{x}{\mu_1} + \frac{1-x}{\mu_2}\right)^{-1}(\bar{\epsilon} + xa_1 + (1-x)a_2)^2 + x\log(x) + (1-x)\log(1-x) + \frac{x}{2}\log(\mu_1) + \frac{1-x}{2}\log(\mu_2)\right)\right). \quad (24)$$

If we now denote

$$\psi(x) = \frac{\beta}{2} \left( \frac{x}{\mu_1} + \frac{1-x}{\mu_2} \right)^{-1} (\bar{\epsilon} + xa_1 + (1-x)a_2)^2 + x \log(x) + (1-x) \log(1-x) + \frac{x}{2} \log(\mu_1) + \frac{1-x}{2} \log(\mu_2).$$

and replace the sum by an integral, we obtain

$$Z_p \approx N^{1/2} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}}} \int_0^1 \left(\frac{x}{\mu_1} + \frac{1-x}{\mu_2}\right)^{-1/2} \exp(-N\psi(x))dx. \quad (25)$$

Since $N$ is large we can use Laplace method which gives

$$Z_p = \frac{(2\pi)^{\frac{N}{2}}}{\beta^{\frac{N-1}{2}}} \left( \frac{1}{\psi''(x^*)} \right)^{\frac{1}{2}} \left( \frac{x^*}{\mu_1} + \frac{1-x^*}{\mu_2} \right)^{-1/2} \exp(-N\psi(x^*)) \quad (26)$$

where $x^*$ is a minimum of $\psi(x)$.

For large $|\bar{\epsilon}|$, the minimum of $\psi(x)$ is either at $x = 0$ or at $x = 1$. More precisely, when $\psi(x)$ is increasing we have

$$Z_p \approx N^{1/2} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}}} \frac{1}{N\psi'(0)} \mu_2^{(N-1)/2} \exp\left(-N\frac{\beta}{2}\mu_2(\bar{\epsilon} + a_2)^2\right), \quad (27)$$

while for decreasing $\psi(x)$, we obtain

$$Z_p \approx N^{1/2} \frac{(2\pi)^{\frac{N-1}{2}}}{\beta^{\frac{N-1}{2}}} \frac{1}{N\psi'(1)} \mu_1^{(N-1)/2} \exp\left(-N\frac{\beta}{2}\mu_1(\bar{\epsilon} + a_1)^2\right).$$

Here $\psi'(0)$ and $\psi'(1)$ need to be determined by computing the corresponding derivatives with respect to $\bar{\epsilon}$ (at constant $\beta$). One can see that we recover the thermodynamics of the single phases.

For the intermediate values of $\bar{\epsilon}$, we can write

$$\bar{\sigma} = \frac{1}{\beta} \frac{\partial \psi(x^*)}{\partial \bar{\epsilon}}|_{\beta} = \left( \frac{x^*}{\mu_1} + \frac{1-x^*}{\mu_2} \right)^{-1} (\bar{\epsilon} + x^*a_1 + (1-x^*)a_2) + \frac{1}{\beta} \frac{\partial \psi}{\partial x^*} \frac{\partial x^*}{\partial \bar{\epsilon}}.$$

Denote

$$\Psi = \frac{\partial \psi}{\partial x} = \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \frac{\beta}{(x/\mu_1 + (1-x)/\mu_2)^2} (\epsilon + xa_1 + (1-x)a_2)^2 + 2 \frac{\beta}{(x/\mu_1 + (1-x)/\mu_2)} (\epsilon + xa_1 + (1-x)a_2) + \log(x) - \log(1-x) + \frac{1}{2}\log(\mu_1) - \frac{1}{2}\log(\mu_2). \quad (28)$$
If $x^*(\bar{\epsilon})$ is the solution of $\Psi(x, \bar{\epsilon}) = 0$, we obtain

$$ \bar{\sigma} = \left( \frac{x^*}{\mu_1} + \frac{1-x^*}{\mu_2} \right)^{-1} (\bar{\epsilon} + x^*a_1 + (1-x^*)a_2). \quad (29) $$

Using (29) we can rewrite equation for $x^*$ in the following way:

$$ \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \log(x^*) - \log(1-x^*) + \frac{1}{2}\log(\mu_1) - \frac{1}{2}\log(\mu_2) = 0. \quad (30) $$

From here we obtain

$$ x^* = \frac{1}{1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left(\frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma}\right)}. $$

This allows us to write the expression for the function $\bar{\epsilon}(\bar{\sigma})$ in the closed form

$$ \bar{\epsilon} = \bar{\sigma} \left( \frac{1}{\mu_1 \left( 1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right) \right)} + \frac{1}{\mu_2 \left( 1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right) \right)} \right) \\ - \left( \frac{1}{1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right)} + \frac{\sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right)}{1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right)} a_1 + \frac{\sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right)}{1 + \sqrt{\frac{\mu_1}{\mu_2}} \exp\left( \frac{\mu_1 - \mu_2}{2\mu_1\mu_2} \beta\bar{\sigma}^2 + 2\beta\bar{\sigma} \right)} a_2 \right). \quad (31) $$

In the case of symmetric double well energy with equal moduli $\mu_1 = \mu_2 = 1$ and $\alpha_1 = -1, \alpha_2 = 1$ we obtain

$$ x^* = \frac{1}{1 + \exp(2\beta\bar{\sigma})}, $$

and

$$ \bar{\sigma} = \bar{\epsilon} + \frac{2}{1 + \exp(2\beta\bar{\sigma})} - 1. $$

We now illustrate the obtained formulas. For the case of equal moduli, we take $a_1 = 1$ and $a_2 = -1$, and for the case of different moduli, we take $a_1 = 1/\sqrt{\mu_1}$ and $a_2 = -1/\sqrt{\mu_2}$. In all cases, the discontinuity point for the stress-strain relationship is at $\epsilon = 0$. In Fig. 8, we plot the isothermal stress-strain curves at different values of $\beta$. At large $\beta$ (low temperatures), our results are close to those obtained through global energy minimization. At finite temperatures the isotherms exhibit finite slope in the transition region, which is a manifestation of “rubber elasticity”. The finite slope is expected because in 1D the Ising model is known not to show phase transition in the thermodynamic limit.

### 3.5 Adiabats in the thermodynamic limit

The general expression for the entropy is given in (7). For large $N$, we can write

$$ \log(Z_p) \approx \frac{N}{2} \log(2\pi) - \frac{N-1}{2} \log(\beta) - N\psi(x^*), $$

and

$$ \beta \frac{\partial \log(Z)}{\partial \beta} \approx -\frac{N-1}{2} - N\beta \frac{\partial \psi(x^*)}{\partial \beta}. $$

Noting that $\partial\psi(x^*)/\partial x^* = 0$, we obtain

$$ \frac{\partial \psi(x^*)}{\partial \beta} = \frac{1}{2} \left( \frac{x^*}{\mu_1} + \frac{1-x^*}{\mu_2} \right)^{-1} (\bar{\epsilon} + x^*a_1 + (1-x^*)a_2)^2. $$
<figure><img src="image_10.png" /><figcaption>Fig. 8 Stress–strain isotherms for different values of β. Left $μ_1 = μ_2$. Right $μ_2 = 2μ_1$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 Stress–strain adiabats for different values of S. Left $μ_1 = μ_2$. Right $μ_2 = 2μ_1$</figcaption></figure>

Then the expression for the entropy takes the form

$$S = kN \left( \frac{1}{2} + \log(2\pi) - \log(\beta) - x^* \log(x^*) - (1 - x^*) \log(1 - x^*) - \frac{x^*}{2} \log(\mu_1) - \frac{(1 - x^*)}{2} \log(\mu_2) \right). \quad (32)$$

This formula is illustrated in Fig. 9 where we plot the stress–strain curves $\bar{\sigma}(\bar{\epsilon}, S)$ for different values of entropy. To this end, the equation $S = \text{const}$ is first solved and the value of $x^*$ is expressed via $\beta$. This value is then substituted into the stress–strain relation.

To compare our thermodynamical theory with the results of the direct numerical simulations we would need to know the temperature–strain distribution along an adiabat $\beta(\bar{\epsilon}, S)$. This relation is shown in Fig. 10 for different values of entropy. Notice that in the mixed-phase region the temperature decreases (see also [26]), which can be explained by the high entropy of the corresponding two-phase inhomogeneous configuration. Indeed, one lattice “melts” before it becomes another lattice and it is then expected that the entropy of the “liquid” phase is above the entropies of the two coexisting solid phases.

A simple formula for the maximum temperature drop associated with this “solid to liquid” transition can be obtained in the case $μ_1 = μ_2$. We observe that in this case $x^* = 0.5$ at $\bar{\epsilon} = 0$ independently of $\beta$ and write

$$S = \frac{1}{2} + \log(2\pi) - \log(\beta^t) - x^* \log(x^*) - (1 - x^*) \log(1 - x^*) = \frac{1}{2} + \log(2\pi) - \log(\beta^t) + \log(2),$$
<figure><img src="image_12.png" /><figcaption>Fig. 10 Temperature–strain relation along the adiabats. Left $\mu_1 = \mu_2$. Right $\mu_2 = 2\mu_1$. On the right plot one can easily check that for the limiting values $\beta_1^2\mu_1 = \beta_2^2\mu_2$</figcaption></figure>

where $\beta^t$ corresponds to the value of temperature at $\bar{\epsilon} = 0$. Once the chain reaches the state where all the springs are in one phase, then $x^* = 0$ or $x^* = 1$, and

$$S = \frac{1}{2} + \log(2\pi) - \log(\beta^t) - x^* \log(x^*) - (1 - x^*) \log(1 - x^*) = \frac{1}{2} + \log(2\pi) - \log(\beta^{\pm}),$$

where $\beta^{\pm}$ characterize the temperature in the single phase states. Since the entropy of the chain does not change along the adiabat, we obtain $\beta^t/\beta^{\pm} = 2$. In the general case of asymmetric wells, $\mu_1 \neq \mu_2$, we can obtain $\beta^t/\beta^{\pm} = \frac{1}{1-x^*}(1/x^* - 1)^{x^*}$, where $x^*$ must be computed at $\bar{\epsilon} = 0$.

## 3.6 Thermodynamic theory versus direct numerical simulations

We can now compare the thermodynamic theory with the results of our direct numerical simulations. We recall that the theory is based on the plausible but heuristic assumption that the motion induced by the irreversible self-heating is ergodic. We also recall that the obtained explicit formulas are approximate in the sense that we implicitly deal with a double-valued energy obtained by the continuation of each of the two parabolic wells beyond the spinodal point into the whole real axis. The last assumption can be easily dropped at the expense of analytical transparency, but we found that our simplified theory is sufficient to show the validity of the thermodynamic approach.

The first problem is to specify the specific adiabat $\beta(\bar{\epsilon}, S)$ resulting from the irreversible self-heating of a ‘cold’ chain. As we have already noted earlier, this transition is irreversible and is accompanied by the entropy increase from $S = 0$ to, say, $S = S_0$. To find the value of $S_0$, we need to know the value of temperature at one value of $\bar{\epsilon}$. Observe that at the spinodal instability threshold the energy which transforms from macro to micro level is known exactly (see Fig. 1). More specifically, half of this (available) energy is transformed into the kinetic energy and half into the elastic energy. The kinetic energy per mode is equal to $kT$ from where the temperature and therefore the value of $\beta$ can be easily computed.

Our Fig. 12 and Fig. 11 compare the adiabats $\bar{\sigma}(\bar{\epsilon})$ and $\beta(\bar{\epsilon})$ obtained using our analytical model with the corresponding results of the direct numerical simulations for the quasi-statically driven chain with 30 springs. In simulations, the numerical results were smoothed using window averaging. We have tested various window sizes. The window size away from the mixed-phase region is taken to be five times larger than that in the mixed-phase region. To compute the stress–strain relation from the statistical model, we obtained the value of $S_0$, generated associated adiabat $\beta(\bar{\epsilon}, S_0)$ and substituted it into the stress strain relation $\bar{\sigma}(\beta, \bar{\epsilon})$.

The stress–strain curves presented in Fig. 11 show a very good agreement between the statistical model and the direct numerical simulations. The temperature–strain distributions $\beta(\bar{\epsilon}, S_0)$ shown in Fig. 12 exhibit a slight mismatch in mixed-phase region, which we associate with the finite averaging window.

The overall agreement between the theory and the numerical experiment can be viewed as a strong argument in favor of the assumption of ergodicity achieved during spinodal instability. It remains to be understood
<figure><img src="image_13.png" /><figcaption>Fig. 11 Comparison of the stress–strain curves obtained in direct numerical simulations (dotted line) with the adiabats computed in the statistical model (solid line). $\bar{\sigma}$ versus $\bar{\epsilon}$. Left $\mu_1 = \mu_2$. Right $\mu_2 = 2\mu_1$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12 Comparison of the direct numerical simulations (dotted line) with the statistical model (solid line). $\beta$ versus $\bar{\epsilon}$. Left $\mu_1 = \mu_2$. Right $\mu_2 = 2\mu_1$</figcaption></figure>

at the formal level why the available potential energy at the instability point is exactly sufficient to thermalize the whole system.

## 4 Miscellaneous

In this section we discuss two issues tangentially related to the main subject of the paper. One concerns the dependence of the results on the type of the loading device and the other one concerns the role of small damping.

We begin with the comparison of the thermodynamical models for a chain with $N$ elements loaded in soft and hard devices. Numerical experiments and statistical modeling for the case of soft device was done in [26] for symmetric energy, and here we generalize these results for the case of asymmetric energy wells.

To analyze the chain in a soft device we need to modify the energy by adding the work of the loading device. The controlling parameter now is the applied stress $\bar{\sigma}$ and the configurational statistical sum takes the form

$$Z_p(\bar{\sigma}, \beta) = \int_{\epsilon_i \in R} \exp\left(-\beta\left(\sum_{i=1}^{N} \frac{\mu_i}{2}(\epsilon_i + \alpha_i)^2 - \bar{\sigma}\sum_{i=1}^{N} \epsilon_i\right)\right) d\epsilon_1, \dots, d\epsilon_N. \quad (33)$$
<figure><img src="image_15.png" /><figcaption>Fig. 13 Stress-strain isotherms for a chain with 30 springs driven in hard and soft devices. Left figure high temperature. Right figure low temperature</figcaption></figure>

One can see that in the case of soft device we deal with noninteracting elements (vs elements interacting through the mean field in the case of hard device). The statistical sum is then a product of the contributions due to independent elements, and we do not even need our simplifying assumption to obtain an analytic result

$$Z_p(\bar{\sigma}, \beta) = \left( \frac{\sqrt{\pi/2}}{\sqrt{\beta\mu_1}} \exp(\bar{\sigma}\beta) \exp\left(\frac{\beta\bar{\sigma}^2}{2\mu_1}\right) \text{ERFC}\left(-\sqrt{\beta\mu_1}(1+\sigma/\mu_1)\right) \right. \\ \left. + \frac{\sqrt{\pi/2}}{\sqrt{\beta\mu_2}} \exp(-\sigma\beta) \exp\left(\frac{\beta\sigma^2}{2\mu_2}\right) (1 - \text{ERFC}\left(-\sqrt{\beta\mu_2}(1+\sigma/\mu_2)\right)) \right)^N, \quad (34)$$

where ERFC is the complementary error function. The stress-strain isotherms for the case of soft device are compared with their hard device analogs (computed using our analytical model) in Fig. 13. One can see that at small N and large β ≫ 1/μ the predictions of the two thermodynamical models can be rather different which is well known for the purely mechanical problem β → ∞ where in the soft device the system ‘snaps’ and in the hard device it ‘pops’ [18]. As one can expect, in the thermodynamical limit the mode of loading becomes irrelevant.

The next issue which is appropriate to discuss in this Section concerns the role of small dissipation. We recall that the discrete system studied in the previous sections was Hamiltonian and therefore after initial ‘self-heating’ which erased the memory of the special preparation of the system, the effective macroscopic behavior under quasi-static cyclic loading was reversible. In many realistic systems with bi- or multi-stability at the micro-scale such reversibility, which can be associated with rubber-type behavior, is replaced by rate independent hysteresis.

To establish a link between the behavior of the undamped system and the overdamped behavior of a hysteretic system [19] we consider a bi-stable FPU system with finite visco-elasticity at the micro-scale(different type of dissipation was introduced in [1,26]). More precisely, we consider springs with viscosity D and write the system of governing equations in the form

$$\rho h \frac{d^2 x_i}{dt^2} - D \frac{d}{dt} \left( \frac{x_{i+1} - 2x_i + x_{i-1}}{h} \right) = E \sigma \left( \frac{x_{i+1} - x_i}{h} \right) - E \sigma \left( \frac{x_i - x_{i-1}}{h} \right) \quad (35) \\ x_0 = 0, \quad x_N(t) = x_N^0 + Vt.$$

Next we introduce nondimensional variables

$$M = V^2 \rho / E, \quad W = DV / (EL), \quad \delta = h / L = 1 / N$$
<figure><img src="image_16.png" /><figcaption>Fig. 14 Kinetic energy versus time for a damped chain during one cycle of loading and unloading</figcaption></figure>

and rewrite the system (35) in the dimensionless form

$$ M\delta \frac{d^2 x_i}{dt^2} - W \frac{d}{dt} \left( \frac{x_{i+1} - 2x_i + x_{i-1}}{\delta} \right) = \sigma \left( \frac{x_{i+1} - x_i}{\delta} \right) - \sigma \left( \frac{x_i - x_{i-1}}{\delta} \right) \quad (36) \\ x_0 = 0, \quad x_N(t) = x_N^0 + t. $$

In the previous sections we studied the limit $W = 0, M \to 0$, while the asymptotics considered in [19] corresponded to the limit $M = 0, W \to 0$. It is then natural to perform numerical simulations for the case when both $M$ and $W$ are small but finite. We assume that $\delta$ is finite, leaving for future studies the comprehensive analysis of the triple asymptotic structure for the singular point $M = 0, W = 0, \delta = 0$.

In Fig. 14 we show the time dependence of the kinetic energy for the initially ‘cold’ chain with 14 springs loaded quasi-statically in a hard device. We show one loading path (solid line) and one unloading path (dotted line). One can see that the elements start transforming from one phase to another at the spinodal point. The transition takes place in the form of avalanches accompanied by the emergence of oscillations. The oscillations quickly disappear due to dissipation before the arrival of the next avalanche. During unloading similar behavior of kinetic energy can be observed again with three major transformation events. In Fig. 15, we present the corresponding stress–strain diagram. Here also the solid line corresponds to elongation, and the dotted line corresponds to contraction. The overall behavior of the system is clearly hysteretic which is due to finite dissipation.

As the value of the parameter $W$ increases, the individual avalanches get smaller and occur in shorter time intervals. The picture asymptotically approaches the one predicted for the overdamped system in [19].

## 5 Conclusions

Our highly idealized model suggests that even in quasi-static conditions the homogenized description of the mechanical structures with internal instability (internal buckling) must be necessarily *thermodynamical* and that the canonical distribution may already be adequate for systems with relatively small number of degrees of freedom.

As we showed, in order to reach the thermodynamic level of description for our bi-stable FPU chain one needs to consider a double limit: $\delta \to 0$ (discrete to continuum) and $v \to 0$ (quasi-static driving). The naive long-wave homogenization of the FPU chain suggests that the limiting problem will be of the form

$$ \sigma(\epsilon)_x = 0, \quad (37) $$
<figure><img src="image_17.png" /><figcaption>Fig. 15 Stress–strain relation for a damped chain during one cycle of loading and unloading</figcaption></figure>

where $\sigma(\epsilon) = e'(\epsilon)$ is the reversible stress strain relation for the elastic material. A more careful line of reasoning going back to von Neumann [23] suggests that the limiting equations will be of thermo-elasticity

$$ \sigma(\epsilon, T)_x = 0, \quad S(\epsilon, T)_t = 0 \qquad (38) $$

where $\sigma(\epsilon, T) = f_c(\epsilon, T)$ is the stress, $S(\epsilon, T) = -f_T(\epsilon, T)$ is the entropy, $T$ is the temperature and $f(\epsilon, T)$ is the free energy. In this paper we derive an explicit relation for the $f(\epsilon, T)$ for the bi-stable FPU chain. Our results suggest that the choice between (37) and (38) depends on the convexity of the function $e(\epsilon)$: if it is convex, we expect (37) to be true for the case of a quasi-static driving of initially ‘cold’ chain, while if it is non-convex, one eventually would have to use (38).

Despite its simplicity our model shows both reversible heat exchange associated with release and absorption of the latent heat and irreversible heating due to energy dissipation. It is also a first example when a system can self-thermalize due to microscopic instability. One interesting extension of the model would be a chain with bi-stable nearest and linear ferromagnetic next to nearest neighbor interactions (e.g. [22]). Such model is expected to exhibit more realistic transformation mechanism with localized phase boundaries separating the two phases. Introducing nonlinearity inside single phases will contribute to thermalization and show a rudiments of diffusive heat transfer. On the mathematical side, the formal justification of the use of canonical measure for the simplest bi-stable FPU model remains a challenge.

**Acknowledgments** We would like to acknowledge useful discussions with Pierre Collet, Michael Ortiz and Henning Struchtrup.

## References

1. Balk, A.M., Cherkaev, A.V., Slepyan, L.I.: Dynamics of chains with non-monotone stress–strain relations. I. Model and numerical experiments. J. Mech. Phys. Solids 49, 131–148 (2001)

2. Berdichevsky, V.: Thermodynamics of Chaos and Order. Addison Wesley, New York (1997)

3. Blanc, F.L.X., Le Bris, C., Patz, C.: Finite-temperature coarse-graining of one-dimensional models: mathematical analysis and computational approaches. J. Nonlinear Sci. 20(2), 241–275 (2010)

4. Campbell, D.K., Rosenau, P., Zaslavsky, G.M.: The Fermi–Pasta–Ulam problem—the first fifty years. Chaos 15, 015101 (2005)

5. Cardin, F., Favretti, M.: On the Helmholtz-Boltzmann thermodynamics of mechanical systems. Contin. Mech. Thermo-dyn. 16, 15–29 (2004)

6. da Fonseca, C., Petronilho, J.: Explicit inverse of a tridiagonal k-toeplitz matrix. Numer. Math. 100, 457–482 (2005)

7. Dreyer, W., Kunik, M.: Cold, thermal and oscillator closure of the atomic chain. J. Phys. A 33, 2097–2129 (2000)

8. Fedelich, B., Zanzotto, G.: Hysteresis in discrete systems of possibly interacting elements with a double-well energy. J. Nonlinear Sci. 2, 319–342 (1992)

9. Fermi, E., Pasta, J., Ulam, S.: Studies of Nonlinear Problems. Los Alamos Technical Reports, Document LA-1940 (1955)
10. Ford, J.: The Fermi–Pasta–Ulam problem: paradox turns discovery. Physics Reports, pp. 271–310 (1992)

11. Gallavotti, E.G.: The Fermi–Pasta–Ulam Problem. Springer, New York (2007)

12. Giada Basile, S.O., Spohn, H.: Energy transport in stochastically perturbed lattice dynamics. Arch. Ration. Mech. Anal. **195**(1), 171–203 (2010)

13. Jackson, E.: Perspectives of Nonlinear Dynamics. Cambridge University Press, Cambridge (1991)

14. Kastner, O.: Molecular-dynamics of a 2D model of the shape memory effect, part II: thermodynamics of a small system. J. Mech. Phys. Solids **18**, 63–81 (2006)

15. Müller, I., Seelecke, S.: Thermodynamic aspects of shape memory alloys. Math. Comput. Model. **34**, 1307–1355 (2001) (Topics in the mathematical modelling of smart materials)

16. Muller, I., Villaggio, P.: A model for an elastoplastic body. Arch. Ration. Mech. Anal. **65**, 25–46 (1977)

17. Ngan, S., Truskinovsky, L.: Thermoelastic aspects of nucleation in solids. J. Mech. Phys. Solids. **50**, 1193–1229 (2002)

18. Puglisi, G., Truskinovsky, L.: Mechanics of a discrete chain with bi-stable elements. J. Mech. Phys. Solids. **48**, 1–27 (2000)

19. Puglisi, G., Truskinovsky, L.: Thermodynamics of rate independent plasticity. J. Mech. Phys. Solids **53**, 655–679 (2005)

20. Stuart, A., Humphries, A.: Dynamical Systems and Numerical Analysis. Cambridge University Press, Cambridge (1996)

21. Truskinovsky, L.: Nucleation and growth in elastodynamics. In: Duxbury, P.M., Pence, T. (eds.) Dynamics of Crystal Surfaces and Interfaces, pp. 185–197. Plenum, New York (1997)

22. Truskinovsky, L., Vainchtein, A.: About the origin of the nucleation peak in transformational plasticity. J. Mech. Phys. Solids **52**, 1421–1446 (2004)

23. von Neumann, J.: Proposal and analysis of a numerical method for the treatment of hydrodynamical shock problem, 1944. In: John von Neumann Collected Works, vol. 6, no. 26. Macmillan, New York (1963)

24. Weinan, E., Ming, P.: Cauchy-born rule and the stability of crystalline solids: dynamic problems. Acta Math. Appl. Sin. Engl. Ser. **23**(4), 529–550 (2007)

25. Weissert, T.: The Genesis of Simulations in Dynamics: Pursuing the Fermi–Pasta–Ulam Problem. Springer, New York (1997)

26. Yoshikawa, T.: Transition from mechanical equilibrium to thermal equilibrium in a chain of unstable springs. Phys. Lett. A **307**, 13–21 (2003)