Introduction

Serrated stress-strain curves have been observed in refractory high-entropy alloy (RHEA)1,2,3. The material is a recently developed alloy family with promising mechanical performance at high temperatures4. Possible applications include jet engines5 and nuclear reactors6. Relatively little is known about the slip statistics of serrations, whose uncertainty undermines the material’s reliability and hinders further applications of RHEA. Here we model their dependence on experimental tuning parameters including temperature and strain rate and compare the results to experiments.

Our simple mean-field theory (MFT)7,8 predicts the scaling behavior of the underlying slip statistics. The large slip avalanches with system-spanning runaways can be attributed to a dynamical weakening effect in the model, which is linked to the solute-dislocation interactions in RHEA. We present the dependence of this weakening on the experimental tuning parameters, such as temperature and strain rate. The new understanding is useful for controlling the solute strengthening mechanisms in RHEAs and optimizing the material’s mechanical properties for high-temperature applications.

A stress drop in the stress-strain curve corresponds to a slip avalanche resulting from the collective movement of dislocations9. Slip avalanches in HEAs have been extensively studied, as discussed in review papers10,11. Recent studies have also explored the avalanche behavior of HEAs at elevated temperatures12,13,14,15. Furthermore, several studies have employed atomistic simulations to determine the universal mechanisms behind slip avalanches and serrated flows in HEAs16,17.

While most past studies have focused on CoCrFeMnNi HEAs or related materials, studies on RHEAs, including HfNbTaTiZr HEA, remain relatively rare. It is still unclear how deformation behaviors at high temperatures differ between these two types of materials. Fortunately, more research has been conducted on RHEAs, as summarized by Bamisaye et al.18. Our study on the universal slip statistics aims to shed light on the underlying mechanisms behind serrated flows in RHEAs. We believe our work will contribute to the development of a unified theory of slip avalanches, not only in RHEAs but also in other related materials.

Mean-field theory of slip avalanches

A simple slip theory that goes beyond depinning by adding weakening to the local thresholds8 provides predictions for the various aspects of slip statistics including size distribution and scaling laws, and temporal avalanche profiles that are also called “shapes”. Without any weakening, the model predictions correspond to the universality class of the mean field depinning phase transition in far-from-equilibrium disordered systems19. The weakening adds another relevant parameter to the underlying critical point with its associated critical exponents. Also, it changes the nature of the depinning transition; while it is continuous at zero weakening, it becomes a first order transition at nonzero weakening. The discontinuity of the first order transition can be seen as macroscopically large slip avalanches with system-spanning runaways.

In the MFT, we regard dislocations as weak spots in the solid. We consider a system that is composed of \(\:N\) sites that can slip when a dislocation passes by. Each site (labeled \(i = 1, \ldots ,N\)) is subject to an accumulated stress, \(\:{\tau}_{i}=\frac{J}{N-1}\sum_{j=1}^{N}\left({u}_{j}-{u}_{i}\right)+K({v}_{app}t-{u}_{i})\), where \(\:{u}_{i}\) is the accumulated displacement, \(\:J\) is the mean-field elastic coupling coefficient, \(\:K\) is the stiffness parameter and \(\:{v}_{app}\) is the applied displacement rate at the sample boundary. At first, the external loading adiabatically increases with time \(\:t\) until a weak spot fails. The failure occurs when its stress \(\:{\tau}_{i}\ge\:{\tau}_{s}\) where \(\:{\tau}_{s}\) is the failure threshold. The failure induces a slip that increases \(\:{u}_{i}\) by a small amount, \(\:{\Delta}{u}_{i}=\frac{{\tau}_{s}-{\tau}_{a}}{J+K}\) which effectively reduces the stress from \(\:{\tau}_{s}\) to an arrest stress \(\:{\tau}_{a}\). In the meantime, the dynamic weakening induced by the failure temporarily lowers the local threshold of each the slipping site until the end of the avalanche, \(\:{\tau}_{s}\to{\tau}_{d}={\tau}_{s}-\epsilon(\tau_s-\tau_a)\), where \(\:\epsilon>0\) is the weakening parameter. Furthermore, the failure also destabilizes the other weak spots because the released stress is redistributed uniformly over the entire system, \(\:{\tau}_{j}\to\:{\tau}_{j}\:+\frac{J}{J+K}\frac{{\tau}_{s}-{\tau}_{a}}{N-1}\), via the mean-field coupling. A chain reaction like this will then generate a cascade of failure events that resembles an avalanche. By iterating both the loading and failure steps, we can simulate the avalanches and compute the corresponding statistics.

