II. MD SIMULATIONS

A. Simulation procedure

To simulate tension and nanoindentation with a 4 nm spherical indenter applying MD, the (111) slab models were constructed for the Al and Pd monolayers, and a four-layer model was constructed for the Al/Pd multilayer [Figs. 2(a)–2(c)]. The lengths of the slabs in the x- and y- directions were determined to fit the incommensurate nature of the Al and Pd lattice constants and to minimize the strains imposed on the interface due to lattice constant mismatch. Specifically, the dimensions are 240, 240, and 160 Å for the Al slab, and 233.4, 233.4, and 155.6 Å for the Pd slab in the x-, y-, and z-directions, respectively. The dimensions of the Al/Pd four-layer models are 280.7, 238.3, and 137.52 Å in the x-, y-, and z-directions, respectively. In the Al/Pd four-layer model, the thickness of the Al layer is 35.07 Å and that of the Pd layer is 33.69 Å, which results in a total height of 137.52 Å for the four-layer model. In the case of tensile loading, periodic boundary conditions were applied on the sides of the slabs in the direction of tension and sink conditions were applied on the other sides; therefore bulk material properties were obtained. This is different from extension of a single layer where beside the direction of tension, periodic boundary condition is also customary applied in the direction perpendicular to the tension in the plane. In the case of indentation loading, the boundary conditions are: periodic conditions in the plane (parallel to the surface) and sink condition out of plane (perpendicular to the surface).

FIG. 2.

Initial atomistic configurations of the (a) Al, (b) Pd (111) slabs, and (c) Al/Pd four-layer model used in the nanoindentation simulations. The red and gray circles represent the Al and Pd atoms, respectively.

In the MD simulations, the open source code LAMMPS was used. ^{Reference Plimpton13}

The accuracy of MD simulations is greatly dependent on the interatomic potential models. ^{Reference Shen and Chen14} In this study, a newly developed long-range empirical potential (LREP) model with a cut-off distance larger than at least the fourth neighbor distance was adopted, to describe Al–Al, Pd–Pd, and Al–Pd interatomic interactions. The LREP is an extension of the Finnis–Sinclare potentials, ^{Reference Finnis and Sinclare15} which has the same form as the embedded-atom-method ^{Reference Daw and Baskes16} except that dissimilar elements can contribute differently to the total electron density at an atomic site. This method has been successfully applied to body-centered cubic (bcc), face-centered cubic (fcc), and hexagonal close-packed (hcp) metal systems, as well as to binary and ternary fcc metal systems. ^{Reference Dai, Kong and Li17–Reference Li, Dai, Liang, Tai, Kong and Liu20} More details about the LREP models for Al–Al, Pd–Pd, and Al–Pd interatomic interactions can be found elsewhere. ^{Reference Kong and Shen21}

In the MD simulation of nanoindentation, the energy of the system was minimized using the conjugate gradient (CG) algorithm for the structure relaxation. The convergence tolerance for the energy minimization was set to 10⁻⁹, meaning that energy minimization stops when the energy change between two successive iterations divided by the energy magnitude is less than 10⁻⁹. After minimization, a sphere indenter of a radius of 40 Å with a contact stiffness of 246 eV/Ų was indented into the specimen using a displacement controlled loading. A smaller contact stiffness of 24.6 eV/Ų was also tested and the obtained force–displacement curve was found similar to that with the larger force constant. For the purpose of simulation efficiency, the contact stiffness of 246 eV/Ų was used in the following nanoindentation simulations.

B. Simulation results

