Molecular dynamics simulation of nanoindentation of Fe3C and tetrahedral Fe4C

Study of nanomechanical response of iron carbides is important because presence of iron carbides greatly in ﬂ uences the performance and longevity of steel components. This work contributes to the literature by exploring nanoindentation of Fe 3 C and tetrahedral-Fe 4 C using molecular dynamics simulation. The chemical interactions of iron and carbon were described through an analytical bond order inter-atomic potential (ABOP) energy function. The indentations were performed at an indentation speed of 50 m/s and a repeat trial was performed at 5 m/s. Load – displacement ( P – h ) curve for both these carbides showed residual indentation depth and maximum indentation depth ( h f / h max ) ratio to be higher than 0.7 i.e. a circumstance where Oliver and Pharr method was not appropriate to be applied to evaluate the material properties. Alternate evaluation revealed Fe 3 C to be much harder than Fe 4 C. Gibbs free energy of formation and radial distribution function, coupled with state of the average local temperature and von Mises stresses indicate the formation of a new phase of iron-carbide. Formation of this newer phase was found to be due to deviatoric strain rather than the high temperature induced in the substrate during nanoindentation. & 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).


Introduction
Iron carbides are prime candidates for sensors, magnets, catalysts, and alternative raw materials used for steel production and for carbon nanotube growth by the pyrolysis of organometallic precursors [1]. Use of iron-carbides also find applications in iron based Fischer-Tropsch catalyst process [2], geophysics [3] and synthesis of nanograined iron carbide crystals through the principal of severe plastic deformation techniques such as equal channel angular pressing (ECAP) and high pressure torsion (or compression shear) [4]. Some of the major phases of iron carbide identified till date are shown in Table 1, which also highlights the constituents and arrangement of these phases and their respective mechanical properties.
It is apparent from Table 1 that a small change in the microstructure, especially the presence of Fe 3 C (cementite), can lead to significant changes in the bulk properties of steels. Also, low carbon concentration is known to strengthen bulk iron while manufacturing steels, but a higher amount of carbon beyond a critical measure and the subsequent formation of cementite-phase iron carbide can have the opposite effect.
Cementite (6.67% carbon) is one of the main iron carbides among all the steel alloys and is known to be hard and brittle. A scratch test showed that the track length of a groove made on pearlite phase was smaller than that made on the ferrite phase [6]. This difference was attributed to the relative higher hardness of pearlite as it contains lamellar cementite (ten times harder than ferrite). In accordance with this experimental observation, finite element analysis (FEM) also confirmed that ductile ferrite offers a larger chip thickness than pearlite [7]. Thus, cementite tends to influence the mechanical, structural, and thermal properties of steel. However, even though it is known that cementite influences the properties of steel during service life of the components at high pressures and temperatures such as in turbine blades and reactor vessels, its mechanical behaviour at a fundamental atomic level is not known [8]. An improved understanding of the nanomechanical response of cementite at nanoscale can help improve the design of newer class of steels.
Chemically and structurally, cohenite (found in meteorites) and cementite (found in steels) are the same phases of iron carbide, having an orthorhombic unit cell which belongs to the space group P nma. The unit cell of cementite contains twelve iron atoms and four carbon atoms, with the positions of the carbon atoms lying interstitially in a lattice of nearly close packed iron atoms [9]. Cementite has a Poisson 0 s ratio of about 0.36 and its elastic constants are reported to be C 11 ¼ 388 GPa, C 12 ¼156 GPa, C 13 ¼164 GPa, C 22 ¼ 345 GPa, C 23 ¼162 GPa, C 33 ¼322 GPa, C 44 ¼15 GPa, C 55 ¼134 GPa, and C 66 ¼134 GPa [10]. Other aspects known about cementite are that it decomposes between 1100 and 1200 K [11] into liquid state. Beyond the eutectic point at 1420 K, liquid consisting of Fe and C solidifies to form austenite and cementite. While cementite is known to be chemically metastable with a positive heat of formation, its melting point has been theoretically estimated to be around 1425 K [10].
During a ball milling test, amorphous Fe 3 C (which has a crystallization temperature of 593 K) was suggested to undergo oxidation [12]. It was proposed that carbon atoms from Fe 3 C are released (in the temperature range of 773-1273 K) to form either CO or CO 2 . Another study reports that the low-nickel metastable austenitic stainless steel undergoes plastic deformation induced phase transformation from γ austenite (fcc) to α 0 -martensite (bcc) and ε-martensite (hcp). It was observed that the formation of ε-martensite occurs due to randomly spaced overlapping stacking faults while α 0 -martensite forms at shear band intersections [13].
In this study, the focus is to understand the behaviour of nanocrystalline Fe 3 C under high temperature and pressure. Several experimental studies have been undertaken to understand the outcome of Fe 3 C under high temperature and high pressure conditions but at macroscale. X-ray emission spectroscopy in a diamond anvil cell has shown that crystalline cementite undergoes a magnetic transition at a pressure of about 25 GPa [14]. This transition is predicted to be a second order phase transition which does not involve any structural changes during the process.
In another study on a polycrystalline cementite, high pressure torsion (HPT) led to the dissolution of polycrystalline Fe 3 C in the form of carbon rich layers at the interface of ferrite grain boundaries at a relatively low pressure of 5 GPa. From several fatigue tests, Ferro et al. [15] concluded that anchoring of dislocations due to carbon interstitials is unstable and thus dislocation nucleation is independent of the interstitial sites. During compression tests performed at 300 K, no phase transformation of Fe 3 C was observed upto a pressure of 30 GPa by Li et al. [16]. Later, Scott et al. [9] used Synchrotron-based X Ray diffraction in conjunction with laser heating to reveal that cementite did not undergo phase transitions at higher limits of 1500 K and 73 GPa, thus confirming the previous findings of Li et al. Rouquette et al. [17] also performed similar studies to reveal that Fe 3 C is stable up to 70 GPa and temperatures as high as 3400 K.
Since the melting point of Fe 3 C is lower (1425 K) than 3400 K, therefore the fact that it retained its phase at such a high temperature needs to be explored better. Therefore, the phenomenon of Fe 3 C retaining its phase at temperatures higher than its melting point was one of the key motivations of the current work. Nanoindentation of a material is an appropriate tool to simultaneously generate the conditions of high temperature and high hydrostatic pressure. While the studies reviewed above provide some experimental insights, no theoretical study is ostensibly evident in the literature which would suggest the nanomechanical response of cementite under high hydrostatic pressure and high temperature conditions [9,16]. This study aims to bridge this gap by examining the nanoindentation [18] behaviour of Fe 3 C and investigating if the reported macroscale experimental studies on compression behaviour of cementite can be extrapolated to nanoscale behaviour of cementite during nanoindentation. To compare and contrast the nanomechanical response of cementite with another iron carbide, tetrahedral Fe 4 C was chosen as the reference material as it was found to be more stable than other iron carbides such as FeC, Fe 7 C 3 , Fe 5 C 2 and Fe 3 C in that order [8].

