## Abstract

We perform atomistic simulations of dislocation nucleation in two-dimensional (2D) and three-dimensional (3D) defect-free hexagonal crystals during nanoindentation with circular (2D) or spherical (3D) indenters. The incipient embryo structure in the critical eigenmode of the mesoregions is analyzed to study homogeneous dislocation nucleation. The critical eigenmode or dislocation embryo is found to be localized along a line (or plane in 3D) of atoms with a lateral extent, $\xi $, at some depth, $Y*$, below the surface. The lowest energy eigenmode for mesoregions of varying radius, r_{meso}, centered on the localized region of the critical eigenmode is computed. The energy of the lowest eigenmode, λ* _{meso}*, decays very rapidly with increasing r

_{meso}and λ

*≈ 0 for $rmeso\u2273\xi $. The analysis of a mesoscale region in the material can reveal the presence of incipient instability even for $rmeso\u2272\xi $ but gives reasonable estimate for the energy and spatial extent of the critical mode only for $rmeso\u2273\xi $. When the mesoregion is not centered at the localized region, we show that the mesoregion should contain a critical part of the embryo (and not only the center of embryo) to reveal instability. This scenario indicates that homogeneous dislocation nucleation is a quasilocal phenomenon. Also, the critical eigenmode for the mesoscale region reveals instability much sooner than the full system eigenmode. We use mesoscale analysis to verify the scaling laws shown previously by Garg and Maloney in 2D [2016, “Universal Scaling Laws for Homogeneous Dissociation Nucleation During Nano-Indentation,” J. Mech. Phys. Solids, 95, pp. 742–754.] for the size, ξ, and depth from the surface, Y*, of the dislocation embryo with respect to indenter radius, R, in full 3D simulations.*

_{meso}## 1 Introduction

In perfect crystals the onset of plasticity is indicated by the events of dislocation nucleation. Instabilities govern the critical behavior of materials during various mechanical processes. 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]. Applications such as chemical mechanical polishing (CMP) [4–6] have nanoscale interconnects in which dislocation nucleation has been shown as a stress relaxation mechanism [7,8]. In microscale crystals, dislocations are typically produced by pre-existing defects in the material. However, in recent years, tools such as nanoindentation and the atomic force microscopy (AFM) have allowed access to mechanical properties in nanoscale samples [9–20] where one can obtain practically defect-free crystals. In these scenarios, one is typically interested in how temperature, loading rate, or the sample or probe size affect the mechanical response.

Miller and Rodney (MR) showed that homogeneous dislocation nucleation (HDN) is a mesoscale process involving tens of atoms forming an *embryo* underneath the indenter [21]. The importance of the nonlocal nature of this problem has also been emphasized in the work of Delph et al. for the development of the Wallace criterion for the prediction of cavitation and crack growth problems in face-centered cubic (FCC) solids [22,23]. Garg and Acharya also proposed a dislocation nucleation criterion based on the kinematic structure of the theory of field dislocation mechanics [24]. For dislocation nucleation, MR developed a predictive criterion based on the mesoscale atomic acoustic tensor [25]. According to this criterion, HDN is indicated by loss of positive definiteness of the mesoscale acoustic tensor. The mesoscale criterion is similar to a previously proposed Λ criterion [25–30]. Λ is the minimum eigenvalue of the local acoustic tensor, calculated on an atom-by-atom basis [25–28]. On the other hand, MR’s criterion is based on the acoustic tensor of a cluster of atoms.

For the non-local analysis described by MR, one has to judiciously choose the meso-region for calculation of the acoustic tensor. MR did not discuss the appropriate procedure to choose the meso-region. The size and location of the meso-region are needed to perform mesoscale analysis to predict dislocation nucleation in nanoindentation simulations. In this work we show that even if the meso-region size is much smaller than the embryo, it can still reveal the instability. However, the meso-region captures the full spatial extent of the embryo only if the meso-region is bigger than the embryo. It has been shown previously by Garg Maloney [31] that the embryo grows nonlinearly with the radius of curvature of the indenter tip. The length of the embryo, $\xi $, and the depth of the embryo from the surface, $Y*$, exhibit nonlinear scaling with indenter radius, $R:\xi \u223cR1/2$ and $Y*\u223cR3/4$. In this work, we use these scaling laws [31] to scale the size of the meso-region, $rmeso$, with the indenter tip radius to capture the defect.

