# Inertia-induced power law scaling in martensites

Oguz Umut Salman, Alphonse Finel, Lev Truskinovsky

► To cite this version:

Oguz Umut Salman, Alphonse Finel, Lev Truskinovsky. Inertia-induced power law scaling in martensites. Mathematics and Mechanics of Solids, 2025, 30 (9), pp.1896 - 1934. 10.1177/10812865251361074. hal-05390223

HAL Id: hal-05390223

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

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

# Inertia-induced power law scaling in martensites

Mathematics and Mechanics of Solids
I–39

© The Author(s) 2025
Article reuse guidelines:
sagepub.com/journals-permissions
DOI: 10.1177/10812865251361074
journals.sagepub.com/home/mms

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

Sage

Oğuz Umut Salman

Laboratoire des Sciences des Procédés et des Matériaux (LSPM), CNRS, Université Sorbonne Paris Nord,
UPR 3407, Villetaneuse, France; Laboratoire d'Etudes des Microstructures, ONERA, Chatillon, France

Alphonse Finel

Université Paris-Saclay, ONERA, CNRS, Laboratoire d'Etude des Microstructures (LEM), 92322, Châtillon,
France

Lev Truskinovsky

Physique et Mécanique des Milieux Hétérogènes (PMMH), ESPCI Paris, PSL University, Sorbonne
Université, Université de Paris, CNRS, Paris, France

Received 11 March 2025; accepted 16 April 2025

Abstract

While martensites subjected to quasi-static deformation are known to exhibit power law distributed acoustic emission in a broad range of scales the origin of the observed scaling behavior and the mechanism for self-organization towards criticality remains obscure. Here, we argue that the power-law structure of intermittent fluctuations can be at least partially attributed to inertia. We build on the insight that inertial dynamics, evidenced by acoustic emission, can become an important factor if the underlying mechanical system is only marginally stable. We first illustrate the possibility of inertia-induced *heavy-tailed* avalanche size distributions using a prototypical example of a discrete chain with bi-stable springs. We then explore the effects of inertia in fully realistic two- and three-dimensional continuum models of elastic phase transitions. In particular, we demonstrate that a three-dimensional model of this type can produce not only qualitative but also quantitative agreement with experiment.

Keywords

Criticality, power laws, martensitic phase transitions, inertia

I. Introduction

Elastic crystals undergoing diffusionless (martensitic) phase transformations are known to belong to
the class of systems exhibiting avalanche-type dynamics under quasi-static driving [1–5]. They can be
then viewed as inherently nonlinear complex systems with intermittent response resulting from cascades

<footnote>**Corresponding author:**
Lev Truskinovsky, Physique et Mécanique des Milieux Hétérogènes (PMMH), ESPCI Paris, PSL University, Sorbonne Université, Université de Paris, CNRS, F-75005 Paris, France.</footnote>

<footnote>Email: lev.truskinovsky@espci.fr</footnote>
of internal instabilities [6–14]. Typical examples of such systems, usually characterized by threshold-
type nonlinearity, long-range interactions, and self-induced disorder, include terrestrial earthquakes,
Barkhausen noise, and amorphous plasticity, to mention just a few [15–21]. A salient feature of all these
systems is that the observable quantities, representing the underlying collective behavior, exhibit heavy-
tailed distributions. If no fine tuning is needed to reach the scaling regime, the emerging *quasi-criticality* is
interpreted either in the framework of self-organized criticality (SOC) [22–29] or within a closely related
paradigm of marginal stability [30–33].

In this paper, we address the origin of intermittent response in quasi-statically driven martensites.
We recall that martensitic transformations are structural, athermal, shear-dominated displacive phase
transitions [34,35]. The presence of lattice mismatch between the low-symmetry martensite and the high-
symmetry austenite phases leads to long-range elastic interactions which conspire with elastic energy
non-convexity in producing complex multiscale microstructures [36–38].

It has long been known that mechanically (or thermally) driven martensitic phase transitions are
accompanied by broadly distributed bursts (avalanches) which are audible. More recently, it has become
clear that the intermittent transformation advance generates detectable acoustic emission (AE) which
carries important information about the mechanism of martensitic transformations. Thus, it reveals that
between consecutive avalanches, the transformation is suppressed while the system deforms purely elas-
tically. Such elastic branches terminate when the system reaches consecutive instability thresholds. It
was established empirically that when the driving is quasi-static, which prevents the overlap of indi-
vidual avalanches, the amplitudes and durations of AE hits (induced by individual elastic instability
events) exhibit power law size distribution. The values of the corresponding exponents were shown to
depend only on crystal symmetry, which points the universality of the underlying non-equilibrium steady
states [1–3,5,14,39–49].

The exact origin of power law distributed avalanches in martensites is still a subject of debate. For
instance, it has been argued that the main factor behind the observed scaling behavior is quenched dis-
order [50]. To corroborate this idea, a driven zero-temperature Random Field Ising Model (RFIM) with
short range interactions of ferromagnetic type was studied extensively [51–53]; same idea is behind paral-
lel studies of the Random Bond Ising Model [54], the Diluted Ising Model [55] and the Random Anisotropy
Ising Model [56]. A shortcoming of these models is that a power law emerges there only at a critical level
of quenched disorder [50–53,57], which points towards 'tuned' rather than 'generic' scaling behavior. It
has been argued, however, that the near-critical domain in the corresponding parameter space may be so
large that it emerges in experiment as an extended criticality [50,51,58].

A different, but closely related, interpretation of the generic nature of scaling in martensites was
proposed in [59], where it was linked to the presence of a limited dislocational activity taking place con-
currently with the martensitic transformation. In this interpretation, an external tuning of disorder is
replaced by self-tuning of the annealed inhomogeneity which takes the form of a co-evolving dislocation
distribution. Based on the analysis of the corresponding stylized model, it was shown that the implied
self-induced disorder is correlated which facilitates the self-organization of the system towards criticality
[4,59–61]. This idea was corroborated to some extent by the observations that the power law behavior
in martensites emerges only after cyclic loading (known as 'training') which ensures that the necessary
level of self-organization has been reached [39,62]. Despite its plausibility, the universality of such an
explanation has been also challenged with reference to some experiments where power law distributed
avalanches were observed without training, already during the first loading cycle and within samples
exhibiting no clear signs of plastic slip [63].

Yet another proposed interpretation of scale invariance in martensites, still relying on the umbrella
RFIM framework, is linking it to the antiferromagnetic nature of elastic long-range interactions. It
suggests that the crucial ingredient missing in the standard RFIM model, is the anisotropy of elastic
interactions which renders corresponding interaction kernels sign indefinite [64,65]. It was shown that
the presence of such kernels changes the nature of interactions fundamentally and can drastically modify
the spatial structure of avalanches. In particular, it was shown that in a model setting of this type, a robust
self-organization towards scaling regime becomes possible [57]. This idea was recently reinforced by an
extensive numerical study of a comprehensive long-range version of the RFIM [66]. However, given that
in the proposed models the role of disorder remains opaque while the robust nature of the scaling is not
sufficiently manifested, the conclusion about the 'generic' nature of the associated criticality remains
still debatable.

In this paper, we complement those studies by focusing on the possibility that *inertia* may also be one of the factors responsible for generic (extended) scaling behavior during martensitic transformations. We elaborate on the evidence that none of the above models takes into account the fact that martensitic phase boundaries move at speeds close to sound velocity, not only producing lattice scale waves but also generating elastic shocks which may affect the transformation [67–76]. The implied scattering and focusing of elastic waves suggests rich dynamics, which is revealed by the detected intense AE [77]. The underlying dynamic activity is not dissimilar to the one generated during terrestrial earthquakes [78,79]. In this respect, of particular importance to our study is the ability of seismic waves to nucleate secondary slips, which is the phenomenon known as *dynamic triggering* of earthquakes [80–82]. This analogy is a reminder that during martensitic transformations, the mechanical system should be considered as fundamentally *underdamped*.

An important role of inertia in the process of self-organization toward criticality and scaling has long been realized in the studies of avalanches in sandpiles. For instance, effects of inertia were engaged to explain the failure to achieve full criticality in physical experiments with actual sand and to justify skewed nature of the associated avalanche shapes [83–89]. However, in these studies, the underdamped nature of the system was modeled only *indirectly*, by accounting for inertia-induced threshold weakening. Inertial effects were also implied in the closely related models of critical behavior caused by kinetic softening, where the underdamped nature of the model was represented by the non-monotonicity of the flow curve [90–94]. A well-known example of such studies is the Burridge–Knopoff model of earthquakes, where inertia is hidden under the phenomenological assumption that dynamic friction is lower than static friction [95]. Behind all these indirect representations of inertia is the idea that dynamically generated acoustic waves allow the system to jump over barriers. The implied positive feedback was shown to be responsible for stick-slip mechanical behavior which is ultimately behind the observed intermittency [92,93,96–98]. An interesting feature of the velocity-weakening-type models of inertial effects is the prediction of the prevalence of large avalanches known as supercriticality.

A crucial importance of the direct account for inertia in the modeling of martensitic transformations has been long realized [67, 70, 99, 100]. Thus, already in an early studies of a one-dimensional continuum model of a prototypical martensitic transformation, it was shown that the formation of twin microstructure can be interpreted as a purely inertial effect which disappears in an overdamped setting, see, for instance, [101]. Inertial effects were then shown to be behind the development of intermittency in a closely related fully dynamic Frenkel–Kontorova model where power law scaling was shown to emerge without any quenched or annealed disorder [102]. The possibility of intermittency and scaling in a molecular dynamics model of a martensitic transformation was demonstrated by Ding et al. [103]; however, the implicitly present inertial effects were not specifically identified as a factor contributing to self-organization towards criticality.

The first continuum model dedicated directly to the study of inertia-induced criticality in martensites was proposed in [104] and [105], see also its subsequent developments in [106]. The authors considered a square-to-rectangle transition in the setting of Ginzburg–Landau model with kinetic energy taken into account. Numerical experiments with such a model, imitating cooling and heating runs, showed the intermittent nature of time series describing the volume fraction of the emerging phase. It was also shown that energy is dissipating in the form of bursts (avalanches) whose statistical distribution follows a power law. To overcome some technical difficulties, the authors had to use a phenomenological truncated elastic kernel, so the elastic interactions in this model were effectively of a short-range type. However, even in such a reduced model, the presence of inertia was shown to ensure that the robust scaling regime can be reached without any involvement of quenched disorder. Notwithstanding its pioneering nature, this work left unanswered the questions about the role of the degree of underdamping, the effect of antiferromagnetic elastic interactions and, the potential importance of elastic incompatibility [107].

In view of its role played in what follows, we recall that the implied *incompatibility* concerns the structure of the stretch tensor *U* that maps the undeformed austenite lattice to one of the equivalent variants of the undeformed martensite lattice. Specifically martensite and austenite phases are considered to be geometrically compatible if there exists a rotation *R* such that *RU* − *I* is a rank-one tensor. Here, *I* is a unit tensor. As it is well known [38], in the case of compatibility, there exist stress-free coherent
interfaces between austenite and martensite phase that are very mobile, allowing the resulting 'thermoelastic' transformation to proceed almost reversibly with only a minimal hysteresis [108–111]. Observe that such an algebraic compatibility of the energy wells should not be confused with the differential elastic compatibility constraint imposed on the strain tensor [106, 112, 113].

The goal of this paper, where we limit our attention to athermal systems and, to avoid spurious time scales, focus exclusively on the case of quasi-static driving, is two-fold. First, we develop an intentionally oversimplified model of an inertial martensitic transformation that can be used to demonstrate the crucial role of the underdamping for reaching the scaling regime. This is achieved by the systematic study of the *crossover* between the limiting overdamped and undamped regimes, neither of which by itself is compatible with power law scaling. Second, to complement the insights obtained from the study of a prototypical model, we perform large-scale numerical experiments with realistic continuum 2D and 3D models simulating some basic martensitic transformations. In particular, this allows us to juxtapose the results for models with compatible and incompatible energy wells. To compare numerical predictions with experimental results, we use in our 3D modeling the parameters that match the available data for single crystals of $Fe_{68.8}Pd_{31.2}$ undergoing a martensitic cubic-to-tetragonal transition. We show that even a basic continuum model of this transformation accounting for both inertia and incompatibility delivers an excellent numerical reproduction of the power law statistics of avalanches with computed exponent values being very close to the ones measured in AE experiments previously conducted on the corresponding crystals [63].

Our important finding on the theoretical side is that a simple 1D discrete snap-spring model is already capable of generating heavy-tailed avalanche size distributions given that the value of a non-dimensional parameter, quantifying the relative strength of inertia versus damping, is chosen in an appropriate range. Interestingly, the ability of such a model to support intermittent fluctuations is lost if one takes a straightforward continuum limit. The problem persists if we move from 1D to 2D continuum setting, while maintaining elastic compatibility of the energy wells. The heavy-tailed distribution is recovered in the continuum setting only if we advance further towards a 3D model while also exacerbating the effect of incompatibility between the energy wells by accounting for a volumetric effect of the transformation. These observations suggests that the abundance of metastable energy wells, emerging in our 1D model due to its *discreteness*, can be viewed as an oversimplified reproduction of the actual complexity of the energy landscape in a 3D continuum model with elastic *incompatibility*. We argue that both of these models effectively describe the same mechanism ensuring 'generic' scale invariance.

The rest of the paper is organized as follows. Our 1D lattice model is introduced in Section II. We analyze separately damped, undamped, and underdamped regimes and show that only in the underdamped regime, the system exhibits intermittency and scaling. We compute numerically the corresponding power law exponents and confirm the validity of the known scaling relations. Section III contains the analysis of a 1D continuum analog of our discrete model obtained by taking a formal continuum (thermodynamic) limit. This model is shown to exhibit neither intermittency nor scaling. In Section IV, we consider a more realistic 2D continuum model and, in this way, limit our attention to the case of compatible energy wells. Here, again it turns out that the desired self-organization to criticality is not achieved. Still, we are able to present in the 2D case a revealing comparison of the transformation mechanisms in overdamped and underdamped regimes showing a remarkable difference in the complexity of the microstructures reachable by these two different types of dynamics. Section V is dedicated to the numerical study of a fully realistic 3D model of the cubic-to-tetragonal transition with incompatible energy wells. We show that in this case, the system finally succeeds in reaching of the scaling regime. We verify that the computed power law exponents both respect the known scaling relations and match the available experimental data. Our conclusions are summarized in the Section VI, where we also mention some open problems.

## 2. One-dimensional discrete model

Consider a 1D chain composed of mass particles linked by elastic springs. Suppose that each particle interacts with nearest (NN) and next to nearest neighbors (NNNs) on each side. Denoted by $u_i$ the horizontal displacement of a particle with index $i=0, ..., N$. We can then write the total energy of the chain in
the form

$$ \mathcal{E} = \mathcal{H} + \mathcal{F}, \qquad (1) $$

where the first term is the kinetic energy

$$ \mathcal{H} = \frac{1}{2} \rho \epsilon \sum_{i=2}^{N-3} \dot{u}_i^2, \qquad (2) $$

and the superimposed dot denotes partial time derivative, $\epsilon$ is the reference inter-particle distance, $\rho$ is the reference mass density. The second term in (1) is the elastic energy

$$ \mathcal{F} = \epsilon \sum_{i=2}^{N-3} (\phi(e_i) + \phi_1(e_i, e_{i-1})). \qquad (3) $$

Here, we introduced the elastic strain

$$ e_i = \frac{u_{i+1} - u_i}{\epsilon}. \qquad (4) $$

The functions $\phi$ and $\phi_1$ introduce NN and NNN elastic interactions, respectively.

To model the simplest elastic (martensitic) phase transition, we assume that the NN potential is bi-stable so that the function $\phi(e)$ has a double-well structure with two energy wells representing different elastic phases. For analytical simplicity, we use in what follows the simplest piecewise quadratic function