First, the stress–strain curves (σ–ε) of Al and Pd, as shown in Figs. 3(a) and 3(b), were obtained by stretching Al and Pd slabs in two directions in the MD simulations. Then, indentation was performed on the Al and Pd slabs and the Al/Pd four-layer model. It can be seen in Fig. 4 that during the entire loading process, the force–displacement curve of the Al/Pd four-layer model is located above that of the Al slab but below that of the Pd slab. It is also found that in the unloading process, the depths where the applied force becomes zero are 9.54, 10.34, and 11.09 Å for the Al slab, Pd slab, and Al/Pd four-layer model, respectively, implying that plastic deformation is much profound in the Al/Pd four-layer model so that the least elastic recovery can be obtained. In indentation, the hardness of a specimen can be estimated from the maximum applied force divided by the residual indent area after unloading. The simulated hardness values are 13.06, 25.07, and 16.46 GPa for the Al slab, the Pd slab, and the Al/Pd four-layer model, respectively. The simulations suggest that the hardness of Al/Pd multilayer is between those of Al and Pd monolayer, which is not consistent with the experiments with the Berkovich indenter, where the measured hardness of Al/Pd multilayer is higher than those of Al and Pd monolayers [Fig. 1(a)]. ^{Reference Dayal, Savvides and Hoffman4} Also, the simulated hardness of Al/Pd multilayer and Al and Pd monolayers is much higher than the experimentally measured hardness values. Explanation for this discrepancy can be drawn from factors not considered in the present MD simulations—for example, solid phase transition processes at the interface between the layers during film deposition and different defects in the measured and simulated samples. Also, recent first-principles calculations of the Al/Pd interface energy suggest that the resulting intermetallics formed in the interface due to large negative enthalpy of mixing between Al and Pd can make huge positive contribution to the hardness enhancement of the Al/Pd thin film. ^{Reference Kong and Shen22} It was also found that the interface between Al and Pd is transparent to dislocations, thus the image forces may also play some role.

FIG. 3.

Tensile stress–strain curve σ–ε of: (a) Al layer and (b) Pd layer.

FIG. 4.

Reaction force P on the indenter versus the penetration depth h of the indenter.

In the present simulation, due to the computational capability, defects such as vacancies, dislocations or grain boundaries were not introduced in the initial configuration. The atoms in the initial Al and Pd slabs and Al/Pd four-layer model were perfectly in their fcc crystalline position. After minimization, due to the mismatch between the lattice constants of Al and Pd, some dislocations will raise in the interface of the Al/Pd four-layer model. The different amount of the defects which leads to the discrepancy between the calculated and measured hardness values can be clearly observed from the atomistic configurations in nanoindentation. Snapshots of the atomistic configurations corresponding to similar penetration depths of 17.14, 17.15, and 17.09 Å for the Al and Pd slabs and the Al/Pd four-layer model are shown in Figs. 5(a)–5(c). The atoms are colored according to their centro-symmetry parameter values and only atoms away from their perfect fcc positions are shown. The dislocation loops can be clearly seen in the figures. A noticeable phenomenon in the Al/Pd four-layer model is that the amount of dislocations formed during the nanoindentation is much larger than that in the Al and Pd slabs. Moreover, the generated dislocations are networked in the Al/Pd four-layer model. It can be found in the Al/Pd four-layer model that all dislocations nucleated under nanoindentation are confined in the first Pd and second Al layers, and there are no new dislocations formed in the third Pd and fourth Al layers except for the initial dislocations in the interface due to lattice constant mismatch. It is known that hardness is conventionally regarded as the resistance to plastic deformation and its value is limited by the smallest resistance to plastic deformation. Since there are a large amount of defects in the measured samples, the calculated hardness from the present nanoindentation simulations cannot be directly compared with the measured one. Instead, the calculated values give the upper bound to the hardness values.

FIG. 5.

Snapshots of the atomistic configurations of the (a) Al, (b) Pd (111) slabs, and (c) Al/Pd four-layer model under nanoindentation. The penetration depths are 17.14, 17.15, and 17.09 Å for the Al and Pd slabs and Al/Pd four-layer model, respectively. Atoms are colored according to their centro-symmetry parameter values and are not shown if they are in perfect fcc lattices. The dislocation loops are clearly seen.

III. MATERIAL PROPERTIES OF Al AND Pd

