The Bain versus Nishiyama–Wassermann path in the martensitic transformation of Fe

Using atomistic simulations and an embedded-atom interatomic potential, which is capable of describing the martensite (fcc → bcc) transition in Fe, we compare the Bain and Nishiyama–Wassermann transformation paths. We calculate the minimum-energy paths for these two transformations at 0 K. For finite temperature, we study the evolution of the free energy along the two paths, calculated via the method of metric scaling, which shows only small differences. However, the variation of the elements of the stress tensor, and of the atomic volume, show clear differences: the Bain path leads to by a factor of five higher compressive pressures compared to the Nishiyama–Wassermann path. This means that the Bain path requires additional work applied on the system in order to accomplish it, and gives atomistic evidence that the martensite transformation rather follows the Nishiyama–Wassermann path in reality.


3
an Fe solid is decisive [11,12]. The philosophy behind our approach is that in such a classical study the empirical interatomic interaction potential incorporates in an effective sense these interactions.

Transformation paths and interatomic potential
As is well known, at 0 GPa, iron has three allotropes [13,14], of which two are of interest in the present study. α-iron (ferrite) with a bcc crystal structure is stable at room temperature and ambient pressure up to 1185 K. Between the Curie temperature of 1043 and 1185 K, the non-ferromagnetic phase is sometimes termed the β-phase. When the temperature is increased beyond 1185 K, it transforms to the fcc structure, called γ -iron or austenite. The reverse change from γ -to α-iron can occur in two different ways depending on the cooling speed [15]. At a sufficiently slow rate, the atomic motion is uncoordinated and the transformation is diffusive. But if the fcc phase is rapidly cooled, there is not enough time to allow for diffusion and the transformation occurs by a cooperative shear mechanism; this latter transformation is called the martensitic transformation of Fe.
From a structural point of view, the martensitic transformation is characterized by the orientational relationship between the two phases, which can be characterized by an associated transformation path. The most intuitive transformation path from the fcc to the bcc crystal structure was proposed by Bain in 1924 [1] as noted above. After defining a bcc unit cell in the fcc crystal, see figure 1(a), we contract the z-axis by about 21% and expand the xand y-axis by about 12%. In spite of its simplicity, the orientational relationship between the fcc (parent) and the bcc (martensite) phases, is not observed in experiment. Closer to what really happens is the Nishiyama-Wassermann path [4,16], which satisfies the orientational relationship We may visualize the entire transformation in two steps, see figure 1(b). In the first step, we fix a {111} fcc plane and the two neighbouring {111} fcc planes are shifted in opposite directions with a common uniform shear of 19 • 28 in the [112] fcc direction. The second step consists of an orthogonal expansion/contraction in order to obtain the desired atomic density.
It is not easy to capture the allotropy of Fe with the help of an empirical interatomic potential in a reliable model to use in molecular dynamics simulations. Usually the models, mainly in the embedded-atom method (EAM) approach, are fitted to low temperature properties at ambient pressure, keeping in mind a defined phase for a particular problem. Only a few interatomic potentials have shown the capability to reproduce the solid-solid phase transition of Fe [17], such as the Meyer-Entel potential [18] or the potential developed by Müller et al which is of the bond-order type [19]. We have used here the Meyer-Entel potential because it is relatively simple, it is fast enough to simulate systems with millions of atoms for long time intervals of the order of nanoseconds using modest computational resources, and it has already produced good results for pure iron and nickel ferrous alloys [10,18], [20]- [22].

