Thermal and structural properties of the martensitic transformations in Fe7Pd3 shape memory alloys: an ab initio-based molecular dynamics study

Ferromagnetic shape memory alloys, including the Fe7Pd3 system, constitute an upcoming class of functional materials, whose atomic-scale physical foundations are still insufficiently understood. The present work employs molecular dynamics simulations, based on ab initio derived embedded atom method potentials, to study martensitic transformations and twin variant reorientation. We address thermal and stress induced austenite-martensite transitions, twinning, as well as twin boundary mobility. While the predicted thermal properties are in accordance with experimental observations, we explore the detailed crystallography underlying transformation as well as twin boundary motion.


Introduction
The ferromagnetic shape memory effect (FSME) opens up possibilities for efficiently working magnetomechanic transducer devices, as they promise to achieve strains of several percent, which exceeds traditional magnetostriction by several orders of magnitude [1]. Two different manifestations of this effect need to be distinguished. One possibility of realization is based on magnetically induced reorientation of variants. Alternatively, the critical temperature for a thermally induced solid-solid phase transition can be shifted relative to the ambient temperature, yielding magnetically induced martensite-austenite transformation [2,3]. While originally demonstrated for the Ni 2 MnGa alloy by Ullakko et al [1], the Fe 7 Pd 3 ferromagnetic shape memory alloy promises advantages, particularly in the field of biomedical applications, due to demonstrated biocompatibility and -functionalizeability [4][5][6][7][8].
However, towards the usage of the FSME in Fe 7 Pd 3 for advanced applications there are still ambiguities concerning the exact details of the incorporated diffusionless (first order) phase transitions. The phase diagram was originally explored long before the discovery of the FSME while studying Invar alloys [9,10]. It is instructive for the present study to point out a few peculiarities on the martensitic and austentic phases of the Fe-Pd system. Crystallographically, the austenite phase is related to the Fe-γ phase and corresponds to a face-centered-cubic (fcc) lattice, whereas the martensitic phases can occur in the three lattices: face-centered-tetragonal (fct), bodycentered-tetragonal (bct), and body-centered-cubic (bcc) [11]. The distinction between the crystallographically equivalent lattices fct and bct is based on different variant axes ratios that occur in the tetragonal phase, viz c/a fct ≈0.94 and c/a bct ≈0.717.
To achieve a better understanding of the metastability and the involved transformations between these Fe-Pd phases of homogeneous chemical disorder of Fe and Pd atoms, MD simulations were performed. Those employed customized embedded atom method (EAM) potentials, previously derived to study the impact of lattice imperfections due to ion beam treatment, nanoindentation as well as application of shock waves on transformation properties [12][13][14]. These potentials are targeted to reliably describe the energetics of disordered Fe 7 Pd 3 alloys, as predicted by density functional theory (DFT) based ab initio calculations, correctly from the fcc Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
to the bcc lattice along the Bain path [15], besides ordered and elemental phases. Entel and Gruner et al [16], who performed complementary DFT studies on large Fe-Pd systems obtained Bain path energy landscapes, consistent with our results.
Building upon this, completely based on first principles computer simulations and concominant to the work done by Entel et al [17][18][19] who studied martensitic transformations in Fe and different ferrous alloys with classical MD, we investigated further properties concerning the martensitic transformations on atomic scales. In doing so, we first investigated thermal properties, viz the transition temperature, also to show the reliability of the EAM potentials. Afterwards, in order to achieve stress induced martensitic phase transitions, a Bain type of transformation was externally imposed on the system. Within this context, it is necessary to point out that experimental studies-depending on the boundary conditions-were reported to comply with these results [20] or yield different transformation paths, including the Nishiyama-Wassermann transformation [21]. As a result of the imposed tetragonal deformation, there was the formation of differently orientated martensitic variants observed, accompanied by the formation of internal interfaces. As one of the prerequisite to the FSME, we focused on detailed analysis of these interfaces and especially their mobility, in order to understanding better the physical processes underlying this effect.