Molecular dynamics simulation of nanoindentation
Compared to a conventional experimental procedure such as Vickers indentation and measurement of the imprint diagonals, instrumented indentation provides a more accurate assessment through recording of the load-displacement data, which enables experimenters to study discrete events such as phase transformations, shear instability and nucleation of dislocation [19] besides measuring the indentation energy and the elastic modulus of the material. Moreover, experimental studies are often cumbersome to be carried out with a high degree of precision especially at the scale of several nanometres in the presence of oxide layers and size effects. Furthermore, it has also been pointed that Kick 0 s Table 1 Various iron carbide alloys [5].  [13]. Molecular dynamics (MD) simulation has been developed and widely used as an alternative intermediate tool to bridge the gap of first principles, experimental and finite element methods. It is an appropriate tool to understand the atomistic tribology of simultaneously occurring processes, the foremost of which are structural transformations in the material, nucleation and propagation of dislocations and chemistry involved during the mechanics of the process. Therefore, to reveal a more phenomenological understanding of the process of nanoindentation of Fe-C system, this paper adopts MD simulation.

MD simulation model
In this work, the "Large-scale atomic/molecular massively parallel simulator" (LAMMPS) [29] was used to perform a series of MD simulations. VMD [30] and Ovito [31] were used to visualize and analyze the atomistic simulation data. The general aspects of simulation of a nanoindenation are quite similar to simulation of nanometric cutting which are described elsewhere [20,21] except for the modifications made in the boundary conditions, shown in Fig. 1.
The simulation model shown in Fig. 1 assumes periodic boundary conditions (pbc) along the Z direction of the model while thermostatic layers and fixed layers of atoms were considered along the X and Y axes. Atoms in the central region directly affected by the interaction of the indenter and the workpiece were allowed to follow Newtonian dynamics (LAMMPS NVE dynamics), while atoms in a thin boundary layer were subjected to a thermostat (LAMMPS NVT dynamics) to dissipate the generated heat in the finite simulation volume. Use of such an ensemble ensures that the deformation mechanics is not affected by any artificial dynamics and is in accordance with the MD simulations performed and reported in the literature [22,23]. In this simulation model, atoms in the indenter were kept fixed (indenter was assumed to be an infinitely rigid body since the intent was to understand the mechanics of the substrate and not of the indenter). This assumption once again is in accordance with previously performed simulations [24-26].

