A Novel Approach to Analysis of Complex Crystallization Behavior in Zr-Based Bulk Metallic Glass by Non-isothermal Kinetics Studies

The complex crystallization behavior of the Zr40Ti15Cu10Ni10Be25 bulk metallic glass (BMG) produced by suction-casting method was studied with the non-isothermal DSC measurements with the heating rate from 5 to 40 K/min. Three exothermic phenomena were observed for the investigated material. The novel evaluation procedure for qualitative and quantitative analysis of intricate crystallization kinetics for Zr-based BMGs is proposed. The unusual deconvolution of the DSC curves based on a Gaussian function and a two-phase exponential decay function allowed for separate, detailed analysis of overlapped peaks. The activation energies for each crystallization stage were studied based on overall (Kissinger) and local (Kissinger–Akahira–Sunose) procedures. The KAS method applied separately for both low and high heating rates showed a significant difference in local activation energies. Finally, the local Avrami exponent evaluation revealed that the first two stages of crystallization are diffusion-controlled with mainly increasing nucleation rate, whereas the third crystallization is more growth-dominated.


I. INTRODUCTION
THE discovery that amorphous internal structures in metallic alloys can be formed through rapid cooling from the liquid state, first observed in 1960, [1] sparked a significant research effort into these materials. This led to the identification of various new compositions of multifunctional metallic glasses with high Glass Forming Ability (GFA) [2] and a better understanding of their physical and mechanical properties. [3][4][5][6][7][8][9][10][11] Additionally, new production methods were developed as a result of this research. [12,13] With the advancement of optimized chemical compositions and empirical rules proposed by Inoue, [14][15][16] it is now possible to produce Bulk Metallic Glasses (BMGs) [17] of considerable size, with diameters ranging from several millimeters to centimeters. [18,19] However, it is important to note that BMGs are not thermodynamically stable (they are metastable after production [20][21][22] ), and will undergo crystallization as their internal energy is greater than that of the crystalline phase. [2,23] The crystallization of glassy alloys takes place at elevated temperatures, when the energy barrier related to the motion of the atoms is overcome and the possibility of displacement of atoms in the alloy is unblocked. [2,22,23] A thorough understanding of the processes and kinetics of crystallization is crucial for both technological and scientific purposes. On the one hand, the ongoing crystallization is a limiting factor for the usage of BMGs at elevated temperatures. On the other hand, the presence of a Supercooled Liquid Region (SLR) between the glass transition temperature and the onset crystallization temperature, [24] in which mechanical and thermodynamic parameters are different from the glass state, [25,26] enables the possibility of the thermoplastic forming of BMGs. [27][28][29][30] In most cases, it is necessary to avoid crystallization during processing, which requires a detailed description of crystallization kinetics. From a scientific standpoint, the study of crystallization in BMGs offers valuable insights into nucleation and growth mechanisms in undercooled liquids.
BMGs composed of multiple elements with a zirconium base are a highly promising class of metallic alloys due to their exceptional properties. They are characterized by high glass forming ability, [9,31,32] enhanced mechanical properties, significant corrosion resistance, favorable biological properties, [33][34][35][36][37][38] and wide SLRs providing suitability for thermoplastic processing. [9,27,30] Despite the unflagging popularity, their crystallization behavior has not been fully explored and often involves complex multi-stage process. [39][40][41][42][43] The crystallization processes are usually described on the basis of Kissinger method, [44] which provides the average crystallization activation energy, as well as Kissinger-Akahira-Sunose (KAS) method, [45,46] which gives the local activation energy (as a function of the crystallized fraction) as an output. Additionally, the Avrami exponent, [47][48][49] which offers insight into the kinetics and mechanism of crystallization, is also frequently analyzed. [50][51][52][53] However, it is important to note that all of these methods require clearly defined peaks in the differential scanning calorimetry (DSC) curve to produce reliable results, which may be challenging to obtain in the case of complex multi-stage crystallization processes.
The aim of this paper is to produce the Zr 40 Ti 15 Cu 10-Ni 10 Be 25 BMG exhibiting complex multi-stage crystallization process characterized by overlapped peaks and describe its crystallization process through the proposed novel non-isothermal DSC data processing procedure with an unprecedented peak separation method. This procedure includes determining the average and local activation energies by applying Kissinger and KAS methods, as well as subsequent calculations of Avrami exponent in the frames of the modified Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation.

