Molecular Dynamics Simulation of the Nanoindentation Process in Cr / CrN and ( Cr / CrN ) 2 Thin Films

Molecular dynamics (MD) simulations were carries out for studying the influence of nanoindentation in the atomistic deformation mechanisms of Cr/CrN and (Cr/CrN)2 coatings with BCC and FCC crystalline structures for Cr and CrN, respectively. The Morse potential was employed in order to determine the atomic interaction forces of the Cr-Cr and Cr-N atoms. A non-deformable potential solid sphere was implemented for determining the role of the nanoindenter. The OliverParr method (OP) was used to obtain the hardness and elastic modulus of the Cr/CrN and (Cr/CrN)2 layers, resulting in values of 18 and 20 GPa for Cr/CrN and (Cr/CrN)2, respectively. The Cheng method was used for correcting the hardness values obtained by the OP method. The Cheng correction showed higher hardness values since it avoids the influence of the scale effect. Regarding the elasticity modulus, Cr/CrN and (Cr/CrN)2 exhibited values of 217.86 GPa and 258.9 GPa, respectively. Simulations of the temperature influence on the hardness were carried out over a range of 300-1000 K. Results indicate that the hardness decreased as a function of the temperature. 4618 S. Amaya-Roncancio et al.


Introduction
In the last decades, thin films have been widely studied because of their several technological applications as hard protective coatings, magnetic storage media and micro/nanoelectromechanical systems (MEMS/NEMS) [1,2].For instance, in the case of hard protective coatings, the mechanical behavior of films during practical applications is different than the performance of the bulk material.For this reason, analysis of thin films mechanical properties is required in order to improve their understanding and reliability [3].The most used technique for studying mechanical properties is nanoidentation.This method allows measuring thin film properties such as the elastic modulus and hardness.In the literature, it is possible to find summaries of nanoindentation method [4,5] and many works presenting nanoindentation experiments for determining mechanical properties of metallic thin films [6][7][8]; nevertheless, deeply studies of the deformation behavior mechanisms are still needed for understanding well the nanoindentation experiments and simulations.Depending on the purpose of the simulation, the studies are concentrated in different features of the nanoindentation procedure.Each numerical model has to be adapted to specific system.Generally, the comparison between results obtained from different simulation models cannot be done directly [9].For instance, finite elements is the most used method that allows a direct comparison with experimental results because of their computational efficiency.This method is useful for determining important parameters of the nanoindentation process; furthermore, it allows obtaining several material parameters using dimensional analysis combined with continuum numerical models.Finite elements is also very useful for microscale studies [10]; nevertheless, when atomistic studies are required for deeper understanding of the plastic deformations, models including smaller scales have to be employed.Because of the decrease in the scales, more computationally expensive processes for the simulations have to be developed and more powerful computational equipment are required; it implies a greater quantity of calculations compared with continuum methods [11].When atomistic details of plastic deformation and the dislocation processes are in the interest of the simulations, numerical models for atomistic simulations have to be developed at the lowest scale: Normally, models at atomic level are employed for systems with sized in the order of 100 nm.On the other hand, one of the most important advantages of the molecular dynamics method is the possibility of identifying the main physical variable tendencies of a certain problem, giving relevant qualitative information in order to understand the physics behind the complex phenomena observing in experimental results; one of these phenomena are the size effect, difficult to study by coarser models.The molecular dynamics (MD) method is considered as one of the most powerful computer simulation method for describing atomic-scale systems.For instance, systems as multilayered metal and ceramic thin films simulated at atomic level are very useful for obtained important details regarding to mechanical behavior of each coating and interfaces, not observed employing continuum models [12].MD simulations of the nanoindentation process of single crystalline thin films have been conducted according to the literature [13].Simulations of nanocrystalline metals was also observed in multilayered metal thin films with thickness in the order of few nanometers [14,15].For instance, Peng et al. [16] used an MD simulation to obtain a relationship between the nanoindentation process parameters of an aluminum thin film on a silicon substrate.Their results indicate that the film hardness increases with the speed of indentation.C.L. Liu et al. [17] simulated the nanoindentation of diamond and gold films in order to analyze the influence of load, loading rate and temperature on the indentation.Fang et al. [18] studied the deformation effect in the nanoindentation process of Al/Ni multilayer films; results exhibited an increase in the depth of the indentation as a function of the maximum load, energy plastic and adhesion.In the literature, no works were found that presented simulations of Cr and CrN as bilayers and multilayers to understand their mechanical properties by accounting for the different crystalline structures of both layers and the mismatch at the interface.On the other hand, experimentally, Cr and CrN materials have reached an important interest as thin films since they can be used for several technological applications [19,20] The aim of this work is to study the mechanical properties of Cr/CrN and (Cr/CrN)2 coatings using molecular dynamics simulations combined with the Oliver-Parr method and the Cheng correction.The Morse interatomic potential will be used for describing the interaction between the atoms of the system.