Potential energy function
The potential energy function to describe the atomic interactions governs the accuracy of a molecular dynamic simulation which in turn defines the reliability of the simulation results. Indeed, the potential energy function in molecular dynamics can be considered analogous to the constitutive equations used in the FEM simulation model in that accurate description of material defined by the function provide better and meaningful simulation results. The analytical bond order potential (ABOP) proposed by Henriksson et al. [11] (which is a modified version of Tersoff [27] potential energy function) describes Fe 3 C and tetrahedral Fe 4 C appropriately, and was used to specify the interactions between Fe and C atoms in this simulation model. The potential function was implemented in LAMMPS in the following manner: where V is the site potential energy. The sub-function V ij describes the energy between two atoms i and j, where i and j denote the atoms of the system, and r ij is the distance between these two atoms.

Table 2
Potential function parameters used in this study [11].
The Morse type terms and bond orders are given by where V R represents a repulsive pair potential, V A represents an attractive pair potential, b ij is the bond order term and θ ijk is the bond angle between the bonds i-j and i-k. From its respective reference, the potential function parameters of ABOP used in the current simulation are listed in Table 2 1

.4. Simulation parameters
The lattice parameters used in the simulation considering using the periodic boundary conditions were a¼ 0.5086 nm, b¼ 0.6521 nm and c¼0.445 nm for Fe 3 C and a¼ b¼c¼0.374 nm for tetrahedral Fe 4 C. These parameter settings give the zero pressure cell volume of 0.1476 nm 3 for Fe 3 C and 0.052 nm 3 for tetrahedral Fe 4 C, in close agreement with previously reported values [1,9]. The lattice structures obtained using these parameters (unit cell) are shown in Table 3, which also highlights the respective positions of iron and carbon atoms. Further details and geometry specifications of the simulation are shown in Table 4.
It may be noted here that usually nanoindentation experiments are executed at finite speeds of few mm/sec in contrast to the speed of 50 m/s used in the current simulation making it somewhat the case of a nanoimpact rather than a nanoindentation. This procedure is however in line with the inception of such studies [22,28,29], and has been adopted so as to mitigate the current limitations of computational time. We take the conventional wisdom to use high speed in the current simulation despite knowing that it might cause some spurious effects. To minimize these effects, sufficiently larger size substrate is used. The choice of substrate needed a large scale molecular dynamics simulation on a high power computing system and the system situated at Daresbury, UK was utilized for the same.
Normally, in the hardness field, pyramidal indenters, such as Vicker 0 s, Knoop, and various conical indenters, are classified as "sharp" indenters while spherical indenters are referred to as Table 3 Unit cell of Fe 3 C and tetrahedral Fe 4 C (Pink colour represents iron atoms while cyan colour represents carbon atoms) with a,b,c as the equilibrium lattice parameters.