II. MATERIALS AND METHODS
High purity Zr, Ti, Cu, Ni and Be elements were melted in an arc-melting furnace under a protective argon atmosphere to produce the Zr 40 Ti 15 Cu 10 Ni 10 Be 25 alloy. The ingot was subsequently remelted five times to ensure compositional homogenization. The cylindrical rod with a diameter of 3 mm was produced by a suction casting process into a water-cooled copper mold. The fabricated BMG rod was then cut to the shape of discs 1.2 mm thick with the help of a wire electrical discharge machine (WEDM). The low current (about 0.5 A) and water cooling were applied during cutting the rod to reduce the influence of heat on the material structure changes. Finally, the samples were thoroughly cleaned from the residual oxides. The final mass of each prepared disc was close to 60 mg and the mass differences between each sample did not exceed 5 pct. The produced as-cast Zr-based BMG discs were subjected to non-isothermal differential scanning calorimetry (DSC, Netzsch STA 449 F1 Jupiter) measurements in a nitrogen atmosphere using an Al 2 O 3 crucible. The DSC investigations were performed from room temperature (RT) up to 1073 K with a heating rate of 5, 10, 15, 20, 30, and 40 K/min. The temperature and sensitivity calibrations of the DSC device were performed for N 2 atmosphere according to the manufacturer's procedures using high purity reference elements.
The collected non-isothermal DSC data was then subjected to the newly proposed analysis procedure, summarized as a flowchart in Figure 1, to describe the multi-stage crystallization process of the Zr-based BMGs on the example of Zr 40 Ti 15 Cu 10 Ni 10 Be 25 alloy. The recorded DSC curves required the separation of the crystallization stages. To preserve the physical nature of the observed complex crystallization processes and achieve the lowest fitting error, the deconvolution method based on the novel approach with the fitting of Gaussian function and a two-phase exponential decay function convolution was proposed to separate the overlapping peaks. This provides the mathematical description of each peak separately, which is a necessary prerequisite for further use in Kissinger and KAS methods as well as Avrami exponent calculations. Additionally, the improvement in the accuracy of KAS method by the division for the low and high heating rates regions is suggested. More detailed description of the aforementioned methods and procedure stages, together with authors modifications (Figure 1) is provided in Section III of this paper.

A. Non-isothermal DSC Crystallization Investigation
The parts of the recorded DSC curves with different heating rates for the Zr 40 Ti 15 Cu 10 Ni 10 Be 25 alloy in the temperature range 625 to 850 K (close to the crystallization temperatures) are shown in Figure 2. The exemplary glass transition temperature (T g ), observed as an endothermic baseline shift for the curve recorded with a heating rate of 20 K/min, is marked as an insert in Figure 2. It is well visible that the investigated alloy is characterized by three deep exothermic crystallization peaks. Observed crystallization behavior is typical for the amorphous Zr-based bulk metallic glasses. In the literature, the first peak is commonly attributed to nanoscale phase separation by spinodal decomposition with primary crystallization of one phase, whereas the second peak corresponds to secondary crystallization. [54][55][56][57] However, Martin et al. suggest a different approach without prior phase separation. [58] The third peak on the DSC curve ( Figure 2) is associated with the formation of high temperature crystalline phases. [59] The shifting of all peaks towards higher temperatures with the increasing heating rate ( Figure 2) is a result of the required activation energy to overcome the crystallization energy barrier. [44] The first two peaks are significantly overlapped when the heating rate exceeds 10 K/ min. It implies that the second crystallization stage starts before the first one is fully completed. Moreover, the third peak cannot be distinguished for the highest heating rate (40 K/min) as its intensity is extremely low, and it merges with the tailing part of the second crystallization stage.
Overlapping peaks are inherently problematic in analysis. In some research, these peaks are evaluated combined as one, [40,60] but it involves inaccuracies in the obtained results. In another popular approach, the Gaussian [61][62][63][64][65][66][67][68] or Lorentzian [60,69] distribution fitting are used to separate individual peaks and describe each peak shape. However, their application not always accurately describes the peaks shape through cumulative curve [63,64,66,68] or allows to perform the fitting of two subsequent peaks to obtain the cumulative curve, which corresponds to the recorded curve. [65] In extreme cases, it is even not possible to satisfactorily fit those functions to the recorded curves. [60] The main problem is related to the inherent symmetry of Gaussian and Lorentzian distribution, which often does not correspond to the characteristics of the peaks shape. In this paper, a convolution of a Gaussian function and a two-phase exponential decay function into one function was proposed to describe the crystallization processes for each peak observed in Figure 2 for the fabricated Zr-based BMG. The proposed convolution allows the function to describe the asymmetric peaks with the pronounced tailing part which are often obtained in DSC measurements. The proposed function is given as: where f(x) is the two-phase exponential decay function: and h(x) is the Gaussian function part: Fig. 1-Proposed novel non-isothermal DSC data processing procedure; b is the heating rate, T p is the peak temperature, E a is the crystallization activation energy, T a is the temperature of the given conversion fraction, a is the crystallized fraction, E al is the local activation energy, n(a) is the local Avrami exponent.