In the present study, the Young's modulus (E) and yield stress (σ_y) of bulk Al and Pd were derived from the MD simulations of tension [Figs. 6(a) and 6(b)]. The tensile σ–ε curves were fitted with smooth curves and the Young's moduli and yield strengths were determined in two directions, that is in the plane of the film and perpendicular to it. It appears that both Al and Pd are stiffer and stronger in the transverse direction than in the plane of the films. The properties of Al are: in the plane of the film E = 63.6 GPa and σ_y = 2.3 GPa, and perpendicular to it E = 71.9 GPa and σ_y = 2.6 GPa. Pd appears stronger and stiffer: in the plane of the film E = 160.3 GPa and σ_y = 5.9 GPa, and perpendicular to it E = 184.4 GPa and σ_y = 8.6 GPa. A material formulation capable of accounting for directional yield strength is the Hill's plasticity material model, according to which one direction is assumed as the reference direction and the properties in the other direction are scaled by a constant. The reference direction here was chosen to be the $\left[ {\bar 101} \right]$ direction and the stress values in the [111] direction were scaled 1.12 and 1.45 times for Al and Pd, respectively. For both materials, the Poisson's ratio was chosen to be 0.3, which is a common value for metals. It can be seen in Figs. 6(a) and 6(b) that after the peak point, the stresses in the MD simulations decrease with the decrease in the strain while they slightly increase in the FEA simulations. This is so because the material model in the FEA simulations cannot accept softening, for which reason a small tangent modulus was introduced.

FIG. 6.

Tensile stress–strain curves σ–ε used in the FEA simulations: (a) Al and (b) Pd.

IV. VALIDATION OF THE MATERIAL FORMULATION USING FEA

Validation of the material properties of Al and Pd and the material formulation of Hill's plasticity was performed by comparing the results of indentation on Al and Pd monolayers and the Al/Pd four-layer film carried out with MD, as detailed above and illustrated in Fig. 4, and FEA software. The FEA model was developed using the general purpose software ANSYS: the size of the Al monolayer was taken to be 16 × 24 nm and the size of the Pd monolayer was taken to be 15.56 × 23.34 nm; the spherical indenter was of 4 nm radius; the multilayered film consisted of two Al layers each of them 3.4 nm thick and two Pd layers each of them 3.5 nm thick, the length of the multilayered film was 23.34 nm.

The assumption that Al and Pd are transversely isotropic in the film plane simplified the simulations and reduced the 3D problem of indentation to an axisymmetric problem. Thus, only half of the film cross-section was modeled and the indenter was simulated as one line and hence assumed to be rigid. The cross-section was meshed with a mapped mesh of 2D quadrilateral elements—the areas of Al and Pd were meshed with 2D four-node linear structural solid elements (PLANE 182) having axisymmetric stress behavior, and the indenter was modeled with 2D rigid target elements (TARGET169). The contact between the indenter and the solid was modeled with 2D contact elements (CONTA171). The friction coefficient at the contact was assumed negligible, at 0.04. Fixed translational boundary conditions were imposed the bottom and on the sides of the films. In the multilayered film, the interface between the Al and Pd nanolayers was assumed perfect. Figures 7(a) and 7(b) show the FEA models of the monolayered and multilayered films. 4842 elements were used to simulate the Al and Pd monolayers [Fig. 7(a)] with the smallest element being 0.69 × 0.46 (Å) and 5242 elements were used to simulate the four-layered Al/Pd film [Fig. 7(b)] with the smallest element being 0.67 × 0.098 (Å).

FIG. 7.

The FEA model in indentation on MD specimens of: (a) Al (Pd) monolayers and (b) Al/Pd multilayers.

The indentation loading was imposed in terms of a vertical displacement applied on the indenter downward. On unloading, the indenter was moved upward until the contact was lost, which corresponds to a contact force of zero. The indentation force was determined as the reaction force on the rigid indenter. The maximum indentation depth imposed on the specimens was 2 nm in the Pd monolayer and Al/Pd multilayer, and 1.9 nm in the Al monolayer. A large number of loading steps were specified, in both loading and unloading, to obtain smooth P–h curves.

The P–h curves from the FEA simulations were compared with those from the MD simulations and were found to closely fit, which validates the material properties and the material formulation [Figs. 8(a)–8(c)]. Noticeably, the MD simulations of indentation of the multilayer predict softening at deep indentations while the FEA hardening. A possible explanation is the material formulation in the FEA model, which cannot resemble softening, as mentioned above in Sec. III.