Material
Number of atoms in one unit cell Unit cell/Lattice structure Top view at (1 0 0)  Readers are requested to refer the web based version of this article for correct interpretation of the colour legends. ''blunt' 0 indenters [30]. Even the sharp indenters have some finite edge radius and therefore, to give a finite radius to the indenter such as those used in practical cases (despite being referred to as extremely sharp), a spherical shaped rigid indenter was used during the simulation.
To quantify the plastic deformation of the atoms, an algorithm proposed by Shimizu et al. [31] was chosen. Atomic local shear strain (von Mises strain) for each atom i was computed by comparing the two atomic configurations (during indentation and the starting configuration) using OVITO. This was accomplished by using a local deformation matrix J i to first compute the local Lagrangian strain matrix as which was then used to compute the von Mises shear strain for each atom i. The von Mises shear strain η i mises , recognized as a good measure of the local inelastic deformation, can be expressed as follows:

Load-displacement (P-h) curve
During displacement controlled indentations (depth-sensing approach), an indenter is pressed against a surface up to a specific depth in the substrate. The data of the load applied to the indenter and the displacement of the indenter is then plotted on the ordinate axis and abscissa axis respectively. Oliver and Pharr [32] method is used to calculate the projected contact area from the analysis of the load-displacement (P-h) curve and the depth. This is done by drawing an imaginary line by typically following the power law on the top 1/3 rd part along the unloading curve. Oliver and Pharr method thus enables to estimate the contact area directly from the P-h curve without assessment of the actual contact area of the indenter curve to consequently estimate material properties such as hardness and elastic modulus. The contact area evaluation is therefore a key parameter in assessment of material properties during nanoindentation.
In a subsequent work, Bolshakov and Pharr [33] asserted that during nanoindentation of sink-in materials or in cases where pile up of material is predominant (specifically when the ratio of residual indentation depth and maximum indentation depth (h f / h max ) exceeds the value of 0.7), the material does not work harden considerably. In such cases, the contact area deduced from the conventional approach of Oliver and Pharr underestimates the true contact area to the tune of up to 60%. Consequently, the conventional method of Oliver and Pharr method overestimates the value of the indentation and elastic modulus by the same margin. Therefore, pile up material must be accounted for obtaining correct estimates of the hardness and elastic modulus. In the past, Christopher et al. [34] used MD simulation to highlight pileup of material during indentation of iron and silver specimens. This was also the case observed in this work (discussed in details later). In what follows, the P-h curves of nanocrystalline Fe 3 C and tetrahedral Fe 4 C obtained from the simulation are detailed.
The P-h plot obtained for nanoindentation of Fe 3 C and Fe 4 C from the current simulation results are shown in Fig. 2(a) and (b) respectively. It can be seen from these plots that the h f /h max ratio in both cases were observed to be higher than 0.7 i.e. 0.76 and 0.84 for Fe 3 C and Fe 4 C respectively. This suggests that the conventional approach of Oliver and Pharr will give overestimated values of hardness from these plots. To mitigate this problem, another parameter referred to as "contact atoms" proposed by Shangda and Fujiu [35] was used.
Shangda and Fujiu 0 s method is extremely simple and can be easily applied to MD simulation based exploration of nanoindentation. As per this method, which aligns with our simulation as the diamond tip was considered a rigid body, the projection of all contact atoms will approximate to a circle. One can easily work out the immediate contact area through a simple expression of Area¼ π(Rþ r 0 ) 2 where R is the radius of the indenter and r 0 is the cutoff radius between the indenter and the substrate material (0.1477 nm here). The reason for addition of r 0 with R is that every peripheral atom on the indenter will have an action range of upto r 0 which will be repelled by the atoms of the specimen being indented. Therefore, when this parameter is added with the indenter radius, it gives the overall radius for calculation of the total contact. This method directly explains that a change in the radius of the tip results in a proportionate change in the contact area between the indenter and the substrate.
Some other unique features associated with the P-h curve were also observed. Firstly, the indenter was observed to experience the cohesive forces even before reaching the surface of the substrate. This is shown in the plot as jump out due to attraction. These jump outs led to a sudden increase in the local temperature, which is  later. From the simulation video, cohesion between the atomic species of the stylus and the substrate was observed to be responsible for these forces. Secondly, unlike Fe 4 C, where some noticeable deviations were observed, the forces plotted during nanoindentation of Fe 3 C showed only a monotonic increase until reaching peak load of 306.88 nN at the indentation depth of 1.5 nm (2.57-1.068 nm). The value of the peak load (upto the same depth of indentation) for Fe 4 C was observed to be far lower than that for Fe 3 C suggesting that Fe 3 C is much harder than Fe 4 C. Moreover, the plot for Fe 4 C is monotonic in the negative area. So while loading, the plot is monotonic for Fe 3 C, but while unloading the plot is monotonic for Fe 4 C. Also, only in Fe 3 C a sharp kink was observed during unloading and not in Fe 4 C. Another feature of these P-h curves was that the unloading curve of Fe 4 C was found to follow the power law whereas Fe 3 C showed a deviation. This deviation also accounts for the difference in the measurement of residual depth (h f ) in the two carbides. Table 5 shows a comprehensive list of all the relevant results for both the simulations.
The average loading forces were estimated to be 180 nN for Fe 3 C and 103 nN for Fe 4 C. Consequently, the nanoindentation hardness (average forces/contact area) for the two carbides Fe 3 C and Fe 4 C were obtained as 23.16 GPa and 13.26 GPa respectively. The magnitude of peak compressive stress was observed to be upto 10 times the magnitude of peak shear stress signifying that nanoindentation is a compression dominated process rather than a shear dominated process. Compression can certainly bring a significant change in the bond length which is primarily responsible for the transformation in the crystal structure of the material as has been shown in the last section.
The negative depth hysteresis observed from the P-h plot has not received much of attention in the past and thus it was of particular interest. Similar to the forces observed during loading stage before reaching to the substrate surface due to cohesion, a strong cohesion between the substrate atoms and the diamond indenter was seen during retraction of the indenter. The work performed by the cohesive energy was against the direction of the retraction force on the indenter until the ultimate separation which led to the observed depth hysteresis. Consequently, a few atoms attached to the indenter undergo a process of tension being pulled by the indenter from one end and by the substrate from the other end during unloading. This results in necking and formation of the piping like structure before the ultimate separation. These atoms from the substrate can be observed to cling with the atoms of the indenter on one end while others resulted in the pile-up on the substrate as shown in Fig. 3(a) and (b) respectively. Overall, a high degree of pile-up was observed in the form of protruding atoms at such a small depth of indentation and consequently lends support to the inapplicability of the Oliver and Pharr method for the current investigation. Fig. 4(a) and (b) show the snapshots from the MD simulation which were processed using OVITO to highlight the volumetric strain rate in Fe 3 C and Fe 4 C respectively at peak loading. The cutoff radius of 3.0 Å and 4.0 Å were found good for atomic strain analysis of Fe 3 C and Fe 4 C respectively. Fe 3 C showed a higher volumetric strain rate of upto 0.3 while Fe 4 C showed a peak strain of 0.25.

