Introduction

It has long been recognized that materials are stronger when smaller1,2. The reason is a size-dependent competition between deformation and fracture that has attracted an intense scientific attention over the past 20 years. Metals have been mostly studied, but strikingly, size-induced toughening was also reported in originally brittle materials such as ceramics and semi-conductors. For instance, silicon and silicides, carbides and more generally oxides have shown unprecedented mechanical properties when probed at the nanoscale3,4,5,6,7.

Size effects have been investigated in both poly- and single-crystalline systems. With decreasing dimensions, bulk polycrystals exhibit the well-known Hall-Petch hardening effect down to a grain size of a few tens of nm, below which softening occurs. Similarly, reducing the size of single-crystal pillars or wires induces strengthening due to two size-dependent contributions: dislocation line tension and dislocation nucleation, both originally observed in metals using transmission or scanning electron microscopy (TEM and SEM), molecular dynamics (MD) and dislocation dynamics (DD) simulations8,9,10,11,12,13. Several mechanical studies performed on metal micropillars (originally containing dislocations) have shown that the size effect relies on preventing the dislocation multiplication process when samples are probed at the microscale (see e.g. ref. 14 and references therein). A similar smaller is stronger trend was observed in ceramics and semi-conductors5,15,16,17. By way of contrast, the strength of nanosamples (size below c.a. 100 nm) such as nanowires and nanoparticles obeys a different rule. Usually, nano-objects are assumed to be dislocation-free (or scarce) as they benefit from soft fabrication routes (e.g. molecular epitaxy or dewetting) as well as image forces that drag dislocations towards the sample surfaces, in contrast with FIB milled samples that contain larger amounts of subsurface defects. Assuming an initial lack of defects, the yield of nano-objects is not due to dislocation multiplication but rather to dislocation nucleation, in a fundamentally different manner compared to macroscopic materials6,11,12,18,19. Consequently, nano-objects are generally able to withstand stresses up to a significant fraction of their ideal strength before plastic relaxation.

In the case of sharp nano-objects, the dislocation nucleation process often initiates from surfaces as commonly observed in atomic-scale simulations (see e.g. ref. 20 for a recent review). Also, minimum energy path (MEP) calculation methods such as the Nudged Elastic Band (NEB) method21,22,23 have been used to compute energy barriers for dislocation nucleation in order to provide reliable mechanical models of incipient plasticity. Indeed, when coupled to the transition state theory (TST), this approach predicts the yield strength of nano-objects on a wide range of temperatures and strain rates that fill the usual gap between classical MD and experimental conditions of deformation13,24,25,26. Nevertheless, these studies are generally restricted to metal nanowires in which only one or a few pre-selected nucleation sites are considered. Depending on the shape and deformation conditions (temperature and strain rate), other nucleation sites might become relevant. Here, we propose a comprehensive study of the site dependence of surface dislocation nucleation and its implication in nano-object mechanics.

In this study, we focus on dislocation nucleation at the surface of magnesium oxide (MgO) nanoparticles. MgO is a well-known refractory material often used as a reference ceramic for mechanical studies27,28,29,30,31,32. MgO has a B1 crystalline structure (a binary face-centered cubic) with \(\frac{1}{2}\langle 110\rangle\){110} and \(\frac{1}{2}\langle 110\rangle\){100} slip systems, typical of cubic oxides. Bulk MgO critical resolved shear stresses are respectively of about few tens and few hundreds MPa at room temperature, where the lattice friction restricts dislocation glide (see e.g. ref. 33 for more details). With the development of modern techniques such as nanoindentation, multi-scale simulation methods and in situ mechanical tests, MgO has become the object of a renewed interest in the fields of micro- and nanomechanics5,6,32,34,35. Issa and collaborators have investigated the mechanical properties of MgO nanocubes under compression using in situ TEM and MD simulations6. After the elastic loading, the MgO nanocubes yielded in the GPa range by the nucleation from the surface of perfect dislocations in the \(\frac{1}{2}\langle 110\rangle\){110} slip systems. Several nucleation sites (corners, edges and surfaces) were observed in both the experiments and simulations (see examples shown here in Fig. 1 or in Fig. 3 of ref. 6). Dislocation nucleation events initiating from within top and bottom surfaces were only rarely observed and were interpreted as due to localized contacts (presence of defects, roughness, misalignment) between the indenter, the substrate and the sample. No \(\frac{1}{2}\langle 110\rangle\){100} dislocation was observed, likely due to the native cubic orientation of the sample.