An important question we address in this work is related to the location of the center of the meso-region: how far can the meso-region move from the embryo center to detect the instability? If the instability is local, one could naively assume that if the meso-region contains the minimum $\Lambda $ location (i.e., the embryo center), it can predict the instability. However, our results contradict this. We determine the critical section of the embryo that the meso-region should necessarily contain to detect the dislocation embryo.

We use the results for size and location of the mesoregion to develop a systematic procedure to use mesoscale analysis for identifying the dislocation embryo in nanoindentation simulations. Since mesoscale analysis can identify instability much sooner than full scale analysis, this highlights the utility of mesoscale analysis in terms of computational time and resources for atomistic simulations. We use the mesoscale analysis in full 3D nanoindentation simulations of FCC crystals to verify the scaling laws given by Garg and Maloney [31] for the embryo size and location with indenter radius.

The rest of this paper is organized as follows. In Sec. 2, we describe the details of our simulations and give a brief kinematic description of the mechanism of homogeneous dislocation nucleation. Section 3 describes our rationale and results of the mesoscale analysis. In particular, we discuss the significance of mesoregion size by fixing the analysis at the center of the embryo, and then vary the mesoregion spatially and temporally. Section 4 describes the three-dimensional nanoindentation simulations of FCC crystals to verify the scaling laws for the embryo size and location for homogeneous dislocation nucleation. Section 5 concludes with a brief summary and discussion on implications of our results, together with an outline of our future work.

## 2 Modeling and Eigenanalysis

We perform athermal quasistatic nano-indentation simulations for 2D hexagonal Lennard Jones (LJ) crystals via energy minimization dynamics. The large-scale atomic/molecular massively parallel simulator (LAMMPS) molecular dynamics framework [32] is used to perform nonlinear energy minimization with the conjugate-gradient algorithm of the Polak–Ribere type, and a custom python code was developed to perform eigenmode analysis using the sparse matrix routines in SciPy. The indenter moves perpendicular to the nearest neighbor axis in hexagonal lattice for all simulations discussed in this paper. The resulting load–displacement and energy–displacement curves for the indentation process are shown in Fig. 1.

*L*

_{x}was chosen to ensure the results are independent of its value. In the rest of the document, all lengths such as

*L*

_{x},

*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. [21,33]. The interaction potential used to model the interaction between the indenter and the atoms is of the form

*R*is the radius of the indenter and

*r*is the distance of the particle from the indenter center.

*U*(

*x*

_{iα},

*D*), is a function of atomic positions,

*x*

_{iα}, and indenter depth,

*D*. The first derivative of energy with respect to the particle position gives the force,

*F*

_{iα}, on each particle as

*H*

_{iαjβ}, is the second derivative of total potential energy as follows:

The lowest four eigenvalues of the Hessian matrix for the relaxed configuration are calculated at each indenter step to identify the reaction coordinate. The system is driven to instability along a single eigenmode, which shows antiparallel motion of small number of atoms on adjacent crystal planes, resulting in the nucleation of a dislocation dipole as shown in Fig. 3(b). We denote the critical depth, *D*_{c}, as the depth of the indenter and the critical force, *F*_{c}, as the indenter force at which the dislocation dipole nucleates. Then, *δD* = *D* − *D*_{c} is the distance of the indenter from the critical depth, and *δF* = *F* − *F*_{c} is the remaining force of the indenter required for nucleation. We call the eigenmode corresponding to the eigenvalue vanishing as *δD*^{0.5}, the critical mode. In Fig. 3(a), the critical mode is along the dotted line. We analyze the critical mode at *δD* ≈ 10^{−6}, when the lowest mode is the critical mode as shown in Fig. 3(b).

*s*is defined as the distance of each triangle from the triangle corresponding to the maximum Ω along the crystal axis of interest. Figure 4(b) shows Ω(

*s*) profiles for