Volumetric strain
The propensity of the strain field and the shape of the strain field in these two materials were investigated further. Fe 3 C showed a larger strain field exhibiting a hemispherical shape which is the well documented shape of the Hertzian strain field observed during the case of nanoindentation [36] and has also become the basis of exploration of ductility in brittle materials. In the case of Fe 4 C, the strain field was observed to be relatively smaller, exhibiting the shape of a truncated cone.

von Mises stress and temperature
In the past, MD simulation studies have not explored the individual contribution of stresses and temperature in influencing the material 0 s behaviour during contact loading such as nanometric cutting and nanoindentation. As shown earlier, nanoindentation is a compression dominated process whereas nanometric cutting is a shear dominated process [37,38]. However, high temperature is common in both processes. In the current work, von Mises stress and the average local temperature in the nanoindentation zone were plotted on the same plot in an attempt to understand the influence of the stresses and temperature on the behaviour of Fe 3 C and Fe 4 C. These plots are shown in Fig. 5(a) and (b) respectively.
One of the contrasting features observed from the plots was a sudden increase of temperature when the indenter approaches the substrate. This jump was also noticed in the P-h curve and can be attributed to the cohesion between the indenter and the substrate. Consequently, a peak temperature of 550 K and 825 K were observed in Fe 3 C and Fe 4 C respectively due to this cohesion. This suggests that Fe 4 C releases more bond-breaking energy (in the form of heat) than Fe 3 C in breaking and reformation of bonds due to the cohesion between substrate and the indenter. As the indenter moved into the substrate, the temperature and the von Mises stresses were both observed to increase gradually. At peak load, the peak stresses were observed to be about 28 GPa for Fe 3 C and 22 GPa for Fe 4 C. These peak stresses correlate with the P-h plot i.e. higher hardness of Fe 3 C compared to that of Fe 4 C aligns with the observation of higher peak stress required for the deformation of the Fe 3 C substrate. An interesting observation from Fig. 5 is that the peaks of stress and temperature do not coincide with each other i.e. at peak loading, temperature lags the peak stress and during unloading, stress lags the temperature. Therefore, the peak temperature at peak stress reached only up to a maximum of 450 K.