Fig. 1: Surface dislocation nucleation in MgO nanocubes and slip systems.
figure 1

a MgO nanocube compression modeled using MD. A perfect \(\frac{1}{2}\langle 110\rangle\){110} dislocation nucleates from a top corner (blue line), just below the indenter. The red arrow refers to the Burgers vector direction. b Surface dislocation nucleation in a MgO nanocube observed during in situ TEM nanocompression test at room temperature (courtesy of Pr. Masenelli-Varlot). c \(\frac{1}{2}\langle 110\rangle\){110} slip systems in MgO.

Full size image

Results and discussion

Modeling approach

Here, surface dislocation nucleation in MgO nanocubes is investigated using a multi-step modeling approach based on the NEB method. We focus on the competition between potential nucleation sites under a given constant load and not on the strain-rate dependence of the process. A representative applied strain of 11% was chosen based on classical MD simulations performed at strain rate of 108 s−1 and at temperatures of 5 and 300 K (see e.g. Fig. 1a). MD simulations show dislocation nucleation at a critical applied strain of about 11–12% corresponding to a compressive strength in the 35–40 GPa range, depending on the nanocube size and temperature. NEB calculations are then performed on the nanocubes at the MD critical strain in order to map the spatial distribution of energy barriers for dislocation nucleation on the nanocube surfaces. The study considers (i) sharp and defect-free MgO nanocubes and (ii) surface dislocation nucleation events from lateral surfaces only (including edges and corners), in accordance with MD and previous experimental tests6. Due to the geometric and crystallographic symmetries of the nanocube, we only need to consider nucleation from one side surface and in one slip system to obtain a complete description of all the potentially activated slip systems i.e. \(\frac{1}{2}[101](\bar{1}01)\), \(\frac{1}{2}[10\bar{1}](101)\), \(\frac{1}{2}[011](0\bar{1}1)\) and \(\frac{1}{2}[01\bar{1}](011)\) (the last two, \(\frac{1}{2}[110](\bar{1}10)\) and \(\frac{1}{2}[1\bar{1}0](110)\), having a zero Schmid factor) from all the four side surfaces (see Fig. 1c).

In the first part of the NEB protocol described in Fig. 2a, a 〈100〉-oriented MgO nanowire (later used to carve out cubic-shaped nanoparticles), with a square cross-section (Lx = Lz = 7.5 nm) and an aspect ratio of 2 (Ly = 15 nm), is constrained between two planar force fields parallel to the z-axis to model both the substrate and a flat punch indenter. The indenter is displaced against the top of the nanowire up to an applied strain of ε = 11%, equivalent to the critical strain computed in the MD simulations. Free boundary conditions (FBCs) are used in the two other directions and the system is relaxed to produce the initial NEB configuration. To construct the final configuration, the wire is then sheared by \(\frac{\sqrt{2}}{2}{a}_{0}\) along [101], within a (\(\bar{1}\)01) slip plane characterised by the z coordinate hp i.e. the height at which the (\(\bar{1}\)01) slip plane intersects the nanowire side surface. At this stage, a preliminary NEB calculation is performed between the elastically compressed wire and its sheared counterpart. The saddle configuration along the MEP is characterised by the nucleation of a \(\frac{1}{2}\)[101](\(\bar{1}\)01) dislocation half-loop, emitted from the (L, L, hp) coordinates (Fig. 2a). The dislocation nucleates at the wire half-length due to symmetry. While the critical shape of the half-loop is close to semicircular, it later elongates along the edge character direction. This behaviour is attributed to the line tension, which is minimum in the screw direction, and to the Peierls lattice friction, which is known to restrict the mobility of \(\frac{1}{2}\langle 110\rangle\){110} screw dislocations in MgO30,33,36.

Fig. 2: Modeling surface dislocation nucleation in MgO nanocubes using a multi-step NEB approach.
figure 2