Zero temperature: minimum-energy path (MEP)
In figure 2, we investigate in detail the potential energy landscape of the Meyer-Entel potential at 0 K and explore the details of the Bain and the Nishiyama-Wassermann path in it. The Bain path evolves along a tetragonal distortion of the lattice, so it can be uniquely characterized by two linear scaling parameters. We choose s x y to describe the linear deformation in the xand the y-direction (cf figure 1(a)) and s z in the z-direction; in each case s = 0 characterizes the fcc, and s = 1 the bcc structure. Figure 2(a) demonstrates that bcc is indeed the ground state structure in this potential, and fcc forms a (metastable) local minimum. We indicate the MEP in this figure; it is defined as the unique path which connects two local minima and runs orthogonal to the equipotential contours of the potential energy surface. It runs through the saddle point and is commonly used in physical chemistry and elsewhere to describe the transformation of reactants to products. We calculate this path by using a variant of the sphere optimization technique [23]. Our analysis demonstrates that the MEP runs close to the diagonal, i.e. s x y ∼ = s z ; note that as a consequence volume does not change linearly along this path (cf also figure 4 below). The saddle point is located at s x y = 0.22, s z = 0.21. Figure 2(b) demonstrates that along this path, also the pressure is quite small; while a 'minimum-pressure path' would look slightly different from the MEP, it would also essentially follow the diagonal s x y = s z . A parameterization of the Nishiyama-Wassermann path is slightly more intricate. While a full characterization would require four parameters, we choose here the non-dimensional volume change s volume (again normalized such that s volume = 0 corresponds to fcc, and s volume = 1 to bcc) instead of three independent changes in linear dimensions of the unit cell. A second parameter is given by the shear s shear parallel to the (111) plane with s shear = 0 for fcc and s shear = 1 (19 • 28 ) for bcc. Figures 2(c) and (d) demonstrate that in this representation, the Nishiyama-Wassermann path runs close to the diagonal, s volume = s shear , and that this path is close to a 'minimum-pressure path'.

