Mechanical Instabilities in Perfect Crystals: From Dislocation Nucleation to Bucklinglike Modes

Often instabilities govern the critical behavior of materials during various mechanical processes. Therefore, it is important to identify the instability mechanisms and understand the mechanical conditions that lead to these instabilities. For the past several decades, dislocation nucleation has been shown to be the instability mechanism that governs the yield strength of crystal films, initially free of defects. Dislocation nucleation from free surfaces has been shown to determine the strength of micro- and nanopillars [1–3]. However, recently Rabkin et al. [4] and Marian and Knapp [5] observed buckling followed by nucleation of dislocations governing the strength of nanopillars. Buckling has also been observed during nanocompression experiments and simulations for nickel crystals [6,7] and silicon nanowires [8]. Marian and Knapp [5] showed transition from volumetric to surface dominated plasticity when the diameter of the wires is less than a critical value, d_c, and provided lower bounds for the yield stress of nanopillars during homogeneous compression. Rabkin et al. [4] observed buckling of these nanopillars only for certain interatomic potentials and associated buckling with FCC–HCP phase transformation. Buckling controls the strength of gold nanowires only when compressed along $〈110〉$ at low temperatures and high strain rates [9]. Their energy barrier calculations could predict the athermal strength fairly well, but only for the cases where buckling was not the critical strength governing phenomenon.

It follows from this short review that buckling is an important instability mechanism preceding nucleation of dislocations during homogeneous compression of gold nanopillars for certain loading directions, crystal geometries, and interatomic potentials. In this work, we focus on nanoindentation of thick crystal films. Presence of strain gradients and high stresses during nanoindentation creates conditions that strongly differ from the homogeneously compressed samples. During nanoindentation experiments of pristine gold films, dislocation nucleation in the bulk has been shown to be the critical instability process [10]. Homogeneous dislocation nucleation events have also been observed in 2D bubble raft experiments [11]. Moreover, dislocation nucleation has also been shown to govern the yield strength of thick pristine crystal films for many different interatomic potentials and loading directions using atomistic simulations [12–15].

In the past decade, many researchers have attempted to formulate a criterion to predict dislocation nucleation [12–15]. Ideally, the criterion should predict both the location and instant of nucleation given the crystal configuration and external conditions. The criteria thus formulated have been both local (e.g., based on local stress or strains) [12–14] and nonlocal (e.g., based on acoustic matrix of cluster of atoms) [15,16] in nature. Dislocation nucleation was shown to be an inherently mesoscale process [15]. The importance of nonlocality was further addressed in a nonlocal instability criterion called the Wallace criterion [17–19]. Specifically, dislocation-dipole nucleation involves an embryo consisting of adjacent planes of atoms slipping against each other to nucleate a dislocation dipole [20–22]. The embryo volume size converges as the number of atoms increases for both the bulk of FCC crystal and for surface nucleation [18].

In this work, we show a new critical eigenmode (or instability mechanism) having a long-wavelength structure similar to buckling in the bulk of the crystal, accompanied by nucleation of several dislocations. This newly observed eigenmode competes with the dislocation-dipole eigenmode as the indenter radius increases. The embryo structure in the buckling eigenmode is sufficiently larger than the embryo in dislocation-dipole eigenmode. The buckling eigenmode exhibits a nontrivial spatial extent and a dominant wavelength. This buckling eigenmode is only observed for the high surface energy crystallographic orientation that has the nearest-neighbor direction coincident with the loading axis. Some questions to consider in this scenario are:

1. (1)

   What is the mechanism of this instability process?

2. (2)

   When does the transition occur?

3. (3)

   How do the nucleation criteria proposed earlier work in this scenario? Can they predict the location and instant of buckling? Can these criteria distinguish between dislocation nucleation and buckling?

We perform eigenmode analysis to understand the bifurcation mechanism. The new instability process is shown to be governed by pitchfork bifurcation similar to classical Euler buckling. The critical axial strain during nanoindentation of crystal films has been compared to the critical axial strain for the instability of homogeneously compressed films. The transition from localized (dislocation-dipole) to delocalized (buckling) eigenmode is shown to occur approximately at 10% below the critical axial strain for homogeneous compression.

