Adaptive hard and tough mechanical response in single-crystal B1 VNx ceramics via control of anion vacancies

1Department of Materials Science and Engineering, Cornell University, Ithaca, NY, 14853, USA 2Department of Materials Science and Engineering, University of California, Los Angeles, CA 90095, USA 3SKF Research and Technology Development Center, 3439 MT Nieuwegein, Netherlands 4Department of Physics, Chemistry and Biology (IFM) Linköping University, SE-581 83, Linköping, Sweden 5Department of Materials Science and the Materials Research Laboratory University of Illinois, 104 South Goodwin, Urbana, IL 61801 6Department of Materials Science, National Taiwan University of Science and Technology, Taipei 10607, Taiwan 7ICAMS, Ruhr-Universität Bochum, D-44780 Bochum, Germany


Introduction
Brittle fracture in ceramics is primarily caused by the limited ability of these materials to dissipate mechanical stresses ahead of a growing crack [1,2]. Accordingly, a strategy adopted to mitigate brittleness involves enhancing hardness while simultaneously hindering or deflecting crack propagation. This is commonly achieved through grain-boundary [3] and nanostructure engineering [4,5]. However, recent experimental results [6][7][8][9] demonstrated that single-crystal NaCl-structure (B1) pseudobinary V0.5Mo0.5N transition-metal (TM) nitride ceramics are intrinsically both hard (~20 GPa) and ductile. This surprising finding confutes the commonly accepted assumption that excellent ductility and high strength (well-correlated to hardness in solids [10,11]) are mutually-exclusive material properties [12,13]. Experiments have also proven that single-crystal B1 V0.5Mo0.5Nx solid solutions become much harder (from 17 to 26 GPa) when the concentration of anion vacancies increases up to 45% [7].
Previous experimental results show that the hardness (H) of B1 Group-VB (i.e. V, Nb, and Ta) nitrides and carbides is enhanced by anion vacancies [14][15][16][17], consistent with H vs. x trends in V0.5Mo0.5Nx [7][8][9]. Nonetheless, although understoichiometric V0.5Mo0.5Nx is considerably less prone to crack than single-crystal B1 TiN and B1 VN, it is more susceptible to fracturing than stoichiometric V0.5Mo0.5N [7]. This raises the question of whether hardness and ductility can be enhanced at the same time by controlling the metal/non-metal compositional ratio in refractory carbonitrides.
In this work, we combine experiments and ab initio molecular dynamics (AIMD) atomistic simulations to demonstrate that control of the N/V stoichiometry in VNx simultaneously enables enhanced hardness, ductility, and toughness. The experiments are conducted on high-quality single-crystalline VNx deposited epitaxially on MgO(001) with nitrogen-to-vanadium ratios x spanning 0.8 -1.0 [18][19][20]. Nanoindentation testing show that the unusual combination of remarkable mechanical properties in understoichiometric VN0.8in preliminary tests, VN0.8 exhibits the best performance among all investigated VNx systemsunambiguously originates from anion vacancies. AIMD simulations of VN and VN0.8 subjected to tensile deformation (up to fracture) and uniform shearing (up to lattice slip) at 300 K are used to unravel key atomistic and electronic mechanisms responsible for the superior mechanical behavior induced by vacancies. The results provide unprecedented insights for simultaneously enhancing hardness, ductility, and toughness in single-crystal TM carbonitride ceramics simply by controlling the lattice stoichiometry.