Despite its simplicity, the slip model has successfully predicted the scaling properties of slip statistics in various deformed solids on different length/time scales, including micro- and nano-pillars of crystalline metals20,21, bulk metallic glasses22,23,24,25,26, high-entropy alloys27,28 and granular materials26,29,30. The reason for the enormous size of the underlying universality class is that the scaling properties of a universality class are independent of the microscopic details. Note that all slip experiments with which the model agrees, display slip localization during avalanches in shear bands or glide planes. For such slip localization the assumptions of the slip model (especially that the slip interactions are positive) are applicable. (If the slip is not localized in any way, then other universality classes might become relevant31,32.) In this paper, we focus on situations where avalanches occur in localized bands or shear planes, so that the MFT with positive couplings is applicable.

Our approach using MFT aims to capture the essential physics of crystal plasticity by focusing on the simplest possible model. While it is true that dislocations can act as both weakening and strengthening sites, our model treats mobile dislocations as weak spots to identify the fundamental mechanisms driving slip behavior. The simplicity of our model is intentional, allowing us to develop intuition and identify the key physical ingredients. Adding more complexity could obscure these basic mechanisms. Our goal is to capture universal aspects of dislocation dynamics, similar to how the Ising model captures fundamental interactions in magnetism. While detailed engineering models include all specific interactions, our approach seeks to highlight the underlying symmetries and interactions applicable across a wide range of systems, providing a general framework for understanding the physics of plasticity.

Besides MFT, the serration statistics associated with the Portevin-Le Chatelier (PLC) effect have been extensively studied, particularly for Al-Mg alloys. Power laws have been observed in stress-drop distributions33,34,35 and in acoustic emission energy distributions36,37. Mäkinen et al. found that the velocity statistics of propagating type-A PLC bands align with a mean-field model of critical avalanche dynamics38. The dynamic strain aging effect explains the temperature and strain rate dependence of the PLC effect39,40,41. This theory attributes negative strain-rate sensitivity to the interaction between diffusing solute atoms and mobile dislocations. Negative strain rate sensitivity is similar to the dynamic weakening effect in MFT, which involves temporary reductions in local threshold. The root causes of both effects are changes in the threshold to slipping for the microscopic weak spots (i.e. the dislocations).

In contrast to the dynamic strain aging theory, the dynamic weakening effect was developed to explain the dynamic phase transitions in many different systems. Using tools from statistical physics we have shown that the results of the dynamic weakening can be applied to systems far beyond just HEAs, to amorphous materials, like bulk metallic glasses25,26,42,43, even to earthquakes44,45,46, considering that seismic faults re-heal their strength during quiescent periods. Here we are testing how far universal aspects of the theory originally developed for earthquakes applies to experiments on HEAs. Although phenomenologically similar, these theories address different phenomena. In our work, we adopt the MFT with the weakening effect, originally developed for earthquakes, to study the universal slip statistics of RHEA, i.e. this paper is linking different structures, different scales, and different fields.

Detrended data and three phases of slip behavior

The experimental serration data used here has previously been published in the paper3. Here we briefly outlined the experimental procedure. Alloy samples of HfNbTaTiZr RHEA were prepared using vacuum arc melting from high-purity elements. These samples underwent three re-melting cycles to ensure uniform composition. The resulting ingots, measuring 75 mm in length, 25 mm in width, and 10 mm in thickness, were cold rolled to create alloy sheets with a 70% reduction in thickness. Subsequently, they were annealed at 950 °C for 1 hour and then rapidly cooled with water quenching. Tensile specimens were cut from these annealed sheets using electrical discharge machining, with dimensions of 19 mm length, 3.0 mm width, and 2 mm thickness. High-temperature tensile tests were conducted using a SHIMADZU AGS-X 100kN universal tester and an ATS series-3210 tube furnace, applying a strain-controlled protocol. Experiments were repeated at least twice for accuracy.