a Illustration of the first NEB step performed on a 2L long MgO nanowire with a square section of edge L constrained by two force fields i.e. the top indenter and the substrate (top image). The saddle point configuration (bottom image) is characterised by a \(\frac{1}{2}\)[101](\(\bar{1}01\)) dislocation half-loop (line and Burgers vector are shown in blue and red, respectively) emitted at (L, L, hp), from the back surface. Only Mg atoms are shown for the sake of clarity (fcc and surface atoms are coloured, respectively, in green and light-grey). The original simulation cell of the nanowire (before loading) is shown using grey-solid lines. The grey-dashed lines refer to a typical nanocube, carved out to compute the site dependence of the surface dislocation nucleation energy barriers during the second NEB step. b Examples of surface dislocation nucleation energy barriers in a nanocube for yi/L ranging from 0 to 0.5 and hp/L = 0.75. c Examples of configurations for dislocation nucleation events from a corner (yi/L = 0, hp/L = 1), from the mid-top edge (yi/L = 0.5, hp/L = 1), the mid-side edge (yi/L = 1, hp = 0.5) and the lower mid-surface region (yi/L = 0.5, hp/L = 0.25).

Full size image

As a second step, the entire NEB chain is sliced along the nanowire length to produce nanocubes of edge L where the nucleation occurs at different positions. Each cube is cut out starting at a distance yi from the original wire lateral surface (see Fig. 2a, c). Varying yi from 0 to L leads to nanocubes with NEB paths where a dislocation nucleates at different lateral positions inside the nanocube. As illustrated in Fig. 2c, the parameter set (yi/L = 0, hp/L = 1) corresponds to a dislocation nucleation event from a top corner of the nanocube while (yi/L = 0.5, hp/L = 1) leads to dislocation nucleation from the middle of the top edge. We identify yi with the initial position of the nucleation centre. Other examples are shown in Fig. 2c where nucleation occurs in the middle (hp/L = 0.5) or in the lower part (hp/L = 0.25) of the nanocube back face. Note that we consider here a single slip system, \(\frac{1}{2}[101](\bar{1}01)\), and nucleation events from one lateral face of the cube. However, as mentioned above, thanks to the symmetries of the sample and of the slip systems, the same calculations would apply to the other faces with the corresponding slip systems. Moreover, for each face, a system and its orthogonal conjugate (for instance \(\frac{1}{2}[101](\bar{1}01)\) and \(\frac{1}{2}[10\bar{1}](101)\), see Fig. 1c) have symmetrical activation energies when counted from the top and bottom surface respectively.

After slicing, the initial NEB paths cut out from the nanowire into nanocubes need to be re-relaxed. Several examples of relaxed MEPs for dislocation nucleation centres along hp/L = 0.75 in the nanocube are illustrated in Fig. 2b. For hp/L = 0.75, the critical energy for surface dislocation nucleation at 0K, ΔE, varies from 0.80 to 1.11 eV and the relaxations of the nine yi/L computed configurations lead to only three different MEPs. For paths initially starting close to a side edge (yi/L equal to 0 and 0.06), the saddle configuration corresponds to a quarter loop emitted with an energy barrier ΔE = 0.80 eV, while a half-loop nucleates for yi/L = 0.50 with a higher energy barrier ΔE = 1.11 eV. In-between, all MEPs approximately converge to a saddle energy of ΔE = 0.98 ± 0.01 eV, characterised by the emission of a half-loop at mid-distance between both aforementioned nucleation sites. During the relaxation, the nucleation centre, which was initially located at yi may move to reach a relaxed position called yc, which is the true physical center of nucleation. A specific protocol was designed to extract yc from the MEPs as described in the Supplementary Methods. This first dataset allows us to draw a preliminary conclusion: several energy paths for dislocation nucleation with close activation energies but different nucleation sites co-exist in nano-objects.

Site dependence of surface dislocation nucleation