Experimental.
The B1 structure of VNx compounds has a large single-phase field that can accommodate compositional variations 0.7 < x ≤ 1, possible through incorporation of anion vacancies. VNx/MgO(001) thin films with x spanning from 0.8 to 1.0 are grown to a thickness of 300 nm under 20 mTorr mixed N2/Ar atmospheres in an ultra-high-vacuum magneticallyunbalanced reactive magnetron sputter-deposition system. Anion sublattice occupancy is controlled by varying the N2 gas fraction between 0.1 and 1.0 and the growth temperature between 430 and 540 °C [19]. Rutherford backscattering spectroscopy is employed to determined nitrogento-vanadium ratios to an accuracy of ±0.05. MgO is selected as substrate material because it is isostructural with VN and exhibits a sufficient film/substrate lattice mismatch (0.9%) to enable film relaxation, a necessity for the investigation of intrinsic properties.
VNx(001) hardness values as a function of x are quantified [21,22] via nanoindentation experiments performed on 300-nm-thick films in a Hysitron TI950 Triboindenter using a sharp Berkovich 142.3° diamond probe (tip radius ~150 nm) calibrated to an epitaxial TiN/MgO(001) sample. Nine indentations, arranged in a 3 × 3 pattern with indents separated by 10 µm, are made in each sample with the maximum tip penetration limited to 10% of the film thickness to avoid spurious substrate effects on measured hardness values. The fracture toughness of VNx samples is assessed via nanoindentation using a cube-corner diamond probe, which is sharper than a Berkovich tip and thus results in higher contact stresses [23]. We performed two different tests: indentations at constant depths of 400 nm (which exceed the film thickness by 100 nm) and indentations as a function of penetration depth between 25 and 400 nm, at steps of 25 nm.
That compositional variations in understoichiometric VNx are accommodated through anion vacancies is concluded based on the evolution of VNx relaxed lattice parameters a0 and mass densities  as a function of x. a0(x) values determined from x-ray diffraction (XRD) high-resolution reciprocal-space-maps decrease from 0.4134 nm (x = 1) to 0.4087 nm (x = 0.8)a trend which is accurately described with first-principles density-functional theory (DFT) results obtained as a function of nitrogen vacancy concentrations [19]. VNx/MgO(001) mass densities, obtained from xray reflectivity scans as well as Rutherford backscattering spectroscopy combined with film thickness measurements [20], decease from 6.06 (x = 0) to 5.98 g·cm-3 (x = 0.24) [20]. VN0.8 containing nitrogen vacancies is expected to have a mass density of 5.9 g·cm-3, in close agreement with the experimentally determined value; V interstitials and VN antisite substitutions would generate larger  values of 7.4 and 6.6 g·cm-3, respectively.

Modeling.
Tensile and shear mechanical testing of VNx systems are modelled using densityfunctional ab initio molecular dynamics [24] simulations at room temperature.
Tensile-testing simulations enable the evaluation of ideal tensile strength T and toughness UT for different strain directions. Deformation is applied orthogonal to (001) and (110) crystal surfaces, for which the energy of formation is lowest in B1-structure ceramics [25]. In AIMD modeling, T corresponds to the maximum zz stress reached during elongation, whereas UT is calculated as the area that underlies the full stress/strain curve. Shear-deformation modeling is used to determine the resistance to change of shape of the materials, which is well-correlated with their hardness [26][27][28]. In addition, these simulations shed light on the mechanisms that govern the activity of {110}〈11 ̅ 0〉 and {111}〈11 ̅ 0〉 slip systems, typically operative in B1-structure nitrides and carbides at room and/or elevated temperatures [29][30][31].
AIMD simulations are carried out with VASP [32] implemented with the projector augmented-wave method [33]. The Perdew-Burke-Ernzerhof (PBE) electronic exchange and correlation approximation [34], Γ-point sampling of the Brillouin zone, and planewave cutoff energies of 300 eV are used in all simulations. The ionic equations of motion are integrated at timesteps of 1 fs, using 10-5 eV/supercell convergence criteria for the total energy during both system equilibration and mechanical testing. Prior to modeling tensile and shear deformation, the supercell equilibrium structural parameters are determined via NPT simulations at 300 K using the Parrinello-Rahman barostat [35] and the Langevin thermostat. Thus, NVT sampling of the configurational space is used to equilibrate the unstrained structures for 5 additional ps using the Nose-Hoover thermostat. In these simulations, it is ensured that the time-averaged |xx|, |yy|, and |zz| stress components are ≲ 0.5 GPa.
Our model supercells are oriented in a convenient manner for investigating the response of Tensile testing is modeled for VNx(001) and VNx(110) supercells with 768 lattice positions (24 layers orthogonal to the strain direction; 32 sites per layer), following the procedure of Ref. [36,37]. Briefly, at each zz strain step (1% of the supercell vertical size), the systems are first coupled to an isokinetic, and then to a Nosé-Hoover, thermostat. The thermostat temperature is set to 300 K, for a total equilibration time of 3 ps. The structures are progressively elongated until they reach their fracture points. VN and VN0.8 simulation boxes used for shear deformation are formed of 576 lattice sites (24 layers parallel to the slip plane; 24 atomic positions per layer). The cells are equilibrated and sequentially deformed (shear strain steps xz=0.02), with procedure analogous to that used for tensile testing. In all VN0.8 supercells, the vacancy distribution on the anion sublattice yields negligible degrees of short-range ordering, which is consistent with experimental analysis [19].
Average tensile zz and shear xz stresses are determined as a function of strain from 500 equilibrated AIMD configurations. We note that, while in previous static first-principles stress/strain calculations thermal disorder was mimicked by introducing arbitrary atomic displacements [27,38], our simulations explicitly include temperature effects, thus providing a realistic description of phase transformations and non-linear elastic effects. The approach is particularly important for modeling crystal structures that are dynamically-stabilized by lattice vibrations, including B1 VN [39].
The simulation results are visualized using the VMD software [40]. Electron-transfer and energy-resolved electron-density maps are calculated for atomic configurations directly extracted from AIMD simulations at 300 K. Electron transfer is obtained by subtracting non-interacting atomic charge densities from the self-consistent electron density of the entire crystal. The electron densities in next-nearest-neighbor d-t2g and 4th-neighbor (across-vacancy) d-eg metallic V states are visualized by resolving the total charge density into the energy intervals [-2 -0 eV] and [-3 --2 eV], respectively, where 0 corresponds to the Fermi level (EF). In our single-crystal B1 VNx(001) samples, understoichiometry is accommodated by vacancies on the anion sublattice [19]. Nanoindentation testing demonstrates that VN0.8 is approximately 20% harder than VN (HVN0.8 = 17.1±0.8 vs. HVN = 14.0±0.8 GPa). This finding is consistent with previous experimental results, which have shown that the hardness of B1 Group-VB nitrides, including VNx, increases with a N content that decreases from 100 to ≈80% [14][15][16][17].