Computational details
For the simulations a massively parallel, classical (group maintained) molecular dynamics (MD) code was used [22]. The simulation cell employed periodic boundary conditions in all spatial directions using the classic minimum-image convention. It is initially composed of 20×20×20 fcc four-atom unit cells, yielding 3.2×10 4 atoms, which initially randomly occupy lattice points following the Fe 7 Pd 3 stoichiometry. The positions and velocities of the constituting particles were integrated using a fifth-order Gear predictor-corrector scheme [23] with a time step of 2.28 fs. The pressure, as well as the temperature of the simulation cell were controlled with the Berendsen [24] barostat (bpc) and thermostat (btc) at coupling constants of m t b = =-· 5 10 Pa fs bpc bpc bpc 1 6 and t = 300 fs btc , respectively.

Thermally induced martensite-austenite transformation
The critical temperature for the first order austenite to martensite phase transition lies within a hysteresis interval. Such a hysteresis is commonly characterized by four temperatures, which are the austenite (martensite) start temperature, where upon heating (cooling) the first austenite (martensite) variants form, and the austenite (martensite) finishing temperature, where upon continuation of the heating (cooling) process the transformation from martensite to austenite (vice versa) is completely accomplished (see e.g. [25] for more details). The reasons for this hysteresis behavior are structural contributions to the martensitic transformation energetics, which are nuclear cooperative movements, the formation of internal interfaces, and the accompanied emergence of internal stress fields [26]. Estimations of the critical austenite to martensite transition temperature for the presently studied Fe 7 Pd 3 system is between 270 and 280 K [5,11] in experiments. Within the framework of the performed MD simulations, however, it was not possible to observe a structural transition from fcc to one of the martensitic phases (and vice versa) by decreasing (increasing) the temperature of the simulation cell. Like in previous studies [27], we attribute this to the absence of large fluctuations, due to the short time scales (of only up to μs) covered by MD, which are definitely too short, compared to experimentally observed rate dependencies for spontaneous nucleation on the scale of minutes. Furthermore, transformations in the system are affected by the thermostat and barostat of the MD-code via the realized thermodynamic ensemble as well as geometric boundary conditions on the simulation box (despite being independent for large enough systems in both cases, viz the thermodynamic limit). This together with the interfering structural contributions that induce the hysteresis, leads to the case that the Gibbs free energy barrier for a transformation to occur cannot be overcome. Geometric boundaries presently include the shape of the cell, which is kept rectangular, which might additionally impede transformation by an increase of effective barriers. For Parinello-Rahman type fluctuating cells [28] this effect is expected to be immanent as well, though to a lesser extent.
The approach of choice to avoid these limitations was the method of creating a multi-phase initial system that can propagate towards the equilibrium state that corresponds to the prespecified temperature (inspired by [29]). As starting point we employed a two-phase cell composed of austenite (γ-disordered Fe 7 Pd 3 on a fcc lattice) and bcc martensite (also disordered) that has already been employed by us and described in detail during the development of the EAM potential [12]. The resulting phases after propagation were analyzed with the common radial distribution function (RDF), as shown in the inset of figure 1. The peaks of the RDF-graphs are rather broad. This accounts to the high disorder and the resulting sampling of the varying inter-atomic distances caused by the different atom sizes of Fe and Pd, the different local atomic environments and therefore potentials in the randomly ordered Fe 7 Pd 3 alloy. At 300 K, which is above the expected critical temperature, the equilibrium crystallographic configuration, as deduced from the value of the lattice parameter a, corresponds to a face centered (fcc) lattice structure, whereas for the lower temperature, 250 K, which is below the expected critical temperature, the equilibrium crystallographic configuration corresponds to a body centered (bcc) lattice structure.
For each temperature taken from an interval of integer values around the expected critical temperature, the initial simulation cell was prepared with 50 uncorrelated initial temperature sets by initializing the atoms with different Boltzmann-distributed velocities. The external stress on the simulation cell was set to zero, allowing the structures to relax according to internal stresses. The simulation was stopped after 30 ps, which was observed to be sufficient time for the cell to propagate towards the resulting phase. After that time the configuration of the simulation cell was determined and the relation between the fcc and bcc configuration was calculated by counting the total occurrences of the fcc structure (figure 1). Due to the best fitting behavior, the fit function was chosen to be a sigmoid function.
In the context of the hysteresis behavior (mentioned above), by interpreting the interval from n(fcc)=0 to n(fcc)=1 as the expansion of the hysteresis curve, one obtains a range of ≈15 K. Interpreting the inflection point of this graph as the critical temperature results in a value of ≈281 K.
Within the context of temperature change induced transitions between austenite and martensite in Fe-Pd, it further deserves mentioning that the specific form of disorder, besides contributions from stress effects, e.g. due to the presence of defects in the phases or at their boundaries, affects the phase transition temperature, particularly for small simulation cells. Thus variations in the phase transition temperature between different studies might be intrinsic manifestations of this scenario.
Nonetheless, the critical temperature determined by us is only about 30 K above the transition temperature stated in recent publications [30,31], showing that the completely ab initio-based EAM potentials are suitable for describing physical properties.