High-temperature tensile tests were performed on HfNbTaTiZr RHEA across different temperatures \(\:T\) and strain rates \(\dot{\gamma }\). For example, we show in Fig. 1a the detrended stress-time curves for different temperatures \(\:\text{T}=473\text{K}\), \(673\;{\text{K}}\) and \(\:873\;{\text{K}}\) at a fixed strain rate \(\dot{\gamma } = 10^{{ - 3}} \,{\text{s}}^{{ - 1}}\). The detrends are obtained by subtracting the moving average from the original data to emphasize the serrated flow that fluctuates around the slow-varying trend. Here, we choose the bin width \(\:{\Delta\:}{t}_{\text{b}\text{i}\text{n}}=2\:\text{s}\text{e}\text{c}\) for the moving average. Note that the resulting serration statistics are insensitive to the bin width if \(\:{\Delta\:}{t}_{\text{b}\text{i}\text{n}}\) is much greater than the serration’s duration, which is less than 1 second in this case. An additional Wiener filter is applied afterward to mitigate the influence of the background noise. The details of the Wiener filtering technique are shown in the Supplemental Information.

We measured slip statistics only during the stationary state when work hardening is negligible. This ensures our data reflects the intrinsic dynamics of weakening and slip processes without the influence of ongoing work hardening, which has been discussed by Hähner and Rizzi40,41. Our rapid weakening effect, driven by dislocation escape and re-pinning, is distinct from the gradual work hardening process. By focusing on the stationary state, we capture the intrinsic behavior of the system’s slip distributions.

We categorize the slip behavior of serrations shown in Fig. 1a into three distinct phases (from low to high temperature): slow-avalanche phase, fast-runaway phase, and quiescent phase. For the quiescent phase at \(\:873\;{\text{K}}\), only small stress drops are observed after detrending the data. In contrast, after detrending both the slow-avalanche phase at \(\:473\;{\text{K}}\) and the fast-runaway phase at \(\:673\;{\text{K}}\) significant serrations are found. To better quantify the observed slip behavior, we examine the statistics of stress drops including size, duration, and maximum velocity. The result is summarized in table 1. Note that the (slip) velocity is defined as the time derivative of stress, which is proportional to the time derivative of the total slip integrated over the slipping area in the RHEA.

Fig. 1
figure 1

(a) Detrended stress-time curves for different temperatures \(\:T=473\;{\text{K}}\)(blue), \(\:673\;{\text{K}}\)(orange) and \(\:873\;{\text{K}}\)(green). The strain rate is fixed at \(\dot{\gamma } = 10^{{ - 3}} \,{\text{s}}^{{ - 1}}\). We categorize the slip behavior into three distinct phases: slow-avalanche phase, fast-runaway phase, and quiescent phase. The detrends are shifted vertically to avoid overlap. (b) The complementary cumulative distribution function (CCDF) of avalanche sizes \(S\). The dashed line indicates the power law with a cutoff. \(\:\text{CCDF}\left(S\right)\sim{\int}_{S}^{\infty}{{S}^{{\prime}}}^{-3/2}{e}^{-\frac{{S}^{{\prime}}}{{S}_{c}}}\:dS{^\prime}\), predicted by the mean-field theory (MFT). The gray shaded area shows the contribution from background noise. (c) Average duration \(\:D\) and (d) Average maximum velocity \(\:{V}_{m}\) against the average size \(\:S\). The dashed lines indicate the scaling laws, \(\:D\sim{S}^{0.5}\) and \(\:{V}_{m}\sim{S}^{0.5}\), predicted by the MFT for small avalanches that do not span the entire system. The dotted lines show the model predictions for the scaling laws for large (i.e., system-spanning) runaways, \(\:D\sim{S}^{0}\) and \(\:{V}_{m}\sim{S}^{1}\).

Table 1 Slip statistics of stress drops for three different phases.

Size distribution and scaling laws

First, let’s look at the size’s statistical distribution. We plot in Fig. 1b the complementary cumulative distribution function (CCDF) of sizes for each phase. With only small slip avalanches, the quiescent phase at \(\:873\;{\text{K}}\) has a probability distribution function (PDF) of avalanche sizes that decays faster than a power law. The fast decay results from the small cutoff Sc, which is almost indistinguishable from the extent of the background noise as shown in gray shaded area. In the following, with “apparent” slip avalanche we refer to avalanches above the noise floor, i.e., avalanches whose size \(\:S\) is much larger than the background noise\(\:\:{S}_{\text{n}\text{o}\text{i}\text{s}\text{e}}\sim0.04\:\text{M}\text{P}\text{a}\). The noise time trace reflects the background stress fluctuations. It was obtained by running a tensile testing experiment without loading any specimen. From this time trace we extracted the corresponding “noise” stress-drop size distribution.