Results and discussion
The increased hardness of VN0.8 is attributed to reduced dislocation mobility. This is also observed indirectly through the incomplete strain relaxation (92%) of VNx layers with x ≥ 0.05 vs. 0 ≤ x ≤ 0.05 (97% nitrogen) [19]. Thus, nitrogen vacancies obstruct dislocation glide (at least for relatively low loads; see the results of AIMD shear deformation below) resulting in higher hardness.
Mechanical testing of single-crystal films demonstrates that VN0.8 is not only harder, but also tougher than the stoichiometric compound. Vacancy-induced toughening in understoichiometric VN0.8 layers is experimentally demonstrated by analyzing nanoindentation load-displacement curves (Fig. 2) and assessing indentation-induced fracturing using scanningelectron microscopy (SEM) and  Scanning probe microscopy (SPM) height isointensity contour maps show that plastic flow in indented VN samples occurs in the shape of a cross (Fig. 3d). Cube-corner indentation produces pronounced material pile-up (red areas in Fig. 3d) along 〈110〉 directions and less pile-up along 〈100〉. Mounds near the faces of the indentation triangle, where stress is concentrated, have heights (≈50 nm) which are more than twice those (≈20 nm) of hillocks near the vertices of the indentation triangle (Fig. 3d). Azimuthal rotations of the indentor with respect to the orientations shown in  GPa, Fig. 4b. Hence, AIMD results show that the ideal tensile strength T of VN0.8 is within ≈10% that of stoichiometric VN. However, VN0.8 exhibits a much higher toughness and resistance to fracture than VN, with the most remarkable differences observed for [001]-strained compounds (Fig. 4a).
At an elongation of 14%, VN(001) cleaves on the (001) plane by sudden and essentially simultaneous breakage of V-N bonds parallel to the loading direction ( Figs. 5a and 4a). In contrast, VN0.8(001) resists fracture up to an elongation of 30% (Figs. 5b and 4a). The crack develops between (001) lattice layers characterized by a relatively low concentration of anion vacancies. In addition, VN0.8 displays a total tensile toughness [48] which is approximately twice that of VN. The simulations reveal that the considerably enhanced resistance to fracture of VN0.8 is due to a transformation toughening mechanism, characterized by buckling of (001) atomic planes, activated by extreme tensile stress (Fig. 5b). Analogous to the results obtained for (001) supercells, VN0.8(110) also displays larger toughness than VN(110) (Fig. 4b) due to changes in bonding geometry activated at high tensile loading (AIMD snapshots not shown). The progressive transformation in the bonding network allows VN0.8 to prevent stress build up, thereby hindering crack nucleation and actively toughening the material. The mechanism is similar to the one observed during AIMD mechanical testing [36] of hard and ductile B1 V0.5Mo0.5Nx alloys [7]. We note that, within the supercell model considered here, the structural transformations identified in VN0.8 compounds are fully reversible upon relaxation; after cleavage on (001) and (110) 6 illustrates the shear stress xz calculated for VN and VN0.8 crystals subjected to [11 ̅ 0] shearing of (110) and (111) crystallographic planes. The peculiar upward bending of xz vs. strain xz curves observed for small (xz < 0.1) shear deformations of VN(110) and VN(111) (Fig. 6a-b) is due to the large third-order elasticity of VN [56] (as elaborated in the following paragraphs).
To elucidate the differences in measured hardness values, it is worth to closely analyze the elastic responses to deformation of B1 VNx crystals. Table 2 summarizes the results of the elastic constants and moduli evaluated by AIMD mechanical testing. The stiffnesses C11, G110, and G111u calculated from the slopes (≤2% deformation) of stress/strain curves presented in Figs. 4a, 6a, and 6b, respectivelyquantify the actual elastic resistances of B1 VNx crystals to 〈001〉 uniaxial tension, {110}〈11 ̅ 0〉 and {111}〈11 ̅ 0〉 shearing at room temperature. Table 2 also lists effective average elastic constants C 11 , C 12 , and C 44 obtained by minimizing the least-square differences between the stress components extracted during AIMD simulations and those predicted according That also implies validity of the equivalences C ij ≡C ij and that actual shear stiffnesses are accurately evaluated using the relations G110≡G 110 ≡(C11-C12)/2 and G111u≡G 111u ≡(C11-C12+C44)/3 [58,59].
Nonlinear elastic effects are characteristic for anharmonic materials, as B1 VN [62], which are dynamically stabilized by lattice vibrations at finite temperatures [39,63]. In the case of {110}〈11 ̅ 0〉 shear deformation, the initial upward curvature of the stress/strain dependence of VN ( Fig. 6a) can be rationalized as a deformation-induced softening of transversal acoustic phonon modes. The {110}〈11 ̅ 0〉 shear deformation effectively corresponds to a pure tetragonal distortion of the cubic lattice (see illustration in figure 4 of [27]). The distortion is energetically facilitated by the B1→tetragonal VN transformation, which is a martensitic transition caused by phonon instabilities of the cubic lattice below 250 K (see schematic representation in figure 5 of [64], and results of electron and x-ray diffraction in figures 6 and 7 of [39]). On the contrary, experimental results indicate that an anion vacancy concentration of 3% thermodynamically stabilizes the B1 structure of VNx at cryogenic temperatures [63].
We propose that the greater shear elastic stiffnesses of VN0.8 explain its experimentallymeasured higher hardness in relation to VN. In addition to being harder than VN, the understoichiometric compound exhibits the remarkable capability to adapt its mechanical response to the loading condition. As the shear strain increases beyond ≈9%, VN0.8 becomes progressively more compliant to change of shape. This facilitates plastic deformation in VN0.8, as evidenced by the smaller ideal shear strengths S and strains S required to activate {110}〈11 ̅ 0〉 and {111}〈11 ̅ 0〉 slip in VN0.8 vs. VN ( Fig. 6 and Table 1).
The results of AIMD simulations allow us understanding the excellent resistance to fracture of VN0.8 in comparison to the brittleness of VN (Fig. 3). The propagation of atomic-scale cracks in, and fracture toughness of, solid crystals depends in a non-trivial manner on the loading condition, crack and lattice geometries, and distribution of stresses within the material [65,66].
Generally speaking, however, crack propagation can be prevented (crack blunting) if the stress concentrated around the crack tip is rapidly dissipated by plastic flow [67]. Crack blunting may occur if the emission of dislocations at an angle to the fracture plane is sufficiently fast, that is, in relation to the speed of crack growth (bond-snapping rate) [67]. According to Schmid's law [ [70,71].
The fact that the tensile/shear strength ratios calculated for VN0.8 (≈0.7) are much larger than those (≈0.5) obtained for VN (see Table 3) is consistent with the superior fracture resistance of the understoichiometric compound (Fig. 3). We note, however, that our experimental results indicate that the toughness of VN0.8 is not uniquely the result of plastic deformation. In fact, lattice slip in VN0.8 occurs gradually, that is, to the extent necessary to prevent stress accumulation. This is demonstrated by absence of sudden pop-ins during experimental stress/strain measurements, Fig.   2, and by more moderate material pile-up subsequent to nanoindentation in comparison to VN, Fig.   3c. The fact that plastic flow (Fig. 6) and local structural transformations (Figs. 4 and 5b) in the understoichiometric compound become operative at elevated stress conditions is a requisite necessary for the simultaneous enhancement of hardness and toughness. The brittleness of VN is reflected by its lower tendency to activate lattice slip, as indicated by considerably larger shear strengths S than in VN0.8, and by the fact that stoichiometric phase undergoes sudden failure beyond its tensile yield points (Fig. 5a).
DFT electronic-structure calculations provide fundamental insights for the electronic mechanisms which enhance hardness and toughness in VN0.8. Energy-resolved electron densities, a technique for visualizing electronic states (see, e.g., Refs. [26,72]), is here employed to demonstrate the enhancement of d-d metallic interactions in the vicinity of anion vacancies. Fig. 7 illustrates the VN0.8 electron densities resolved in energy intervals that primarily correspond to d-t2g (Fig. 7a) or d-eg states (Fig. 7b)  in VNx has been previously suggested for lattices with isotropic vacancy ordering [73]. Due to high metallic nature, the electron density of VN0.8 rearranges considerably when the material is subjected to external loads comparable to its ideal tensile T or shear S strengths. As detailed below, strain-mediated electron transfer both leads to modifications in bonding geometries during elongation and assists lattice slip upon shearing, thus enabling rapid stress dissipation. Fig. 8 illustrates the electron-transfer maps calculated for VN0.8 under tensile deformation.
In Fig. 8a, five adjacent atomic layers are labeled in alphabetic order from "a" to "e". V and N atomic sites along each plane are numbered sequentially from "1" to "4". Note that positions "2d" and "4d" are vacant. In unstrained VN0. 8 and "c" V layers, which is evidenced by zig-zag patterned electron accumulation between all V-V pairs (Fig. 8a, 20% strain).  Fig. 8b) results in increased electron accumulation in the remaining (reinforced) vertical V-V and V-N bonds (note, e.g., Vi -Vii -Ni -Viii -Nii chain of bonded atoms in strained VN0.8). The strainmediated electronic mechanism leads to bond-bending and buckling of (001) layers, which confer on VN0.8 a great fracture resistance.
Our stress/strain results indicate that the initially relatively higher resistance to shearing of VN0.8, Fig. 6, originates from an overall stiffening of nearest-neighbor V-N bonds. The effect is due to a larger fraction of electrons that can be employed in d-eg-p bonding states, analogous to the vacancy-induced strengthening demonstrated for B1 VMoNx alloys (see figure 7a in [7]).
However, for xz > 0.09, the VN0.8 response to shearing changes from hard to compliant, which facilitates lattice slip. A representative example is illustrated in Fig. 9, where electron transfer maps are calculated for a sequence of VN0.8 atomic configurations sampled during {110}〈11 ̅ 0〉 slip. Time progression increases in alphabetic order from Fig. 9a to Fig. 9i. The (110) atomic layer labelled as "D" slides against the (110) layer "C". Lattice slip is assisted by continuous reorganization of d-electron clouds near to Nvac: back and forth transfer of electrons among  d-eg V2-V3 fourthneighbor bonds normal to the slip plane and  d-t2g V1-V3 second-neighbor bonds parallel to the glide direction. In particular, a shortening of the  d-eg V2-V3 bond lifts the V3 atom upward (Fig.   9c-g), thus favoring lateral [1 ̅ 10] translation of layer "D". Presumably, N diffusion also facilitates the crystal glide process, as suggested by the fact that, while the N1 and N3 atoms (see Fig. 9a-d) migrate out of the plane of view, two additional vacancies become visible upon completion of the slip process (compare Fig. 9a with Fig. 9i).

Summary and perspectives
We combine experiments and first-principles simulations to demonstrate that the control of anion vacancy concentrations in epitaxial single-crystal B1 VNx(001)/MgO(001) ceramics allows simultaneous enhancement of both material's hardness and toughness.  [48] UT is the energy density absorbed during deformation: area that underlies a stress/strain curve in Fig. 4.   Table 3. Calculated tensile-to-shear strength ratios, here used as indicators of fracture resistance.      For both {110}〈11 ̅ 0〉 and {111}〈11 ̅ 0〉 slip systems, VN0.8 exhibits an initially stronger resistance to shear deformation (up to ≈9%), followed by more facile activation (lower ideal shear strengths) of layer glide than calculated for the stoichiometric compound.