First-principles phonon-based model and theory of martensitic phase transformation in NiTi shape memory alloy

First-principles Phonon-based Model and Theory of Martensitic Phase Transformation in NiTi Shape Memory Alloy

Journal Pre-proof

First-principles Phonon-based Model and Theory of Martensitic Phase Transformation in NiTi Shape Memory Alloy Pawan Kumar, Umesh V. Waghmare PII: DOI: Reference:

S2589-1529(20)30019-3 https://doi.org/10.1016/j.mtla.2020.100602 MTLA 100602

To appear in:

Materialia

Received date: Accepted date:

23 January 2020 25 January 2020

Please cite this article as: Pawan Kumar, Umesh V. Waghmare, First-principles Phonon-based Model and Theory of Martensitic Phase Transformation in NiTi Shape Memory Alloy, Materialia (2020), doi: https://doi.org/10.1016/j.mtla.2020.100602

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

First-principles Phonon-based Model and Theory of Martensitic Phase Transformation in NiTi Shape Memory Alloy Pawan Kumar Theoretical Sciences Unit and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore 560064 India

Umesh V. Waghmare∗ Theoretical Sciences Unit and School of Advanced Materials, Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore 560064 India

Abstract Remarkable shape memory property of Nitinol (NiTi) originates from a diffusionfree reversible phase transition it exhibits upon cooling from the cubic austenite to a low-symmetry martensite structure. Atomistic mechanisms of this martensitic transformation (MT) and consequent emergence of its various lowsymmetry structures are not understood yet. Starting with first-principles density functional theoretical calculations, we present here a phonon-based model and its statistical mechanical analysis to obtain atomistic insights into martensitic phases and transformation in NiTi. We uncover seven order parameters that are relevant to the MT in NiTi. From Monte Carlo simulations of an effective model Hamiltonian derived to capture its low energy landscape, we determine its soft phonons and establish the cell-doubling M50 mode as the primary order parameter. Using Landau theoretical analysis, we show that relative strengths of its third-order coupling with secondary order parameters (e.g. strain) determine the symmetry of low-T structures emerging at its MT. These couplings can be used as the descriptors of stability of martensitic phases ∗ Corresponding

author Email address: [email protected] (Umesh V. Waghmare)

Preprint submitted to Journal of LATEX Templates

February 4, 2020

that will guide the strategies to improve the shape memory of NiTi through substitutional alloying. Keywords: Shape Memory, NiTi, Phonon, Free energy and Martensitic Phase Transformation

1. Introduction Shape Memory Alloys (SMAs) are of great importance to technologies ranging from medical stents to smart structures in the aerospace industry. They exhibit a shape memory effect (SME), in which a highly deformed material recovers its predetermined shape upon heating [1]. Fundamental to the SME is the martensitic transformation (MT), in which crystal structure of NiTi changes from the high-temperature cubic austenite (B2) structure to a low-temperature monoclinic martensite (B190 ) structure [2] through a non-diffusive first-order phase transition involving atomic displacements (phonons) and large deforma10

tions in the crystal shape (strain) [3]. MT in SMAs has been investigated through experimental [4, 5, 6, 7, 8] and theoretical studies over the past few decades [9, 10, 11, 12, 13]. While Nitinol (NiTi) is the prototypical member of an important class of SMAs, a clear understanding of the precise mechanisms of MT and relative stability of its various low-symmetry structures relevant to its SME is still lacking. Depending on the composition of NiTi-based alloys, their MT involves five crystalline structures: B2 (cubic, space group (SG) P m¯ 3m), R (trigonal, SG P 3) [14], B19 (orthorhombic, SG P mma), B190 (monoclinic, SG P 21 /m) and BCO (base-centered orthorhombic, SG Cmcm) [15], a special case of B190 struc-

20

ture with lattice parameters satisfying certain geometric conditions. At hightemperatures, these alloys have the cubic B2 structure, which deforms spontaneously through one of the three transition pathways (Fig. 36 of Ref. [16]) to low-symmetry structures upon cooling. The transformation from B2 → B190

occurs in quenched NiTi, B2 → R → B190 occurs in N i4 T i3 and B2 → B19 → B190 occurs in ternary T i49.5 N i45.5 Cu5.0 and quaternary T i50 N i44 Cu5 Al alloys

2

[17]. Fundamental understanding of the microscopic coupling governing these transition pathways is essential to design of improved NiTi-based shape memory alloys. Several first-principles studies of NiTi [14, 15, 18, 19, 20] focused on accu30

rate prediction of structural parameters of its low-temperature phases and their electronic and vibrational properties, complementing the experimental results. In earlier theoretical works on NiTi [21, 22], B19 phase was found to be a relevant metastable phase, and it was suggested that the MT occurs in two steps B2 → B19 → B190 . However, this was suggested to be unlikely in a recent work of Strachan et al. [18], who showed that B19 phase is unstable. Seminal first-principles theoretical analysis of Huang et al. [15] showed that the BCO structure is the ground state of NiTi, while it is yet to be observed experimentally. More importantly, BCO structure can not store shape memory because it is connected to B2 structure along non-unique paths of transformation.

40

Thus, theoretical prediction of BCO as the ground state structure of NiTi had been a puzzling result. Haskins et al. [12] resolved this puzzle by determining the free energies of B190 and BCO structure using thermodynamic integration within ab initio molecular dynamics (AIMD). They confirmed that BCO is the ground state structure of NiTi, while B190 structure gains stability at T > 75 K Some of the recent theoretical works used molecular dynamics (MD) [23, 24, 25] to simulate temperature-dependent structural phase transformation in NiTi, with a semiempirical model derived from the second-moment approximation of tight-binding method. With a focus primarily on the transition temperature (TM ), their estimates of TM are in good agreement with

50

experiment. Haskins et al. [12] performed AIMD simulations to estimate free energies along the path of transformation between B2 and B190 structures. Their estimate of the transformation temperature is approximately 500 K, which is about 180 K above experimental results [26]. However, the microscopic picture and interactions governing the MT are not quite clear. The monoclinic B190 structure of NiTi stores its shape memory, as it undergoes a reversible transformation upon heating to the B2 phase involving 3

{0¯11}<011> shuffle and {100}<011> non-basal shear. In the absence of nonbasal shear, B2 transforms into the B19 structure [27, 28]. Elastic moduli of the cubic structure corresponding to {0¯ 11}<011> basal shear and {100}<011> 60

non-basal shear are C 0 =(C11 − C12 )/2 and C44 respectively, and their ratio

is called the anisotropy factor A= C44 /C 0 . Ren and Otsuka [29] presented a Landau-like model to study MT in NiTi taking into account {0¯ 11}<011> shuffle, {100}<011> non-basal shear and {0¯ 11}<011> basal shear as three order

parameters, and determined their contributions in different phases of NiTi. Otsuka et al. [16] proposed that a material with 10 < A < 20 exhibits a B2 to B19 transformation, while a material with the smaller value of A (∼ 2) exhibits B2 to B190 transformation. Anisotropy factor (A) is a physically sensible macroscopic descriptor, as its large value essentially means stabilizing B19 structure as the non-basal shear is suppressed in the low-temperature phase due to high 70

energy (∝ C44 ) cost. They also highlighted that {0¯ 11}<011> shuffle is related to softening of M-point phonons [16]. Though strain and phonon softening provide important information relevant to the MT, it is not clear (i) what are the primary and secondary order parameters, and (ii) which strain-phonon coupling drives the MT to B19 versus B190 structures. Indeed, the pathway of MT depends on a delicate competition between different microscoping couplings between phonons and strain, and a quantitatively accurate microscopic model and its theoretical analysis are needed to understand the physics of MT. Here, we present a detailed ab initio and Landau theoretical analysis of the MT in NiTi with a focus to uncover microscopic mechanisms and couplings.

80

While the martensitic transformation involves twin and planar fault structures deviating from the single crystalline phases, our goal here is to analyze the single crystal structural transformation of the B19 phase to B2 phase with temperature. In Sec. 2, we present details of the methodology used in our first-principles calculations. In Sec. 3, we present T-dependent vibrational free energies of B19, B190 and BCO structures, adding to the resolution of puzzling issue of their relative stability of B190 and BCO structures of NiTi. In Sec. 4, we uncover the relevant order parameters of MT and their relationship with phonons and 4

strain. In Sec. 5, we present an effective Hamiltonian (Hef f ) with parameters determined to capture the low energy landscape of various NiTi structures. In 90

Sec. 6, we present results of Monte Carlo simulations of Hef f to determine temperature and pressure dependent MT, uncovering the soft mode behavior. We present Landau theoretical analysis in Sec. 7 to identify specific couplings in Hef f that stabilize B19 and B190 phases. In Sec. 8, we present results of structural disorder and microstructures relevant to MT in NiTi using Hef f analysis, and finally we summarize our work in Sec. 9.

2. Methods Our first-principles calculations within the density functional theory (DFT) were performed with a plane-wave pseudopotential scheme implemented in the Quantum Espresso (QE) package [30], with a generalized gradient approx100

imation (GGA) and Perdew-Burke-Ernzerhof (PBE) [31] form of exchangecorrelation energy functional. In self-consistent Kohn-Sham (KS) calculations with primitive unit cells, the Brillouin Zone (BZ) integrations were sampled on uniform meshes of 18 × 18 × 18 k-meshes for B2 structure, and 18 × 12 × 12 k-

meshes for B19, B190 and BCO structures, and note that these uniform k-meshes include Γ-point. To check the accuracy of our quantum-mechanical calculations, specially the energy of B1900 (at γ =∼ 101o , a monoclinic structure intermediate to B190 and BCO [18]) relative to B190 (at γ =∼ 97.8o ) and BCO structures (which is within ∼1 meV), we tested the convergence of total energies with respect 110