½3
The final equation after performing the convolution is given as: where limits z 1 and z 2 are given as: and parameters as: y 0 is the y offset, x c is the x center, w is the Gauss width, A 1 is the area 1, t 1 is the decay time 1, A 2 is the area 2, t 2 is the decay time 2.
The above-mentioned approach significantly improves the reflection of the recorded DSC curve shape and reduces the difference between experimental data and the fitted curve. The results obtained from numerical analysis of DSC curves for each heating rate with the corresponding initial DSC curve and the cumulative curve obtained by summation of single peak fits are depicted in Figure 3(a). The strong correlation between the recorded and cumulative curves is well seen. The coefficients of determination (R 2 ) for cumulative and experimental curves are over 0.96 for every heating rate, indicating very good fitting. In comparison, the popular Gauss fitting, despite reporting also very good R 2 for low heating rates (below 20 K/min), reaches only R 2 = 0.88 for 30 K/min curve. Another common fitting with the Lorentz peak function is also characterized by very high R 2 for low heating rates (even 0.99), however, it often misses the true shape of the peaks and distorts the beginning of the first crystallization because of the wide tails. The R 2 of the Lorentz fitting also decreases for the high heating rates reaching 0.92 for 30 K/min. Moreover, as both standard peak fitting functions are symmetric, they have difficulties in fitting the tailing part of the second crystallization, do not reflect the true peak's maximum horizontal coordinate, completely miss the first crystallization for 30 K/min heating rate, and are unable to reasonably fit the 40 K/min. The proposed fitting is fully adaptable to cover all of these conditions with good R 2 , showing superiority to Gaussian and Lorentzian fitting. As a consequence, in this paper, it was possible to determine separately the integrated peak area (Figure 3(b)) for further quantitative crystallization analysis.
On the basis of curves presented in Figures 2 and 3(a), the characteristic temperatures, such as the glass transition temperature (T g ) and the onset crystallization temperature (T x ), along with peak crystallization temperatures (T p ) for each crystallization stage of the produced Zr 40 Ti 15 Cu 10 Ni 10 Be 25 BMG for different heating rates were estimated and summarized in Table I. The width of the SLR was calculated as: On the basis of the obtained results (Table I), it is well visible that each characteristic temperature shifts towards a higher temperature with the increasing heating rate. However, the change in T g is less significant than the change in T x , which results in a widening of the SLR from 73.7 to 99.5 K. The high values of DT x indicate good thermal stability of the investigated alloy.
The obtained values of DT x are in a comparable range as for other Zr-based alloys with similar chemical compositions. [39,42]