*L*= 160 and

*R*= 40. To define the size of the embryo, these Ω(

*s*) profiles are fit to Gaussian functions of the form

*s*) curves upto half their maximum value using a generalized least-square method and calculate

*ξ*, the embryo size.

## 3 Mesoscale Analysis of Incipient Dislocation

In this section, we study the critical eigenmode of the Hessian matrix of a circular mesoregion of radius, *r*_{meso}, of a cluster of atoms in the crystal. For an undeformed configuration, the lowest energy eigenvalue of an undeformed mesoregion scales as $1/rmeso2$, where *r*_{meso} is the size of the mesoregion as shown in the top panel of Fig. 5. This is expected for an isotropic, linear, elastic 2D sheet of atoms in a crystal with fixed boundaries. The critical eigenmode for an undeformed region of radius *r*_{meso} = 8 is a longwavelength mode as shown in the bottom left panel of Fig. 5.

Just before nucleation, the critical eigenmode for a mesoregion containing the embryo is localized as shown in the bottom right panel of Fig. 5. The mesoregion critical eigenmode, defined as the mesomode, looks similar to the full system eigenmode shown in Fig. 3(b). The lowest eigenvalue corresponding to the critical eigenmode falls off faster than any power law as shown by the curve in the top panel of Fig. 5.

A mesoregion of radius as small as six interatomic distances can capture the defect as shown by the mesomode in Fig. 6. Above a critical radius of the meso-region centered at the center of the incipient dislocation, defined by *hot atom*, the structure of the meso-mode saturates. In Fig. 6, the mesomode for *r*_{meso} = 32 has the same spatial structure as the full system eigenmode. The energy of the mesoregion or the lowest eigenvalue initially decreases as the mesoregion radius increases and then plateaus similar to the spatial structure of the mesoregion. Figure 7 is a log–lin plot for the mesoregion eigenvalue versus *r*_{meso} for different *R*. The height of the plateau is not a function of *R*; it is governed by the distance from nucleation, *δD*, for the mesoscale analysis. On the other hand, the structure of the localized mesomode is independent of the proximity to nucleation. Therefore, the rest of this work is focused on the structure of the mesomode versus the eigenvalue associated with the mesoregion.

### 3.1 Mesoregion Centered at the Hot Atom.

We compute the eigenmode corresponding to the lowest energy eigenvalue, defined as mesomode, for various *r*_{meso} with mesoregions centered at the embryo. Figure 6 shows the mesomode corresponding to *r*_{meso} = 6 and *r*_{meso} = 32 for *L* = 160, *R* = 120, and the critical portion of the eigenmode for full system (*r*_{meso} = ∞). The Ω(*s*) profiles corresponding to these modes is shown in Fig. 8. The width of the Ω curve for *r*_{meso} = 8 is smaller than that of the full system. The Ω curves converge when *r*_{meso} is greater than a critical value. As expected, they converge for $rmeso\u2273\xi $; however, mesoregions smaller than *ξ* definitely have information about the incipient instability.

We next measure the spatial extents, *ξ*_{meso}, of the mesoregion eigenmodes for all our systems. *ξ*_{∞} is the spatial extent obtained from the full system eigenmode. Figure 9(a) shows *ξ*_{meso} ∼ *r*_{meso} for different indenter radii. As expected from Fig. 9(a), for small *r*_{meso}, *ξ*_{meso} increases before plateauing at *ξ*_{∞}. The increase of *ξ*_{meso} with *r*_{meso} for mesoregions smaller than *ξ*_{∞} follows a power law. The insets show the collapse of these curves when rescaled by their corresponding plateau value: *ξ*_{∞}. The collapsed rescaled curve in Fig. 9(b) is universal for all indenter-crystal geometries. The plateauing of *ξ*_{meso} observed for *r*_{meso} > *ξ*_{∞} can be reasonably expected from the fact that those mesoregions completely encompass the incipient dislocation.

