Ab Initio Simulations of Temperature Dependent Phase Stability and Martensitic Transitions in NiTi Justin B. Haskins,†,§ Alexander E. Thompson,‡,§ and John W. Lawson∗,¶ AMA Inc., Thermal Materials Protection Branch, NASA Ames Research Center, Moffett 1 Field, California 94035, USA, USRA, Thermal Materials Protection Branch, NASA Ames Research Center, Moffett Field, California 94035, USA, and Thermal Materials Protection Branch, NASA Ames Research Center, Moffett Field, California 94035, USA E-mail: [email protected] ∗To whom correspondence should be addressed †AMA Inc., Thermal Materials Protection Branch, NASA Ames Research Center, Moffett Field, Califor- nia 94035, USA ‡USRA, Thermal Materials Protection Branch, NASA Ames Research Center, Moffett Field, California 94035, USA ¶Thermal Materials Protection Branch, NASA Ames Research Center, Moffett Field, California 94035, USA §Contributed equally to this work 1 Abstract 2 For NiTi based alloys, the shape memory effect is governed by a transition from a 3 low-temperature martensite phase to a high-temperature austenite phase. Despite con- 4 siderable experimental and computational work, basic questions regarding the stability 5 of the phases and the martensitic phase transition remain unclear even for the simple 6 case of binary, equiatomic NiTi. We perform ab initio molecular dynamics simulations 7 to describe the temperature-dependent behavior of NiTi and resolve several of these 8 outstanding issues. Structural correlation functions and finite temperature phonon 9 spectra are evaluated to determine phase stability. In particular, we show that finite 10 temperature, entropic effects stabilize the experimentally observed martensite (B19’) 11 and austenite (B2) phases while destabilizing the theoretically predicted (B33) phase. 12 Free energy computations based on ab initio thermodynamic integration confirm these 13 results and permit estimates of the transition temperature between the phases. In 14 addition to the martensitic phase transition, we predict a new transition between the 15 B33 and B19’ phases. The role of defects in suppressing these phase transformations is 16 discussed. 17 2 Shape memory alloys are materials that after deformation recover their original shape 18 upon heating. They are technologically important for a wide range of applications, including 19 actuators, shape-morphing wings and next generation space suits, among others. Nickel 20 Titanium (nitinol) is perhaps the best-known example in this class of alloys and figures 21 prominently in many commercial applications. The shape memory effect in NiTi is driven 22 by a martensitic phase transition from a low temperature martensite phase (B19’) to a high 23 temperature austenite phase (B2).1,2 Many applications involving shape memory alloys are 24 tied to the specific value of the martensitic phase transition temperature.3 Having the ability 25 to tune this transition temperature, for example, through ternary additions in NiTi-based 26 alloys,3–5 will open the door to significantly more far-reaching applications. However, even 27 for binary, equiatomic NiTi, which is the simplest example in this class of materials, basic 28 questions regarding the stability of the phases and the martensitic phase transition remain 29 unclear. In this paper, we resolve several of these important, outstanding issues. 30 Experimentally, the high temperature austenite phase of NiTi has the cubic B2 (Pm¯3m 31 symmetry) structure. The low temperature martensite phase has the monoclinic B19’ 32 (P2 /m symmetry) structure, with an experimentally determined angle γ of 98◦.6,7 The 33 1 transition temperature between the two phases is reported to be approximately 341 K.8 Rel- 34 evant crystal structures are shown in Figure 1. Considerable computational work has been 35 performed to understand the phases of NiTi and related materials. In particular, density 36 functional theory (DFT) studies9–22 have provided many insights into the energetics and 37 properties of NiTi; but they have also generated new unanswered questions. For example, 38 DFT formation energies for B2 are in good agreement with experiments;23–26 however, B2 is 39 predicted to be dynamically unstable at T = 0, i.e. certain phonons modes have imaginary 40 frequencies.27 Recent attempts using small systems to include finite temperatures effects into 41 B2 stability analyses have given contradictory results.28,29 On the other hand, B19’ at the ex- 42 perimental monoclinic angle γ of 98◦ is dynamically stable at T = 0; however, the computed 43 structure is unstable to shear.12 Huang et al. determined the DFT ground state of NiTi 44 3 at T = 0 to be a new orthorhombic phase (B33) with an angle of γ = 107.3◦ .30 However, 45 the B33 structure has not been observed experimentally in NiTi and its crystal symmetry 46 (Cmcm) is incompatible with the shape memory effect, and therefore cannot represent the 47 martensitic phase of this material. Thus, after considerable computational analysis, we are 48 in the unsatisfying position that the two experimentally observed phases for NiTi have un- 49 determined stability; whereas the only computed stable phase has never been observed and 50 is incompatible with the shape memory effect. 51 To address these discrepancies, we perform high accuracy, ab initio molecular dynamics 52 (AIMD) simulations based on density functional theory combined with extended thermody- 53 namic integration methods to evaluate the stability and relative free energies for the defect- 54 free, single crystal phases (B2, B19’, B33) of NiTi for a range of temperatures up to 900 K. 55 These materials are strongly anharmonic, and therefore, methods based primarily on phonon 56 analysis, even at finite temperatures, will not capture the full behavior. This necessitates 57 high accuracy computations of the free energy. We show that finite temperature, entropic 58 effects resolve many of the controversies derived from previous studies, bringing computa- 59 tion into much closer agreement with experiment. Specifically, we show that entropic effects 60 stabilize both B2 and B19’ while destabilizing B33. Furthermore, the martensitic transition 61 temperature is estimated between these stable phases. In addition, we also identify a new 62 phase transition between B33 and B19’. 63 For B2, B19’, and B33, we consider phase stability from several complementary view- 64 points. For each case, the lattice vectors of the AIMD simulation cells are optimized such 65 that all finite temperature components of the stress tensor are zero. This procedure not only 66 accounts for thermal expansion, but also places the system at a critical point on the free 67 energy surface. Next, we examine deviations of the crystalline structure from ideality during 68 the course of the simulations in these optimized cells. Structural evolution is evaluated quan- 69 titatively with: (1) normalized position correlation functions (NPCFs) 31 and (2) atomic dis- 70 placementscatterdiagrams.32 TheNPCFisproportionalto(cid:80) (cid:104)(r (t−t )−R0)·(r (t )−R0)(cid:105) 71 i i 0 i i 0 i 4 where r (t) are the atomic trajectories from the AIMD simulation, R0 is the ideal reference 72 i i lattice vectors of interest and the brackets are ensemble averages. For long times (t → ∞), 73 vibrational motion becomes uncorrelated, and therefore, NPCF → 0 indicates stabilization 74 with respect to the reference lattice whereas nonzero values indicate the converse. We plot 75 the atomic displacements explicitly relative to the reference structures on scatter diagrams. 76 Significant deviations from zero displacement signal an instability. 77 Both NPCFs as well as atomic displacement scatter plots are shown in Figure 2 for 78 144 atom cells of B33, B19’, and B2. Large cell sizes are required to eliminate finite size 79 effects (see Supplemental Documentation for extensive convergence tests for all computed 80 properties). Figure 2a, b, and c, shows very different behavior for the three phase at different 81 temperatures. For B33, the NPCFs indicate structural stability at lower temperatures, 50 K 82 and 300 K, but instability for T > 300 K. Convergence times at 300 K are almost two 83 orders of magnitude larger than at 50 K. This may indicate the proximity of a stability 84 transition for B33. Interestingly, the B19’ phase loses its T = 0 non-zero shear stress even 85 at low T and maintains its ideal configuration across the full 50 to 600 K temperature range 86 considered. Perhaps most striking is that while the B2 structure is unstable at 50 K, it 87 stabilizes for T > 300 K. Unlike B33, the NPCF convergence rates for B2 increase with 88 increasing temperature. 89 The atomic displacements scatter plots for B2, B19’, and B33 are shown in the overlay 90 plots in Figure 2a, b, and c, respectively. For each case, displacements are provided for 91 50 K (blue circles) and 600 K (red squares). At low temperatures, B33 displacements are 92 negligible; however, at high temperatures, large displacements on the order of 0.5 Å can 93 be seen in the a direction. Displacements in a result from thermally induced motion along 94 the [100](011) stacking fault, which previously has been shown from DFT calculations to 95 be important for the martensitic transition.9,10 At low and high temperatures, the B19’ 96 phase exhibits only minor displacements, ∼ 0.05 Å. The B19’ displacements do not show 97 any particular ordering and can most likely be attributed to vibrational motion, as indicated 98 5 by the loss of correlation in the NPCF. The B2 phase at low temperatures shows large 99 displacements ∼ 0.4 Å from ideality. At higher temperatures, however, these displacements 100 largely vanish, as seen in the tight clustering near the origin. Both the NPCF and the atomic 101 scatter plots indicate that for T > 300 K, the high temperature phase of NiTi is very closely 102 approximated by ideal B2. 103 Phase stability is further investigated by explicit computation of temperature-dependent 104 phonons as derived from the AIMD simulations.33,34 Imaginary phonon modes (represented 105 as negative numbers) indicate the crystal structure is dynamically unstable, i.e. it is not a 106 local minima of the energy. Phonon spectra are shown in Figure 3 at both zero-temperature 107 and at 600 K. It is important to note that these results are very sensitive to cell size, and 108 therefore, using sufficiently large cells is crucial to obtain reliable results. Because of this, 109 144 atom cells were employed for computations on all phases (see Supplemental Figures S3- 110 S8 for convergence study results). The B33 phase, given in Figure 3a, develops imaginary 111 modes at 600 K through the lowering of the TA mode along the Γ → A direction. Phonon 112 dispersions for the B19’ phase shows stability across the full temperature range investigated, 113 as shown in Figure 3b. Most dramatically perhaps, the imaginary modes reported for the 114 B2 T = 0 K phonon dispersion lift and become positive at 300 K, as shown in Figure 3a, 115 indicating stabilization of this phase, consistent with the structure analysis of Figure 2. 116 The stress tensor, structure, and phonon analyses provide a complementary picture of 117 the temperature dependent stability of the three phases that is consistent with experiment. 118 Namely, stable phases at a given temperature exhibit the following properties: all compo- 119 nents of the stress tensor (normal stresses and shears) are on average zero; the NPCF goes 120 to zero in finite time; and all phonon modes are positive. Our results therefore show that 121 finite temperature, entropic effects stabilize the high-temperature B2 phase between 50 and 122 300 K. Similarly, the low-temperature B33 phase is progressively destabilized, fully losing 123 stability between 300 K and 600 K. The B19’ phase, on the other hand, is unstable to shear 124 at T = 0 but exhibits full stability from 50 K up to 700 K. 125 6 To obtain further insights into phase stability as well as transitions between the phases, 126 we next compute the relative free energies of the phases. Vibrational entropy is frequently 127 evaluated via the quasi-harmonic approximation (QHA). However, stability issues at T = 0 128 invalidates this approach for B2 due to the appearance of imaginary phonon modes. Alter- 129 natively, stable, finite temperature phonon spectra can be used with the QHA expressions 130 to obtain entropy estimates. However, for strongly anharmonic materials such as NiTi, this 131 approach will not be a good approximation to the full anharmonic entropy. This strong 132 anharmonicity therefore necessitates the use of high accuracy methods for computing free 133 energies. For this reason, we use two different methods based on thermodynamic integration 134 to compute the free energies. The computations should be exact to within the accuracy of 135 DFT. 136 Our first approach is a generalization of the stress-strain methods developed previously 137 for transition metals.35 Those methods based on Bain path integration are necessarily vol- 138 ume conserving; however, many systems of interest including NiTi do not conserve volume 139 between the phases. We generalized that approach to account for arbitrary volume changes 140 in an exact way. We expect this method to have applicability beyond what is presented 141 here. Our generalized stress-strain method requires a well defined, continuous path in lattice 142 vector space between the two given phases. For NiTi, the monoclinic angle, γ, provides a 143 natural, continuous parameter to connect the three phases of interest shown in Figure 1. 144 In general, multiple paths can be considered; however, the B33→B19’→B2 path was deter- 145 mined to be the best behaved and is equivalent to motion along the (cid:104)100(cid:105){110} generalized 146 stacking faults. Spontaneous motion along this fault was found in the high temperature B33 147 phase during structural stability tests. The B33→B19’ path is largely a transformation in γ- 148 space, as the lattice vectors are of comparable magnitude, while the B19’→B2 path involves 149 non-trivial changes to both γ as well as the lattice vectors. Optimization of the simulation 150 cells to obtain zero stress is required to ensure the obtained free energy differences, which 151 are Helmholtz free energy differences, are equivalent to Gibbs free energy differences. We 152 7 find the internal atomic coordinates for the 144 atom cell for this path to transform contin- 153 uously and that the stresses converge rapidly (< 10 ps) (see Supplemental Figure S9). It 154 should be noted that while free energy differences between stable phases can be rigorously 155 computed, evaluation of free energies differences involving unstable structures is still an area 156 of active investigation. Therefore, free energies involving unstable structures may contain 157 some systematic error as discussed recently.36 158 Our second approach uses the Einstein crystal method to compute free energy differences 159 atisolatedpointsalongthetransformationpaths. Thesecomputationswereusedtocheckthe 160 stress-strainmethodandwereonlyperformedatfreeenergymiminaalongthetransformation 161 path. The reference harmonic free energy is obtained from the force constants associated 162 with the temperature dependent phonon dispersions. This approach overcomes difficulties 163 in using T = 0 K phonon dispersions with imaginary modes. For each stable crystal at a 164 given temperature, thermodynamic integration is performed from the system described by 165 the harmonic reference potential to the one described by DFT. The integration can result in 166 anharmonic contributions to the free energy on the order of 5 meV/atom compared to the 167 harmonic reference free energy. This nontrivial anharmonic contribution to the free energy 168 can shift the transition temperature by as much as 100 K and thus confirms the need for 169 high accuracy or exact methods to study these systems. 170 Free energy results using the generalized stress-strain method at 0, 50, 300, and 600 K 171 are given in Figure 4a as a function of γ. Einstein crystal results are shown as open symbols 172 for validation. Agreement between the methods is excellent (≤ 1 meV/atom). The T = 0 173 curve reproduces previous DFT results, and clearly shows that B2 and B19’ are not energetic 174 minima whereas B33 is a stable minimum, as reported by Huang et al.30 We see however 175 that the free energy surface changes considerably as a function of temperature. Between 176 B33 and B19’, a small but distinguishable barrier develops between the phases for T = 177 50K −300K. Above 600 K, however, the free energy is monotonically decreasing from B33 178 to B19’. Importantly, B19’ develops a clearly defined free energy minima above 50 K. Thus, 179 8 B19’ is entropically stabilized and develops into a separate phase distinct from B33. The 180 B2 phase is unstable to transitions to B19’ until 300 K, above which a free energy barrier 181 develops stabilizing this phase as a local minima. These results are consistent with the 182 structural and phonon analysis. 183 Further detail is provided by Figure 4b where free energy differences relative to the most 184 stable phase are mapped as a function of T and γ. Blue and red represent small and large 185 free energy differences, respectively. The free energies are again derived from the generalized 186 stress-strain method. The free energy map illustrates the stability regions associated with 187 each of the phases: B33, B19’, and B2. White circles indicate stable points of each phase, 188 i.e. all finite temperature stresses are zero and all finite temperature phonons are real and 189 positive. Thus each white circle represents a stable, free energy miminum for that phase and 190 the white lines denote the extent of the stable free energy basins. Free energy and stability 191 results are provided for 50 K as well as between 0 and 900 K in steps of 100 K. The regions 192 of stability for each phase are found to be 0 < T < 300 K for B33, 50 < T < 600 K for B19’, 193 andT > 300KforB2. Interestingly, theB19’angleisshowntobeafunctionoftemperature, 194 ranging from ∼ 100◦ at 50 K to ∼ 98◦ at 600 K. Furthermore, the stable basin of B19’ is 195 relatively shallow suggesting that γ values for this phase might be fairly sensitive to small 196 changes in stress. This could be important since stress fields associated with defects could 197 potentially alter the value of γ quoted here. The free energy results allow us to estimate the 198 phase transition temperatures. In particular from Figure 4, the transition between B33 and 199 B19’ appears to occur between 50 and 300 K. The B19’→B2 free energy path, is uphill until 200 600 K, indicating that a transition occurs between 300 and 600 K. 201 DifferencesintheGibb’sfreeenergy(∆G)betweenthevariousstablephasesasafunction 202 of T are plotted in Figure 5. Vanishing of the free energy difference indicates a phase 203 transition. For ∆G between B33 and B19’, a new phase transition is predicted to ocurrs 204 at 75±26 K. For B19’ to B2, ∆G goes to zero at 500±14 K. The larger error for the low- 205 T transition is a function of the slope of ∆G and the target accuracy of 1 meV/atom. 206 9 The low values for the B33→B19’ transition temperature explains the lack of experimental 207 evidence for B33, despite being energetically favored at 0 K. If B33 becomes unstable at 208 low temperatures, it may be difficult to synthesize and therefore to observe. The B19’→B2 209 transition temperature is roughly 150 K above the experimental value of 341 K.8 However, 210 the methods used to obtain this value, based on ab initio thermodynamic integration, are 211 numerically exact to within the accuracy of DFT. Therefore, we expect it to be a reliable 212 estimate of the martensitic transition temperature for defect-free, single crystal NiTi. The 213 non-trivial difference with experimental values is most likely due to defects that have been 214 shown to suppress transition temperatures in this and related materials.37 This also suggests 215 that improved processing resulting in higher material quality could produce materials with 216 higher measured transition temperatures. In addition, it is also known that the transition 217 temperature is dependent on the heating and cooling rates with slower rates giving higher 218 transition temperatures. Since we use equilibrium methods to estimate this temperature, 219 our results correspond effectively to infinitely slow rates. For that reason, we expect them 220 to be an upper bound for the experimental transition temperature. 221 We have performed a comprehensive computational analysis based on ab initio molecular 222 dynamics of the stability and transitions between the major phases of NiTi: B2, B19’, and 223 B33. Considerable previous computational analysis based mainly on T = 0 DFT resulted 224 in significant discrepancies between experiment and computation. We have shown that by 225 including temperature dependent entropic effects into the computations, many of these dif- 226 ferences can be resolved. We show that B2 and B19’ are stabilized due to these entropic 227 effects whereas B33 is destabilized. These materials are shown to be highly anharmonic. An- 228 harmonic contributions to the free energy can shift the transition temperature by as much 229 as 100K and thus necessitates the need for high accuracy or exact methods to study these 230 systems. Wedevelopangeneralizedstress-strainmethodtoperformsuchcomputations. The 231 phase transition temperature between B2 and B19’ is estimated to be approximately 500 K 232 for defect-free, single crystals which is about 150 K above experimental results. Defects and 233 10