Modeling the Thermal Stability of the $α/ω$ Microstructure in Shocked Zr: Coupling between defect state and phase transformation
Stephen R. Niezgoda Thaddeus Song En Low






l1]T.S.E. Low
l1,l2]S.R. Niezgodacor1 
[cor1]Corresponding Author
niezgoda.6@osu.edu


[l1]Department of Materials Science and Engineering, The Ohio State University, Columbus, OH 43210, USA
[l2]Department of Mechanical and Aerospace Engineering, The Ohio State University, Columbus, OH 43210, USA


Under high pressure, Zr undergoes a transformation from its ambient equilibrium hexagonal close packed α phase to a simple hexagonal ω phase. Subsequent unloading to ambient conditions does not see a full reversal to the α phase, but rather a retainment of significant ω. Previously, the thermal stability of the ω phase was investigated via in-situ synchrotron X-ray diffraction analysis of the isothermal annealing of Zr samples shocked to 8 and 10.5 GPa at temperatures 443, 463, 483, and 503 K. The phase volume fractions were tracked quantitatively and the dislocation densities were tracked semi-quantitatively. Trends included a rapid initial (transient) transformation rate from ω→α followed by a plateau to a new metastable state with lesser retained ω (asymptotic). A significant reduction in dislocation densities in the ω phase was observed prior to initiation of an earnest reverse transformation, leading to the hypothesis that the ω→α transformation from is being hindered by defects in the ω phase. As a continuation of this work, we present a temperature dependent model that couples the removal of dislocations in the ω phase and the reverse transformation via a barrier energy that is associated with the free energy of remaining dislocations. The reduction of dislocations in the ω phase occurs as a sum of glide and climb controlled processes, both of which dictate the transient and asymptotic behavior of the annealing process respectively. 



Zirconium
Omega Phase
High Pressure
Phase Transformation
Model
Annealing







§ INTRODUCTION

At ambient conditions, the group IV transition metals (Zr, Ti, Hf) are thermodynamically stable as a hexagonal close packed (HCP) crystal (P6_3/mmc, c/a=1.593), denoted as the α phase. At high pressure, the α phase in these metals undergoes a transformation to the ω phase, which has a simple hexagonal crystal structure (P6/mmm, c/a=0.623). After shock loading <cit.> or severe plastic deformation by high-pressure torsion (HPT) <cit.>, the ω phase displays a strong hysteresis; with a significant fraction remaining after pressure is released. In the case of ω material yielded from shock impact experiments, samples with retained ω fractions as high as  80% of ω phase have been observed <cit.>. The resulting two-phase microstructure is (meta)stable for years at standard temperature and pressure. 

Over the past 60 years, there have been numerous studies on the α→ω forward transformation, with the majority focusing on the determination of the equilibrium temperature/pressure phase diagram, crystallography, mechanical properties, and electronic structure of the ω phase <cit.>. The significant disagreement amongst these studies attest to the influence of the transformation hysteresis on the determination of basic physical properties of group IV transition metals. As an example, the equilibrium transformation pressure for Zr has been reported to range from 2–7  GPa for the static pressure loading case <cit.>. However, even though the metastability of the ω phase is well documented, there is only a limited body of work on the stability of the α/ω dual phase microstructure. 

The α→ω phase transformation in Zr has been shown to be achievable by several techniques, namely 1) hydrostatic pressure, 2) shock loading, and 3) high pressure torsion methods. It is not clear whether there is a single mechanism for the α→ω transformation or if multiple transformation pathways exist.   Depending on the experimental conditions and loading mechanisms, multiple (α/ω) orientation relationships have been reported. The amount of retained ω, transformation pressures, and the apparent degree of ω stability vary significantly between experiments <cit.>. Hence there is a need for clear need for systematic studies of both the forward transformation and of the reverse ω→α reverse transformation. More specifically, models that link the evolution of the microstructure to thermodynamic state variables, coupled with in-situ experimental observation, are required to characterize the kinetics of the transformation under various conditions in order to provide evidence for or against proposed transformation or reverse-transformation mechanisms <cit.>. 

In this paper, emphasis is placed upon the thermal stability of the ω phase generated from shock experiments. The ω phase produced from shock impact is generally characterized by the  (0 0 0 1)_α∥(1 0 1 1)_ω and [1 0 1 0]_α∥[1 1 2 3]_ω orientation relationship identified by Song & Gray <cit.>.  Brown et al. first investigated the thermal stability of the ω phase by heating shocked samples under a constant temperature ramp from room temperature to 620 K at a rate of 0.05 K/s <cit.>. A subsequent study by Low et al. involved the isothermal annealing of shocked Zr samples over a similar temperature range <cit.>. Both experiments showed a significant decrease in dislocation content prior to the onset of the reverse transformation, leading to postulation that the metastability of the ω phase was highly correlated with the defect state of the microstructure. Furthermore, it was speculated that the dislocation population was somehow arresting the reverse transformation and that given sufficient thermal energy, the defects would gain mobility and annihilate through recovery mechanisms; allowing the reverse transformation to proceed. Here, we present a model based on this hypothesis and demonstrate that model predictions are consistent with the experimentally observed transformation kinetics. Section <ref> reviews and summarizes the previously collected experimental data which informs the current work as well as a brief discussion of  previous modeling efforts. This is followed by section <ref> which provides detail into the model, after which section <ref> provides the results from the model. Further discussion on the results and their implications are provided in section <ref>.




§ EXPERIMENTAL METHODS AND PRIOR STUDIES
 




 §.§ Production of Metastable ω Through Shock Impact