It has been shown earlier [15,20,22] that the Λ field is quasi-local at the dislocation-dipole loop in the traditional case of dislocation-dipole nucleation. In the case of buckling-type mode, Λ field is highly delocalized in the buckled region in the bulk of the crystal. In both cases, dislocation-dipole nucleation and buckling, Λ becomes negative before nucleation. An important result presented here is that the Λ criterion developed by VanVliet et al. [13] is useful in detecting the transition from dislocation-type instability to buckling-type instability. Another criterion, based on the linear stability analysis of field dislocation mechanics (FDM), has been shown to predict dislocation-dipole nucleation [16,22] previously. The dislocation nucleation criterion in FDM is not meant to predict smooth bifurcation modes as shown in Refs. [16,22], and in this work, the FDM-based criterion shows no signature for the buckling-type mode. In this sense, the field dislocation mechanics based criterion can also distinguish dislocation-dipole eigenmode from the buckling eigenmode.

In Sec. 2, we describe the simulation formalism for different crystal orientations. We discuss the key parameters involved in the transition from dislocation-dipole eigenmode to the new long-wavelength buckling eigenmode. In Sec. 3, we study the bifurcation mechanism associated with this new instability process. In Sec, 4, we present the Λ criterion and FDM linear stability based analysis of the delocalized buckling eigenmode. Section 5 presents some concluding remarks.

We press a rigid, frictionless indenter on two different orientations of a 2D hexagonal Lennard–Jones (LJ) crystal. The athermal, quasi-static nanoindentation simulations are performed via energy minimization dynamics. The LAMMPS molecular dynamics framework [23] is used to perform nonlinear energy minimization using the conjugate-gradient algorithm with the Polak–Ribiere modification. Crystallographic orientations O₁ and O₂ are shown in Fig. 1. O₁ has the indenter motion direction perpendicular to the nearest-neighbor direction. O₂ is the highest surface energy orientation having nearest-neighbor direction parallel to the indenter motion direction.

Here, R is the radius of the indenter, and r is the distance of the particle from the indenter center.

The indenter load, F, and the total potential energy, U, increase as the indenter moves toward the crystal until the crystal becomes unstable. For all the systems in O₁ and some systems in O₂, the load drop is accompanied by nucleation of a dislocation dipole. However, O₂ crystal shows significantly larger drop in F followed by nucleation of many dislocations simultaneously when indented by large indenters. F versus D curves for 2D hex lattices of thickness, L_y, pressed by indenters of radius, $R=Ly,2Ly, and 3Ly$⁠, are shown in Fig. 2. For $R=3Ly$⁠, there is an unsystematic increase in the F discontinuity.

where t and c are the first and second derivatives of the bond energy, and $nijα$ is the unit normal pointing from particle i to particle j. Then, $Hiαjβ=−Miαjβ$ for off-diagonal terms and $Hiαiβ=∑jMiαjβ$ for diagonal terms. However, for multibody potentials like EAM potential, the calculation of Hessian matrix is more involved.

We compute the lowest four eigenvalues of the Hessian matrix for the relaxed configurations at each indenter step. In all the cases, we observe that the system is driven to instability along a single eigenmode. The critical eigenmodes for crystal films in O₂ indented by indenters having $R=Ly,2Ly, and 3Ly$ are shown in the left panel of Fig. 3. For a given scalar quantity, we compute Ω as in Ref. [26]. Ω is the transverse gradient of the eigenmode along a particular crystal axis in the deformed configuration. The lattice is first triangulated, and on each triangle, Ω is computed by linear interpolation of the critical mode vectors. The Ω field corresponding to the critical eigenmodes is shown in the right panel of Fig. 3. For R = L_y and $R=2Ly$⁠, the critical eigenmodes are sharply localized along a slip plane. There are adjacent planes of atoms slipping with respect to each other to nucleate a dislocation dipole. For $R=3Ly$⁠, the critical eigenmode is delocalized having a long-wavelength structure similar to buckling. The spatial extent of this delocalized region could be governed by the system size and indenter radius.

The critical axial strain, $D*/Ly$⁠, at which an infinite crystal during homogeneous compression becomes unstable is analytically calculated using the dynamical matrix calculation as done in Refs. [20,22]. $D*$ is the indenter displacement at which the eigenvalue of the dynamical matrix becomes negative. Figure 4 shows the axial strain for various indenter radii during nanoindentation at which the systems become unstable. The critical axial strain for O₁ and O₂ is shown by the solid lines in Fig. 4. The circles in Fig. 4 represent the instability via dislocation nucleation, and crosses represent the instability via buckling. It is interesting that this transition is observed only for O₂ and roughly $10%$ below the critical axial strain at $R=2.75Ly$⁠. O₁, on the other hand, becomes unstable via dislocation-dipole nucleation even for very large indenter radii.