to energy cutoffs used to truncate the plane-wave basis representing the KohnSham wave functions (Fig. S1). While total energies of these structures shifted with increasing the energy cutoff and converged at 100 Ry (Fig. S1), their relative energies do not change above the cutoff energy of 40 Ry. We thus performed the structural optimization calculations at energy cutoffs of 40 and 100 Ry. We relaxed B190 structure to minimize energy with respect to atomic posi-

5

tions keeping lattice parameters fixed at these in Ref.[15] until the force on each atom is less than 2 meV/˚ A. At these lattice parameters of B190 , shear stresses (σxy = −10.5 kbar) are a bit high, and the structure is thus not completely op120

timized. When B190 structure is relaxed with respect to both lattice parameters and atomic positions until the magnitude of stresses are less than 0.1 kbar (and force on each atom is less than 2 meV/˚ A), it transforms into B1900 structure. Thus, it is clear that B190 at experimental lattice parameters is unstable with respect to B1900 (Table S1). Optimization of the BCO structure with respect to its atomic positions and lattice parameters relaxed until the magnitude of stresses are less than 0.1 kbar (and force on each atom is less than 1 meV/˚ A), and it is lower in energy with respect to B1900 structure (Table S1). We reconfirm structural parameters and relative energies of the structures of NiTi using (i) local density approximation (LDA) as parametrized by Perdew-

130

Zunger [32] pseudopotential implemented in QE at energy cutoffs of 50 Ry and 400 Ry to truncate the plane-wave basis set for representing Kohn-Sham wave functions and charge density respectively at the uniform mesh of 20 × 16 × 16 k-points, and (ii) SCAN meta-GGA exchange-correlation functional [33, 34] as implemented in the VASP code [35, 36] along with projector augmented wave (PAW) potentials [37]. Our estimates of structural parameters and energies of the B2, B19, B190 , B1900 and BCO phases of NiTi agree well with earlier theoretical works [18, 38, 39, 40], structural parameters of B2 and B190 structures agree well with experiments [4, 6] as well, and reconfirm that BCO is the most stable of the structures of NiTi (see Table 1, and Table S1 in SI), though their relative

140

energies depend on the flavor choice of (exchange-correlation energy functional) of DFT. We also notice that differences between lattice parameters and relative energies of these phases estimated with calculations at energy cutoffs of 40 and 100 Ry are very small. Therefore, the rest of our first-principles calculations have been performed at 40 Ry energy cutoff on plane wave basis of Kohn-Sham wavefunction. We refer to B1900 as B190 in our rest of analysis, as both the structures are monoclinic and exhibit SME [18]. We determined dynamical matrices and phonon spectra within the frame6

Table 1: Structural parameters (4 atoms per unit cell) and energies of B2, B19, B190 and BCO phases of NiTi relative to B190 structure as a reference, obtained with different computational schemes. DFT estimates of the lattice parameters and energies of equiatomic structures at zero temperature and pressure, while Ref. [4] and Ref. [6] used T i49.5 N i50.5 and Ti-49.2 at.%Ni alloys, respectively. Phase

Space group

Method

a(˚ A)

b(˚ A)

c(˚ A)

γo

P m¯3m

E − EB190 (meV/f.u)

B2

PBE-GGA

3.004

4.248

4.248

90

75.34

B19

B190

BCO

Pmma

P 21 /m

Cmcm

PZ-LDA

2.937

4.153

4.153

90

112.13

SCAN metaGGA

2.966

4.195

4.195

90

156.67 84.00

USPP-GGA [15]

3.009

4.255

4.255

90

Exp [4]

3.015

4.264

4.264

90

PBE-GGA

2.832

4.589

4.168

90

28.34

PZ-LDA

2.644

4.553

4.165

90

19.78

SCAN metaGGA

2.805

4.565

4.097

90

45.86

USPP-GGA [15]

2.776

4.631

4.221

90

24.00

PBE-GGA

2.935

4.733

4.027

100.9

0.00

PZ-LDA

2.846

4.722

3.932

103.8

0.00

SCAN metaGGA

2.907

4.726

3.953

101.9

0.00

USPP-GGA [15]

2.929

4.686

4.048

97.8

0.00

Exp [6]

2.898

4.646

4.108

97.8

PBE-GGA

2.928

4.910

4.006

107.1

-1.35

PZ-LDA

2.851

4.813

3.923

107.1

-2.23

SCAN metaGGA

2.891

4.888

3.937

107.3

-0.40

USPP-GGA [15]

2.940

4.936

3.997

107.0

-16.00

work of DFT linear response (DFT-LR) at q-points on a 6 × 6 × 6 mesh for B2,

and 2 × 2 × 2 mesh in the BZ for B19, B190 and BCO structures implemented in 150

QE. To estimate the contributions to vibrational free energy of stable phonons, we use harmonic approximation and treat them as quantum oscillators:

ph Fstable (S) =

kB T X hωqν ¯ ln[2sinh( )], Nq q,ν 2kB T

(1)

ph where Fstable (S) vibrational free energy of stable phonons of the optimized

7

structure S (=B19, B190 and BCO). ωq,ν is the frequency of ν phonon mode at q wave vector, and Nq is the number of q wave vectors. Since equation (1) does not hold for unstable phonons (ω 2 < 0), we determine their contributions to vibrational free energies by treating them as quantum anharmonic oscillators (QAO). We include fourth-order anharmonic terms by fitting to double-well energy surface obtained from first-principles, and determined energy spectra of 160

the QAO’s using higher-order finite difference method [41]. We estimate the vibrational free energy of unstable modes using:

ph Funstable (S) =

