Thermal and athermal crackling noise in ferroelastic nanostructures

The evolution of ferroelastic microstructures under external shear is determined by large-scale molecular dynamics simulations in two and three dimensions. Ferroelastic pattern formation was found to be almost identical in two and three dimensions, with only the ferroelastic transition temperature changing. The twin patterns generated by shear deformation depend strongly on temperature, with high wall densities nucleating under optimized temperature conditions. The dynamical tweed and mobile kink movement inside the twin walls is continuous and thermally activated at high temperatures, and becomes jerky and athermal at low temperatures. With decreasing temperature, the statistical distributions of dynamical tweed and kinks vary from a Vogel–Fulcher law P(E) ~ exp-(E/(T-TVF)) ?> to an athermal power-law distribution P(E) ~ E-ε ?>. During the yield event, the nucleation of needles and kinks is always jerky, and the energy of the jerks is power-law distributed. Low-temperature yield proceeds via one large avalanche. With increasing temperature, the large avalanche is thermally broken up into a multitude of small segments. The power-law exponents reflect the changes in temperature, even in the athermal regime.


Introduction
The dynamics of ferroelastic domain formation and straininduced de-twinning [1][2][3][4][5][6][7] often produce crackling noise [8,9]. This characteristic sound stems from the displacement discontinuities of propagating domain boundaries, which dissipate energy over a wide range of domain sizes [8,10,11] and frequencies. Crackling noise is common in nature. A piece of paper crackles when crumpled, the Earth emits intermittent noise in earthquakes, and so on. Similar behavior is observed in magnetic systems when the magnetization evolves through steps under an applied magnetic field, which generates so-called Barkhausen noise [12][13][14] when the movement is jerky.
The study of the crackling noise in ferroelastic materials is of wider importance. For geologists, it is highly relevant for the deformation of rocks. Many minerals are ferroelastic, they make up more than 30% of the Earth's crust and most of the Earth's mantle. In earthquakes, fore-shocks and aftershocks are commonly observed [15,16]. The correlations between avalanches and fore-shocks open up the possibility of being able to predict major collapses [17,18]. For materials scientists, the emerging field of 'domain boundary engineering' has generated great interest due to the functionalities of the domain boundary in ferroelastic materials [19][20][21]. The The evolution of ferroelastic microstructures under external shear is determined by large-scale molecular dynamics simulations in two and three dimensions. Ferroelastic pattern formation was found to be almost identical in two and three dimensions, with only the ferroelastic transition temperature changing. The twin patterns generated by shear deformation depend strongly on temperature, with high wall densities nucleating under optimized temperature conditions. The dynamical tweed and mobile kink movement inside the twin walls is continuous and thermally activated at high temperatures, and becomes jerky and athermal at low temperatures. With decreasing temperature, the statistical distributions of dynamical tweed and kinks vary from a Vogel-Fulcher law ~ P E E T T ( ) exp-( /( -)) VF to an athermal power-law distribution ~ ε P E E ( ) -. During the yield event, the nucleation of needles and kinks is always jerky, and the energy of the jerks is power-law distributed. Low-temperature yield proceeds via one large avalanche. With increasing temperature, the large avalanche is thermally broken up into a multitude of small segments. The power-law exponents reflect the changes in temperature, even in the athermal regime.
Key words: crackling noise, avalanche, jerks, needle domain, ferroelastic phase transition (Some figures may appear in colour only in the online journal) 3 Author to whom any correspondence should be addressed.