The data used to develop and validate this model came from annealing studies of shocked Zr performed by Low et al. <cit.> and Brown et al. <cit.>.  The shocked material for those studies came from gas gun experiments performed by Ceretta and collaborators <cit.>. Prior to shock loading, all Zr samples were prepared from a high-purity crystal bar Zr (<100 ppm impurities) which was upset forged, clock rolled, and annealed, resulting in a plate with homogenous and fully recrystallized microstructure, with average grain size of 15-20 μ m <cit.>. The plate exhibited a strong basal texture (>8 times uniform random distribution) nearly aligned with the normal or through-thickness direction (TT) direction of the plate, with prism planes uniformly distributed about the in-plane directions of the plate (transversely isotropic). Shock loading by dynamic impact was performed utilizing a gas-driven 80-mm single stage launcher <cit.>. The impact direction was along the TT direction or parallel to the c-axis of the majority of grains in the sample. The targets were impacted by 2.5 mm thick Zr flyer plates accelerated to velocities of 640 or 835 m/s, resulting in peak compressive stresses of 8 or 10.5 GPa respectively on the Zr samples <cit.>. The total duration of the shock process (plasticity and transformation included) was approximately 0.5 μs <cit.>. The samples were soft-recovered and were further sectioned for X-ray and microscopy. The real-time process of the shock deformation was tracked via velocity interferometer system for any reflector (VISAR) measured wave profiles. The wave profiles and their interpretations are given in figure <ref>. The VISAR traces indicate that samples shocked to 10.5 GPa complete the α→ω transformation over a significantly shorter time duration than samples shocked to 8 GPa. It was also speculated that the ω phase in samples shocked to 10.5 GPa experienced additional plastic deformation, in addition to the higher pressure, resulting in a difference in defect population compared to the 8 GPa shocked samples.  The resulting volume fractions of retained ω phase are ∼ 60% and ∼ 80% for samples shocked to 8 and 10.5 GPa respectively. Complete details of the shock experiments can be found in <cit.>. 





 §.§ In-situ Annealing Experiments

The in-situ annealing experiments of the α/ω samples consist of work by Brown et al. <cit.> and Low et al. <cit.> and were performed via x-ray diffraction (XRD) characterization methods at the 1ID-C beam line at the Advanced Photon Source (APS), Argonne National Laboratory <cit.>. The experimental setup was typical for powder diffraction on a polycrystalline sample. A high energy monochromated incident beam (E = 86 keV) was used to illuminate the samples and the resulting diffraction images were collected on a a two-dimensional detector placed at an appropriate distance to capture at least 5 diffraction rings from each phase<cit.>. The individual (α/ω) phase fractions were tracked quantitatively and the dislocation densities were tracked semi-quantitatively. 

The experimental work by Brown et al. involved the heating of Zr samples shocked to 8 and 10.5 GPa from 300 to 620 K at a constant temperature ramp rate of 0.05 K/s. Significant ω→α transformation was observed within the temperature range of 475 < θ < 550 K  <cit.>.  Low et al. further investigated the transformation under isothermal conditions at temperatures of 443, 463, 483 and 503 K in samples similarly shocked to 8 and 10.5 GPa <cit.> and observed a rapid initial transformation rate. At all temperatures there was a continuous deceleration of the transformation, with the transformation rate approaching zero with significant ω phase remaining. This new metastable state appeared to be asymptotically stable. The higher the temperature applied, the further the extent of the reverse transformation (see figures <ref> and <ref> in section <ref>). Both phases contain high dislocation densities (≈10^15 m^-2) prior to annealing and observe a reduction in dislocation densities during annealing. Furthermore, in the initial stage of annealing, a significant reduction of the ω phase dislocations was observed prior to onset of the reverse transformation, prompting the authors to postulate that the ω→α transformation was arrested by dislocations in the ω phase. The experiments also noted that at temperatures 483 and 503 K, the 10.5 GPa samples observed a more substantial reverse transformation such the α volume fractions surpassed those of the 8 GPa samples, despite the 10.5 GPa samples starting off with a lower α volume fraction. This cross-over trend was also observed when applying different temperature ramp rates from room temperature to 503 K on samples shocked to similar conditions. EBSD characterization of ex-situ annealed samples indicated that the reverse transformation proceeded primarily by growth of existing α laths at 443 K and below and nucleation and growth of new α laths at higher temperatures (463–503 K). 

The annealing experiments raised an important question. The ω phase is very far from equilibrium at atmospheric pressure. The occurrence of a rapid transformation at low homologous temperatures (T_m ≈ 2128 K for Zr) indicated the ω was only weakly kinetically frustrated and that application of relatively small thermal driving forces was sufficient to overcome the transformation barrier. Yet, the asymptotic behavior, in which the transformation plateaus to another metastable equilibrium with retained ω, suggests a more stable arresting mechanism. The question posed by the previous experiments can be phrased briefly as “What is the source for this type of hysteresis and how does it simultaneously explain both the transient and asymptotic behavior of the transformation?” 



 §.§ Prior Modeling and Simulation


Nisoli et al. proposed a model which predicted the asymptotic stability of the ω phase under isothermal annealing, however this model did not address the initial rapid transformation rate during the initial stages of heating <cit.>. Nisoli et al. evolved the dislocation densities based on the Kocks and Mecking formulation <cit.> that relates the dislocation evolution to shear evolution via a phenomenological model:


    dρd γ = c_1√(ρ) - c_2ρ.
 


The source of shear was assumed to be a combination of  thermally activated plasticity and phase transformation (transformation shear). For the plastic shear contribution, Nisoli et al. takes the Orowan equation <cit.> γ̇_p = ρ_m b v, implying that plastic shear occurs by the movement of mobile dislocation populations present in the ω and that in turn contributes to dislocation elimination via the Kocks-Mecking relation listed above. Coupling between the dislocation evolution and the phase transformation was done in a similar manner to the modeling of martensitic transformations in TRIP-steels by Stringfellow et al. <cit.>, in that the phase volume fraction v_ω is related to the shear γ via: 


    v̇_ω∝ - v_ωγ̇.


The reverse transformation rate and the dislocation evolution were described to be mutually coupled, in that occurrence of both processes lead to a net shear which in turn drives further dislocation evolution and phase reversion. One can think of this as a feedback loop in which the input and output continuously drive another. In the context of an asymptotic framework, this description is adequate in that the initial work required to kick-start the microstructural evolution has been done and does not require specification. While the model adequately captures the long-time behavior of the α/ω microstructure, it does not address the initial transient phase of the transformation.