h X X 1 X [ Eqν0 − kB T ln{1 + e−(Eqνl −Eqν0 )/kB T }, Nq q,ν qν

(2)

l=1

ph where Funstable (S) is the contribution of unstale phonon modes to the vibra-

tional free energy of structure S (=B19, B190 ). Eqν0 and Eqνl are the ground state and lth excited state energies of the QAO representing phonon mode ν at q wave vector. Details of this scheme presented in supplementary information (Sec S3 in SI). Using equation (1) and (2), we estimate the total free energy as follows.

ph ph G(S) = E(S) + Fstable (S) + Funstable (S),

170

(3)

where E(S) and G(S) are the total and free energies of the optimized structures S, respectively. 3. Addressing the puzzle of relative stability of B190 and BCO structures We note that the puzzling issue of BCO as the ground state structure of NiTi has been resolved by Haskins et al. [12], who obtained free energies of 8

these structures from AIMD and showing that the B190 structure stabilizes above T > 75 K. We note that their analysis is based on classical statistical mechanics, and we complement it here with quantum statistical mechanics, in which Bose-Einstein statistics of phonons is used. We use equation (3) to 180

determine free energies and relative thermodynamic stability of B2, B19, B190 and BCO structures at low-temperature. We determined their phonon spectra using DFT linear response (DFT-LR), and confirm the results of Ref. [10, 12, 13, 14] that M50 phonon mode is the strongest lattice instability of B2 structure (ω ∼57i cm−1 , Fig. 1a). In B19 structure, one of the transverse acoustic (TA) phonon branches is weakly unstable (ω ∼33i cm−1 ) along Γ → Y direction (Fig. 1b), little different from the results of Ref. [14] which do not show any imaginary frequency along Γ → Y direction. A TA phonon branch in Γ → E

direction of B190 structure exhibits a weak instability (ω ∼8i cm−1 , Fig. 1c)

in our analysis, slight different from the results of Ref. [14, 12, 13] (where no 190

imaginary frequency phonons found along Γ → E direction). We note that

unstable modes are expected in the phonon spectra of B19 and B190 as they are consistent with this instability (against shear strain) of these structures with respect to BCO structure [15]. Our calculated phonon dispersion of the BCO structure, the ground state, exhibits no instability anywhere in the BZ (Fig. 1c), and is in good agreement with Ref. [12, 13]. We now estimate the free energies of B19, B190 and BCO structures using equation (3), and we find ∆G = GB19 − GB190 > 18 meV/f.u. at T up to

temperatures 400 K, and thus conclude that B19 structure is unlikely to occur at low-T. At T = 0 K, ∆G = GBCO − GB190 = -0.10 meV/f.u. (while this 200

is within the error of our computations, we have tested its convergence with respect to q-mesh), revealing nearly equal stability of BCO and B190 structures within our error bars. However, at T > 43 (46) K B190 structure becomes more stable than the BCO structure with (without) including contributions of unstable modes to free energy of B190 structure (Fig. 1d). We note that with inclusion of contributions of unstable modes to free energy leads to stabilization of B190 structure at a slightly lower temperature (Fig. S3 in SI), because weak 9

instabilities in the TA branch of B190 structure contribute strongly to vibrational entropy, stabilizing B190 structure at an even lower temperature (43 K). The reduction in temperature is small because there are very few unstable modes. We confirm the transformation from BCO to B190 structure as it was studied earlier using AIMD [12], and quasi-harmonic approximation (QHA) [13]. Our estimate of transformation temperature from BCO to B190 is lower than that in Ref. [12] because of the quantum approach used here and classical approach used in their work (zero-point motion typically suppress the ordering at low-T), and Ref. [13] (at T = 100 K). However, the BCO phase deforms readily and transforms to B190 phase in the response to stresses, and shape memory effect is thus strongly influenced by stresses and planar twining faults. B2

75

[010]

-1

Frequency (cm )

150

4’

5

2’

1

0

25’

R

[-1 1 1]

X Γ Y

M (d)

B19’ BCO

225 150 75 0 B

100

0

Γ

M

300

Γ

Y

A

Γ

B19

200

15

5’ X

300

∆G (meV/f.u.)

-1

15 15

5’

-75 Γ (c)

(b)

5’

4’

225

-1

Frequency (cm )

(a)300

Frequency (cm )

210

E

Z

Γ

Z Γ

S

U

Γ

T R

Γ

24 18

0.1

12

0

6

-0.1

0 0

∆G = GB19-GB19’ ∆G = GBCO-GB19’

43 K

0

15

100

30

45

60

200 T (K)

300

400

Figure 1: (color online). Phonon spectra along high-symmetry lines of B2 (a), B19 (b), B190 (black) and BCO (red) (c) structures of NiTi obtained using DFT-LR calculations. Unstable modes with imaginary frequencies (ω 2 < 0) are shown with negative values. Difference between Gibbs free energies of B19 and B190 structures (black line) and BCO and B190 structures (red line) (d).

10

4. Identification of order parameters We now identify the order parameters of the MT in NiTi by deriving atomic 220

positions and lattice vectors of its low-symmetry B19 and B190 structures as changes or distortions of the reference B2 structure expressed in terms of its phonons modes and strains respectively. While M50 , M20 and M40 phonons are involved in the transformation of B2 to B190 (Fig. 2a-2c, 2e), M50 phonon alone drives the B2 to B19 structural transformation (Fig. 2a, 2d). Dominant structural instability of M50 mode makes it the order parameter common to both B2 → B19 and B2 → B190 phase transformations. The quantitative estimation

of phonons and strains involved in transformations to B19 and B190 structures of NiTi are listed in Table 2. While we find that {0¯ 11}<011> basal shear does not constitute an order parameter of any of these structures, it can arise when the 230

transformation occurs under applied stress as seen in earlier reports [27, 28]. It is evident from our analysis that Bain strain, pure shear and hydrostatic strains are involved in B2 to B19 (and B190 ) transformations. In B2 to B19 transformation, we find only four order parameters (M50 , s2 , s3 and s4 ), while in B2 to B190 seven order parameters (M50 , M20 , M40 , s2 , s3 , s4 and s5 ) are involved. Thus the couplings of {100}<011> non-basal shear deformation with M20 and M40 phonons are expected to be relevant to the transformation to B190

structure. The dependence of secondary order parameters on the primary order parameter is nonlinear (See Fig. S4); hence we can not combine them into one as a linear combination.

240

5. Construction of effective Hamiltonian Having identified the relevant phonons and strains, we now construct an effective Hamiltonian (Hef f ) that models the low-energy landscape of NiTi with special attention to these phonons, strains and their interactions. Since the BCO structure occurs in a vanishingly small subspace of configurations given by certain geometric conditions on B190 structural parameters (SI, Sec. S2), and their relative stability depends on quantum vibrational energy not captured by 11

Figure 2: (color online). Order parameters of B2 → B19 and B2 → B190 structural transformations. Atomic displacements of M50 (a), M20 (b) and M40 (c) modes at q~ = π/a(011), √ √ shown in 1 × 2 × 2 supercell of B2 structure. While M50 phonon and orthorhombic strain give B19 structure (d), M50 , M20 , M40 modes and monoclinic strain together give B190 structure (e). The planar unit shaded with light red colour at c/2 distance along [0 ¯ 1 1] direction contains one Ti (blue) and one Ni (red) atom on its edges-centers.

Monte Carlo or Molecular dynamics, we will focus our attention on the energy landscape and paths relevant to B2, B190 and B19 phases in our statistical mechanical analysis of the MT. 250

We follow here the scheme of the lattice Wannier functions (LWFs) in Ref. [42], and identify the symmetry invariant subspaces of the relevant degrees of freedom starting with full phonon dispersion of B2 structure obtained along the high-symmetry lines (Γ → X → M → Γ → R → M , see Fig. 1a) in the BZ. Subspaces of acoustic and optical phonons are separated by a gap in frequencies and involve modes dominated by Ni and Ti displacements respectively. Using symmetries and eigenvectors of the zone boundary phonons, we determine Ni-centric {~η } and Ti-centric {~τ } localized LWFs to span these subspaces of phonons. Using the cubic symmetry (of B2 phase), we express Hef f as an explicit

260

symmetry invariant Taylor expansion in ~ηi , ~τi and homogeneous strains (εαβ ),

12

Table 2: Order parameters (M-point phonon modes at q~ = π/a(011) of the cubic B2 structure and strain tensor components) associated with various low-symmetry structures of NiTi, obtained from first-principles calculations.

The amplitudes of phonon eigenmodes are in

unit of lattice constant of B2 structure, and phonon eigenmodes are expressed in terms of atomic displacements eˆph = |T ix , T iy , T iz , N ix , N iy , N iz >.

The strain eigenmodes

s = |ε1 , ε2 , ε3 , ε4 , ε5 , ε6 >, are in the Voight notation. (ε1 = εxx , ε2 = εyy , ε3 = εzz , ε4 = 2εyz , ε5 = 2εzx and ε6 = 2εxy ).

Modes

Eigenmode

Character

B2

B19

B190

M50

|0, 0.33, 0.33, 0, 0.94, 0.94 >

{0¯11}<011> Basal shuffle

0

0.075

-0.095

M20

|0, 0, 0, 1, 0, 0 >

{0¯11}<100>, Ni displacements

0

0

0.066

0

0

-0.096

{0¯11}<011> Basal shear

0

0

0

M

40

s1 s2 s3 s4 s5

|1, 0, 0, 0, 0, 0 >

√1 |0, 1, −1, 0, 0, 0 > 2 √1 | − 2, 1, 1, 0, 0, 0 > 6 √1 |1, 1, 1, 0, 0, 0 > 3

|0, 0, 0, 1, 0, 0 >

√1 |0, 0, 0, 0, 1, 1 2

>

{0¯11}<100>, Ti displacements Bain strain

0

0.072

0.047

Hydrostatic strain

0

0.003

0.017

{010}<001> Pure shear

0

0.099

0.162

{100}<011> Non-basal shear

0

0

-0.199

i indicating the lattice unit cell to which ~ηi and ~τi belong. Hef f consists of four parts:

Ni Ti Hef f = Hef ηi ) + Hef τi ) + Hspc (~ηi , ~τi , εαβ ) + Helastic (εαβ ), f (~ f (~

(4)

Ni Ti where Hef ηi ) and Hef τi ) operate in the subspaces of acoustic and optic f (~ f (~

phonons respectively. Helastic (εαβ ) is the linear and nonlinear elastic energy of homogeneous strain εαβ , and the coupling of εαβ with both sets of LWFs is included in Hspc (~ηi , ~τi , εαβ ). We note that inhomogeneous (spatially varying) α strain field is captured here by the acoustic phonons: εαβ (r) = 12 ( ∂η ∂rβ +

∂ηβ ∂rα ).

Ni Energetics of εαβ (r) captured in Hef f has to be consistent with the energies 270

of homogeneous strain εαβ captured in Helastic . This is achieved by ensuring Ni equivalence between Helastic and Hef f , Hspc in long wavelength limit.

13

Ni 5.1. Effective Hamiltonian of acoustic modes: Hef f Ni Hef ηi ) represents energetics of acoustic phonons (and inhomogeneous strain), f (~

which include all the lattice instabilities in NiTi, i.e. unstable modes (Fig. 1a), This part of effective Hamiltonian is invariant under translational and rotational symmetries. It includes harmonic and anharmonic terms, which relate to linear and nonlinear elastic energy terms in Helastic in the long-wavelength (continuum) limit. To impose the translational symmetry, we express its terms using differences in ~ηi ’s at neighboring sites and their dot products. In the harmonic 280

part, we consider differences between ~ηi ’s up to third NN sites, with a general form permitted by the symmetry of space group P m¯ 3m as follows:

N

Hhar (~η )

=

+

12 X j=1

+

6

1XX [ {A11 |~η1ij |2 + A12 (~η1ij · dˆ1j )2 } 2 i=1 j=1

8 X j=1

{A21 (~η2ij · dˆ2j )2 + A22 (~η2ij · dˆ21j )2 + A23 (~η2ij · dˆ22j )2 } {A31 |~η3ij |2 + A32 (~η3ij · dˆ3j )2 }]

(5)

where ~η1ij , ~η2ij , ~η3ij denote the LWFs differences between ~ηi at site i and its ~ηj first, second and third NNs at site j respectively. dˆ1j , dˆ2j ,dˆ3j denote unit vectors along the directions of first, second and third NN sites j respectively. dˆ21j and dˆ22j are unit vectors perpendicular to dˆ2j . Aij ’s are the harmonic Ni coefficients in Hef f , and determined from force constants of acoustic phonons at

high symmetry q-points (M and R) and ~qΣ =

π 3a (0, 1, 1).

Linear combinations of

these coefficients give the corresponding eigenvalues of acoustic phonons. Since we have six equations (Table 3) and seven coefficients, we use singular value 290

decomposition (SVD) to determine these coefficients, and have listed them in Table 6. Ni In the anharmonic part of Hef ηi ’s up to third f , we considered differences in ~

NN sites in the third-, and fourth-order terms, which are needed to include inhomogeneous strain correctly in the long-wavelength limit (interaction of the

14

N i at a q-point give Table 3: Linear combinations of coefficients in the harmonic terms in Hef f

the eigenvalue of corresponding phonon. M20 and M50 represent TA2 and doubly degenerate unstable (LA and TA1 ) modes at M-points respectively. Σ1 , Σ3 and Σ4 represent LA , TA1 and TA2 modes at Σ-point respectively, and R250 represents triply degenerate acoustic mode at R-point. Acoustic phonon

Linear combination of the coefficients

Eigenvalue of acoustic phonon

M20

4A11 + 4A21 + 4A22

9.22

M50

4A11 + 2A12 + 2A21 + 2A22 + 4A23

-4.19

(eV/f.u)

Σ1 Σ3 Σ4 R250

1 2 A11

+ 14 A12 + A21 + 14 A22 + 12 A23 + 32 A31 + A32

7.24

+ 41 A12 + 14 A21 + A22 + 12 A23 + 32 A31

0.80

+ 12 A21 + 12 A22 + 34 A23 + 23 A31 + 12 A32

1.95

1 2 A11 1 2 A11

6A11 + 2A12 + 8A31 + 83 A32

5.65

first and second NN site is zero for shear strain in third-order terms (See Table 4), and hence we need the third NN site interaction). Due to cubic symmetry, odd-order terms do not contribute to the energy of phonons at high-symmetry q-point, but do contribute in the long-wavelength limit. We approximated interaction in the third-order terms by considering the dominant term in interaction 300

with NN sites allowed by the cubic symmetry. Fourth-order interaction terms include the full cubic anisotropy for the first NNs and isotropic terms for the second and third NNs. We simplified sixth- and eighth-order terms restricting to the isotropic ones in differences between ~ηi ’s up to first NNs. Since our nonlinear elastic energy includes terms up to fourth-order (Sec. 5.3), we did not Ni include fifth- and seventh-order terms in the anharmonic part of Hef f . More-

over, the MT in NiTi is governed by unstable M-point phonons, which are not affected by the fifth- and seventh-order terms. We have tested that the errors associated with such truncation are less than 1% of the energy differences.

15

Hanh (~η )

=

N X i=1

+

6 X j=1

[B31

6 X j=1

(~η1ij · dˆ1j )3 + B32

12 X j=1

(~η2ij · dˆ2j )3 + B33

+

j=1

310

j=1

(~η3ij · dˆ3j )3

4 4 4 2 2 2 2 2 2 {B41 (η1ijx + η1ijy + η1ijz ) + B42 (η1ijx η1ijy + η1ijy η1ijz + η1ijz η1ijx )

+ B43 (~η1ij · dˆ1j )4 + B44 |~η1ij |2 (~η1ij · dˆ1j )2 } + B45 6 X

8 X

{B61 |~η1ij |6 + B81 |~η1ij |8 }]

12 X j=1

|~η2ij |4 + B46

8 X j=1

|~η3ij |4 (6)

Ni where Bij ’s are coefficients of anharmonic terms. Third-order terms of Hef f

coefficients are related to the nonlinear elastic moduli. Thus, these third-order coefficients (B31 , B32 and B33 ) are obtained from the corresponding third-order elastic moduli. Since there are only three coefficients, we use three inequivalent strain modes (hydrostatic strain, Bain strain and shear strain) to calculate these Ni coefficients (See Table 4). Fourth-order terms in Hef f contribute to energy of

phonons at both high-symmetry q-points and long-wavelength limit. We have Ni six fourth-order coefficients in the Hef f , which are calculated by freezing doubly

degenerate unstable M50 (ηy and ηz ), stable M20 (ηx ) acoustic phonons and their linear combinations (Fig. 3). We fitted the total energy of configurations 320

obtained as structural distortions with these phonon modes to a polynomial of 8th order. We used two strain modes (hydrostatic and shear strain modes) to connect with behavior in the long wavelength limit (Table 5). To determine the Ni 4 coefficients of sixth- and eight-order terms in Hef f , (128B61 = −3.7×10 eV/f.u)

and (512B81 = 8.3 × 105 eV/f.u), we fit the polynonial to double-well energy

of M50 mode. We used gnuplot software for fitting these data to polynomials (values of all these coefficients are listed in Table 6), and find that the maximum fitting errors in second-, fourth-, sixth- and eight-order coefficients are ±0.13%, ±0.31%, ±1.58% and ±3.85% respectively. These are negligible and do not affect the energy landscape and its impact on the MT.

16

M-point = π/a(0,1,1)

E (meV/f.u)

60

η=ηx η=ηy η=(ηy+ηz)/√2

40

η=(ηx+ηy)/√2

20 0 0

0.02

0.04 η (aB2 unit)

0.06

0.08

Figure 3: (color online). Total energies of cell-doubling structural distortions (~ η at q~ = π/a(0, 1, 1)). Solid lines represent the fits obtained with the 4th , 6th and 8th anharmonic parameters in Hef f .

Table 4: Strain mode is represented with s = |ε1 , ε2 , ε3 , ε4 , ε5 , ε6 >, in the Voight notation. Hydrostatic, Bain and shear strain modes are represented by shydro = ε|1, 1, 1, 0, 0, 0 >, sBain = ε| − 2, 1, 1, 0, 0, 0 > and sshear = ε|0, 0, 0, 1, 1, 1 > respectively. Coefficients of the N i are linear combinations of the 3rd order elastic moduli. 3rd order terms in Hef f

Strain mode

330

Linear combination of

Linear combination of

3rd order coefficients

3rd order elastic moduli 1 2 C111

shydro

6B31 + 33.94B32 + 41.57B33

sBain

−12B31 + 8.49B32

−C111 + 3C112 − 2C123

sshear

9.24B33

C456

+ 3C112 + C123

Ti 5.2. Effective Hamiltonian of optical phonons: Hef f Ti Hef τi ) models the energetics of optical phonons which are all stable (ω 2 > f (~

0). Hence, we include only harmonic interactions between ~τi ’s up to third nearest neighbor sites with a general form permitted by the symmetry of P m¯ 3m space group [43]:

Ti Hef τ) f (~

=

+

N 6 X 1X ˜ [A01 |~τi |2 + {A˜11 |~τ1ij |2 + A˜12 (~τ1ij · dˆ1j )2 } 2 i=1 j=1

12 X j=1

{A˜21 |~τ2ij |2 + A˜22 (~τ2ij · dˆ2j )2 } + 17

8 X j=1

A˜31 |~τ3ij |2 ]

(7)

Table 5:

The amplitudes of phonon eigenmodes are in unit of lattice constant of

B2 structure, and phonon eigenmodes are expressed in terms of atomic displacements eˆph = |T ix , T iy , T iz , N ix , N iy , N iz >.

The eigenmodes at M-point ( π (011)) are ηx = a

|0, 0, 0, 1, 0, 0 > of M20 phonon, ηy = |0, 0, 0.33, 0, 0.94, 0 > and ηz = |0, 0.33, 0, 0, 0, 0.94 > of doubly degenerate M50 mode. shydro and sshear are two strain modes. Mode

Linear combination of 4th order coefficients

4th order fitting coefficients in in polynomial (eV/f.u)

ηx

32B41 + 64B45

468.68

ηy

32B41 + 16B43 + 16B44 + 64B45

1057.35

+ ηz )

16B41 + 8B42 + 8B43 + 16B44 + 64B45

999.11

+ ηy )