Stress induced phase transition
To study structural transformation and mechanical properties of the Fe 7 Pd 3 shape memory alloy a simulation cell in fcc-austenite configuration was probed against uni-axial tetragonal deformation of the initially cubic simulation cell. This setup is inspired by the Bain path [15], which describes the transformations between fcc and bcc lattices by means of characteristic lattice c/a-axes ratios occurring during such a tetragonal deformation. This simulation setup was modeled by a continuous contraction of the c-axis, which was chosen to be in [001]direction. By varying the magnitude of the strain rate it was found that a strain rate of −10 −7 fs −1 , where z 0 is the height of the initial simulation cell, was low enough for the system to adjust in a quasi-static like form to the deformation induced stresses. This constitutes a compromise between the simulation time accessible in largescale classical MD simulations and the velocity of sound comprising an upper bound. At 2.0 ns, the most critical transformation processes (complete conversion from fcc associated lattice to bcc associated lattice) already occurred.
Three different parameter dependencies were sampled for the simulations, viz., the temperature, the initial temperature set, and the atomic order. The temperature, ranging from 260 to 300 K, was increased in 5 K steps. The atomic order sampling was realized by randomly assigning the atom types to the lattice positions according to the Fe 7 Pd 3 stoichiometry. For both the temperature set and the atomic order there were five different values tested per temperature. It was observed that the simulated transformation process can be divided into three stages that exhibit different dominating structural transformation mechanisms. The stages are marked in figure 2(d), showing the mean radial-distribution-separation-function (D RDF ) over all parameters. The Δ RDF is defined as: where r is the inter-atomic distance, ( ) g r t 0 is the common RSF of the system at a given time t 0 , and δt is the time step between two evaluated simulation cell configurations. For each simulation the time series of RDFs were obtained by computing the RDF of the current state every 10 ps≡δt. The idea of the Δ RDF is that large absolute differences of consecutive RDF give information about major local structural changing events that occur in the simulation cell during the transformation. The extent of the locality can even be adjusted by the choice of the RDF cut-off radius r cut or even by setting the entry radius r 0 to values larger than 0. This procedure allows for analysis of structure changing events localized within defined sphere shells. In this work r 0 was 0 and the cut-off radius was 10 Å, which is approximately three lattice constants.
From visual inspection and the help of the Δ RDF , three different stages can be distinguished: Stage 1. The first stage is governed by formation of isolated tetragonal variants of different orientation within the simulation cell. Due to the random alloy character of the atomic arrangement and the resulting high variation of local inter-atomic potentials, the local stress fields do not only have contributions along the z-axis of the simulation cell, in which the strain is exclusively applied, but also along other directions. In order to accommodate these shear stresses the atoms consequently perform a movement that leads to local shear strains. As a result differently oriented tetragonal variants form. In this stage the lattice structure of the simulation cell appears to be contorted ( figure 3(b)), which is due to the absence of defined interfaces between the loosely distributed neighboring orientation domains.
Upon increasing the strains, the tetragonal variants eventually percolate to form structural more favorable martensitic domain patterns, which means that neighboring orientation domains either align or form interfaces between twins ( figure 3(b)). These structural adjustments lead, due to the cooperative shear movement of the concerned atoms, to peaks in the graph of the Δ RDF . The peak magnitudes depend on the number of atoms that are involved in this process. The repeatedly recognizable angle between neighboring variants of the contorted structure ( figure 3(b)) is α≈92.8°. An evaluation that is based on simple geometric considerations under the assumption of a-c twinning with the formula [32] where a/c is the ration between the long tetragonal axis and the short tetragonal axis, yields that the c/a-axes ratio is approximately 0.94, which is characteristic for the fct lattice structure. The occurring peaks, grouped and termed as (A) in figure 2(d), are therefore assumed to represent the formations of larger domains of fct variants. At the end of the first stage there is another group of major peaks, termed (B), in the Δ RDF -graph ( figure 2(d)). They emerge upon the connection of larger orientation domains across the periodic boundaries of the simulation cell and introduce the second stage, accompanied by local major crystallographic changes.
Stage 2. After neighboring domains have connected to form larger domains, that orient to reduce stresses, the simulation cell is, in principle, divided into three martensite domains with different orientations. There is, on the one hand, the dominating bulk martensite in a body centered structure (evaluated with the help of the adaptive common neighbor analysis algorithm [33]). Embedded inside of this, on the other hand, there are two martensite layers, which expand parallel to the (101) fcc and (¯) 011 fcc plane (figures 3(c) and 4(b)). At the boundaries of these embedded layers of tetragonal variants, a characteristic tilt angle of α≈105°is observed, which according to equation (2) corresponds to a a/c-axes ratio of 0.77. Compared to the a/c-axes ratio of the bct phase (0.717), one can conclude that the occurring structures can be related to a bct-martensite configuration.
Further compression of the simulation cell leads to sharper interfaces of the embedded layers, while their thickness decreases. Stage 3. Eventually the thicknesses of the orientation layers vanishes, causing the embedded tetragonal variants to disappear. For both layers, theses processes occur non-simultaneously, which can also be determined by two distinct peaks in the Δ RDF -graph (figure 2(d), peaks in group (C)). The remaining structure is body centered cubic martensite ( figure 3(d)).
The transformation showed the same behavior of the three phases and the embedded layers in the bctstructure for all the tested parameters (figures 2(a)-(c)). Comparing the mean D RDF over all temperature sets and atomic orders graphs for different temperatures one obtains that the peak intensities decreases for increasing temperatures. This may be caused by the lower mobility of the atoms and the resulting abruptness of the cooperative shear movements. From the first-order phase transition character of the martensitic transformations, it is expected to observe a stress dependency of the transition temperature, as predicted by Clausius-Clapeyron type equations [35]. Experimentally, this temperature-stress connection has already been studied and evaluated [12]. Presently, however, the links between the parameters volume and transformation stresses for different temperatures, required for an analysis with the Clausius-Clapeyron equation, are hard to establish. The reason for this is the disorder of the system. The exact shape of the coincidentally transforming areas and the corresponding stress fields are not accessible for the extend of the necessary accuracy of a global Clausius-Clapeyron type analysis covering all cell.
While the graphs of the mean Δ RDF over all temperatures and atomic orders show no significant dependence on the temperature set, there are larger observable differences between the graphs of the mean Δ RDF over all temperatures and temperature sets. This observation shows that there is a dependency of structural details of the transformation on the atomic order of the simulation cell, which can be understood by considering the disorder effects of the atomic environment. Similar to the behavior described by the Clausius-Clapeyron equation, the transition temperature of an alloyed system can be treated described as a function of disorder. The underlying mechanism to this can be traced back to local shear stress fields, which are strongly correlated to the varying disorder potentials. Hence, as already predicted by the Clausius-Clapeyron equation, the transition temperature, and in this context the critical global transformation stress values, are shifted due to the effective local stress fields. A much more detailed discussion of this phenomenon be found in [12].
Compared to the graphs of the thermodynamic variables, e.g. internal energy, volume or stress, only the spatial stress in the z-direction of the simulation cell and the Δ RDF significantly display the fct-forming events taking place at the beginning at the simulation (figures 1(d)/(A)). This shows that the Δ RDF is a helpful tool to identify and distinguish events that occur during the time evolution of a system, where there are local lattice deformations incorporated. Between the spatially resolved strain variables of the system and the Δ RDF , however, there is a close connection but the advantage of latter is the focus on local structural changes.
The occurrence of different tetragonal variants in the simulation cell that appear in the course of the martensitic conversion, and the fact that they are not aligned with the compressed z-axis let us assume that shearing is the dominant mechanism for the transformation. This is typical for the transformation path proposed by Nishiyama and Wassermann [36,37] (figure 4(a)). There the conversion from a fcc to bcc lattice is accomplished by atomic shear movement along the [¯] || [¯] 112 011 fcc bcc direction. Comparing this to the orientation of defect meshes from the dislocation-analysis-algorithm [38] (DXA) one observes similarities. Given a particular lattice, the DXA identifies all defects of this prespecified structure and returns the defect mesh, which separates the crystal regions that belong to the given lattice from those that are different, such as defects or other crystallographic domains. Also for pure iron, as the main ligand of the here studied alloy, the Nishiyama-Wassermann path was found to be energetically preferred over the Bain path for the martensitic transformation [19]. The occurrence of two different orientations can be explained by the absence of a preferred direction due to cubic symmetry.  [36,37] fcc to bcc transformation path. (b) Defect meshes from the dislocationanalysis-algorithm [38] that indicate the boundaries between the bct matrix and embedded bct layers.