In a different study, Zong et al. performed a molecular dynamics simulation of the reverse transformation. The kinetics of the reverse transformation were found to follow a modified Kohlrausch–Williams–Watts equation, which is  sometimes used to describe “glassy” or frustrated kinetics <cit.>. Zong was also able to extract  a relaxation time constant which could be related to the thermal activation energy <cit.>. In said work, the authors noted that (a) the effective activation energies were higher in samples shocked to higher peak shock pressures (8 vs. 10.5 GPa); (b) the transformation can be modeled as a thermally activated process resulting from heterogenous nucleation from defects (ie. α/ω interface, dislocations, grain boundaries, etc.). The molecular dynamics study of Zong suggested that the reverse transformation via interfacial growth of existing martensitic α laths was less likely to occur than the defect mediated heterogenous nucleation process at all temperatures, which was inconsistent with the ex-situ annealing experiments from Low et al. <cit.>.  This apparent disagreement between molecular dynamics predictions and experimental observation indicates that a more detailed description of the defect states and their role in the reverse transformation is warranted. 



§ MODEL FOR THE EVOLUTION OF THE TWO-PHASE Α/Ω MICROSTRUCTURE
 
In this work we present a temperature dependent model that couples the static recovery of dislocations with that of the ω→α transformation. Instead of coupling both processes via the evolution of shear, the dislocation density reduction is driven by the system's tendency to reduce the overall free energy. This reduction is subsequently coupled with the reverse transformation using an expression that describes the resistance to the reverse transformation in terms of the remaining dislocation content. 



 §.§ Dislocation Density Evolution Model

At the beginning of the isothermal heating experiments, the post-shock defect state of the ω phase contains a high dislocation density arranged in a complex tangled dislocation network <cit.>. Under an annealing process, the removal of dislocation occurs by recombination or elimination reactions, which ultimately leads to the growth and maturation of these dislocation networks <cit.>. The proposed model for the evolution of the dislocation network is based on the framework of Nes <cit.> for the modeling of static recovery in iron, aluminum, and AlMg alloys. Nes considers the reduction of dislocation density through the growth and refinement of dislocation networks; which occurs through the processes of thermally activated glide, thermally activated cross slip, climb and solute drag as potential rate limiting mechanisms.

Following Nes, we assume these reactions have different rates depending on the underlying mechanism of defect annihilation. Given the high purity of the Zr samples we can neglect solute drag. We have also decided not to consider thermally activated cross slip, as we do not have a sufficient understanding of the likely deformation modes in ω and the relative mobilities of different dislocation populations especially in the as-shocked condition. Rather than introducing ad-hoc fitting parameters to account for the relative ease of the different mechanisms, we have decided to simplify the framework and only consider two rate limiting mechanisms:

	
  * glide controlled: based on glide mobility of heavily jogged screw dislocations.
	
  * climb controlled: based on the climb of edge dislocations.


Our supposition is that during the early stage of the annealing process, sufficient energy is provided to activate network growth that is glide controlled. Upon exhaustion of defect populations that eliminate by means of glide controlled processes, the network growth then proceeds by removal of dislocations by a climb controlled process. Note that the rate limiting step for thermally activated glide is the motion of jogs in screw dislocation segments. Therefore, both glide and climb controlled reactions are fundamentally rate-limited by the climb of edge segments. 

In network growth, the average distance between dislocations, r, increases at a rate of 


    drdt = k_1 v


Here, k_1 is a proportionality constant that, ideally, would be characterized experimentally. v is the average velocity of dislocations
which depends on a driving force associated with a decrease in free energy by network growth: 


    F = k_2
    μ b^2r
 

where k_2 is a constant of order unity, μ is the bulk modulus, and b is the burgers vector magnitude. 
 


  §.§.§ Glide Controlled Elimination




To derive an expression for the glide speed v, we consider that in a high purity metal, the rate controlling mechanism is expected to be the climb of jogs with edge character <cit.>. Consider an initially straight screw dislocation that contains jogs of alternating signs separated at a length of l_j along the dislocation, as shown in figure <ref>. Under the driving force F a local vacancy equilibrium will be established at each jog 


    c=c_0 exp±(Fbl_j/rk_Bθ)


where k_B is the Boltzmann constant, θ is the temperature. The ± indicates that jogs of opposite sign will respectively either emit or absorb vacancies during climb. The consequent presence of concentration gradients means that the rate of climb will be controlled by vacancy diffusion between jogs of opposite signs <cit.>. In the case of bulk self diffusion at lower stresses, one can expect the dislocation mobility to be controlled by the slower moving vacancy absorbing jogs. The vacancy producing jogs also move while leaving vacancies behind, but would exert a weaker dragging force on the dislocation <cit.>. Considering this, the glide speed can be expressed as: 


    v = 4π k_3D_sbsinh (Fbl_jk_Bθ ) ≃4π k_3D_sbexp (Fbl_jk_Bθ )


where D_s is the bulk self diffusivity and k_3 is a constant that needs to be determined experimentally. The sinh function has been replaced with the exp function on the assumption that Fbl_j/k_Bθ > 1. Combining the glide speed expression with equation (<ref>), the network growth rate that is climb controlled is given as: 


    ( drdt)_glide  = 4 π k_1 k_3 D_sbexp( Fb l_j0k_Bθ)


l_j0 is the initial average jog spacing of the heavily jogged screw dislocation. For simplicity we have chosen to keep the average jog spacing constant for this model. The physical meaning of this assumption is that the lateral-jog-glide-mobility is low and the frequency of annihilation of opposite signed jogs is also low. This is a reasonable assumption in metals with low stacking faults or in solid solution strengthened alloys where glide is impeded by dislocation-solute interactions. Neither of these cases directly apply to  ω-Zr as the SFE of ω≈ 479 mJ/m^2 based on atomistic modeling <cit.>. While this may be a less than satisfactory assumption, there is not a physically motivated way to estimate the lateral-jog-glide-mobility. Rather than introducing an additional fitting parameter in the model to control this behavior, we felt that the assumption of a constant jog spacing was a better choice based on the limitations on the current state of understanding of ω phase defects. The jog spacing is related to the average dislocation separation by: 


    l_j0 = κ r_0
 