The power law scaling of $\xi meso$ suggests that for $rmeso$*smaller* than $\xi \u221e$, the spatial structure of the lowest meso-mode does reveal the structure of the critical mode. However, for meso-regions of radius of one or two particle spacings, the meso-mode invariably ceases to have the structure typical of an incipient dislocation. Although even small meso-regions host a lowest mode that resembles a dislocation embryo, it is only for $rmeso\u2265\xi \u221e$, that the meso-region gives a reliable characterization of the embryo size. It is in this sense that we define that nucleation is a *quasi-local* phenomenon. We define a quasi-local instability as one for which strictly atom-wise quantities (the atomic level virial or quantities derived from it such as the atomic level acoustic tensor) do not determine the stability of the system, but one for which the analysis of a finite region will asymptotically determine the stability of the system as the spatial extent of the analysis region grows. This scenario, where we can pick up the signature of an incipient dislocation with meso-regions smaller than its spatial extent , but not too small, indicates that homogeneous dislocation nucleation is a *quasi-local* phenomenon.

### 3.2 Mesoscale Analysis Centered Off From the Highest Ω Triangle.

In this section, we analyze the lowest eigenmode for the mesoregion, called *mesomode*, when the mesoregion is not centered at the embryo center. We define *critical mesomode* as the mode that contains the signature of an incipient dislocation. In Fig. 10, the mesomode is shown when the mesoregion is moved spatially along the slip plane. Upto a critical distance, *δ*_{c}, the mesomode is the critical mesomode. At *δ*_{c}, the long wavelength mode has lower energy than the critical mesomode, and the critical mesomode no longer remains the mesomode. The important question we address here is: How *δ*_{c} scales with the radius of mesoregion, *r*_{meso}, and the embryo size, *ξ*?

We perform mesoscale analysis centered at atoms around the embryo for different embryo sizes. In Fig. 11, the size of the mesoregion is fixed, *r*_{meso} = 40. The colors at each atom represent maximum Ω for the mesomode. The mesomode is found to be the critical mesomode when the mesoregion is centered on the red atoms. The area formed by the red atoms depends on *r*_{meso} and *ξ*.

*δ*

_{c}is the critical distance from the center of the embryo upto which the critical mesomode has a lower energy than the long wavelength mode. If we assume that the mesoregion should contain some critical portion of the embryo as shown in Fig. 12,

*δ*

_{c}can be written as in Eq. (5), where

*x*is the critical portion of the embryo.

*r*

_{meso}≫

*x*,

*θ*= 0,

*x*, we move our mesoregion along the slip plane for various

*r*

_{meso}and

*ξ*as shown in Fig. 13(a). For

*r*

_{meso}≫

*ξ*,

*x*is constant with varying

*r*

_{meso}and increases with an increase in

*ξ*(or

*R*). From Fig. 13(b),

*x*is found to be 1.2

*ξ*. This value of

*x*is substituted back in Eq. (6). Then, in Fig. 11,

*δ*

_{c}based on Eq. (6) is shown by the black curves for different

*R*(or

*ξ*). The black curves predict the area formed by red atoms fairly well. This implies if

*r*

_{meso}>

*ξ*, then it must encompass roughly 120$%$ of

*ξ*to be the critcal mesomode.

### 3.3 Mesoscale Analysis Away From Nucleation (Temporally).

Here, we address at what critical indenter depth, *D*_{meso}, the meso-mode becomes localized. *F*_{meso} is the indenter force at *D*_{meso}, the indenter depth at which the meso-mode can detect the incipient dislocation. In other words, at what *D*_{meso} or *F*_{meso} does the energy of the critical meso-mode become lower than the long wavelength meso-mode. *δ*_{Dmeso} is defined as the indenter depth difference between $Dc$ and *D*_{meso}. In other words, higher *δ*_{Dmeso} implies sooner prediction of incipient dislocation by the meso scale analysis. We study the dependence of *δ*_{Dmeso} with *r*_{meso}.

*λ*, of a critical meso-mode decreases with distance to dislocation nucleation,

*, as:*

*δ*D*λ*, is a function of the size of the meso-region,

_{e}*r*

_{meso}, given by

*D*

_{meso}, or critical indenter force,

*F*

