Site dependence of surface dislocation nucleation in ceramic nanoparticles

The extremely elevated strength of nanoceramics under compression arises from the necessity to nucleate highly energetic dislocations from the surface, in samples that are too small to contain pre-existing defects. Here, we investigate the site dependence of surface dislocation nucleation in MgO nanocubes using a combination of molecular dynamics simulations, nudged-elastic-band method calculations and rate theory predictions. Using an original simulation setup, we obtain a complete mapping of the potential dislocation nucleation sites on the surface of the nanoparticle and find that, already at intermediate temperature, not only nanoparticle corners are favorable nucleation sites, but also the edges and even regions on the side surfaces, while other locations are intrinsically unfavorable. Results are discussed in the context of recent in situ TEM experiments, sheding new lights on the deformation mechanisms happening during ceramic nanopowder compaction and sintering processes.


INTRODUCTION
It has long been recognized that materials are stronger when smaller 1,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 nanoscale [3][4][5][6][7] .
Size effects have been investigated in both poly-and singlecrystalline 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) simulations [8][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 semiconductors 5,[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 materials 6,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) method [21][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 deformation 13,[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 studies [27][28][29][30][31][32] . MgO has a B1 crystalline structure (a binary face-centered cubic) with 1 2 h110i{110} and 1 2 h110i{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 nanomechanics 5,6,32,34,35 . Issa and collaborators have investigated the mechanical properties of MgO nanocubes under compression using in situ TEM and MD simulations 6 . After the elastic loading, the MgO nanocubes yielded in the GPa range by the nucleation from the surface of perfect dislocations in the 1 2 h110i{110} slip 1 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 1 2 h110i{100} dislocation was observed, likely due to the native cubic orientation of the sample.

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 strainrate dependence of the process. A representative applied strain of 11% was chosen based on classical MD simulations performed at strain rate of 10 8 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 tests 6 . 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. 1 2 ½101ð101Þ, 1 2 ½101ð101Þ, 1 2 ½011ð011Þ and 1 2 ½011ð011Þ (the last two, 1 2 ½110ð110Þ and 1 2 ½110ð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 (L x = L z = 7.5 nm) and an aspect ratio of 2 (L y = 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 ffiffi 2 p 2 a 0 along [101], within a (101) slip plane characterised by the z coordinate h p i.e. the height at which the (101) 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 1 2 [101](101) dislocation half-loop, emitted from the (L, L, h p ) coordinates (Fig. 2a). The dislocation nucleates at the wire halflength 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 1 2 h110i{110} screw dislocations in MgO 30,33,36 .
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 y i from the original wire lateral surface (see Fig. 2a, c). Varying y i 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 (y i / L = 0, h p /L = 1) corresponds to a dislocation nucleation event from a top corner of the nanocube while (y i /L = 0.5, h p /L = 1) leads to dislocation nucleation from the middle of the top edge. We identify y i with the initial position of the nucleation centre. Other examples are shown in Fig. 2c where nucleation occurs in the middle (h p /L = 0.5) or in the lower part (h p /L = 0.25) of the nanocube back face. Note that we consider here a single slip system, 1 2 ½101ð101Þ, 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 1 2 ½101ð101Þ and 1 2 ½101ð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 h p /L = 0.75 in the nanocube are illustrated in Fig. 2b. For h p /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 y i /L computed configurations lead to only three different MEPs. For paths initially starting close to a side edge (y i /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 y i /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 middistance between both aforementioned nucleation sites. During the relaxation, the nucleation centre, which was initially located at y i may move to reach a relaxed position called y c , which is the true physical center of nucleation. A specific protocol was designed to extract y c 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, y i and y c , as well as the critical loop radius, r c , 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 r c = 1 nm (Fig. 3b). The rest of the top edge (0.18 < y i /L < 0.82) shows the emission of half loops with larger activation energies (0.81 ≤ ΔE ≤ 0.82 eV) and smaller radii (0.66 ≤ r c ≤ 0.77 nm).
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 (y i /L = 0 and 1) and mid-domains (y i /L = 0.5). Actually, ΔE values at mid-height (h p /L = 0.50) are amongst the highest computed (1.02 ≤ ΔE ≤ 1.12 eV), with r c 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 observed 6 . Overall, r cut 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 r c , 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 h p /L = 0.25 and 0.75. However, this symmetry fails in the lower part of the nanocube i.e. for h p /L < 0.15, where the NEB calculations do not converge due to the lack of available free volume (h p /L < 0.15 corresponds to less than four atomic rows on the side surface).  symmetry between the top and bottom surfaces of the nanocube observed in Fig. 1.
Except close to y i /L = 0 and y i /L = 0.5, the relaxed position of the dislocation nucleation centre y c shifts significantly away from its original position y i during the NEB relaxation, as shown Fig. 3c. In particular, y c /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 (y i ) and relaxed (y c ) centre positions in Fig. 4a, b, respectively. While the nanocube top edge exhibits only two preferential nucleation domains (at corners and mid-edge), lower h p /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 where ν 0 is the attempt frequency, ΔF the activation free energy and k B T the thermal energy scale. Here we use the approximation of a homogeneous surface disordering temperature ΔS = ΔE/T m , where ΔS is the entropy and T m 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 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.10 11 /s (as in ref. 24 ) and T m = 1550 K i.e. half the melting temperature of MgO 37 . We have checked that the results presented below do not depend significantly of our estimation of T m . The log 10 (k S1 + k S2 ) 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 10 7 -10 8 /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 experiments 6 . 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 (y c /L = 0, h p /L = 0.75), likewise in the middle of the top edge (y c /L = 0.5, h p /L = 1) where the cumulative rate is in the 1 /s range. The rest of the nucleation sites is irrelevant (k S1 + k S2 ≤ 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 compaction 38,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 multistep 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 r c 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 nanocompression 6 for which the multi-step NEB method will have to be applied at lower critical strains. Activation energy ΔE, dislocation radius r c 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.

Softwares and interatomic potentials
The LAMMPS code 40 is used for both NEB and MD simulations. Atomic structures are designed using ATOMSK 41 and visualised using OVITO 42 . Atomic interactions are described using the rigid ion model of Ball and Grimes 43 successfully employed to describe MgO surface diffusion 43 , elastic constants and dislocation properties 44 as well as nanocube compression tests using classical MD 6 . 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 Hardy 45 , 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 10 8 /s using a timestep of 2 fs and the Nosé-Hoover thermostat 46 . 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 Fig. 4 Mapping of the surface dislocation nucleation energy and rate for an MgO nanocube at ε = 11%. a, b show the 0 K activation energy ΔE as function of the initial and relaxed nucleation centre positions, respectively y i /L and y c /L. Both plots consider only one slip system S1 on each face. c, d show the cumulative dislocation nucleation rate log 10 (k S1 + k S2 ) 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. 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 on 22 . Finally, the last NEB calculation is carried out down to a force norm of 10 −2 eV/Å.

DATA AVAILABILITY
The datasets generated and/or analysed during the current study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
MD and NEB codes are part of the LAMMPS package available at https://lammps. sandia.gov. Homemade routines generated during the current study are available from the corresponding author upon reasonable request.