The excess number of large stress drops is visible as a bump in the PDFs or a steep drop-off in the CCDFs of both, the slow-avalanche phase at \(\:473\;{\text{K}} ,\) and the fast-runaway phase at \(\:673\;{\text{K}}\). And for both cases, the stress drops seem to be consistent with the MFT-predicted distribution, that is, \(\:\text{C}\text{C}\text{D}\text{F}\left(S\right)\sim{\int}_{S}^{\infty}{{S}^{{\prime}}}^{-\frac{3}{2}}{e}^{-\frac{{S}^{{\prime}}}{{S}_{c}}}d{S}^{{\prime}}\) where \(\:{S}_{c}\) is a size cutoff. According to the MFT, the exact value of the cutoff \(\:{S}_{c}\) is nonuniversal and material-specific, but its scaling with experimental parameters is universal and independent of the microscopic details. For example, the cutoff is predicted to scale with the stiffness of the machine as \(\:{S}_{c}\sim{K}^{-2}\). In the Supplemental Information, we show fits of different distribution functions to the data, including the MFT-predicted power law with a cutoff, a purely exponential distribution, and a Gaussian distribution. The result shows that the data agrees with the MFT although we cannot rule out the possibility of other distributions. Having more data in the future to decrease the size of the statistical errorbars might allow us to narrow the number of possible distributions.

Besides the size distribution, we also notice that the slope of stress drops at \(\:673\;{\text{K}}\) are much steeper than those of \(\:473\;{\text{K}}\), even though their overall sizes are comparable. To visualize this qualitative distinction, we plot the duration \(\:D\) and maximum velocity \(\:{V}_{m}\) against size \(\:S\) in Fig. 1c and d to highlight the significant difference in scaling relations. Note that each dot in Fig. 1c and d represents the average \(\:D\) or \(\:{V}_{m}\) over the avalanches within a logarithmic size bin, and the error bar denotes the 2.5%- and the 97.5%-percentiles of \(\:D\) or \(\:{V}_{m}\). As we expect, the slow-avalanche phase at \(\:473\;{\text{K}}\) has the slip statistics exhibiting the predicted scaling laws, \(\:D\sim{S}^{0.5}\) and \(\:{V}_{m}\sim{S}^{0.5}\), which is consistent with the avalanche statistics derived from the MFT for the avalanches in the power-law scaling regime.

The fast-runaway phase at \(\:673\;{\text{K}}\) has different scaling laws, that is, \(\:D\sim{S}^{0}\) and \(\:{V}_{m}\sim{S}^{1}\) which notably deviate from those of the slow-avalanche phase. According to the MFT, this deviation reflects the runaway nature of the avalanches, which is caused by the dynamic weakening during plastic deformation8. In fact, the weakening effect has been already observed in study of CoCrFeMnNi HEAs28,47 and has been attributed to the solute-dislocation interaction in the alloy. Note that these scaling relations have previously also been used to identify the runaway slip avalanches in compressed bulk metallic glasses and sheared granular materials26,42.

The scaling laws allow us to rule out the possibility that large slips are merely due to finite size effects or overlapping effects. The former occurs when the finite system size limits the slip avalanches, while the latter results from the merging of temporally overlapping small events into synthetic large events. Their scaling laws, (with \(\:D\sim{S}^{0.5}\) and \(\:{V}_{m}\sim{S}^{0.5}\) for finite size effects, and\(\:\:D\sim{S}^{1}\) and \(\:{V}_{m}\sim{S}^{0}\) for overlapping effects) are qualitatively different from those of runaways characterized by \(\:D\sim{S}^{0}\) and \(\:{V}_{m}\sim{S}^{1}\) as seen here.

So far, the empirical avalanche statistics match the prediction by the MFT. Yet, we should point out that there is a small discrepancy between the experimental data and the MFT: the runaway avalanches in empirical data have a relative broad distribution in size, which deviates from the MFT prediction that the runaways should have a characteristic size with small variation. This discrepancy is interesting but not unexpected because the characteristic size of runaways is expected to depend on the system size, which is a material-specific detail and could vary with time as the dislocation structure evolves during experiments for example.

