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

We perform atomistic simulations of nanoindentation on Lennard–Jones 2D hexagonal crystals. In this work, we find a new spatially extended buckling-like mode of instability, which competes with the previously known instability governed by dislocation-dipole nucleation. The geometrical parameters governing these instabilities are the lattice constant, a, the radius of curvature of the indenter, R, and the thickness of the indenter layer, Ly. Whereas dislocation nucleation is a saddle-node bifurcation governed by R/a, the buckling-like instability is a pitchfork bifurcation (like classical Euler buckling) governed by R=Ly. The two modes of instability exhibit strikingly different behaviors after the onset of instability. The dislocation nucleation mode results in a stable final configuration containing a surface step and a stable dislocation at some depth beneath the surface, while the buckling modes are always followed immediately by subsequent nucleation of many dislocation dipoles. We show that this subsequent dislocation nucleation is also observed immediately after buckling in free standing rods, but only for rods which are of sufficiently wide aspect ratio, while thinner rods exhibit stable buckling followed only later by dislocation nucleation in the buckled state. Finally, we study the utility of several recently proposed local and quasi-local stability criteria in detecting the buckling mode. We find that the so-called K criterion, based on the stability of a representative homogeneously deformed lattice, is surprisingly useful in detecting the transition from dislocation-type instability to buckling-type instability. [DOI: 10.1115/1.4034564]