_{meso}, the energy of the critical meso-mode becomes equal to the long-wavelength mode energy. This gives the analytical relation between

*δ*

_{Dmeso}and

*r*

_{meso}:

*δ*

_{Dmeso}and

*F*

_{meso}versus

*r*

_{meso}curves from simulations as shown in Fig. 14. The power law derived in Eq. (10) is also observed in simulations.

*versus*

*δ*_{D}*r*

_{meso}curves collapse when

*r*

_{meso}is scaled by the intrinsic embryo size, $\xi $, and

*F*

_{meso}versus

*r*

_{meso}collapse when

*r*

_{meso}is scaled by $\xi 1.5$. This can be derived from the fact that hardness,

*F*

_{meso}scaled by

*R*, is the critical quantity equivalent to

*δ*

_{Dmeso}and $\xi $ increases as a square root of

*R*. From these curves it can be seen that as the system size approaches infinity,

*δ*

_{Dmeso}goes to zero. In other words, the meso-scale analysis predicts nucleation much earlier than the full eigenmode analysis.

## 4 Scaling Laws in Full Three-Dimensional Simulations Using Mesoscale Analysis

We perform computational nano-indentation on a FCC lattice using a spherical indenter of radius *R* as described in Sec. 2. At each indenter step, the Hessian matrix is computed for the FCC lattice, similar to the two-dimensional simulations. Diagonalization of the Hessian matrix for a three-dimensional system can be computationally very expensive. Some of the systems considered for this work contain million atoms approximately. We use mesoscale analysis described in Sec. 3 to compute the critical mesomode. First, we measure the embryo size, *ξ*, and the embryo location, *Y**, for a smaller system size, *L*, and a smaller indenter radius, *R*, using full system kinematic analysis. Then, the results for the smaller *R* are extrapolated to calculate *ξ* and *Y** for the required *R*, assuming *ξ* ∝ *R*^{0.5} and *Y** ∝ *R*^{0.5}. These scaling laws for *ξ* and *Y** were shown previously by Garg and Maloney [30] that the embryo grows in a nonlinear way with the radius of curvature of the indenter tip. The length of the embryo, *ξ*, and depth of the embryo from the surface, *Y**, exhibit nontrivial scaling with indenter radius, *R*: *ξ* ∼ *R*^{1/2} and *Y** ∼ *R*^{3/4}. In this section, they are shown to hold true for the FCC lattice in full three-dimensional simulations.

It was observed in Secs. 3.1 and 3.2 that mesoregions, *r*_{meso} ≈ 6, can capture the embryo. Also, it was shown in Sec. 3.3 that the mesoregion not centered at the embryo can also capture the incipient dislocation. To search for the embryo, the mesoscale analysis was performed only once in each spherical zone of a radius of six atomic units. Since the mesoregions used were much smaller than full system size (*r*_{meso} ≈ 6), and the analysis was performed at a fewer locations, locating the embryo was computationally cheaper than full-scale eigenmode analysis. Once the embryo was located, we used a mesoregion of radius *r*_{meso} ≈ 1.5*ξ* centered at the embryo to capture the structure of the critical eigenmode. We utilized the result obtained in Sec. 3.3 that *r*_{meso} ≈ *ξ* is needed to capture the structure of the embryo. It was verified that increasing *r*_{meso} did not affect the results significantly.

The critical mesomode for *L* = 65 and *R* = 65 is shown in Fig. 15. At the onset of instability, there are two planes slipping with respect to each other resulting in the nucleation of a dislocation loop. The critical slip plane is (1 1 1) and the slip direction is 〈1 2 1〉. As shown, the embryo structure located on (1 1 1) is hexagonal. We triangulate the (1 1 1) plane containing the embryo using Delaunay triangulation. Each triangle and the vertex atom on the adjacent plane form a tetrahedron. The critical mesomode is interpolated on these tetrahedrons using linear finite element shape functions. Using this interpolation, Ω is defined as the curl of the critical mesomode resolved in the direction perpendicular to the slip plane, as done in Sec. 3.2 for the two-dimensional simulations.