where κ is a proportionality constant, and r_0 is the initial average separation between dislocations. κ can be interpreted as a proxy of the extent of deformation experienced during the shock – a lower value of κ indicates a higher extent of prior deformation. A higher proportion of mobility limiting jogs indicate relatively more interactions during plasticity of the shock deformation. We would expect κ to be lower in shocked samples than in quasi-static deformed samples. 



  §.§.§ Climb Controlled Elimination

In the case of the annihilation of two edge dislocations, the movements of the dislocations prior to annihilation can be conceptually decomposed into that of a pure glide of length l_g followed by a pure climb of length l_c. By treating the overall glide/climb movement as a summation of elementary activations (each of climb step b), the activation length can be written as l_a = ξ b, where ξ = 1+l_g/l_c. Again assuming a vacancy diffusion controlled climb for bulk diffusion, the average dislocation speed for climb control can be expressed as <cit.>: 


    v = 2 πξ D_s b ln ( 2rb) [ exp ( Fb^2k_B θ ) - 1 ]


which upon combination with equation <ref> results in the following network growth rate:


    (drdt )_climb   = 2 π k_1 ξ D_s b ln ( 2rb) [ exp ( Fb^2k_B θ ) - 1 ]
       = 2 πk̃ D_s b ln ( 2rb) [ exp ( Fb^2k_B θ ) - 1 ]
 

It is worth mentioning that the determination reasonable values of ξ is non-trivial. Again, to reduce the number of fitting parameters in the absence of experimental values ξ and k_1 have been combined into k̃. 



  §.§.§ Dislocation Evolution

The total growth rate of the dislocation networks is then just a sum of the aforementioned contributions: 


    drdt =  (drdt )_glide +  (drdt )_climb


Approximating the dislocation density as ρ≈ (1/r^2), we are able to  evolve the dislocation density with respect to time via: 


    dρdt = -2ρ^3/2drdt
 

from which the dislocation density can be used to calculate the microstrains via:


    ε = √(ρρ_0)ε_0


where ρ_0 and ε_0 are the initial dislocation density and experimentally measured microstrain respectively. 




 §.§ Phase Transformation Model

The rate at which the ω phase volume fraction f_ω varies with time is given as a double exponential rate equation:  


    ḟ_ω = ḟ_0 f_ωexp (  α_1Δ GR θ )exp ( β_1 B(ρ)k_B θ )


where ḟ_0 is a reference transformation rate and R is the gas constant. Conceptually, one can think of the first exponential term being that of the driving force for the ω→α transformation due to the Gibbs free energy differential (Δ G)  at the current temperature and pressure. This driving term is assumed to be independent of the current dislocation state of the ω phase. α_1 is a unitless parameter which increases non-linearly with temperature and can be interpreted as an indicator for the specification of the ω→α transformation mechanism (ie. growth of existing α laths vs nucleation and growth of new α plates). 
Δ G is evaluated using an equation of state (EOS) for Zr detailed by Greeff <cit.>. The EOS for Zr provides the Helmholtz free energy, H^eos, via: 

    H^eos(V,θ) = ϕ_0(V) + H^ion(V,θ) + H^el(V,θ)
 
where ϕ_0 is the static lattice energy, H^ion is the ion motion free energy and H^el is the electronic excitation free energy. According to Greeff, ϕ_0 dominates the pressure or volumetric contribution to the Helmholtz free energy expression, while H^ion has the strongest temperature dependence. For the analytical expressions and calibrated parameters for these individual contributions, and the reader is referred to the work of Greeff <cit.>. Consequently, the Gibbs free energy for each phase can be calculated as a function of specific-volume and temperature via:

    G = H^eos + Pv

where P is the pressure and v is the specific volume of each corresponding phase. 

The second exponential term in Eq. <ref> describes the barrier energy, due to the defect state of the microstructure, that resists the ω→α transformation under heating. The barrier energy is related to the energy required to mobilize the arresting dislocations, subsequently allowing elimination reactions to occur. Not all dislocations will have the same mobility or be equivalent in their propensity to arrest the transformation. Here we make the assumption that effectiveness of a dislocation at resisting the transformation is directly related to the energy requirement for mobilizing that defect. Experimentally, Low et al. observed a large decrease in the ω phase dislocation density before any significant transformation occurred. Here, we postulate that the highly mobile defects exert little to no arresting influence and due to their high mobility are removed early during the recovery process. By this logic, the barrier energy or the resistance to the ω→α transformation should increase as the dislocation density is reduced as the remaining dislocation population are less mobile and more effective at retarding the transformation. The barrier energy is evaluated as follows:


	
    B(ρ)    = E_climb + E_glideρ
    
    	E_climb   = Fb^2
    
    	E_glide   = Fbl_j

	
 

The expressions for the individual energy contributions are from the exponential arguments of the individual dislocation network growth rate equations (eq. (<ref>) and (<ref>)). Note that the energy terms in equation <ref> are not constant, as the driving force F varies inversely proportional to the average dislocation separation. Taking into account that the barrier energy B represents an averaged activation energy for mobilization of individual dislocations, the β_1 parameter can be thought to be a scaling of this energy to represent a reference defect population providing this resistive energy, and is on the order of 10^-1ρ_ref≈ 10^14.



 §.§ Model Parameters
 
Table <ref> is a summary of the parameters used for the dislocation density evolution model. The parameters listed here are calibrated to the fitting of the annealing of 8 GPa samples.