Structural transformations
The radial distribution function (RDF) g(r) is the primary linkage between macroscopic thermodynamic properties of the  [35] π Â ((2.85/2) þ 0.147) 2 7.77 nm 2 π Â ((2.85/2) þ0.147) 2 ¼7.77 nm 2 Average hardness (GPa) 23. 16 13.26 Overall peak temperature during entire simulation run (during loading) (K) 725 825 Average peak temperature in the indentation zone (K) at peak stress 450 420 Average peak von Mises stress in the indentation zone (GPa) 28 22 Average peak compressive stress (s yy ) in the indentation zone (GPa) 40 28.8 Average peak shear stress (τ xy ) in the indentation zone (GPa) 4 4.5 material and the intermolecular interactions. In the past, g(r) has been used as a very effective measure to indicate structural transformations in the materials especially during contact loading [39,40]. Similar to these papers, g(r) was used in this work to monitor the variation in the interatomic bond length between carbon and iron atoms of Fe 3 C and Fe 4 C, before, upon, and after the loading. The results obtained from the MD simulation during nanoindentation of Fe 3 C and Fe 4 C during peak loading and upon unloading are compared in the graphs with the equilibrium structure shown in Figs. 6(a) and (b). It can be seen from Fig. 6(a) that before loading, Fe 3 C exhibited a crystalline peak at an interatomic bond length of 0.195 nm with a corresponding second peak at 0.275 nm. At peak loading, these peaks didn 0 t show any noticeable changes except in the count of the number of atoms. However, an additional peak was found to form at about 0.175 nm. Also, the long range peaks (beyond 0.3 nm) of Fe 3 C appeared to have become disordered in contrast to their periodicity before the indentation. A similar phenomenon was observed during the indentation of Fe 4 C as shown in Fig. 6(b). As evident from Fig. 6(b), the initial configuration of Fe 4 C started from a peak interatomic bond length of 0.175 nm. This interatomic bond length was found to change to 0.195 nm upon peak loading and remained unchanged even upon complete unloading. The pattern of the high order peaks becoming disordered was similar to patterns observed during indentation of Fe 3 C beyond 0.3 nm.