Finite temperature
We use this information when performing finite-temperature simulations. In these, it is not the 'minimum energy' but a 'minimum free-energy' path which is decisive for the transformation. In view of the considerable increase in computational cost for free-energy calculations (see below), we calculate the free energies not in the two-dimensional parameter space as indicated in figure 2, but only along the diagonal; our 0 K calculations may be used to justify this procedure.
We have calculated the free-energy difference along the Bain and Nishiyama-Wassermann transformation paths using the metric-scaling technique [17]. We start with the Helmholtz free-energy difference dA for an isothermal process given by dA −δW. ( By including the general expression for mechanical work δW as a result of stress-induced infinitesimal strains d i j we have where σ i j is the (i, j) element of the stress tensor. Now we consider the fcc → bcc transformation. We can perform the change as reversibly as possible by subdividing the transformation path into small steps. The final free-energy difference will be given by A = k dA k . By considering the reverse process, we can find bounds for A as We perform our finite temperature molecular-dynamics simulations using the LAMMPS code [24]. Our simulation box contains 12288 atoms in an fcc crystal configuration. Our cubic simulation box employs periodic boundary conditions, whereby the size and shape of the cell are allowed to change during the simulations. Before starting the simulation, we let the system relax for 40 ps at the desired temperature using the Nosé-Hoover thermostat.
In the following, we describe the simulation technique for the Nishiyama-Wassermann transformation; for the Bain path, we proceed in an analogous way. We change the shear angle by steps of 1% of the total shear angle 19 • 28 . After each step, the system is equilibrated with the Nosé-Hoover thermostat during 20 ps. In all simulations, simultaneously to shearing, we linearly scale the edges of the simulation volume in order to obtain a smooth transition from the fcc to the bcc density. Thus, after 2000 ps we obtain the bcc structure at its equilibrium density. We have also calculated the free-energy differences in the reverse direction, from bcc to fcc, in order to obtain the second bound to the total free-energy difference, cf equation (5). Figure 3 displays the free-energy differences at four temperatures for the two transformation paths. The error, as obtained from the bounds to the free energy, equation (5), is about 1%. This plot demonstrates that the transition temperature is around 500 K, in close agreement with the results obtained via thermodynamic integration [17]. Since the transformation must be independent of the path, the difference between the free energies of the bcc phase as evaluated by following the two paths is an indication of the systematic accuracy of our calculations. Note that within this accuracy, also the transformation barriers for the two paths coincide. This is astonishing and indicates that for processes occurring in thermodynamic equilibrium, the two paths are equivalent.
However, in reality the transformation does not occur in equilibrium, but in a strained and stressed environment as well as in the presence of other martensitic variants. Therefore, we also calculated the variation of the elements of the stress tensor, σ i j , along the transformation paths. From these we calculate the total pressure given by The corresponding data are shown in figure 4. They reveal a distinct difference between the two paths. Both the shear components and the total pressure are considerably larger along the Bain path than along the Nishiyama-Wassermann path. For the Bain path, the compressive pressure is higher by a factor of 5 than for the Nishiyama-Wassermann path. For the Nishiyama-Wassermann path, the diagonal elements of the stress tensor partially compensate, leading to a small hydrostatic pressure, while for the Bain path, the diagonal elements add up to a sizeable pressure. This behaviour can be rationalized by considering the evolution of the atomic volume in figure 4(d), which shows a considerably smoother variation for the Nishiyama-Wassermann path.
The transition temperature T 0 is by definition that temperature at which the free energies of the austenitic and martensitic phases cross as a function of T (whereby the free energy of each phase is artificially extrapolated to the corresponding other phase)-or it is approximately the T 0 of Landau theory of first-order transition. In experiment, the martensitic transformation in alloys usually involves a martensite start (M s ) and final (M f ) temperature (and an austenite start (A s ) and final (A f ) temperature for the reverse transformation); however, the corresponding transformations in elemental iron occur over a much narrower temperature range than in the alloys. In the simulations of alloys one may obtain a similar hysteresis behaviour if a sufficient number of defects is retained in the simulations-although the magnitude of the hysteresis is often unrealistically large because one does not simulate the real composition of the alloy. In the case of elemental iron the transformation in the supercells used in the simulations occurs more or less spontaneously, so that the relationship between the transition temperature and the experimentally determined temperatures may be identified tentatively as T 0 = (M s + A s )/2 [25].  is σ x z = σ yz = 0. Subfigure (c) compares the evolution of the total pressure along the two paths, and subfigure (d) the evolution of the atomic volume.

Conclusions
We conclude that our simulations offer an atomistic explanation of the question why in elemental iron, the Nishiyama transformation path is preferred over the Bain path. While the free energy itself-and hence the work done along the path-is very similar along both paths, the shear stresses and in particular the compressive pressure is considerably higher along the Bain path. This means that the Bain path requires additional work applied on the system in order to accomplish it. This finding gives atomistic evidence that the martensitic transformation in iron rather follows the Nishiyama-Wassermann path in any real experiment. We note that the use of boundary conditions and our restricted system size may affect the fluctuations and hence the evaluation of statistical averages at finite temperature. However, since we deal with effectively one single martensitic variant in a simulation box with periodic boundary conditions and use long equilibration times before starting the actual simulation, we do not expect qualitative changes for larger system sizes concerning the explanation offered in this paper of why the Nishiyama path is favoured over the Bain path.
In summary, our analysis indicates the optimum, i.e. minimum-energy, path which the system takes during transformation. The optimum path for the Bain transformation deviates from the usual (constant-volume) path, while for the Nishiyama-Wassermann transformation, we find a simultaneous volume expansion combined with shear. Using these paths, we have calculated the free energy for each of the transformation paths corresponding to the work done along each path. While the work, i, j σ i j d i j , is similar along the two paths, the stresses, σ i j , are not. This explains why elemental iron will follow the Nishiyama path when undergoing a martensitic transformation. This conclusion holds as long as we consider the formation of a single martensitic variant. In the real highly stressed system with multiple martensitic variants competing with each other, the situation may perhaps change. Such a situation is beyond the scope of this work.
With respect to the Kurdjumov-Sachs transformation [5] mentioned in the introduction, we note that this is similar but not identical to the Nishiyama-Wassermann relationship. Instead of equation (2) see figure 1(c). The similarity follows from the observation that both transformation mechanisms consist of a main shear of 19 • 28 in [112] fcc direction in the (111) fcc plane and that in both cases, the close-packed planes of austenite and martensite are parallel to each other. However, the close-packed directions of austenite and martensite are parallel in the case of Kurdjumov-Sachs while they deviate slightly (by 5 • 16 ) from each other in the Nishiyama-Wassermann case (see, e.g. the corresponding figures in [26,27]). Because of the similarity and only slight deviation between Kurdjumov-Sachs and Nishiyama-Wassermann, we therefore expect that the energetics calculated here for the Nishiyama-Wassermann path will not change too drastically when considering the Kurdjumov-Sachs path. An additional remark concerns the concept of an order parameter when dealing with the formation of martensite. Here, again, for the Bain transformation, the tetragonal distortion, i.e. the ratio c/a of the fcc-bcc transformation has often been used in statistical models. A corresponding order parameter may be used for describing the Nishiyama transformation. However, at finite temperatures, in particular in the non-equilibrium situation of a highly stressed system, not in thermal equilibrium during the transformation, the use of a conserved order-parameter concept may not always be useful. Molecular dynamics or Monte Carlo simulations may yield reasonable results also for non-equilibrium systems. Nonetheless, the characterization of the transformation by a parameter c/a or by a two-component parameter comprising shear and shuffle as for the bcc-hcp transformation used in [28], is still a very useful concept.