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

HAL
open science

# 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. 2025. hal-05148156

HAL Id: hal-05148156

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

Preprint submitted on 7 Jul 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inertia-induced power-law scaling in martensites

O. U. Salman¹,², A. Finel², and L. Truskinovsky³

¹Laboratoire des Sciences des Procédés et des Matériaux (LSPM), CNRS,
Université Sorbonne Paris Nord, UPR 3407, F-93430, Villetaneuse, France

²Laboratoire d'Etudes des Microstructures, ONERA,
29 Avenue de la Division Leclerc, 92322 Chatillon, France and

³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

(Dated: July 7, 2025)

Martensites subjected to quasistatic deformation are known to exhibit power law distributed acoustic emission in a broad range of scales, however, the origin of the observed scaling behavior and the mechanism of self-organization towards apparent criticality remains obscure. Here we argue that the power law structure of intermittent fluctuations can be at least partially interpreted as an effect of 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 can produce not only qualitative but also quantitative agreement with experiment.

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

# 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 are then examples of inherently nonlinear complex systems with intermittent response resulting from cascades 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 usually interpreted either in the framework of the concept 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 nonconvexity in producing complex multiscale microstructures [36–38].

It has been long 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 an important information about the mechanism of martensitic transformations. Thus, it reveals that between consecutive avalanches the transformation is suppressed while the system deforms purely elastically. 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 individual avalanches, the amplitudes and durations of AE hits (induced by individual elastic instability events) exhibit power law behavior. The value of the corresponding exponents were shown to depend only on crystal symmetry which points towards the universality of the underlying nonequilibrium 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 disorder [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 the parallel 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 concurrently 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 the 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 makes the 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 the present 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 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 acoustic emission [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 towards criticality and scaling has been long 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 response which is ultimately behind the observed intermittency [92, 93, 96–98]. An interesting salient 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 also long realized [67, 70, 99, 100]. Thus, already in an early study 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 [101]. Inertial effects were then shown to be behind the development of intrmittency 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 in [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, 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 ad 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]. We note 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 the present 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 twofold. 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₆₈.₈Pd₃₁.₂ undergoing a martensitic cubic-to-tetragonal transition. We show that even a very 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 exponent values very close to the ones measured in acoustic emission (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 vs damping, is chosen in an appropriate range. Interestingly, the ability of such a model to generate 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 prototypical 1D lattice model is introduced in Section II. We analyze separately damped, undamped and intermediate 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. Our 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. Our 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 manages to reach 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 final Section VI where we also formulate some open questions.

## II. ONE DIMENSIONAL DISCRETE MODEL

Consider a one dimensional chain composed of mass particles linked by elastic springs. Suppose that each par-
ticle interacts with nearest (NN) and next to nearest neighbors (NNN) on each side. Denote 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 nearest neighbor (NN) and next to nearest neighbor (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 piece wise 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

$$ \delta \mathcal{L} / \delta u = \partial \mathcal{L} / \partial u - \partial / \partial t (\partial \mathcal{L} / \partial \dot{u}) \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 = ϵN and times using the viscous time scale τ = γ/κ, 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) - \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 vs 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 the system (12) under the assumption that the two non-dimensional parameters V and δ are small but finite. The limit V → 0 corresponds to quasi-static driving and in our numerical experiments we use the value V ~ 10⁻⁵. In the limit δ → 0 we obtain continuum model, however, in all our numerical experiments we use a small but finite value δ ~ 10⁻³. 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⁻⁴.

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 ~ 10 ~ 10⁵. 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 FFT method. The numerical algorithm is discussed in more detail in Appendix A.

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, ..., N and N = T/Δ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 any more. 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, \quad (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 duration of same avalanche

$$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 the 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 V_{k+s}, \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 both, the auto-correlation function

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

We'll 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].

## A. 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 quasistatically. The details of our numerical implementation can be found in Appendix A.

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 Fig. 1. In the horizontal yielding segments of the
<figure><img src="image_3.png" /><figcaption>FIG. 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_4.png" /><figcaption>FIG. 2: Two consecutive individual avalanches in the overdamped regime which is shown in Fig. 1. Parameters: $N = 1000$ and $C = 80$.</figcaption></figure>

hysteresis cycle the system exhibits a succession of almost equal size stress drops. The typical individual avalanches, responsible for such drops, are shown in Fig. 2. The avalanches are separated by the quiescent intervals and each avalanche has roughly the same shape with an exponentially decaying tail. While 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 Fig. 1 is the absence of the nucleation peak predicted theoretically for 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), 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 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-a-vis the nucleation peak in a purely viscoelastic overdamped system.

<figure><img src="image_5.png" /><figcaption>FIG. 3: Power spectrum of the signal $V(t)$ in the overdamped regime which is shown in Fig. 1. Here $C = 80$ and $N = 1000$.</figcaption></figure>

One can see that transformational yield in Fig. 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 Onsagerian dissipative potential operative at the level of individual mass points [122, 125].

The power spectrum (22) computed in this regime is shown in Fig. 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.

<figure><img src="image_6.png" /><figcaption>FIG. 4: Cumulative probability distribution representing statistics of avalanches in the overdamped regime which is shown in Fig. 1. Here $C = 80$ and $N = 1000$. Inset: several consecutive avalanches occurred in a short time interval.</figcaption></figure>

To explain these observations we turn to Fig. 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 (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 Fig. 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 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 nonzero inertia in the system.

Finally, in Fig. 4 we show the cumulative probability distribution 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-a-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].

## B. 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.
<figure><img src="image_7.png" /><figcaption>FIG. 5: Evolution of an undamped discrete chain subjected to quasistatic loading in a hard device: the first loading cycle. Here $N = 1000$ and $C = 100000$. Inset: a zoom on the kinetic energy evolution in the pre-transition time interval.</figcaption></figure>

The response shown in Fig. 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 in (31) that $\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 Fig. 5 the originally 'cold' system first deforms homogeneously while maintaining an affine configuration. One can see, 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.

<figure><img src="image_8.png" /><figcaption>FIG. 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>

As loading continues beyond the ensuing massive instability, we observe a quasi-deterministic average stress-strain response with superimposed chaotic fluctuations, see Fig. 6. In accordance with [126] all modes acquire exactly 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 Fig. 6, at $C = 100000$ 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}. \quad (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.

The corresponding power spectrum (22), which is not accessible by purely thermodynamic analysis, is shown in Fig. 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.

<figure><img src="image_9.png" /><figcaption>FIG. 7: Power spectrum of the signal $V(t)$ for the undamped discrete chain subjected to quasistatic cyclic loading in a hard device, which is shown in Fig. 6. Here $C = 100000$ and $N = 1000$.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>FIG. 8: Cumulative probability distribution for the magnitude of the effective energy dissipation during individual avalanches in the undamped discrete chain subjected to quasistatic cyclic loading in a hard device, which is shown in Fig. 6. Inset shows a fragment of the actual time series. Here $C = 100000$ and $N = 1000$.</figcaption></figure>

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 no underdamped.

### C. 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 Fig. 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].

According to Fig. 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
<figure><img src="image_11.png" /><figcaption>FIG. 9: Stress-strain response of a quasistatically driven chain with inertia and dissipation balanced at $C = 8000$. Here $N = 1000$.</figcaption></figure>

at an almost constant level, we observe a broad, heavy-tailed distribution of stress drops.

<figure><img src="image_12.png" /><figcaption>FIG. 10: Power spectrum of the underdamped system shown in Fig. 9. Here $C = 8000$ and $N = 1000$.</figcaption></figure>

The computed power spectrum is shown in Fig. 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 because at these scales dissipation is negligible and the structure of avalanches is shaped by dynamics as much as by dissipation.

Finally, we show in Fig. 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 Fig. 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)$$

now with the exponent $\kappa \approx 2.2$, see Fig. 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 Fig. 12(b). The revealed ubiquity of power law behaviors corroborates the idea that the system develops in this range robust scale free correlations.
<figure><img src="image_13.png" /><figcaption>FIG. 11: Cumulative probability distribution in log-log scale of the energy dissipated by avalanches in the underdamed system loaded quasistatically in a hard device and shown in Fig. 9. Here C = 8000 and N = 1000. Inset: a fragment of the actual time series.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 12: Statistical signatures of avalanche behavior for the system shown in Fig. 9: (a) Cumulative probability distribution of avalanche durations. (b) Cumulative probability distribution of avalanche sizes. (c) Cloud plot showing joint distribution of avalanche sizes and durations and revealing strong correlation between these two quantities.</figcaption></figure>

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 Fig. 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 [114, 115, 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 Fig. 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 cumulative probability distribution (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 Fig. 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 cumulative probability distributions shown in Fig. 12(a)(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 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.

## D. Regime diagram

The goal of this section is to consolidate and summarize the results obtained so far. Our Fig. 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 Fig. 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 vs. 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.

<figure><img src="image_15.png" /><figcaption>FIG. 13: Regime diagram in the space of parameters showing the stress-strain response of the quasistatically 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>

According to Fig. 13, there are three main regimes with qualitatively different behavior. As we have seen, in the 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 Fig. 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.

### E. 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 [50, 130-132] as 'sweeping' of an instability. If we interpret transformational yield as such an instability, we see that 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 so on. Various prototypical models of the implied feedback allowing the system to mainten near marginal stability, and safeguarding the 'generic' nature of the scaling behavior, were studied in [57, 124].

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 self-organized criticality (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 fronts [125, 140].

### III. 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 (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 $\mathcal{E} = \mathcal{K} + \mathcal{F}$, where the (dimensionless) kinetic energy is

$$ \mathcal{K} = \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 (12) is that in (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 one-dimensional 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.

## IV. TWO DIMENSIONAL CONTINUUM MODEL

The goal of this section is to analyze a straightforward 2D analog of (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 meaningless example of a martensitic phase transition in 2D involves square and rectangle phases [72, 112, 148]. To describe the corresponding transformation, first turn to dimensional variables and again write the total energy in the form $\mathcal{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, \quad (48) $$
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. \quad (49) $$

Note first that the elastic energy density in (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). \quad (50) $$

As it is clear from (49), we assumed that the strain component $e_2$ is the primary order parameter. The first two harmonic terms in (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 (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 = 0$ and the corresponding term is mentioned here only because it will reappear in the three-dimensional model. We further assumed that the regularizing gradient energy term in (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} e_i^2(\mathbf{r}, t) d^2x, \quad (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}. \quad (52) $$

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

$$ \tau = \frac{t}{t_0}. \quad (53) $$

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}. \quad (54) $$

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

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

where $\tau$ is the dimensionless temperature. In Fig. 14 we show that, for instance at $\tau = 1$, which is the point of a 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 A, 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.
<figure><img src="image_16.png" /><figcaption>FIG. 14: The elastic energy density $\phi(e_2)$ at different values of the temperature $\tau$.</figcaption></figure>

Consider first an almost *undemped* regime with $\tilde{\rho} = 1$ and $\tilde{\gamma}_2 = \tilde{\gamma}_1 = \tilde{\gamma}_3 = 0$. This choice corresponds to the limit $C \to \infty$ in the 1D model. In Fig. 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 were 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 Fig. 15(a). Our Fig. 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_17.png" /><figcaption>FIG. 15: Time evolution of a pre-existent martensitic embryo in an almost 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 Fig. 14 for the identification of particular variants.</figcaption></figure>

As we see in Fig. 15(a) the transformation begins around the infinitesimal embryo in a form of a single variant of the martensite effectively chosen by the type of the inserted imperfection. Already in Fig. 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 Fig. 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 Fig. 15(d). As such multi-scale microstructure formation process reaches the boundaries of the domain, austenite phase largely disappears, see Fig. 15(e). Eventually the microstructure stabilizes by coarsening and all independently developing transformed zones finally merge into a single multiscale texture, see Fig. 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 exclusively with 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 *overdamped* regime with $\tilde{\rho} = 1$ and $\gamma_2 = 1$, $\gamma_1 = \gamma_3 = 0.5$. As we see in Fig. 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 Fig. 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 Fig. 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
<figure><img src="image_18.png" /><figcaption>FIG. 16: Time evolution of a pre-existent martensitic embryo in an underamped 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 Fig. 14 for the identification of particular variants.</figcaption></figure>

pattern is regular and the whole pattern shows much less complexity than the multi-scale texture observed in the undamped regime.

<figure><img src="image_19.png" /><figcaption>FIG. 17: Post nucleation pattern in a 2D almost fully damped continuum model, see Fig. 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 = 0$ in Eq. 44). In this limit we obtain time-dependent Ginzburg-Landau model with viscous damping, see for instance [72, 101]. In this limit the breakdown of an unstable austenite phase results a predictable emergence of a coarse highly regular twin microstructure, see Fig. 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. Moreover, 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 multi-scale 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.

V. 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^3 x \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}-type and <110>-type directions:

$$ 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 (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 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.

<figure><img src="image_20.png" /><figcaption>FIG. 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>
A. 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}, U_2 = \begin{pmatrix} a^* & 0 & 0 \\ 0 & c^* & 0 \\ 0 & 0 & a^* \end{pmatrix}, U_3 = \begin{pmatrix} a^* & 0 & 0 \\ 0 & a^* & 0 \\ 0 & 0 & c^* \end{pmatrix}, \qquad (61) $$

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 (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 [113, 156, 157]. 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 \simeq 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 Fig. 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_0 f_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} \text{ 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 \text{ kg m}^{-3}$ which suggests that the value $\tilde{\rho} \sim 1$ can be used as representing a typical underdamped regime.
## B. 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 A2 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 $\dot{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}, \quad (63)$$

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

<figure><img src="image_21.png" /><figcaption>FIG. 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.</figcaption></figure>

The unfolding of the transformation processes is illustrated in Fig. 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 Fig. 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 (variants 1 and 2 shown in green and red colors, respectively), see Fig. 19(b). The texture of internally twinned martensite laminates continue to grow dynamically while progressively generating more and more complex hierarchical pattern, see Fig. 19(c,d). Finally, the austenite is completely transformed into the variants of martensite which self-organize to form a multi-scale structure, see Fig. 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, as we have seen, 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_22.png" /><figcaption>FIG. 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>

### C. Statistics of avalanches

Consider next the issues of intermittency and scaling in the 3D model. In Fig. 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 Fig. 20(a,b,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 Fig. 20 (a). The system size avalanche, however, is formed itself by a large number of small bursts representing pre- and after-shocks Fig. 20 (b,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.

<figure><img src="image_23.png" /><figcaption>FIG. 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>

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. The avalanche size *S* is calculated as in 1D using Eqs. 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, combining 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 Fig. 21. Over four decades it is clearly of power
<figure><img src="image_24.png" /><figcaption>FIG. 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_25.png" /><figcaption>FIG. 23: Correlations between avalanche energies *E* and avalanche amplitudes in 3D underdamped continuum model.</figcaption></figure>

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 Fig. 22 (a,b). here we show 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)$$

which exponents $\tau = 2.4$ and $\kappa = 2.15$. Finally, in Fig. 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.

<table><caption>TABLE I: 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>

Given that we fixed the simulation parameters using experimental measurements available in the literature for the cubic to tetragonal martensitic transformation in Fe<sub>68.8</sub>Pd<sub>31.2</sub>, 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 I, 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₆₈.₈Pd₃₁.₂ 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 linked to the underdamped nature of dynamics.

## VI. CONCLUSIONS

The goal of this study was to identify the origin of the power law acoustic emission (AE) in martensites subjected to quasistatic loading. Despite numerous proposals regarding the possible mechanism 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 systematically the possibility that the implied 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 observed 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 all 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 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 nearest neighbors (NN) are characterized by a nonconvex potential, harmonic next to nearest neighbor (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 Random Field Ising Model (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 as 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 intermittent 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 completely 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 the present study would be to take into account the possibility 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 associated 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 phenomena 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 the present 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 to 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.

# VII. ACKNOWLEDGMENTS

The authors are grateful to Francisco Perez Reche for valuable comments and suggestions. 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.

# Appendix A: Spatial and temporal discretisation schemes

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

## 1. 1D discrete problem

To solve (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^{-iq_n j \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 (A1) $$

$$ \frac{d\hat{u}(q)}{dt} = y_2 \qquad (A2) $$

and rewrite our second order equation as a first-order system

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

$$ \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 (A4) $$

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 $\hat{\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 (A5) $$

$$ 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 (A6) $$

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

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}'(q)}{\delta} + L(q) \left(y_{1,n} + \frac{\Delta t}{2} k_{1,1}\right) \right] \qquad (A8) $$

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

$$ 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}'(q)}{\delta} + L(q) \left(y_{1,n} + \frac{\Delta t}{2} k_{2,1}\right) \right] \qquad (A10) $$

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

$$ 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}'(q)}{\delta} + L(q) \left(y_{1,n} + \Delta t k_{3,1}\right) \right]. \qquad (A12) $$

Here, $\hat{\phi}'(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}'(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 (A13) $$

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

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

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

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 in [151].

## 2. 1D 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 (A17) $$

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 (A18) $$

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

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} - dt f_t^{t+1} + \frac{dt^2 f_{tt}^{t+1}}{2}, \qquad (A20) $$

$$ f^{t-1} \approx f^{t+1} - 2dt f_t^{t+1} + 2dt^2 f_{tt}^{t+1}. \qquad (A21) $$

We can also write

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

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 (A23) $$

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 (A24) $$
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 (A25)$$

Then, we obtain the equation

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

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 (A27)$$

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 Fast Fourier Transform (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 (A28)$$

The nonlinear terms appearing in the RHS of Eq. A28 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.

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

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

[3] Antoni Planes, Lluís Mañosa, and Eduard Vives. Acoustic emission in martensitic transformations. J. Alloys Compd., 577:S699-S704, November 2013.

[4] Francisco J. Perez-Resche. Modelling avalanches in martensites. In Avalanches in functional materials and geophysics, pages 99-136. 2017.

[5] X. Balandraud, N. Barrera, P. Biscari, M. Grédiac, and G. Zanzotto. Strain intermittency in shape-memory alloys. Phys. Rev. B, 91(17):174111, 2015.

[6] J. P. Sethna, M. K. Bierbaum, K. A. Dahmen, C. P. Goodrich, J. R. Greer, L. X. Hayden, and S. Zapperi. Deformation of crystals: Connections with statistical physics. Annu. Rev. Mater. Res., 47(1):217-246, 2017.

[7] M. J. Alava, P. K. Nukala, and S. Zapperi. Statistical models of fracture. Adv. Phys., 55(3-4):349-476, 2006.

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

[9] M. C. Miguel, A. Vespignani, S. Zapperi, J. Weiss, and J. R. Grasso. Intermittent dislocation flow in viscoplastic deformation. Nature, 410:667-671, 2001.

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

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

[12] Oguz Umut Salman and Lev Truskinovsky. Minimal integer automaton behind crystal plasticity. Phys. Rev. Lett., 106 (17):175503, April 2011.

[13] P Zhang, O U Salman, J Weiss, and L Truskinovsky. Variety of scaling behaviors in nanocrystalline plasticity. Phys Rev E, 102(2-1):023006, August 2020.

[14] O. U. Salman, A. Finel, R. Delville, and D. Schryvers. The role of phase compatibility in martensite. J. Appl. Phys, 111 (10), 2012.
[15] D. S. Fisher. Sliding charge-density waves as a dynamic critical phenomenon. Phys. Rev. B, 31:1396, 1985.

[16] A. Rosso, J. P. Sethna, and M. Wyart. 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, pages 277–305. 2022.

[17] M. LeBlanc, J. T. Uhl, and K. A. Dahmen. Avalanches in strained bulk metallic glasses. Phys. Rev. Lett., 112:155501, 2014.

[18] D. Denisov, K. Lorincz, J. Uhl, K. Dahmen, and P. Schall. Universality of slip avalanches in flowing granular matter. Nature Communications, 7:10641, 2016.

[19] C. Liu, E. E. Ferrero, F. Puosi, J.-L. Barrat, and K. Martens. Driving rate dependence of avalanche statistics and shapes at the yielding transition. Phys. Rev. Lett., 116:065501, 2016.

[20] P. Y. Chan, G. Tsekenis, J. Dantzig, K. A. Dahmen, and N. Goldenfeld. Plasticity and dislocation dynamics in a phase field crystal model. Phys. Rev. Lett., 105:015502, 2010.

[21] N. Friedman, A. T. Jennings, G. Tsekenis, J.-Y. Kim, M. Tao, J. T. Uhl, J. R. Greer, and K. A. Dahmen. Statistics of dislocation slip avalanches in nanosized single crystals show tuned critical behavior. Phys. Rev. Lett., 109:095507, 2012.

[22] Per Bak, Chao Tang, and Kurt Wiesenfeld. Self-organized criticality: An explanation of the 1/f noise. Phys. Rev. Lett., 59(4):381–384, 1987.

[23] R. Dickman, M. A. Muñoz, A. Vespignani, and S. Zapperi. Paths to self-organized criticality. Brazilian Journal of Physics, 30:27–41, 2000.

[24] H. J. Jensen. *Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems*, volume 10. Cambridge University Press, 1998.

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

[26] Gunnar Pruessner. *Self-Organised Criticality: Theory, Models and Characterisation*. Cambridge University Press, 2012.

[27] D. Sornette. Critical Phenomena in Natural Sciences: Chaos, Fractals, Self-Organization and Disorder: Concepts and Tools. Springer Science and Business Media, 2006.

[28] C. Gros. Self-organized criticality. In *Complex and Adaptive Dynamical Systems: A Comprehensive Introduction*, pages 203–239. Springer International Publishing, 2024.

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

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

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

[32] M. Ovaska, A. Lehtinen, M. J. Alava, L. Laurson, and S. Zapperi. Excitation spectra in crystal plasticity. Phys. Rev. Lett., 119(26):265501, 2017.

[33] S. Franz and S. Spigler. Mean-field avalanches in jammed spheres. Phys. Rev. E, 95(2):022139, 2017.

[34] Zenji Nishiyama. Martensitic transformation. Elsevier, 2012.

[35] AL Roitburd. Martensitic transformation as a typical phase transformation in solids. In *Solid state physics*, volume 33, pages 317–390. Elsevier, 1978.

[36] Richard D James and Kevin F Hane. Martensitic transformations and shape-memory materials. Acta Mater, 48(1): 197–222, 2000.

[37] Armen G Khachaturyan. *Theory of structural transformations in solids*. Courier Corporation, 2013.

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

[39] Francisco-José Pérez-Reche, Marcelo Stipcich, Eduard Vives, Lluís Mañosa, Antoni Planes, and Michel Morin. Kinetics of martensitic transitions in cu-al-mn under thermal cycling: Analysis at multiple length scales. Phys. Rev. B, 69:064101, Feb 2004.

[40] U. Chandni, A. Ghosh, H. S. Vijaya, and S. Mohan. Criticality of tuning in athermal phase transitions. Phys. Rev. Lett., 102(2):025701, 2009.

[41] M. L. Rosinberg and E. Vives. Metastability, hysteresis, avalanches, and acoustic emission: Martensitic transitions in functional materials. In *Disorder and Strain-Induced Complexity in Functional Materials*, pages 249–272. Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.

[42] E. Vives, D. E. Soto Parra, A. Planes, L. Mañosa, R. Romero, R. S. Edwards, and S. Dixon. Acoustic emission avalanches in martensitic transitions: new perspectives for the problem of source location. Solid State Phenom., 172:144–149, 2011.

[43] M. C. Gallardo, J. Manchado, F. J. Romero, J. Del Cerro, E. K. Salje, A. Planes, and M. Stipcich. 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, 81(17):174102, 2010.

[44] D. L. Beke, L. Daróczi, L. Z. Tóth, M. K. Bolgár, N. M. Samy, and A. Hudák. Acoustic emissions during structural changes in shape memory alloys. Metals, 9(1):58, 2019.

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

[46] Eduard Vives, Daniel Soto-Parra, Lluís Mañosa, Ricardo Romero, and Antoni Planes. Driving-induced crossover in the avalanche criticality of martensitic transitions. Phys. Rev. B, 80:180101, Nov 2009.

[47] Genís Torrents, Xavier Illa, Eduard Vives, and Antoni Planes. Geometrical model for martensitic phase transitions: Understanding criticality and weak universality during microstructure growth. Phys. Rev. E, 95(1-1):013001, January 2017.
[48] D L Beke, M K Bolgár, L Z Tóth, and L Daróczi. On the asymmetry of the forward and reverse martensitic transformations in shape memory alloys. J. Alloys Compd., 741:106-115, April 2018.

[49] Yintao Song, Xian Chen, Vivekanand Dabade, Thomas W Shield, and Richard D James. Enhanced reversibility and unusual microstructure of a phase-transforming material, *Nature*, 502(7469):85–88, October 2013.

[50] J. P. Sethna, K. A. Dahmen, and C. R. Myers. Crackling noise. *Nature*, 410:242, 2001.

[51] James P. Sethna, Karin Dahmen, Sivan Kartha, James A. Krumhansl, Bruce W. Roberts, and Joel D. Shore. Hysteresis and hierarchies: Dynamics of disorder-driven first-order phase transformations. Phys. Rev. Lett., 70(21):3347-3350, May 1993.

[52] Olga Perković, Karin Dahmen, and James P. Sethna. Avalanches, barkhausen noise, and plain old criticality. Phys. Rev. Lett., 75(24):4528-4531, Dec 1995.

[53] Karin Dahmen and James P. Sethna. Hysteresis, avalanches, and disorder-induced critical scaling: A renormalization-group approach. Phys. Rev. B, 53(22):14872-14905, Jun 1996.

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

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

[56] Eduard Vives and Antoni Planes. Hysteresis and avalanches in the random anisotropy ising model. Phys. Rev. B, 63(13):134431, Mar 2001.

[57] Francisco-José Pérez-Reche, Lev Truskinovsky, and Giovanni Zanzotto. Driving-induced crossover: From classical criticality to self-organized criticality. Phys. Rev. Lett., 101(23):27, December 2008.

[58] J. P. Sethna, K. A. Dahmen, and O. Perkovic. Random-field Ising models of hysteresis. Academic Press, New York, 2006.

[59] Francisco-José Pérez-Reche, Lev Truskinovsky, and Giovanni Zanzotto. Training-induced criticality in martensites. Phys. Rev. Lett., 99(7), August 2007.

[60] F. J. Perez-Reche, C. Triguero, G. Zanzotto, and L. Truskinovsky. Origin of scale-free intermittency in structural first-order phase transitions. Phys. Rev. B, 94(14):144102, 2016.

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

[62] Lluís Carrillo, Lluís Mañosa, Jordi Ortín, Antoni Planes, and Eduard Vives. Experimental evidence for universality of acoustic emission avalanche distributions during structural transitions. Phys. Rev. Lett., 81(9):1889-1892, August 1998.

[63] Erell Bonnot, Ricardo Romero, Lluís Mañosa, Eduard Vives, and Antoni Planes. Elastocaloric effect associated with the martensitic transition in shape-memory alloys. Phys. Rev. Lett., 100(12):125901, March 2008.

[64] J. D. Eshelby. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the Royal Society of London A, 241:376, 1957.

[65] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet. Elastic consequences of a single plastic event: A step towards the microscopic modeling of the flow of yield stress fluids. Europhysics Journal E, 15:371, 2004.

[66] S. Rossi, G. Biroli, M. Ozawa, and G. Tarjus. Far-from-equilibrium criticality in the random-field ising model with eshelby interactions. Phys. Rev. B, 108(22):L220202, 2023.

[67] PC Clapp. How would we recognize a martensitic transformation if it bumped into us on a dark & austy night? Le Journal de Physique IV, 5(C8):C8-11, 1995.

[68] Yucong Miao and Joost J Vlassak. Explosive martensitic transformation of supercooled austenite in cuzr-based thin-film shape memory alloys. Acta Mater., 200:162-170, 2020.

[69] Stefan Schwabe, Klara Lünser, Daniel Schmidt, Kornelius Nielsch, Peter Gaal, and Sebastian Fähler. What is the speed limit of martensitic transformations? Sci. Technol. Adv. Mater., 23(1):633-641, 2022.

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

[71] Lev Truskinovsky. Nucleation and growth in elastodynamics. In Dynamics of Crystal Surfaces and Interfaces, pages 185-197. Springer, 2002.

[72] T. Lookman, S. R. Shenoy, K. Ø. Rasmussen, A. Saxena, and A. R. Bishop. Ferroelastic dynamics and strain compatibility. Phys. Rev. B, 67(2):024114, 2003.

[73] Lev Truskinovsky and Anna Vainchtein. Dynamics of martensitic phase boundaries: discreteness, dissipation and inertia. Contin. Mech. Thermodyn., 20(2):97-122, May 2008.

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

[75] J-Y Cho, AV Idesman, VI Levitas, and Taehyo Park. Finite element simulations of dynamics of multivariant martensitic phase transitions based on ginzburg-landau theory. International Journal of Solids and Structures, 49(14):1973-1992, 2012.

[76] Xiaoying Liu, Daniel Schneider, Martin Reder, Paul W Hoffrogge, and Britta Nestler. Modeling of martensitic phase transformation accounting for inertia effects. Int. J. Mech. Sci, 278:109443, 2024.

[77] Eduard Vives, Jordi Ortín, Lluís Mañosa, Ismael Ràfols, Ramon Pérez-Magrané, and Antoni Planes. Distributions of avalanches in martensitic transformations. Phys. Rev. Lett., 72(11):1694-1697, 1994.

[78] Ll Mañosa, A Planes, D Rouby, M Morin, P Fleischmann, and J L Macqueron. Acoustic emission field during thermoelastic martensitic transformations. Appl. Phys. Lett., 54(25):2574-2576, June 1989.

[79] Jordi Baró, Álvaro Corral, Xavier Illa, Antoni Planes, Ekhard KH Salje, Wilfried Schranz, Daniel E Soto-Parra, and Eduard Vives. Statistical similarity between the compression of a porous material and earthquakes. Phys. Rev. Lett.
110(8):088702, 2013.

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

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

[82] GG Kocharyan. Nucleation and evolution of sliding in continental fault zones under the action of natural and man-made factors: A state-of-the-art review. *Izvestiya, Physics of the Solid Earth*, 57:439–473, 2021.

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

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

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

[86] J. M. Carlson and J. S. Langer. Properties of earthquakes generated by fault dynamics. *Phys. Rev. Lett.*, 62:2632, 1989.

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

[88] G. A. Held, D. H. Solina, D. T. Keane, W. J. Haag, P. M. Horn, and G. Grinstein. Experimental study of critical-mass fluctuations in an evolving sandpile. *Phys. Rev. Lett.*, 65:1120, 1990.

[89] H. M. Jaeger, C.-h. Liu, and S. R. Nagel. Relaxation at the angle of repose. *Phys. Rev. Lett.*, 62:40, 1989.

[90] M C Marchetti. Depinning and plasticity of driven disordered lattices. In M Carmen Miguel and Miguel Rubi, editors, *Jamming, Yielding, and Irreversible Deformation in Condensed Matter*, Lecture Notes in Physics, pages 137–157. Springer-Verlag, 2006.

[91] D. V. Denisov, K. A. Lorincz, W. J. Wright, T. C. Hufnagel, A. Nawano, X. Gu, J. T. Uhl, K. A. Dahmen, and P. Schall. Universal slip dynamics in metallic glasses and granular matter – linking frictional weakening with inertial effects. *Scientific Reports*, 7:43376, 2017.

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

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

[94] T. W. de Geus and M. Wyart. Short-range depinning in the presence of velocity-weakening. *arXiv preprint*, 2024.

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

[96] K. M. Salerno, C. E. Maloney, and M. O. Robbins. Avalanches in strained amorphous solids: does inertia destroy critical behavior? *Phys. Rev. Lett.*, 109:105703, 2012.

[97] K. M. Salerno and M. O. Robbins. Effect of inertia on sheared disordered solids: Critical scaling of avalanches in two and three dimensions. *Phys Rev E*, 88:062206, 2013.

[98] K. Karimi, E. E. Ferrero, and J.-L. Barrat. Inertia and universality of avalanche statistics: the case of slowly deformed amorphous solids. *Phys Rev E*, 95:013003, 2017.

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

[100] ACE Reid and RJ Gooding. Pattern formation in a 2d elastic solid. *Physica A*, 239(1-3):1–10, 1997.

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

[102] F. J. Elmer. Avalanches in the weakly driven frenkel-kontorova model. *Phys Rev E*, 50(6):4470, 1994.

[103] X Ding, T Lookman, Z Zhao, A Saxena, J Sun, and EKH Salje. Dynamically strained ferroelastics: Statistical behavior in elastic and plastic regimes. *Phys. Rev. B*, 87(9):094109, 2013.

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

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

[106] Arya Paul, Jayee Bhattacharya, Surajit Sengupta, and Madan Rao. Non-affine deformation in microstructure selection in solids ii: Elastoplastic theory for the dynamics of solid state transformations. *J. Phys.: Condens. Matter*, 20(36):365211, 2008.

[107] John M Ball and RD James. Incompatible sets of gradients and metastability. *Arch. Ration. Mech. Anal.*, 218:1363–1416, 2015.

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

[109] Zhiyong Zhang, Richard D James, and Stefan Müller. Energy barriers and hysteresis in martensitic phase transformations. *Acta Mater.*, 57(15):4332–4352, 2009.

[110] KC Atli, BE Franco, I Karaman, D Gaydosh, and RD Noebe. Influence of crystallographic compatibility on residual strain of tini based shape memory alloys during thermo-mechanical cycling. *Mater. Sci. Eng. A*, 574:9–16, 2013.

[111] Yuanchao Yang, Yangyang Xu, Yumei Zhou, Xiangdong Ding, Jun Sun, Turab Lookman, and Dezhen Xue. Nonthermoe-lastic martensitic features in ideal martensites due to volume effects. *Phys. Rev. B*, 108(2):024102, 2023.

[112] S. R. Shenoy, T. Lookman, A. Saxena, and A. R. Bishop. Martensitic textures: Multiscale consequences of elastic compatibility. *Phys. Rev. B*, 60(18):R12537–R12541, 1999.

[113] Rajeev Ahluwalia, Turab Lookman, and Avadh Saxena. Dynamic strain loading of cubic to tetragonal martensites. *Acta Mater.*, 54(8):2109–2120, 2006.
[114] János Kertész and László B Kiss. The noise spectrum in the model of self-organised criticality. *Journal of Physics A: Mathematical and General*, 23(9):L433, 1990.

[115] Matthew C. Kuntz and James P. Sethna. Noise in disordered systems: The power spectrum and dynamic exponents in avalanche models. Phys. Rev. B, 62:11699-11708, Nov 2000.

[116] Lasse Laurson, Mikko J Alava, and Stefano Zapperi. Power spectra of self-organized critical sandpiles. J. Stat. Mech., 2005(11):L11001, 2005.

[117] Lasse Laurson and Mikko J Alava. 1/f noise and avalanche scaling in plastic deformation. Phys. Rev. E, 74(6):066106, 2006.

[118] Manoj Kumar Nandi, Alessandro Sarracino, Hans J Herrmann, and Lucilla de Arcangelis. Scaling of avalanche shape and activity power spectrum in neuronal networks. Phys. Rev. E, 106(2):024304, 2022.

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

[120] Steven M Kay. Fundamentals of statistical signal processing: estimation theory. Prentice-Hall, Inc., 1993.

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

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

[123] Alexander Mielke and Lev Truskinovsky. From discrete visco-elasticity to continuum rate-independent plasticity: Rigorous results. Arch. Ration. Mech. Anal., 203(2):577-619, September 2011.

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

[125] Lev Truskinovsky and Anna Vainchtein. Kinetics of martensitic phase transitions: Lattice model. SIAM J. Appl. Math., 66(2):533-553, 2005.

[126] Yalchin R Efendiev and Lev Truskinovsky. Thermalization of a driven bi-stable FPU chain. Continuum Mech. Thermodyn., 22(6-8):679-698, September 2010.

[127] Bruce J West and Michael F Shlesinger. On the ubiquity of 1/f noise. Int. J. Mod. Phys. B, 03(06):795-819, June 1989.

[128] E Milotti. 1/f noise: a pedagogical review. arXiv: Classical Physics, April 2002.

[129] Stefano Polizzi, Francisco-José Pérez-Reche, Alain Arneodo, and Françoise Argoul. Power-law and log-normal avalanche size statistics in random growth processes. Phys Rev E, 104(5):L052101, 2021.

[130] Didier Sornette. Sweeping of an instability: an alternative to self-organized criticality to get powerlaws without parameter tuning. Journal de Physique I, 4(2):209-221, 1994.

[131] Didier Sornette. Power laws without parameter tuning: An alternative to self-organized criticality. Phys. Rev. Lett., 72 (14):2306, 1994.

[132] Eric Vivès, J Ortín, Ll Mañosa, I Ràfols, and A Planes. Distribution of acoustic emission avalanches in martensitic transformations. Le Journal de Physique IV, 5(C2):C2-59, 1995.

[133] J S Urbach, R C Madison, and J T Markert. Interface depinning, self-organized criticality, and the barkhausen effect. Phys. Rev. Lett., 75(2):276-279, July 1995.

[134] Gianfranco Durin, Giorgio Bertotti, and Alessandro Magni. Fractals, scaling and the question of self-organized criticality in magnetization processes. Fractals, 03(02):351-370, June 1995.

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

[136] Oguz Umut Salman, Benson Muite, and Alphonse Finel. Origin of stabilization of macrotwin boundaries in martensites. Eur. Phys. J. B, 92:1-9, 2019.

[137] J M Carlson, J S Langer, B E Shaw, and C Tang. Intrinsic properties of a burridge-knopoff model of an earthquake fault. Phys. Rev. A, 44(2):884-897, July 1991.

[138] Jean Schmittbuhl, Jean-Pierre Vilotte, and Stéphane Roux. Velocity weakening friction: A renormalization approach. J. Geophys. Res., 101(B6):13911-13917, 1996.

[139] F J Elmer. Is self-organized criticality possible in dry friction? In Physics of Sliding Friction, pages 433-447. Springer Netherlands, Dordrecht, 1996.

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

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

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

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

[144] Kaushik Bhattacharya, Nikan B Firoozye, Richard D James, and Robert V Kohn. Restrictions on microstructure. Proc. R. Soc. Edinb. A, 124(5):843-878, 1994.

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

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

[147] Yury Grabovsky and Lev Truskinovsky. Rigidity-induced critical points. Phys Rev E, 110(6):064114, 2024.

[148] S H Curnoe and A E Jacobs. Twin wall of proper cubic-tetragonal ferroelastics. Phys. Rev. B: Condens. Matter Mater. Phys., 62(18):R11925-R11928, 2000.

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

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

[151] O U Salman. Modeling of spatio-temporal dynamics and patterning mechanisms of martensites by phase-field and Lagrangian methods. PhD thesis, Université Pierre et Marie Curie, 2009.

[152] JM Ball, Philip J Holmes, RD James, RL Pego, and PJ Swart. On the dynamics of fine structure. J. Nonlinear Sci., 1: 17-70, 1991.

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

[154] S. H. Curnoe and A. E. Jacobs. Time evolution of tetragonal-orthorhombic ferroelastics. Phys. Rev. B, 64(6):064101, 2001.

[155] A. E. Jacobs, S. H. Curnoe, and R. C. Desai. Simulations of cubic-tetragonal ferroelastics. Phys. Rev. B, 68(22):224104, 2003.

[156] S Muto, R Oshima, and F E Fujita. Elastic softening and elastic strain energy consideration in the fcc-fct transformation of fe-pd alloys. Acta Metall. Mater., 38(4):685-694, 1990.

[157] S Kartha, J A Krumhansl, J P Sethna, and L K Wickham. Disorder-driven pretransitional tweed pattern in martensitic transformations. Phys. Rev. B, 52(2):803-822, July 1995.

[158] M Sato, B H Grier, S M Shapiro, and H Miyajima. Effect of magnetic ordering on the lattice dynamics of FCC Fe1-xPdx. J. Phys., 12(10):2117-2129, October 1982.

[159] R Oshima, M Sugiyama, and F E Fujita. Tweed structures associated with fcc-fct transformations in fe-pd alloys. Metall. Trans., 19(4):803-810, April 1988.

[160] L R Testardi and T B Bateman. Lattice instability of high-transition-temperature superconductors. II. single-Crystal V3Si results. Phys. Rev., 154(2):402-410, February 1967.

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

[162] Mario Pitteri and Giovanni Zanzotto. Continuum models for phase transitions and twinning in crystals. Chapman and Hall/CRC, 2002.

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

[164] Mahmood Mamivand, Mohsen Asle Zaeem, and Haitham El Kadiri. A review on phase field modeling of martensitic phase transformation. Comput. Mater. Sci., 77:304-311, September 2013.

[165] Y Ni, Y M Jin, and A G Khachaturyan. The transformation sequences in the cubic→tetragonal decomposition. Acta Mater, 55(14):4903-4914, August 2007.

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

[167] Tetsu Ichitsubo, Katsushi Tanaka, Masahiro Koiwa, and Yoshihiro Yamazaki. Kinetics of cubic to tetragonal transformation under external field by the time-dependent ginzburg-landau approach. Phys. Rev. B, 62:5435-5441, Sep 2000.

[168] F E Hildebrand and C Miehe. A phase field model for the formation and evolution of martensitic laminate microstructure at finite strains. Philos. Mag. (Abingdon), 92(34):4250-4290, December 2012.

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

[170] Aaron Clauset, Cosma Rohilla Shalizi, and M E J Newman. Power-law distributions in empirical data. SIAM Rev. Soc. Ind. Appl. Math., 51(4):661-703, November 2009.

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

[172] Randall J. LeVeque. Finite Difference Methods for Ordinary and Partial Differential Equations, Steady State and Time Dependent Problems. SIAM, 2007. ISBN 978-0-898716-29-0.