The value of the exponential pre-factors (k_1, k_3, k̃) from equations (<ref>)and (<ref>) are not readily obtainable from experimentation. Even though they have a clear physical significance there is no way to directly estimate the parameters or constrain them with physically motivated bounds. Since these terms affect the outcome of the calculation of the network growth rate in a similar manner to the self diffusivity D_s by acting as an exponential pre-factor, we can consider the product of the these parameters in Eq. <ref> to be an effective diffusivity and verify that it is consistent with the limited experimental observations. Zr is known to have vanishingly low self diffusivities at low homologous temperatures (≈10^-23 at 800 K) <cit.>. However, it has been reported that precipitation of the ω phase has been correlated with an anamolous increase in diffusivity <cit.>. For this work we have chosen, somewhat arbitrarily, for the diffusivity at 443 K to be ≈10^-22.  By keeping k_1, k_3, k̃ as temperature independent parameters, we can extract the trend in diffusivities with increasing temperature.  

In previous work, Brown et al. assumed the dislocation density as being proportional to the microstrains reported from the XRD analysis <cit.>. For the α phase, this can be made quantitative through the application of line profile analysis. Line profile analysis requires detailed knowledge of the available slip systems and likely dislocations present <cit.>; data which is not currently available for ω. For the ω phase, a problem arises in that we have quantitative values of the microstrain but no direct measurement of the dislocation density within the ω phase through conventional line profile analysis <cit.>. This means that the initial ω phase dislocation density needs to be estimated. We chose the 443 K sample shocked to 8 GPa as a baseline reference sample. Based on a measured α phase dislocation density of  ρ^α_0 ≈ 9× 10^14 m^-2 <cit.>, we chose to set the initial ω dislocation density to be ≈ 10^15 m^-2 for this sample. By taking the initial values of the microstrain and dislocation density of the 443 K sample shocked to 8 GPa as reference microstrain ε_ref and reference dislocation density ρ_ref, respectively, the initial dislocation densities in the other samples can be calculated as:  

    ρ_0 =  (ε_0ε_ref )^2 ρ_ref




In table <ref>, a list of parameters related to the transformation model is listed. 





§ RESULTS
 


 §.§ Calibration to Isothermal Annealing Experiments
 
Figure <ref> illustrates model fits to the isothermal annealing of the 8 GPa sample. Using the parameters provided in section <ref> as a baseline, the remaining parameters that remain unknown are the self diffusivities of the ω phase that influences the dislocation removal rate and the α_1 parameter of the transformation rate equation (equation <ref>). Due to any lack of experimental reports for these parameters, the temperature dependent trends of these two parameters are treated as model outcomes and potential validation points should the diffusivities be measured or calculated in the future. The extracted values for these two parameters are provided in table <ref>. For the microstrains illustrated in figure <ref>, the fits for the higher temperatures showed good correspondence with the samples annealed at higher temperatures (463-503 K) but showed a mild deviation in the curvature at 443K. The fits for the ω were also satisfactory and able to capture the initial rapid transformation and asymptotic behavior. 








The results of exercising the model on the 10.5 GPa shock samples is shown in figure <ref>. The parameters calibrated to the 8 GPa samples were retained, with only the initial condition (e.g. volume fractions and ω dislocation density) varied. As can be seen from the figure, the predictions of the ω phase fractions are reasonable. However, the model over-predicts the reduction of microstrain for all temperatures except 443K. Low et al. observed that the evolution of microstrain (and by extension dislocation density) of the 10.5 GPa samples was fundamentally different than the 8 GPa samples, leading to a speculation that the defect populations of the two samples were qualitatively different likely due to more extensive ω phase plasticity during the shock event <cit.>. The model predictions are consistent with this hypothesis and suggest that within the 10.5 GPa samples there is a population of defects that are not readily annealed out at these low homologous temperatures but is also not actively resisting the ω→α transformation. This point is discussed further in section <ref>.








 §.§ Effect of Jog Spacing


The poor predictions for the ω volume fraction for the 483K 10.5 GPa sample (see figure <ref>) can be sensibly and simplistically remedied by a small change in the value of κ from 0.023 to 0.0215 as shown in figure <ref>. The variation of κ has a significant effect on the initial transient behavior of the ω volume fraction. As κ decreases, the initial transformation rate increases. From a physical understanding, the adjustment of the κ parameter is equivalent to varying the jog spacing on screw dislocations, which concerns the dislocation population that eliminates by glide controlled mechanisms. In a more deformed sample, one would expect a higher frequency of jogs along its dislocations. Again, this is consistent with the notion that the 10.5 GPa samples have undergone more severe plastic than the 8 GPa samples. As the degree of deformation is increased, the heterogeneity is also expected to increase. The in-situ annealing samples were small 3mm disks cut from the large impact specimen <cit.>. It is not surprising to find spatial variation in defect state in the shocked samples and that the 8 GPa samples are more homogeneous than the 10.5 GPa samples. From this analysis it would seem that the jog spacing, κ, could serve as a good indicator of the plastic history of the samples.






 §.§ Validation: Linear Temperature Ramp

As a validation of the framework, the model was exercised on the constant temperature ramp data extracted from Brown et al. <cit.>. As described above all model parameters except the diffusivities and α_1 were considered to be temperature independent. In order to match/simulate the heating from room temperature to 620 K, the values for the diffusivities and α_1 had to be extrapolated beyond the range given in Table <ref> to higher and lower temperatures. The self diffusivities were fit to a two term exponential equation. For the α_1 parameter, a piecewise linear interpolation has been applied for the temperature range of 443-503 K. For temperatures in the range (503 < θ≤ 620 K) , this fit was continued by linear extrapolation to the point (θ, α_1) = (620, 24). For temperatures less than 443 K, a linear extrapolation to the origin was chosen. Remember that the α_1 parameter acts as a modifier to the bulk Gibbs free energy based to account for potentially different transformation mechanisms over different temperature ranges. 