Note that several studies have investigated slow and runaway avalanches in various dislocation models, as reviewed by Papanikolaou et al.48. However, none has explicitly related these two regimes within a unified framework. For instance, the dynamic aging effect has been considered in continuum models40,41,49,50, where dislocation motion is slowed when arrested by a solute atmosphere. Additionally, forest dislocations can also impede dislocation motion by acting as obstacles in slip planes51,52. On the other hand, runaway-like dislocation avalanches have been observed in nanopillar compression tests53 and attributed to the breakaway of dislocation pileups54,55. In our study, we model the transition between these two avalanche regimes using a MFT, providing an intuitive and quantitative understanding of the weakening effect that drives this non-equilibrium phase transition in RHEAs.

Scaling collapse of avalanche shapes and runaway avalanches

To highlight the difference between the slow-avalanche phase and fast-runaway phase, we show the size-dependent avalanche shapes in Fig. 2 for \(\:473\;{\text{K}}\) and \(\:673\;{\text{K}}\). Each curve corresponds to the average shape \(\left\langle {v\left( t \right)|S} \right\rangle\) of avalanches within a given size bin \(\:S\). And the shaded regions show the 95%-confidence interval evaluated with the bootstrap method56. (Here, we resample data with replacement within a given size bin and then estimate the corresponding confidence interval of the average shape.) To demonstrate the avalanche’s scale invariance, we have rescaled the shapes rescaled by \(\:{S}^{1/2}\). The inset plots show the original shapes before the rescaling. Predicted by the MFT, the rescaled shapes of the small avalanches in Fig. 2a collapse into the predicted universal curve, \(\left\langle {v\left( t \right){\text{|}}S} \right\rangle \sim Ate^{{ - \frac{{At^{2} }}{{2S}}}}\) (the dashed line). Note that \(\:A\simeq70\:\text{M}\text{P}\text{a}/{\text{s}}^{2}\) is a non-universal scaling factor.

Fig. 2
figure 2

Temporal avalanche profiles (“shapes”) for small avalanches of given size \(\:S\) rescaled by \(\:{S}^{1/2}\) for (a) the slow-avalanche phase at 473 K and (b) the fast-runaway phase at 673 K. Insets show the original shapes before the rescaling. Each curve corresponds to the average shape \(\left\langle {v\left( t \right)|S} \right\rangle\) of avalanches within a given size bin \(\:S\). The shaded regions indicate the 95%-confidence interval evaluated with the bootstrap method. Predicted by the mean-field theory, the rescaled shapes in (a) collapse into a universal curve, \(\left\langle {v\left( t \right){\text{|}}S} \right\rangle \sim Ate^{{ - \frac{{At^{2} }}{{2S}}}}\), indicated by the dashed line. In contrast, the sharp peaks in (b) reflect the rapid nucleation of runaways, whose shapes are predicted not to collapse in the same way as those of the small avalanches. This model prediction is also confirmed here by the experimental data.

Deviating from the universal curve, the emergence of sharp peaks in Fig. 2b reflects the nucleation of runaway and its rapid expansion. The mechanism behind runaways is dynamic weakening which temporarily lowers the material’s strength. At a microscopic scale, the reduction of material’s strength is caused by the escaping of depinning dislocations from solute segregation, which originally strengthens the material by pinning the dislocations. This weakening will effectively boost the avalanche expansion by facilitating the subsequent slip events. It forms a positive feedback loop because the subsequent slip events weaken the material’s strength even further when more dislocations have escaped from the solute segregation. As a result, an avalanche with sufficient large size will become unstoppable if its cumulated weakening effect can surpass the dissipation. The avalanche ends up being a runaway quickly spanning over the entire system.

In summary, our results show two distinct slip behaviors in slowly deformed RHEAs: the slow-avalanche and fast-runaway phases. In Fig. 1b, we illustrate that these two phases produce large slip avalanches significantly exceeding the background noise. Figure 1c and d demonstrate that the scaling relations between size, duration, and maximum velocity align with MFT predictions. Specifically, the slow-avalanche phase features duration and maximum velocity scaling with \(\:S^{1/2}\), while the fast-runaway phase shows duration saturation and maximum velocity scaling with \(\:S\).