16B41 + 8B42 + 4B43 + 8B44 + 64B45

√1 (ηy 2 √1 (ηx 2

shydro

6B41 + 6B43 + 6B44 + 192B45 + 648B46

sshear

0.75B41 + 0.375B42 + 30.75B45 + 168B46

637.52 1 8 C1111

+ C1112 + 34 C1122 + 32 C1123 =815.68 1 8 C4444

+ 34 C4455 =112.34

where ~τi denotes Ti-centric LWFs at site i, ~τ1ij , ~τ2ij , ~τ3ij denote the LWFs differences between ~τi and ~τj at its first, second and third NN sites j respectively. Ti A˜ij ’s are the coefficients of harmonic interactions in Hef τ ), which are deterf (~

mined from linear fits of the frequencies of optical phonons at high symmetry q-points (Sec. S5 in SI), and are listed in Table 6. Table 6: Coefficients of harmonic and anharmonic terms in effective Hamiltonian in unit of eV/f.u..

Coeff.

Values

Coeff.

Values

Coeff.

Values

A11

-2.6

A˜11

-4.5

B41

11.3

A12

3.3

A˜12

20.8

B42

-1.6

A21

4.6

A˜21

-1.2

B43

-16.9

A22

0.3

A˜22

1.7

B44

53.7

A23

-2.6

A˜31

-0.3

B45

1.7

A31

0.8

B31

0.8

B46

0.3

A32

3.2

B32

-6.9

B61

-286.1

A˜01

75.6

B33

-3.7

B81

1614.3

18

340

5.3. Hamiltonian of homogeneous strain: linear and nonlinear elastic energy Large deformations (homogeneous strain) are involved during MT in NiTi as a result of its soft second-order elastic moduli. We considered non-linear elasticity (Helastic (εαβ )) with terms upto fourth order allowed by the cubic symmetry [44, 45, 46, 47, 48, 49]. Using the Voigt notation (ε1 = εxx , ε2 = εyy , ε3 = εzz , ε4 = 2εyz , ε5 = 2εzx and ε6 = 2εxy ), elastic energy is expressed as:

Helastic (ε)

= + + +

N C1 (ε1 + ε2 + ε3 ) N {C11 (ε21 + ε22 + ε23 ) + 2C12 (ε1 ε2 + ε2 ε3 + ε3 ε1 ) + C44 (ε24 + ε25 + ε26 )} 2 N {C111 (ε31 + ε32 + ε33 ) + 3C112 (ε21 (ε2 + ε3 ) + ε22 (ε3 + ε1 ) + ε23 (ε1 + ε2 )) 6 6C123 ε1 ε2 ε3 + 3C144 (ε1 ε24 + ε2 ε25 + ε3 ε26 )

+ 3C155 (ε1 (ε25 + ε26 ) + ε2 (ε26 + ε24 ) + ε3 (ε24 + ε25 )) + 6C456 ε4 ε5 ε6 } N {C1111 (ε41 + ε42 + ε43 ) + 4C1112 (ε31 (ε2 + ε3 ) + ε32 (ε3 + ε1 ) + ε33 (ε1 + ε2 )) + 24 + 6C1122 (ε21 ε22 + ε22 ε23 + ε23 ε21 ) + 12C1123 ε1 ε2 ε3 (ε1 + ε2 + ε3 ) + 6C1144 (ε21 ε24 + ε22 ε25 + ε23 ε26 ) + 6C1155 (ε21 (ε25 + ε26 ) + ε22 (ε26 + ε24 ) + ε23 (ε24 + ε25 )) +

12C1255 (ε1 ε2 (ε24 + ε25 ) + ε2 ε3 (ε25 + ε26 ) + ε3 ε1 (ε26 + ε24 ))

+ 12C1266 (ε1 ε2 ε26 + ε2 ε3 ε24 + ε3 ε1 ε25 ) + 24C1456 ε4 ε5 ε6 (ε1 + ε2 + ε3 ) + C4444 (ε44 + ε45 + ε46 ) + 6C4455 (ε24 ε25 + ε25 ε26 + ε26 ε24 )}

where N is the number of unit cells, Cij , Cijk and Cijkl are second-, thirdand fourth-order elastic moduli, C1 represents a pressure term. Our estimates of these compliances are listed in Table 7. The detailed description of the 350

calculations of elastic moduli has been given in SI (Sec. S5). While third- and fourth-order elastic moduli of NiTi have not been determined theoretically or experimentally to the best of our knowledge, our estimates of the linear elastic moduli are in good agreement with earlier theoretical [9, 50] and experimental works [51] (Table S4 in SI). 5.4. Coupling of homogeneous strain with phonons Homogeneous strain and its coupling with LWFs (~η 0 s and ~τ 0 s) play a crucial role in MT. In Hspc (~ηi , ~τi , εαβ ), we include couplings between homogeneous

19

(8)

Table 7: Elastic moludi of Helastic (ε) in GPa, determined from first-principles.

Coeff.

Values

Coeff.

Values

Coeff.

Values

C1

-1.1

C144

-205.2

C1144

2081.3

C11

174.2

C155

-200.9

C1155

1700.1

C12

154.9

C456

-201.1

C1255

2173.7

C44

48.3

C1111

15514.7

C1266

457.6

C111

-463.8

C1112

-1806.4

C1456

616.8

C112

-544.4

C1122

4731.7

C4444

2144.5

C123

-385.8

C1123

759.5

C4455

527.9

strains and (i) quadratic terms in ~ηi ’s (differences between ~ηi at site i and ~ηj at its first NN sites j, see equation 9), and (ii) with products of Ti and Ni-centric 360

LWFs (see equation 10). The former captures the coupling between primary order parameter, M50 mode and strain, while the latter gives mixing between acoustic and optical phonons due to strain distortions. We find that the thirdorder strain-phonon coupling terms are not sufficient to capture the details of low energy phases, and hence include fourth-order coupling between axial strain and Ni-centric LWFs (see equation 11).

3 Hspc (~η , ε)

=

N X 6 X {g31 (ε1 + ε2 + ε3 )|~η1ij |2 i=1 j=1

2 2 2 + g32 (ε1 η1ijx + ε2 η1ijy + ε3 η1ijz )

+ g33 (ε1 + ε2 + ε3 )(~η1ij · dˆ1j )2 2 2 2 + g34 (ε1 η1ijx d21jx + ε2 η1ijy d21jy + ε3 η1ijz d21jz ) 2 2 2 2 2 2 + g35 (ε1 (η1ijy + η1ijz )d21jx + ε2 (η1ijz + η1ijx )d21jy + ε3 (η1ijx + η1ijy )d21jz )

+ g36 (ε4 η1ijy η1ijz + ε5 η1ijz η1ijx + ε6 η1ijx η1ijy ) + g37 (ε4 η1ijy η1ijz d21jx + ε5 η1ijz η1ijx d21jy + ε6 η1ijx η1ijy d21jz )}

20

(9)

3 Hspc (~η , ~τ , ε)

=

N 8 X X [h31 (ε1 + ε2 + ε3 ){ηix (τjy djy + τjz djz )djx + c.p} i=1

+

j=1

8 X h32 {ε1 ηix (τjy djy + τjz djz )djx + c.p} j=1

+

+

+

+

+

4 Hspc (~η , ε)

=

h33 {ηix

8 X

h34 {ε4 ηix h35 {ε4

(ε2 τjy djy + ε3 τjz djz )djx + c.p}

j=1

8 X

τjx djy djz + c.p}

j=1

8 X

(ηiy τjy + ηiz τjz )djy djz + c.p}

j=1

h36 {ε4 ηix h37 {ηix

8 X

(τjy djz + τjz djy )djx + c.p}

j=1

8 X

(ε5 τjy + ε6 τjz )djy djz + c.p}]

(10)

j=1

N X 6 X {g41 (ε21 + ε22 + ε23 )|~η1ij |2 i=1 j=1

2 2 2 + g42 (ε21 η1ijx + ε22 η1ijy + ε23 η1ijz )

2 2 2 + g43 (ε1 (ε2 + ε3 )η1ijx + ε2 (ε3 + ε1 )η1ijy + ε3 (ε1 + ε2 )η1ijz ) 2 2 2 + g44 (ε1 ε2 η1ijz + ε2 ε3 η1ijx + ε3 ε1 η1ijy )

(11)

where gij ’s denote couplings of strain with Ni-centric LWFs, and hij ’s are 370

the third-order coupling that cause strain-induced mixing between acoustic and optic phonons. These coefficients were determined from calculations of structures obtained by freezing M-point phonon eigenvectors at the different values of strain. Values of the third- and fourth-order coupling coefficients (See Table 8) are obtained as first and second derivatives of harmonic coefficients with respect to strain respectively (Sec. S5 in SI).

6. Monte Carlo Simulations We now analyze Hef f with MC simulations on a periodic system containing L×L×L (L = 16) unit cells of B2 structure. We used a single-flip update within 21

Table 8: Coefficients of strain phonon couplings part of effective Hamiltonian in unit of eV/f.u..

Coeff.

Values

Coeff.

Values

Coeff.

Values

g31

-2.8

g37

20.1

h32

33.8

g32

-7.9

g41

14.7

h33

48.3

g33

20.2

g42

130.5

h34

10.5

g34

-25.7

g43

-30.5

h35

-8.4

g35

-11.2

g44

-10.4

h36

25.8

g36

-18.0

h31

-22.6

h37

-16.4

Metropolis scheme and adjusted the step-size of configurational {~ηi , ~τi } updates 380

to maintain the acceptance ratio of ∼ 0.5 near the transformation temperature. In each Monte Carlo sweep (MCS), we picked ηi , τi randomly (totally 2L3

updates) and homogeneous strain variables L times. Thus, each MCS involves (2L3 + L) attempts of updating configurations. To identify different phases from the configurations sampled in MC simulations at a given temperature, we accumulated absolute values of Fourier components of ηi and τi , and obtained averages of (M50 , M20 and M40 modes), and averages of strain (<ε>) components at each temperature. To assess the possibility of hysteresis, we approach the transformation from high (low)-temperatures by cooling (heating) the system. In cooling simulations, we start from 200 K tem390

perature, equilibrate the system to the cubic phase, and reduce the temperature in steps of 5 K (and 1 K near the transformation temperature), down to a lowtemperature of 100 K. In heating simulations, we start from T=100 K taking B190 structure as the initial configuration and increase temperature in steps of 5 K (and 1 K near the transformation temperature) to a high-temperature of 200 K (we checked that our results do not depend on the choice of initial configuration). At each temperature, we used 10000 MCS for thermal equilibration and 50000 MCS for thermodynamic averaging. Average values of homogeneous strain components ε (Fig. 4a), M50 , M20 and M40 modes (Fig. 4b) vanish at high-temperatures, and change discon22

164 K Cooling 163 K 162 K

180

0.02

B19’

0 100

(d) B19’

0.03

0.06

M5’ (aB2 unit)

0.09

0.12

140 160 T (K)

B2 180

200

Cooling Heating

6

4 3 2

4

1 0

120

5

2

B2

50 0 -0.03

0.04

200

150 100

0.06

ω x 10 (M5’)

200

140 160 T (K)

0.08

Cooling

Cooling

<ε> Histogram (arbitrary unit)

250

120

Heating

-0.2 100

M5’ M2’ M4’

Heating

0 -0.1

(c)

B2

<|O|> (aB2 unit)

B19’

0.1

(b) 0.1

ε1 ε2 ε3 ε4 ε5 ε6

C4 (M5’)

(a) 0.2

0 100

3 2 1

100

120

150

T (K)

200

140 160 T (K)

180

200

Figure 4: (color online). Averages of components of strain tensor (a), absolute values of M50 , M20 and M40 modes (b) as a function of temperature at ambient pressure. Histograms of M50 mode obtained from configurations sampled during cooling (c) at TM and its two temperatures in its close vicinity. Square of frequency of M50 phonon and its Binder cumulant (d).

23

400

tinuously to non-zero values at low-temperatures. In cooling simulations, this discontinuity marks the martensitic transformation at TM = 163 K. In heating simulations, our model transforms to austenite phase at TA = 169 K. It is clear (Fig. 4a and 4b) that NiTi transforms from B2 to B190 structure directly ruling out the intermediate B19 structure. Results of our MC simulations reveal a sharp discontinuity in order parameters at TM . Even more remarkable is that they change only weakly below TM . Histograms of the primary order parameter obtained from cooling (Fig. 4c) simulations reveal (i) relatively narrow distributions and (ii) a sharp jump across TM . Bimodal nature of the histogram at TM confirms coexistence of B2 and B190 , and the first-order character of the

410

MT. In order to further corroborate these results, we present frequency (ω) and Binder cumulant (C4 ) of the primary order parameter (M50 mode) obtained in cooling simulations:

ω2

∝

C4

=

kB T (hO2 i − hOi2 )−1 N hO4 i , hO2 i2

(12) (13)

where N = L3 is the number of unit cells used in simulations and O is ~q = π/a(011) Fourier component of ηy . Softening in M50 phonon frequency and its T-coefficients (Fig. 4d) below and above TM , and a discontinuity in C4 (inset, Fig. 4d) of M50 mode at T = 163 K confirm that M50 mode is the primary order parameter of MT in NiTi. Our estimate of the equilibrium transition temperature (T0 = 12 (TM +TA )) [52] is underestimated with respect to 420

experiment (TM = 318 K, TA = 320 K and T0 = 319) [26]. This is similar to the efficacy of model Hamiltonian schemes in prediction of ferroelectrics structural transitions [43]. This is probably because of the GGA-DFT errors in energies and structural parameters, and needs further investigation. To assess the sensitivity of our estimates of T0 to DFT-errors in lattice constants, we simulated Hef f of NiTi at -5 GPa and 5 GPa pressures. At P =24

(a) 0.1

(b)

Cooling Heating

0.1

0.02 0 120

B2

B19’

160

200 T (K)

<|O|> (aB2 unit)

0.04

M5’ M2’ M4’

Cooling

0.06

Heating

<|O|> (aB2 unit)

0.08

0.05

0

240

M5’ M2’ M4’

B19’

B19

50

100 T (K)

B2

150

200

Figure 5: Absolute values of M50 , M20 and M40 modes as a function of temperature at P = -5GPa (a) and at P = 5GPa (b).

5 GPa, transformation from B2 to B190 structure (Fig. 5a) occurs at a notably higher temperature with hysteresis over a much wider range of temperature: TM ∼ 155 K, TA ∼ 225 K and T0 = 190 K. At P = 5 GPa, in heating (we start from T=10 K taking B190 structure as the initial configuration) simulations, 430

the MT occurs in two steps (i) B190 → B19 at 60 K, and (ii) B19 → B2 at 155 K (Fig. 5b), with a weaker hysteresis over ∆T = 5 K range. While in cooling simulations, B2 structure transforms into B19 structure at 155 K, and remains in this phase up on further cooling down to 10 K (B19 structure does not transform into B190 structure). This is probably because at P = 5 GPa, the energy difference between B19 and B190 structure is small (Table S6 in SI) and the energy barrier is higher, needing longer simulation for a transition to occur. Thus, B19 phase can arise as a metastable phase due to inhomogeneous pressure and stress fields in NiTi near the transformation temperature. Secondly, T0 as well as the range of temperature of hysteresis in MT increase with negative

440

pressure (i.e. increase in volume of the unit cell of B2 structure). We compared our pressure-induced martensitic transformation with the earlier experimental work [53], which has been done on Ni-rich (50.375 at.% Ni) composition of NiTi. Figure 6 of Ref. [53] shows that TM increases with increasing uniaxial stress, which is unlikely in our work probably because we have applied hydrostatic pressure on equiatomic NiTi.

25

We now discuss comparison between our results with the work in Ref. [12] in which martensitic transition temperature was estimated using free energies obtained using thermodynamic integration within AIMD. While TM predicted by our Hef f based analysis is underestimated with respect to experiment, free450

energy based analysis [12] gives an overestimated TM . As both the works are fundamentally based on density functional theory treated within the same approximation, it would appear that this discrepancy is likely due to further approximations or truncations in construction of the effective Hamiltonian. Indeed, this comparison uncovers important aspects of first-order structural phase transitions, which are typically fluctuation-driven first-order phase transition in which strain-phonon coupling plays a crucial role, as was shown in similar phase transitions in perovskite ferroelectrics [43, 54]. If fluctuations in strain (εαβ ) are not allowed in a simulation, called as clamped lattice analysis, the same effective Hamiltonian gives a second order phase transition at much higher Tc . To

460

verify this, we simulated our Hef f of NiTi with strain clamped at the lattice parameters of B2 and B190 ’s structures and found that (a) the thermal hysteresis vanishes, confirming a second order transition, and (b) Tc ’s are much higher (200 K for zero strains of the cubic lattice, and higher than 500 K for strains clamped at the B190 ’s lattice at P=0 GPa) (See Fig. S5). We note that free energy estimation in Ref. [12] is based on thermodynamic integration along a path in the structural (strain) space, and lattice constants are kept fixed in MD simulation at each point along this path. As fluctuations in strain are thus not included, the estimated TM ’s from the resulting free energies are higher. More work is necessary to obtain deeper understanding of the martensitic transitions

470

in this context. The hysteresis found in the vicinity of martensitic transition in our heating and cooling simulations is essentially due to first-order character of the martensitic transition, in a way similar to that found in simulations of ferroelectric phase transitions [43]. The origin of this has been identified in fluctuations generated by the strain-phonon coupling. While there is no time-scale in Monte Carlo simulations, the hysteresis does arise in our quasi-static simulations, in 26

which temperature is varied in very small steps, and the system is equilibrated at each temperature. Such simulated hysteresis should apply to experimental measurements carried out at very slow rates of heating and cooling. As the 480

systems simulated in our work (8192 atoms) are much too small, their phase transition occurs homogeneously, i.e. without nucleation and growth of domains (twins). Hence, it is not possible for such simulations to be directly related to Ms and Mf temperatures observed in experiments. The simulated phase transition occurs marking a sharp, global change in structure throughout the system. Indeed, our work forms a starting point of multi-scale modeling and simulations [55] that will allow us to treat much larger systems.

7. Landau Theoretical Analysis of Hef f To determine the couplings that are responsible for the observed MT in NiTi, we now present a Landau theory obtained by projecting the Hef f into 490

the subspace of phonons at M-point π/a(011) and strain degrees of freedom, Landau Ti Ni {ηx , ηy , τx } (equation 14), (ii) τi ) projected on Hph ηi ) and Hef (i) Hef f (~ f (~

Landau Hspc (~ηi , ~τi , εαβ ) projected on Hspc {ηx , ηy , τx , ε1 , ε2 , ε4 , ε5 } (equation 15), and Landau {ε1 , ε2 , ε4 , ε5 } (equation 12 in SI). To (iii) Helastic(εαβ ) projected on Hstr

make Landau theoretical analysis simple, we exclude the high frequency M50 optical phonon mode, which has a negligible contribution to the energetics of the relevant structures of NiTi (Table S7 and S8). While this simplification is made to facilitate analytical treatment, we include all the terms in our exact numerical analysis results of which are shown in Fig. 6. In the projected subspace, ηy and ηz span the subspace of doubly degenerate unstable M50 phonon 500

modes, ηx and τx span the subspace of M20 and M40 phonon modes respectively. We restrict to structures with the symmetry of orthorhombic and monoclinic phases of NiTi with ηz = ηy , ε3 = ε2 and ε6 = ε5 , and write Hamiltonian in these subspace:

Landau Hph (ηx , ηy , τx )

= a1 ηx2 + a2 ηy2 + a3 τx2 + a4 ηx4 + a5 ηy4 + a6 ηx2 ηy2 27

+

a7 (ηx2 + 2ηy2 )3 + a8 (ηx2 + 2ηy2 )4

Landau Hspc (ηx , ηy , τx , ε1 , ε2 , ε4 , ε5 )

=

(14)

(g1 ε1 + g2 ε2 )ηx2 + (g3 ε1 + g4 ε2 + g5 ε4 )ηy2

+ g6 ε5 ηx ηy + g7 ε5 τx ηy + g8 ε4 ηx τx +

(g9 ε21 + g10 ε22 + g11 ε1 ε2 )ηx2

+

(g12 ε21 + g13 ε22 + g14 ε1 ε2 )ηy2

(15)

where ai and gi are the coefficients of phonon and strain-phonon coupling terms in the Landau energy function respectively (their values are listed in Table 9). The instability of M50 mode (a2 =-8.4 eV/f.u. < 0) corresponds to {0¯ 11} <0 1 1> shuffle, whose coupling with ε1 (g3 =116.4 eV/f.u.) is positive, and with ε2 (g4 =-188.6 eV/f.u.) is negative. Therefore, ε1 must be negative and ε3 = ε2 510

positive to reduce the energy of system, consistent with Bain strain distortion that lowers the symmetry of cubic structure of NiTi to tetragonal one. Coupling of ηy with ε4 ({010}<0 0 1> pure shear, g5 =-143.7 eV/f.u.), transforms B2 structure to B19 structure. Though ηx and τx represent stable phonons at Mpoint and linear elastic moduli are positive definite [56], the non-basal shear strain ε6 = ε5 couples strongly with ηx , τx and ηy (g6 = -126.3 eV/f.u. and g7 =137.4 eV/f.u.) stabilizing B190 structure relative to B19 structure through structural distortions involving M20 and M40 modes. Rest of the terms in Landau energy function play a supporting role in giving correct details of energetics and lattice parameters of B19 and B190 structures. From exact analysis including

520

all the terms, we determined the phase diagram of stability of B19 and B190 structures (see Fig. 6) in the plane of these g-couplings, which can thus be used as descriptors of stability of these phases. If the values of these couplings (g6 and g7 ) were weak, B2 structure of NiTi would transform through its MT to B19 structure, not B190 . The primary order parameter (phonon mode M50 ) drives the MT, M20 and M40 modes and strains constitute the secondary order parameters of the MT in NiTi that determine stabilization of B19 vs B190 phases. 28

1.5

g/g7

1.2

B19’

0.9 0.6

B19

0.3 0 0

0.3

0.6 0.9 g’/g6

1.2

1.5

Figure 6: Phase diagram showing the relative stability criteria between B19 and B190 structures of NiTi as a function of third-order strain-phonon coupling coefficients g and g 0 .

Our analysis suggests that the three order parameters [29] are not sufficient to study MT in NiTi.

8. Structural disorder and microstructures relevant to MT in NiTi 530

In addition to bulk crystalline phases, a martensitic transformation invariably involves planar faults and microstructures (e.g. twinning and stacking faults). We now demonstrate the efficacy of Hef f in capturing the physics of phases with structural disorder and microstructures in NiTi. We present stability analysis of (i) average austenitic structures reported in Refs. [57, 58], and (ii) (100)[011]B2 twinning and stacking faults that occur in the martensitic B190 structure. 8.1. Stablity of average austenitic structures To simulate the reported average austenitic structures, we start with B2 structure and relax its structural distortions in the subspace of unstable M50

540

phonons using Hef f in MC simulation at T=10 K, constraning ~τi0 s and ε to be zero. We find a structure resulting from a combination of two unstable M50 phonons at distinct M-points, π/a(011) and π/a(101), which break the 29

Table 9: Coefficients of terms in Landau energy function are linear combinations of the parameters in Hef f . Landau coeff.

Hef f coeff.

Values (eV/f.u)

a1

4A11 + 4A21 + 4A22

9.2

a2

8A11 + 4A12 + 4A21 + 4A22 + 8A23

-8.4

a3

1 ˜ 2 A01

+ 4A˜11 + 8A˜21 + 4A˜22

17.2

a4

32B41 + 64B45

468.9

a5

64B41 + 32B42 + 32B43 + 64B44 + 256B45

4.0×103

a6

64B42 + 64B44 + 256B45

2.1×103

a7

128B61

-3.7×104

a8

512B81

8.3×105

g1

8g31 + 8g32

-85.8

g2

16g31 + 8g35

-134.2

g3

16g31 + 8g33

116.4

g4

32g31 + 16g32 + 16g33 + 8g34 + 8g35

-188.6

g5

8g36

-143.7

g6

16g36 + 8g37

-126.3

g7

16 3 h36

137.4

g8

8 3 h34

28.0

g9

8g41 + 8g42

1.2×103

g10

16g41 + 8g43

-9.1

g11

16g44

-167.0

g12

16g41

235.0

g13

32g41 + 16g42 + 16g44

2.4×103

g14

16g43 + 16g44

-655.3

translational symmetry of B2 structure, transforming it to a stable hexagonal austenite structure reported in Refs. [57, 58]. This hexagonal structure (see Fig. 7a) has a periodic unit cell with 24 atoms (12 Ni and 12 Ti atoms) and √ √ ahex = bhex = 2 2a and chex = 3a (a is the lattice parameter of B2 structure) along [0¯ 11]B2 , [¯ 101]B2 and [111]B2 directions respectively. It is lower in the energy by 9.8 meV /f.u. than the B2 structure, and higher in the energy by 71 meV /f.u. than the B190 structure. To stabilize this austenitic structure, Ni 550

and Ti atoms displace from their high symmetry positions in the B2 structure

30

along the eigenvectors of unstable M50 phonons such that the average structure remains B2. This is one of the many possible stable austenitic structures, and similar to the hexagonal structure of Refs. [57, 58]. The main difference between this hexagonal structure and that of Refs. [57, 58] is the size of the respective periodic unit cell. The hexagonal structure reported in Refs. [57, 58] has 54√ √ atoms unit cell, where ahex = bhex = 3 2a and chex = 3a. Thus, Hef f does capture the physics and stability of average austenite structure like the one reported in Refs. [57, 58]. 8.2. Microstructure of NiTi 560

To examine microstructures relevant to the MT in NiTi, we simulated (100)[011]B2 twinning and stacking faults in its martensitic B190 structure using Hef f . Separated by a (100) plane of the B2 structure, we introduced two domains (See Fig. 7b-7d) in a L × L × L system. Domains (I, x ≤ L/2) and (II, x > L/2) contain states of B190 structure with distinct orientation, charecterized by or-

der parameters. In the twinned structure, M50 and non-basal shear (ε5 and ε6 ) order parameters have opposite values in the two domains. On the other hand, in stacking faulted structure, M50 , M20 and M40 have opposite values in the two domains though non-basal shear strains are the same. Hydrostatic, Bain and pure shear strains (secondary order parameters) are same in both domains 570

of twinned and stacking faulted structures. We constructed the twinning and stacking faults in B190 structure of NiTi using the optimized values of order parameters obtained by Hef f ( see Table S8 in SI). Our estimates of the twinning and stacking fault energies (using equation (13) in SI) are 96.3 mJ/m2 and 217.4 mJ/m2 , respectively, confirming that twinning faults are more likely to occur than stacking faults in NiTi. Our MC simulations at T= 10 K initialized with the structures containing twinning and stacking faults reveal that both the structures remain stable in their respective forms, and that the twinned structure is lower in energy than the stacking faulted structure (see Table S9 in SI). While our Hef f -based estimates of twinning and stacking faults of NiTi are

580

bit different in comparison with their DFT values [59], their relative stability is 31

Hexagonal austenitic structure (a)

(1 1 1)B2 plane view

Single Crystal B19' Ni

(b)

Ti [-1 0 1]B2

[0 -1 1]B2

M5'(+)

ε5=ε6= ε

M5'(+)

ε5=ε6= ε

Stacking fault B19'

Twinned B19' (c)

(d)

Domain (II)

M5'(-)

ε5=ε6=-ε

Domain (I)

M5'(+)

ε5=ε6= ε

M5'(-)

ε5=ε6= ε

M5'(+)

ε5=ε6= ε

[1 0 0]B2

[0 1 1]B2

Figure 7: (Color online) Analysis based on Hef f gives stable atomic structure of the hexagonal austenitic phase (a). Optimized atomic structures of single-crystal (b), twinned B190 (c), and stacking faulted (d) B190 phase of NiTi obtained from MC simulations of Hef f at 10 K temperature. Blue boxes in (b), (c) and (d) represent B190 structural unit cell. M50 (+) is the primary order parameter, and M50 (−) represents its opposite value. ε5 and ε6 are the components of non-basal shear strain, which are the same in both domains of the stacking faulted structures, while have opposite values in the two domains of twinned structure of NiTi.

captured qualitatively well by our effective Hamiltonian. In the phonon-based theoretical framework of our Hef f , a ferroelastic domain wall separating domains of B190 phases with distinct order parameter orientation gives the planar twinning faults are known to be relevant to the MT in NiTi.

9. Summary In summary, we have presented a microscopic picture of martensitic phase transformation in NiTi in terms of its phonons and their interactions derived from first-principles and captured in an effective Hamiltonian. We have shown that vibrational entropy of soft modes stabilizes the monoclinic martensite B190

32

590

structure over its BCO structure at T>43 K. Through Monte Carlo simulations of an effective Hamiltonian derived to capture its low-energy landscape of NiTi, we (a) determine its soft modes and establish the cell-doubling M50 phonon of the cubic phase as the primary order parameter of MT, and (b) there are SIX other secondary order parameters that are relevant to the MT in NiTi. We show that pressure can introduce an intermediate phase during the MT, thanks to the interesting physics of seven coupled order parameters. Using Landau theoretical analysis, we show that relative strengths of the third-order coupling between primary and secondary order parameters including strain determine the specific symmetry of low-T structures emerging from its MT. These can be

600

used as first-principles descriptors in designing materials with improved shape memory properties through substitutional alloying in NiTi. We finally show that our Hef f qualitatively captures stabilities of other stable austenitic structures, twinning and stacking faults in NiTi, which are relevant to its MT.

10. Acknowledgments We acknowledge useful discussions with K. M. Rabe, D. Banerjee and U. Ramamurty. U. V. Waghmare acknowledges funding from a JC Bose National Fellowship of DST, Govt of India and Sheikh Saqr Fellowship. Pawan Kumar acknowledges funding from the India-Korea Virtual Network Center in Computational Materials Science supported by DST, Govt of India.

610

References [1] K. Otsuka, T. Kaseshita, Science and technology of shape-memory alloys: New developments, Mater. Res. Bull. 27 (2002) 91–100. [2] G. M. Michal, R. Sinclair, The structure of TiNi martensite, Acta Cryst. B 37 (1981) 1803–1807. [3] P. A. Lindgard, O. G. Mouritsen, Theory and model for martensitic transformations, Phys. Rev. Lett. 57 (1986) 2458–2461. 33

[4] E. Goo, R. Sinclair, The B2 to R transformation in Ti50 Ni47 Fe3 and Ti49.5 Ni50.5 alloys, Acta Metallurgica 33 (1985) 1717–1723. [5] Y. Wang, A. G. Khachaturyan, Three-dimensional field model and com620

puter modeling of martensitic transformations, Acta Materialia 45 (1997) 759–773. [6] Y. Kudoh, M. Tokonami, S. Miyazaki, K. Otsuka, Crystal structure of the martensite in Ti-49.2 at.% Ni alloy analyzed by the single crystal X-ray diffraction method, Acta Metallurgica 33 (1985) 2049–2056. [7] M. J. Marcinkowski, A. S. Sastri, D. Koskimaki, Martensitic behaviour in the equi-atomic NiTi alloy, Philosophical Magazine 18 (1968) 945–958. [8] F. E. Wang, S. Pickart, H. A. Alperin, Mechanism of the NiTi martensitic transformation and the crystal structures of NiTi-II and NiTi-III phases, Journal of Applied Physics 42 (1972) 97–112.

630

[9] Y. Y. Ye, C. T. Chan, K. M. Ho, Structural and electronic properties of the martensitic alloys NiTi, TiPd, and TiPt, Phys. Rev. B 56 (1997) 3678–3689. [10] X. Huang, C. Bungaro, V. Godlevsky, K. M. Rabe, Lattice instabilities of cubic NiTi from first principles, Phys. Rev. B 65 (2001) 014108–015112. [11] S. Kadkhodaei, A. van de Walle, First-principles calculations of thermal properties of the mechanically unstable phases of the PiTi and NiTi shape memory alloys, Acta Materialia 147 (2018) 296 – 303. [12] J. B. Haskins, A. E. Thompson, J. W. Lawson, Ab initio simulations of phase stability and martensitic transitions in NiTi, Phys. Rev. B 94 (2016) 214110–214119.

640

[13] Z.-Y. Zeng, C.-E. Hu, L.-C. Cai, X.-R. Chen, F.-Q. Jing, Lattice dynamics and phase transition of NiTi alloy, Solid State Communications 149 (2009) 2164–2168.

34

[14] K. Parlinski, M. Parlinska-Wojtan, Lattice dynamics of NiTi austenite, martensite, and R phase, Phys. Rev. B 66 (2002) 064307–064314. [15] X. Huang, G. J. Ackland, K. M. Rabe, Crystal structures and shapememory behaviour of NiTi, Nature Materials 2 (2003) 307–311. [16] K. Otsuka, Physical metallurgy of NiTi-based shape memory alloys, Progress in Materials Science 50 (2005) 511–678. [17] T. Fukuda, 650

T. Saburi,

T. Chihara,

Y. Tsuzuki,

Mechanism of

B2 − B19 − B19 transformation in shape memory NiTiCu alloys, Materials Transactions JIM 36 (1995) 1244–1248. [18] K. G. Vishnu, A. Strachan, Phase stability and transformations in NiTi from density functional theory calculations, Acta Materialia 58 (2010) 745– 752. [19] M. Sanati, R. C. Albers, F. J. Pinski, Electronic and crystal structure of NiTi martensite, Phys. Rev. B 58 (1998) 13590–13593. [20] M.-X. Wagner, W. Windl, Lattice stability, elastic constants and macroscopic moduli of NiTi martensites from first principles, Acta Materialia 56 (2008) 6232–6245.

660

[21] S. Kibey, H. Sehitoglu, D. D. Johnson, Energy landscape for martensitic phase transformation in shape memory NiTi, Acta Materialia 57 (2009) 1624–1629. [22] J. M. Zhang, G. Y. Guo, Microscopic theory of the shape memory effect in NiTi, Phys. Rev. Lett. 78 (1997) 4789–4792. [23] R. Mirzaeifar, K. Gall, T. Zhu, A. Yavari, R. DesRoches, Structural transformations in NiTi shape memory alloy nanowires, Journal of Applied Physics 115 (2014) 194307–194314. [24] Y. Zhong, K. Gall, T. Zhu, Atomistic study of nanotwins in NiTi shape memory alloys, Journal of Applied Physics 110 (2011) 033532–033540. 35

670

[25] D. Mutter, P. Nielaba, Simulation of structural phase transitions in NiTi, Phys Rev. B 82 (2010) 224201–224209. [26] R. R. Pezarini, H. S. Bernab´e, F. Sato, L. C. Malacarne, N. G. C. Astrath, J. H. Rohling, A. N. Medina, R. D. dos Reis, F. C. G. Gandra, On the use of photothermal techniques to study NiTi phase transitions, Mater. Res. Express 1 (2014) 026502–026510. [27] R. F. Hehemann, G. D. Sandrock, Relations between the premartensitic instability and the martensite structure in TiNi, Scripta Metallurgica 5 (1971) 801–805. [28] K. Otsuka, X. Ren, Martensitic transformations in nonferrous shape mem-

680

ory alloys, Mater. Sci. Eng. A 273 (1999) 89–105. [29] X. Ren, K. Otsuka, The role of softening in elastic constant C44 in martensitic transformation, Scripta Materialia 38 (11) (1998) 1669 – 1675. [30] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. D. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, R. M. Wentzcovitch, Quantum espresso: a modular and open-source software project for

690

quantum simulations of materials, J. Phys. Condens. Matter 21 (2009) 395502–395520. [31] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (1996) 3865–3868. [32] J. P. Perdew, A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23 (1981) 5048– 5079.

36

[33] J. Sun, A. Ruzsinszky, J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115 (2015) 036402–36407. 700

[34] J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. V. Waghmare, X. Wu, M. L. Klein, J. P. Perdew, Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional, Nature Chemistry 8 (2015) 831–836. [35] G. Kresse, J. Furthm¨ uller, Efficient iterative schemes for ab initio totalenergy calculations using a plane-wave basis set, Phys. Rev. B 54 (1996) 11169–11186. [36] G. Kresse, J. Furthm¨ uller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1996) 15–50.

710

[37] G. Kresse, D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59 (1999) 1758–1775. [38] W.-S. Ko, B. Grabowski, J. Neugebauer, Development and application of a ni-ti interatomic potential with high predictive accuracy of the martensitic phase transition, Phys. Rev. B 92 (2015) 134107–134128. [39] D. Holec, M. Fri´ ak, A. Dlouh´ y, J. Neugebauer, Ab initio study of pressure stabilized niti allotropes: Pressure-induced transformations and hysteresis loops, Phys. Rev. B 84 (2011) 224119–224126. [40] J. Lee, Y. Ikeda, I. Tanaka, Solution effect on improved structural compatibility of niti-based alloys by systematic first-principles calculations, Journal

720

of Applied Physics 125 (2019) 055106–055123. [41] J. R. Chelikowsky, N. Troullier, K. Wu, Y. Saad, Higher-order finitedifference pseudopotential method: An application to diatomic molecules, Phys. Rev. B 50 (1994) 11355–11364.

37

[42] K. M. Rabe, U. V. Waghmare, Localized basis for effective lattice hamiltonians: Lattice wannier functions, Phys Rev. B 52 (1995) 13236–13246. [43] U. V. Waghmare, K. M. Rabe, Ab initio statistical mechanics of the ferroelectric phase transition in PbTiO3 , Phys Rev. B 55 (1997) 6161–6173. [44] H. Wang, M. Li, Ab initio calculations of second-, third-, and fourth-order elastic constants for single crystals, Phys Rev. B 79 (209) 224102–224111. 730

[45] X. Li, S. ch¨ onecker, J. Zhao, L. Vitos, B. Johansson, Elastic anharmonicity of bcc Fe and Fe − based random alloys from first-principles calculations, Phys. Rev. B 95 (2017) 024203–024213. [46] M. Lopuszy´ nski, J. A. Majewski, Ab initio calculations of third-order elastic constants and related properties for selected semiconductors, Phys. Rev. B 76 (2007) 045202–045210. [47] D. Holec, M. Fri´ ak, J. Neugebauer, P. H. Mayrhofer, Trends in the elastic response of binary early transition metal nitrides, Phys. Rev. B 85 (2012) 064101–064109. [48] M. D. Jong, I. Winter, D. C. Chrzan, M. Asta, Ideal strength and ductility

740

in metals from second- and third-order elastic constants, Phys. Rev. B 96 (2017) 014105–014116. [49] J. Zhao, J. M. Winey, Y. M. Gupta, First-principles calculations of secondand third-order elastic constants for single crystals of arbitrary symmetry, Phys. Rev. B 75 (2007) 094105–09411. [50] N. Hatcher, O. Y. Kontsevoi, A. J. Freeman, Role of elastic and shear stabilities in the martensitic transformation path of NiTi, Phys. Rev. B 80 (2009) 144203–144220. [51] O. Mercier, K. N. Melton, G. Gremaud, J. Hagi, Singlecrystal elastic constants of the equiatomic NiTi alloy near the martensitic transformation,

750

Journal of Applied Physics 51 (1980) 1833–1834. 38

[52] R. Salzbrenner, M. Cohen, On the thermodynamics of thermoelastic martensitic transformations, Acta Metallurgica 27 (1979) 739 – 748. [53] G. J. Pataky, E. Ertekin, H. Sehitoglu, Elastocaloric cooling potential of NiTi,Ni2 FeGa, and CoNiAl, Acta Materialia 96 (2015) 420 – 427. [54] K. M. Rabe, U. V. Waghmare, Strain coupling in perovskite structural transitions: A first principles approach, Ferroelectrics 194 (1997) 119–134. [55] U. V. Waghmare, First-principles theory, coarse-grained models, and simulations of ferroelectrics, Accounts of Chemical Research 47 (2014) 3242– 3249. 760

[56] F. Mouhat, F. X. Coudert, Necessary and sufficient elastic stability conditions in various crystal systems, Phys Rev. B 90 (2014) 224104–224107. [57] N. A. Zarkevich, D. D. Johnson, Stable atomic structure of NiTi austenite, Phys. Rev. B 90 (2014) 060102–060105 (R). [58] N. A. Zarkevich, D. D. Johnson, Shape-memory transformations of NiTi: Minimum-energy pathways between austenite, martensites, and kinetically limited intermediate states, Phys. Rev. Lett. 113 (2014) 265701–265705. [59] T. Ezaz, H. Sehitoglu, H. Maier, Energetics of twinning in martensitic NiTi, Acta Materialia 59 (2011) 5893 – 5904.

39

Declaration of interests The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: U. V. Waghmare acknowledges funding from a JC Bose National Fellowship of DST, Govt of India and Sheikh Saqr Fellowship. Pawan Kumar acknowledges funding from the India-Korea Virtual Network Center in Computational Materials Science supported by DST, Govt of India.