The model predictions vs. the Brown et al. data is displayed in figure <ref>. The figure visually confirms an excellent consistency between the model and the experimental values. With the exception of κ, all other model parameters remained fixed at the values calibrated for the 8GPa isothermal annealing case. κ=0.217 and κ=0.210 was used  for the 8 and 10.5 GPa constant heating ramp samples respectively. The changes represent a difference of approximately 6% to our the isothermal 8 GPa values. Reiterating that the κ value represents the jog spacing on screw dislocations, a 6% difference in value is quite reasonable due to the annealing samples being obtained from different regions of the larger shock loaded Zr plate (Also the annealing samples were removed at different times, more than one year separated the sample prep for the Brown at al. and Low et al. studies). Differences of a few percent in the initial volume fraction of ω phase also points to the samples from both experiments seeing slightly different deformation conditions as well. 








§ DISCUSSION
 


 §.§ Glide vs. Climb Controlled Dislocation Removal
 
One of the key features of the current model is the concept of elimination of dislocation populations by two modes, more specifically the glide and climb controlled removal processes. Nisoli et al. had previously modeled the reverse transformation during annealing in the limit of long times <cit.> and an overview of that work was presented in section <ref>. Nisoli et al. starts with a Kocks Mecking formulation in which the dislocation density is evolved in relation to shear, and considers two sources for shear: i) thermally activated dislocation motion and ii) transformation shear. Nisoli further makes the assumption, based on experimentally observed transformation rates, that the transformation rate decreases to zero as time goes to infinity and all the shear can be assumed to come from thermally activation of dislocations at asymptotic times. By considering a Taylor expansion of the Kock Mecking equation around the point ρ(t)=ρ(∞),  Nisoli et al. eventually derives an expression for the dislocation evolution at long times, as shown: 

    ρ = ρ_0(1+t/τ)^1/ψ

where is a τ is a critical timescale which scales with the temperature dependent mobility of the dislocation population. This equation provides a very crucial and insightful link, as it is of the same form that Nes derives for climb controlled dislocation elimination in annealing processes <cit.>. Nes also notes that the climb controlled elimination processes dominates when glide contribution to the dislocation removal is exhausted, therefore providing a key to the transient vs asymptiotic behavior of the dislocation removal trends. 

A more detailed picture of the dislocation reduction process is shown in figure <ref>, which illustrates the microstrain evolution in the ω phase if only one of the two dislocation reduction processes (glide and climb) were included in the model. It can be observed that at lower temperatures (443 K), the glide process dominates strongly in the initial stage but eventually tapers off, in turn letting the climb process dictate dislocation reduction trends at long times. At higher temperatures (503 K), the contribution of climb appears to be more significant and consistently dominates the dislocation removal trends, whereas the red line depicting glide controlled populations sees a rapid exhaustion (figure <ref>). Naturally, the rate of the process rate limited by climb follows an Arrhenius trend with respect to temperature and dominates more at higher temperatures. At sufficient temperatures, if the process is completely climb controlled, standard Avrami kinetics for the transformation should be recovered. Nisoli was able to demonstrate that this indeed occurs for dislocation evolution processes with forms given by equation <ref>. Our proposed model and the model of Nisoli et al. both give the same predictions for long times and high temperatures. The key difference is that our model ascribes the behavior to a specific dislocation mechanisms, climb controlled recovery, as opposed to generic thermally activated dislocation motion. 

Since the dislocation reduction trends are coupled to the reverse transformation, the transient behavior of the reverse transformation is associated to the rapid exhaustion of dislocations that eliminate by glide controlled processes. This effect is more obvious at lower temperatures, when the climb controlled processes are more passive and insufficient to drive a rapid elimination of dislocations in the initial stages of annealing. This is also shown in the work of Nisoli et al., in which fits to the ω decay at low temperatures were unable to capture the initial rapid descent behavior.

It is also important to consider the range of applicability of the current model given our choice to only consider glide and climb controlled recovery mechanisms. Of particular interest is determining if the current framework can be used to assess the ω→α reversion that occurs during the immediate unloading following the shock event. Both dislocation removal mechanisms in the model are essentially controlled by climb of edge segments. Here the assumption is that the time-scale on which individual climb events occurs is short compared to the annealing time and we can consider the rate of occurrence in a statistical sense relative to the total dislocation density. In contrast the shock process provides high strain rates over a period on the order of a  μ s. It is unlikely that dislocation removal via climb leads to the reversion during the shock loading due to the longer time scale under which climb occurs. The current framework works under a longer time length scale following the shock, in which the pinning effect of the ω phase stems from the interactions of dislocation populations during the shock. In contrast, the type of plasticity events activated prior to, during, and after phase transformation may vary according to the applied peak shock pressure and duration of shock, implying multiple dislocation interaction possibilities which result in energetically varying post-shock dislocation populations. Unfortunately we do not have a understanding of the material state and defect state during the shock event to speculate on the nature of the immediate ω→α reversion. Modeling the simultaneous reversion would require further knowledge on possible dislocation interactions at high strain rates and how their interaction products eventually pin the ω phase.  









 §.§ Distinguishing Peak Shock Effects (κ parameter)
 
Referencing figure <ref>, it was noted that the model was able to generally capture the ω phase evolution in 10.5 GPa samples despite overpredicting the reduction in microstrain. The simulation results were obtained by directly using parameters calibrated to the 8 GPa samples. We also conjectured that this observational mismatch could be attributed to the difference in defect state in the 10.5 GPa samples. In order to justify this claim, we revisit the VISAR profiles of the shock deformation in figure <ref>, and note that samples shocked to 10.5 GPa complete the α→ω transformation in a shorter period of time compared to the samples shocked to 8 GPa. This transformation period is then followed by a period of constant velocity,  indicating continued deformation of the sample at peak pressure during this period. Due to the earlier completion of the transformation, the ω phase in the 10.5 GPa samples undergoes a longer period under deformation which consequently affords more time for dislocation interactions. Additionally, the exaggerated pressure loading could also activate different slip modes and act to further diversify the interaction possibilities. Ultimately, this could imply that the ω phase dislocation populations may be distinct, based on the duration of deformation following the α→ω transformation. With the exception of the 443K sample, all of the 10.6 GPa samples exhibit the same reduction of microstructure during annealing. This leads us to speculate that there are two distinct populations of defects present: i) the arresting dislocations which are recoverable through the climb and glide processes in our model and ii) dislocations that are not readily removed at these low temperatures but do not exert a significant resisting force on the transformation. 