OPEN ACCESS
Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. domain boundaries can be superconducting [22][23][24], multiferroic [25,26] and so on. Such properties can be exclusively contained in the twin boundaries and do not exist in the bulk. The application of unique properties requires not only the functionality of the domain boundaries, but also a high density of domain boundaries. Domain boundaries are usually created by fast-temperature quench or low-temperature shear strain [27,28]. Quench and low-temperature shear lead to dynamic processes of domain evolution and can be optimized to maximize the density of domain boundaries. Ferroelastic materials also play an important role in multiferroic devices [29,30] as templates for twin structures.
Experimentally, 'crackling noise' is generated by rapid changes in microstructures [1][2][3][4][5][6][7]30], and is often measured via the acoustic emission under shear deformation or temperature changes [5]. Another experimental technique is the direct calorimetric observation of heat fluxes, thermal expansion and the non-smooth temperature evolution of the elastic moduli [1,31]. The energy spectra are often a superposition of jerky and smooth components, depending on the time evolution of twin patterns. Energy jerks often stem from abrupt pinning/de-pinning events (e.g. domain walls interact with external defects or are jammed by other twin walls).
Smooth behavior relates to the unpinned movement of twin walls. Both acoustic emission and heat-flux measurements typically show the jerky energy spectra, with power-law scaling of the jerk distribution. The power-law exponents, ε, are found in the range between 1.3 and 2.3 [1,[31][32][33]. Nevertheless, it is not easy to obtain a direct link between microstructural changes and energy emissions. Simplified theoretical works, such as Bak's sand-pile model [34][35][36][37], the random field Ising model [11,38] and renormalization group methods [39], provided a first attempt to understand avalanche behavior and estimate the energy exponents, but they cannot help decipher the underlying physical processes of changing microstructures.
During recent years, large-scale molecular dynamics (MD) simulations became possible and have contributed to the analysis of avalanches and jerks. The evolution of domain structures during shear deformation has already been simulated in some detail [40,41], and it was shown that pinning/ de-pinning events often exist without external defects but are a consequence of domain jamming. The main result was that athermal avalanches dominate at low temperatures, while thermally activated ones were found at high temperatures (Vogel-Fulcher behavior) [40,41].
In this paper, we report on how high-resolution statistical analysis of such models confirms the previous conclusions and sheds new light on the underlying mechanisms. For this purpose, we use almost the same model parameters as in [40] but increase the atomic mass to reduce the temperature of the ferroelastic phase transition. This allows us to scale the temperature relative to the transition temperature, T C . Twin pattern nucleation was found predominantly at the center of the sample, so boundary effects (open boundary conditions) are not relevant in our simulations. We find that the previous restriction to two-dimensional models is not significant because we observe almost the same transition behavior and pattern formation when we extend the model to three dimensions.
This paper is organized as follows. In section 2, we briefly introduce the model and methods for our MD simulations. Section 3 focuses on the results of the MD simulations. We first check the difference two-dimensional and three-dimensional model. The origin of the power law and exponential behaviors is shown from the analysis of the simulated jerk excitations. Then we compare the evolution of the microstructure under applied shear and the jerk energy statistics in different deformation stages where different mechanisms dominate the dynamical behavior. In section 4, we discuss the temperature effect on jerky behavior for different strain regimes.

Computer simulation
Our model is constructed to be generic for all ferroelastic materials by choosing a typical spontaneous strain that leads to ferroelastic twinning. This twinning and the mobility of the twin walls define ferroelasticity. The large-scale MD simulations are performed based on the interatomic potentials, by which the elementary steps leading to the needle domain movement and kink propagation are well reproduced. Forcefield simulations are not adopted since many atomic details would be averaged and cannot be analyzed. The potential consists of three interactions [42][43][44] in a monoatomic lattice, which are (1) harmonic nearest-neighbor interactions (elastic springs), = − U r r ( ) 20( 1) 2 , (the black springs in figure 1); (2) double well potentials between the next-nearest-neighbors (the red springs in the diagonal directions in figure 1 ; and (3) fourth-order interactions (springs) between the third-nearest-neighbors, = − − U r r ( ) ( 2) 4 (long green lines), where r is the distance vector (figure 1). The double well potentials between the next-nearest-neighbors were designed using inspiration from The interatomic potentials involve nearest-neighbor (black springs), next-nearest-neighbor (red springs) and third-nearest-neighbor (green dashed arrows) interactions. The black and green springs are harmonic and set the length scale of the interfaces, while the red springs represent the double well potentials (Landau springs) described in [42][43][44].
the Landau potentials to form a 4° shear angle. The shear angles of ferroelastic oxide materials are typically below 4°, while many metallic martensites have larger shear angles. We constructed the model so that the shear angle was fixed to 4°, which is a good compromise for metallic as well as oxide materials. The model potential mimics the ferroelastic phase transition in SrTiO 3 [45]. For three-dimensional models, we add potentials between adjacent layers: the interaction between the first-nearest-neighbors is the same as that in the plane, with harmonic springs adopted for the nextnearest-neighbors, = − U r r ( ) 10( 2 ) 2 . These interactions lead to a monoclinic structure that allows strain-compatible twinning.
We use this potential to simulate the shear instability in two and three dimensions with more than 1 million particles. Compared with our previous works [40,41], we increased the mass of particles from 1 to 1000 to rescale the simulation temperature to approach the phase transition temperature, T C [46]. The initial configuration has a two-dimensional sandwich configuration with two pre-existing horizontal domain boundaries. We then anneal each configuration at a given temperature for 5 × 10 5 time steps. After this relaxation, two buffer layers at the top and bottom of the two-dimensional sheet are sheared by the prescribed shear strain (hard boundary conditions). The shear is performed over 2 × 10 7 time steps, which stabilizes one domain orientation (-4°) and destabilizes the other (+4°). Free boundary conditions are adopted to allow the needle domains to nucleate from the free surface. The LAMMPS computer   code was used with an NVT ensemble, and the temperature of the sample was held constant by the Nosé-Hoover thermostat.

Results
The domain structures after molecular dynamic relaxations were analyzed at temperatures in two and three dimensions. For a three-dimensional model, the sample of 1.3 million particles contained 32 layers perpendicular to the plane, in which strain-compatible twinning is allowed by symmetry. The macroscopic shear angle and the averaged shear angles (proportional to the spontaneous strain) show the same decrease with increasing temperature (their squared values decrease linearly with T in figure 2(a)). The extrapolated transition temperature is T C = 130 K for the three-dimensional model. To evaluate the effect of the thickness of the sample, we reduced the sample to one layer. The temperature evolution of the strain is very similar to the thick sample, and the transition temperature, T C , is around 160 K.    Time renormalization relates all the time scales to the phonon frequency of our model. The phonon frequencies were measured using the dissipation-fluctuation theorem by the Fourier transform of the displacement patterns in two and three dimensions ( figure 2(b)). Phonon resonances appear at 0.077 ps -1 , so the phonon time is 13 ps. The phonon time does not depend significantly on the dimension of the sample. As the results in two and three dimensions are closely related, we make use of the higher number of particles in the plane of the twin boundaries and focus on two-dimensional simulations.
The evolution of the potential energy (Pe) with strain is shown for constant strain rate in figure 3. We measure the energy of the jerks in units of = J e e ( ) (dPe/d ) 2 . The derivative of the potential energy Pe with respect to strain e, ( = v e dPe/d ), reflects the velocity of the energy change. The energy of the jerks can be visualized as the kinetic energy (~v 2 ) of the dynamic process [31,47]. The energy spectra are evaluated with high time resolution from the potential energy of the system measured at every simulation step. In the spectrum, jerks start at strain e i when the J(e) signal crosses a fixed threshold and finishes at e j , when the signal remains below the threshold (red line in figure 4). The thresholds are set to be J 0.1 , where J is the average over all jerk energies in any given spectrum. The integral of the normalized energy signal (J J / ) for the duration ( ~ e e i j ) of the event describes the jerk energy E dissipated by jerks. The probability functions P(E) are the numbers of jerks in an energy interval between E and + E E d . The waiting time between adjacent jerks, using the same threshold, is calculated by converting + e e i j 1 into units of time. We analyze the jerk distribution by amplitude filtering. We define the energy windows and count the number of jerk peaks inside each window ( figure 4). The distribution contains cut-off energies as upper and lower window energies. In power-law distributions, the relative change in the probability ε = P E P E E E d ( )/ ( ) -d / is proportional to the relative change in energy and does not depend on any intrinsic energy scale. For thermally active processes, the probability is ~ P E E k T T ( ) exp-( / ( -)) VF , where T VF denotes the threshold of the thermal activation. The different statistics are clearly seen in figure 5, where the intervals in the plastic regime at  The raw data of the jerk spectra are shown for two temperatures, and two different structural regimes are collected in figure 5. The microstructural regimes are defined in figure 3, where the change in the potential energy of the sample is shown as a function of the applied shear strain. The initial elastic regime (I) leads to a regime with some nucleation of kinks in the interfaces between the sandwich layers (II), and finally to the yield regime (III) where the main pattern formation takes place. After yield, de-twinning occurs first for unstable needle domains, and then de-twinning proceeds essentially via the propagation of kinks in the twin walls (IV) until a single domain sample is recovered. Snapshots of the structural states in the various regimes are shown in figures 6 and 7.
The jerks in the various regimes differ in amplitude and profile. In figure 8, we compare the jerks in the elastic, the yield, and the plastic regime. Figure 8(a) shows that the amplitude of jerks is highest in the yield regime (black). The weak jerks are found in the elastic and plastic regime, where virtually all excitations are phononic. A surprising difference between the yield regime and the other regimes is that the individual jerks are smooth functions of the applied strain outside the yield regime. During yield, all the jerks are rugged with 'sub-jerks' inside them. This observation places emphasis on the scale invariance of jerks during yield: jerks represent avalanches, whereby each avalanche contains sub-avalanches that follow each other without the main avalanche coming to rest. This effect is even more obvious in the case of low-temperature avalanches in the yield regime. In figure 9, we show that the jerk distribution at low temperatures (0.5 and 1 K) has a clear onset for small strains and remains active for very large strain intervals. This clearly demonstrates that the individual jerks are part of a large, collective avalanche and not a collection of uncorrelated jerks. Such collective jerks also occur in the early elastic regime, which contains the low-temperature tweed dynamics. The jerks become smooth and uncorrelated at higher temperatures, when the elastic regime is dominated by both tweed dynamics and some front propagation of the twin boundaries. Such individual, uncorrelated jerks dominate at high temperatures at all conditions outside the yield regime.
The integrated distributions for the given regimes are summarized in figure 10. Jerks related to the extended avalanches are always power-law distributed, with an energy exponent ε =1.30 at T=10 K and ε =1.7 at T=1 K. The waiting time correlation follows a power law with an exponent of 2.0 at 1 K and a larger exponent at 10 K. The yield behavior is hence essentially athermal and scale-invariant, but shows different inter-jerk correlations with power-law distributed waiting times between the jerks.
Increasing the temperature from 1 to 10 K does not fundamentally change the avalanche mechanism (although the energy exponent changes slightly). The main difference is that the large avalanche at 1 K is broken up at 10 K into smaller segments. The increase in the energy exponent at higher temperatures appears empirically to be related to this breakup. The waiting time distributions show strong correlations between the individual jerks, similar to the results for the collapse of porous materials and earthquakes. The exponent of the waiting time power-law distribution is 2.4.
The effect of temperature is strong outside the yield regime. Thermal excitations, and hence the appearance of an energy scale, are very pronounced in the elastic and the plastic regimes, where we find that exponential distributions apply for jerks and for the waiting times between jerks. The only power-law distribution occurs just before the yield regime with a characteristically larger energy exponent and a strong exponential tail at high energies. At low temperatures, we find power-law distributions for jerks in all regimes besides the near-yield regime, where stretched exponentials characterize the crossover between the power law and the exponential regime.

Discussion
The discovery of a temperature-generated crossover between thermal and athermal regimes [40] is confirmed by our simulations. The origin of the crossover is related to the strain regimes outside the yield regime. These elastic and plastic regimes are thermally activated at high temperatures, while such activations are frozen out at sufficiently low temperatures. The spectra of the jerks show quasi-continuous profiles for each jerk, which relates to the continuous movement of kinks and needle domains under the applied strain field. These movements are thermally activated at sufficiently high temperatures. At low temperatures, we find that domain boundary freezing is in close agreement with the experimental observations [48,49]. The Vogel-Fulcher temperature is similar to that observed during a domain freezing process.
Domain freezing is much less important during yield events with energies much larger than the thermal energy. Domain nucleation during yield is scale-invariant and not thermally activated. The initial nucleation is a runaway process at low temperatures (0.5 and 1 K), so it only exhausts itself when the final number of domains has been produced. This process is hence geometrically one big avalanche, while all the other events are jerks with very weak correlations. It is therefore crucial to distinguish between jerks related to avalanches (during the yield event) and all other jerks, which can be thermally induced and do not lead to extended avalanches.
We now analyze the jerk profiles to identify a way to distinguish between the two types of jerk. Averaged jerk profiles were defined by summing over the interpolated profiles where the number of data points per jerk varies dramatically. We interpolated linearly between the data points within each jerk profile. The normalized profiles for the various regimes and temperatures are shown in figure 11. Subtle changes in the profiles show parabolic distributions in the yield regime (large-scale avalanches), while Gaussian distributions occur for cases where the excitation is mainly related to phonons in the plastic regime.
The temperature dependence of the yield effect causes a slight change in the energy exponent, but also has a very strong effect on the resulting microstructure (figure 7) . Sparse microstructures are found at low temperatures (front propagation) and high temperatures (thermal excitations), while the densest microstructure occurs at 1 and 10 K. The energy exponent changes accordingly from 1.88 (0.5 K), 1.66 (1 K) to mean field values near 1.35 at higher temperatures.