B. General Crystallization Activation Energy
On the basis of the determined peak crystallization temperatures and their shifting phenomenon, it is possible to calculate the crystallization activation energy using the Kissinger equation [44,45] : where b is the constant heating rate, T p is the peak temperature of crystallization (K), A is the pre-exponential factor, E a is the activation energy, R is the gas constant.
The specific form of this equation (linear relationship) allows presenting it on a graph, where the left side of the equation relates to the value of the ordinate axis, whereas the inversed crystallization peak temperature is presented on the abscissa. The logarithmic term remaining on the right is a constant value. For each crystallization stage, the T p recorded for a given b was used to calculate the points presented in Figure 4. These points for each crystallization stage were fitted with a    characteristic straight line. The slope of the fitted lines corresponds to the ÀE a /R ratio in Eq. [8].
The fitting parameters along with the calculated E a are listed in Table II. The good linearity of data points is confirmed by coefficients of determination (R 2 ), which exceed a value of 0.995 for each linear fit shown in Figure 4. The first crystallization possesses the lowest E a = 178 kJ/mol and the last one (third crystallization) is characterized by a 2.25 times higher value (402 kJ/mol). The calculated crystallization activation energies for the first two crystallization stages are comparable to those of the other Zr-based alloys with a similar composition reported by Cheng [39] and Yang. [70] On the other hand, it is worth noticing that for the third stage of crystallization the E a for the investigated alloy is about 100 kJ/mol higher than in abovementioned alloys with close compositions. The observed difference in E a for the third crystallization is related to minor changes in Zr-based BMG compositions.
Taking into account the thermal stability of metallic glasses, represented by the amorphous to crystalline transition, the activation energy (first stage of crystallization) for the studied alloy (178 ± 6 kJ/mol) is slightly higher than reported for the (Ti 41 Zr 25 Be 28-Fe 6 ) 93 Cu 7 (155.62 ± 11.58 kJ/mol). [71] However, in vast amount of cases the E a of the studied Zr 40 Ti 15 Cu 10-Ni 10 Be 25 alloy is noticeably lower in comparison with common BMGs, such as Zr 54 Al 10.2 Ni 9.4 Cu 26.4 (340.1 ± 31.7 kJ/mol), [60] Ti 20 Zr 20 Cu 20 Ni 20 Be 20 (286.9 ± 19.6 kJ/mol) [72] or Ti 40 Zr 25 Ni 8 Cu 9 Be 18 (294 kJ/mol). [73] This contrast suggests a lower energy barrier for crystallization in the investigated alloy despite having a wider SLR than the referenced alloys. It is consistent with previous findings that alloys with similar compositions undergo structural decomposition at the nanoscale, which facilitates the crystallization at low temperatures after sufficient time. [54,55] C. Local Crystallization Activation Energy Despite the analytical and comparative convenience, the Kissinger method has inevitable drawbacks. Among others, it reports only one average, temperature independent, activation energy hence missing the complexity of the thermodynamic processes or assumes that maximal reaction rate occurs at the same conversion factor a for all b. This subsequently leads to its low accuracy and oversimplification of the result. [45,74] To account for those problems, the isoconversional Kissinger-Akahira-Sunose (KAS) method, which is based on the reaction conversion in temperature, was proposed.
For proper use of the KAS method, the reaction conversion (crystallized fraction a) as a function of temperature should be determined for each crystallization stage. The previously performed deconvolution of overlapped peaks makes it possible to analyze each peak separately. It should be noted that in the numerical analysis, the integrated peak area between visible deviations from the baseline is proportional to a. Normalized summated fitted peak areas from Figure 3(b) compared to the integrated and normalized overlapped first and second crystallization recorded DSC curves are presented in Figure 5. In this procedure, the normalization was performed to 1, which indicates the end of both processes. The characteristic inflection in Figure 5 divides the curves for first and second crystallization. It is well seen that the curves corresponding to fitted and recorded DSC curves in general overlap with each other for every heating rate. The slight differences (typically less than 0.04 a) in the bottom parts of curves originate from the first crystallization and are caused by the necessary baseline selection during peak integration. The performed fitting inherently have a flat baseline and does not take into account the vertical shift of DSC signal between the T g and onset of the first crystallization. However, the overall good coincidence proves that the proposed deconvolution procedure is valid for the analysis of each peak separately.
Considering the proposed deconvolution method, Figure 6 shows the crystallized fraction as a function of temperature for each crystallization process separately. The typical sigmoidal curves are visible for each crystallization stage, however, the significant asymmetry is observable for the second stage ( Figure 6(b)) due to the different crystallization rates in the initial stages and the end of the crystallization process.
The sigmoidal curves, presented in Figure 6, were used for the KAS method to determine the local crystallization activation energy (E al ) in the process extent. The KAS method is analogous to Kissinger approach, however in this case the peak temperature (T p ) is replaced by the temperature of the given conversion fraction (T a ) for different heating rates [45,46] : On the basis of the Eq. [9] the results from numerical analysis of crystallization processes for every stage of crystallization are shown in Figures 7(a) through (c). All calculations were conducted for a from 0.05 to 0.95 with a 0.05 step. For each a a set of points with different b (b = 5, 10, 15, 20, 30, 40 K/min) was calculated as a base for linear fitting. However, it is clearly seen that the points do not follow a linear dependence, especially at the beginning and the end of the crystallization process. For each case, except for the initial stages of third crystallization, the points turn downward with increasing heating rate, so when they are shifted to a higher temperature. What is interesting, this effect was not pronounced in Kissinger plot. The deviation from linear dependence visible in KAS plots (Figures 7(a) through (c)) is higher for the latest stages of crystallizations (high a) so also for higher temperatures (T). This behavior is consistent with the theory [74] which predicts that the activation energy decreases with increasing temperature. The decrease in E a results in the curving down non-linearity in Kissinger plot as b rises. Similar non-linearity was obtained for other Zr-based BMGs even in Kissinger analysis for a broad b range. [75] As a consequence, to increase the accuracy, the plots in (Figures 7(a) through (c)) were divided into two regions reflecting low (below 15 K/min) and high (above 15 K/min) heating rates. It allows to determine the activation energy dependence in lower and higher temperatures separately. In both ranges, linear fitting was performed for    each a to fulfill the method recommendations. [45] As indicated in Tables III and IV, in most cases, the R 2 coefficients are above 0.99. From the slopes of approximated linear fittings, following Eq. [9], where E al is a part of the slope factor, it is possible to calculate the local crystallization activation energy in a given crystallization conversion.   The coefficients of determination (R 2 ) for each fitting are also presented in the table presented independently for each crystallization step. It is visible that the crystallization activation energy fluctuates with the crystallized fraction, which was not taken into account in Kissinger analysis (Table II). Moreover, there is also a noticeable difference between high and low b ranges. For first and second crystallization, the local activations energies are smaller in high b rage (when shifted to higher temperatures) what is consistent with theoretical predictions [74] and previous studies on different Zr-based alloys. [75] In contrast, the abnormal behavior was observed in the initial part of the third crystallization (Figure 7(f)), which is characterized by significant uncertainties, what can be caused by a complex microstructure of material exposed to previous crystallizations. According to ICTAC Kinetics Committee recommendations by Vyazovkin et al. [45] the E al (a) can be approximated to the average value E a when the difference between the minimum and maximum values of E al is less than 20 to 30 pct of this average. For the first crystallization (Figure 7(d)) the average E al = 185 ± 3 kJ/mol in the low b range is comparable to 178 ± 6 kJ/mol obtained in Kissinger analysis, whereas in the high b range the average E al = 154 ± 5 kJ/mol is noticeably lower. For both ranges, the E al (a) slightly increases with the crystallized fraction and decreases before the end of crystallization (for a ‡ 0.8). This behavior is similar to that reported by Zhuang [75] for the first crystallization in other Zr-based BMG. However, Cheng et al. [39] reported quite different behavior for a similar composition. The difference is connected with the lack of peaks deconvolution and merging the first two neighboring crystallizations.
For the second crystallization (Figure 7(e)), the average E al = 280 ± 18 kJ/mol (calculated according to ICTACS recommendations) in the low b range is comparable to 267 ± 9 kJ/mol obtained from Kissinger analysis and is again significantly lower for the high b range. However, there is a difference in the E al (a) behavior, as for the low b range the E al (a) increases with crystallized volume fraction, whereas for the high b range E al (a) decreases.
In the case of the third crystallization (Figure 7(f))) the average E al = 369 ± 7 kJ/mol in the low b range is rather constant and slightly lower than for Kissinger analysis value (402 ± 16 kJ/mol). In the high b region (15 to 30 K/min), as the crystallized fraction raises, the E al decreases steeply below the values for the low b region, and the value obtained from Kissinger analysis. It should be mentioned that Kissinger analysis reports appropriate average values of activation energies for low hating rates (b £ 15 K/min), however the double linear fitting methodology proposed by authors and shown in Figures 7(a) through (c), revealed that reduction of activation energy for high heating rates (b ‡ 15 K/min) and for higher temperatures is missed in Kissinger analysis.
Generally, it could be stated that for the Zr 40 Ti 15-Cu 10 Ni 10 Be 25 BMG in the low b region, the E al increases moderately or remains fairly constant. The increase can be associated with multiple ongoing processes with different activation energies. In the case of primary crystallization, the change in E al is also connected with an increasing difficulty in the movement of atoms in the remaining amorphous matrix. As nanocrystallization progresses, smaller glassy volumes are available in the structure, which limits the resulting grains sizes. Moreover, the driving force for crystallization decreases as the local diffusion coefficients in the amorphous residual phase decreases. [76] The other factor is the equilibration of the local chemical composition and reduction in concentration gradient. [76][77][78] In combination this results in slowing down crystallization kinetics and an increase in activation energy. The fairly constant value of E al , reported for the latest crystallization stage (Figure 7(f)), is a result of complex crystallization behavior and high temperature reducing activation energy. [74] For the high b range, the crystallization temperatures are shifted to higher values and the temperature reducing effect on activation energy is so significant that E al for each stage of crystallizations decreases with reaction conversion. This could explain the discrepancy between E al (a) behavior in the low and high b range for the second and third crystallization (Figures 7(e) and (f)). The only deviation is visible in the high b range for first crystallization, which takes place in the lowest temperatures. In this case, the activation energy decrease is not so pronounced in comparison with the increasing atoms movement difficulty. Still, there is also a characteristic turning point (at a = 0.7) and a decrease in E al at the end of the process at higher temperature.