Furthermore, Fig. 2 highlights that avalanches in the slow-avalanche phase can be collapsed in a scaling collapse, where the rescaled temporal avalanche shapes collapse onto universal scaling function predicted by MFT. The amount by which the individual shapes have to be rescaled is also predicted by MFT and consistent with the data. In contrast, the system-spanning events in the fast-runaway phase display a Gaussian peak, which strongly deviates from the universal curve. This Gaussian shape is also predicted by the MFT to indicate the rapid nucleation of runaway avalanches. Scaling collapses constitute strong, state of the art tests for model predictions, because they test not only scaling exponents, but entire scaling functions. The presented observations thus support the applicability of MFT to our data, while also highlighting the unique characteristics of the fast-runaway phase.

We demonstrate that the empirical data aligns with the power-law size distribution predicted by MFT. However, it’s essential to note that solely observing a power-law size distribution does not definitively confirm consistency with MFT. Therefore, the paper delves much deeper into various slip statistics, encompassing statistical properties of distributions of size, duration, velocity, and temporal avalanche shapes. The scaling collapses are a much stronger test than just trying to fit power laws to distributions. The strength of the statistical tests results from the combination of scaling collapses with distribution tests, and getting a picture that consistently fits the MFT for all six statistical measures that we test in this paper. By thoroughly analyzing these statistics, we can assert with greater confidence that the slip avalanches of RHEA indeed are consistent with the MFT predictions.

Phase diagram of slip behavior

Now that the connection between runaway avalanches and dynamic weakening is established, it will be interesting to see how the slip behaviors of the RHEA depend on experimental conditions. As a summary, we present the constructed nonequilibrium phase diagram in Fig. 3 to summarize the crucial roles of solute-dislocation interaction in the slip behaviors of RHEA at different temperatures \(\:T\) and strain rates \(\dot{\gamma }\). The result is also listed in Table 2. Notice that the \(\:x\)-axis in Fig. 3 is the reciprocal temperature, \(\:1000/T\).

Fig. 3
figure 3

Nonequilibrium phase diagram of slip behaviors for different temperatures \(\:T\) and strain rates \(\dot{\gamma }\). Note that the \(x\)-axis is the reciprocal temperature, \(\:1000/T\). The black dots, the red crosses and the blue squares denote the quiescent phase, the fast-runaway phase, and the slow-avalanche phase, respectively. Different color of contours in the background indicates the ratio, \(\:\frac{\varDelta{t}_\text{Slip}}{\varDelta{t}_\text{Diffusion}}\sim\frac{\rho{D}_{T}}{\dot{\gamma\:}}\) where \(\rho \simeq 3 \times 10^{{15}} \:{\text{m}}^{{ - 2}}\) is the empirical dislocation density and \(\:{D}_{T}={D}_{0}{e}^{-\frac{Q}{RT}}\) is the temperature-dependent diffusion coefficient with \(D_{0} \simeq 3.77 \times 10^{8} \:{\text{m}}^{2} /{\text{s}}\) and \(Q \simeq 42.1\:{\text{kcal}}/{\text{mol}}\). The right-hand side is the ratio of the typical slip time scale, \(\:\varDelta{t}_\text{slip}\sim\frac{{b}^{2}\rho}{\dot{\gamma}}\), to the typical diffusion time scale \(\:\varDelta{t}_\text{Diffusion}\sim\frac{{b}^{2}}{{D}_{T}}\), where \(\:b\) is the magnitude of a Burger’s vector. Reddish color shows the diffusion-dominated regime and bluish color shows the slip-dominated regime. The slow-avalanche phase and fast-runaway phase are in an intermediate area when two time scales are roughly equal to each other, \(\:\varDelta{t}_\text{Slip}\approx\varDelta{t}_\text{Diffusion}\).

As we expect, the fast-runaway phase (red crosses) appears at a higher temperature than the slow-avalanche phase (blue squares), when high temperature accelerates solute diffusion, facilitates solute segregation, and thus enhances the associated weakening effect. For either high temperature \(\:T\ge\:873\;{\text{K}}\) and/or low strain rate \(\:\dot{\gamma\:}\le\:{10}^{-4}\:\text{s}^{-1}\), the quiescent phase (black dots) has no serrations observed in the detrended data. It implies that the serrated flow is effectively suppressed by thermally assisted creep, which has already been shown by a recent study with the depinning theory in finite dimensions57.