Possible new phase of iron carbide
It is interesting to note that Scott et al. [9] in their X-ray diffraction studies on cementite noted a discontinuity in the trend of the lattice parameters at about 30 GPa but did not offer any explanation of this variation. During the nanoindentation simulation of Fe 3 C in the current work, peak von Mises stresses were fairly close (28 GPa) to the experimental pressure of 30 GPa where a discontinuity in the lattice parameter was observed by Scott et al. [9]. It appears that there exists an identifiable newer phase of iron carbide exhibiting a peak interatomic bond length of 0.195 nm with a second order peak at about 0.275 nm with other long range order peaks being highly disordered as evident from Figs. 6(a) and (b). This observation aligns with the fact that this newer phase could also be an outcome of the magnetic transition, if 25 GPa of stress was observed to be the ultimate limit for such transitions [14].
In order to assert whether the occurrence of this possible new phase of iron carbide was a consequence of the peak pressure or peak temperature, Gibb 0 s free energy of formation of Fe 3 C was plotted using the various thermodynamic equations reported in the literature (Fig. 7). The transition in Gibb 0 s free energy can be seen to occur at 1281 K. This means that a reaction to decompose Fe 3 C from its original phase at a temperature lower than 1281 K will not be spontaneous and will be disfavoured. The peak temperature observed during the indentation for both Fe 3 C and for Fe 4 C (o500 K) was far smaller in comparison to temperature of 1281 K leading to thermodynamically favourable decomposition of Fe 3 C. This finding suggests that the formation of the new phase of iron carbide is likely to be an outcome of the high pressure phase transition rather than high temperature transition. The former brings volumetric strain induced changes in the interatomic bond length of iron carbides (Fe 3 C and Fe 4 C).
While this study and the reported simulation results provide some basic inputs on a newer phase of iron carbide, its crystal structure is yet to be fully discovered and can better be examined experimentally using an X-ray diffraction study or some similar methods.