The full mapping of the surface dislocation nucleation activation energy is shown in Figs. 3 and 4. To characterize these events, we use four quantities, their activation energy ΔE, the initial and relaxed positions of the dislocation nucleation centre, yi and yc, as well as the critical loop radius, rc, which is computed in the activated state using a dedicated protocol described in the Supplementary Fig. 1. We see in Figs. 3a and 4a, b that even for a simple cubic nanoparticle shape, the spatial distribution of ΔE is particularly heterogeneous, including several subdomains. Lower activation energies are observed at the top edge (corners and mid-region) of the nanocube, below the indenter, which by symmetry is equivalent to the bottom edge with a nucleation in the orthogonal slip system. Both top corners show the lowest energy levels with 0.33 ≤ ΔE ≤ 0.36 eV and are characterised by the nucleation of quarter loops with a critical radius of about rc = 1 nm (Fig. 3b). The rest of the top edge (0.18 < yi/L < 0.82) shows the emission of half loops with larger activation energies (0.81 ≤ ΔE ≤ 0.82 eV) and smaller radii (0.66 ≤ rc ≤ 0.77 nm).

Fig. 3: Side-surface site dependence of dislocation nucleation in a MgO nanocube at ε = 11%.
figure 3

a 0 K activation energy ΔE, b critical dislocation radius rc and c relaxed position of the nucleation centre yc/L, as a function of the initial position yi/L.

Full size image
Fig. 4: Mapping of the surface dislocation nucleation energy and rate for an MgO nanocube at ε = 11%.
figure 4

a, b show the 0 K activation energy ΔE as function of the initial and relaxed nucleation centre positions, respectively yi/L and yc/L. Both plots consider only one slip system S1 on each face. c, d show the cumulative dislocation nucleation rate log10(kS1 + kS2) calculated using Eq. (1), respectively, at T = 300 and 1000 K. Both orthogonal slip systems S1 and S2 are considered on each face. The black grid shows the NEB calculation mesh. Forbidden sites are marked as hatched domains.

Full size image

The activation energy increases when moving away from the indenter, as shown by the ΔE interpolation on the side surfaces of the sample (Fig. 4a), especially at side edges (yi/L = 0 and 1) and mid-domains (yi/L = 0.5). Actually, ΔE values at mid-height (hp/L = 0.50) are amongst the highest computed (1.02 ≤ ΔE ≤ 1.12 eV), with rc between 0.83 and 0.91 nm. As a matter of fact, the energy cost at 0 K for surface dislocation nucleation from the middle of a side edge (ΔE = 1.12 eV) is more than three times larger than from a corner. This corroborates in situ TEM and MD results in which surface dislocation nucleation from the corners of MgO nanocubes (as shown Fig. 1) is the most frequently observed6. Overall, rcut ranges from 0.65 to 1.25 nm with extrema values located at the mid-surface and side edge domains. While ΔE mostly scales with the dislocation line length for half loops and large values of rc, quarter loops show a less regular behaviour that confirms the influence of the lateral free surfaces on the nucleation process. Besides, the ΔE distribution is mostly symmetric from the top to the bottom side of the nanocube i.e. ΔE = 0.80 ± 0.01 eV (side edge) and 1.10 ± 0.01 eV (mid-domain), for symmetric hp/L = 0.25 and 0.75. However, this symmetry fails in the lower part of the nanocube i.e. for hp/L < 0.15, where the NEB calculations do not converge due to the lack of available free volume (hp/L < 0.15 corresponds to less than four atomic rows on the side surface). ΔE data lines for hp/L = 0.15 and hp/L = 0.90 compare well with each other, whereas hp/L = 0.10 is the starting point of ΔE divergence. The few converged configurations computed at hp/L = 0.10 confirm the hardening of surface dislocation nucleation in \(\frac{1}{2}[101](\bar{1}01)\) when localized near the bottom of the nanocube (ΔE = 1.3 eV for yi/L = 0.5 and hp/L = 0.10). Note that, while dislocation nucleation in 1/2[101](\(\bar{1}\)01) is not possible near the bottom surface, nucleation is likely in the orthogonal \(\frac{1}{2}[101](\bar{1}01)\) system, whose activation energies are obtained by mirror symmetry with respect to hp/L = 0.5, thus recovering the symmetry between the top and bottom surfaces of the nanocube observed in Fig. 1.