$$ \phi(e_i) = \frac{\kappa}{2}(e_i - d)^2, \qquad (5) $$

where $d=0$ for $e_i < e^c$ and $d=a$ for $e_i > e^c$. Here $\kappa$ is the elastic modulus, $a$ is the transformation strain and $e^c$ is the critical strain.

To capture the ferromagnetic nature of the RFIM-type short-range interactions in the simplest form, we assume that the NNN potential $\phi_1$ is of the form

$$ \phi_1(e_i, e_{i-1}) = \frac{\mu}{2\epsilon^2}(e_i - e_{i-1})^2, \qquad (6) $$

where $\mu > 0$ is the second-order elastic modulus.

Finally, to account for environmental friction, we introduce the standard quadratic Rayleigh function

$$ \mathcal{R} = \frac{\gamma}{2} \sum_{i=2}^{N-3} \dot{u}_i^2, \qquad (7) $$

where $\gamma$ is the effective viscosity coefficient.

The dynamic equations describing the resulting finite dimensional discrete mechanical system can be written the form

$$ \frac{\delta \mathcal{L}}{\delta u_i} = \frac{\delta \mathcal{R}}{\delta \dot{u}_i} \qquad (8) $$

where

$$ \mathcal{L} = \mathcal{H} - \mathcal{F} \qquad (9) $$

is the Lagrangian and

$$ \frac{\delta \mathcal{L}}{\delta u} = \frac{\partial \mathcal{L}}{\partial u} - \frac{\partial}{\partial t} \left( \frac{\partial \mathcal{L}}{\partial \dot{u}} \right) \qquad (10) $$

is the variational derivative. We assume that the system is loaded in a hard device and therefore set the boundary conditions in the form

$$ u_0 = u_1 = 0, \dot{u}_{N-1} = \dot{u}_{N-2} = v, \qquad (11) $$
where $v$ is the parameter characterizing the loading rate.

It will be convenient to write the resulting system of equations in a dimensionless form. If we normalize lengths using the system size $L = \epsilon N$ and times using the viscous time scale $\tau = \gamma/\kappa$, we obtain the system

$$C\delta\ddot{u}_i = (\tilde{\phi}'((\tilde{u}_{i+1} - \tilde{u}_i)/\delta) - \tilde{\phi}'((\tilde{u}_i - \tilde{u}_{i-1})/\delta) - \dot{\tilde{u}}_i + (G/\delta^3)(\tilde{u}_{i+2} + \tilde{u}_{i-2} - 4\tilde{u}_{i+1} - 4\tilde{u}_{i-1} + 6\tilde{u}_i)). \quad (12)$$

where $\tilde{u}_i = u_i/L$, $\tilde{t} = t/\tau$ and $\tilde{\phi} = \phi/\kappa$. The dimensionless parameter characterizing the degree of discreteness is

$$\delta = \frac{1}{N}. \quad (13)$$

The relative role of inertia versus damping is measured by the non-dimensional parameter

$$C = \frac{\rho\kappa L^2}{\gamma^2}. \quad (14)$$

The effect of NNN interactions scales with dimensionless parameter

$$G = \frac{\mu}{\kappa}. \quad (15)$$

Yet another dimensionless ratio entering the problem through the boundary conditions characterizes the rate of loading

$$V = \frac{v\gamma}{\kappa L}. \quad (16)$$

In what follows, we consider the behavior of system (12) under the assumption that the two non-dimensional parameters $V$ and $\delta$ are small but finite. The limit $V \to 0$ corresponds to quasi-static driving, and in our numerical experiments, we use the value $V \sim 10^{-5}$. In the limit $\delta \to 0$, we obtain continuum model; however, in all our numerical experiments, we use a small but finite value $\delta \sim 10^{-3}$. The parameter $G$, characterizing the strength of ferromagnetic short-range interactions, is kept sufficiently small to ensure that phase boundaries are localized; specifically throughout the paper, we use the value $G = 10^{-4}$.

In contrast, the parameter $C$, playing the role analogous to a Reynolds number in the theory of turbulence as it characterizes the ratio of inertia to dissipation, is varied in a broad range $C \sim 10 - 10^5$. In this way, we are able to cover different regimes, from fully damped to fully undamped. The goal is to obtain a phase (regime) diagram and search for the range where the system may exhibit power law scaling. To increase statistics, we studied cyclic loading and, in each of the regimes, performed hundreds of cycles.

To solve the equations of motion (12) numerically, we carried out a time-stepping approach using a fourth-order explicit predictor-corrector algorithm. For spatial discretization, we used the fast Fourier transform (FFT) method. The numerical algorithm is discussed in more detail in Appendix 1.

In an attempt to reproduce the results of AE experiments [63], we recorded at each avalanche the associated dimensionless energy dissipation. We started with discretization of the time domain $[0, T]$ introducing the points $t_a = a\Delta t$, where $a = 0, 1, 2, \dots, N$ and $N = T/\Delta t$. As the system is being loaded, an avalanche begins when at least one NN spring starts changing its energy well, and we suppose that this takes place at the time $t = t_i$. During the avalanche, the energy dissipation $\mathcal{R}(t_a)$ is recorded at each time step $t = t_a$ between the time $t = t_i$ and the avalanche termination time $t = t_f$ when the energy wells are not switched anymore. We then compute for such an avalanche the total dissipation

$$E = \sum_{t=t_i}^{t_f} \mathcal{R}(t_a) \Delta t, \quad (17)$$

where the summation is over the corresponding discrete time points $t_a$. In what follows, we associated dissipation only with the fluctuating part of the displacement field

$$w_i(t) = u_i(t) - \frac{i+1}{N} u_{i-1}(t). \quad (18)$$
While it may be natural to adopt as a measure of avalanche size the total transformation strain, one can show [63] that an almost equivalent result is obtained if we use instead a closely related quantity with the same dimensionality

$$S = \sum_{t=t_i}^{t_f} V(t_a) \Delta t, \qquad (19)$$

where

$$V(t) = \mathcal{R}(t)^{1/2}. \qquad (20)$$

In addition to the two macroscopic variables E and S, which are computed for each avalanche, we also record the third macroscopic observable T defined as the avalanche duration

$$T = t_f - t_i. \qquad (21)$$

In this way, we can check the validity of various general scaling relations linking fluctuations in the observables E, S and T.

To complement the avalanche-related information, we also studied a different type of observable which is not tied to intermittent structure of the AE signal. Specifically, we recorded the power spectrum of the time series V(t) defined by the formula [114–119]:

$$PS(f) = \left| \sum_{k=1}^{n} V_k e^{-i2\pi fk} \right|^2, \qquad (22)$$

where $V_k = V(t_k)$ is an abstract dimensionless discrete signal with $t_k = k/n$ and $i = 1, ..., n$.

It is clear that by analyzing function (22), one can identify the presence of correlations in the time series even in the absence of intermittency and avalanches. Indeed, consider an auto-correlation function of the discrete signal $V_k$

$$C(s) = \frac{1}{(N-s)\langle V_k^2 \rangle} \sum_{k=1}^{n-s} (V_k - \langle V_k \rangle)(V_{k+s} - \langle V_k \rangle), \qquad (23)$$

where

$$\langle V_k \rangle = \frac{1}{n} \sum_{k=1}^{n} V_k. \qquad (24)$$

While in the case of a fully uncorrelated white noise $C(s) = 0$, in the presence of short-range correlations, one can expect an exponentially decaying behavior $C(s) \sim \exp(-s/\tau)$ (Lorentzian spectrum). Of particular interest to us will be the time series exhibiting long-range correlations. In this case, one can expect a slower, power law decay behavior for the auto-correlation function

$$C(s) \sim s^{-\gamma}. \qquad (25)$$

We will be using the fact that the corresponding power spectrum is

$$PS(f) \sim 1/f^{\alpha}. \qquad (26)$$

with $\alpha = 1 - \gamma$; in the limiting case $\alpha = 1$ the auto-correlation function is known to exhibit a slower logarithmic decay [120].
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Hysteretic stress-strain response of an overdamped chain subjected to quasi-static cyclic loading in a hard device. Parameters: $N = 1000$ and $C = 80$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Two consecutive individual avalanches in the overdamped regime which is shown in Figure 1. Parameters: $N = 1000$ and $C = 80$.</figcaption></figure>

## 2.1. Overdamped regimes

The regime where dissipation completely overcomes inertia corresponds to the damped limit $C \to 0$. This regime in our 1D discrete setting has been studied analytically before, see for instance [121–124]. It was shown that under cyclic quasi-static loading, the system exhibits rate independent hysteresis with individual springs changing phase (energy well) sequentially, one after another. Dissipation can be then represented as a sequence of periodically spaced identical events (trivial avalanches). Here, we present for completeness the results of numerical experiments in the overdamped regime which corroborate these analytical findings.

In our numerical experiments, we loaded the system from a homogeneous state where springs were in the first energy well until all springs have switched into the second energy well and the system reached again a homogeneous configuration. We then unloaded the chain following the same protocol with both loading and unloading performed quasi-statically. The details of our numerical implementation can be found in Appendix 1.

We found that in the range $C \sim 20 – 80$, the dissipation is still sufficiently strong to deliver the expected over-damped hysteretic response. The computed strain-stress relation is shown in Figure 1. In the horizontal yielding segments of the hysteresis cycle, the system exhibits a succession of almost equal size stress drops. The typical individual avalanches, responsible for such drops, are shown in Figure 2. The avalanches are separated by the quiescent intervals, and each avalanche has roughly the same shape with an exponentially decaying tail. In the studied overdamped regimes with $C \sim 20 – 80$, only a very small number of springs switch from one energy well to another during each avalanche.

An intriguing feature of Figure 1 is the absence of the nucleation peak predicted theoretically for an NNN system in [121]. Such a peak would correspond to a collective nucleation event when many springs change phase simultaneously. It allows subsequent transformation to proceed through a growth of a nucleus and implies spring-by-spring advance of a phase boundary. There are several reasons why the situation here may be different. First of all, we use “extra hard device” boundary conditions (11),
<figure><img src="image_6.png" /><figcaption>Figure 3. Power spectrum of the signal V(t) in the overdamped regime which is shown in Figure 1. Here C = 80 and N = 1000.</figcaption></figure>

which prevent the formation of a nucleus at the boundary of the domain. Then, instead of viscoelastic dissipation, which allows one to close the problem in terms of strain variables only, we use environmental viscous friction formulated in terms of displacements. Note that the latter represent nonlocal functions of strains. Finally, we performed our numerical experiments at small but finite inertia, which apparently also contributed to the inhibition of collective nucleation. As we show in what follows, the reduction of friction leads to re-emergence of the nucleation peak; however, it will then represent a different physical phenomenon vis-à-vis the nucleation peak in a purely viscoelastic overdamped system.

One can see that transformational yield in Figure 1 takes the form of a sequence of similar instability events, each one representing a transition from a marginally stable state to the nearest meta-stable state. Then, in the limit $\delta \to 0$, the yielding emerges as a collection of infinitely many events that are all infinitely small. One can say that in this limit all of these events merge together, and the system can be viewed as remaining all the time in a marginally stable state. The emerging coarse grained macroscopic dissipative potential is known to be a (non-Onsagerian) homogeneous function of degree one which replaces the classical quadratic Onsagerian dissipative potential operative at the level of individual mass points [122,125].

The power spectrum (22) computed in this regime is shown in Figure 3. At high frequencies, it exhibits a range of the power law behavior

$$PS(f) \sim 1/f^2. \qquad (27)$$

There is an almost flat segment at small frequencies with perhaps a single characteristic frequency expressed in a more pronounced way.

To explain these observations, we turn to Figure 2 showing the structure of two consecutive avalanches. Observe that each avalanche can be approximated by an exponentially decaying function

$$R(t) = R_0 e^{-\lambda t}, \qquad (28)$$

where $\lambda$ is the decay rate which appear to be the same for both avalanches. The Fourier transform of equation (28) is given by the formula

$$\hat{R}(f) = \frac{R_0}{\lambda + if} \qquad (29)$$

and the corresponding power spectrum takes the form

$$PS(f) = \frac{R_0^2}{\lambda^2 + f^2}. \qquad (30)$$

One can see that both regimes visible in Figure 3 are captured by this formula, which suggests that the computed power spectrum is a representation of a sum of almost identical individual avalanches of
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Cumulative probability distribution representing statistics of avalanches in the overdamped regime which is shown in Figure 1. Here, $C = 80$ and $N = 1000$. Inset: several consecutive avalanches occurred in a short time interval.</figcaption></figure>

exponential shape. Note also that the transition between high- and low-frequency regimes is positioned around the frequency characterizing the decay rate of the exponential tail. The presence of a characteristic frequency in the low-frequency range may be the sign of the superimposed 'ringing' in the system due to small but non-zero inertia in the system.

Finally, in Figure 4, we show the cumulative probability distribution (CPD) of avalanche sizes represented by the variable $E$. One can see that the distribution is localized on a small interval of avalanche sizes with an exponential cut off tail. This suggests that in the overdamped regime, avalanches are over-correlated vis-à-vis a potential power law scaling regime with individual avalanches being both almost equidistant and having almost the same size. We can conclude that in the ovedamped regime small inertia is almost invisible, playing the role of a quenched quasi-Gaussian disorder, see also [124].

## 2.2. Underdamped regimes

The fully undamped limit $C \to \infty$ with neglected NNN interactions ($G=0$) was studied analytically and numerically in [126]. In this case, there is no hysteresis and in cyclic quasi-static loading the system exhibits, after a one-cycle transient, a reversible behavior reminiscent of entropic elasticity. This means that already during the first cycle, the system fully thermalizes after the system departs from the homogeneous elastic branch. In other words, a state of equipartition is reached in the course of the development of the associated spinodal instability. After reaching the state of thermal equilibrium, the system behaves under quasi-static driving as an adiabatically loaded thermoelastic body. Our numerical experiments fully corroborate these analytical and numerical findings.

The response shown in Figure 5, where $C = 100000$, is very close to the behavior of the system with $C = \infty$. The latter is described by the equation