We call the three diagonal axes of the hexagonal embryo as *n*_{1}, *n*_{2}, and *n*_{3}, as shown in Fig. 16. Ω versus *s* is plotted in Fig. 16 to obtain the embryo size (as done in Sec. 3.2). Interestingly, the Ω versus *s* curves along the three diagonal axes collapse, implying that the embryo is regular hexagon. The collapse of these curves was observed for all indenter radii. The embryo size, Ω, can be calculated from these Ω versus *s* curves by fitting Gaussian profiles as done in Sec. 3.2.

*Y**, and the embryo size,

*ξ*(Fig. 17).

*R*limit, but there was a systematic indentation size effect observed where $Fc/R$ increases at small R. $\xi /R0.5$ also exhibited the same type of indentation size effect as $Fc/R$ where $\xi /R0.5$ increases at small

*R*for all systems. Similar indentation size effect is also observed in the three-dimensional simulations at small

*R*.

## 5 Discussion and Summary

In this work, we studied the nucleation of dislocation dipoles in the bulk of perfect 2D and 3D crystals subjected to nanoindentation with a circular (and spherical in 3D), atomistic indenter under athermal, quasistatic conditions. We performed mesoscale analysis of configurations at the stability threshold, showing that small, nonlocal regions of the crystal, centered at the embryo, contain significant information about an incipient nucleation event. However, unlike previous work that utilized the minimum eigenvalue of the mesoregions as the main analysis tool [21–23], we focused our attention on the spatial structure of the lowest mesoregion eigenmode. We found that the relation between *ξ*_{meso}, the embryonic size inferred from the mesoregion, and *r*_{meso}, the size of the mesoregion itself, is universal. The lowest mesomode and eigenvalue were found to provide excellent estimates of the structure and energy of the true critical mode but only for mesoregions larger than *r*_{meso} > 1.5*ξ*. We also showed that the mesoregions not centered at the embryo center also reveal the presence of embryo, provided they encompass 1.2*ξ*. This scenario leads us to think of homogeneous dislocation nucleation as quasilocal: full information about the nature of the embryo can only be obtained by analyzing sufficiently large regions; however, its existence can be inferred by examining regions much much smaller than its intrinsic size.

We have also understood the effects of the proximity to nucleation, *δD*, on mesoscale analysis. Mesoscale analysis detects the embryo presence much before that the full system kinematic analysis. Mesoscale analysis can be used to make atomistic simulations computationally efficient, especially for 3D simulations where full system kinematic analysis could be challenging or, in some-cases impossible. We showed an initial study for full 3D FCC crystal nano-indentation simulations using the mesoscale analysis. The critical mesomode was used to calculate the embryo size, *ξ*, and the embryo location, *Y**. *ξ* and *Y** vary as the square-root of indenter radius for 3D FCC lattice. These scaling laws are consistent with the two-dimensional hexagonal lattice simulations for various interatomic potentials and crystal orientations. There are two main points to make regarding the overall utility of our results. As a practical consideration, it is dramatically more numerically expedient to employ a meso-scale eigenmode analysis as opposed to a global one. Miller and Rodney first proposed this approach, but understanding its limits is therefore important. Second, we address a fundamental question about the nature of the critical mode and how it evolves as one approaches the limit of stability. We find it likely that, in the future, when attempts to derive a theory for, e.g., activation rates for thermally assisted nucleation in a crystal in a state of inhomogeneous strain, the meso-scale properties we study here and the associated notion of locality/non-locality should be important.

The ultimate use of the analysis presented here would be to inform coarser-grained models, which do not explicitly take into account the atomic degrees of freedom, about the creation of new dislocations out of the void. For example, in field dislocation mechanics [36,37], one introduces a continuous field to represent the dislocation density. It is our hope that a criterion for dislocation nucleation based on a mesoscale analysis like we presented here could serve as a guide for the introduction of atomistic details at the dislocation embryo in concurrent multiscale schemes built on field theories like field dislocation mechanics.

## Acknowledgment

We thank Amit Acharya for useful discussions at various stages during this work. This material is based upon the work supported by the National Science Foundation (Award Number CMMI-1100245).