During the early stages of shock the α phase material undergoes plastic deformation by slip and twinning <cit.>, generating a high dislocation density. Following the shear transformation from α→ω a portion of these formerly mobile α dislocations may survive in the ω in sessile high-energy configurations. We hypothesis that these defects strongly contribute to the arresting force. Given the high energy configuration and high density of these defects, they will rapidly annihilated when they are mobilized. In contrast, we speculate that the second group of dislocations are native to the ω slip systems and were generated during the longer deformation period during the 10.5 GPa shock.  Kumar and Beyerlein recently computed the generalized stacking fault surfaces for ω-Zr at room temperatures and reported that the favorable slip modes for the ω phase consist of prismatic <c>, prismatic-II <101̅0> and pyramidal-II <c+a>, which are different from the conventional slip modes in the α phase <cit.>. As more experimental and computational studies aimed at exploring the fundamental deformation mechanisms of the ω are performed, we will hopefully be able to verify this hypothesis either by reanalysis of the Low et al. data or by additional in-situ annealing studies. 



In section <ref> we also illustrated the possibility of varying the κ parameter to achieve better fits in the 10.5 GPa ω phase fractions, as shown in figure <ref>. From a physical understanding, we equated the act of decreasing the κ parameter to decreasing the jog spacing in screw dislocations that are removed in a glide controlled manner. Decreasing κ effectively describes a more deformed system, as one would expect smaller jog separations in a highly deformed state. From equations <ref> and <ref>, we can see that decreasing the value of κ for the 10.5 GPa sample would also decrease the dislocation removal rate, given a constant temperature. One would then logically expect a consequent decrease in the reverse transformation. However, this is not the case as we see in figure <ref>, in which a lower κ value results in a higher initial transformation rate. In order to explain this, a more detailed look at the effect of κ on the dislocation evolution and its subsequent effect on the reverse transformation is required. This requires that, we first examine (i) the barrier energy coupling both processes and (ii) the energetic implications of varying κ. 



  §.§.§ The Barrier Energy

In previous work, it was observed that the dislocation density in the ω phase decreased considerably before the ω→α transformation occurred in earnest <cit.>. It was thought that reduction to a certain critical value of defect level in the ω phase had to be achieved prior to initiation of significant ω→α transformation. While sufficient in explaining the rapid initial transformation, this hypothesis was inadequate in explaining the progression of the reverse transformation to new metastable levels depending on the temperature applied. Instead of a critical barrier, this model presents a continuous barrier energy that evolves in tandem with the dislocation density. This barrier energy acts to resist further reverse transformation of the ω back into α. If we associate the immobile dislocations with the arresting of the ω phase, and that the mobilization and subsequent removal of these dislocations allow for the reverse transformation, we can then alternatively describe the resistive barrier energy as being related to the energy required to mobilize these dislocations. 

This mobilizing energy can be thought to follow a distribution, such that dislocations with lower mobilizing energy requirements are favorably set in motion for elimination reactions. Subsequent removal of dislocations with lower mobilizing energy will mean that the remaining dislocation populations are those with higher mobilizing energy and consequently offer a stronger resistance to the ω→α transformation. In the current model, this mobilizing energy for individual dislocations is related to the dislocation energy in a straightforward manner given in equation <ref>. The increase of the average mobilizing energy with respect to the amount of dislocations remaining is given in figure <ref>. The observation of such an increasing trend falls in line with observed reduction trends in the microstrains (figures <ref> and <ref>), in that a steep reduction occurs in the beginning of the annealing process followed by a much slower reduction in the asymptotic region. The steep reduction in the initial stage comes from the rapid removal of large dislocation populations that are unstable due to their low mobilizing energies (glide controlled), following which they are exhausted and the dislocation removal is proceeded by the removal of more stable dislocation populations (climb controlled). At high homologous temperatures, there is sufficient energy to overcome all barriers and the transformation goes to completion. However, at lower temperatures, a balance is produced between the available thermal energy and the energy necessary to remove dislocations arresting the transformation. In the application to the transformation rate equation (eq. <ref>), we note that the barrier energy is normalized by the current dislocation density and represents an averaged resistance of a single dislocation at a point in time. Hence, the β_1 parameter is applied to scale this barrier energy to represent a collective contribution from an effective dislocation population. In the current implementation, we have chosen β_1 = 1.9× 10^14 which is ∼ 20% of the reference dislocation density ρ_ref=1 × 10^15. 






  §.§.§ Varying the κ Parameter