Introduction
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][2][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 h110i 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][13][14][15].
In the past decade, many researchers have attempted to formulate a criterion to predict dislocation nucleation [12][13][14][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][13][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][18][19]. Specifically, dislocation-dipole nucleation involves an embryo consisting of adjacent planes of atoms slipping against each other to nucleate a dislocation dipole [20][21][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 1 the nearest-neighbor direction coincident with the loading axis. Some questions to consider in this scenario are: (1) What is the mechanism of this instability process? (2) When does the transition occur? (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 K field is quasilocal at the dislocation-dipole loop in the traditional case of dislocation-dipole nucleation. In the case of buckling-type mode, K field is highly delocalized in the buckled region in the bulk of the crystal. In both cases, dislocation-dipole nucleation and buckling, K becomes negative before nucleation. An important result presented here is that the K 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 dislocationdipole 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 longwavelength buckling eigenmode. In Sec. 3, we study the bifurcation mechanism associated with this new instability process. In Sec, 4, we present the K criterion and FDM linear stability based analysis of the delocalized buckling eigenmode. Section 5 presents some concluding remarks.

Transition From Localized (Dislocation Nucleation) to Long-Wavelength (Buckling) Instability
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 1 and O 2 are shown in Fig. 1. O 1 has the indenter motion direction perpendicular to the nearest-neighbor direction. O 2 is the highest surface energy orientation having nearest-neighbor direction parallel to the indenter motion direction.
We configure 2D hexagonal crystals as shown in Fig. 1 with periodic boundaries on the sides, rigid base at the bottom, and a circular indenter on top of the crystal. A large enough L x was chosen to ensure that the results are independent of its value. In the rest of the paper, all the lengths such as Lx, L, R, and C are measured in units of the lattice constant, a. Energy and forces are measured in LJ units. We use a stiff, featureless, harmonic, repulsive, cylindrical indenter for all our simulations as in Refs. [14,15]. The interaction potential used to model the interaction between the indenter and the atoms is of the form 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 1 and some systems in O 2 , the load drop is accompanied by nucleation of a dislocation dipole. However, O 2 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 ¼ L y ; 2L y ; and 3L y , are shown in Fig. 2 The total potential energy, Uðx ia ; DÞ, is a function of atomic positions, x ia , and indenter depth, D. The first derivative of energy with respect to the particle position gives the force, F ia , on each particle as Latin characters are used to index particle number, and Greek characters to index Cartesian components. Then, we use the Lanczos algorithm, as implemented in the SciPy toolkit [24], to compute the lowest four eigenvalues of the Hessian matrix for the relaxed configuration at each indenter step. The Hessian matrix, H iajb , is the second derivative of total potential energy as The forces induced by an infinitesimal external indenter motion, N ia , must be balanced by the internal atomic rearrangements as shown in Ref. [21] Here, _ x¼ : ðdx=dDÞ. We refer to the derivative of the particle positions with respect to indenter depth as velocities since D plays the role of time in quasi-static indentation. We may solve Eq. (4) at each indentation step to compute the atomic velocities. The analytical expression of H iajb can be simply derived for pair potentials such as LJ potential and Morse Potential using the below equation [25] M where t and c are the first and second derivatives of the bond energy, and n ija is the unit normal pointing from particle i to particle j. Then, H iajb ¼ ÀM iajb for off-diagonal terms and H iaib ¼ P j M iajb 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 2 indented by indenters having R ¼ L y ; 2L y ; and 3L y are shown in the left panel of Fig. 3. For a given scalar quantity, we compute X as in Ref. [26]. X 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, X is computed by linear interpolation of the critical mode vectors. The X field corresponding to the critical eigenmodes is shown in the right panel of Fig. 3. For R ¼ L y and R ¼ 2L y , 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 ¼ 3L y , 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 Ã =L y , 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 1 and O 2 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 2 and roughly 10% below the critical axial strain at R ¼ 2:75L y . O 1 , on the other hand, becomes unstable via dislocation-dipole nucleation even for very large indenter radii.

Bifurcation Mechanism
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 X field corresponding to the eigenmode is shown for R ¼ L y ; 2L y ; and 3L y . 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 ¼ L y ; 2L y ; and 3L y , O 2 . For R ¼ L y and 2L y , a single eigenmode descends through the spectrum according to dD 0:5 scaling. This square-root scaling corresponds to the saddle-node bifurcation as shown by Hasan and Maloney [21]. However, for R ¼ 3L y , the critical eigenvalue vanishes linearly with dD. 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 bucklingtype 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 L x =L y ¼ 0:2. The sequence of structure of a rod with L x =L y ¼ 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 L x =L y ¼ 0:33, dislocation nucleation occurs immediately after Transactions of the ASME buckling as shown in Fig. 7(b), and for L x =L y > 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.

Nucleation Criterion
Many researchers [12][13][14][15] have tried to formulate a criterion to predict homogeneous dislocation nucleation. Li and coworkers [12,13] developed a criterion known as the K criterion based on Hill's analysis of the stability of plane waves in a homogeneously deformed crystal. K is related to the dynamical matrix which gives the vibrational frequencies of phonons of a given wavevector. The dynamical matrix, D il , for the homogeneously deformed lattice at the ith particle can be computed as For a 2D hexagonal lattice, D il is a 2 Â 2 matrix indexed by l; . R ij is the displacement vector defined from particle, i, to the neighboring particle, j, in the homogeneously deformed crystal. H l contains the elements of the Hessian matrix for a homogeneously deformed crystal. k andk are, respectively, the magnitude and direction of the wavevector k. Then For particle i, K is the minimum eigenvalue of K il over all the values of k. k 2 K min is then the lowest squared phonon frequency for the phonon of given wavevector. A negative value of K min indicates an unstable phonon mode [20,22]. Miller and Rodney showed that K becomes negative before the actual instability. However, K 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 K cannot distinguish between dislocation-dipole eigenmode and phase-boundary nucleation. A following question is: How does K change as the critical eigenmode transforms in localized to delocalized eigenmode as the indenter radius increases. Figure 8 shows the K field for R ¼ L y ; 2L y ; and 3L y for L y ¼ 120. As observed earlier [20,22], K is minimum (and negative) at the localized embryo for R ¼ L y and, therefore, predicts the location of dislocation nucleation. However, for R ¼ 2L y , K is negative at the embryo and at the center line or the indentation axis. For R ¼ 3L y , K is negative in a huge region surrounding the diffused embryo. However, K becomes negative before the actual instability both for dislocation-dipole and long-wavelength buckling instability. K field changes from a local field for dislocationdipole instability to a delocalized field for buckling-type instability, and therefore, it can distinguish between dislocation-dipole instability and buckling instability.
Recently, it was shown that the linear stability analysis of the evolution equation of the dislocation density tensor in field dislocation mechanics can predict both the location and instant of dislocation nucleation. The result holds for many interatomic potentials and crystallographic orientations in 2D and 3D nanoindentation simulations. The evolution equation of the dislocation density field is given by where the time derivative corresponds to the spatial representation of the a field. Here, a is the dislocation density tensor, v is the material velocity vector, V is the dislocation velocity vector relative to the material, and s is a dislocation nucleation rate tensor. For any physically reasonable constitutive equation for the dislocation velocity V, it may be assumed that V ¼ 0 if a ¼ 0. As shown in Refs. [16,22], the governing equation for linear stability analysis becomes In terms of components with respect to a rectangular Cartesian coordinate system, Eq. (9) can be written as where e jrs is a component of the third-order alternating tensor. The notation da ¼ a is used. Equation (10) can be rewritten as where a is the perturbation field of the dislocation density tensor. Then, the components of a ij are considered as a 9 Â 1 vector A and ðv j;m d ir À v k;k d ir d jm Þ as a 9 Â 9 array denoted by N to write Eq. (11) as The maximum of the real parts of the eigenvalues, eig realÀpart , of N at any point is given by the value of the field g at that point For dislocation nucleation, g increases by orders of magnitude just before nucleation. In Fig. 9, g field is shown for the dislocationdipole and buckling eigenmode. g is two orders of magnitude higher for the dislocation-dipole mode as compared to the longwavelength buckling mode. The linear stability analysis can distinguish between dislocation and buckling instability.

Summary and Future Directions
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 sequence of rod configuration as it is compressed: (b) configuration just before buckling at "red cross" in (a), (c) configuration just after buckling, and (d) the buckled configuration emits dislocations from edges when compressed further at "black cross" in (a) as shown. Fig. 7 Potential energy, U, as a function of yy : (a) L x =L y 5 0.2, (b) L x =L y 5 0.33, and (c) L x =L y 5 0.5. As the crystal film thickness increases, the buckled configuration emits dislocations sooner. For large enough L x =L y as in (c), the buckling accompanies nucleation of dislocations.  7)) criterion for (a) R 5 L y , (b) R 5 2 L y , and (c) R 5 3 L y . The arrows in the figures represent the critical eigenmode. Fig. 9 Critical eigenmodes just before instability for buckling and dislocation-dipole nucleation as shown in (a) R 5 3L y and (b) R 5 L y . Corresponding g (given by Eq. (13)) calculated using linear stability of FDM, for LJ crystal, O 2 : (c) R 5 3L y and (d) R 5 L y . g is orders of magnitudes higher for dislocation-dipole nucleation as compared to the buckling.
The transition from localized to delocalized eigenmode is observed only for certain loading directions. Orientation, O 1 , crystal film fails by the long-wavelength mode when compressed homogeneously. However, O 1 crystal film shows dislocationdipole nucleation during nanoindentation even when it is compressed with a large indenter radius (R ¼ 15L y ), and the critical axial strain is very close to the critical axial strain for homogeneous loading. The highest surface energy orientation, O 2 , 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 K criterion was useful for predicting dislocation-dipole nucleation, it is also useful for predicting buckling. We find that the K field is localized at the embryonic loop for dislocation-dipole nucleation; however, for buckling-type instability, K field is delocalized in the long-wavelength region. For both dislocation-dipole nucleation and buckling, K 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 longwavelength mode spans entire crystal length. However, our preliminary results showed that the size of the buckled region in O 2 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.