Except close to yi/L = 0 and yi/L = 0.5, the relaxed position of the dislocation nucleation centre yc shifts significantly away from its original position yi during the NEB relaxation, as shown Fig. 3c. In particular, yc/L curves exhibit constant plateaus and steps that emphasize forbidden and preferential nucleation sites (see Supplementary Information for details on the definition of the forbidden nucleation sites). For comparison, ΔE maps are shown as a function of the initial (yi) and relaxed (yc) centre positions in Fig. 4a, b, respectively. While the nanocube top edge exhibits only two preferential nucleation domains (at corners and mid-edge), lower hp/L are characterised by additional nucleation sites, separated by forbidden regions where dislocation nucleation is unfavorable.

Dislocation nucleation rate and temperature

The probabilistic nature of the dislocation nucleation process at finite temperature (i.e. when the applied load is lower than the athermal strength) can be rationalised combining the previously discussed NEB results with the TST. We will use here a simple approximation, which, however, already clearly illustrates the competition between nucleation sites at finite temperatures. As a first approximation, the nucleation rate per site k at a given ε and T can be written as

$$k={\nu }_{0}{e}^{-\frac{{{\Delta }}F(\varepsilon,T)}{{k}_{B}T}},$$
(1)

where ν0 is the attempt frequency, ΔF the activation free energy and kBT the thermal energy scale. Here we use the approximation of a homogeneous surface disordering temperature ΔS = ΔE/Tm, where ΔS is the entropy and Tm the surface disordering temperature, as done in ref. 24. Thus, the free energy ΔF introduced in Eq. (1) can be expressed using the 0K activation energy ΔE as

$${{\Delta }}F(\varepsilon,T)={{\Delta }}E(\varepsilon )(1-\frac{T}{{T}_{m}}).$$
(2)

The spatial distribution of the dislocation nucleation rate is calculated combining Eqs. (1) and (2) as well as ΔE computed at ε = 11%. ν0 is estimated as 3.1011 /s (as in ref. 24) and Tm = 1550 K i.e. half the melting temperature of MgO37. We have checked that the results presented below do not depend significantly of our estimation of Tm. The log10(kS1 + kS2) rate that accounts for dislocation nucleation in both orthogonal slip systems S1 and S2 is shown in Figs. 4c, d, at T = 300 and 1000 K, respectively. While both figures show qualitatively similar patterns, the nucleation rate differs by several orders of magnitude depending on the temperature and site location. Both figures cover very different ranges of nucleation rates: 16 orders of magnitude at 300 K compared to 1 at 1000 K. At room temperature, dislocation nucleation is more likely to appear from the corners with a rate in the 107–108 /s range compatible with the MD timescale (a few hundred of ns). It is also in good agreement with in situ TEM observations (see e.g. Fig. 1b) even if MgO nanocubes yield at a lower compression stress in experiments6. However, the nucleation rate quickly decreases along the side edge but is still compatible with experimental timescales i.e. few events per seconds are expected for sites at (yc/L = 0, hp/L = 0.75), likewise in the middle of the top edge (yc/L = 0.5, hp/L = 1) where the cumulative rate is in the 1 /s range. The rest of the nucleation sites is irrelevant (kS1 + kS2 ≤ 10−3 /s) at room temperature. Increasing the temperature reduces the gap between nucleation sites. At intermediate T = 1000 K (less than 30% of MgO melting temperature), the nucleation rate varies within the same order of magnitude, whatever the site location (Fig. 4d). Thus, while corners remain the first choice for dislocation nucleation at low temperature, side edges, mid-top edge as well as mid-surface domains become quickly relevant nucleation sites already at a relatively low temperature for a ceramic. These results confirm the importance of a comprehensive description of the surface dislocation nucleation sites in nanoparticles. As a matter of fact, they have major implications in the processing of bulk nanostructured ceramics. Compaction and sintering are usually performed at higher temperature (about 2000 K for MgO powders) and it is usually considered that nanopowders only reorganize spatially during compaction38,39: here we confirm they should plastically deform. High temperature applications such as coating and thermal protection for which the mechanical history of grains is of prime importance to ensure functionality are also concerned.

