Thermalization of a driven bi-stable FPU chain

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.


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, . . ., N the coordinates of the particles and by e( ) the elastic energy of the springs.Then where σ ( ) = e ( ).We assume zero initial velocities and apply the following time-dependent boundary conditions: The system (1) depends on the two dimensionless parameters.The first parameter, measuring the degree of discreteness of the system, is δ = 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 /( √ E/ρ), where ρ 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 δ → 0 (thermodynamic limit) and v → 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.

Numerical studies
Assume that δ = 1 and denote by j = x j+1 − x j the strain in the jth spring, j = 0, . . ., N − 1.To make analytical computations possible, we assume that the function e( ) is bi-quadratic where the spinodal strain is a 0 = (1 − κ)/(1 + κ), κ = √ μ 2 /μ 1 (see Fig. 1).While one of the parameters μ 1 or μ 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 < a 0 belongs to Phase I and the interval > a 0 -to Phase II.The spinodal region is represented in this model by a single-point = 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).
Fig. 1 Stress-strain and energy-strain relations for a bi-linear snap-spring.Available energy is released during spinodal instability 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 The general motion can then be represented in the form where ω 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 ω k (t) is the solution of the linear system 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 μ = 1 (otherwise we can re-scale the time) we obtain that the eigenvalues of A are λ l = −2 + 2 cos lπ N +1 and the corresponding eigenvectors are Once the eigenvectors and eigenvalues are known, the matrix A can be represented as where D is a diagonal matrix consisting of eigenvalues and the columns of U are corresponding eigenvectors.
The the system can be represented as a set of noninteracting harmonic oscillators d dt 2 φ = Dφ, where φ = U ω.The energy of each mode (each oscillator) is then Consider now the symmetric bi-quadratic energy with μ 1 = μ 2 = μ.In this case the definition of the modes can be formally extended into the regime where the springs change phase and where oscillators are  2 Total energy of the modes for the system of five springs during quasistatic stretching of an initially 'cold' chain 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 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 where k is a dimensional (Boltzmann) constant and τ is the averaging window.One can see that after a finite time the system achieves the state of equipartition with 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.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 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 μ 1 = μ 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 where E μ i is the total energy of the modes in the single phase stage of deformation with stiffness μ i , (i = 1, 2).In terms of the temperature of the particles this relation translates into 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.

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]).To define the ensemble we fix the total strain = N i=1 i and denote the mean strain by ¯ = /N .We introduce a canonical measure by assuming that a microstate with the energy E = N i=1 , where β = (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 and the integration is carried over all micro-configurations characterized by N −1 strains 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 where E p = N i=1 e( i ) is the total potential energy of the system.After the statistical sum is known, the free energy can be computed as follows: where The expression for the entropy takes the form

Single phase system
To illustrate the above formulas we first consider quadratic potential with e( ) = 1 2 μ 2 .In this case the problem reduces to the computation of a Gaussian integral As is common in such calculations we introduce the new variables r 0 = 0, r k = k i=1 i , and then use the known formulas to obtain The average potential energy takes the form E p = (N −1) 2β + 1 2 N μ¯ 2 .At large N the energy density E p /N splits into the macroscopic elastic energy 1  2 μ 2 and the microscopic energy 1/(2β) = 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 = (∂ E/∂ T ) = k N = 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 One can see that the stress-strain relation at constant temperature remains purely mechanical.The expression for the finite size entropy takes the form 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 which agrees with the relation between the limiting temperatures found in our numerical experiments.

Two-phase system
For the two-phase system with a bi-quadratic energy we can write Here the discrete parameter α which locates the bottoms of the energy wells takes only two values −a 1 and −a 2 .Similarly, μ i takes values μ 1 or μ 2 depending on the value of i .Note that both μ i and α i depend on 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 i = 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: To compute the integrals, we can again introduce the variables r k In terms of these new variables, the function Z p can be written in the form To proceed further, we denote Then ( 13) can be rewritten as where (•, •) denotes scalar product and r = (r 1 , . . ., r N −1 ).Next, we use again the well known property of Gaussian integrals exp to obtain where P −1 N −1,N −1 is the last entry (last column and last row) of the inverse of P(μ).Now, noting that r N = + i α i , the expression for Z p can be rewritten as Every term that involves P depends on μ i and consequently on α i .Thus det(P) can be evaluated explicitly det(P) = N i=1 D/μ i , where D = μ 1 . ..μN . Te term P −1 N −1,N −1 can be written in the closed form (see [6]) where , we obtain the following recursive relation We then note that and by induction argument show that Using these expressions we finally obtain 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 μ i and can be rewritten in terms of the number of springs in the first phase p.First we observe that Then, we can write the following explicit expression capturing the size effect:

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 β → ∞ 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 β one can approximate the partition function Z p by where the matrix A is defined by and x 0 is the point of the minimum of P(x).Now, suppose for simplicity that we are dealing with equal moduli case and assume that E( 1 , . . ., N ) has minimum E min .In the neighborhood of the minimum E has quadratic behavior with therefore, by using the Laplace formula ( 21) we obtain By using the explicit expression for the ground state configurations [8,18], we can rewrite (22) in the form In terms of the average strain = /N we obtain The corresponding stress-strain relation reads where p is a piece wise constant function of (see [8,18] for more detail).

Thermodynamic limit
To study the limiting case N → ∞, we introduce the new variable x = p/N and write If we now denote and replace the sum by an integral, we obtain Since N is large we can use Laplace method which gives where x * is a minimum of ψ(x).
For large | |, the minimum of ψ(x) is either at x = 0 or at x = 1.More precisely, when ψ(x) is increasing we have while for decreasing ψ(x), we obtain Here ψ (0) and ψ (1) need to be determined by computing the corresponding derivatives with respect to (at constant β).One can see that we recover the thermodynamics of the single phases.
For the intermediate values of , we can write Using (29) we can rewrite equation for x * in the following way: From here we obtain .
This allows us to write the expression for the function (σ ) in the closed form In the case of symmetric double well energy with equal moduli μ 1 = μ 2 = 1 and α 1 = −1, α 2 = 1 we obtain and 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/ √ μ 1 and a 2 = −1/ √ μ 2 .In all cases, the discontinuity point for the stress-strain relationship is at = 0.In Fig. 8, we plot the isothermal stress-strain curves at different values of β.At large β (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.

Adiabats in the thermodynamic limit
The general expression for the entropy is given in (7).For large N , we can write and Noting that ∂ψ(x * )/∂ x * = 0, we obtain Then the expression for the entropy takes the form This formula is illustrated in Fig. 9 where we plot the stress-strain curves σ ( , S) for different values of entropy.To this end, the equation S = const is first solved and the value of x * is expressed via β.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 β( , 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 = 0 independently of β and write On the right plot one can easily check that for the limiting values where β t corresponds to the value of temperature at = 0. Once the chain reaches the state where all the springs are in one phase, then x * = 0 or x * = 1, and where β ± characterize the temperature in the single phase states.Since the entropy of the chain does not change along the adiabat, we obtain β t /β ± = 2.In the general case of asymmetric wells, μ 1 = μ 2 , we can obtain x * , where x * must be computed at = 0.

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 β( , 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 .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 β can be easily computed.
Our Fig. 12 and Fig. 11 compare the adiabats σ ( ) and β( ) 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 β( , S 0 ) and substituted it into the stress strain relation σ (β, ).
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 β( , 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 at the formal level why the available potential energy at the instability point is exactly sufficient to thermalize the whole system.

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 σ and the configurational statistical sum takes the form 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 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 Next we introduce nondimensional variables In the previous sections we studied the limit W = 0, M → 0, while the asymptotics considered in [19] corresponded to the limit M = 0, W → 0. It is then natural to perform numerical simulations for the case when both M and W are small but finite.We assume that δ is finite, leaving for future studies the comprehensive analysis of the triple asymptotic structure for the singular point M = 0, W = 0, δ = 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].

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: δ → 0 (discrete to continuum) and v → 0 (quasi-static driving).The naive long-wave homogenization of the FPU chain suggests that the limiting problem will be of the form σ (ε) x = 0, (37) where σ (ε) = e (ε) 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 σ (ε, T ) x = 0, S(ε, T ) t = 0 (38) where σ (ε, T ) = f ε (ε, T ) is the stress, S(ε, T ) = − f T (ε, T ) is the entropy, T is the temperature and f (ε, T ) is the free energy.In this paper we derive an explicit relation for the f (ε, 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(ε): 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 absorbtion 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.
Fig.2Total energy of the modes for the system of five springs during quasistatic stretching of an initially 'cold' chain

Fig. 3 Fig. 4
Fig.3The motion of the chain with 80 particles after spinodal instability

Fig. 5
Fig. 5 Total energy of linear modes versus time for the system with μ 1 = μ 2 = 1

)Fig. 13
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

Fig. 14
Fig. 14 Kinetic energy versus time for a damped chain during one cycle of loading and unloading

Fig. 15
Fig. 15 Stress-strain relation for a damped chain during one cycle of loading and unloading