Simulations Model
Systems of one bilayer and two bilayers of Cr/CrN and (Cr/CrN)2 were built.In this design, the sample size in the x-and y-directions was 86.7 Å.The thickness of each layer (z-direction) is remained fixed at 44 Å for the Cr/CrN and 22 Å for (Cr/CrN)2.In the first part of the simulations, the temperature was kept at 300 K. Performing simulations under these conditions is important for comparing the results of hardness with those obtained before for Cr and CrN [21].The system was composed of 115926 atoms, which is a sufficient number for describing the behavior of the system under study [18].In the second part of the simulations, the temperature was varied between 300 K and 1000 K to observe the influence of the temperature on the mechanical properties of Cr/CrN and (Cr/CrN)2.The Cr and CrN samples exhibit BCC and FCC crystalline structures, respectively.This structure is oriented in the plane (100) and placed on a substrate of silicon oriented in the [100] direction.Both the Cr and CrN layers are built with the same dimensions and crystal orientation.The indenter has a radius of 3 nm and is located on the surface of the sample at a position of 43.35 Å in xand y-axes.The sample is thermally stable to the temperature of work using the velocity scaling method [22].During program execution, the sample temperature was controlled using the stochastic collision method [23].The boundary conditions in the sample are fixed at the edges (x-y plane) and are free at the surface.The indenter was modeled as a rigid sphere following the potential described by the equation: (  ) = ( −   ) 3  (1 where K is a constant related to the effective stiffness of the indenter (10 eV/ Å 3 ) and R is the radius of the indenter (3 nm).The value of K is similar to that used in the research presented by Lilleodden [24].The indenter stiffness values below 1 eV/ Å 3 are not suitable for the deformation of a sample [25].A similar behavior occurs when varying the radius of the sphere: for values less than 1 nm, the deformation of the sample is not reliable because this size is comparable to the size of the unit cell of the nanoindenter and the sample, thus affecting the results of the load-unload curves and producing noise [24].The Morse interatomic potential [18] was used to describe the interaction between the atoms of the system, as described by the equation: where Dij is the cohesion energy, which is the energy value for the equilibrium point between two bonded atoms; rij is the instantaneous distance between atom i and atom j; ro is the equilibrium distance; and  is the elastic modulus.The parameters used in the simulations, including the lattice parameters, are listed in Table 1.For a pair potential as, removing an atom from the lattice involves breaking bonds.The cohesive energy of the lattice comes from adding energy of those bonds.Hence, the cohesive energy is equal to the vacancy formation energy; nevertheless, in metals, the vacancy formation energy can be typically one third of the cohesive energy.For instance, in the case of metallic bonding, the valence electrons are highly delocalized.Such a situation is for example most closely realized in the alkali metals, which readily give away their only weakly bound s electron in the valence shell.In these cases, is clear that a contribution arising from delocalized bonding cannot be described by a sum of pair potentials, because they are only adequate where there is a negligible distortion of the atomic electron density when the tom is added to the solid.If the binding arose only out of pairwise bonds with the nearest neighbor, the cohesive energy would be proportional to the coordination number that is not true in metals with highly delocalized electrons [26].For solving this issue, the potential of each atom is calculated as a sum over the contribution not only on the nearest neighbors but also on the all lattice atoms; The total energy E(r) of a crystal whose atoms are at room temperature can be obtained by summing equation ( 2) over the entire atoms of the crystal.This is done by choosing one atom in the lattice as an origin and calculating its interaction with other atoms in the crystal [27].The result obtained from such summation should be multiplied by N/2, where N is the atomic number of the material.The functional form of the total energy E(r) of a crystal in terms of the Morse potential, from this procedure is given as: The Lorentz-Berthelot mixing rule was used to estimate the interatomic Morse potential for materials A and B with parameters Dij, ro, and  for a mixed pair of atoms using the following formulas [28]  − = (    ) 1/2 (4) =   − (2/ − ) where DA-B, roA-B and A-B are the cohesion energy, equilibrium distance and lattice constant, respectively, for materials A and B (A and B are Cr, N or Si, depending on the case).The equations of motion were integrated using the Verlet algorithm with a dynamic time step of 1.92 fs.During the nanoindentation process, the rigid sphere is located at 7.2 Å above the surface of the sample.The indentation forces were obtained by summing the z-axis force of each individual atom that is affected in the sample.The dynamic time scale is calculated using the equation: Where  = √ 2  ⁄ is the dynamics time of the system, in which m, x and E are the physical quantity of mass, distance and energy, respectively.The variables µ,  and  are the scaling constants for the mass, distance and energy, respectively [23].Dynamic time was used to find an adequate t that can be used to perform the time integration.The dynamic time must be chosen carefully because if it is very large, the indenter will move toward the sample slowly and the sample deformation will not correspond to the deformation produced by the applied force; because the simulation time step is very short (1 fs), the random movement of the atoms cause the force to be unstable during the dump step [16].In different research studies using molecular dynamics [18], the use of femtosecond (fs) time scales is recommended.This femtosecond order of magnitude of time is obtained in a semiempirical way; furthermore, it is suitable for avoiding system collapse, thereby allowing the extraction of reliable data from the simulation.In this case, the timestep t=2.5 fs and 100 steps of dynamic time were used, having as a total simulation time of 250 fs, which proved to be adequate and reproducible, according to the values reported in the literature [18].
A reduced elastic modulus (E * ) is obtained from the following equation, which takes into account the combination of the elastic effect of the tip and the film [5]: where  = 1 is included for a spherical indenter, S is the slope of the initial part of the unloading curve, and Ac is the projected contact area at a maximum load.The elastic modulus of the sample (Es) is obtained from the expression [29]: where Ei is the elastic modulus of the indenter, and i and s are the Poisson's ratios of the sample and indenter, respectively.Ac is a function of indenter radius, R, and the contact depth, hc, as:   = ℎ  (11) where hc may be expressed as: where hmax is the maximum depth, Pmax the maximum load, and Smax is the slope of the unloading curve at the maximum load.Material hardness is defined as the local resistance to plastic deformation.Thus, the hardness, H, is determined from the maximum indentation load Pmax divided by the projected contact area, i.e.:

Results and Discussion
Figures 1 (a) and (b) show simulations of the load-unload curves during the nanoindentation process for Cr/CrN and (Cr/CrN)2, respectively.The process was performed 10 times while varying the indentation depth.The indentation depth was determined by measuring the penetration distance of the indenter into de sample as the force is applied [30].According to this figure, at the beginning of the indentation process, the local plasticity in the film is observes at very shallow indentation depth (i.e. the indentation size effect).When the indentation approaches to certain value, plastic deformation rises in the substrate and interfaces and the process enter to the next region.In this region, the local plasticity of the substrate appears at the indentation depth close to interface (i.e. the interaction of indentation size and substrate effect).[31].
Figures 1 (a) and (b) show simulations of the load-unload curves during the nanoindentation process for Cr/CrN and (Cr/CrN)2, respectively.The process was performed 10 times while varying the indentation depth.The indentation depth was determined by measuring the penetration distance of the indenter into de sample as the force is applied [30].
According to this figure, at the beginning of the indentation process, the local plasticity in the film is observes at very shallow indentation depth (i.e. the indentation size effect).When the indentation approaches to certain value, plastic deformation rises in the substrate and interfaces and the process enter to the next region.In this region, the local plasticity of the substrate appears at the indentation depth close to interface (i.e. the interaction of indentation size and substrate effect).[31].

S. Amaya-Roncancio et al.
These simulations were carried at a constant temperature of 300 K.According to the literature, as indentation depth increases, the temperature in the indentation region rises and the areas affected by the temperature can expand to subsurface.The bigger the indentation depth, the more effect there is on subsurface layer.But in general, the temperature is still below to the material phase transition.Moreover, as the indentation depth increases, atoms will be random-close-packed and the motion of atoms will be close to stagnation, decreasing the atomic vibrations normally produced by the temperature.Then, the free volume becomes smaller and atomic motion is obstructed [32].For this reason, is suitable to perform this kind of simulations as a constant temperature, as is done my many authors [4,33].From figure 1, the hardness for both systems was extracted using the Oliver-Parr and Cheng methods, as shown in Figures 2 (a) and (b), respectively.As the indentation depth continues to increase, the load curve starts to go up until it reaches a maximum depth.After reaching the maximum depth, the indenter begins to unload and returns to its original position.When the nanoindenter draws back, the plastic-deformed region undergoes partial elastic recovery, indicating the irreversibility of the plastic behavior during nanoindentation that is manifested in the load-displacement curves as a hysteretic response [34].According to this figure, the hardness increases as a function of the nanoindentation depth.In the case of Figure 2(a) (OP method), the hardness reaches values of 18 GPa and 20 GPa for Cr/CrN and (Cr/CrN)2, respectively.These values are similar to those reported by experimental studies.For example, Kot el al. [35] presented a study of multilayers of alternate Cr/CrN coatings deposited on AISI 301 steel using a pulsed Nd:YAG laser, which used micro indentation to measure both the hardness and the Young's modulus as a function of the bilayer period.The hardness reported for the (Cr/CrN)2 of approximately 17 GPa is not far from the value obtained in our work (20 GPa).Romero et al. [36] reported a study of polycrystalline Cr/CrN multilayers deposited onto silicon substrates using r.f.magnetron sputtering.The bilayer period values varied from the nanometric range (2-20 nm) to higher values.Nanoindentation measurements indicated multilayer hardness values ranging from 15.7 to 28 GPa.The same authors [37] presented other work of Cr/CrN multilayers between 15 and 23 GPa, depending on the layers thickness.Nevertheless, in the literature, there are few works presenting simulations of Cr/CrN, although similar works can be found for other materials.For example, Te-Hua Fang et al. [18] showed the nanoindentation response of multilayered thin films using molecular dynamics simulations.According to this work, when the indentation depth increased, the maximum load also increased.This behavior was also observed in our work, as is shown in Figures 1 and 2. Ping-Peng et al. [16] presented a three-dimensional molecular dynamic (MD) simulation of the nanoindentation of aluminum thin film on a silicon substrate.The results indicated that the hardness of the film increased with the loading rate.According to the results presented in Figure 2(b), the hardness obtained using the Cheng method [38] is higher than those calculated using the OP method, presenting values of 28 GPa and 34 GPa for Cr/CrN and (Cr/CrN)2, respectively.The values obtained from the correlation indicate the existence of scale effects in the system of multilayers.According to Cheng [39] and Bolshakov [40], the error can be attributed to the stress in the interlayer caused by the variations in the lattice parameter and the deformations and dislocations generated in the material during stabilization of its structure.
It is important to analyze curves of figure 2 compared with experiments.In simulations presented here, the samples can be considered as single crystals, because they have only one orientation.According to Liu et.al. [41], at the atomic depth indentation the hardness of the single crystal is higher than the hardness values measures from nanocrystallines samples, that is in opposition to the classical Hall-Petch effect observed at bulk materials.Then, systems with the layer thicknesses in the nanoscale for multilayers fallow the inverse Hall-Petch effect that establishes an inverse relationship between hardness and thickness.
Instance, Jang et al [42] studied the indentation process with and without gran boundaries.They observed lower hardness in the case of grain boundaries which is attributed to the interaction between dislocations and grain boundaries.As the indentation depth increases, in real materials more defects are found restricting the increase of the hardness.In the case of the Cr/CrN and (Cr/CrN)2 bilayers, a particular phenomenon in the interface occurs in which a structural change is produced.In the case of the interlayers, epitaxial stress plays an important role.The mechanical properties depend strongly on the type of interlayer, the crystalline structure and the lattice parameter.According to Sergey [43], the interlayers of multilayered materials can be classified in three categories: coherent, semi-coherent and incoherent.In the first case, the materials involved have similar crystalline structures and lattice parameters; in the second case, semi-coherent interfaces exist, and the crystalline structure is the same but the lattice parameter is notably different; finally, in the case of an incoherent interface, both, the crystallographic structure and the lattice parameter are different.Incoherent interfaces exhibit no continuous slip plane and thus are strong barriers to slip transmission.For example, Cr/CrN belongs to the last group because Cr has a BCC structure with a lattice parameter of 2.89 Å and CrN has a FCC structure with a lattice parameter of 4.18 Å.At very low scales, the interlayers can be considered to be grain boundaries in the internal structure, as presented in Figure 3 [44].It was found that several interfaces cause a greater hardness, especially for the case of incoherent interfaces, where compressive and tensile stresses in the interlayer are generated.The sliding process in the material is abruptly interrupted by the change in the crystalline structure, as shown in Figure 4.
In the case of the Cr/CrN bilayers, the incoherent-type interlayer generates a deformation around the zone where the crystalline structure changes, as shown in Figure 3.Such incoherence produces tensions in the atoms around the interlayer [45].This behavior causes the force exerted by the indenter to be rapidly dissipated by the material in different directions to the load application, as shown in Figure 4.
To reach an important deformation in the sample, a greater work of the indenter was required, thus causing the load-unload curve to exhibit greater values for the case of multilayers.As a consequence, the hardness increased for the multilayer case.The deformation and sliding caused by the tensions at the interlayers of the Cr/CrN and (Cr/CrN)2 samples are identified in Figure 4.The dislocation generations are not clearly observed.This may be due to the characteristics observing during other molecular dynamics simulations [4,46] that also showed that the deformation of indentation is highly localized and the stored elastic energy is not sufficient to create this kind of defects in the atomic scale.For the case of Cr/CrN and (Cr/CrN)2 thin films, the following equation is used, which is based on the Oliver-Pharr theory for the nanoindentation process [43], according to eq. (10) In general, to obtain the elasticity modulus of the multilayers of different materials, the Poisson ratio of the upper layer (CrN in this case) is used.Note that while the Poisson ratio of Cr is 0.2 [46], the value for CrN is 0.24 [47], indicating a very small difference in determining the Poisson ratio of the multilayer system.
For the particular case presented in this work, values of the elasticity modulus of the layers were obtained.In the case of Cr/CrN, the elasticity modulus presented values not greater than 217.86 GPa, and in the case of (Cr/CrN)2, a value of 258.9 GPa was observed.In both cases, the elasticity modulus is greater than that obtained for CrN, which, in this work, is 186.89GPa.
Having taken into account that the elasticity modulus is directly related to the deformation of the material, a material with a high value of the Young's modulus will experience less deformation (strain) for a given force (stress) than one with a lower value of the Young's modulus.In this case, as the number of bilayers was increased, the elasticity modulus also increased.This increase can be attributed to the change in the crystalline structure of the material.The changes in the directions of deformations caused by the indenter penetration cause the lower layers to undergo elastic deformations that help the material to recover its initial shape, contrary to the effect that can occur in monolayer thin films, which can suffer inelastic deformations.Because the elasticity modulus is the opposition to the deformation, the samples of Cr/CrN and (Cr/CrN)2 exhibit higher elasticity modulus values [48].To analyze the effect of temperature on the multilayer system hardness, simulations including variation in the temperature from 300 K to 1000 K were performed while keeping the other parameters fixed.The range of temperature was lower that the melting points of materials present in the simulations (2180 K, 2043 K and 1687 K for Cr, CrN and Si respectively).The values of the hardness obtained from the simulations for the Cr/CrN and (Cr/CrN)2 systems are shown in Figure 5, which reveals a tendency for the hardness to decrease with increasing temperature. of the incoherent interlayers, the changes in the temperature cause an increase in the number of dislocations; this behavior occurs because the interlayer incoherence was favored, according to Figure 6, which causes instability in the crystalline structure, thus causing the hardness to depend on the changes in the positions of the atoms caused by the instantaneous kinetic energy of each particles that is provided by the temperature.The result is a decrease in the hardness that depends on the vacancies and structural dislocations caused by the temperature as well as the deformations produced for stabilizing the tensions in the interlayer.In the literature, no results of an annealing process of Cr/CrN multilayers were found; nevertheless, a work carries out by Barshilia et al [54] shows the hardness evolution depending on the annealing temperature for TiN/CrN multilayers.The decrease of hardness as a function of the temperature was observed.This behavior was attributed to interdiffusion at the interfaces.An interesting point is the similitude with result presented in our work; the hardness decrease from about 36 to 28 GPa when temperature variyes between RT and 970 K.
As a final remark, it is necessary to clarify that most of results have to be considered as a qualitative behavior of the Cr/CrN multilayers.For obtaining more realistic results, it would be necessary to include other elements; for instance, to use the embedded atom methods, that has into account the effect of the electron density of the atoms included in the system.This is interesting especially when metals are present.Other important factor that affects the results is the presence of imperfections as vacancies, dislocations, grain boundaries among others.However, results obtaining in this work can, in a qualitative way, reproduce the load-unload curves, the variation of the hardness depending on the indentation depth and the annealing temperature.It makes these results interesting, considering that no reports for simulations this specific system is found in the literature.

Conclusions
A computational model based on the molecular dynamics method was developed to simulate the nanoindentation process of solid materials; this model was used to reproduce the mechanical behavior of Cr/CrN and (Cr/CrN)2 thin-film systems.The hardness was determined using both the Oliver Parr and Cheng correction methods.The Cheng correction resulted in higher hardness values because it avoids the influence of the scale effect.Regarding the elastic modulus, Cr/CrN and (Cr/CrN)2 exhibited simulated values of 217.86 GPa and 258.9 GPa, respectively.This work also presented a study of the influence of temperature on the systems; according to the results, the hardness decreases with temperature because of the instabilities in the crystalline structure, thereby causing the hardness to depend on the changes in the positions of the atoms caused by the instantaneous kinetic energy of each particle that is provided by the temperature.

Figure 2 .
Figure 2. Hardness versus the nanoindentation depth for Cr/CrN and (Cr/CrN)2 obtained using the (a) OP theory and (b) Cheng theory.Source: The authors.

Figure 3 .
Figure 3. Deformations generated by the tensions in the interface before the indentation for (a) Cr/CrN and (b) (Cr/CrN)2.Source: The authors.

Figure 4 .
Figure 4. Deformations generated by the tensions in the interface after the indentation for (a) Cr/CrN and (b) (Cr/CrN)2.Source: The authors.

Figure 5 .
Figure 5. Hardness as a function of the temperature for (a) Cr/CrN and (b) (Cr/CrN)2 (lines are for the eye guide).Source: The authors.

Table 1 .
Morse potential parameters for Cr, N and Si and their interactions.