The lowest four eigenvalues of the Hessian matrix for the relaxed configurations are computed at each indenter step to identify the reaction coordinate. In all the cases, we observe that the system is driven to instability along a single eigenmode.

In the left panel of Fig. 5, the Ω field corresponding to the eigenmode is shown for $R=Ly,2Ly, and 3Ly$⁠. In the right panel of Fig. 5, we plot four lowest eigenvalues of the Hessian matrix as we approach the instability for three indenter radii, $R=Ly,2Ly, and 3Ly$⁠, O₂. For R = L_y and $2Ly$⁠, a single eigenmode descends through the spectrum according to $δD0.5$ scaling. This square-root scaling corresponds to the saddle-node bifurcation as shown by Hasan and Maloney [21]. However, for $R=3Ly$⁠, the critical eigenvalue vanishes linearly with δD. The linear scaling corresponds to the pitchfork bifurcation, also known for the classical Euler buckling. There are two degenerate eigenmodes descending through the spectrum.

The two mechanisms for instability result in strikingly different behaviors after the onset of instability. The dislocation-dipole nucleation is followed by a configuration with a surface step and a stable dislocation in the bulk of a crystal, whereas the buckling-type instability is immediately followed by nucleation of several dislocation dipoles. The nucleation of several dislocation dipoles in buckling-type instability results in the drop in F versus D curves as shown in Fig. 2. The nucleation of several dislocations immediately following buckling is also observed for thick free standing atomistic rods as discussed below.

Classical Euler buckling of a continuum rod is a continuous pitchfork bifurcation. The discontinuity in F versus D plots as shown in Fig. 2 is not expected for Euler buckling of a rod. In this work, for the buckling in the bulk of the crystal, buckling occurs in the absence of free surface and in the presence of strain gradients. The discontinuity in F versus D curves is also observed when a crystal film is compressed homogeneously. This implies that strain gradient most likely is not the primary reason for the discontinuity. We compare homogeneously compressing a continuum rod to an atomistic crystal film (or atomistic rods) to understand the discontinuity in F versus D plots. There can be two possible main differences: (a) A crystal film has an atomistic length scale whereas a continuum rod does not, and (b) a crystal film has a finite thickness as compared to the negligible thickness of a continuum rod.

Here, we perform an experiment to understand the discontinuity in F versus D plots for homogeneously compressed atomistic films. We homogeneously compress LJ 2D rods with varying aspect ratios. Potential energy, U, as a function of strain, $ϵyy$⁠, is shown in Fig. 6(a) for a rod with $Lx/Ly=0.2$⁠. The sequence of structure of a rod with $Lx/Ly=0.2$ during compression is shown in Figs. 6(b)–6(d). The rod $buckles$ at 12% $ϵyy$ as shown in Fig. 6(a), and the corresponding buckled configuration is shown in Fig. 6(c). When the rod buckles, U versus $ϵyy$ curve is continuous. At 14% strain, the buckled rod emits dislocations from the highest curvature (strain gradient) locations. The emission of dislocations is accompanied by the discontinuity in U versus D curve as shown in Fig. 6(a). The buckled configuration having dislocations at 14% strain is shown in Fig. 6(d).

The preceding discussion shows that a thin LJ 2D rod under homogeneous compression first becomes unstable due to buckling and then after further compression, buckling is followed by nucleation of dislocations. In Fig. 7, we show that the dislocation nucleation is immediately observed after buckling for thick free standing rods. The red cross in Fig. 7 corresponds to the instant of buckling, and black cross corresponds to the instant of dislocation emission. As the aspect ratio (or the film thickness) increases to $Lx/Ly=0.33$⁠, dislocation nucleation occurs immediately after buckling as shown in Fig. 7(b), and for $Lx/Ly>0.33$⁠, buckling is immediately followed by dislocation emission as shown in Fig. 7(c). Consequently, the crystal film having thickness greater than a critical value emits dislocations as soon as it buckles resulting in discontinuity in the F versus D curves.

For particle i, Λ is the minimum eigenvalue of $Λiμν$ over all the values of $k$⁠. $k2Λmin$ is then the lowest squared phonon frequency for the phonon of given wavevector. A negative value of Λ_min indicates an unstable phonon mode [20,22]. Miller and Rodney showed that Λ becomes negative before the actual instability. However, Λ_min corresponds to the location of dislocation embryo and, therefore, can be used for predicting the location of dislocation nucleation for many different interatomic potentials and orientations as shown by Garg and Maloney [20,22]. In the paper by Garg et al. [16,22], it was shown that Λ cannot distinguish between dislocation-dipole eigenmode and phase-boundary nucleation. A following question is: How does Λ change as the critical eigenmode transforms in localized to delocalized eigenmode as the indenter radius increases.