FIG. 8.

P–h curves in indentation simulated using MD and FEA on: (a) Al monolayer, (b) Pd monolayer, and (c) Al/Pd multilayer.

V. FEA ON Al AND Pd MONOLAYERS, AND ON Al/Pd MULTILAYERS

With the validated material properties of Al and Pd, indentation could be simulated on Al and Pd monolayers and Al/Pd multilayers of a geometry used for experimental nanoindentation. Thus, the effect of indenter shape and indentation depth on the hardness values can be ascertained and compared with experiments. A large number of FEA simulations were carried out with a spherical indenter of 100 nm radius—the estimated tip radius of a diamond Berkovich indenter. The thickness of the Al and Pd monolayers was chosen as 500 nm, sufficiently large to avoid any boundary effects. The layer thickness in the multilayered films was varied from 2.25 to 43 nm. The maximum indentation depth (h _max) was varied from 3 to 54 nm, to cater for a large range of shallow, intermediate, and deep indentations. It is a practical rule of thumb that h _max should be less than 10% of the total film thickness to minimize the substrate effect on the measured mechanical properties. The present study imposed a more stringent rule—the stresses at the boundaries were outputted and when any boundary effect was noticed, the film thickness was increased. The width of the films was also checked against the stresses at the boundary and increased when required. An additional requirement for the multilayered films is that h _max should be greater than the thickness of one layer, to account for the multilayered nature of the film. In this study, h _max was set to at least 1.25 times the layer thickness (Fig. 9). For all those reasons, and to decrease the computer time required in the simulations, the size of the FEA models varied with h _max. The FEA models were simplified to axisymmetric and the same element types were used as in the FEA study in Sec. IV. The convergence of the FEA results was checked by carrying out simulations on models with different element sizes.

FIG. 9.

Schematic of a multilayered film indented with a spherical indenter.

The present study sought to answer the following questions “how is indentation hardness influenced by: (i) layer order, (ii) layer thickness, (iii) indentation depth, and (iv) the way the contact radius is measured?” Although most of the FEA simulations were carried out for multilayers having Pd as the top layer and Al as the bottom layer, to answer question (i), some simulations were performed with the inverse order—Al was placed at the top and Pd at the bottom. To answer question (ii), multilayered films of various layer thickness were studied, i.e., 2.25, 3.75, 5.5, 7.5, 10, 14, 15, 20, 25, 29, 36, 39, and 43 nm. Question (iii) above was answered by indenting the films to various depths, spanning from 3 to 54 nm. To reveal if hardness at a particular depth was layer-thickness dependent, some simulations were conducted for two or three different film structures, at the same maximum indentation depth. For example, the h _max of 18.75 nm was achieved in two ways—by indenting a multilayer of a layer thickness d = 7.5 nm to h _max = 2.5 × d and by indenting a multilayer of d = 10 nm thick layers to a depth of h _max =1.875 × d. An additional consideration was the number of the layers—a sufficient number of layers were taken to avoid any boundary effect: 12, 16, 24, 28, 32, and 48 layers. Apart from Al/Pd multilayered films, Al and Pd monolayers were also indented to the same h _max. Question (iv) above was answered by determining the contact radius in three ways, that is, (i) by Field and Swain's formula, (ii) at the point of inflection of the imprint at the maximum indentation depth, and (iii) at the point of inflection of the imprint after unloading.

It needs to be noted that, although care was taken to ensure computational accuracy of the results, the main objective of the study was not the accuracy of the hardness determination but rather qualitative outcomes. Also, the material properties of Pd and Al were extracted from MD simulations and may not necessarily be the material properties of experimentally deposited films. It is well recognized that the elastic modulus is essentially scale-independent; however, the material strength at the nanoscale is significantly higher than the one at the macroscale. Nevertheless, using MD to provide information which otherwise is unavailable at the macroscale proves a valuable tool in obtaining a deeper insight into the material behavior of thin films.