$$ \delta\ddot{u}_i = (\tilde{\phi}'((\tilde{u}_{i+1} - \tilde{u}_i)/\delta) - \tilde{\phi}'((\tilde{u}_i - \tilde{u}_{i-1})/\delta) + (G/\delta^3)(\tilde{u}_{i+2} + \tilde{u}_{i-2} - 4\tilde{u}_{i+1} - 4\tilde{u}_{i-1} + 6\tilde{u}_i), \quad (31) $$

where the time scale is chosen differently than in the case $C \neq 0$. Specifically, we assumed that the time scale $\tau$ and the dimensionless loading rate $V$ are now defined by the formulas $\tau = \sqrt{\rho/\kappa L}$ and $V = v\sqrt{\rho/\kappa}$. We again started with a homogeneous state where all springs were in a single well and then loaded the system quasi-statically (and effectively isentropically) in a hard device. According to Figure 5, the originally 'cold' system first deforms homogeneously while maintaining an affine configuration. One can see,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of an undamped discrete chain subjected to quasi-static loading in a hard device: the first loading cycle. Here $N = 1000$ and $C = 100,000$. Inset: a zoom on the kinetic energy evolution in the pre-transition time interval.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Stress-strain response of an undamped discrete chain subjected to quasistatic cyclic loading in a hard device. Here $N = 1000$ and $C = 100000$.</figcaption></figure>

however, the presence of superimposed small elastic 'ringing' because our loading is still characterized by small but finite rate. When the system reaches the elastic instability limit, the elastic homogeneous state breaks down giving rise to a complex dynamical regime.

As loading continues beyond the ensuing massive instability, we observe a quasi-deterministic average stress-strain response with superimposed chaotic fluctuations, see Figure 6. In accordance with [126] all modes acquire the same energy which allows one to introduce the idea of temperature. The subsequent quasi-thermoelastic, hysteresis-free averaged behavior proceeds in full agreement with the formulas of classical thermodynamics describing the corresponding adiabatic processes. See [126] for a fully explicit analysis in the case of bi-quadratic NN potential and $G=0$ when both free energy and entropy can be computed analytically.

According to Figure 6, at $C = 100,000$, the overall behavior of the system is practically indistinguishable from the limiting case $C = \infty$. While in the limiting Hamiltonian regime, the actual dissipation vanishes, in the analysis of the underdamped regime, we can still resort to the direct study of the quantity

$$V(t) = \sqrt{\sum_{i} (\dot{u}_{i})^{2}}. \qquad (32)$$

where summation is over all the mass points. Therefore, the statistical study of environmental dissipation in the underdamped system is somewhat similar to the study of temperature fluctuations in the limiting undamped limit.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Power spectrum of the signal $V(t)$ for the undamped discrete chain subjected to quasi-static cyclic loading in a hard device, which is shown in Figure 6. Here $C = 100000$ and $N = 1000$.</figcaption></figure>

The corresponding power spectrum (22), which is not accessible by purely thermodynamic analysis, is shown in Figure 7. We observe fully uncorrelated equilibrium fluctuations at small frequencies producing a flat part of the spectrum $PS(f) \sim 1/f^\alpha$ with $\alpha = 0$. This indicates that long-wave collective modes are effectively subjected to white noise. At larger frequencies, the power law behavior persists with exponent stabilizing at the value $\alpha = 2$. Such behavior is indicative of the underlying Brownian motion of individual mass points which can be interpreted as a Brown noise. It suggests that these points are subjected to fully uncorrelated random forces which should be expected in the state of thermal equilibrium. More generally, the Lorentz-like broad structure of the computed power spectrum is indicative of an exponential decay of correlations in the time domain which in turn points to fast relaxation times characteristic of strongly chaotic systems involving a wide range of frequencies.

Even though we largely lose correlations and intermittency in the underdamped regime, some analog of avalanche distribution can be still constructed from the obtained time series by appropriate thresholding. In Figure 8, we show the distribution of the implied (quasi) avalanches where the quantity $V^2(t)$ was used as the analog of dissipation. According to Figure 8, the resulting distribution is close to Gaussian with an almost flat behavior at small event sizes. We only show here the distribution beyond the lower bound $x_{min}$ determined by minimizing the Kolmogorov-Smirnov distance. At intermediate range, we observe a small value of the exponent $\beta \approx 0.23$, which suggests that inertia leads to over-correlation among quasi-avalanches. At larger event sizes, we observe a characteristic exponential cut off reflecting the finite size of the system.

To summarize, while we have seen that in the overdamped regimes, the dynamics is close to being fully deterministic and regular, which suggests very limited complexity, in the underdamped regimes, the dynamics is close to being fully chaotic, which again means that complexity is minimal. As we show below, more interesting structure of correlations emerges in the intermediate regimes where the system in neither overdamped nor underdamped.

## 2.3. Intermediate regimes

At finite values of the parameter $C$, the direct integration of the system (12) subjected to quasi-static driving produces intermittent dynamics with avalanches of widely different sizes. Moreover, as we show below, in a well-defined interval of the values of $C$, the distribution of avalanche sizes exhibits a range of power law behavior indicating scale-free internal organization of dynamic microstructure.

A characteristic stress-strain response of the system at $C = 8000$ under cyclic loading over 20 cycles is shown in Figure 9. In each cycle, we observe recurrent system size characteristic event when initially homogeneous (affine) configuration is breaking down due to a massive elastic instability. Using the language of plasticity theory, one can say that the system undergoes brittle yielding [13,57].
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Cumulative probability distribution for the magnitude of the effective energy dissipation during individual avalanches in the undamped discrete chain subjected to quasi-static cyclic loading in a hard device, which is shown in Figure 6. Inset shows a fragment of the actual time series. Here $C = 100,000$ and $N = 1000$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Stress-strain response of a quasi-statically driven chain with inertia and dissipation balanced at $C = 8000$. Here $N = 1000$.</figcaption></figure>

According to Figure 9, after the first system-size avalanche, the subsequent avalanches are all of smaller size, and in each cycle, the system apparently reaches a steady-yield regime. While the average stress in this regime is maintained at an almost constant level, we observe a broad, heavy-tailed distribution of stress drops.

The computed power spectrum is shown in Figure 10, where one can identify three regimes corresponding to small, intermediate, and large frequencies. In the range of small frequencies, the spectrum is flat, indicating the absence of persistent system size collective correlations. The characteristic times associated with such frequencies, which would have characterized the durations of the corresponding correlated events, are much bigger than the durations of the largest avalanches. It suggests that this range has nothing to do with intermittency, and the fluctuational response is reminiscent of what we have seen in the underdamped regime. In the range of large frequencies, we again see the Lorenzian behavior, with $1/f^2$-type decay, which is similar to what we have seen in the overdamped regime. This is apparently due to the dominant role played in this range by the time scale of viscous dissipation, and it reflects the relaxational behavior of individual masses at the end of each avalanche.

The most interesting behavior is observed in the range of intermediate frequencies, where we see a power law decay with a nontrivial exponent $\alpha = 1$ characteristic of a so-called 1/f noise [114–116,119,127, 128]. The associated frequencies correspond roughly to the durations of avalanches whose magnitudes are inside the ‘inertial’ (power law) range, see below. Note that a similar 1/f noise has been also recorded in the closely related discrete models of crystal plasticity [117,124]. If small avalanches die off in the intermediate regime, in the range of frequencies where $PS(f) \sim 1/f$, such a decay is practically invisible
<figure><img src="image_13.png" /><figcaption>Figure 10. Power spectrum of the underdamed system shown in Figure 9. Here C = 8000 and N = 1000.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Cumulative probability distribution in log-log scale of the energy dissipated by avalanches in the underdamed system loaded quasi-statically in a hard device and shown in Figure 9. Here, C = 8000 and N = 1000. Inset: a fragment of the actual time series.</figcaption></figure>

because at these scales, dissipation is almost negligible, and the structure of avalanches is shaped by dynamics as much as by dissipation.

Finally, we show in Figure 11 the associated statistics of avalanche sizes represented through distribution of the dissipated energies E. Here, we zoomed into a robust range of power law behavior

$$P(E) \sim E^{-\beta} \qquad (33)$$

stretching over two decades with exponent $\beta \approx 2.1$. In the inset, we show the structure of the corresponding intermittent signal. Our Figure 12 compares this distribution with the distribution of the avalanche magnitudes S which has again a power law form

$$P(S) \sim S^{-\kappa}, \qquad (34)$$
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Statistical signatures of avalanche behavior for the system shown in Figure 9: (a) Cumulative probability distribution of avalanche durations, (b) cumulative probability distribution of avalanche sizes, and (c) cloud plot showing joint distribution of avalanche sizes and durations and revealing strong correlation between these two quantities.</figcaption></figure>

now with the exponent $\kappa \approx 2.2$, see Figure 12(a), and the distribution of avalanche durations also exhibiting a power law behavior

$$P(T) \sim T^{-\tau}, \qquad (35)$$

now with the exponent $\tau \approx 2.5$, see Figure 12(b). The revealed ubiquity of power law behaviors corroborates the idea that the system develops in this range robust scale-free correlations.

It has been shown in [115] that there exists a relation between the exponents representing the power spectrum and the exponents characterizing the power law tail of the probability distribution for avalanche sizes of the form

$$PS(f) = f^{-(3-\kappa)/\xi}. \qquad (36)$$

Here, the parameter $\xi$ can be found from another scaling relation

$$\langle S \rangle \sim \langle T \rangle^{1/\xi}, \qquad (37)$$

and our cloud plot juxtaposing distributions of sizes and durations, shown in Figure 12(c), suggests that in our case $1/\xi \approx 1.2$. Given that we also obtained the values $\kappa \approx 2.2$ and $\alpha \approx 1$, the implied universal relation

$$\alpha = \frac{3 - \kappa}{\xi} \qquad (38)$$

is respected. We also note that the exponents of the power laws for durations $\tau$ and energies $\beta$ can be related to $\kappa$ and $\xi$ through another two scaling relations discussed in [115], [114], [118]:

$$\tau = 1 + \frac{\kappa - 1}{\xi} \qquad (39)$$

and

$$\beta = 1 + \frac{\kappa - 1}{2 - \xi}. \qquad (40)$$

Using the computed values of the exponents, we again find an agreement which points to consistency of our numerical results.

We highlight the emergence in Figure 12 of a feature which can be interpreted as a supercritical behavior in the range of large events: it may not look as a conventional 'bump' because instead of PDF, we present here a CPD. This feature reveals the presence of recurrent large 'brittle' nucleation events that occur in each cycle, which are clearly visible in the stress-strain curves shown in Figure 9. We reiterate that supercritical 'bumps' have been found to be characteristic of underdamped behavior in many other systems [92,93,96-98]. In our representation, these large events present themselves as almost a single data point because in the CPDs shown in Figure 12(a)-(c), we do not perform any binning.

Finally, as a word of caution, we stress that our three numerically observed fluctuation patterns still need to be subjected to careful statistical analysis as it is usually difficult to distinguish between different
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Regime diagram in the space of parameters showing the stress-strain response of the quasi-statically driven chain in underdamped, overdamped and intermediate regimes. Here $M = v^2\rho/\kappa$ and $V = \gamma v/(L\kappa)$. The boundaries between different regimes are delineated approximately in accordance with multiple numerical tests conducted at different values of the parameter $C = M/V^2$. The particular regimes a,b and c are discussed in the above in more detail.</figcaption></figure>

types of heavy-tailed distributions on finite intervals and from noisy data. In particular, power law and log-normal distributions are often confused with the latter sometimes interpreted erroneously as power laws with exponential decay [129]. We plan to address this technical question in a separate study.

## 2.4. Regime diagram

The goal of this section is to consolidate and summarize the results obtained so far. Figure 13 presents the regime diagram illustrating qualitatively different responses of the system in the space of dimensionless parameters

$$M = v^2 \rho / \kappa, \quad V = \gamma v / (L\kappa). \qquad (41)$$

Note that we have chosen both parameters $M$ and $V$ to be dependent on the rate of loading $v$. In this way, we can highlight that our assumption of quasi-static loading implies a double asymptotics

$$M \to 0, \quad V \to 0. \qquad (42)$$

Therefore, in Figure 13, we are first of all interested in the behavior of the system around the origin of the parameter space. In view of non-commutative nature of the limits (42), the behavior of the system turns out to be dependent on the asymptotic path controlled by the value of the parameter

$$C = M/V^2. \qquad (43)$$

In other words, qualitatively different fluctuational patterns can be expected to take place as one changes the relative strength of inertia versus the strength of dissipation. But this is exactly what we have seen in our numerical experiments which were all conducted in a broad range of values of $C$ even though we reported above the results only for selected values of this parameter.

According to Figure 13, there are three main regimes with qualitatively different behavior. As we have seen, in two of them, an overdamped/damped regime (small $C$) and an underdamped/undamped regime (large $C$), the behavior of the system is relatively simple and well understood. It is either uncorrelated
(disorder) or overcorrelated (order) behavior and therefore in none of these regimes the system shows *complexity exemplified by intermittency and scaling*. In the intermediate regime, which is identified in Figure 13 by the word ‘scaling’, the inertia and the dissipation are balanced in the sense that none of them dominates. Our numerical experiments showed that distributions with a robust scaling range over two decades are observed over a finite range of values of the parameter *C*. The characteristic feature of the scaling regime is a delicate interplay between *inertia, dissipation, and discreteness* allowing the system to self-organize towards scale-free non-equilibrium steady state.

## 2.5. Extended criticality

While the interval of parameters *C* where scaling was observed is relatively narrow and while the power law response associated with this interval extends only along a finite range of event sizes, here we argue that one can still interpret the observed behavior as self-organization towards ‘generic’ quasi-critical regime.

One reason to think that we deal here with extended criticality is the near *marginality* of the underlying mechanical system. We have seen that this is definitely the case for the overdamped regime, where the system in a large *N* limit exhibits a quasi-plastic yield on a stress–strain plateau. It is clear that the yielding takes place because in the overdamped regime, the system explores only a part of the energy landscape where the barriers, separating different metastable states, are extremely small. Then, in the scaling regime, the system can be thought as being on one side close to barrierless marginal stability limit and, on the other side, being able to easily cross the remaining energy barriers due to underdamping. Apparently, it is the implied delicate balance between near marginality with its abundant elastic instabilities and the ‘fluidity’ of navigating among small energy barriers which allows the system to self-organize toward a dynamical critical state. While the resulting yielding takes place around a fixed average stress, the ‘tuning’ is only apparent because stress is not a control parameter in the hard device loading experiments.

In this perspective, the emergence of ‘generic’ scale invariance can be also associated with what has been coined in [130], see also [131], [132], and [50], as ‘sweeping’ of an instability. If we interpret transformational yield as such an instability, we see that the threshold is indeed continuously swept due to the presence of a feedback provided by long-range elastic interactions. In particular, a hard device loading ensures that the threshold is overshot with subsequent stress drop bringing the system back into a stable state. Then the system is destabilized again and so on. Various prototypical models of the implied feedback allowing the system to maintain near marginal stability, and safeguarding the ‘generic’ nature of the scaling behavior, were proposed in [124] and [57].

The fact that the observed dynamic activity with a broad distribution of scales takes place generically, without any special tuning of parameters, can be also interpreted as SOC. In this respect, it may be appropriate to mention that our account for NNN interactions is a poor man’s attempt to imitate the effect of demagnetizing forces in the theory of Barkhausen effect, which is crucial for reaching a robust critical regime [52,133,134]. Our NNN interactions can be also viewed as imitating strain incompatibility between austenite and martensite which is another important factor of reaching robust criticality [5,49]. Specifically, the presence of NNN interactions contributes to the creation of metastability by preventing most of the avalanches from sweeping over the whole system. The classical analogs of the implied metastability are the stable critical nuclei in isochoric thermodynamic systems [135], see also [136].

Yet another factor which may explain the ‘generic’ nature of scaling in our model is the presence of inertia-induced dynamic disorder. While the quenched disorder in RFIM had to be tuned to reach the critical state, in our model, the effectively annealed dynamic disorder is both self-induced and self-tuning. One can then argue that the implied feedback loops regulate the level of such disorder driving the system to criticality without external fine tuning.

Finally, we recall that bringing inertia into the model leads to softening of the kinetic relation which regulates propagation of phase boundaries [73,125]. As we have already mentioned, such softening is an important source of intermittent stick-slip behavior and a factor often used to explain the robust emergence of ‘generic’ scaling [137–139]. While in our oversimplified model, the interfaces between phases are not explicitly tracked, the kinetic softening is still present, as it can be shown through the analysis of radiative damping by propagating phase transition fronts [125,140].
### 3. One-dimensional continuum model

Given that in a discrete model we were able to reach the regime of extended scale-free behavior by adjusting the value of a single parameter C, it is tempting to argue that scaling in martensites is entirely due to an interplay between inertia and dissipation. It is then natural to check if this conjecture survives in the corresponding continuum model. To this end, we need to perform in equation (12) a formal limit $\delta \to 0$. The resulting dimensionless partial differential equation is of the form:

$$ C\ddot{u} = \partial_x\phi'(\partial_x u) + \partial_{xx}\dot{u} - G\partial_{xxxx}u. \quad (44) $$

Here, $u(x, t)$ is the continuum displacement, $\phi$ is our double-well potential, and the only dimensionless parameters are the ratio of inertia over dissipation C and the measure of nonlocality G. The corresponding total energy is again $E = \mathcal{H} + \mathcal{F}$, where the (dimensionless) kinetic energy is

$$ \mathcal{H} = \int \frac{C}{2} \dot{u}^2(x, t) dx. \quad (45) $$

The (dimensionless) elastic energy has the form:

$$ \mathcal{F} = \int \left( \phi(e) + \frac{G}{2} (\partial_x e)^2 \right) dx, \quad (46) $$

where we introduced the continuum strain $e = \partial_x u$ which plays in this case the role of the classical Landau–Ginzburg order parameter. The only difference from equation (12) is that in equation (44), instead of environmental friction, we used the (dimensionless) viscoelastic Rayleigh dissipative potential

$$ \mathcal{R} = \int \sum_{i}^{3} \frac{1}{2} \dot{e}^2(\mathbf{r}, t) d^2x. \quad (47) $$

This choice is dictated by the necessity to eventually extend the continuum model from 1D to 3D, where the environmental viscosity would not make much sense.

The 1D continuum model (44) has been studied extensively in both statics and dynamics, see for instance [73] and the references therein. Here, we used (44) to perform numerical experiments with a system subjected to quasi-static cyclic loading.

Our main conclusion is that, independently of the value of parameter C, the model (44) does not produce either intermittency or scaling. The reason is a rather limited ability of such continuum model to support pinning of phase boundaries and, more generally, to generate elastic metastability which is necessary for the emergence of the recurrent quiescence and the associated intermittency. One can say that with discreteness lost, the effective energy landscape becomes too simple which prevents the emergence of (near) marginality and prevents the system from reaching sufficiently complex non-equilibrium steady state.

It may sound counter intuitive, but in this respect, the finite dimensional discrete system is much richer than its straightforward infinite dimensional analog. Indeed, since springs can transform only after reaching individual thresholds, the effective energy landscape of a discrete system has an exponentially large number of local minima and, while this number diverges in the continuum limit, the minima themselves become progressively more shallow and eventually disappear [121,141–143]. In particular, the continuum model (44) loses the capability of trapping phase boundaries which compromises the hysteretic nature of martensitic transformations.

In the next two sections, we study whether the problem can be remedied by moving from an over schematic 1D continuum model, towards its 2D and 3D continuum analogs which can already represent realistic martensitic transformations.
## 4. Two-dimensional continuum model

The goal of this section is to analyze a straightforward 2D analog of equation (44). To focus exclusively on the effect of higher dimensionality, we have chosen to intentionally simplify the targeted martensitic transformation by neglecting its volumetric effect. This makes the energy wells, describing the high-symmetry austenite and the low-symmetry martensite, automatically rank one connected and therefore, kinematically compatible. By neglecting this important source of metastability in the system, we naturally place our 2D model between the automatically compatible 1D model and the generically incompatible 3D model [107,144–147].

A simplistic but nevertheless meaningful example of a martensitic phase transition in 2D involves square and rectangle phases [72,112,148]. To describe the corresponding transformation, we first turn to dimensional variables and again write the total energy in the form $E = \mathcal{K} + \mathcal{F}$, where the kinetic energy has the standard form

$$ \mathcal{K} = \int \frac{\rho}{2} \dot{u}^2(\mathbf{r}, t) d^2x, \qquad (48) $$

and where $\rho$ is the mass density and $u(\mathbf{r}, t)$ is the displacement vector. The elastic energy must be specialized to reflect the chosen crystal symmetry and we write it in the form

$$ \mathcal{F} = \int \left( \frac{A_1}{2} e_1^2 + \frac{A_3}{2} e_3^2 + \phi(e_2, \tau) + K e_1 e_2^2 + \frac{G}{2} |\nabla e_2|^2 \right) d^2x. \qquad (49) $$

Note first that the elastic energy density in equation (49) is presented as a function of three components of the linear elastic strain tensor

$$ e_1 = \left( \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \right), e_2 = \left( \frac{\partial u_x}{\partial x} - \frac{\partial u_y}{\partial y} \right), e_3 = \frac{1}{2} \left( \frac{\partial u_y}{\partial x} + \frac{\partial u_x}{\partial y} \right). \qquad (50) $$

As it is clear from equation (49), we assumed that the strain component $e_2$ is the primary order parameter. The first two harmonic terms in equation (49) with the coefficients $A_1$ and $A_2$ describe the classical physically linear elasticity in a solid with rectangular symmetry. The cubic term in equation (49) is needed to potentially introduce a nonzero volumetric effect of the phase transition and therefore the third order elastic modulus $K$ can be viewed as controlling the strength of the corresponding nonlinear shear-dilatation coupling [149,150]. However, as we have already mentioned, in our numerical experiments we set $K$ to be very small and the corresponding term is mentioned here only because it will reappear in the 3D model. We further assumed that the regularizing gradient energy term in equation (49), mimicking NNN interactions in the discrete model, includes only the gradient of the order parameter $e_2$ with the coefficient $G$ bringing into the continuum model a finite internal length scale. To complete the model, we introduce viscoelastic Rayleigh dissipative potential:

$$ \mathcal{R} = \int \sum_{i}^{3} \frac{\gamma_i}{2} \dot{e}_i^2(\mathbf{r}, t) d^2x, \qquad (51) $$

where $\gamma_i$ are the corresponding effective viscosity coefficients.

The numerical implementation of the equation of motion requires spatial and temporal discretization. If we use as a spatial scale, the grid size $d_0$, the dimensionless spatial coordinates $\tilde{x}$ and displacements $\tilde{u}$ become

$$ \tilde{x} = \frac{x}{d_0}, \quad \tilde{u} = \frac{u}{d_0}. \qquad (52) $$

By introducing the characteristic time scale $t_0$, we similarly normalize time

$$ \tau = \frac{t}{t_0}. \qquad (53) $$
<figure><img src="image_17.png" /><figcaption>Figure 14. The elastic energy density $\phi(e_2)$ at different values of the temperature $\tau$.</figcaption></figure>

We also introduce the energy density scale $f_0$ which allows us to write the remaining dimensionless parameters of the model in the form:

$$ \tilde{\rho} = \frac{\rho d_0^2}{t_0^2 f_0}, \quad \tilde{G} = \frac{G}{d_0^2 f_0}, \quad \tilde{\gamma} = \frac{\gamma}{t_0 f_0}, \quad \tilde{A}_i = \frac{A_i}{f_0}. \qquad (54) $$

Finally, the minimal Landau-type energy density function $\phi$ capturing the schematics of square-to-rectangle transition is written in the form

$$ \phi(e_2, \tau) = \frac{\tau}{2}e_2^2 - e_2^4 + \frac{e_2^6}{2}, \qquad (55) $$

where $\tau$ is the dimensionless temperature. In Figure 14, we show that, for instance at $\tau = 1$, which is the point of an equilibrium first-order transition, the energy density (55) exhibits three equivalent minima, describing two variants of the low symmetry martensite (rectangular phase) phase and one minimum describing the high symmetry austenite (square phase); below this temperature, the austenite stops being the ground state.

In our numerical experiments, the dynamic equations were discretized on a rigid grid of size 512 × 512 with periodic boundary conditions and solved using a Fourier pseudo-spectral spatial scheme with the corresponding time marching temporal algorithm detailed in Appendix 1, see also [151]. Since we study the 2D model only, to provide qualitative illustrations of the effects of inertial dynamics, we used generic values of dimensionless parameters. For instance, throughout this study we assumed that $\tilde{G} = 5$, other parameters are specified below. A physically meaningful calibration of the model is postponed till the next section where we consider a fully realistic 3D model of the same basic type.

Consider first an underdamped regime with $\tilde{\rho} = 1$ and $\tilde{\gamma}_2 = \tilde{\gamma}_1 = \tilde{\gamma}_3 = 0.5$. This choice corresponds to the limit of intermediate values of $C$ in the 1D model. In Figure 15, we illustrate the results of our numerical experiments describing dynamic nucleation of a martensite. At time $t = 0$, a homogeneous configuration of austenite (square) phase was equilibrated at $\tau = 1$. To destabilize this configuration, we decrease the dimensionless temperature $\tau$ in $10^{-4}$ increments till the transformation begins.

To break geometrical degeneracy and ensure controlled nucleation of the martensite, we placed a small martensitic embryo in the middle of the square phase. Specifically, the embryo was installed by initiating the value of the order parameter to 1 in small circular domains in the middle of the computational domain, see red regions in Figure 15(a). Our Figure 15(a)–(f) illustrate different stages of the evolution of the transformation starting from the initial nucleation event. Green color represents the undeformed austenite phase where $e_2 = 0$. Red and blue colors correspond to two variants of the emerging martensite.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Time evolution of a pre-existent martensitic embryo in an underdamped 2D continuum model. Panels (a)–(f) illustrate different stages of the evolution. Green color represents the austenite, red and blue colors represent the martensite, see Figure 14 for the identification of particular variants.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Time evolution of a pre-existent martensitic embryo in an almost undamped 2D continuum model. Panels (a)–(f) illustrate different stages of the evolution. Green color represents the austenite, red and blue colors represent the martensite, see Figure 14 for the identification of particular variants.</figcaption></figure>

As we see in Figure 15(a), the transformation begins around the infinitesimal embryo in the form of a single variant of the martensite effectively chosen by the type of the inserted imperfection. Already in Figure 15(b), we see that the homogeneous martensitic nucleus evolves into a complex multiscale two-variant twinned microstructure which grows dynamically in an apparently self-similar way till it hits the boundaries of the domain, see Figure 15(c). As the transformation process unfolds, interaction of internally generated elastic waves produces multiple new micro-twinned patches appearing almost spontaneously all over the computational domain, see Figure 15(d). As such multiscale microstructure formation process reaches the boundaries of the domain, austenite phase largely disappears, see Figure 15(e). Eventually, the microstructure stabilizes by coarsening and all independently developing transformed zones finally merge into a single multiscale texture, see Figure 15(f). The resulting complex variant mixture contains differently oriented but geometrically compatible martensite laminates forming an intricate hierarchical pattern. In view of the geometric compatibility of all participating energy wells, the formation of microstructure of such complexity can be viewed as associated with an interplay between small dissipation and inertial dynamics. It implies, generation of elastic waves by the transformation events, interaction of these waves with inhomogeneities, their reflection from the boundaries and, most importantly, their self focusing which triggers secondary transformation events.

Consider next the almost undamped regime with $\tilde{\rho}=1$ and $\gamma_2 = 1$, $\gamma_1 = \gamma_3 \approx 0$. As we see in Figure 16(a), the transformation begins again in a form of a single variant of the martensite effectively chosen by the nature of the inserted imperfection. Already in Figure 16(b), we see that the growing martensitic nucleus splits into two variants and then extends as a lenticular domain. The latter contains two twinned martensitic variants forming fully compatible interfaces between themselves and with the austenite, see Figure 16(c). As this initial band hits the boundaries of the domain, the growth process continues along the direction perpendicular to its boundaries involving thickening of the two already formed martensitic variants. The presence of self-generated dynamic activity compromises such purely coarsening process producing secondary nucleation events and generating sequential appearance of zipping martensitic variants with alternating strains, somewhat similar to what has been predicted in [101]. The resulting fine poly-twin pattern, which emerges only after a very long time, is regular, and the whole pattern shows much less complexity than the multiscale texture observed in the underdamped regime.
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Post-nucleation pattern in a 2D almost fully damped continuum model, see Figure 14 for the identification of the colors with particular variants of the martensitic phase.</figcaption></figure>

Finally, consider the *damped* regime with $\tilde{\gamma}_2 = 1$, $\tilde{\gamma}_1 = \tilde{\gamma}_3 = 0.5$ and $\tilde{\rho} \approx 0$. This is equivalent to the assumption that $C \approx 0$ in equation (44). In this limit, we obtain time-dependent Ginzburg–Landau model with viscous damping, see for instance [72,101]. Then, the breakdown of an unstable austenite phase results a predictable emergence of a coarse highly regular twin microstructure, see Figure 17. The whole process can be described as an elastic spinodal decomposition with subsequent coarsening driven by weak interaction of the twin boundaries with the possibility of being stuck in one of the kinetically induced locking configurations [152,153].

The examples presented in this section show unambiguously that suppression of inertia in this type of problems blocks the emergence of complexity. Also, as our more exhaustive study have shown, the metastability problem detected in the 1D continuum model persists in the 2D compatible problem. Specifically, we have checked that, independently of the degree of over- or under-damping the ensuing dynamics does not support either intermittency or scaling. Still, the relative transparency of the associated 2D problem allowed us to illustrate in some detail the dramatic difference between the processes of microstructure formation in overdamped and underdamped models. In particular, we could show that bringing inertia into the model contributes to multiscale nature of the emerging microstructure through dynamic creation of virtual nucleation sites. However, while this continuum model generates a succession of progressively finer scales, we did not see any evidence that it supports robust self-organization towards criticality.

**5. Three-dimensional continuum model**

In this section, we show that the desired complexity of the energy landscape, comparable with what we have seen in our discrete model and, therefore, compatible with scaling and extended criticality, can be recovered if we account in our continuum modeling for the incompatibility of the energy wells. Using 3D instead of 2D framework is beneficial in this respect as it allows one to perform a comparison with actual physical experiments. It also considerably increases the set of accessible data which ensures the detection of scaling regime with higher statistical certainty.

Specifically, we consider, in this section, a realistic continuum model of a cubic-to-tetragonal martensitic transition in 3D which, in particular, takes into account a nonzero volumetric effect of the transformation. The goal is to reproduce in our numerical experiments the result of [63], where the authors investigated the acoustic activity which was measured experimentally as a single crystal of $Fe_{68.8}Pd_{31.2}$ was quasi-statically driven through the cubic-to-tetragonal phase transition. As we
have already mentioned, the analysis of statistical distributions of dissipated energy and avalanche durations performed in [63], revealed a power law statistics and produced particular values of the exponents.

While the kinetic energy in the desired 3D continuum model can be again taken in the form (48), the elastic free energy must be adjusted to capture the specifics of a particular martensitic transformation. Specifically, to describe a generic cubic-to-tetragonal martensitic transition in 3D, we use the following expression for the Ginzburg-Landau type energy [113,148,154,155]

$$ \mathcal{F} = \int (\phi_L + \phi_G)d^3x \qquad (56) $$

where the energy densities describing local and gradient contributions are given by the formulas:

$$ \phi_L = A_2(e_2^2 + e_3^2) + A_3(e_4^2 + e_5^2 + e_6^2) + A_4 e_3(e_3^2 - 3e_2^2) + A_6(e_2^2 + e_3^2)^2 + A_1(e_1 - K(e_2^2 + e_3^2))^2, \qquad (57) $$

and

$$ \phi_G = \frac{G}{2}(|\nabla e_2|^2 + |\nabla e_3|^2). \qquad (58) $$

One can see that we use here two coupled primary-order parameters represented by two deviatoric components of the linear elastic strain tensor $\mathbf{e} = (1/2)(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ and describing shear deformations in $110 < \bar{1}10 >$-type shear deformations:

$$ e_2 = \frac{1}{\sqrt{2}}(\epsilon_{xx} - \epsilon_{yy}), \quad e_3 = \frac{1}{\sqrt{6}}(\epsilon_{xx} + \epsilon_{yy} - 2\epsilon_{zz}). \qquad (59) $$

The remaining non-order parameter components of the strain tensor are:

$$ \begin{aligned} e_1 &= \frac{1}{\sqrt{3}}(\epsilon_{xx} + \epsilon_{yy} + \epsilon_{zz}), \\ e_4 &= \epsilon_{xy} + \epsilon_{yx}, e_5 = \epsilon_{xz} + \epsilon_{zx}, e_6 = \epsilon_{yz} + \epsilon_{zy}. \end{aligned} \qquad (60) $$

While the coefficients $A_i$ in equation (56) with $i=1, 2, 3$ are the classical linear elastic moduli of a tetragonal phase, the coefficients $A_4$ and $A_6$ describe the coupling between the order parameters which is necessary to destabilize the cubic phase in favor of the tetragonal phase [113,148,154,155]. As we have already seen in the 2D case, adding of the coupling coefficient $K$ allows one to capture shear-induced volumetric effect of the transformation. The chosen form of the energy density (57) guarantees that there are in general three (geometrically compatible) energy wells corresponding to symmetry-related tetragonal variants of the martensitic phase, and one (geometrically incompatible) energy well corresponding to a higher symmetry austenite phase.

## 5.1. Calibration of the model

Fixing parameters of the energy density (56) requires the knowledge of such experimentally measured quantities as the homogeneous transformation strain, the elastic constants, and the interfacial energy. Fortunately, all of these parameters have been experimentally measured for FePd alloys close to the transition temperature [113].

Specifically, it is known that $a = 3.725Å$ and $c = 3.795Å$ are the lattice constants of the equilibrium martensite (tetragonal phase) at stress-free state while $a_0 = 3.756Å$ is the lattice constant of the unstrained austenite (cubic phase). The corresponding stretch tensors mapping the cubic lattice to three symmetry-related variants of the tetragonal phase can be then written in the form:

$$ U_1 = \begin{pmatrix} c^* & 0 & 0 \\ 0 & a^* & 0 \\ 0 & 0 & a^* \end{pmatrix}, \quad U_2 = \begin{pmatrix} a^* & 0 & 0 \\ 0 & c^* & 0 \\ 0 & 0 & a^* \end{pmatrix}, \quad U_3 = \begin{pmatrix} a^* & 0 & 0 \\ 0 & a^* & 0 \\ 0 & 0 & c^* \end{pmatrix}, \qquad (61) $$
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Typical energy landscape (contour plot) in the space of primary-order parameters adopted in the 3D continuum model. Three minima correspond to three symmetry related variants of the stable martensite (tetragonal phase).</figcaption></figure>

where $c^* = c/a_0$ and $a^* = a/a_0$. In view of these measurements we know, for instance, that the equilibrium value of the nonzero primary order parameter for, say variant 3, is $e_3^0 = 0.0152$ while the volume change associated with the transformation is $e_1^0 = 0.0072$. From equation (61) one can see that even without the volumetric effect taken into consideration, austenitic and martensitic wells are not geometrically compatible.

The elastic moduli of FePd close to the transition temperature have been also experimentally measured, see the discussion in [156], [157] and [113]. In particular, it was observed that this martensitic transformation is accompanied by softening of the deviatoric elastic modulus of the austenite $C'^A = (C_{11}^A - C_{12}^A)/2$, while the elastic constants $C''^A = (C_{11}^A + C_{12}^A + 2C_{44}^A)/2$ and $C_{44}^A$ vary only slightly, which is all in accordance with our model. Since the elastic moduli of the martensitic phase are largely unknown, it is usually assumed that $C'^M \approx 2C_{44}^M = 2C_{44}^A$, were it is implied that the elastic constants of the martensite are expressed in the undeformed reference state of the austenite.

If we normalize the elastic moduli in each of the phases using the same energy density scale $f_0 = 20.4$ GPa taken from [158], and introduce the dimensionless constants

$$ \tilde{A}_i = A_i/f_0, \quad \tilde{K} = K/f_0, $$

we can use the experimental data from [156] to obtain the following numerical values for the dimensionless parameters of the energy density (57): $\tilde{A}_2 = 1$, $\tilde{A}_4 = -131.57$, $\tilde{A}_6 = 4328.2$, $\tilde{A}_1 = 14$, $\tilde{A}_3 = 2$, and $\tilde{K} = 31.11$. A typical contour plot of the ensuing elastic potential is shown in Figure 18.

Turning to the gradient term (58), we observe that the adopted expression is the simplest one involving only primary-order parameters $e_2$ and $e_3$. For simplicity, we treated both of these parameters equally and introduced a single coefficient $G$ bringing into the problem an internal length scale $d_0$. To fix this parameter, we first introduce the dimensionless coordinates $\tilde{x} = x/d_0$ and displacements $\tilde{u} = u/d_0$. We then use the known dimensional value $G = 3.15 \times 10^{-8} J/m$ obtained from microstructural data in [159] and assume arbitrarily that $G/(d_0^2 f_0) = 5$ which sets the value of the internal length scale at $d_0 = 1.81$nm.
To deal with dissipative dynamics, we again adopt for the Rayleigh potential the simplest viscoelastic form:

$$ \mathcal{R} = \int \sum_{i} \frac{\gamma_i}{2} \dot{e}_i^2(\mathbf{r}, t) d^2x. \qquad (62) $$

where $\gamma_i$ are the associated generalized viscosity coefficients. They can be nondimensionalized using the characteristic time scale $t_0$ of visco-elastic relaxation which we set to be $t_0 \simeq 1$ ps to ensure that in the numerical experiments the duration of realistic avalanches in dimensionless time scale $\tau = t/t_0$ is of order 1. We further assume that the non-dimensional viscosity parameters $\tilde{\gamma}_i = \gamma_i/(t_0f_0)$, describing relaxation of the primary-order parameters $\tilde{\gamma}_2$ and $\tilde{\gamma}_3$, are both equal to 1. This would mean that the corresponding dimensional viscosity coefficients are of the order $\gamma \simeq 21 \times 10^{-3}$ Ns/m$^2$, which is close to the measured damping parameter for V$_3$S alloy also undergoing a cubic-to-tetragonal transition [160]. The remaining viscosity coefficients were chosen to be much smaller: $\tilde{\gamma}_1 = \tilde{\gamma}_4 = \tilde{\gamma}_5 = \tilde{\gamma}_6 = 0.1$.

Finally, we need to fix the degree of underdamping. It is clear that the scale of inertial effects in our model is characterized by the dimensionless parameter

$$ \tilde{\rho} = \frac{\rho d_0^2}{t_0^2 f_0}, $$

which plays the role analogous to the parameter C in the discrete model since $\rho$ represents the strength of inertia while the time scale $t_0$ characterizes the strength of dissipation. Following [113], we use the estimate $\rho \sim 10^4$ kg m$^{-3}$ which suggests that the value $\tilde{\rho} \sim 1$ can be used as representing a typical underdamped regime.

## 5.2. Numerical experiments

Evolution of the system in such regime was studied numerically by solving the governing equations with a time step $\Delta t = 0.01$ and using periodic boundary conditions, see Appendix 1 for details. The simulations were carried out using a grid of size of $256 \times 256 \times 256$, which corresponds to a domain size of $0.46\mu m \times 0.46\mu m \times 0.46\mu m$, where we used the adopted value of $d_0$. Note that the domain, used in the targeted experimental studies of avalanche-type dynamics, was much larger, of order of mm, which suggests that our experiments are most probably affected by the system size and any quantitative comparison with experiment must incorporate the anticipated size effect.

In our numerical experiments, we simulated thermally, rather than mechanically, driven transformation while employing a very small cooling rate. Specifically, we effectively decreased the temperature by reducing the value of $\tilde{A}_2$ using the increments $10^{-5}$ every 1000 time steps. We started at $t=0$ where the state was homogeneous, and the displacements were absent, $\tilde{u}_i(\tilde{\mathbf{r}}) = 0$. To break the degeneracy of such a state, we placed a single defect in the middle of the computational domain, which created a controlled nucleation site. The implied perturbation was included into the energy density through a linear term $U_l(\tilde{\mathbf{r}}) = -s(\tilde{\mathbf{r}})(e_2 + e_3)$, where $s(\tilde{\mathbf{r}})$ is the external stress field acting on primary-order parameters and mimicking a configuration of balanced force couples. One can see that the latter induce a locally tetragonal distortion. Specifically, following [161], we assumed that

$$ s(\tilde{\mathbf{r}}) = \frac{s_0 e^{-|\tilde{\mathbf{r}} - \tilde{\mathbf{r}}_0|^2}}{\zeta^2}, \qquad (63) $$

where $\tilde{\mathbf{r}}_0$ is the defect coordinate and where we set $s_0 = 3$ and $\zeta = 4$.

The unfolding of the transformation processes is illustrated in Figure 19, where we show the evolution of the system from the nucleation moment till the final state is reached under constant cooling rate. As we see in Figure 19(a), the homogeneous austenite (shown as transparent) transforms first into a single variant of martensite (variant 3 shown in blue color) around the defect which obviously biases this particular choice. Due to complex interaction of emitted waves, other variants start to pop up almost immediately
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Time evolution of the system after a nucleation event in the underdamped 3D continuum system. Austenite is presented by transparent gray. Colors represent different martensitic variants. (a) $t = 150$, (b) $t = 300$, (c) $t = 3000$, (d) $t = 4000$ and (e) $t = 16,000$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Time evolution of the dissipated energy after the nucleation event in 3D underdamped continuum model. Panels (a)–(c) represent the same signal at different temporal scales.</figcaption></figure>

(variants 1 and 2 shown in green and red colors, respectively), see Figure 19(b). The texture of internally twinned martensite laminates continue to grow dynamically while progressively generating more and more complex hierarchical pattern, see Figure 19(c) and (d). Finally, the austenite is completely transformed into the variants of martensite which self-organize to form a multi-scale structure, see Figure 19(e). Note that each pair of martensitic variants in the ensuing mixture state is geometrically compatible and therefore forming what is known as compound twins [162].

However, due to incompatibility of the austenitic and martensitic energy wells, the formation of a single variant of martensite inside an austenite matrix is hardly possible since it would have required a forbiddingly large energy cost. Instead, right after the initial instability, we observe the emergence of a complex combination of multi-variant lamellae, incorporating all symmetry-related energy wells of martensitic phase. The system effectively relies on elastic screening and compensation to minimize globally the effect of the formation of energetically expensive boundaries separating martensite and austenite. The purely elastic drive towards the formation of microstructure is, of course, moderated by the gradient terms in the energy which limit the spatial scale of the emerging microstructure.

To show that the complexity of such microstructure depends crucially on the presence of inertial terms in the governing equations, we also performed numerical experiments under the assumption of an overdamped dynamics with $\tilde{\rho}=0$. In this case, we essentially deal with the classical Ginzburg–Landau model, which ends up generating considerably simpler equilibrium patterns even though microstructures do form as compatible coexistence of martensitic variants with the austenite matrix is still an issue [155, 163–169]. The overdamped model obviously misses the possibility of self-focusing of elastic waves and therefore does not allow for secondary nucleation events which in the underdamped case contribute in a crucial way to the development of multiple length scales. In other words, without bringing inertia into the Ginzburg–Landau model, one completely misses the effect of *dynamic triggering* of microstructural complexity.
<figure><img src="image_24.png" /><figcaption>Figure 21. Probability density function of avalanche energies in the 3D underdamped continuum model. Inset: final configuration with 3 martensitic variants represented by different colors.</figcaption></figure>

### 5.3. Statistics of avalanches

Consider next the issues of intermittency and scaling in the 3D model. In Figure 20, we show the time evolution of the normalized dissipated energy $\mathcal{R}$ during the whole cooling process beyond the initial nucleation stage. We use several levels of magnification in Figure 20(a)–(c) to show that the structure of the time series remains basically the same. It suggests that the corresponding dynamics is intermittent exhibiting broad distribution of scales.

The biggest burst occurs at the very early stage of the phase transformation when the transition is just initiated. As we have seen, at this stage, a complex texture of variously oriented variants of martensite appears almost instantly inside a significant part of the volume of the austenite phase, see Figure 20 (a). The system size avalanche, however, is shaped itself by a large number of small bursts representing pre- and after-shocks, see Figure 20(b) and (c). As a large number of localized transformation events occur at almost the same time, the total dissipated energy is computed by integration over the whole sample which makes some superposition of bursts inevitable.

The statistical structure of the observed fluctuations is characterized using the distribution of our three main observables: $S$, $E$, and $T$. More specifically, we interpreted the number of time steps during an avalanche as its duration $T$. The total dissipated energy $E$ is defined again as the discrete sum of incremental values of the energy dissipation over the duration of each avalanche. We also computed the maximum of the amplitude for each of the avalanches shown in Fig. 20 and denoted by $S$, because its distribution turned out to be statistically equivalent to the distribution of avalanche amplitudes, calculated in 1D using equations (19) and (20) within the duration $T$ of an avalanche. As in the 1D model, we used an irrelevant threshold to identify individual avalanches and ended up processing approximately 1800 avalanche events. Statistical distributions were constructed using logarithmic binning, with the lower bound $x_{min}$ determined by minimizing the Kolmogorov–Smirnov distance. Scaling exponents were calculated via fitting to data above $x_{min}$, see [170] for details.

By dropping particularly small- and large-size events, representing the effects of viscosity, gradient regularization, and post-processing as well as the system size avalanches at the beginning of each new loading cycle, we obtain the probability density distribution of the dissipated energy shown in Figure 21. Over four decades, it is clearly of power law type

$$P(E) \sim E^{-\beta} \qquad (64)$$

with exponent $\beta = 1.6$. The attendant statistical distributions of avalanche durations $T$ and avalanche amplitudes $S$ are shown in Figure 22(a) and (b). Here, we see a broader range of event sizes with again almost perfect power law-type ranges characterized by the distributions

$$P(T) \sim T^{-\tau}, \quad P(S) \sim S^{-\kappa}, \qquad (65)$$
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Statistics of avalanches in 3D underdamped continuum model: (a) Distribution of avalanche amplitudes. (b) Distribution of avalanche durations.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 23.** Correlations between avalanche energies *E* and avalanche amplitudes in 3D underdamped continuum model.</figcaption></figure>

<table><caption>**Table 1.** Comparison of the power law exponents measured in experiments [63] with those computed in the 3D underdamped continuum model.</caption><thead><tr><th>κ<sup>exp</sup></th><th>τ<sup>exp</sup></th><th>β<sup>exp</sup></th><th>z<sup>exp</sup></th></tr></thead><tbody><tr><td>2.26 ± 0.1</td><td>-</td><td>1.64 ± 0.1</td><td>1.97 ± 0.4</td></tr><tr><td>κ</td><td>τ</td><td>β</td><td>z</td></tr><tr><td>2.1 ± 0.1</td><td>2.4 ± 0.1</td><td>1.6 ± 0.05</td><td>1.8 ± 0.3</td></tr></tbody></table>

which exponents $τ = 2.4$ and $κ = 2.15$. Finally, in Figure 23, we show the distributions of avalanche energies *E* plotted against the distribution of avalanche amplitudes *S* which shows a presence of a correlation of the form

$$E \sim S^z \qquad (66)$$

with exponent $z = 1.8 \pm 0.3$; the fact, that $z \neq 2$ is a signature of complexity developing in this system at both temporal and spatial levels.

Given that we fixed the simulation parameters using experimental measurements available in the literature for the cubic to tetragonal martensitic transformation in Fe₆₈.₈Pd₃₁.₂, we can compare our
numerical results with experimental data obtained from the measurements of acoustic activity in such crystals subjected to cyclic thermal driving, see [63] for details. In experiments, avalanches amplitudes $S$ and avalanche durations $T$ were extracted directly from the AE signal while the avalanche energies $E$ were obtained by integrating the square of such signal. In our Table 1, the numerical values of the exponents found in experiments and simulations are juxtaposed.

One can see that the computed exponents characterizing avalanche amplitudes ($\kappa$) and avalanche energies ($\beta$) are in excellent agreement with experiment. The experimental value of the exponent characterizing energy-amplitude correlations ($z$) is within the error bars of our computed value. Since in experiments apparently no power law behavior was identified for avalanche durations, the picture is now completed by the results of our simulations, where small but finite interval of power law behavior was found.

To summarize, despite all the drastic simplifications made in the modeling of martensitic transformation in $Fe_{68.8}Pd_{31.2}$ crystals, we seem to have succeeded to capture surprisingly well the values of experimentally measured exponents. This can be seen as a reflection of the scale-free (critical) nature of the studied phenomenon, which makes the values of the exponents insensitive to the details of the modeling scheme [62,171]. Given that our model accounts for neither particular quenched disorder, nor dislocational activity or any other type of annealed disorder, and knowing that intermittency and scaling are not captured in the corresponding overdamped model, we may conclude that the emergence of an extended 'inertial range' in our numerical experiments is indeed intimately linked to the underdamped nature of dynamics.

## 6. Conclusion

The goal of this study was to identify the origin of the power law AE in martensites subjected to quasi-static loading. Despite numerous proposals regarding the possible mechanisms of observed scaling behavior, the definitive model accounting for self-organization towards 'generic' criticality in martensites has been missing. Here, by following some earlier insights, we explored the possibility that the extended scaling regime can be, at least partially, interpreted as an effect of a constructive interplay between inertia and elastic incompatibility.

Our interest in the effects of inertia is rooted in the idea that dynamics should play an important role during steady-state transformational plasticity because the underlying mechanical system is close to being only marginally stable. Indeed, it is known that the associated quasi-plastic yielding operates in the part of the energetic landscape where the barriers, separating generic metastable states, are extremely small. When inertia is incorporated into such a model, randomly generated elastic waves would not be severely obstructed and scattered by these barriers and can therefore interact sufficiently freely. In particular, they can self-focus creating virtual nucleation sites which can potentially facilitate growth of a new phase simultaneously at many scales. In a different perspective, the emerging wave activity can be perceived as representing annealed self-induced disorder that can self-tune internally allowing the system to reach the scaling regime. In all these processes, the background role of elastic incompatibility reduces to shaping the energy landscape towards coexistence of sufficiently rich metastability with almost marginality.

While the detailed mathematical structure of the implied feedback mechanisms still remains unclear, in this paper, we provided compelling evidence that in a model accounting for inertial effects and at the same time capable of exhibiting both metastability and marginality, the desired regime of intermittency and scaling can be reached. One can say that in such systems inertia serves as a mechanism of self-organization towards an extended dynamical critical state.

We started with a study of an elementary example of a 1D chain of mass points connected by springs. While the crucial assumption, allowing one to model in this way elastic phase transitions, was that the interactions between the NNs are characterized by a nonconvex potential, harmonic NNN interactions were also taken into account. Such an extension of the minimal NN model was supposed to mimic the ferromagnetic interactions, and therefore the developed model can be viewed as a soft-spring version of a 1D RFIM. The main novelty is in the focus on underdamped dynamics instead of a more conventional overdamped dynamics. We showed that in this setting one can reach the scale-free behavior regime without quenched disorder and without tuning. The incorporation of NNN interactions into a 1D model
is a ‘poor man’s’ attempt to account for strain incompatibility between austenite and martensite energy
wells. In the same spirit, the assumption of discreteness, responsible for the abundance of energy wells in
such model, can be viewed as a way to imitate the complexity of actual energy landscape in transforming
martensites.

We then explored the effects of underdamping in more realistic 2D and 3D continuum models. Here,
the desired balance between metastability and marginality can be attributed to elastic incompatibility
which generates long-range elastic interactions. While overdamped Ginzburg-Landau-type continuum
models do not exhibit either intermittency or scaling in the absence of quenched disorder, we showed
that the problem can be remedied by the account of inertia. Specifically, using a physically realistic 3D
continuum model, we were able to show that the energies, the amplitudes, and the durations of inter-
mittent transformation-induced avalanches all exhibit the expected power law behavior. The computed
exponents were found to be in good agreement with those found in experimental studies which suggests
that the model captures adequately the associated inertia-dominated universality class.

This, of course, does not exclude the relevance of the alternative mechanisms of self-organization
towards criticality in martensites including, for instance, the one relying on dislocational activity. In
this perspective, an interesting generalization of this study would be to take into account the possibil-
ity of inertia-induced dislocation nucleation. This would allow one to clarify the role of plastic activity
in thermoelastic martensitic transformations. As we have already mentioned, some of the martensitic
transformations exhibit scaling only after extensive cyclic loading, and the role of inertia in the associ-
ated self-organization processes still remains poorly understood. Inertia may be also playing an important
role in dynamic phase transitions, say inside shock waves, which would be another class of problems to be
addressed using the tools developed in this paper. Those are of course severely rate-dependent phenom-
ena, and the emergence in such problems of additional, loading dependent, time scales can be expected
to compromise scaling at least to some extent.

The most immediate extension of this work would be to investigate the role of inertia in martensitic
transitions differing by their crystallographic symmetry and the degree of incompatibility of the energy
wells. This can help to explain why some of these phase transitions do not exhibit scale-free AE, which
remains a puzzling question.

Finally, we mention that the partition of the energy of wave motion excited by the transformation into
elastic radiation, measured in AE experiments, and thermal heating, detected by infrared cameras, is still
not well understood due to nonlinearity involved in the energy transfer from long to short waves and in the
process of the eventual thermalization of lattice scale waves. Studying these phenomena inside a single
model can reveal additional relations between inertia and dissipation and contribute to the understanding
of how a particular degree of underdamping, ensuring scaling, and criticality, can be actually achieved
in real systems.

Authors' note

The paper is dedicated to Marcelo Epstein on the occasion of his 80th birthday.

Acknowledgements

The authors are grateful to Francisco Perez Reche for valuable comments and suggestions.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: O. U. S. was supported by the grants ANR-18-CE42-0017-03, ANR-19-CE08-0010-01, ANR-20-CE91-0010 and MRTN-CT-2004-505226. A. F. was supported by the grants MRTN-CT-2004-505226. L. T. was supported by the grants ANR-17-CE08-0047-02, ANR-21-CE08-MESOCRYSP and ERC-H2020-MSCA-RISE-2020-101008140.
References

[1] Planes, A, and Vives, E. Avalanche criticality in thermal-driven martensitic transitions: the asymmetry of the forward and reverse transitions in shape-memory materials. *J Phys Condens Matter* 2017; 29(33): 334001.

[2] Porta, M, Castán, T, Saxena, A, et al. Influence of the number of orientational domains on avalanche criticality in ferroelastic transitions. *Phys Rev E* 2019; 100(6): 062115.

[3] Planes, A, Mañosa, L, and Vives, E. Acoustic emission in martensitic transformations. *J Alloys Compd* 2013; 577: S699-S704.

[4] Perez-Reche, FJ. Modelling avalanches in martensites. In: Planes, A, Saxena, A, and Salje, EKH (eds) *Avalanches in functional materials and geophysics*. Cham: Springer, 2017, pp. 99-136.

[5] Balandraud, X, Barrera, N, Biscari, P, et al. Strain intermittency in shape-memory alloys. *Phys Rev B* 2015; 91(17): 174111.

[6] Sethna, JP, Bierbaum, MK, Dahmen, KA, et al. Deformation of crystals: connections with statistical physics. *Annu Rev Mater Res* 2017; 47(1): 217-246.

[7] Alava, MJ, Nukala, PK, and Zapperi, S. Statistical models of fracture. *Adv Phys* 2006; 55(3-4): 349-476.

[8] Alava, MJ, Laurson, L, and Zapperi, S. Crackling noise in plasticity. *Eur Phys J Spec Top* 2014; 223(11): 2353-2367.

[9] Miguel, MC, Vespignani, A, Zapperi, S, et al. Intermittent dislocation flow in viscoplastic deformation. *Nature* 2001; 410: 667-671.

[10] Zaiser, M. Scale invariance in plastic flow of crystalline solids. *Adv Phys* 2006; 55(1-2): 185-245.

[11] Borja da Rocha, H, and Truskinovsky, L. Mean field fracture in disordered solids: statistics of fluctuations. *J Mech Phys Solids* 2022; 158: 104646.

[12] Salman, OU, and Truskinovsky, L. Minimal integer automaton behind crystal plasticity. *Phys Rev Lett* 2011; 106(17): 175503.

[13] Zhang, P, Salman, OU, Weiss, J, et al. Variety of scaling behaviors in nanocrystalline plasticity. *Phys Rev E* 2020; 102(2-1): 023006.

[14] Salman, OU, Finel, A, Delville, R, et al. The role of phase compatibility in martensite. *J Appl Phys* 2012; 111(10): 103517.

[15] Fisher, DS. Sliding charge-density waves as a dynamic critical phenomenon. *Phys Rev B* 1985; 31: 1396.

[16] Rosso, A, Sethna, JP, and Wyart, M. Avalanches and deformation in glasses and disordered systems. In: *Emergent dynamics in glasses and disordered systems: correlations and avalanches*, Spin Glass Theory and Far Beyond: Replica Symmetry Breaking After 40 Years. Singapore: World Scientific Publishing, 2022, pp. 277-305.

[17] LeBlanc, M, Uhl, JT, and Dahmen, KA. Avalanches in strained bulk metallic glasses. *Phys Rev Lett* 2014; 112: 155501.

[18] Denisov, D, Lorincz, K, Uhl, J, et al. Universality of slip avalanches in flowing granular matter. *Nat Commun* 2016; 7: 10641.

[19] Liu, C, Ferrero, EE, Puosi, F, et al. Driving rate dependence of avalanche statistics and shapes at the yielding transition. *Phys Rev Lett* 2016; 116: 065501.

[20] Chan, PY, Tsekenis, G, Dantzig, J, et al. Plasticity and dislocation dynamics in a phase field crystal model. *Phys Rev Lett* 2010; 105: 015502.

[21] Friedman, N, Jennings, AT, Tsekenis, G, et al. Statistics of dislocation slip avalanches in nanosized single crystals show tuned critical behavior. *Phys Rev Lett* 2012; 109: 095507.

[22] Bak, P, Tang, C, and Wiesenfeld, K. Self-organized criticality: an explanation of the 1/f noise. *Phys Rev Lett* 1987; 59(4): 381-384.

[23] Dickman, R, Muñoz, MA, Vespignani, A, et al. Paths to self-organized criticality. *Braz J Phys* 2000; 30: 27-41.

[24] Jensen, HJ. *Self-organized criticality: emergent complex behavior in physical and biological systems*, vol. 10. Cambridge: Cambridge University Press, 1998.

[25] Dhar, D. Theoretical studies of self-organized criticality. *Physica A* 2006; 369(1): 29-70.

[26] Pruessner, G. *Self-organised criticality: theory, models and characterisation*. Cambridge: Cambridge University Press, 2012.

[27] Sornette, D. *Critical phenomena in natural sciences: chaos, fractals, self-organization and disorder: concepts and tools*. Berlin and Heidelberg: Springer, 2006.

[28] Gros, C. Self-organized criticality. In: Gros, C (ed.) *Complex and adaptive dynamical systems: a comprehensive introduction*. Cham: Springer, 2024, pp. 203-239.

[29] Tadić, B, and Melnik, R. Self-organised critical dynamics as a key to fundamental features of complexity in physical, biological, and social networks. *Dynamics* 2021; 1(2): 181-197.

[30] Müller, M, and Wyart, M. Marginal stability in structural, Spin, and electron glasses. *Annu Rev Condens Matter Phys* 2015; 6: 177.

[31] Shang, B, Guan, P, and Barrat, JL. Elastic avalanches reveal marginal behavior in amorphous solids. *Proc Natl Acad Sci* 2020; 117(1): 86-92.

[32] Ovaska, M, Lehtinen, A, Alava, MJ, et al. Excitation spectra in crystal plasticity. *Phys Rev Lett* 2017; 119(26): 265501.

[33] Franz, S, and Spigler, S. Mean-field avalanches in jammed spheres. *Phys Rev E* 2017; 95(2): 022139.
[34] Nishiyama, Z. *Martensitic transformation*. Amsterdam: Elsevier, 2012.

[35] Roitburd, A. Martensitic transformation as a typical phase transformation in solids. *Solid State phys* 1978; 33: 317-390.

[36] James, RD, and Hane, KF. Martensitic transformations and shape-memory materials. *Acta Mater* 2000; 48(1): 197-222.

[37] Khachaturyan, AG. *Theory of structural transformations in solids*. North Chelmsford, MA: Courier Corporation, 2013.

[38] Bhattacharya, K. *Microstructure of martensite*. Oxford: Oxford University Press, 2003.

[39] Pérez-Reche, FJ, Stipcich, M, Vives, E, et al. Kinetics of martensitic transitions in Cu-Al-Mn under thermal cycling: analysis at multiple length scales. *Phys Rev B* 2004; 69: 064101.

[40] Chandni, U, Ghosh, A, Vijaya, HS, et al. Criticality of tuning in athermal phase transitions. *Phys Rev Lett* 2009; 102(2): 025701.

[41] Rosinberg, ML, and Vives, E. Metastability, hysteresis, avalanches, and acoustic emission: martensitic transitions in functional materials. In: Planes, A, Saxena, A, Fukuda, T, et al. (eds) *Disorder and strain-induced complexity in functional materials*. Berlin and Heidelberg: Springer, 2011, pp. 249-272.

[42] Vives, E, Soto Parra, DE, Planes, A, et al. Acoustic emission avalanches in martensitic transitions: new perspectives for the problem of source location. *Solid State Phenom* 2011; 172: 144-149.

[43] Gallardo, MC, Manchado, J, Romero, FJ, et al. Avalanche criticality in the martensitic transition of Cu 67.64 Zn 16.71 Al 15.65 shape-memory alloy: a calorimetric and acoustic emission study. *Phys Rev B* 2010; 81(17): 174102.

[44] Beke, DL, Daróczi, L, Tóth, LZ, et al. Acoustic emissions during structural changes in shape memory alloys. *Metals* 2019; 9(1): 58.

[45] Illa, X, Winkelmayer, P, and Vives, E. Local strain variability and force fluctuations during the martensitic transition under different driving mechanisms. *Phys Rev B* 2015; 92(18): 184107.

[46] Vives, E, Soto-Parra, D, Mañosa, L, et al. Driving-induced crossover in the avalanche criticality of martensitic transitions. *Phys Rev B* 2009; 80: 180101.

[47] Torrents, G, Illa, X, Vives, E, et al. Geometrical model for martensitic phase transitions: understanding criticality and weak universality during microstructure growth. *Phys Rev E* 2017; 95(1-1): 013001.

[48] Beke, DL, Bolgár, MK, Tóth, LZ, et al. On the asymmetry of the forward and reverse martensitic transformations in shape memory alloys. *J Alloys Compd* 2018; 741: 106-115.

[49] Song, Y, Chen, X, Dabade, V, et al. Enhanced reversibility and unusual microstructure of a phase-transforming material. *Nature* 2013; 502(7469): 85-88.

[50] Sethna, JP, Dahmen, KA, and Myers, CR. Crackling noise. *Nature* 2001; 410: 242.

[51] Sethna, JP, Dahmen, K, Kartha, S, et al. Hysteresis and hierarchies: dynamics of disorder-driven first-order phase transformations. *Phys Rev Lett* 1993; 70(21): 3347-3350.

[52] Perković, O, Dahmen, K, and Sethna, JP. Avalanches, Barkhausen noise, and plain old criticality. *Phys Rev Lett* 1995; 75(24): 4528-4531.

[53] Dahmen, K, and Sethna, JP. Hysteresis, avalanches, and disorder-induced critical scaling: a renormalization-group approach. *Phys Rev B* 1996; 53(22): 14872-14905.

[54] Vives, E, and Planes, A. Avalanches in a fluctuationless first-order phase transition in a random-bond ising model. *Phys Rev B* 1994; 50(6): 3839-3848.

[55] Tadic, B. Nonuniversal scaling behavior of Barkhausen noise. *Phys Rev Lett* 1996; 77(18): 3843-3846.

[56] Vives, E, and Planes, A. Hysteresis and avalanches in the random anisotropy ising model. *Phys Rev B* 2001; 63(13): 134431.

[57] Pérez-Reche, FJ, Truskinovsky, L, and Zanzotto, G. Driving-induced crossover: from classical criticality to self-organized criticality. *Phys Rev Lett* 2008; 101(23): 27.

[58] Sethna, JP, Dahmen, KA, and Perkovic, O. *Random-field ising models of hysteresis*. New York: Academic Press, 2006.

[59] Pérez-Reche, FJ, Truskinovsky, L, and Zanzotto, G. Training-induced criticality in martensites. *Phys Rev Lett* 2007; 99(7): 075501.

[60] Perez-Reche, FJ, Triguero, C, Zanzotto, G, et al. Origin of scale-free intermittency in structural first-order phase transitions. *Phys Rev B* 2016; 94(14): 144102.

[61] Perez-Reche, FJ, Truskinovsky, L, and Zanzotto, G. Martensitic transformations: from continuum mechanics to spin models and automata. *Contin Mech Thermodyn* 2009; 21: 17-26.

[62] Carrillo, L, Mañosa, L, Ortín, J, et al. Experimental evidence for universality of acoustic emission avalanche distributions during structural transitions. *Phys Rev Lett* 1998; 81(9): 1889-1892.

[63] Bonnot, E, Romero, R, Mañosa, L, et al. Elastocaloric effect associated with the martensitic transition in shape-memory alloys. *Phys Rev Lett* 2008; 100(12): 125901.

[64] Eshelby, JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. *Proc R Soc Lon A* 1957; 241: 376.

[65] Picard, G, Ajdari, A, Lequeux, F, et al. Elastic consequences of a single plastic event: a step towards the microscopic modeling of the flow of yield stress fluids. *Eur J E* 2004; 15: 371.

[66] Rossi, S, Biroli, G, Ozawa, M, et al. Far-from-equilibrium criticality in the random-field ising model with Eshelby interactions. *Phys Rev B* 2023; 108(22): L220202.
[67] Clapp, P. How would we recognize a martensitic transformation if it bumped into us on a dark & austy night? *J Phys IV* 1995; 5(C8): C8-11.

[68] Miao, Y, and Vlassak, JJ. Explosive martensitic transformation of supercooled austenite in CuZr-based thin-film shape memory alloys. *Acta Mater* 2020; 200: 162–170.

[69] Schwabe, S, Lünser, K, Schmidt, D, et al. What is the speed limit of martensitic transformations? *Sci Technol Adv Mater* 2022; 23(1): 633–641.

[70] Lin, J, and Pence, TJ. Pulse attenuation by kinetically active phase boundary scattering during displacive phase transformations. *J Mech Phys Solids* 1998; 46(7): 1183–1211.

[71] Truskinovsky, L. Nucleation and growth in elastodynamics. In: Duxbury, PM, and Pence, TJ (eds) *Dynamics of crystal surfaces and interfaces*. Boston, MA: Springer, 2002, pp. 185–197.

[72] Lookman, T, Shenoy, SR, Rasmussen, KO, et al. Ferroelastic dynamics and strain compatibility. *Phys Rev B* 2003; 67(2): 024114.

[73] Truskinovsky, L, and Vainchtein, A. Dynamics of martensitic phase boundaries: discreteness, dissipation and inertia. *Contin Mech Thermodyn* 2008; 20(2): 97–122.

[74] Steinbach, I, and Shchyglo, O. Phase-field modelling of microstructure evolution in solids: perspectives and challenges. *Curr Opin Solid State Mater Sci* 2011; 15(3): 87–92.

[75] Cho, JY, Idesman, A, Levitas, V, et al. Finite element simulations of dynamics of multivariant martensitic phase transitions based on Ginzburg–Landau theory. *Int J Solids Struct* 2012; 49(14): 1973–1992.

[76] Liu, X, Schneider, D, Reder, M, et al. Modeling of martensitic phase transformation accounting for inertia effects. *Int J Mech Sci* 2024; 278: 109443.

[77] Vives, E, Ortín, J, Mañosa, L, et al. Distributions of avalanches in martensitic transformations. *Phys Rev Lett* 1994; 72(11): 1694–1697.

[78] Mañosa, L, Planes, A, Rouby, D, et al. Acoustic emission field during thermoelastic martensitic transformations. *Appl Phys Lett* 1989; 54(25): 2574–2576.

[79] Baró, J, Corral, Á, Illa, X, et al. Statistical similarity between the compression of a porous material and earthquakes. *Phys Rev Lett* 2013; 110(8): 088702.

[80] Gomberg, J, and Johnson, P. Dynamic triggering of earthquakes. *Nature* 2005; 437(7060): 830–830.

[81] Johnson, PA, and Jia, X. Nonlinear dynamics, granular media and dynamic earthquake triggering. *Nature* 2005; 437(7060): 871–874.

[82] Kocharyan, G. Nucleation and evolution of sliding in continental fault zones under the action of natural and man-made factors: a state-of-the-art review. *Izv Phys Solid Earth* 2021; 57: 439–473.

[83] Khfifi, M, and Loulidi, M. Scaling properties of a rice-pile model: inertia and friction effects. *Phys Rev E* 2008; 78: 051117.

[84] Maimon, R, and Schwarz, JM. Continuous depinning transition with an unusual hysteresis effect. *Phys Rev Lett* 2004; 92: 255502.

[85] Papanikolaou, S. Shearing a glass and the role of pinning delay in models of interface depinning. *Phys Rev E* 2016; 93: 032610.

[86] Carlson, JM, and Langer, JS. Properties of earthquakes generated by fault dynamics. *Phys Rev Lett* 1989; 62: 2632.

[87] Prado, CPC, and Olami, Z. Inertia and break of self-organized criticality in sandpile cellular-automata models. *Phys Rev A* 1992; 45: 665.

[88] Held, GA, Solina, DH, Keane, DT, et al. Experimental study of critical-mass fluctuations in an evolving sandpile. *Phys Rev Lett* 1990; 65: 1120.

[89] Jaeger, HM, Liu, Ch, and Nagel, SR. Relaxation at the angle of repose. *Phys Rev Lett* 1989; 62: 40.

[90] Marchetti, MC. Depinning and plasticity of driven disordered lattices. In: Carmen Miguel, M, and Rubi, M (eds) *Jamming, yielding, and irreversible deformation in condensed matter*, Lecture Notes in Physics. Berlin and Heidelberg: Springer-Verlag, 2006, pp. 137–157.

[91] Denisov, DV, Lorincz, KA, Wright, WJ, et al. Universal slip dynamics in metallic glasses and granular matter – linking frictional weakening with inertial effects. *Scientific Reports* 2017; 7: 43376.

[92] Nicolas, A, Barrat, JL, and Rottler, J. Effects of inertia on the steady-shear rheology of disordered solids. *Phys Rev Lett* 2016; 116: 058303.

[93] Karimi, K, and Barrat, JL. Role of inertia in the rheology of amorphous systems: a finite-element-based elastoplastic model. *Phys Rev E* 2016; 93(2): 022904.

[94] de Geus, TW, and Wyart, M. Short-range depinning in the presence of velocity-weakening. arXiv [preprint], 2024. DOI: 10.48550/arXiv.2411.06732.

[95] Clancy, I, and Corcoran, D. Criticality in the Burridge-Knopoff model. *Phys Rev E* 2005; 71(4): 046124.

[96] Salerno, KM, Maloney, CE, and Robbins, MO. Avalanches in strained amorphous solids: does inertia destroy critical behavior? *Phys Rev Lett* 2012; 109: 105703.

[97] Salerno, KM, and Robbins, MO. Effect of inertia on sheared disordered solids: critical scaling of avalanches in two and three dimensions. *Phys Rev E* 2013; 88: 062206.
[98] Karimi, K, Ferrero, EE, and Barrat, JL. Inertia and universality of avalanche statistics: the case of slowly deformed amorphous solids. *Phys Rev E* 2017; 95: 013003.

[99] Truskinovsky, L. Transition to detonation in dynamic phase changes. *Arch Ration Mech Anal* 1994; 125: 375-397.

[100] Reid, A, and Gooding, R. Pattern formation in a 2D elastic solid. *Physica A* 1997; 239(1-3): 1-10.

[101] Bales, G, and Gooding, R. Interfacial dynamics at a first-order phase transition involving strain: dynamical twin formation. *Phys Rev Lett* 1991; 67(24): 3412-3415.

[102] Elmer, FJ. Avalanches in the weakly driven Frenkel-Kontorova model. *Phys Rev E* 1994; 50(6): 4470.

[103] Ding, X, Lookman, T, Zhao, Z, et al. Dynamically strained ferroelastics: statistical behavior in elastic and plastic regimes. *Phys Rev B* 2013; 87(9): 094109.

[104] Ahluwalia, R, and Ananthakrishna, G. Power-law statistics for avalanches in a martensitic transformation. *Phys Rev Lett* 2001; 86(18): 4076-4079.

[105] Sreekala, S, Ahluwalia, R, and Ananthakrishna, G. Precursors and power-law statistics of acoustic emission and shape memory effect in martensites. *Phys Rev B* 2004; 70: 224105.

[106] Paul, A, Bhattacharya, J, Sengupta, S, et al. Non-affine deformation in microstructure selection in solids II: elastoplastic theory for thedynamics of solid state transformations. *J Phys Condens Matter* 2008; 20(36): 365211.

[107] Ball, JM, and James, R. Incompatible sets of gradients and metastability. *Arch Ration Mech Anal* 2015; 218: 1363-1416.

[108] Bouville, M, and Ahluwalia, R. Effect of lattice-mismatch-induced strains on coupled diffusive and displacive phase transformations. *Phys Rev B* 2007; 75(5): 054110.

[109] Zhang, Z, James, RD, and Müller, S. Energy barriers and hysteresis in martensitic phase transformations. *Acta Mater* 2009; 57(15): 4332-4352.

[110] Atli, K, Franco, B, Karaman, I, et al. Influence of crystallographic compatibility on residual strain of TiNi based shape memory alloys during thermo-mechanical cycling. *Mater Sci Eng A* 2013; 574: 9-16.

[111] Yang, Y, Xu, Y, Zhou, Y, et al. Nonthermoelastic martensitic features in ideal martensites due to volume effects. *Phys Rev B* 2023; 108(2): 024102.

[112] Shenoy, SR, Lookman, T, Saxena, A, et al. Martensitic textures: multiscale consequences of elastic compatibility. *Phys Rev B* 1999; 60(18): R12537-R12541.

[113] Ahluwalia, R, Lookman, T, and Saxena, A. Dynamic strain loading of cubic to tetragonal martensites. *Acta Mater* 2006; 54(8): 2109-2120.

[114] Kertész, J, and Kiss, L. The noise spectrum in the model of self-organised criticality. *J Phys Math General* 1990; 23(9): L433.

[115] Kuntz, MC, and Sethna, JP. Noise in disordered systems: the power spectrum and dynamic exponents in avalanche models. *Phys Rev B* 2000; 62: 11699-11708.

[116] Laurson, L, Alava, MJ, and Zapperi, S. Power spectra of self-organized critical sandpiles. *J Stat Mech* 2005; 2005(11): L11001.

[117] Laurson, L, and Alava, MJ. 1/f noise and avalanche scaling in plastic deformation. *Phys Rev E* 2006; 74(6): 066106.

[118] Nandi, MK, Sarracino, A, Herrmann, HJ, et al. Scaling of avalanche shape and activity power spectrum in neuronal networks. *Phys Rev E* 2022; 106(2): 024304.

[119] Ruseckas, J, and Kaulakys, B. Scaling properties of signals as origin of 1/f noise. *J Stat Mech* 2014; 2014(6): P06005.

[120] Kay, SM. *Fundamentals of statistical signal processing: estimation theory*. Hoboken, NJ: Prentice-Hall, 1993.

[121] Truskinovsky, L, and Vainchtein, A. The origin of nucleation peak in transformational plasticity. *J Mech Phys Solids* 2004; 52(6): 1421-1446.

[122] Puglisi, G, and Truskinovsky, L. Thermodynamics of rate-independent plasticity. *J Mech Phys Solids* 2005; 53(3): 655-679.

[123] Mielke, A, and Truskinovsky, L. From discrete visco-elasticity to continuum rate-independent plasticity: rigorous results. *Arch Ration Mech Anal* 2011; 203(2): 577-619.

[124] Salman, O, and Truskinovsky, L. On the critical nature of plastic flow: one and two dimensional models. *Int J Eng Sci* 2012; 59: 219-254.

[125] Truskinovsky, L, and Vainchtein, A. Kinetics of martensitic phase transitions: lattice model. *SIAM J Appl Math* 2005; 66(2): 533-553.

[126] Efendiev, YR, and Truskinovsky, L. Thermalization of a driven bi-stable FPU chain. *Continuum Mech Thermodyn* 2010; 22(6-8): 679-698.

[127] West, BJ, and Shlesinger, MF. On the ubiquity of 1/f noise. *Int J Mod Phys B* 1989; 3(6): 795-819.

[128] Milotti, E. 1/f noise: a pedagogical review. 2002. DOI: 10.48550/arXiv.physics/0204033.

[129] Polizzi, S, Pérez-Reche, FJ, Arneodo, A, et al. Power-law and log-normal avalanche size statistics in random growth processes. *Phys Rev E* 2021; 104(5): L052101.

[130] Sornette, D. Sweeping of an instability: an alternative to self-organized criticality to get powerlaws without parameter tuning. *J Phys* 1994; 4(2): 209-221.

[131] Sornette, D. Power laws without parameter tuning: an alternative to self-organized criticality. *Phys Rev Lett* 1994; 72(14): 2306.
[132] Vivès, E, Ortín, J, Mañosa, L, et al. Distribution of acoustic emission avalanches in martensitic transformations. *J Phys IV* 1995; 5(C2): C2-59.

[133] Urbach, JS, Madison, RC, and Markert, JT. Interface depinning, self-organized criticality, and the Barkhausen effect. *Phys Rev Lett* 1995; 75(2): 276-279.

[134] Durin, G, Bertotti, G, and Magni, A. Fractals, scaling and the question of self-organized criticality in magnetization processes. *Fractals* 1995; 3(2): 351-370.

[135] Schmelzer, J, and Ulbricht, H. Thermodynamics of finite systems and the kinetics of first-order phase transitions. *J Colloid Interface Sci* 1987; 117(2): 325-338.

[136] Salman, OU, Muite, B, and Finel, A. Origin of stabilization of macrotwin boundaries in martensites. *Eur Phys J B* 2019; 92: 1-9.

[137] Carlson, JM, Langer, JS, Shaw, BE, et al. Intrinsic properties of a Burridge-Knopoff model of an earthquake fault. *Phys Rev A* 1991; 44(2): 884-897.

[138] Schmittbuhl, J, Vilotte, JP, and Roux, S. Velocity weakening friction: a renormalization approach. *J Geophys Res* 1996; 101(B6): 13911-13917.

[139] Elmer, FJ. Is self-organized criticality possible in dry friction? In: Persson, BNJ, and Tosatti, E (eds) *Physics of sliding friction*. Dordrecht: Springer, 1996, pp. 433-447.

[140] Gorbushin, N, Vainchtein, A, and Truskinovsky, L. Transition fronts and their universality classes. *Phys Rev E* 2022; 106(2): 024210.

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

[142] Puglisi, G, and Truskinovsky, L. A mechanism of transformational plasticity. *Contin Mech Thermodynamics* 2002; 14: 437-457.

[143] Puglisi, G, and Truskinovsky, L. Rate independent hysteresis in a bi-stable chain. *J Mech Phys Solids* 2002; 50(2): 165-187.

[144] Bhattacharya, K, Firoozye, NB, James, RD, et al. Restrictions on microstructure. *Proc R Soc Edinb A* 1994; 124(5): 843-878.

[145] Smyshlyaev, V, and Willis, J. On the relation of a three-well energy. *Proc R Soc Lond A* 1999; 455(1983): 779-814.

[146] Grabovsky, Y, and Truskinovsky, L. A class of nonlinear elasticity problems with no local but many global minimizers. *J Elast* 2023; 154(1): 147-171.

[147] Grabovsky, Y, and Truskinovsky, L. Rigidity-induced critical points. *Phys Rev E* 2024; 110(6): 064114.

[148] Curnoe, SH, and Jacobs, AE. Twin wall of proper cubic-tetragonal ferroelastics. *Phys Rev B Condens Matter* Mater Phys 2000; 62(18): R11925-R11928.

[149] Budiansky, B, and Truskinovsky, L. On the mechanics of stress-induced phase transformation in zirconia. *J Mech Phys Solids* 1993; 41(9): 1445-1459.

[150] Fadda, G, Truskinovsky, L, and Zanzotto, G. Unified Landau description of the tetragonal, orthorhombic, and monoclinic phases of zirconia. *Phys Rev B* 2002; 66(17): 174107.

[151] Salman, OU. *Modeling of spatio-temporal dynamics and patterning mechanisms of martensites by phase-field and lagrangian methods*. PhD Thesis, Université Pierre et Marie Curie, Paris, 2009.

[152] Ball, J, Holmes, PJ, James, R, et al. On the dynamics of fine structure. *J Nonlinear Sci* 1991; 1: 17-70.

[153] Friesecke, G, and McLeod, J. Dynamic stability of non-minimizing phase mixtures. *Proc R Soc Lond A* 1997; 453(1966): 2427-2436.

[154] Curnoe, SH, and Jacobs, AE. Time evolution of tetragonal-orthorhombic ferroelastics. *Phys Rev B* 2001; 64(6): 064101.

[155] Jacobs, AE, Curnoe, SH, and Desai, RC. Simulations of cubic-tetragonal ferroelastics. *Phys Rev B* 2003; 68(22): 224104.

[156] Muto, S, Oshima, R, and Fujita, FE. Elastic softening and elastic strain energy consideration in the fcc—fct transformation of Fe-Pd alloys. *Acta Metall Mater* 1990; 38(4): 685-694.

[157] Kartha, S, Krumhansl, JA, Sethna, JP, et al. Disorder-driven pretransitional tweed pattern in martensitic transformations. *Phys Rev B* 1995; 52(2): 803-822.

[158] Sato, M, Grier, BH, Shapiro, SM, et al. Effect of magnetic ordering on the lattice dynamics of FCC Fe1-xPdx. *J Phys* 1982; 12(10): 2117-2129.

[159] Oshima, R, Sugiyama, M, and Fujita, FE. Tweed structures associated with Fcc-Fct transformations in Fe-Pd alloys. *Metall Trans* 1988; 19(4): 803-810.

[160] Testardi, LR, and Bateman, TB. Lattice instability of high-transition-temperature superconductors. II. Single-Crystal V3Si results. *Phys Rev* 1967; 154(2): 402-410.

[161] Rodney, D, Bouar, YL, and Finel, A. Phase field methods and dislocations. *Acta Mater* 2001; 51: 17-30.

[162] Pitteri, M, and Zanzotto, G. *Continuum models for phase transitions and twinning in crystals*. Boca Raton, FL: Chapman and Hall/CRC, 2002.

[163] Artemev, A, Jin, Y, and Khachaturyan, AG. Three-dimensional phase field model and simulation of cubic → tetragonal martensitic transformation in polycrystals. *Philos Mag A* 2002; 82(6): 1249-1270.

[164] Mamivand, M, Zaeem, MA, and El Kadiri, H. A review on phase field modeling of martensitic phase transformation. *Comput Mater Sci* 2013; 77: 304-311.
[165] Ni, Y, Jin, YM, and Khachaturyan, AG. The transformation sequences in the cubic→tetragonal decomposition. *Acta Mater* 2007; 55(14): 4903–4914.

[166] Turteltaub, S, and Suiker, ASJ. A multiscale thermomechanical model for cubic to tetragonal martensitic phase transformations. *Int J Solids Struct* 2006; 43(14–15): 4509–4545.

[167] Ichitsubo, T, Tanaka, K, Koiwa, M, et al. Kinetics of cubic to tetragonal transformation under external field by the time-dependent Ginzburg-Landau approach. *Phys Rev B* 2000; 62: 5435–5441.

[168] Hildebrand, FE, and Miehe, C. A phase field model for the formation and evolution of martensitic laminate microstructure at finite strains. *Philos Mag* 2012; 92(34): 4250–4290.

[169] Idesman, A, Cho, J, and Levitas, V. Finite element modeling of dynamics of martensitic phase transitions. *Appl Phys Lett* 2008; 93: 043102.

[170] Clauset, A, Shalizi, CR, and Newman, MEJ. Power-law distributions in empirical data. *SIAM Rev Soc Ind Appl Math* 2009; 51(4): 661–703.

[171] Vives, E, Goicoechea, J, Ortín, J, et al. Universality in models for disorder-induced phase transitions. *Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics* 1995; 52(1): R5–R8.

[172] LeVeque, RJ. *Finite difference methods for ordinary and partial differential equations, steady state and time dependent problems*. Philadelphia, PA: SIAM, 2007.

# Appendix I

## Spatial and temporal discretisation schemes

Here, we provide some technical details behind our numerical implementation of discrete and continuum models.

*One-dimensional discrete problem.* To solve equation (12) with periodic boundary conditions, we used discrete Fourier transform. For a system of $N$ nodes, we can write

$$ \hat{u}(q_n) = \sum_{j=0}^{N-1} \tilde{u}_j e^{-iqnj\delta}, \quad q_n = \frac{2\pi n}{L}, $$

where $i = \sqrt{-1}$. Combining all terms, we arrive at the Fourier space representation:

$$ C\delta \frac{d^2\hat{u}(q)}{dt^2} - \frac{d\hat{u}(q)}{dt} = -2i \sin(q\delta) \frac{\hat{\phi}'}{\delta} + \frac{4G}{\delta^3} (1 - \cos(q\delta))^2 \hat{u}(q) $$

The nonlinear function $\tilde{\phi}'$ is first evaluated at each grid point using the current values of $\tilde{u}_i$. Once it is computed across the entire spatial domain, we apply the Fourier transform to obtain its spectral representation.

If we define a mode-dependent operator $L(q) = \frac{4E}{\delta^3}(1 - \cos(q\delta))^2$, we can rewrite our equation in the form

$$ C\delta \frac{d^2\hat{u}(q)}{dt^2} - \frac{d\hat{u}(q)}{dt} = -2i \sin(q\delta) \frac{\hat{\phi}'(q)}{\delta} + L(q)\hat{u}(q) $$

To integrate this equation in time, we employ a standard fourth-order Runge-Kutta RK4 algorithm. We first define new variables

$$ \hat{u}(q) = y_1 \qquad (67) $$

$$ \frac{d\hat{u}(q)}{dt} = y_2 \qquad (68) $$
and rewrite our second-order equation as a first-order system

$$ \frac{dy_1}{dt} = y_2 \qquad (69) $$

$$ \frac{dy_2}{dt} = \frac{1}{C\delta} \left[ y_2 - 2i \sin(q\delta) \frac{\hat{\phi}'(q)}{\delta} + L(q)y_1 \right] \qquad (70) $$

For each time step from $t_n$ to $t_{n+1} = t_n + \Delta t$, the RK4 algorithm proceeds as follows:

1. Compute the nonlinear term $\tilde{\phi}'$ in real space using the current $\tilde{u}_i$ values, then Fourier transform it to obtain $\hat{\phi}'(q)$

2. Execute the four RK4 steps:

$$ k_{1,1} = y_{2,n} \qquad (71) $$

$$ k_{1,2} = \frac{1}{C\delta} \left[ y_{2,n} - 2i \sin(q\delta) \frac{\hat{\phi}'(q)}{\delta} + L(q)y_{1,n} \right] \qquad (72) $$

$$ k_{2,1} = y_{2,n} + \frac{\Delta t}{2} k_{1,2} \qquad (73) $$

To compute $k_{2,2}$ and perform subsequent steps, we need to update the nonlinear term as follows:

- Transform $y_1 + \frac{\Delta t}{2}k_{1,1}$ back to real space

- Compute the nonlinear term $\tilde{\phi}'$ with these updated values

- Fourier transform it to obtain the update of $\hat{\phi}'(q)$

$$ k_{2,2} = \frac{1}{C\delta} \left[ \left(y_{2,n} + \frac{\Delta t}{2} k_{1,2} \right) - 2i \sin(q\delta) \frac{\hat{\phi}'(2)(q)}{\delta} + L(q) \left(y_{1,n} + \frac{\Delta t}{2} k_{1,1} \right) \right] \qquad (74) $$

$$ k_{3,1} = y_{2,n} + \frac{\Delta t}{2} k_{2,2} \qquad (75) $$

$$ k_{3,2} = \frac{1}{C\delta} \left[ \left(y_{2,n} + \frac{\Delta t}{2} k_{2,2} \right) - 2i \sin(q\delta) \frac{\hat{\phi}'(3)(q)}{\delta} + L(q) \left(y_{1,n} + \frac{\Delta t}{2} k_{2,1} \right) \right] \qquad (76) $$

$$ k_{4,1} = y_{2,n} + \Delta t k_{3,2} \qquad (77) $$

$$ k_{4,2} = \frac{1}{C\delta} \left[ \left(y_{2,n} + \Delta t k_{3,2} \right) - 2i \sin(q\delta) \frac{\hat{\phi}'(4)(q)}{\delta} + L(q) \left(y_{1,n} + \Delta t k_{3,1} \right) \right]. \qquad (78) $$

Here, $\hat{\phi}'(i)(q)$ denotes the Fourier transform of the nonlinear force $\phi'$ evaluated at the $i$th intermediate estimate of the displacement. Specifically, for each stage $i=2, 3, 4$, we compute:

* Transform $y_1^{(i)}$ back to real space via inverse FFT,

* Evaluate $\tilde{\phi}'$ pointwise in real space using the updated displacement,

* Then apply FFT to obtain $\hat{\phi}'(i)(q)$ in spectral form.

This ensures that the nonlinearity is updated consistently at each intermediate Runge-Kutta stage.

3. Update the solution:

$$ y_{1,n+1} = y_{1,n} + \frac{\Delta t}{6} (k_{1,1} + 2k_{2,1} + 2k_{3,1} + k_{4,1}) \qquad (79) $$

$$ y_{2,n+1} = y_{2,n} + \frac{\Delta t}{6} (k_{1,2} + 2k_{2,2} + 2k_{3,2} + k_{4,2}) \qquad (80) $$
4. Recover the Fourier coefficients:

$$ \hat{u}_{n+1}(q) = y_{1,n+1} \qquad (81) $$

$$ \left. \frac{d\hat{u}(q)}{dt} \right|_{n+1} = y_{2,n+1} \qquad (82) $$

This computational scheme provides fourth-order temporal accuracy and does not introduce numerical dissipation, preserving the spectral properties of our original equation. To implement non-periodic (fixed) boundary conditions in this setup, we used the ghost points technique which is described by Salman [151].

**One-dimensional continuum model.** To overcome the stiffness of the time-dependent problem, we had to use an implicit-explicit time-marching scheme with high temporal accuracy. The more straightforward explicit methods require less complex structure but are not suitable in our case because they require small time step $dt$ which does not allow one to reach steady state in a reasonable computational time. Below we illustrate the method using the simplest one-dimensional framework. Exactly the same type of equation is also solved in our 2D and 3D numerical experiments, and the corresponding generalization is straightforward.

In the interest of analytical transparency, we present the analysis for the following simplified quasi-linear differential equation which has the minimal required property of being second-order in time and fourth-order in space

$$ C\ddot{u} = \partial_x\phi'(\partial_x u) + \partial_{xx}\dot{u} - G\partial_{xxxx}u. \qquad (83) $$

Here, it is implied that the continuum system is discretized at $N$ nodes and that we deal again with periodic boundary conditions. The nonlinear term $f(\partial_x u)$ is taken to be the same as in our model 1D equation. At least second-order accuracy in time is required [172], and therefore, we choose second-order approximations to discretize time

$$ \ddot{u} \approx \frac{2u^{t+1} - 5u^t + 4u^{t-1} - u^{t-2}}{dt^2}, \qquad (84) $$

$$ \dot{u} \approx \frac{-3u^t + 6u^{t-1} + u^{t-2}}{6dt}. \qquad (85) $$

The remaining terms are approximated around the time step $t+1$. Thus, the Taylor expansion of the nonlinear terms $f^t$ and $f^{t-1}$ gives

$$ f^t \approx f^{t+1} - dtf^t_{t+1} + \frac{dt^2 f^t_{tt+1}}{2}, \qquad (86) $$

$$ f^{t-1} \approx f^{t+1} - 2dtf^t_{t+1} + 2dt^2f^t_{tt+1}. \qquad (87) $$

We can also write

$$ f^{t+1} \approx 2f^t - f^{t-1} + O(dt^2). \qquad (88) $$

The remaining linear terms are computed implicitly at time $(t+1)$. This allows us to express $u_i^{t+1}$ in terms of $u_i^t$, $u_i^{t-1}$ and $u_i^{t-2}$

$$ \left( \frac{2C}{dt^2} - \frac{2\partial_{xx}}{6dt} + G\partial_{xxxx} \right) u^{t+1} = \partial_x[2f_i^t - f^{t-1}] + \rho \frac{5u^t - 4u^{t-1} + u^{t-2}}{dt^2} - \partial_{xx} \frac{-3u_i^t + 6u^{t-1} + u^{t-2}}{6dt}. \qquad (89) $$

This equation can be rewritten compactly if we introduce the linear operator

$$ H = \left( \frac{2C}{dt^2} - \frac{2\partial_{xx}}{6dt} + G\partial_{xxxx} \right) \qquad (90) $$
and the non linear function $h_i^t$

$$h_i^t = \partial_x[2f^t - f^{t-1}] + C \frac{5u_i^t - 4u_i^{t-1} + u_i^{t-2}}{dt^2} - \partial_{xx} \frac{-3u^t + 6u^{t-1} + u^{t-2}}{6dt}. \quad (91)$$

Then, we obtain the equation

$$Hu^{t+1} = h^t \quad (92)$$

Observe first that dealing with spatial derivatives of fourth order is numerically challenging if they are computed in real space and therefore we transform our equation again into Fourier space

$$\hat{H}_{q_n} \hat{u}_{q_n}^{t+1} = \hat{h}_{q_n}^t \quad (93)$$

Here it is implied that the wave vector $q$ is quantized and takes the following discrete values:

$$q_n = \frac{2\pi n}{L}, \quad n = 0, 1, 2, \dots, N-1,$$

where $n$ is the mode index running from 0 to $N-1$, $L$ is the physical size of the system and $\Delta x = L/N$ is the spacing between nodes. The Fourier representation of the operator

$$\hat{H}_{q_n} = \frac{2C}{dt^2} + \frac{2q_n^2}{6dt} + Gq_n^4, \quad n = 0, 1, 2, \dots, N-1$$

is particularly suitable for numerical implementations using the FFT algorithm, as the corresponding operations in Fourier space for a diagonal operator can be performed efficiently. In particular, the operator $\hat{H}_{q_n}$ can be easily inverted. We can then write the solution of the problem in the form

$$\hat{u}_{q_n}^{t+1} = \hat{H}_{q_n}^{-1} \hat{h}_{q_n}^t. \quad (94)$$

The nonlinear terms appearing in the RHS of equation (94) are first calculated in real space and transformed into Fourier space, and their derivative is calculated as $\partial_x[2f^t - f^{t-1}] = \mathcal{F}^{-1}(iq_n\mathcal{F}([2f^t - f^{t-1}]))$, where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote direct and inverse Fourier transforms, respectively.