Figure 8 shows the Λ field for $R=Ly,2Ly, and 3Ly$ for L_y = 120. As observed earlier [20,22], Λ is minimum (and negative) at the localized embryo for R = L_y and, therefore, predicts the location of dislocation nucleation. However, for $R=2Ly$⁠, Λ is negative at the embryo and at the center line or the indentation axis. For $R=3Ly$⁠, Λ is negative in a huge region surrounding the diffused embryo. However, Λ becomes negative before the actual instability both for dislocation-dipole and long-wavelength buckling instability. Λ field changes from a local field for dislocation-dipole instability to a delocalized field for buckling-type instability, and therefore, it can distinguish between dislocation-dipole instability and buckling instability.

For dislocation nucleation, η increases by orders of magnitude just before nucleation. In Fig. 9, η field is shown for the dislocation-dipole and buckling eigenmode. η is two orders of magnitude higher for the dislocation-dipole mode as compared to the long-wavelength buckling mode. The linear stability analysis can distinguish between dislocation and buckling instability.

In this paper, we have shown a new instability mechanism for thick crystal films during nanoindentation. The commonly observed dislocation-dipole eigenmode competes with the long-wavelength buckling eigenmode as the indenter radius increases. Unusually, buckling has been shown to occur in the absence of free surfaces and in the presence of strain gradients. As the critical eigenmode transitions from the dislocation-dipole nucleation to buckling-type instability, the bifurcation mechanism changes from the saddle-node bifurcation to the pitchfork bifurcation. After the dislocation-dipole nucleation, one dislocation forms a surface step and the other one is a stable dislocation at some depth beneath the surface. The buckling-type instability is followed by nucleation of several dislocations. The subsequent nucleation of several dislocations is also shown to occur after buckling in free standing rods. When a thin rod is compressed, it buckles and then, during further compression the buckled rod emits dislocations from the highest strain gradient regions. For sufficiently wide aspect ratio rods, buckling is immediately followed by nucleation of several dislocations.

The transition from localized to delocalized eigenmode is observed only for certain loading directions. Orientation, O₁, crystal film fails by the long-wavelength mode when compressed homogeneously. However, O₁ crystal film shows dislocation-dipole nucleation during nanoindentation even when it is compressed with a large indenter radius (⁠$R=15Ly$⁠), and the critical axial strain is very close to the critical axial strain for homogeneous loading. The highest surface energy orientation, O₂, crystal film buckles during nanoindentation approximately 10% below the critical axial strain for homogeneous compression.

We study the effectiveness of several dislocation-dipole nucleation criteria previously proposed for detecting buckling. To the extent that the local Λ criterion was useful for predicting dislocation-dipole nucleation, it is also useful for predicting buckling. We find that the Λ field is localized at the embryonic loop for dislocation-dipole nucleation; however, for buckling-type instability, Λ field is delocalized in the long-wavelength region. For both dislocation-dipole nucleation and buckling, Λ becomes negative before nucleation, and therefore, it cannot be used to predict the instant of nucleation. The dislocation nucleation criterion in FDM cannot predict smooth bifurcation modes as shown in Refs. [20,22] and, therefore, cannot predict buckling. This specificity is its power in the sense that this analysis can be used to differentiate dislocation-dipole eigenmode from the buckling-type eigenmode.

The scaling of the wavelength of the buckling-type eigenmode with system size, L_y, needs to be understood. For homogeneous compression, the size of the buckled region or the long-wavelength mode spans entire crystal length. However, our preliminary results showed that the size of the buckled region in O₂ crystal film does not increase significantly when the indenter radius is doubled. Resolving this paradox may help in constructing appropriate framework to predict this instability.

Finally, the results shown here are for 2D atomistic simulations. This instability process still needs to be seen in 3D simulations and experiments at finite temperature. For dislocation-dipole nucleation, the energy barrier decreases with increase in temperature. Buckling is a geometrical instability, and the critical stress lacks a thermally activated component. Therefore, the energy barrier for buckling instability might not be temperature dependent. In that case, it would be difficult to observe buckling at higher temperatures.

This material is based upon the work supported by the National Science Foundation under Award No. CMMI-1100245.