We also compare our phase diagram of slip behaviors with the traditional phase diagram of the PLC effect, which has been categorized into three different types (type-A, type-B, and type-C)3,58. Typically observed during a tensile test, type-A bands exhibit small, temporally irregular stress drops during continuous propagation across a tensile specimen. In contrast, type-B bands propagate intermittently at roughly equal intervals in space, with larger, irregular stress drops, while type-C bands display larger, temporally periodic stress drops and occur at spatially random nucleation sites throughout a specimen. By comparing these two phase diagrams, we found that both the slow-avalanche and fast-runaway phases correspond to the type-A PLC effect. However, the fast-runaway phase is near the dynamic phase transition from type-A to type-B or type-C. We suspect that the fast-runaway phase may be related to strain localization occurring in space when a continuously propagating type-A band transitions to randomly nucleating type-B or type-C bands.

This correspondence also helps explain why, in Fig. 1, the \(\:{S}^{1/2}\)-scaling of \(\:D\) and \(\:{V}_{m}\) at \(\:473\;{\text{K}}\) extends to sizes that exceed (by one decade) those observed at \(\:673 \;{\text{K}}\). This difference arises from the difference in effective slip-system size facilitating the slips for the respective type of PLC bands. At \(\:473 \;{\text{K}}\), the slow-avalanche phase has a larger effective system size because of propagating type-A PLC bands, which encounter more dislocations as they travel through the material, making the system size equal to the total number of dislocations in the specimen. In contrast, at \(\:673 \;{\text{K}}\), the fast-runaway phase is characterized by localized type-B or type-C bands, limiting the effective system size for the dynamics to the dislocations within the localized region.

In MFT, the system size refers to the total number of weak spots accessible during avalanches, depending on the dislocation density and whether the PLC band is propagating or localized. For non-propagating type-B or type-C bands, the system size is confined to the localized band. Conversely, propagating type-A bands increase the system size as they move and encounter more dislocations. This is supported by a study on Al alloys38, which showed that avalanche sizes follow a power-law distribution with a cutoff equal to the specimen length. Thus, the difference in effective system sizes for the dynamics between \(\:473 \;{\text{K}}\) and \(\:673 \;{\text{K}}\) explains the observed scaling behavior in Fig. 1.

Upon comparing our phase diagram with traditional PLC phase diagrams in Table 2, we observed that both the slow-avalanche and fast-runaway phases correspond to the type-A band. However, the fast-runaway phase is closer to the boundary between type-A and type-B or C bands. This suggests that the runaway phenomenon and associated weakening effect may significantly influence the dynamical phase transition of PLC effects. Despite the limitation of a small dataset, our results distinctly illustrate the qualitative differences between these phases. Future experiments could complement our findings by adding more data points to the phase diagram. This would provide further insight into the relationship between microscopic slip dynamics and the macroscopic PLC phenomenon.

Notice that a reversed strain-rate dependence of PLC band types is observed at the highest tested temperature in Table 2. Specifically, at 873 K, the deformation bands transition from type-A to type-B and then to type-C as the strain rate increases. We attribute this inverse trend to thermal softening3 and the resulting overlapping effect of avalanches at high temperatures. Unlike conventional cases where type-A bands dominate at high strain rates, the high dislocation density in RHEAs promotes random strain localization, leading to type-B band formation at \(\:\dot{{\gamma}}={10}^{-3}\:\text{s}^{-1}\) and type-C band formation at \(\:\dot{{\gamma}}={10}^{-2}\:\text{s}^{-1}\). A detailed discussion and scaling analysis are provided in the Supplementary Information. Our findings highlight the critical role of avalanche overlapping in this process: as the strain rate increases, slip avalanches occur more frequently and begin to merge, forming larger events. This overlapping effect, confirmed by the scaling behavior of avalanche statistics, leads to the emergence of type-B and type-C bands at higher strain rates. These results not only clarify the deformation behavior of refractory HEAs but also highlight the differences in dominant deformation mechanisms between body-centered cubic and face-centered cubic HEAs at elevated temperatures.

Table 2 Slip behaviors and the corresponding Portevin-Le Chatelier (PLC) bands at different temperatures \(T\) and strain rates \(\dot\gamma\). The slip behaviors are categorized into three phases: the quiescent phase, slow-avalanche phase, and fast-runaway phase. The PLC bands are categorized into three different types: type-A, type-B and type-C PLC bands3. Note that blank entries mean no tests have been conducted under those conditions.