D. Local Avrami Exponent
The crystallization kinetics of processes involving crystals nucleation and growth is usually studied within the framework of the JMAK equation. [47][48][49]79,80] The JMAK approach illustrates the crystallized fraction over time and describes the sigmoidal crystallization curves. The model tends to give a fairly good approximation of crystallization kinetics, even if the transformation departs from its initial requirements. [81][82][83] The JMAK model can be expressed by the equation: where n is the Avrami exponent, a t ð Þ is the crystallized fraction vs time, t is the time, s is the incubation time, K is the reaction rate constant.
What is more, the physically meaningful n value describes the kinetics of crystallization and its mechanism [84,85] by the component values [86,87] : where a is the nucleation kinetics index, p is the growth mechanism index, b is the morphology index dimensionality of the growth).
The characteristic values of the indexes in Eq. [11] are listed in Table V. The n values are usually obtained through the linear fitting on a double logarithmic plot. [90] However, this approach provides only a single averaged value, which can be inaccurate considering more complex crystallization processes. Therefore, the concept of local Avrami exponent as a function of crystallized fraction (n(a)) was introduced. [90] Although the original equation was established for crystallization under isothermal conditions, it is possible to extend it to non-isothermal regimes through the Nakamura isokinetic approximation. [91,92] On that basis, Bla´zquez et al. proposed the direct extension of the JMAK equation (Eq. [10]) to non-isothermal conditions [93][94][95] : ½12 Equation 12 allows for the determination of the local Avrami exponent n as a function of a using just one selected heating rate b and raw estimation of E a . Despite Eq. [12] being derived for a constant value of activation energy, it was reported by Bla´zquez et al. [93] that the results can be more accurate while using a varying (local) one.
The second approach for determining the local Avrami exponent n(a) in the extent of crystallization was introduced by Lu et al. [96] by a more straightforward equation: This approach inherently uses the local activation energy concept.
In the presented studies, Eqs. [12] and [13] were used to evaluate the local Avrami exponent on the basis of DSC measurements, calculated crystallized fractions obtained from peak deconvolution (Figure 6), and already determined local activation energies (Tables III  and IV). To ensure the function continuity, the discrete local activation energy values were fitted with the 5th order polynomial function (Figures 7(d) through (f)). Results of n(a) calculations applying both methods for each crystallization stage and heating rate are shown in Figure 8. The values obtained for the low crystallized fractions (a £ 0.1) have been reported to be exceptionally dependent on the DSC baseline selection, [95] as well as numerical instabilities during integration and differentiation. Therefore, this region of uncertainty was marked on the plots with a red color and is excluded from discussion.
It is visible that results obtained from two different methods (described by Eqs. [12] and [13]) are almost perfectly identical, although the Bla´zquez method (Figures 8(a) through (c)) tends to report marginally lower values than the Lu method (Figures 8(d) through (f)). For each crystallization stage, the n(a) behavior is similar, starting from a steep exponential decrease with a nearly linear decline above a % 0.3. Analogous local Avrami exponent evolution characteristic is commonly seen in many kinds of metallic glasses. [62,[95][96][97][98][99] The n(a) dependence in Figure 8 illustrates the variable main crystallization mechanism and slowing down kinetics of crystallization in its every stage. The differences in n(a) values for varied heating rates are significant at the beginning of each crystallization stage (a ‡ 0.1), however they converge in the middle of the crystallization process (a = 0.5) what fulfills the isokinetic (independent of the heating rate) behavior in a proper way. It is worth noticing that in each graph in Figure 8 there is no clear dependence between n(a) and applied heating rate b.
According to Table V and Eq. [11], it should be noted that different crystallization mechanisms (e.g., interface or diffusion-controlled growth) can result in the same Avrami exponent value. [83] Therefore, un unambiguous process description is not possible using only this parameter, especially when it varies and runs through characteristic values and regions responsible for different crystallization mechanisms. [100] In the case of BMG, it is well known that primary crystallization in metallic glasses is mainly driven by volume diffusion [72,78,97,99,[101][102][103] with 3-dimensional growth. [83] Therefore, such a mechanism was assumed in this paper for the first crystallization stage. However, it should be bear in mind that the diffusion-controlled growth may not be the sole instance in the initial phase of the transformation where also the interface-controlled growth and transient stage can occur [76,104] affecting the n(a) analysis. Nevertheless, as mentioned earlier, the initial stages of crystallizations were excluded from analysis anyway. Going to the analysis of first crystallization, the initial values of n(a) for a £ 0.2 are high-in  the range above 8-which indicates the increasing nucleation rate. With the progress of crystallization, n(a) values converge to 2.5 for a % 0.9. Therefore, according to Table V, it is evident that during the substantial majority of the first stage of crystallization, the crystallization mechanism is diffusion-controlled growth with an increasing nucleation rate, which starts to decrease at the end of the crystallization process (a ‡ 0.9). The high share of nucleation is related to the creation of nanocrystalline structure typical for Zr-based alloys with similar compositions. The high value of n(a) at the beginning of the crystallization process shows that the dynamics of the nucleation rate increase is the highest just after crystallization begins. It is consistent with previous findings for similar compositions of BMGs, where the initial spinodal decomposition is followed by rapid onset of nucleation in the amorphous phase. [54,55] The full analysis of second crystallization kinetics is far more difficult as there are no reported data regarding the crystallization mechanism, and it does not result explicitly from the local Avrami characteristics. Moreover, the crystallization process and a number of subsequent exothermic events during the DSC measurement can vary between alloys with different chemical compositions. [59,101] Nevertheless, it can be expected that this crystallization stage is also diffusion controlled. The initial structure (nanocrystals uniformly embedded in amorphous matrix) is complex and different from the resulting one which, in similar alloys, contains different intermetallic phases. [58,105,106] However, interface-controlled growth is characteristic for reactions that do not affect the local atomic composition of the internal structure. [88,89] Though, in the present case, the local atomic arrangement is likely to change during the precipitation of phases with a different composition than the amorphous matrix and nanocrystals, what accounts for diffusion-controlled growth. Moreover, a detailed analysis of n(a) curves in Figures 8(b) and (e) reveals that the local Avrami exponent reaches values below 1 just near the end of the second stage of crystallization which, according to data gathered in Table V, are inaccessible for interface-controlled growth. Consequently, it can be stated that through a major part of the second crystallization (up to a % 0.7) the process is governed by diffusion-controlled growth with an increasing nucleation rate. Initial, remarkably high values of n(a) result from the rapid onset of the crystallization, which can be connected with the preexisting nanocrystalline structure. Further, the crystallization proceeds with the decreasing nucleation rate up to a % 0.85, when it reaches the zero nucleation rate. Then the process is governed only by growth with decreasing dimensionality (Table V). This behavior can be explained by the equilibration of chemical composition, reduction of diffusion coefficients, [76] and impingement, preventing the spherical growth in favor of lower dimensionality growth in non-transformed volumes. For both first and second crystallization, the increasing nucleation rate during most of the process can be connected with a decrease of the amorphous matrix stability with increasing temperature.
In the case of third crystallization, it is not possible to carry out detailed analysis, as it was done for first and second crystallization, due to the creation of unidentified phases before and after. [58,106] However, some conclusions can be drawn from the n(a) behavior presented in Figures 8(c) and (f). Comparing all crystallization stages, it is worth noticing that n(a) exhibits much lower values for third one for a £ 0.7. According to Table V, it indicates a more growth-oriented crystallization process. The majority of the transformation in Figures 8(c) and (f) is represented by the n value in the range 1.5 to 3. These values point at zero nucleation rate and decreasing dimensionality of the growth in the case of interface-controlled growth, or initially increasing and then decreasing to zero nucleation rate in the case of diffusion-controlled growth.

IV. CONCLUSIONS
The crystallization behavior in the studied Zr 40 Ti 15 Cu 10 Ni 10 Be 25 BMG exhibits three exothermic events observed in non-isothermal DSC measurements. The first two peaks significantly overlap (for DSC curves recorded with high heating rates), what shows the complexity of the crystallization process in the investigated alloy and requires a sophisticated novel approach to analysis. Evaluation of overlapped peak is inherently troublesome, however the proposed non-isothermal DSC curve deconvolution based on a Gaussian function and a two-phase exponential decay function allows for their accurate independent analysis. The suggested approach provides crystallized fractions which are in good agreement with the ones obtained from the raw curves and can be successfully used for other BMGs characterized by complex crystallization behavior.
The crystallization activation energies for first, second, and third crystallization obtained by commonly used Kissinger method are 178 ± 6, 267 ± 9, and 402 ± 16 kJ/mol, respectively. The first value, corresponding to primary crystallization, observed in the investigated alloy, is smaller than that for other Zr-based BMG with narrower SLR. Another approach utilizing more advanced Kissinger-Akahira-Sunose isoconversional method, with a proposed in this paper unprecedented division for the low and high heating rates region, showed that activation energy dependence is more complex than a single value obtained from the Kissinger method, and can be expressed as a function of crystallized fraction. Kissinger analysis reports average values of activation energies which are in good agreement with KAS method only for low heating rates. However, their reduction for high heating rates and higher temperatures is missed.
The calculated local Avrami exponent is the highest at the early stage of each crystallization, and then decreases, indicating slowing down crystallization kinetics. Based on proposed analysis utilizing the Avrami exponent components, it is seen that the first two stages of crystallization are diffusion-controlled with mainly increasing nucleation rate (Avrami exponent higher than 2.5). However, at the end of the second crystallization stage, the nucleation rate reaches zero what is followed by a decrease in the dimensionality of the growth. The third crystallization stage described by lower values of Avrami exponent is more growth-dominated.