Hyperbolic metamaterials with extreme mechanical hardness
Arrigo Calzolari, Alessandra Catellani, Marco Buongiorno Nardelli, Marco Fornari
HHyperbolic metamaterials with extreme mechanical hardness
Arrigo Calzolari,*, ∗ Alessandra Catellani, Marco Buongiorno Nardelli,, and Marco Fornari CNR-NANO Research Center S3, Via Campi 213/a, 41125 Modena, Italy Department of Physics, University of North Texas, Denton, TX 76203, USACenter for Autonomous Materials Design,Duke University, Durham, NC 27708, USA Department of Physics and Science of Advanced Materials Program,Central Michigan University, Mt. Pleasant, MI 48859 USACenter for Autonomous Materials Design,Duke University, Durham, NC 27708, USA (Dated: February 23, 2021)
Abstract
Hyperbolic metamaterials (HMMs) are highly anisotropic optical materials that behave as metals or asdielectrics depending on the direction of propagation of light. They are becoming essential for a plethoraof applications, ranging from aerospace to automotive, from wireless to medical and IoT. These applicationsoften work in harsh environments or may sustain remarkable external stresses. This calls for materials thatshow enhanced optical properties as well as tailorable mechanical properties. Depending on their specificuse, both hard and ultrasoft materials could be required, although the combination with optical hyperbolicresponse is rarely addressed. Here, we demonstrate the possibility to combine optical hyperbolicity and tunablemechanical properties in the same (meta)material, focusing on the case of extreme mechanical hardness.Using high-throughput calculations from first principles and effective medium theory, we explored a largeclass of layered materials with hyperbolic optical activity in the near-IR and visible range, and we identified areduced number of ultrasoft and hard HMMs among more than 1800 combinations of transition metal rocksaltcrystals. Once validated by the experiments, this new class of metamaterials may foster previously unexploredoptical/mechanical applications. a r X i v : . [ phy s i c s . op ti c s ] F e b NTRODUCTION
The ability to engineer new materials by integrating different components at the nanoscale hasrevolutionized science and technology in the last 100 years. The transistor [1], the giant magneto-resistance [2], and the pioneering work of Veselago and Pendry [3, 4] in the area of negative refractionhave shown unexpected physics and opened the doors to new technologies [5, 6]. Hyperbolic opti-cal metamaterials (HMMs) fall within this scientific endeavor. They are highly anisotropic opticalmaterials which behave as metals or as dielectrics depending on the direction of propagation of theelectromagnetic field. Their properties are reflected in the sign of the diagonal terms of the dielectrictensor that is controlled during deposition by engineering the composition of the multilayer systemin the subwavelength regime [7]. Due to their hyperboloidal isofrequency surface, HMMs have beenproposed for a broad range of applications working from THz to UV-visible frequencies, includingnegative refraction [8, 9], optical cavities [10], biosensing [11], and waveguides [12]. The hyperbolicbehavior of materials is particularly relevant in the visible range for application such as superlenses[13], cloaking [14], sub-wavelength imaging [15, 16], and perfect-absorbers for IoT [17]. In addition,the potential for highly directional propagation of electromagnetic modes localized at subwavelengthscales, e.g. volume plasmon polaritons [18, 19], opens a new route into several nanophotonic applica-tions as fluorescence engineering [20], super-Planckian thermal emission [21], subsurface sensing [22],cryptography [23], and Dyakonov plasmons [24].The outstanding optical value of HMMs for applications cannot be disentangled by their mechanicalproperties, especially when the devices are expected to perform at extreme conditions such as hightemperature, high pressure, under large tensile strain, and so on. For example, highly efficient antennasand radars are necessary for improved wireless communication, space vehicle navigation, and GPSsatellites. Antenna technology has gained great advantages from the integration with HMMs materials[25] that enhance the dipole emission to free space, and provide a broadband impedance matching[26, 27]. However, in space conditions these devices need to withstand uniquely harsh environments:strength, hardness and stiffness become unavoidable properties to resist to the forces and the collisionsthese systems are exposed to. Similar arguments hold for HMM-based devices for microwave [28] oroptical [14] cloaking exploitable in naval or military applications as well as for radar scanning systemsfor automotive [29]. For all these applications the possibility to use hard (or ultrahard ) materials couldimprove the mechanical resistance of the optical device. On the opposite side, ultrasoft metamaterialscould be of great interest for applications, e.g., in the medical industry [30], where the ability tomanipulate electromagnetic waves and provide diagnostic images can be coupled to the extremelyflexibility of soft matter and the ability to encapsulate and transport molecular systems (e.g., drugs)in biological systems. On a general ground, the combination of hyperbolic optical properties andtailorable mechanical properties would allow to optimize the characteristics (e.g. endurance) of theoptical devices, opening to previously unexplored solutions. In order to achieve this goal, the conceptof optical metamaterials has to be extended to the mechanical realm [31] by considering the effect ofthe multilayer geometry on the elastic properties and mechanical hardness.Unfortunately, both ultrasoft and hard materials are rare and the ones known so far such as diamond,cubic BN, metal borides and carbides (WB , ReB , WC, TiC), and ceramics (Si N , Al O ), do notexhibit the special optical properties of HMMs [32]. This calls for new classes of artificial materials.By synergistically optimizing the optical properties and the mechanical responses, we have designedHMMs with extreme ( ultrasoft and hard ) mechanical hardness starting from simple transition-metalrocksalt crystals, such as TiN, ZrN, which are practical for the realization of hyperbolic superstructures[33]. Our results are based on high-throughput first principles and effective medium theory calculationsof the optical and mechanical properties of periodic superlattices and provide an efficient strategy to fineengineering the optical and mechanical properties. We demonstrate the effectiveness of the approach2 a)(b) e || e ⊥ h h Sc Ti V Cr Mn Fe Co Ni Cu ZnY Zr Nb Mo Tc Ru Rh Pd Ag CdLa Hf Ta W Re Os Ir Pt Au Hg
B C N U [ ] [ ] [ ] z Figure 1. a) Scheme of the superlattice models used to design the metamaterials. Red arrows indicate thecomponents of the dielectric function along the direction parallel ( (cid:15) (cid:107) ) and perpendicular ( (cid:15) ⊥ ) to the opticalaxis z . High symmetry directions of cubic lattice are aligned to cartesian axes. b) Table of TMX rocksaltcompounds obtained combining the 30 transition metals (TM) with B, C, and N (X). Gray sections indicatethat the material is mechanically unstable (U). by identifying several ultrasoft and hard HMMs among more than 1800 rocksalt combinations. DESCRIPTORS AND SELECTION CRITERIA
Single TMX crystals.
We considered superlattices composed of the alternating layers of two rocksaltTMX materials along the cubic [001] direction (Figure 1a). Ninety TMX compounds were obtainedusing all possible transition metals (TM) with the non-metal elements X=(B, C, N). All structural,electronic, optical and elastic properties of single TMXs were evaluated from first principles, by usingapproaches based on density functional theory (DFT). The optimized lattice parameter ( a ), theformation energy (∆ U ) with respect to the stable phases of the single components, and the complexdielectric function (ˆ (cid:15) = (cid:15) r + i(cid:15) i ) were calculated with the Quantum Espresso code [34]. The ElaSticcode [35] has been used to obtain the elastic stiffness coefficients ( c ( (cid:96) ) ij ) of the (cid:96) -th TMX crystal; seeSec. 5 and Sec. S1 of Supporting Information (SI) for the full description of the method.About one third of the systems (28/90) resulted to be thermodynamically and mechanically unstablehaving both formation energies ∆ U > c <
0. MoN has negative formationenergy but it has been discarded for its mechanical instability ( c <
0) in the cubic phase, in agreementwith previous theoretical calculations [36]. The remaining 62 stable compounds formed the building3locks for the superlattices. A summary of the rocksalt compounds with negative formation energycan be found in Figure 1(b).Carbide and nitride compounds of the groups 4 and 5 (e.g. TiC, TiN, VN, etc.) are stable andvery well characterized from both experimental [37, 38] and theoretical [39–42] perspectives. Theseare hard refractory materials exploited in a large range of applications from mechanics to electronicsand optics. The carbide and nitride compounds with the remaining transition metals are more rare(e.g. RuC, AuC) and partially unexplored (e.g. FeC, MnN). Monoboride compounds are also quiterare, being the diboride structures the most common ones. Most monoborides favor the orthorhombicor tetragonal phases [43], although the metastable NaCl structures have been predicted [44].All bulk components are metallic, except ScN, YN, and LaN that exhibit small indirect bandgaps( E g < . eV ). The analysis of the main electronic and optical properties of the stable TMX crystals isreported in Sec. S2 of SI. The stiffness coefficients, the elastic moduli, and the mechanical parametersof all 62 crystals are summarized in Sec. S3 of SI, in very good agreement with previous results (seeSI for references). TMX superlattices.
By combining all possible pairs of stable TMX rocksalts, we constructed 1891binary superlattices composed of alternating layers of different materials along the [001] direction(Figure 1a). The fractional volume of the (cid:96) th layer is defined as f (cid:96) = h (cid:96) /h ( (cid:96) = 1 , h = h + h isthe volume of the superlattice unit. The optical and elastic properties of the superlattices were obtainedwithin the effective medium theory (EMT) [45]. Experimentally, the thickness of the individual layersranges from 5 to 200 nm. It is worth noticing, that although the effective medium theory used inthis work depends only on the volume ratio and not on the physical thickness of the single layers,confinement and interface effects in ultrathin layers could have important consequences. However,previous experimental [33] and theoretical results [46] demonstrated that EMT holds in the case ofmetallic layers (e.g. TiN) down to the shortest 5 nm thickness. We have tested the validity of theEMT for selected configurations, as discussed in Sec. S4 of SI.HMMs are anisotropic uniaxial materials that derive their name from the shape of the isofrequencysurface ω ( k ) where k , and ω are the wavevector and the frequency of the radiation. In non-magneticuniaxial materials, the dielectric tensor can be written in terms of two diagonal parameters: (cid:15) (cid:107) and (cid:15) ⊥ , where the subscripts (cid:107) and ⊥ indicate components parallel and perpendicular to the anisotropyaxis (Figure 1). If one of the two parameters is negative, the isofrequency surface opens up into ahyperboloid. The choice (cid:15) ⊥ > (cid:15) (cid:107) < (cid:15) ⊥ < (cid:15) (cid:107) > (cid:15) is given by˜ (cid:15) (cid:107) = (cid:15) (cid:15) f (cid:15) + f (cid:15) (1)˜ (cid:15) ⊥ = f (cid:15) + f (cid:15) , where (cid:15) and (cid:15) are the complex dielectric functions of the two constituents, respectively, and f (cid:96) arethe volume fraction of the (cid:96) th layer within the superlattice.The elastic properties of superlattices are also derived from the bulk constituents along the linesproposed by Grimsditch and Nizzoli [48]. This approach allows to determine the total stress (¯ σ ) andstrain (¯ η ) tensor of the superlattice in terms of the stress and strain of the two individual layers:¯ σ = f σ + f σ = ¯ C ¯ η, (2)where ¯ C = { ¯ c ij } is the effective elastic matrix relative to the whole superlattice (see SI, Sec. S4,for the complete theoretical description). This approach has been profitably used to study severalsuperlattices and multilayer compounds, such as TiN/WN [49], GaN/AlN [50], Si/Ge [51], and Ni/Mo452]. When all the effective stiffness (¯ c ij ) and compliance (¯ s ij ) parameters are known, all the derivedquantities, such as the elastic moduli, can be obtained by substituting the effective parameters inthe original formula. The effective bulk modulus ( B ) and shear modulus ( G ) for the superlatticescan be evaluated within the Voigt–Reuss–Hill approximation [53], see Sec. S1 of SI for the completeformulation. In a similar way, the Young modulus ( E y ) and the Poisson ratio ( ν ) as well as othermechanical properties of superlattices can be straightforwardly derived from the effective moduli inthe Hill notation B = B H and G = G H .For example, the development of plastic deformations in crystals or the determination of micro-cracks in materials are related to the elastic anisotropy [53], which affects the mobility of dislocations.Here, we adopted the universal anisotropic index A U proposed by Ranganathan and Ostoja-Starzewski[54], defined as: A U = 5 G V G R + B V B R − ≥ , (3)where the subscripts V and R identify the Voigt and Reuss expressions, respectively (Sec. S1, SI). A U is zero for isotropic materials, while a deviation from zero indicates anisotropic mechanical properties.Another relevant parameter is the Pugh modulus G/B , which is involved in the strain of fractures andin the ductility of solids. The competition between plastic flow and brittle fracture can be quantifiedin terms of the solidity index S [55], defined as S = GB and running from 0 and 1. The index is zero fora liquid and reaches its maximum value for diamond. S = 0 .
23 is assumed as dividing point, betweenductile (
S < .
23) and brittle (
S > .
23) materials, which means that materials with lower solidityindex show higher plasticity. We used the Vickers model to estimate the hardness of materials [56].The Vickers scale measures the indentation resistance which results from the combination of resistanceto plastic flow, phase transformation, and fracturing. The
Vickers hardness H V is estimated by usingthe Tian empirical approach [57]: H V = 0 . (cid:16) GB (cid:17) . G . . (4)The comparison with alternative formulations of H V is reported in Table S4 of SI. Materials can beroughly classified in term of their hardness scale [58]. Superhard materials are defined as those withVickers hardness greater than 40 GPa and ultrahard materials as having a hardness exceeding 80 GPa.Very low Vickers hardness ( H V < . ultrasoft materials .A material is hard if it resists to indentation, i.e. to plastic deformations. This involves thecapability to contrast the nucleation and the motion of dislocations that may generate fractures.Thus, in the indentation hardness tests, one particular role is played by the fracture toughness ( K IC ),which measures the resistance of a material against crack propagation [59]. Fracture toughness iscritical for brittle materials, where the easy formation of cracks may cause to the failure of the sample. K IC can be estimated combining the bulk and shear moduli as [60]: K IC = V / G (cid:16) BG (cid:17) / , (5)where V is the volume per atom. Notably, mechanical properties such as hardness, and fracturetoughness are qualitative rather than quantitative characteristics and they depend not only on theproperties of a material, but also on the measurement method and interpretation of the results [61].We initially considered superlattices with 50%-50% of constituent composition (i.e. f = f = 0 . a of 4%. Thisthreshold allows to include systems such as TiN/NbN (∆ a = 3 . a = 3 . a = 1 . a = 4 . a = 6 . E ∈ [0 . − . hyperbolic if the hyperbolic condition (cid:15) ⊥ · (cid:15) (cid:107) < H V < . H V > . RESULTS AND DISCUSSION
We focus on the 17 superlattices that simultaneously fulfill the hyperbolicity and hardness conditionsdescribed above. A few of them, such as HfN/ScN [71, 72], ZrN/ScN [73], TiC/TiN [74, 75], andHfC/HfN [76], have been realized and experimentally characterized, although not studied for hyperbolicoptical applications.The list of the resulting compounds and their main elastic properties are summarized in Table I.Only two out of the seventeen systems are ultrasoft materials and both include RhB as a constituent.The remaining systems are hard materials each including at least one nitride component (mostly ScN,HfN). More generally, the most recurrent TM elements are from groups 4 and 5, such as Zr, Nb, Hf,Ta. TMX compounds from groups 7-12 (except ZnN) do not match the selection criteria: most of themhave pure metallic character in the energy range under investigation. All the selected superlattices donot exhibit magnetic ordering. 6 able I. Elastic properties of hyperbolic superlattices. Elastic moduli (
B, G, E y ), Young parameters (E iy ) areexpressed in GPa. Percent lattice mismatch (∆ a ), Poisson ratio ( ν ), universal elastic anisotropy (A U ) areadimensional parameters.System ∆ a B G ν E y E [100] y E [001] y E [110] y A U RhB/RuC 1.3 310.0 50.1 0.42 142.6 400.9 389.5 83.3 8.20RhB/ZnN 2.2 231.5 46.3 0.41 130.3 280.3 276.2 144.8 3.84MoB/HfN 1.4 270.9 160.7 0.25 402.4 560.8 560.7 355.6 0.60NbB/HfN 1.9 249.1 154.1 0.24 383.2 532.8 532.1 338.4 0.60NbB/ZrN 0.2 241.2 152.3 0.24 377.4 506.1 506.3 339.2 0.47TaB/HfN 1.6 257.6 160.1 0.24 397.9 553.9 554.0 351.3 0.61WC/ScN 2.7 273.2 167.2 0.25 416.7 520.0 473.6 404.7 0.19TaB/ZrN 0.2 249.4 158.3 0.24 392.1 527.2 526.4 352.1 0.48NbC/HfN 0.6 281.1 175.4 0.24 435.5 574.2 573.4 402.6 0.41NbC/ZrN 2.4 272.3 173.7 0.24 429.6 547.4 543.4 402.4 0.31HfC/HfN 2.7 248.9 164.7 0.23 404.7 502.7 498.6 383.6 0.25HfN/ScN 0.1 234.2 159.2 0.22 389.3 453.0 436.0 382.6 0.10ScN/WB 0.8 248.0 166.3 0.23 407.7 473.5 448.4 399.4 0.09ZrN/ScN 1.9 226.8 157.3 0.22 383.2 426.0 416.5 379.8 0.05TiC/TiN 2.0 259.4 180.5 0.22 439.6 489.6 486.5 426.6 0.06NbC/ScN 0.5 248.6 183.0 0.20 440.8 474.1 448.6 437.0 0.02TaC/ScN 0.8 259.0 193.5 0.20 464.8 524.3 478.3 457.1 0.05Figure 2. Elastic properties of hyperbolic superlattices: (a) bulk B and shear G moduli, (b) universal elasticanisotropy A U . Full (empty) symbols correspond to superlattice (single TMX) systems. Inset in panel (b)zooms on the lower A U values. The effective bulk and shear moduli of the hyperbolic superlattices are shown in Figure 2 (filledsymbols), along with the corresponding values of the single rocksalt components (empty symbols).7 igure 3. 2D plots of angular Young moduli for selected hyperbolic superlattices, projected on the (a) (001),and (b) (1¯10) planes.
Except for RhB/RuC system ( B = 310 GPa), the bulk modulus of all other compounds is in the range B ∈ [200 − G <
100 GPa) and hard material ( G ∼
150 GPa). This confirmsthe predominant role of shear coefficient in the definitions of hardness (Eq. 4). Notably, the valueresulting for the superlattice is always included between the values of the corresponding componentmaterials.Figure 2b shows the elastic anisotropy trend: the larger is the deviation from zero, the stronger isthe anisotropy. A U index is larger than 3.8 for RhB/RuC and RhB/ZnN; while A U is less than 0.6for the remaining systems. This indicates that ultrasoft materials are more anisotropic than the hardones. The directional dependence of Young modulus is used to characterize the elastic anisotropy alongcrystallographic direction (Sec. S1, SI), as it describes the ability of a material to resist to a directionaldeformation within the range applicable in linear elasticity. Figure 3 shows the 2D Young plots forselected hyperbolic superlattices, projected on the (001) (panel a), and (1¯10) plane (panel b). SystemsRhB/RuC and RhB/ZnN have a lobate (flower-like) and anisotropic behavior, especially along the[001] stacking direction. The other systems have more symmetrical character, being almost circular(i.e. isotropic) on the (001) plane (Figure 3a), rather all curves are prolate along the growth direction[001]. This points to the indirect interplay between elastic anisotropy and mechanical properties [77],through the formation of dislocations in the material [78]. The presence of dislocations increases theelastic anisotropy and reduces the elastic moduli of the material [79].The mechanical properties of hyperbolic materials summarized in Table II. Figure 4a shows thesolidity index S (i.e. the Pugh modulus) that is used to discriminate the brittle vs ductile behaviorof superlattices (full symbols), with respect to their rocksalt components (empty symbols). Again,we observe a net difference between systems RhB/RuC and RhB/ZnN which are ductile ( S < . S coefficient is also a measure of the inhomogeneity of the electron density (the so-called covalence ) and8 able II. Mechanical and hyperbolic optical properties of hyperbolic superlattices. Vickers hardness (H V ) areexpressed in GPa, fracture toughness ( K IC ) is expressed in MPa · m / . Solidity index ( S ) is adimensional.System S H V K IC typeRhB/RuC 0.12 1.9 1.8 IIRhB/ZnN 0.15 2.2 1.5 IMoB/HfN 0.44 18.5 3.1 INbB/HfN 0.46 18.9 3.0 INbB/ZrN 0.47 19.1 2.9 ITaB/HfN 0.47 19.5 3.1 IWC/ScN 0.46 19.7 3.2 IITaB/ZrN 0.48 19.8 3.0 INbC/HfN 0.47 20.9 3.3 INbC/ZrN 0.48 21.3 3.3 IHfC/HfN 0.50 21.3 3.1 IHfN/ScN 0.51 21.5 2.9 I+IIWB/ScN 0.50 21.8 3.0 IIZrN/ScN 0.52 21.8 2.8 I+IITiC/TiN 0.52 24.1 3.2 INbC/ScN 0.55 26.0 3.2 IITaC/ScN 0.56 27.5 3.4 IIFigure 4. Mechanical properties of hyperbolic superlattices: (a) solidity index S , (b) fracture toughness K IC ,(c) Vickers hardness H V . Full (empty) symbols correspond to superlattice (single TMX) systems. is used to discriminate the interatomic bonding character of a crystal. The Pugh modulus increases inparallel with the directionality of the bonds: for metallic bonding we obtain S < .
3, while S ≈ . S > . S ≈ .
15) are typical of metallic systems, while the values of the othersamples ( S ≈ .
5) indicate an ionic bonding character.Since indentation resistance in hardness measurements is related to both the resistance to compres-sion (i.e. bulk modulus) and to the resistance to shape changes (i.e. shear modulus), it is not surprisingthat the Pugh modulus correlates also with both fracture toughness and the Vickers hardness, as shown9n panels 4b and c, respectively. Pugh modulus is negative correlation with fracture toughness: duringdeformation bonds break and reform resulting in displacement of atoms and slipping of atomic planes.Materials with low fracture toughness usually exhibit high ductility and yield critical strains. On thecontrary, high Pugh modulus, which results from directional bonds, increases the shear modulus andlimits the motion of dislocations that may generate fractures, thereby increases the material hardness.Thus, RhB/RuC and RhB/ZnN are simultaneously ductile, ultrasoft, and have the lowest fracturetoughness, while all other systems are brittle, hard, with higher fracture toughness.The calculated H V values for superlattices are in agreement with the reduced set of available ex-perimental data. For example, the measured microhardness for the TiC/TiN superlattice, H V =30.5-31.5 GPa [74], well matches the theoretical value, H V =24.1 GPa, obtained with the present method.Again, the calculated hardness of the superlattices is intermediate with respect to the single compo-nents (Figure 4c). This means from one hand that the formation of multilayers does not waste theaverage hardness of the components, as it may happen in some alloyed mixtures [81]. On the otherhand, layering cannot be assumed as a general strategy to increase the strength of hard materialsto obtain super- or ultra-hard ones. Indeed, we do not observe the increase of hardness ( hardening effect) reported for similar systems such as TiN/NbN superlattices, for which a large range and higherhardness values have been measured: H V = [17 . − .
0] GPa [82, 83]. The calculated hardness for thesame system is H V =16.2 GPa (see Table S3, SI) and close to the lower bound experimental limit. Theorigin of such large variability and of the hardening process may depend on several aspects [69, 84, 85],such as crystalline micro-granurality, grain boundaries, interface widths, cohesive forces and stressat the interfaces. A leading role in this process is also played by dislocations and their dynamics.Many mechanisms have been proposed to explain the dislocation effects on hardening of superlattices,including the barrier to dislocation motion provided by layers with differing shear moduli [86]; theinhibition of dislocation motion due to layer coherency strains [87], or the misfit dislocation arrays atthe interfaces [88]. This analysis goes beyond the aim of the present work and the calculated valuescan be safely assumed as lower limits of superlattice hardness.The hyperbolic character of the superlattices and their specific operational window is shown inFigure 5a, where we focus in the energy range E ∈ [0 . − .
5] eV, which includes near-IR, visible,and near-UV range that are the most interesting for optical applications. RhB/RuC and RhB/ZnNsuperlattices are hyperbolic in the visible-UV part of the spectrum, all other systems exhibit a type-Icharacter at higher energies or a type-II character at lower energies.We investigate in more details the optical properties of HfN/ScN, that we assume as the referencefor this class of TMX multilayers. HfN/ScN has been experimentally grown [71, 72] and behavesas both type-I and type-II material, depending on the energy of the incoming radiation. Figure 5bshows the real and imaginary part of both parallel ˜ (cid:15) (cid:107) and perpendicular ˜ (cid:15) ⊥ components of the effectivedielectric function (Eq. 2); see Fig. S2 in SI for the details of the constitutive dielectric functions.We define E I and E II , the energies at which the system changes its optical behavior from regular(e.g. metal) to hyperbolic type-I ( E I ) or type-II ( E II ), respectively. The parallel component (bluelines) has the behavior typical of semiconductors: the real part (straight line) is always positive for E < E II = 1 .
89 eV and reaches the dielectric constant value ( (cid:15) = 21.8) in the limit for E → E = 1 .
76 eV, which mainly derives from interbandoptical transitions in the ScN part. The perpendicular component (red lines) has, instead, a metalliccharacter: for
E < E I = 2 .
14 eV the real part is negative and diverges for E →
0, with a typicalDrude-like behavior that mostly comes from the intraband transitions of HfN component. Thus, for
E < E II Re [˜ (cid:15) (cid:107) ] is positive and Re [˜ (cid:15) ⊥ ] is negative, which corresponds to a type-II system. The oppositehappens for E > E I , where Re [˜ (cid:15) (cid:107) ] < Re [˜ (cid:15) ⊥ ] >
0, and the system has a type-I character. In theintermediate range between E II and E I , both components are negative and the system is metallic.In order to quantify the hyperbolic character of the multilayers, we define two figures of merit10 a) (c)(d)(b) Figure 5. (a) Energy distribution of hyperbolic superlattices. Gray shaded area indicates the visible range.(b-d) Hyperbolic optical properties of HfN/ScN superlattice: (b) real and imaginary part of the effectivedielectric functions ˜ (cid:15) (cid:107) and ˜ (cid:15) ⊥ ; (c) parallel and perpendicular quality factors Q j ; and (d) strength of dielectricanisotropy ∆˜ (cid:15) . Shaded areas in panel (b-d) indicate the energy regions with type-II ( E < E II , red area) andtype-I ( E > E I , blue area) hyperbolic behavior. Vertical dotted lines in panel (b) mark selected energies E i ,discussed in the angular analysis of the dielectric function (cid:15) ϕ (see Figure 6). [89, 90], namely the quality factor Q j = − Re [˜ (cid:15) j ] Im [˜ (cid:15) j ] , with j = (cid:107) , ⊥ ; and the strength of dielectric anisotropy∆˜ (cid:15) = Re [˜ (cid:15) (cid:107) − ˜ (cid:15) ⊥ ]. The quality factor Q accounts for the energy losses only in the direction in which thematerial has a metallic character. Systems with the quality factor Q ≈ E ∼ Q ⊥ =2.8, this qualifies HfN/ScN as a hyperbolic material with reduced energy loss in thenear-IR range. At higher energies ( E > . (cid:15) factors quantifies the metal-dielectric anisotropy along the direction paralleland perpendicular to the optical axis. Systems with ∆˜ (cid:15) >
20 are considered good hyperbolic materials[91, 92]. In the case of HfN/ScN, this condition is fulfilled almost in the entire range
E < E II , wherethe system has a type-II behavior.The energy range, the hyperbolic type, and the gain factors strictly depend on the composition(i.e. the crossover energies, and optical transitions in each TMX layer) and on the thickness ratio(e.g. f i ) of the multilayer, as shown in Figure S8 of SI, where the gain factors of selected superlatticesare shown. Ultrasoft systems RhB/RuC and RhB/ZnN have low quality factors ( Q < .
5) and low11ielectric anisotropy (∆˜ (cid:15) <
20) for the largest part of the energy spectrum under consideration. Wecan conclude that RhB/RuC and RhB/ZnN conjugate rare qualities such as ultrasoft hardness andhyperbolic light dispersion, but are plagued by high energy dissipation because of interband transitionsfrom non-metal sp to transition-metal d bands. Hard superlattices HfN/ScN, ZrN/ScN and TaC/ScNhave instead the highest gain factors that correspond to smaller energy loss. Depending on the specificneed of an application - such as the energy working range or the dissipation tolerance - it is possible toselect the optimal TMX pair and to fine tune the chemical and structural details of the superlattice.The strength of dielectric anisotropy also affects the dispersion relation of the radiation that maypropagate along the medium, as (cid:15) (cid:107) and (cid:15) ⊥ define the geometric characteristics of the hyperbolicisosurface equation: k + k (cid:15) (cid:107) + k (cid:15) ⊥ = k ⊥ (cid:15) (cid:107) + k (cid:107) (cid:15) ⊥ = ω c . (6) Figure 6. 2D k-dispersion plot for dielectric (air and ScN) and hyperbolic (HfN/ScN) materials at (a) E = 0 . E = 2 .
50 eV. Labels (R) and (C) indicate if only real part or the complete complex dielectricfunction is considered in Eq.6.
The k-dispersions for HfN/ScN case are shown in Figure 6, along with the plots for representativedielectric systems (namely air and ScN), included for comparison. We considered two radiation energies E = 0 .
83 eV (panel a) and E = 2 .
50 eV (panel b) that are representative for type-II and type-I character of HfN/ScN superlattice. The k wavevector follows a truly hyperbolic surface (i.e. openform) only for the ideal case with zero absorption processes, i.e. when only the real part of the dielectricfunction is considered in Eq. 6. In the case of absorbing materials, the isosurface equation assumes aform whose final dispersion crucially depends on the imaginary part of the dielectric function. In thecase of air, the dielectric function reduces in all direction to the dielectric constant (cid:15) r = 1 .
0, and thek-dispersion relation is a circumference for both energies. The case of ScN at E = 0 .
83 eV is similar:the imaginary part of the dielectric function is almost zero. At higher energy ( E = 2 .
50 eV) ScN isoptically active and the imaginary part plays a crucial role ( (cid:15) r = 11 . (cid:15) i = 9 . E = 0 .
83 eV (type-IImaterial) the inclusion of the imaginary part has minor effect on the k-dispersion. In this case thehyperboloid isosurfaces reduce to a hyperbola (Figure 6a). At E = 2 .
50 eV (type-I material) thehyperbola is open if only the real part of the dielectric function is considered. The inclusion of theimaginary part (i.e. dissipation) gives rise to a close loop (panel b).12ince electromagnetic radiations propagate through dielectric systems ( (cid:15) r > (cid:15) r < (cid:15) ⊥ · (cid:15) (cid:107) < volume plasmon polaritons (VPPs) [19, 92]. In anisotropic materials, the electric-field vector E andthe electric-displacement vector D are not usually parallel. As a consequence, the Poynting vector S and the direction of the wavefront k are no more parallel (Figure 7a). (a) z Sk 𝜑 ! 𝜗 ! (b)(c) 𝜑 ! 𝜗 ! Figure 7. (a) Vector diagram for transverse magnetic propagating waves (TM) and hyperbolic dispersionisosurface corresponding to a type-II HMM; (b) Real part of the effective function (cid:15) ϕ and (d) ϑ ( ϕ ) anglebetween the extraordinary wave and the optical axis, at different energies ( E i ) of the incoming electric field,for the HfN/ScN superlattice. Energies ( E i ) are defined in Figure 5b. Vertical dashed lines identify the criticalangles ϕ c and ϑ c , respectively. The coupling between the Poynting vector and the wavevector can be investigated through theangular description of the dielectric function [92]:1 (cid:15) ϕ = sin ϕ ˜ (cid:15) (cid:107) + cos ϕ ˜ (cid:15) ⊥ , (7)where ϕ is the angle between the wavevector of the radiation k and the optical axis z (Figure 7a).Since ˜ (cid:15) (cid:107) and ˜ (cid:15) ⊥ are functions of the energy of the incoming radiation, (cid:15) ϕ is also a function of theenergy. Figure 7b shows the real part of (cid:15) ϕ of the HfN/ScN system, for a set of selected energies, E i ,marked as vertical lines in Figure 5b. E I and E II are the energies that delimit the type-I and type-IIranges, respectively and are both in the visible range. E = 0 .
83 eV ( λ = 1 . µ m) is characteristic13nergy (wavelength) for telecommunications, E = 1 .
29 eV is a representative energy in the near IR, E = 1 .
66 eV and E = 1 .
73 eV are selected energies in the visible range. For all energies, except E I , Re [ (cid:15) ϕ ] has the same general trend, being negative at low angles and positive for ϕ → π/
2, theopposite holds for E I . For ϕ = 0 ( ϕ = π/ (cid:15) ϕ = ˜ (cid:15) ⊥ ( (cid:15) ϕ = ˜ (cid:15) (cid:107) ) that is negative (positive) for a type-II(type-I) HMM. As the energy increases from E to E I , Re [ (cid:15) ϕ ] we observe a reduction of the intensityespecially of the negative part, along with a shift of the maxima/minima towards higher angles.The condition Re [ (cid:15) ϕ ] = 0 determines the angular boundary between the metallic and dielectricresponse of the metamaterial, which is the analogous of plasmon excitation in planar interfaces. Thecritical angle ϕ c is the direction of radiation at which there can be the excitation of a travellingplasmon-polariton wave in the system. If we consider the interface with another dielectric medium,the wavevector at the interface follows the Snell’s law. Thus, the angle θ between the extraordinarywave (i.e. the Poynting vector) and the optical axis is: ϑ ( ϕ ) = arctan (cid:18) Re [ ˜ (cid:15) (cid:107) ˜ (cid:15) ⊥ ] tan ( ϕ ) (cid:19) . (8)At the critical angle ϕ c , the quantity θ c represents the angle of propagation of the VPP along themetamaterial [94]. This means that VPP propagates throughout the volume of the HMM along a cone with axis coincident to the optical axis and aperture θ c (Figure 7a). The trend for ϕ c follows the oneobserved for Re [ (cid:15) ϕ ], as shown in Figure 7b. The energy dependent ϑ ( ϕ ) and the critical angles ϑ c areshown in panel (c). We observe a change in the concavity of the ϑ ( ϕ ) as a function of the incomingradiation energy (e.g. E and E ). For E ≤ E , ϑ c is negative implying an anomalous refraction for theextraordinary wave. The refraction is anomalous in the sense that the transmitted wave is a backwardwave with the Poynting vector and the wavevector that are in opposite lateral directions. This is nota case of double negative refraction as originally formalized for left-handed materials by Veselago [3],i.e. negative dielectric permittivity (cid:15) < µ <
0. Rather, it hasbeen shown that it possible to obtain backward waves with uniaxial anisotropic media, under certainconditions and when only one parameter has a negative value [93]. In the present case, µ (cid:107) , µ ⊥ , (cid:15) (cid:107) > (cid:15) ⊥ <
0. This allows for backward waves only for transverse magnetic mode radiation [93]. For
E > E , ϑ is positive, even though (as for E and E II ) the system as the same type-II character. Thisis due to the ˜ (cid:15) (cid:107) / ˜ (cid:15) ⊥ fraction in Eq. 8. For E > E I both ϕ c and ϑ c are quite zero and positive, and (cid:15) ϕ (cid:39) ˜ (cid:15) ⊥ is small and positive (Figure 5b). The excited transverse magnetic mode waves can thus travelalmost parallel to the optical axis as through a low dielectric medium. This analysis opens up thepossibility to use hard materials also for negative refraction applications and cloaking in the visiblerange.Finally, we investigated the effect of the filling factor f (cid:96) in the results presented above. For theseventeen selected superlattices we calculated the mechanical properties and the hyperbolic behaviorfor f ranging from 0 to 1 ( f = 1 − f ). In the case, e.g. of RhB/RuC pair, f = 0 corresponds to theRuC single crystal and f = 1 corresponds to the RhB one, the same for all the other systems. Theresults for the calculated Vickers hardness are shown in Figure 8a. The volume ratio f (cid:96) enters in thedefinition of both the effective stiffness coefficients of the superlattices and the effective parallel andperpendicular dielectric functions. Thus, the selective criteria for the hardness or hyperbolic behaviordo not necessarily hold for all the f (cid:96) values. For example, in the case of RhB/ZnN system, all thecriteria are fulfilled only in the range f ∈ [0 . − . a) (b) Figure 8. (a) Hardness of hyperbolic superlattices as a function of the filling factor f . (b) Energy distributionof hyperbolic response for HfN/ScN, as a function of the filling factor f . the energy distribution of the hyperbolic response for HfN/ScN, as a function of the filling factor f .The largest energy range for hyperbolic optical responses is close to f = 0 .
5, while both type-I andtype-II ranges reduce as the layer composition moves from 50-50%. In other cases, as for RhB/RuC,NbB/HfN, NbB/ZrN, TaB/HfN and TaB/ZrN, the hyperbolic character changes with the volumecomposition: RhB/RuC is type-I for f < .
45 and type-II for f ≥ .
45; NbB/HfN is type-I in thevis-UV region for the entire range f ∈ [0 . − .
60] and it is also type-II in the near-IR part of thespectrum for f ∈ [0 . − . f (cid:96) = 0 .
5, fulfill the query for different volume ratio, for example at f = 0 .
75 AgC/ScN and CuC/RuCsuperlattices are ultrasoft materials ( H V =2.0, and 2.3, GPa) with type I+II and type-II hyperboliccharacter, respectively. We conclude that, within the limits of the effective medium approach, thevolume ratio is a critical degree of freedom that can be tailored to reach the specific mechanical andoptical properties needed for a specific application. CONCLUSION
The realization of HMM relies on the fabrication of multilayer structures where the choice of thecomponent materials provides the strategy to fine engineering the optical and mechanical properties.We used high-throughput techniques based on density functional and effective medium theory todesign HMMs with selected mechanical properties, by combining the properties of simple TMX rocksaltcrystals. The simultaneous requirement of structural lattice match, extraordinary mechanical hardness(both solid-state ultrasoft and hard materials) and hyperbolic optical response restrict the wide rangeof possible combinations to 17 superlattices, out of 1891 combinations, most of which include nitridecomponents. A few compounds include quite exotic crystal structures such as RuC and ZnN, ormonoborides whose stability has to be confirmed by experiments. The majority of the systems arecomposed of transition metal nitrides and carbides (e.g. HfN, TiN, ZrN, NbC, TiC, TaC) that areeasily grown and whose stability is largely demonstrated. A few of the proposed superlattices, such as15fN/ScN, ZrN/ScN, seem to be particularly promising: they have been experimentally realized, andare expected to be simultaneously mechanically hard and optically hyperbolic over a large range offrequencies. This calls for further experimental validations and characterizations. The choice of thematerial layers and their relative thickness can be fine tuned to fulfill ad hoc mechanical and opticalrequirements to optimize the performance of the final optoelectronic application. This opens up waysfor the design and growth of hyperbolic materials with superior mechanical properties that can beused in extreme conditions as for aerospace or security applications.
METHODOLOGY AND COMPUTATIONAL DETAILS
Electronic structure calculations were carried out by using a first-principles total-energy-and-forcesapproach based on density functional theory (DFT), as implemented in the QUANTUM ESPRESSO(QE) suite of codes [34]. The Perdew, Burke, and Ernzerhof (PBE) [95] generalized gradient approxi-mation was adopted for parametrization of the exchange-correlation functional. Ionic potentials weredescribed by normconserving Trouillier-Martin pseudopotentials. Single particle Khon-Sham orbitalsare expanded in plane waves up to a kinetic energy cutoff of 150.0 Ry. A uniform (24 × × F m ¯3 m , number 225), with two inequivalent atoms in the primitive cell. Theequilibrium lattice constance (a ) of each system is obtained through a variable cell optimization [96].See Sec. S1 of SI for further computational details.In the case of semiconducting compounds (ScN, Yn, LaN), the electronic properties have beencalculated ab initio by using a recent pseudohybrid Hubbard implementation of DFT+U, namelyACBN0 [97], that profitably corrects the energy bandgap [97, 98] as well as the dielectric and vibrationalproperties of semiconductors [99]. The optimized values for the studied compounds are U N (ScN)=3.1eV, U Sc (ScN)=0.1 eV, U N (YN)=3.10 eV, U Y (YN)=0.1 eV, U N (LaN)=3.3 eV, U La (ScN)=0.0 eV. Thespin degrees of freedom for magnetic compounds (MnB, LaB, CrC, MnC, FeC, CrN, MnN, FeN) havebeen described within the local spin-density approximation (LSDA). Preliminary tests confirmed thatspin-orbit coupling (SOC) contributions give negligible corrections to the ground state of single crystalsin the fcc structure.The optical properties of rocksalt crystals are determined using the epsilon.x code, contained in theQE package, which implements a band-to-band independent-particle (IP) formulation of the frequency-dependent Drude-Lorentz model for the dielectric function ˆ (cid:15) ( ω ), where both intraband (Drude-like)and interband (Lorentz-like) contributions are explicitly considered [100]:ˆ (cid:15) ( ω ) = 1 − (cid:88) k ,n f n,n k ω p ω + iηω + (cid:88) k ,n (cid:54) = n (cid:48) f n,n (cid:48) k ω p ω k ,n,n (cid:48) − ω − i Γ ω , (9)where ω p is the bulk plasma frequency; (cid:126) ω k ,n,n (cid:48) = E k ,n − E k ,n (cid:48) is the vertical band-to-band transitionenergy between occupied and empty Bloch states calculated at the DFT level; { k , n } and { k , n (cid:48) } arethe k-point and band quantum numbers, respectively. η, Γ → + are the Drude-like and Lorentz-likerelaxation terms, while f n,n k and f n,n (cid:48) k are the corresponding oscillator strengths within the dipoleapproximation (see S.I. for further details).Despite the simplicity, the capability of the present approach in simulating the optical properties ofplasmonic materials has been extensively proved, e.g., in Refs [100–102]. In the particular case of TiN,the accuracy of this approach has been tested in Refs. [103] and [104], where the dielectric function ofTiN bulk calculated with the Drude-Lorentz approach is compared to one based on Time-DependentDensity Functional Perturbation Theory and to the experimental data [105, 106]. A similar Drude-Lorentz approach has been profitably used by other authors for simulating the dielectric function of16etallic rocksalts (e.g. ZrN, HfN, TaC, WC) [107], in very good agreement with the present results. Inthe case of ScN (see SI, Sec. S2), the comparison with experimental results [108] confirms the accuracyof the present approach also in the case of small gap semiconductors.Second order elastic stiffness (c ij ) and compliance (s ij ) coefficients have been calculated by meansof a polynomial fitting of the strain-stress relation of deformed crystals, as implemented in the ElaS-tic simulation package [35]. For each rocksalt crystal, we considered 11 distorted structures with aLagrangian strain value between − η max and + η max , with η max = 1 × − . The symmetry-dependentdeformation types follows the universal linear-independent coupling strains (ULICS) classification,proposed in Ref. [109]. The energies and stresses of the distorted structures are calculated ab initioby using the QE code (see S.I. for further details).High throughput DFT, optical, and ElaStic calculations have been run by using the automaticworkflows implemented in the AFLOW π infrastructure [110]. Supporting Information
Supporting Information is available from the Wiley Online Library or from the author. The Sup-porting Information file includes: the description of the DFT implementation of the optical and elasticproperties of crystals (Sec. S1); the summary of the main structural, electronic and optical propertiesof single TMX systems (Sec. S2); the complete list of the elastic and mechanical properties of singlerocksalt crystals (Sec. S3); the detailed description of the mean theory approach for the evaluationof the elastic properties in superlattices (Sec. S4); the complete list of the elastic and mechanicalproperties of all TMX-based superlattices (Sec. S5); and the description of the gain factors (Sec. S6)discussed in the text.
Acknowledgements
We acknowledge the support of the High Performance Computing Center at the University of NorthTexas and the Texas Advanced Computing Center at the University of Texas, Austin. The authorsthank Rita Stacchezzini for graphical help as well as Andrew Supka and Sharad Mahatara for technicalsupport and discussions. ∗ Email: [email protected][1] J. Bardeen, Nobel Lecture , 318 (1956).[2] P. A. Grunberg, Nobel Lecture , 92 (2007).[3] V. Veselago, Sov. Phys. Usp. , 517 (21967).[4] S. A. Ramakrishna, Rep. Progr. Phys. , 449 (2005).[5] J. B. Pendry, Phys. Rev. Lett. , 3966 (2000).[6] N. Engheta, IEEE Ant. Wireless Prop. Lett. , 10 (2002).[7] L. Ferrari, C. Wu, D. Lepage, X. Zhang, and Z. Liu, Prog. Quant. Elect. , 1 (2015).[8] J. B. Pendry and D. R. S. D. Schurig, Science , 1780 (2006).[9] S. Bang, S. So, and J. Rho, Sci. Rep. , 1 (2019).[10] X. Yang, J. Yao, J. Rho, X. Yin, and X. Zhang, Nature Photonics. , 450 (2012).[11] A. V. Kabashin, P. Evans, S. Pastkovsky, W. Hendren, G. A. Wurtz, R. Atkinson, R. Pollard, V. A.Podolskiy, and A. V. Zayats, Nature Mater. , 867 (2009).[12] A. A. Govyadinov and V. A. Podolskiy, Phys. Rev. B , 155108 (2006).[13] J. Rho, Z. Ye, Y. Xiong, X. Yin, Z. Liu, H. Choi, G. Bartal, and X. Zhang, Nature Commun. , 143(2010).[14] W. Cai, U. Chettiar, A. Kildishev, and V. Shalaev, Nature Photonics , 224 (2007).[15] Z. Jacob, L. V. Alekseyev, and E. Narimanov, Opt. Expr. , 8247 (2006).
16] A. Al`u, A. Salandrino, and N. Engheta, Optics Express , 1557 (2006).[17] M. Amiri, F. Tofigh, N. Shariati, J. Lipman, and M. Abolhasan, IEEE Internet of Things Journal , 1(2020).[18] I. Avrutsky, I. Salakhutdinov, J. Elser, and V. Podolskiy, Phys. Rev. B , 241402 (2007).[19] S. V. Zhukovsky, O. Kidwai, and J. E. Sipe, Opt. Exp. , 14982 (2013).[20] H. N. S. Krishnamoorthy, Z. Jacob, E. Narimanov, I. Kretzschmar, and V. M. Menon, Science ,205 (2012).[21] Y. Guo, C. L. Cortes, S. Molesky, and Z. Jacob, Appl. Phys. Lett. , 131106 (2012).[22] T. Taubner, D. Korobkin, Y. Urzhumov, G. Shvets, and R. Hillenbrand, Science , 1595 (2006).[23] X. Liu, J. Wang, L. Tang, L. Xie, and Y. Ying, Adv. Funct. Mater. , 5515 (2016).[24] Z. Jacob and E. E. Narimanov, Appl. Phys. Lett. , 1595 (2008).[25] F. Monticone and A. Al`u, Rep. Progr. Phys. , 036401 (2017).[26] I. M. Yusupov and D. S. Filonov, IEEE Conf. Paper. , 169 (2020).[27] C. A. Valagiannopoulos, M. S. Mirmoosa, I. S. Nefedov, S. A. Tretyakov, and C. R. Simovski, J. Appl.Phys. , 163106 (2014).[28] D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith,Science , 977 (2006).[29] R. Liu, C. Ji, Z. Zhao, and T. Zhou, Engineering , 179 (2015).[30] K. V. Sreekanth, Y. Alapan, M. ElKabbash, E. Ilker, M. Hinczewski, U. A. Gurkan, A. De Luca, andG. Strangi, Nature Materials , 621 (2016).[31] F. M. Gao and L. H. Gao, J. Superhard Mater. , 148 (2010).[32] M. T. Yeung, R. Mohammadi, and R. B. Kaner, Ann. Rev. Mater. Res. , 465 (2016).[33] G. V. Naik, B. Saha, J. Liu, S. M. Saber, E. A. Stach, J. M. K. Irudayaraj, T. D. Sands, V. M. Shalaev,and A. Boltasseva, Proc. Nat. Ac. Sci. USA , 7546 (2014).[34] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. Chiarotti,M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann,C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini,A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov,P. Umari, and R. M. Wentzcovitch, J. Phys.: Cond. Matt. , 395502 (2009).[35] R. Golesorkhtabar, P. Pavone, J. Spitaler, P. Puschnig, and C. Draxl, Comp. Phys. Commun. ,1861 (2013).[36] Z. Wu, X.-J. Chen, V. V. Struzhkin, and R. E. Cohen, Phys. Rev. B , 214103 (2005).[37] X. Chen, V. V. Struzhkin, Z. Wu, M. Somayazulu, J. Qian, S. Kung, A. N. Christensen, Y. Zhao, R. E.Cohen, H.-k. Mao, and R. J. Hemley, Proc. Nat. Ac. Sci. USA , 3198 (2005).[38] C. Y. Allison, F. A. Modine, and R. H. French, Phys. Re. B , 2573 (1987).[39] M. Iuga, G. Steinle-Neumann, and J. Meinhardt, Eur. Phys. J. B , 127 (2007).[40] J. Chang, G.-P. Zhao, X.-L. Zhou, K. Liu, and L.-Y. Lu, J. Appl. Phys. , 083519 (2012).[41] M. G. Quesne, A. Roldan, N. H. de Leeuw, and C. R. A. Catlow, Phys. Chem. Chem. Phys. , 6905(2018).[42] D. A. Papaconstantopoulos, W. E. Pickett, B. M. Klein, and L. L. Boyer, Phys. Rev. B , 752 (1985).[43] F. Cardarelli, Materials Handbook , A Concise Desktop Reference (Springer Science & Business Media,2008).[44] P. Mohn and D. G. Pettifor, J. Phys. C: Sol. St. Phys. , 2829 (2000).[45] A. H. Sihvola, Electromagnetic mixing formulas and applications (IET., London, UK, 1999).[46] A. Catellani and A. Calzolari, Opt. Mater. Exp. , 1459 (2019).[47] A. Poddubny, I. Iorsh, P. Belov, and Y. Kivshar, Nature Photonics , 948 (2013).
48] M. Grimsditch and F. Nizzoli, Phys. Rev. B , 5891 (1986).[49] J. Buchinger, N. Koutn´a, Z. Chen, Z. Zhang, P. H. Mayrhofer, D. Holec, and M. Bartosik, ActaMaterialia , 18 (2019).[50] A. Mohamed, M. Dutta, and M. Stroscio, IEEE Electron Device Letters , 1175 (2019).[51] C. Prieto, A. d. Bernab´e, R. Casta˜ner, A. Mu˜noz-Mart´ın, R. J. Jim´enez-Riobo´o, M. Garc´ıa-Hern´andez,and A. d. Andr´es, J. Physics: Cond. Matt. , 2931 (2000).[52] F. Martin, C. Jaouen, J. Pacaud, G. Abadias, P. Djemia, and F. Ganot, Phys. Rev. B , 33 (2005).[53] C. M. Kube, AIP Advances , 095209 (2016).[54] S. I. Ranganathan and M. Ostoja-Starzewski, Phys. Rev. Lett. , 055504 (2008).[55] A. H. Cottrell, Introduction to the modern theory of metals (Institute of Metals, London, UK, 1988).[56] A. C. Fischer-Cripps,
Nanoindentation (Springer, 2004).[57] Y. Tian, B. Xu, and Z. Zhao, Int. J. Refract. Metals Hard Mater. , 93 (2012).[58] V. Kanyanta, Microstructure-Property Correlations for Hard, Superhard, and Ultrahard Materials (Springer, Cham, 2016).[59] Y. Feng, T. Zhang, and R. Yang, J. Am. Ceram. Soc. , 332 (2011).[60] H. Niu, S. Niu, and A. R. Oganov, J. Appl. Phys. , 065105 (2019).[61] V. V. Brazhkin and V. L. Solozhenko, J. Appl. Phys. , 130901 (2019).[62] M. Shinn, L. Hultman, and S. Barnett, J. Mater.Res. , 901?911 (1992).[63] J. O. Kim, J. D. Achenbach, M. Shinn, and S. A. Barnett, J. Mater. Res. , 2248 (2011).[64] U. Helmersson, S. Todorova, S. A. Barnett, J. E. Sundgren, L. C. Markert, and J. E. Greene, J. Appl.Phys. , 481 (1987).[65] J. Musil, Surf. Coat. Tech. , 322 (2000).[66] X. Zhao, Y. Zhuo, S. Liu, Y. Zhou, C. Zhao, C. Wang, and Q. Yang, Surf. Coat. Tech. , 200 (2016).[67] Q. Yang and L. R. Zhao, J. Vac. Sci. Technol. A , 558 (2003).[68] B. Saha, G. V. Naik, S. Saber, C. Akatay, E. A. Stach, V. M. Shalaev, A. Boltasseva, and T. D. Sands,Phys. Rev. B , 125420 (2014).[69] X. Chu and S. A. Barnett, J. Appli. Phys. , 4403 (1995).[70] X. Chu, M. S. Wong, W. D. Sproul, and S. A. Barnett, J. Mater. Res. , 2500 (1999).[71] S. Chakraborty, H. Uchiyama, M. Garbrecht, V. Bhatia, A. I. K. Pillai, J. P. Feser, D. T. Adroja,S. Langridge, and B. Saha, Appl. Phys. Lett. , 1 (2020).[72] M. Garbrecht, B. Saha, J. L. Schroeder, L. Hultman, and T. D. Sands, Sci. Rep. , 1 (2017).[73] M. Garbrecht, J. L. Schroeder, L. Hultman, J. Birch, B. Saha, and T. D. Sands, J. Mater. Sci. , 8250(2016).[74] M. Azadi, A. S. Rouhaghdam, and S. Ahangarani, Strength of Materials , 279 (2016).[75] Y. Zhao, G. Lin, J. Xiao, C. Dong, and L. Wen, Vaccum , 1 (2010).[76] J. Keem, in T. S. Sudarshan (ed.), Surface modification Technologies: An Engineers Guide (MarcelDekker, New York, US, 1990).[77] Y. Zhang, D. Liu, Y. Shi, Z. Sun, L. Wu, and Y. Gao, App. Phys. Exp. , 015501 (2019).[78] C. A. Sweeney, W. Vorster, S. B. Leen, E. Sakurada, P. E. McHugh, and F. P. E. Dunne, J. Mech.Phys. Sol. , 1224 (2013).[79] S. Groh, B. Devincre, L. P. Kubin, A. Roos, F. Feyel, and J. L. Chaboche, Phil. Mag. Lett. , 303(2010).[80] J. Haines, J. L´eger, and G. Bocquillon, Annu. Rev. Mater. Res. , 1 (2001).[81] M. Shinn and S. A. Barnett, Appl. Phys. Lett. , 61 (1994).[82] J. Musil, Surf. Coatings Tech. , 322 (2000).[83] M. Shinn, L. Hultman, and S. Barnett, J. Mater. Res. , 901 (1992).
84] S. Zhang, D. Sun, Y. Fu, and H. Du, Surf. Coat. Tech. , 113 (2003).[85] C. Subramanian, K. S. Wear, and 1993, Wear , 85 (1993).[86] J. S. Koehler, Phys. Rev. B , 547 (1970).[87] J. W. Cahn, Acta Met. , 1275 (1963).[88] K. Yoshii, H. Takagi, M. Umeno, and H. Kawabe, Metall. Trans. A , 1273 (1984).[89] A. J. Hoffman, L. Alekseyev, S. S. Howard, K. J. Franz, D. Wasserman, V. A. Podolskiy, E. E. Nari-manov, D. L. Sivco, and C. Gmachl, Nature Materials , 946 (2007).[90] V. P. Drachev, V. A. Podolskiy, and A. V. Kildishev, Optics Express , 15048 (2013).[91] K. Korzeb, M. Gajc, and D. A. Pawlak, Optics Express , 25406 (2015).[92] S. Ishii, A. V. Kildishev, E. Narimanov, V. M. Shalaev, and V. P. Drachev, Laser & Phot. Rev. , 265(2013).[93] I. V. Lindell, S. A. Tretyakov, K. I. Nikoskinen, and S. Ilvonen, Microwave Opt. Tech. Lett. , 129(2001).[94] S. V. Zhukovsky, O. Kidwai, and J. E. Sipe, Optics Express , 14982 (2013).[95] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996).[96] R. Wentzcovitch, J. Martins, and G. Price, Phys. Rev. Lett. , 3947 (1993).[97] L. A. Agapito, S. Curtarolo, and M. Buongiorno Nardelli, Phys. Rev. X , 011006 (2015).[98] P. Gopal, M. Fornari, S. Curtarolo, L. A. Agapito, L. S. I. Liyanage, and M. Buongiorno Nardelli,Phys. Rev. B , 245202 (2015).[99] A. Calzolari and M. Buongiorno Nardelli, Sci. Rep. , 2999 (2013).[100] A. Calzolari, A. Ruini, and A. Catellani, ACS Photonics , 703 (2014).[101] A. Catellani, P. D’Amico, and A. Calzolari, Phys. Rev. Mater. , 015201 (2020).[102] M. Eaton, A. Catellani, , and A. Calzolari, Opt. Express , 5342 (2018).[103] A. Catellani and A. Calzolari, Phys. Rev. B , 115145 (2017).[104] D. Shah, A. Catellani, H. Reddy, N. Kinsey, V. Shalaev, A. Boltasseva, and A. Calzolari, ACS Photonics , 2816 (2018).[105] J. Pfl¨uger, J. Fink, W. Weber, K. P. Bohnen, and G. Crecelius, Phys. Rev. B , 1155 (1984).[106] A. A. Herzing, U. Guler, X. Zhou, A. Boltasseva, V. Shalaev, and T. B. Norris, Appl. Phys. Lett. ,171107 (2016).[107] M. Kumar, N. Umezawa, S. Ishii, and T. Nagao, ACS Photonics , 43 (2015).[108] D. Gall, M. St¨adele, K. J¨arrendahl, I. Petrov, P. Desjardins, R. T. Haasch, T.-Y. Lee, and J. E. Greene,Phys. Rev. B , 125119 (2001).[109] R. Yu, J. Zhu, and H. Q. Ye, Comp. Phys. Commun. , 671 (2010).[110] A. R. Supka, T. E. Lyons, L. Liyanage, P. D’Amico, R. A. R. Al Orabi, S. Mahatara, P. Gopal, C. Toher,D. Ceresoli, A. Calzolari, S. Curtarolo, M. Buongiorno Nardelli, and M. Fornari, Comput. Mater. Sci. , 76 (2017)., 76 (2017).