From figure <ref>, decreasing the value of κ was shown to increase the initial transformation rate, although inspection of equations <ref> and <ref> indicate that it should lower the removal rate of glide dislocation populations and subsequently reduce the initial transformation rate. The explanation for this observation is that the behavior of the reverse transformation is not purely reliant on dislocation removal, but rather the energetics of the removal process. In the current model, the barrier energy resisting the ω→α transformation is described as being related to the mobilizing energy requirement of the dislocation populations. By decreasing the κ parameter, we are also effectively increasing the mobilizing energy requirement for the glide dislocation populations, and consequently slowing down the dislocation recovery process. Equivalently, we are widening the energy distribution solely for the glide populations while keeping the same dislocation content so that the glide population now includes dislocations with higher energy requirements for mobilization. However, this energy increase does not shift the glide populations to be on par with the climb populations, they are still exhausted before climb controlled recovery. Since the mobilization energy requirement of a dislocation is directly related to its resistive strength against the transformation, in the case of the inevitable exhaustion of glide populations, increasing their energetic requirements increases their resistive strengths against the transformation, and that removal of a population of `stronger' resistive dislocations results in a faster initial reverse transformation. To simplify this argument further, varying κ effectively determines how strongly the glide population arrests the reverse transformation. We see this reflected in figure <ref>, in which decreasing the jog spacing κ results in a higher initial rate of transformation, with minor effects on the asymptotic trend.

In the isothermal annealing experiments, given sufficiently high temperatures (>483 K), a cross-over trend was observed when comparing the α phase fractions of both the 8 and 10.5 GPa samples. More specifically, even though the 10.5 GPa samples start off with a lower amount of α phase than the 8 GPa samples, the growth of α in the 10.5 GPa is rapid enough that it catches up and surpasses the α content in the 8 GPa samples. Figure <ref> depicts the cross-over trend observed at 483 K when comparing samples shocked different pressures. Note that the fit applied to the 10.5 GPa curve here is the same one that was used in figure <ref> with κ=0.0215. Low et al. speculated that the transformation rate was directly related to the current ω fraction so that similar cross-over behavior would be observed anytime two samples with different initial ω content were simultaneously annealed given enough time. In the model analysis, we presented an alternate view that the transformation rate was more related to the types of defects present and not the current ω fraction.  In order to highlight the difference in these views we ran an additional simulation where both the 8 and 10.5 GPa samples were assigned equal values of κ=0.023 and annealed at 483 K to longer times of t=100,000 s and found no cross-over points. Figure <ref> illustrates the results of this hypothetical simulation. Notice that at asymptotic times, both 8 and 10.5 GPa samples have achieved a plateau in the reverse transformation, and are unlikely to cross-over at further times. Again, this result points to the difference in glide controlled dislocation population as the microstructural feature which directly controls the initial transient transformation rate. This argument is further solidified by the observation that at asymptotic times, which we previously established to be dictated by climb controlled removal of dislocations, both samples behave identically and are metastable to long times with no observable trends pointing to a cross-over. 












 §.§ Nucleation vs. Interfacial Growth of α
 
As discussed earlier, the model is able to be extrapolated across a range of temperatures to reproduce Brown et al.'s annealing experiments as shown in figure <ref>. Upon close observation of the simulation prediction in figure <ref>, the volume fraction curves have abrupt changes in the slope (which are especially noticeable in the beginning of the transformation). This is related to the linear piecewise fitting that has been applied to the α_1 parameter as depicted in figure <ref>. Here, we attempt to ascribe a meaning to the α_1 parameter. Referring to figure <ref>, if we consider only the temperature range from the isothermal experiments, and proceed along the x-axis in increasing temperature, we would observe a rapid initial rise of the α_1 parameter (stage 1) followed by a very noticeable decrease in the slope at 463 K<θ<483 K (stage 2). Further increasing the temperature sees another rapid increase in the α_1 parameter (stage 3), although slightly less rapid as the initial increase. 

In Low et al. the authors compared the EBSD scans of samples in the as shocked state, post annealing at 443 K, and post annealing at 463 K. In the sample annealed at 443 K, the growth of the α phase seemed to proceed by interfacial growth of existing α lath as the number of laths in the annealed condition did not increase but the laths were visibly thicker <cit.>. In the sample that was annealed at 463 K, the volume fraction of the α phase was not only higher but the number of α laths increased and the average lath was significantly thinner suggesting nucleation and growth of new laths. Here, we suggest that the temperature trends of the α_1 parameter reflect the transition of the reverse transformation from being interfacial growth dominant to being nucleation dominant, and that this transition occurs in stage 2 of the α_1 vs. temperature fitting. Although the current model supports the hypothesis of the reverse transformation being closely related to the removal of defects, much work is left to be done in terms of understanding the mechanisms by which the reverse transformation occurs, eg. the type of defects at the interface and their removal mechanisms which lead to growth of current interfaces or nucleation of a new α lath. For the purposes of this model the α_1 parameter then acts as a modifier to the bulk Gibbs free energy difference to account for the variation in specific transformation mechanisms active at a given temperature.  



§ CONCLUSION

In prior work, in-situ XRD experiments were undertaken to investigate the thermal stability of the α/ω microstructure generated from shock loading pure Zr. Provided sufficient thermal energy, the material exhibited a rapid initial transformation rate back to α followed by a plateau to a metastable state with a lower amount of retained ω. Similar trends were observed in the dislocation densities (associated with the measured microstrains) in the ω phase. However, comparing both reductions simultaneously revealed that a significant reduction of dislocation densities in the ω phase preceded the reverse transformation, leading to postulation that the dislocation populations within the ω phase were responsible for the metastability of the ω phase at ambient conditions. Building atop this hypothesis, a model that couples the dislocation removal to the reverse transformation has been developed. Both phenomena are coupled together via a barrier energy that is a function of the defect state of the microstructure i.e. the remaining dislocation population that effectively resists the reverse transformation. The removal of dislocations favors of defects that have a lower energetic requirement for mobilization. As more dislocations with low mobilizing energies are removed, the remaining dislocation population are much more stable and harder to mobilize, effectively providing a stronger resistance to the reverse transformation, thus the occurrence of a plateau to metastable states with retained α under heating.

The dislocation removal in the ω phase has been attributed to two distinct modes of occurrence, namely the glide and climb controlled process. Glide controlled removal of dislocations has been demonstrated to dictate or aid significantly in the initial, rapid reduction of dislocations. Upon exhaustion of dislocations population associated with glide elimination, the asymptotic rate removal then becomes climb controlled. Due to the coupling between dislocation reduction and the reverse transformation, the association of the dislocation reduction with two simultaneously occurring modes (which differ in contribution strength depending on the annealing stage) allows for a better and sensible description of the transient and asymptotic plateauing of the reverse transformation.

The model also adequately describes the effects of peak shock pressures (difference in extent of deformation) on the annealing responses. This is done via the κ parameter from equation <ref> which relates the jog separation on screw dislocations to the average dislocation separation distance. A lower κ value reflects a more deformed state of dislocations. Doing so allowed for verification of cross-over trends observed in annealing the experiments, in which the 10.5 GPa samples ended up with lower amounts of retained ω compared to 8 GPa samples under similar isothermal environments, despite the 10.5 GPa samples starting off with higher amounts of ω. 


















































































































































































§ REFERENCES