In conclusion, the site dependence of surface dislocation nucleation in nanoparticles was analyzed by means of a multi-step NEB modeling approach, applied here to MgO nanocubes under compression at ε = 11%, close to the MD critical strain. In this study, we proposed an original method to probe multiple nucleation sites at the surface of a nanoparticle and provide a full mapping of the sample surface that shows wide spatial heterogeneity in terms of dislocation nucleation activation energy and rate when combined with a rate theory. Parts of guessed nucleation sites drift during the NEB relaxation and forbidden sites appear at the top edges and side surfaces. Moreover, while dislocation nucleation from the corners is the most favourable process at low temperature, sites at side and top edges, as well as mid-surface domains become active already at intermediate temperatures. In the simulation conditions of deformation (ε = 11%; σ = 38.6 GPa), dislocation nucleation critical radii rc are about 1 ± 0.25 nm. These results are in quantitative agreement with MD simulations and point out the importance of a holistic approach when applied to the study of potential dislocation nucleation sites in nano-objects. They could have a major impact for our understanding of bulk nanostructured ceramic processes such as compaction and sintering. Moreover, they provide a first hint at the interpretation of in situ TEM nanomechanical tests of MgO nanocompression6 for which the multi-step NEB method will have to be applied at lower critical strains. Activation energy ΔE, dislocation radius rc and spatial distribution of favourable nucleation sites can be used as reliable inputs for higher scale simulations as DD and crystal plasticity. We have here only discussed the site dependence of dislocation nucleation for a given applied load. We are currently applying the same methodology at varying applied loads in order to extend our study to the effect of strain rates.

Methods

Softwares and interatomic potentials

The LAMMPS code40 is used for both NEB and MD simulations. Atomic structures are designed using ATOMSK41 and visualised using OVITO42. Atomic interactions are described using the rigid ion model of Ball and Grimes43 successfully employed to describe MgO surface diffusion43, elastic constants and dislocation properties44 as well as nanocube compression tests using classical MD6. Cut-off and summation radii of 8 and 10 Å are used for short- and long-range interactions, respectively, the latter being computed using the multi-level summation method of Hardy45, with a relative error on the forces of about 10−5.

Nanocompression simulation method

Indenter and substrate are modelled using two harmonic force fields with a force constant of about 1000 eV/Å3 that characterises the stiffness of the indenter and bottom surface (see LAMMPS fix indent documentation for more details). This value is chosen high enough to limit the penetration of the indenter and substrate into the sample. In the MD, the top indenter is moved at a constant displacement rate leading to an engineering strain rate of 108 /s using a timestep of 2 fs and the Nosé-Hoover thermostat46. The reader can refer to ref. 6 for more details about the MD protocol.

Multi-step NEB protocol

In the following, we present the technical details of the multi-step NEB method. First, an elastically deformed (ε = 11%) MgO nanowire is minimized down to a force norm of 10−4 eV/Å, at 0 K temperature. This configuration is referred as the perfect strained nanowire (PSN) in the following. From the PSN, a sheared nanowire is generated as described in the text. 48 nanowire configurations are linearly interpolated between the PSN and its sheared counterparts. Among these configurations, we select the first one, which relaxes toward the sheared nanowire; it is then used as the end configuration of a a free-end NEB calculation using 12 replicas down to a force norm of 0.5 eV/Å (see LAMMPS fix neb documentation for more information). Both parallel and perpendicular nudging forces of 10 eV/Å are used. The target energy of the free end is set to the one of the initial configuration with a spring constant of 20 Å−1. Then, the NEB relaxation for the nanocubes is carried out in two steps. In the first step, we use a parallel spring constant of 0.5 eV/Å and a perpendicular one of 3 eV/Å. The target energy of the free end is set to the energy of the first replica (i.e. the perfect strained cube) with a spring constant of 1 Å−1. During that stage, we prevent intermediate replicas from having an energy lower than the first replica by removing gradient contribution for any replica having an energy lower than that of the first replica. This preliminary relaxation is ended after 10,000 relaxation steps or when the force norm becomes lower than 0.1 eV/Å. In the second step, the energy constraint on intermediate replicas is removed, a parallel spring constant of 0.1 eV/Å and a perpendicular one of 3 eV/Å are used. The same settings for the free end are kept and the climbing NEB is turned on22. Finally, the last NEB calculation is carried out down to a force norm of 10−2 eV/Å.