In summary, the phase diagram suggests that there exists an optimal temperature range (for fixed strain rate) where the weakening effect is significant, and runaway avalanches emerge as a result. Here, we attribute such “weakening regime” to the intricate balance between dislocation slip and atomic diffusion. If the dislocation slip is much faster than the atomic diffusion, the weakening effect would be suppressed because dislocations slip away from pinning sites before solute atmosphere forms. On the other hand, fast atomic diffusion and/or slow dislocation slip inhibit the weakening effect because the solute atmosphere could now migrate with slow moving dislocations.

The “weakening regime”, therefore, corresponds to an intermediate regime that dislocations stay at pinning sites long enough for solute atmosphere to form, yet slip frequently/fast enough to escape from the atmosphere. To justify the argument, we estimate the typical time scales of atomic diffusion and dislocation slip for the RHEA and compare their values at different temperature \(\:T\) and strain rate \(\:\dot{\gamma\:}\). Note that the similar time scale analysis has already been utilized in the recent study on CoCrFeMnNi HEA47.

According to the Orowan’s equation59, the dislocation-mediated plasticity has a strain rate that can be written as \(\:\dot{\gamma}=\rho b{v}_{D}\) where \(\:\rho\:\) is the dislocation density, \(\:b\:\)is the magnitude of a Burger’s vector and \(\:{v}_{D}\) is the average dislocation velocity. Therefore, the typical time scale for a dislocation to slip across one Burger’s vector is given by \(\:\Delta{t}_{\text{s}\text{l}\text{i}\text{p}}\sim\frac{b}{{v}_{D}}=\frac{{b}^{2}\rho\:}{\dot{\gamma}}\), where \(\:\rho\simeq3\times{10}^{15}\:{\text{m}}^{-2}\) for the RHEA60. On the other hand, the diffusion time scale can be estimated by considering the time required for an atom to diffuse to a neighboring lattice site, that is, \(\:\Delta{t}_{\text{D}\text{i}\text{f}\text{f}\text{u}\text{s}\text{i}\text{o}\text{n}}\sim\frac{{b}^{2}}{{D}_{T}}\) where \(\:{D}_{T}={D}_{0}{e}^{-\frac{Q}{RT}}\) is the temperature-dependent diffusion coefficient with \(D_{0} \simeq 3.77 \times 10^{{ - 8}} \:{\text{m}}^{2} /{\text{s}}\) and \(\:Q\simeq42.1\:\text{k}\text{c}\text{a}\text{l}/\text{m}\text{o}\text{l}\) for the RHEA61,62.

The ratio between these two time scales is then equal to

$$\frac{{\Delta t_{{{\text{Slip}}}} }}{{\Delta t_{{{\text{Diffusion}}}} }}\sim \frac{{\rho D_{0} e^{{ - \frac{Q}{{RT}}}} }}{{\dot{\gamma }}}$$

which is a function of \(\:T\) and \(\dot{\gamma }\). The contours in Fig. 3 represent different values of the ratio. Reddish color indicates the diffusion-dominated regime at high \(\:T\) and low \(\dot{\gamma }\), while bluish color indicates the slip-dominated regime at low \(\:T\) and high \(\dot{\gamma }\). As we expect, the fast-runaway phase is located at an intermediate regime where the time scales are comparable to each other. There is no runaway avalanche in the diffusion-dominated regime and the slip-dominated regime where the dynamic weakening is effectively suppressed.

Conclusion

Our result establishes the connection between slip behavior of serration flows and the dynamic weakening effect which is associated with the strong solute-dislocation interaction in HfNbTaTiZr RHEA. In comparison with the recent MFT study on CoCrFeMnNi HEA28 in which researchers focused solely on the size distribution of avalanches, here we examine many more aspects of slip statistics including not only the size distributions but also the scaling relations between size, duration and maximum velocity and the scaling collapse of rescaled avalanche shapes. The observed scaling exponents are all consistent with the critical values predicted by the MFT. Our comprehensive analysis reveals the qualitative differences in scaling laws between the slow-avalanche phase and the fast-runaway phase, which are sometimes hard to distinguish based on the size distribution alone. As a result, we conclude that the large serrations in the fast-runaway phase are indeed caused by the weakening effect in the HfNbTaTiZr RHEA. Furthermore, comparing the slip statistics with traditional PLC diagrams shows that the fast-runaway phase occurs near the type-A/B and type-A/C boundaries, suggesting that the weakening effect may strongly influence the dynamic phase transition of PLC effects. Future research on high-entropy materials can adapt the nonequilibrium phase diagram constructed in this paper to prevent undesirable serrations and thus improve the material’s reliability and controllability.