Sensitivity of MD results to the velocity of indenter
A constraint of MD simulation is that the computational time and resources available currently do not permit scaling of the simulation parameters to the experimental scale. Not only the length scales are restricted to up to few nanometres (or few million of atoms) but also the velocity of indentation or velocity of cutting such as the one used in the current simulation (50 m/s) against experimental speed of upto $ 1 m/s in nanometric cutting or 0.1-10 mm/s during nanoindentation is a factor that might lead to incorrect deductions. A higher rate of loading and unloading such as the one used in this work (50 m/s) may also cause a higher strain rate in comparison to the experiments. It is therefore necessary to assess the effect of such strain rates on the sensitivity of the results obtained from the MD simulation. Additional trials were therefore done at different indentation speeds (slower speeds), the results of which are listed in Table 6. Fig. 8 shows the P-h plots obtained from the nanoindentation trials of Fe 3 C and Fe 4 C, using the same geometry as in earlier trials, but at a different indentation speed (5 m/s) and retraction speed (20 m/s). The number of crests and troughs exhibited by the loading curve was observed to be much higher during the indentation made Hertzian strain field of hemispherical shape can be seen to surround the indentation zone with strains getting reduced gradually towards the periphery.(b) Top view of the substrate and side view to highlight local von Mises strain at peak loading of the indenter during indentation made in Fe 4 C. A classical Hertzian strain field of truncated cone shape can be seen to surround the indentation zone with strains getting reduced gradually towards the periphery. at slower speed. From Table 6, the peak stresses and the peak temperature were observed to be similar at different indentation speeds (variation in the temperature was unexpected and this is an area which will require further research). Table 6 however also shows that the plastic depth of indentation and its ratio to the maximum indentation depth both were same confirming that the proposed new phase is not an outcome of the high speed of the indentation used during the simulation. Relying on the accuracy of the ABOP potential function and by neglecting the anisotropy of the specimen [24], these findings thus seem to assert that the MD    simulation results obtained above are insensitive to the speed of indentation of 50 m/s. However, to supplement the above findings, further analysis was done to investigate the strain rate sensitivity (only during loading). Stress (s) during the indentation can be calculated as the ratio of the instantaneous load and the projected area (hardness) while the indentation strain rate _ ε(instantaneous descent rate of the indenter divided by that depth) and m (strain rate sensitivity) are two other important parameters expressed as follows [42]: where h i is the displacement of the indenter at ith time step, dh/ dt is the velocity of the indenter (since here it is a constant displacement indentation), A is the projected area (2π Â R Â h i for spherical indenter) and m is the strain rate sensitivity. Strain rate sensitivity is an important indicator to assert if a large variation in the indentation speed will lead to significant changes in the outcome of the indentation results. Fig. 9 and Fig. 10 show the plots of the log stress vs. log strain rate obtained from the indentation data (Table 6). These plots were fitted with a linear trend line to estimate the slope of each plot. The value of the slope (strain rate sensitivity) obtained from Fig. 9 for Fe 3 C was in the range of À 1.16 to À 1.2 whereas that for Fe 4 C (Fig. 10) was observed to vary between À 1.08 and À 1.589 across the speed of indentation of 5-50 m/s. The variations in the magnitude of strain rate sensitivity are marginal, and therefore signify that the MD simulation results obtained in this work are insensitive to the speed of the indenter.

Conclusions
Despite a rich literature providing a depth of knowledge and understanding on the properties of steel and iron carbides, there are several gaps in the current pool of knowledge of high-pressure surface science of iron carbides. For example, the nanomechanical response of the Fe 3 C (cementite) and Fe 4 C phase during nanoindentation has not been theoretically investigated. MD simulation is used in this work to explore this new phenomenon. The following conclusions can be drawn based on the aforementioned discussions: 1. P-h curve for both the simulated materials revealed that the ratio of residual indentation depth and maximum indentation depth (h f /h max ) exceeded the value of 0.7. Consequently, the conventional approach of Oliver and Pharr method was inapplicable to evaluate the material properties. A modified method was used to calculate the contact area of the atoms which revealed Fe 3 C to be much harder than Fe 4 C. 2. Atoms from the substrate were found to cling to the indenter upon unloading causing a negative depth hysteresis in the P-h curve. This was due to the strong cohesive dynamics of the atoms of iron and carbon. Thus, studying the phenomena of cohesive dynamics was realized to be one of the major strengths of MD simulation method used in this work.

The measure of von Mises strain and the von Mises stresses
were both found to align with the P-h curve suggesting that Fe 3 C requires higher energy to deform than does Fe 4 C. Interestingly, the propensity and shape of the Hertzian field in these two materials were found to be significantly different i.e. Fe 3 C showed a larger strain field exhibiting a hemispherical shape while Fe 4 C showed a relatively smaller strain field exhibiting the shape of a truncated cone. 4. Gibbs free energy of formation and radial distribution function coupled with state of the temperature and stresses present the possibility of formation of a newer phase of iron-carbide during nanoindentation of either Fe 3 C or Fe 4 C. It was found that the formation of this newer phase is an outcome of the deviatoric strain rather than the high temperature induced in the substrate during nanoindentation. Based on the previous literature, this new phase is proposed to have lost the magnetic properties which were inherent in the pristine material. Fig. 9. Comparison of log (stress) and log (strain rate) for Fe 3 C obtained from the MD simulation with respect to change in the indentation speed.