Characteristics of the twin boundaries
The twin boundaries were nucleated by applying stress with the constant-strain mode (as described above). At a strain rate of 10 −7 fs −1 , the effect of elongation along the [001]-direction of the simulation cell was modeled. The expansions of the x-and y-dimension of simulation cell were capable of adjusting to occurring stresses so that upon elongation in the z-direction there was an observable compression of the x-y plane. As known from the Fe 7 Pd 3 alloy the preferred tetragonal structure is the one with a short c axis and two long a axes. In the course of the deformation process two far-ordered domains of differently oriented tetragonal variants form at a critical strain in z-direction of approximately 1.3 (inset of figure 5), both with longer a axes, one parallel to the z direction of the simulation cell and short c axes pointing either towards the x-z plane or the x-y plane ( figure 6(a)). At the boundaries, the tetragonal variants coincide under an angle that varies with increasing strain (figure 5). The optimum of the structure that can also be found as a minimum in the curves of the mean potential energy per atom (inset of figure 5) is achieved for an angle that corresponds to a c/a-axes ratio of 0.717, which is the characteristic ratio of the bct structure. At this optimum angle, the simulation was stopped and the configuration with the far-order twin-pattern was used to prepare the simulation for investigating the mobility  of the twin boundaries. The common neighbor analysis also revealed a body-centered configuration for the atoms in the simulation cell.
The setup therefore was to allow the simulation cell to relax, according to internal stresses caused by the twin boundaries. In the course of the relaxation, the boundaries of the inner orientation domain moved towards each other (figures 6(a) and (b)) until they eventually approached and vanished. The duration of this observed process shows a dependency on the coupling constant of the Berendsen barostat m t b = -· bpc bpc bpc 1 : strong coupling, viz. small time constants τ bpc , leads to a faster interface movement ( figure 7). However, it was observed, that neither the temperature, in the range around the working environment of desired applications from 24to 340 K, nor coupling to the temperature control had influence on the mechanisms described below.
In order to investigate the nature of the interface movement, the velocities of the atoms were analyzed ( figure 6). From the cell geometry one finds that the interfaces move along a vector  v i in the x-y plane with either both positive (interface 1) or both negative values (interface 2) for the components v x and v y . The determined average values, on the other hand, show that there is a more pronounced component of the velocity in positive y direction and negative x direction, both inverting signs as the interfaces begin to merge.
From this, one can conclude several properties. At first, one finds that the interface mobility can be described by the well studied mechanism of the propagation of partial dislocations, moving perpendicular to the normal of the twin boundaries (see e.g. [39]). During the merging process, resulting in the vanishing of the interfaces, an annihilation of the antithetically oriented propagating partial dislocations occurs together with a complete annihilation of the two opposite Burgers vectors associated with the twin dislocations. Also during this process the atomic velocities caused by the partial dislocations movement are reflected, which is substantiated by the sign inversion of the spatial components of the mean velocity over all atoms. After this annihilation, the components of the velocities do still not vanish completely. This behavior can be attributed either to be an artifact from the Savitzky Golay smoothing procedure of the highly noisy velocity data, to the release of energy stored in the partial dislocations or to the influence of the highly strained simulation cell by which the simulation was initialized.
From the occurring imbalance of the values of the mean velocities over all atoms in x and y direction as the interfaces approach each other before the merging process (figure 6), one can conclude, that the interface movements are of different intensity, as their contributions do not cancel each other. This leads to the assumption that one interface propagates more actively compared to the other. By tracking the distance between the two interfaces over the time one obtains an approximation for the highest achievable interface mobility, viz for the strongest coupling to the barostat, of ≈250K500 m s −1 , taken into account the non-equal participation behavior, as one of the interfaces move slightly faster. The value for this velocity is about an order magnitude smaller than the velocity of sound in iron or palladium. The observed behavior of an accelerated motion of the twin boundaries that progresses to constant velocity ( figure 7(b)) can be attributed to the phenomenon of dislocation drag, which was exhaustively discussed for Fe 7 Pd 3 systems in [40].
Further characterization of the interfaces was performed by analyzing two quantities of the graphs of the mean internal energy per atom ( figure 7(a)). The first one is the energy difference between the structures with twin boundary interfaces and without interfaces, after they completely vanished. By normalization of this energy to the interface area one obtains a value of 0.205 J m −2 , which is of the same order of magnitude as boundary energy values of other materials that were studied in the framework of EAM potential based MD [41,42]. For twin boundaries occurring in the face-centered crystallographic regime, much smaller interface energies and a much higher interface mobility are expected. However, due to the absence of clearly defined interfaces in the simulation data, no more detailed evaluation of those quantities were possible.
Nanoindentation measurements on acoustic emission in Fe 7 Pd 3 shape memory alloys, were previously performed [43], implying that stress induced reorientation of twin variants can occur abruptly and can be correlated to the measurement of phonons. These results agree with the here observed avalanche-like transformation events that occurred upon deformation and the observation of phononic oscillation due to the release of interface energy.

Conclusion
The good accordance of the estimated martensitic phase transition temperature with experimental values corroborates the reliability of the EAM potentials and strongly includes that they are a reasonable model to investigate properties of the studied disordered Fe 7 Pd 3 alloys on an atomistic scale by means of MD simulations. The here employed radial-distribution-separation-function, presently represents a complementary tool to evaluate time series of MD data with respect to local structure changing events. Further on, the Fe 7 Pd 3 martensitic transformation dynamics revealed that, despite a provided tetragonal transformation path, the system accomplished the fcc to bcc conversion by cooperative shear movements along directions that are closely related to the Nishiyama-Wassermann [36] path. An analysis of the twin boundary movement shows their high mobility with velocities of up to a few hundred ms −1 . This velocity can be understood as a limiting factor for the operating-frequency of the FSME effect. An interface energy for the twin boundaries was determined by a phonon related analysis of oscillations that occurred in the graphs of the mean potential energy per atom. This energy, which is in accordance to previously found interface energy values for iron [44], can serve as a prediction for nano-indention measurements.