1. Introduction

In the last few decades the scientific community has focused significant effort on the design and fabrication of artificial nanostructures, with the intent of circumventing limitations imposed by natural materials, like absorption, and provide novel functionalities by harnessing processes that may occur only at the nanoscale, while simultaneously reducing the size of optical devices. Among the plethora of proposed novel materials and nanostructures, those operating near their epsilon-near-zero (ENZ) region have been shown to possess a set of peculiar properties [1–3]. These materials may be exploited to manipulate radiation direction and polarization [4–8], but may also serve as platforms to study nonlinear optical interactions both in the continuous wave and in the pulsed regimes [9–11] thanks to their ability to enhance local fields under specific circumstances, i.e., TM-polarized field incident at oblique incidence [3]. Among the explored nonlinear processes [12–19] that have been studied both theoretically and experimentally [20–29], second and third harmonic generation stand out. The ENZ condition is known to manifest itself in any natural material. For example, GaAs, Si, and GaP have their zero crossing for the real part of the permittivity in the UV, metals like Au, Ag and Cu in the visible, and highly-doped oxides like ITO and AZO in the infrared [30]. However, it is possible to design artificial nanostructures and tailor the frequency response at will [31] or to compensate losses [21,32,33], which are thought to represent the main limitation for field enhancement in ENZ materials [34]. Still, loss compensation techniques based on the inclusion of active materials might be challenging for fabrication and not practical for certain applications.

An alternative route to achieving high local field enhancement without resorting to loss compensation techniques is to exploit the effective anisotropy of artificial nanostructures via the realization of longitudinal epsilon-near-zero, or LENZ, materials [35]. Indeed, it has been shown that LENZ possess exceptional abilities that span from perfect light bending [36,37] and angular filtering and polarization control [38], to coherent perfect absorption [39] and control of leaky wave radiation [40]. Moreover, since the effective anisotropy tends to enhance the local field intensity and improve tolerances in terms of acceptance angles and a weaker sensitivity to material thickness [35], LENZ turn out to be good candidates for the investigation of nonlinear processes.

Recently, LENZ have been studied in the context of second order nonlinear optical processes arising from symmetry breaking at each interface of a layered structure. It was shown that a Si/Dy:CdO (dysprosium-doped Cadmium Oxide) multilayer can boost the enhancement of second harmonic generation (SHG) beyond their isotropic counterparts [41]. Here we start by investigating third harmonic generation (THG) with the intent to further examine the role of anisotropy on nonlinear processes in ENZ media. We simulate an effective anisotropic medium and demonstrate that, regardless of the thickness of the medium, THG can be enhanced dramatically by properly increasing the material’s degree of anisotropy, and demonstrate that nanostructures with metallic inclusions can further boost nonlinear interactions with respect to nanostructures that show a similar degree of anisotropy but are composed of dielectrics and semiconductors. We then analyze a realistic scenario, where a metal/conducting-oxide multilayer stack composed of Au and Dy:CdO is used to test the impact of anisotropy on the same third order nonlinear process. We compare the response of the multilayer with and without the effects of nonlocality, and introduce the contribution of hot electrons, two processes that can be triggered in both metals and conducting oxides, and compete to either blueshift or redshift the resonant linear and nonlinear responses of the structure. We demonstrate that hot electrons contributions can further boost the nonlinear response by several orders of magnitude and, therefore, their inclusion in the theoretical investigation is crucial to make accurate plans and predictions for experiments. Finally, we demonstrate that the high degree of anisotropy achieved in a multilayer environment also enhances second order nonlinear processes even in the absence of bulk nonlinearities.

2. Improving nonlinear phenomena through anisotropy

Anisotropic ENZ or LENZ have been shown to increase SHG conversion efficiencies in the absence of dipole-allowed quadratic nonlinearity [41]. This improvement is possible thanks to the high local fields that can be achieved in anisotropic ENZ [35] even without the aid of geometric resonances. Although it is clear from previous investigations that anisotropy facilitates nonlinear optical processes in the presence of the ENZ condition, to date it is still not clear if nanostructures with a low degree of anisotropy can perform as well as structures with a high degree of anisotropy, and whether or not the improvement associated with the anisotropy also applies to bulk nonlinearities.

In order to shed light on these two aspects we first investigate THG in a homogenous slab of material having the dielectric permittivity of Dy:CdO [42]. Then we artificially modify the transverse components of the dielectric permittivity (ɛ_z= ɛ_Dy:CdO, ɛ_x = ɛ_y = ɛ_Dy:CdO +Δɛ) to simulate a Dy:CdO-like material and introduce an arbitrary degree of anisotropy. For simplicity, in this section we assume the material is isotropic at the third harmonic frequency. The geometry of the structure under consideration is depicted in Fig. 1: a slab of thickness d is illuminated by a TM-polarized field incident at a variable angle ϑ_i.

 

Fig. 1. Sketch of the structure under investigation: a slab of material of thickness d is illuminated by a TM-polarized pump wavevector k_FF at an angle ϑ_i. Forward and backward third harmonic signals are generated with wavevector ± k_TH, respectively.

Download Full Size | PDF

The fundamental frequency (FF) corresponds to λ_FF = 2147 nm, where the real part of the dielectric permittivity of Dy:CdO displays a zero crossing. Pump irradiance is set to 1 GW/cm². In this preliminary analysis we solve both linear and nonlinear electromagnetic problems by using a frequency-domain, finite-element solver (COMSOL Multiphysics). Since we are simulating a highly-doped oxide, we assume the material possesses an isotropic, bulk, dispersion-less and near-instantaneous $\chi _{}^{(3)} = {10^{ - 20}}\textrm{ }{\textrm{m}^\textrm{2}}\textrm{/}{\textrm{V}^\textrm{2}}$ at both FF and third harmonic (TH) frequency. This value is compatible with those measured for other doped oxides, such as ITO and AZO [17,20]. No dipole-allowed second-order nonlinearity is considered. The TM-polarized field is assumed to be a superposition of two monochromatic signals at the FF and TH harmonic frequencies.

The TH signal is initially null. Therefore, when solving Maxwell’s equations, the time-independent, nonlinear complex current densities at FF and TH may be written as [43]:

In order to understand the impact of anisotropy we begin our investigation by calculating the field intensity enhancement, FIE =$|{{\left\langle {{E_z}} \right\rangle } / {{E_0}{|^2}}}$, where $\left\langle {{E_z}} \right\rangle $ is the average value of the component of the electric field in the direction of propagation (see Fig. 1), evaluated inside the slab. E₀ is the amplitude of the incident electric field. In Fig. 2 we show FIE maps for different slab thicknesses, respectively for d = 20 nm [Fig. 2(a)], d = 60 nm [Fig. 2(c)], and d = 100 nm [Fig. 2(e)]. FIE is calculated as a function of the angle of incidence and degree of anisotropy Δɛ. Specifically, the larger the value of |Δɛ| is, the larger the difference will be between longitudinal (direction of propagation) and transverse permittivities. As one departs from the isotropic condition (Δɛ = 0), FIE increases for a certain range of angles of incidence, regardless of slab thickness. We note that the sign of the anisotropy also appears to impact FIE, which attains slightly larger values when Δɛ is positive. The plots also confirm that the tolerance to the angle of incidence improves as the thickness of the slab is increased, as indicated in Ref. [35]. Although FIE improvement is generally a good guideline when attempting to boost nonlinear optical processes, the maps of total (forward plus backward) THG for thicker slabs [Figs. 2(d) and (f)] show a somewhat different trend: while Figs. 2(c) and (e), show that FIE reaches higher values when $\Delta \varepsilon \gg 0$, the THG maps in Figs. 2(d) and (f) reveal that total third harmonic conversion efficiency generally improves when the degree of anisotropy increases, but achieves larger values when $\Delta \varepsilon \ll 0$.

 

Fig. 2. FIE for a slab of material with variable anisotropy and angle of incidence for thicknesses (a) d = 20 nm, (c) d = 60 nm, and (e) d = 100 nm; Total third harmonic generated signal (forward plus backward) for a slab of material with variable anisotropy and angle of incidence and thickness (b) d = 20 nm, (d) d = 60 nm, and (f) d = 100 nm. All maps have been calculated using the finite-element solver COMSOL Multiphysics.

Download Full Size | PDF

The reason behind the discrepancy between the FIE and THG trends in thicker slabs resides in the definition of FIE, which considers the average value of the longitudinal electric field inside the slab. In fact, the maximum value of the longitudinal electric field component E_z (not shown), which occurs at the first interface between vacuum and the slab, is largest for large negative anisotropy ($\Delta \varepsilon \ll 0$), and produces a stronger TH signal. Put another way, while a large degree of anisotropy is desirable in order to improve FIE, the sign of the anisotropy may also be relevant to determine the efficiency of the nonlinear optical process, especially in thicker slabs.

3. Third-harmonic generation from Au/Dy:CdO multilayer stack

A careful look at Fig. 2 reveals that relatively high |Δɛ| are generally needed to significantly improve the efficiency of nonlinear processes compared to the isotropic case (Δɛ = 0). For example, a 100nm thick slab shows a maximum THG efficiency occurring near 65°, when $\Delta\varepsilon \, \cong \, -100$. At first sight, the problem of assembling a practical device with these characteristics may appear to be daunting. The requirement is to achieve high anisotropy while preserving a zero-crossing in the real part of the effective permittivity in the longitudinal direction. However, if we aim to realize an anisotropic material with $\Delta \varepsilon \ll 0$, one may resort to the inclusion of metals because they naturally display large, negative permittivities in practically all spectral ranges of interest. For this reason we direct our attention to structures similar to that depicted in Fig. 3: a five periods metal-conducting oxide stack composed of Au [30] and Dy:CdO [42], each 10 nm thick (a = b = 10 nm), for a total thickness d = 100 nm, comparable in size to the slab of Figs. 2(e) and (f). The multilayer is situated on top of a SiO₂ semi-infinite substrate, and it is illuminated by a TM-polarized field incident at a variable angle ϑ_i.

 

Fig. 3. Sketch of the multilayer structure: five periods of Dy:CdO (a=10 nm) and Au (b=10 nm) are alternated to obtain an anisotropic response. The structure is illuminated by a TM-polarized pump wavevector k_FF at an angle ϑ_i. Only backward third harmonic signal with wavevector k_TH is monitored.

Download Full Size | PDF

The effective permittivities of the stack calculated at λ_FF = 2147 nm (the crossing point of Dy:CdO) are ɛ_z,FF= 2.1×10⁻⁷ + 0.34i, ɛ_x,FF = ɛ_y,FF = −101.34 + 11.34i at the FF, and ɛ_z,TH= 12.46 + 0.51i, ɛ_x,TH = ɛ_y,TH = −6.52 + 0.92i at the TH. We stress that the solution presented here is by no means unique and similar values for the longitudinal and transverse permittivities of the multilayer may be obtained by considering a different set of materials, as shown in Appendix A. Additionally, while in the previous section the anisotropy at the TH frequency was ignored, a realistic nanostructure will indeed display anisotropic behavior also at the TH as the calculated effective permittivity values indicate. To boot, in our preliminary analysis we also assumed a purely real Δɛ. The introduction of metals in the multilayer stack adds a realistic degree of damping in the system, as revealed by the imaginary part of the effective permittivity in the transverse direction.

Metals like gold and silver typically require one free and two bound electron contributions in order to accurately reproduce the local dielectric constant over a broad range that includes UV wavelengths, where interband transitions significantly contribute to the dielectric response. Conducting oxides, on the other hand, typically require one Drude and one Lorentz electron species. Therefore, gold data were modeled by using a combination of a Drude (${\omega _{p,f}} = 1.27 \times {10^{16}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\gamma _f} = 9.43 \times {10^{13}}\textrm{ }{\textrm{s}^{ - 1}}$) and two Lorentz oscillators (${\omega _{p,b1}} = 7.73 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\omega _{0,b1}} = 5.18 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\gamma _{b1}} = 2.26 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\omega _{p,b2}} = 1.12 \times {10^{16}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\omega _{0,b2}} = 8.48 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\gamma _{b2}} = 2.26 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$) that provide a good fit for data down to approximately 200nm [30]. On the other hand, Dy:CdO [42] data were fit with a single Drude oscillator (${\omega _{p,f}} = 1.99 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\gamma _f} = 2.73 \times {10^{13}}\textrm{ }{\textrm{s}^{ - 1}}$) for free electrons, and one Lorentz oscillator for bound electrons with a resonance far from the ENZ condition (${\omega _{p,b1}} = 8.87 \times {10^{16}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\omega _{0,b1}} = 4.33 \times {10^{16}}\textrm{ }{\textrm{s}^{ - 1}}$, ${\gamma _{b1}} = 1.88 \times {10^{15}}\textrm{ }{\textrm{s}^{ - 1}}$). As it was done for the single slab, we evaluate both linear and nonlinear responses of the structure so that we may keep track of how the two regimes are related. In order to model a realistic stack we resort to a detailed, microscopic hydrodynamic description of light-matter interactions that has previously been used to simulate wave propagation in an ITO nanolayer [44–47]. We also simulated both linear and nonlinear electromagnetic problems by means of a frequency-domain, finite-element solver (COMSOL Multiphysics). This will help us verify that the results are robust and consistent in the continuous wave and pulsed regimes. Using the finite element solver, gold and Dy:CdO layers are modeled with a near-instantaneous, isotropic nonlinear response with $\chi _{}^{(3)} = {10^{ - 20}}\textrm{ }{\textrm{m}^\textrm{2}}\textrm{/}{\textrm{V}^\textrm{2}}$. On the other hand, the hydrodynamic approach allows for a realistic description of both linear and nonlinear dispersions without a preconceived notion of the magnitude of $\chi _{}^{(3)}$. Both approaches yield qualitatively similar results. We stress that, while gold has been reported to have a nonlinear susceptibility of the order of 10⁻¹⁶ m²/V² in the visible range, its magnitude generally depends on pulse duration and may decrease substantially away from the plasma frequency of the material [17], a fact that can easily be verified using the Drude-Lorentz oscillator approach. In addition, the fields are strongly localized inside the CdO layers, a circumstance that greatly diminishes the nonlinear response of gold.

Pump irradiance is assumed once again to be 1 GW/cm². FIE and THG maps as functions of incident wavelength and angle of incidence are shown in Figs. 4(a) and 5(a), respectively. In the multilayer, FIE =$|{{\left\langle {{E_z}} \right\rangle } / {{E_0}{|^2}}}$ is obtained by evaluating the average value of the component of the electric field in the direction of propagation $\left\langle {{E_z}} \right\rangle $ inside all layers of the stack, keeping in mind, however, that the field is localized mostly inside the Dy:CdO layers. Both FIE and THG peak in the vicinity of the zero-crossing for the effective permittivity in the propagation direction. The peaks are also relatively wide as a function of the angle of incidence, in contrast with the response of isotropic ENZ materials [41]. We note that the TH conversion efficiency shown in Fig. 5(a) refers only to the TH reflected signal. The transmitted TH signal is two orders of magnitude smaller than the reflected component and is, therefore, ignored. Although the peak values shown in Fig. 4(a) and Fig. 5(a) are similar to the results presented in the previous section for the 100 nm slab having similar values of longitudinal and transverse permittivities, [see white dots in Figs. 2(e) and (f)], we are now analyzing a nanostructure that presents a non-homogeneous electric field distribution. Therefore, the concept of effective permittivity falls short of an adequate description of nonlinear interactions, and may only be used as a guideline to design multilayers. Care should be exercised when predicting the nonlinear response from the structure using effective medium approaches. In the multilayer scenario, three important aspects need to be taken into account for both linear and nonlinear simulations: (i) the thickness of each layer is 10 nm, so that nonlocal contributions should be taken into account in both linear and nonlinear regimes; (ii) it is known that hot carrier contributions may be triggered in both Au and Dy:CdO, significantly altering the linear and nonlinear responses [45]; (iii) a second-order response that arises from symmetry breaking at the interfaces in not negligible is this kind of structure [41], may lead to cascaded THG [29], and must be included in the calculations [45,47,48].

 

Fig. 4. FIE for the multilayer structure sketched in Fig. 3, calculated as a function of incident wavelength and angle of incidence in the (a) local, (b) nonlocal ($\chi _{hot,FF/TH}^{(3)} = 0$) and (c) nonlocal scenario where hot carrier contributions are included. (b) and (c) reveal how nonlocality blueshifts the resonance, while the nonlinearity due to hot electrons imparts a redshift of the same resonant feature. All maps have been calculated using the finite-element solver COMSOL Multiphysics. The hydrodynamic model yields similar results.

Download Full Size | PDF

 

Fig. 5. TH conversion efficiency maps calculated as functions of pump wavelength and angle of incidence in: (a) local scenario [Eqs. (1) and (2) are implemented], (b) nonlocal scenario [Eqs. (3) and (4) are implemented, without hot carrier contributions, i.e. $\chi _{hot,FF/TH}^{(3)} = 0$] and (c) nonlocal scenario where hot carrier contributions are included - Eqs. (3) and (4) implemented as shown. (d) sections of TH conversion efficiency maps assuming ϑ_i= 65° (blue line, circle markers - local scenario), ϑ_i= 69° (red line, square markers - nonlocal scenario) and ϑ_i= 71° (green line, triangle markers - nonlocal scenario assuming hot electrons contributions). (e) sections of THG conversion efficiency for ϑ_i= 65°, for different pulse durations, as indicated, calculated in the nonlocal scenario assuming hot electrons contributions. Maps in (a), (b) and (c) and sections in (d) have been calculated using the finite-element solver COMSOL Multiphysics. Curves in (e) have been calculated by means of the microscopic, hydrodynamic approach.

Download Full Size | PDF

Therefore, in addition to the nonlinear contributions that arise from bound electrons, which peak in a different wavelength range [45] and are taken into account directly into Maxwell’s equations, the nonlinearities associated with Eqs. (1) and (2) are now related directly to hot electrons, and are treated in the Drude portion of the simulation. Adding the usual nonlocal terms to Eqs. (1) and (2), with the understanding that we are now dealing with free electrons, we obtain the following novel equations of motion for the nonlinear, free currents:

The introduction of nonlocal and hot electron contributions alters significantly both linear and nonlinear responses of the nanostructure. In what follows, we will discuss and compare three different regimes: (i) local regime – where only Eqs. (1) and (2) are implemented, (ii) nonlocal regime, where Eqs. (3) and (4) are implemented but we assume $\chi _{hot}^{(3)} = 0$ and (iii) nonlocal + he regime, where Eqs. (3) and (4) are implemented as shown and $\chi _{hot}^{(3)} = {10^{ - 18}}$ m²/V².

We start by examining the average FIE in the nonlocal regime and with the model that includes both nonlocal and hot electron contributions, respectively [Figs. 4(b) and (c)]. When comparing the maps in Fig. 4, two macroscopic differences may be discerned: (i) the nonlocality blue-shifts the structure’ spectral response [50,51]; (ii) FIE values drop for all wavelengths and angles of incidence. Moreover, while hot carrier contributions do not seem to alter the maximum FIE that may be achieved in the multilayer, significant distortion of both spectral and angular features are noticeable when hot electrons contributions are introduced [Fig. 4(c)]. Note that the degree of distortion is proportional to pump irradiance value (not shown) [45].

The impact of nonlocality and hot carriers on the third order nonlinear process are reported in Fig. 5. We report a full comparison of TH conversion efficiency as a function of pump wavelength and angle of incidence in the local [Fig. 5(a)], nonlocal [Fig. 5(b)] and nonlocal – with hot electron contributions [Fig. 5(c)] regimes. In Fig. 5(d) we show three sections of the maps taken at the angle of incidence where maximum conversion efficiency is registered. The competing effects of nonlocal and hot carrier contributions are quite noticeable here as well. By abandoning the local model to include nonlocal contributions [$\chi _{hot,FF/TH}^{(3)} = 0$, compare Figs. 5(a) and (b)] we record a significant blue-shift (∼ 70 nm) of the third harmonic peak, and a slight shift in the angular response (maximum conversion efficiency moves from ϑ_i= 65° to ϑ_i= 69°). Furthermore, since the introduction of nonlocal effects hampers electric field enhancement [see Figs. 4(a) and (b) for FIE calculations], overall conversion efficiency decreases by one order of magnitude, to 10⁻⁸. On the other hand, the introduction of hot electron contributions [Fig. 5(c)] pushes the spectral resonance toward its original “local” spectral position and causes distortions in the angular response that reflects what happens to the FIE [Fig. 4(c)]. However, while the FIE in the nonlocal scenario that includes the hot electron contributions [Fig. 4(c)] is comparable with the FIE in the absence of hot electron contributions [Fig. 4(b)], TH conversion efficiency is now of order 10⁻⁵, surpassing the simulated local and nonlocal scenarios [Figs. 5(a) and (b)]. One way to understand the large change in conversion efficiency is as follows. There are two sources of THG: the ENZ condition, which is resonant and dominates near the crossing point, and a non-resonant contribution due to bound electrons, which dominates far from the ENZ condition. Once the hot electron nonlinearity is triggered, newly generated TH photon flood the system [45] and contribute to increase conversion efficiency. The analysis presented in Fig. 5 therefore confirms that treating the multilayer as an effective medium or failing to consider these additional contributions in the nonlinear dynamics will return an inadequate prediction of the nanostructure linear and nonlinear response. The competition and interaction of nonlocal effects and hot electron contributions lead in fact to non-trivial results that cannot be neglected when either designing or interpreting experimental results.

In Fig. 5(e) we show THG conversion efficiency for three different pulse durations, calculated using the hydrodynamic model discussed above. Longer pulses can resolve the resonance, with local fields reaching larger values and thus yielding improved conversion efficiencies and further resonance shifts. We note that the degree of distortion induced in the nonlocal scenario that includes hot electron contributions [Fig. 5(c) and green curve, triangle markers in Fig. 5(d)] strongly depends on pump irradiance and can lead to bistability at relatively low pump intensities. The onset of bistability is strongly suggested in Fig. 5(e), for 150fs incident pulses.

 

Fig. 6. (a) Pump absorption spectra calculated for low pump irradiance value (1 KW/cm²) in the local (red line, square markers) and nonlocal (blue line, circle markers) approximation for a 60 fs incident pulse and nonlocal approximation for a 120 fs incident pulse (black line, triangle markers); (b) Pump absorption spectra calculated for high pump irradiance value (1 GW/cm²) in the local (black line, triangle markers), nonlocal (blue line, circle markers) and local scenario with hot electron contributions (red line, square markers) with all pulses 60fs long. Angle of incidence is ϑ_i=65° for all plots. All curves have been calculated by means of the microscopic, hydrodynamic approach.

Download Full Size | PDF

Finally, we monitor linear and nonlinear absorption of the multilayer stack in the three regimes. Results are shown in Figs. 6(a) and (b), respectively. Figure 6(a) shows the linear pump absorption profiles calculated in the local and nonlocal approximations, for 60 fs and 120 fs incident pulses and low pump irradiance value (1 KW/cm²) so as not to trigger nonlinear effects. Remarkably, the nonlocal terms, whose strength is predetermined by the Fermi energy and effective free electron mass, blueshift the resonance by nearly 200 nm. The effect of pulse duration may also be ascertained from the figure: the 120 fs pulse couples more strongly with the structures, leading to larger local fields, stronger light-matter coupling, and more pump absorption. In Fig. 6(b) we show the effect of increasing input power density to 1 GW/cm², and simultaneously triggering hot electron contributions. We monitor pump absorption within the multilayer in the following scenarios: local, with (red line, square markers) and without hot electron contributions (black line, triangle markers); and nonlocal, with the hot electron nonlinearity included (blue line, circle markers) [see Fig. 6(b)]. The inclusion of the hot electron contribution in Fig. 6(b) mitigates the blueshift displayed in the linear regime [Fig. 6(a)], broadening, red-shifting, and smearing the resonance as the plasma frequency decreases (red line, square markers). The significance and complexity of this dynamical, time-dependent interplay cannot be overstated, and it impacts the fundamental and its harmonics. All absorption curves in Fig. 6 were calculated assuming ϑ_i=65°.

As mentioned earlier, and as previously reported [38], the second order nonlinear response in these structures is also significant. Additional effects occur at the metal/doped-oxide interfaces due to the large field confinement and discontinuities in the free electron density [47,52]. In Fig. 7 we limit ourselves to depicting the predicted, reflected SHG conversion efficiency spectrum for ϑ_i= 65° incident angle, and different pulse durations (50 fs pulses – blue, solid line, 100 fs pulses – red, dashed line, 150 fs pulses – green, dotted line). SHG undergoes shifts similar to THG as pulse duration is increased, but the maximum conversion efficiency is not altered significantly because SHG is more sensitive to field discontinuities than field localization. Finally, we compare the SH signal generated by the anisotropic, metal/conducting-oxide multilayer stack with the SH signal generated by a single, 20 nm-thick Dy:CdO layer set in a Kretschmann configuration [20], assuming a MgO prism and a 150 fs pulse incident at ϑ_i= 36° (black line, diamond markers). The figure suggests that at least one order of magnitude and 100nm separate the efficiency of the SH signals generated under the two different circumstances.

 

Fig. 7. (a) SHG conversion efficiency for the multilayer in Fig. 3 with ϑ_i= 65°, for 50 fs pulse (blue, solid line), 100fs pulse (red, dashed line) and 150 fs (green, dotted line); Curves are calculated in the nonlocal + he scenario (b) SHG conversion efficiency for a 20 nm thick Dy:CdO layer in Kretschmann configuration (black line, diamond markers) assuming ϑ_i= 36° and 150 fs pulse. Curves in (a) have been calculated by means of the microscopic, hydrodynamic approach.

Download Full Size | PDF

4. Conclusions

In conclusion, we have attempted to clarify the role anisotropy plays in third order nonlinear processes. We demonstrated that a high degree of anisotropy can enhance nonlinear processes thanks to the high local field values that can be achieved. We also found that the sign of the anisotropy plays a different role in FIE and TH conversion efficiency when thicker slabs are investigated, suggesting that it may possible to exploit metals to achieve a large degree of anisotropy. Finally, we demonstrated that while an effective medium approach may be useful as a general guideline to design appropriate nanostructures, care should be exercised when non-homogeneous field distributions are achieved, as in the proposed multilayer. The introduction of nonlocal effects and hot electron contributions in fact reveals how these two effects compete and interact, leading not only to significant spectral and angular shifts, but also to the possible suppression/enhancement of the generated signal that depends on pump excitation conditions.

Appendix A

In Section 2 we showed that strongly anisotropic structures can improve FIE and, therefore, boost harmonic generation. A multilayer structure composed of a conducting oxide and a metal provides a good combination of permittivity values that allows to achieve high anisotropy while preserving a zero-crossing in the real part of the effective permittivity in the longitudinal direction. Although the inclusion of metals in the nanostructure might be necessary to achieve simultaneously a negative transverse permittivity and high anisotropy, the choice of materials that can be included in the multilayer stack is wide and depends also on the operating frequency range one aims to exploit. Here we report and compare the longitudinal and transverse effective permittivities of four multilayer stacks composed of different metals and different highly doped conducting oxides. More specifically Fig. 8 shows the real and imaginary parts of (a) the longitudinal ɛ_z and (b) the transverse effective permittivities ɛ_x, ɛ_y of a multilayer stack composed of Au and Dy:CdO (i.e., the same structure discussed in Section 3). Figures 9, 10 and 11 show the same quantities for a Ag/Dy:CdO, Au/ITO and Ag/ITO multilayers, respectively. All effective permittivities are calculated assuming an infinite number of layers each 10nm thick. The plots show that the zero-crossing point for the longitudinal component of the dielectric permittivity can be chosen by including a material having its plasma frequency in another frequency range (e.g., Dy:CdO shows its plasma frequency around 2.1 µm, while ITO shows its plasma frequency around 1.2 µm). Metals, on the other hand, impact the overall damping and the value of |Δɛ| at the operating frequency. For example, by including silver in the multilayer we generally obtain a lower |Δɛ| with respect to structures that include gold [compare Fig. 8(b) with Fig. 9(b) and Fig. 10(b) with Fig. 11(b)].

 

Fig. 8. Real and imaginary part for the (a) longitudinal and (b) transverse effective permittivities for a Au/Dy:CdO multilayer. The thickness of each layer is assumed equal to 10 nm.

Download Full Size | PDF

 

Fig. 9. Real and imaginary part for the (a) longitudinal and (b) transverse effective permittivities for a Ag/Dy:CdO multilayer. The thickness of each layer is assumed equal to 10 nm.

Download Full Size | PDF

 

Fig. 10. Real and imaginary part for the (a) longitudinal and (b) transverse effective permittivities for a Au/ITO multilayer. The thickness of each layer is assumed equal to 10 nm.

Download Full Size | PDF

 

Fig. 11. Real and imaginary part for the (a) longitudinal and (b) transverse effective permittivities for a Ag/ITO multilayer. The thickness of each layer is assumed equal to 10 nm.

Download Full Size | PDF

Appendix B

To solve the electromagnetic problem at the fundamental and harmonic frequencies we use two different, independent methods. Here we briefly describe the main characteristics of these implementations and the set of equations that have been used to solve the electromagnetic problem at the fundamental, second and third harmonic frequencies.

Hydrodynamic approach

The hydrodynamic approach is implemented by means of a fast Fourier transform pulse propagation method [47]. In this context Dy:CdO is fully described by a combined Lorentz-Drude model as follows:

(B1)$${\varepsilon _{Dy:CdO}}(\omega ) = 1 - \frac{{\omega _{p,b,Dy:CdO}^2}}{{{\omega ^2} - \omega _{0,b,Dy:CdO}^2 + i{\gamma _{b,Dy:CdO}}\omega }} - \frac{{\omega _{p,Dy:CdO}^2}}{{{\omega ^2} + i{\gamma _{f,Dy:CdO}}\omega }}. $$

Here, $\omega _{p,b,Dy:CdO}^{}$ is the bound electron plasma frequency, defined similarly to the free electron counterpart defined in Section 3; $\omega _{0,b,Dy:CdO}^{}$ is the resonance frequency; and ${\gamma _{b,Dy:CdO}}$ the bound electron damping coefficient. The real part of the total dielectric constant of Dy:CdO crosses zero in the infrared. The Lorentzian portion is peaked in the UV range. The dielectric constant of the metal is similar to the function depicted in Eq. (B1), except that it contains parameters appropriate for gold and an additional Lorentzian function to describe one free electron and two bound electron species. The dielectric functions are never specified in the model. Instead, we fit the measured, local dielectric constants to determine damping coefficients, effective masses, and densities to be inserted in dynamical equations of motion. The resulting, coupled material equations of motion that describe separate polarization components produced by free and bound electron may be written as follows:

(B2)$$\begin{aligned} \ddot{\textbf{P}}_f^{} &+ {\gamma _f}{{\dot{\textbf{P}}}_f} = \frac{{{n_{0,f}}{e^2}}}{{m_f^\ast ({T_e})}}\textbf{E} - \frac{{{e^{}}}}{{m_f^\ast ({T_e})}}\textbf{E}({\nabla \bullet {\textbf{P}_f}} )+ \frac{{{e^{}}}}{{m_f^\ast ({T_e})}}{{\dot{\textbf{P}}}_f} \times \textbf{H}\\ &+ \frac{{3{E_F}}}{{5m_f^\ast ({T_e})}}({\nabla ({\nabla \bullet {\textbf{P}_f}} )+ {\nabla^2}{\textbf{P}_f}} )- \frac{1}{{{n_{0,f}}e}}[{({\nabla \bullet {{\dot{\textbf{P}}}_f}} ){{\dot{\textbf{P}}}_f} + ({{{\dot{\textbf{P}}}_f} \bullet \nabla } ){{\dot{\textbf{P}}}_f}} ]\end{aligned}, $$

(B3)$${\ddot{\textbf{P}}_b} + {\gamma _b}{\dot{\textbf{P}}_b} + \omega _{0,b}^2{\textbf{P}_b} + {\textbf{P}_{b,NL}} = \frac{{{n_{0,b}}{e^2}}}{{m_b^\ast }}\textbf{E} + \frac{e}{{m_b^\ast }}{(}{\textbf{P}_b} \bullet \nabla {)}\textbf{E} + \frac{e}{{m_b^\ast c}}{\dot{\textbf{P}}_b} \times \textbf{H}. $$

Equations (B2) and (B3) describe the dynamics of free and bound electrons, respectively; T_e is the free electron temperature; ${\textbf{P}_{b,NL}} = \alpha {\textbf{P}_b}{\textbf{P}_b} - \beta ({\textbf{P}_b} \bullet {\textbf{P}_b}){\textbf{P}_b} + \ldots .$ is the bulk crystal’s nonlinear polarization; $m_b^\ast $ is the bound electron mass; ${n_{0,b}}$ is the bound electron density. α = 0 for centrosymmetric media like Dy:CdO or noble metals. The third order nonlinearity $\beta ({\textbf{P}_b} \bullet {\textbf{P}_b}){\textbf{P}_b}$ is assumed to be isotropic. We note that for up to third harmonic generation, each of Eqs. (B2) and (B3) splits into three separate equations that describe envelope function without ever introducing the slowly envelope approximation. Equations (B2) and (B3) preserve linear and nonlinear dispersions and can describe harmonic generation from any type of surface. The free electron component in Dy:CdO display an effective mass $m_f^\ast $ that is a function of temperature. If the electron temperatures do not exceed a few thousand degrees Kelvin, then one may write:

(B4)$$m_f^\ast ({T_e}) \approx m_0^\ast{+} \alpha {K_B}{T_e} = m_0^\ast{+} \alpha {K_B}\Lambda \int {\int {\textbf{J} \bullet \textbf{E}\,\textrm{d}{\textbf{r}^3}} } dt, $$

where K_B is Boltzmann’s constant; α and Λ are constants of proportionality; $m_0^\ast $ is the electron’s rest mass for no applied field; $\Lambda \int {\int {\textbf{J} \bullet \textbf{E}\,\textrm{d}{\textbf{r}^3}} } dt$ represents absorption; and $\textbf{J} = {\dot{\textbf{P}}_f}$ is the free current density. For practical purposes, Eq. (B4) may be simplified so that $\textbf{J} = {\sigma _0}\textbf{E}$, where ${\sigma _0}$ is a constant to be determined. Then, using a simplified temperature dependent expression for the effective mass, the plasma frequency takes the following form:

(B5)$$\frac{{{n_{0,f}}{e^2}}}{{m_f^\ast ({T_e})}}\textbf{E} \approx \,\,\frac{{{n_{0,f}}{e^2}}}{{m_0^\ast }}\textbf{E} - \,\tilde{\Lambda }({\textbf{E} \bullet \textbf{E}} )\textbf{E} + \,\,{\tilde{\Lambda }^2}{({\textbf{E} \bullet \textbf{E}} )^2}\textbf{E} - {\tilde{\Lambda }^3}{({\textbf{E} \bullet \textbf{E}} )^3}\textbf{E} + \ldots .$$

The coefficient $\tilde{\Lambda }$ thus determines the spatio-temporal dynamics of the redshift impressed upon the plasma frequency. In the expansion we retain terms up to the 7^th power of the field because the large local fields and slow convergence require one to do so, with consequential lowest order contributions to odd harmonics (in our case FF and TH). As a final step, the total polarization is written as the vector sum of free and bound electron contributions, ${\textbf{P}_{Total}} = {\textbf{P}_f} + {\textbf{P}_b}$, and is inserted into Maxwell’s equations, which take the following form:

(B6)$$\begin{aligned} &\nabla \times \textbf{E} ={-} \frac{{\partial \textbf{H}}}{{\partial t}}\\ &\nabla \times \textbf{H} = \frac{{\partial \textbf{E}}}{{\partial t}} + \,\frac{{\partial {\textbf{P}_{Total}}}}{{\partial t}} \end{aligned}. $$

Finite element method

The solution for the fundamental frequency field is found numerically by solving Maxwell’s equations as in Eqs. (B6). We then use the solution obtained at the fundamental frequency to evaluate the third order nonlinear response. For this purpose, we use the same finite-element solver and write the nonlinear polarization P_Total as the sum of the linear polarization P_L and the nonlinear polarization P_NL, where the time derivative of P_NL (J_NL) in the local regime is written as Eqs. (1) and (2) for the fundamental and third harmonic frequency, respectively. The expression of J_NL then is modified as described in Eqs. (3) and (4) when nonlocal and hot electrons contributions are considered. The solution for second harmonic frequency is also found numerically by solving the weak form of the inhomogeneous Helmholtz equation. To account for second order nonlinear effects we evaluate the effective second-order response that arises from symmetry breaking at the interface, modeled by means of second harmonic (SH) current density sources that takes into account volume and surface contributions from magnetic dipoles (Lorentz force) and convective nonlinear sources [47,48]. In particular, SH current density sources are introduced as the superposition of two terms: a volume current, ${\textbf{J}_{\textrm{vol}}}$, and a surface current, ${\textbf{J}_{\textrm{surf}}}$. These currents can then be linked to the FF electric field and to the free electron hydrodynamic parameters as follows [49,53,54]:

(B7)$$\hat{\textbf{n}} \cdot {\textbf{J}_{\textrm{surf}}} = i\frac{{{n_{0,f}}{e^3}}}{{2m_{0,f}^{{\ast} 2}}}\frac{{3 + {\varepsilon _{\textrm{FF}}}}}{{{{({\omega + \textrm{i}{\gamma_f}} )}^2}({2\omega + \textrm{i}{\gamma_f}} )}}E_{\textrm{FF},\textrm{n}}^2, $$

(B8)$$\hat{\textbf{t}} \cdot {\textbf{J}_{\textrm{surf}}} = i\frac{{2{n_{0,f}}{e^3}}}{{m_{0,f}^{{\ast} 2}}}\frac{1}{{{{({\omega + i{\gamma_f}} )}^2}({2\omega + i{\gamma_f}} )}}E_{\textrm{FF},\textrm{n}}^{}E_{\textrm{FF},\textrm{t}} , $$

(B9)$$\begin{aligned}{\textbf{J}_{\textrm{vol}}} &= \frac{{{n_{0,f}}{e^3}}}{{m_{0,f}^{{\ast} 2}}}\frac{1}{{\omega ({\omega + i{\gamma_f}} )({2\omega + i{\gamma_f}} )}}\\ &\left[ {\frac{{{\gamma_f}}}{{\omega + i{\gamma_f}}}({{\textbf{E}_{\textrm{FF}}} \cdot \nabla } ){\textbf{E}_{\textrm{FF}}} - \frac{i}{2}\nabla ({{\textbf{E}_{\textrm{FF}}} \cdot {\textbf{E}_{\textrm{FF}}}} )} \right] \end{aligned}, $$

where ${n_{0,f}} = {{{\varepsilon _0}m_{0,f}^\ast \omega _{p,f}^2} / {{e^2}}}$ is the free electrons density, the effective electron mass is assumed to be $m_{0,f}^\ast $= $m_e^{}$, e is the elementary charge. $\omega _{p,f}^{}$ and ${\gamma _f}$ are the free-electrons plasma frequency and the electron gas collision frequency, respectively. ${{\varepsilon}_{\textrm{FF}}}$ is the relative permittivity at the FF, $\omega$ is the angular frequency of the FF field, ${\textbf{E}_{\textrm{FF}}}$ is the FF electric field phasor, and $\hat{\textbf{n}}$ and $\hat{\textbf{t}}$ are unit vectors pointing in directions outward normal and tangential to the slab surface, respectively. Moreover, $E_{\textrm{FF},\textrm{n}}^{}$ and $E_{\textrm{FF},\textrm{t}}^{}$ are the normal and tangential components of the FF electric field in the local boundary coordinate system defined by $\hat{\textbf{n}}$ and $\hat{\textbf{t}}$, respectively, and are evaluated inside the slab region.

Funding

Army Research Laboratory (W911NF-20-2-0078).

Acknowledgments

Research was sponsored by the Army Research Laboratory and was accomplished under Cooperative Agreement Number W911NF-20-2-0078.

Disclosures

The authors declare no conflicts of interest.

References

1. R. W. Ziolkowski, “Propagation in and scattering from a matched metamaterial having a zero index of refraction,” Phys. Rev. E 70(4), 046608 (2004). [CrossRef]  

2. M. Silveirinha and N. Engheta, “Tunneling of Electromagnetic Energy through Subwavelength Channels and Bends using ɛ-Near-Zero Materials,” Phys. Rev. Lett. 97(15), 157403 (2006). [CrossRef]  

3. S. Campione, D. de Ceglia, M. A. Vincenti, M. Scalora, and F. Capolino, “Electric field enhancement in ɛ-near-zero slabs under TM-polarized oblique incidence,” Phys. Rev. B 87(3), 035120 (2013). [CrossRef]  

4. S. Enoch, G. Tayeb, P. Sabouroux, N. Guérin, and P. Vincent, “A Metamaterial for Directive Emission,” Phys. Rev. Lett. 89(21), 213902 (2002). [CrossRef]  

5. G. Lovat, P. Burghignoli, F. Capolino, D. R. Jackson, and D. R. Wilton, “Analysis of directive radiation from a line source in a metamaterial slab with low permittivity,” Antennas and Propagation,” IEEE Trans. Antennas Propag. 54(3), 1017–1030 (2006). [CrossRef]  

6. D. Rocco, M. A. Vincenti, and C. De Angelis, “Boosting Second Harmonic Radiation from AlGaAs Nanoantennas with Epsilon-Near-Zero Materials,” Appl. Sci. 8(11), 2212 (2018). [CrossRef]  

7. A. Alù, M. G. Silveirinha, A. Salandrino, and N. Engheta, “Epsilon-near-zero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern,” Phys. Rev. B 75(15), 155410 (2007). [CrossRef]  

8. P. Ginzburg, F. J. R. Fortuño, G. A. Wurtz, W. Dickson, A. Murphy, F. Morgan, R. J. Pollard, I. Iorsh, A. Atrashchenko, P. A. Belov, Y. S. Kivshar, A. Nevet, G. Ankonina, M. Orenstein, and A. V. Zayats, “Manipulating polarization of light with ultrathin epsilon-near-zero metamaterials,” Opt. Express 21(12), 14907–14917 (2013). [CrossRef]  

9. J. Wu, B. A. Malomed, H. Y. Fu, and Q. Li, “Self-interaction of ultrashort pulses in an epsilon-near-zero nonlinear material at the telecom wavelength,” Opt. Express 27(26), 37298–37307 (2019). [CrossRef]  

10. A. Ciattoni, C. Rizza, A. Marini, A. D. Falco, D. Faccio, and M. Scalora, “Enhanced nonlinear effects in pulse propagation through epsilon-near-zero media,” Laser Photonics Rev. 10(3), 517–525 (2016). [CrossRef]  

11. T. Zhai and X. Zhang, “Epsilon-near-zero metamaterials for tailoring ultrashort pulses,” Appl. Phys. B 113(2), 185–189 (2013). [CrossRef]  

12. C. Rizza, A. Ciattoni, and E. Palange, “Two-peaked and flat-top perfect bright solitons in nonlinear metamaterials with epsilon near zero,” Phys. Rev. A 83(5), 053805 (2011). [CrossRef]  

13. A. Ciattoni, C. Rizza, and E. Palange, “Multistability at arbitrary low optical intensities in a metal-dielectric layered structure,” Opt. Express 19(1), 283–288 (2011). [CrossRef]  

14. A. Ciattoni, C. Rizza, and E. Palange, “Extreme nonlinear electrodynamics in metamaterials with very small linear dielectric permittivity,” Phys. Rev. A 81(4), 043839 (2010). [CrossRef]  

15. A. Ciattoni, C. Rizza, and E. Palange, “Transmissivity directional hysteresis of a nonlinear metamaterial slab with very small linear permittivity,” Opt. Lett. 35(13), 2130–2132 (2010). [CrossRef]  

16. D. de Ceglia, S. Campione, M. A. Vincenti, F. Capolino, and M. Scalora, “Low-damping epsilon-near-zero slabs: Nonlinear and nonlocal optical properties,” Phys. Rev. B 87(15), 155140 (2013). [CrossRef]  

17. O. Reshef, I. De Leon, M. Z. Alam, and R. W. Boyd, “Nonlinear optical effects in epsilon-near-zero media,” Nat. Rev. Mater. 4(8), 535–551 (2019). [CrossRef]  

18. M. Z. Alam, S. A. Schulz, J. Upham, I. De Leon, and R. W. Boyd, “Large optical nonlinearity of nanoantennas coupled to an epsilon-near-zero material,” Nat. Photonics 12(2), 79–83 (2018). [CrossRef]  

19. M. Z. Alam, I. De Leon, and R. W. Boyd, “Large optical nonlinearity of indium tin oxide in its epsilon-near-zero region,” Science 352(6287), 795–797 (2016). [CrossRef]  

20. T. S. Luk, D. de Ceglia, S. Liu, G. A. Keeler, R. P. Prasankumar, M. A. Vincenti, M. Scalora, M. B. Sinclair, and S. Campione, “Enhanced third harmonic generation from the epsilon-near-zero modes of ultrathin films,” Appl. Phys. Lett. 106(15), 151103 (2015). [CrossRef]  

21. M. A. Vincenti, S. Campione, D. de Ceglia, F. Capolino, and M. Scalora, “Gain-assisted harmonic generation in near-zero permittivity metamaterials made of plasmonic nanoshells,” New J. Phys. 14(10), 103016 (2012). [CrossRef]  

22. M. A. Vincenti, D. de Ceglia, A. Ciattoni, and M. Scalora, “Singularity-driven second- and third-harmonic generation at ɛ-near-zero crossing points,” Phys. Rev. A 84(6), 063826 (2011). [CrossRef]  

23. A. Ciattoni and E. Spinozzi, “Efficient second-harmonic generation in micrometer-thick slabs with indefinite permittivity,” Phys. Rev. A 85(4), 043806 (2012). [CrossRef]  

24. A. Capretti, Y. Wang, N. Engheta, and L. Dal Negro, “Comparative Study of Second-Harmonic Generation from Epsilon-Near-Zero Indium Tin Oxide and Titanium Nitride Nanolayers Excited in the Near-Infrared Spectral Range,” ACS Photonics 2(11), 1584–1591 (2015). [CrossRef]  

25. X. Wen, G. Li, C. Gu, J. Zhao, S. Wang, C. Jiang, S. Palomba, C. Martijn de Sterke, and Q. Xiong, “Doubly Enhanced Second Harmonic Generation through Structural and Epsilon-near-Zero Resonances in TiN Nanostructures,” ACS Photonics 5(6), 2087–2093 (2018). [CrossRef]  

26. V. Roppo, M. A. Vincenti, D. de Ceglia, and M. Scalora, “Deep-subwavelength waveguiding via inhomogeneous second-harmonic generation,” Opt. Lett. 37(15), 3093–3095 (2012). [CrossRef]  

27. M. A. Vincenti, D. de Ceglia, J. W. Haus, and M. Scalora, “Harmonic generation in multiresonant plasma films,” Phys. Rev. A 88(4), 043812 (2013). [CrossRef]  

28. A. Capretti, Y. Wang, N. Engheta, and L. Dal Negro, “Enhanced third-harmonic generation in Si-compatible epsilon-near-zero indium tin oxide nanolayers,” Opt. Lett. 40(7), 1500–1503 (2015). [CrossRef]  

29. D. de Ceglia, M. A. Vincenti, N. Akozbek, M. J. Bloemer, and M. Scalora, “Nested plasmonic resonances: extraordinary enhancement of linear and nonlinear interactions,” Opt. Express 25(4), 3980–3990 (2017). [CrossRef]  

30. E. D. Palik and G. Ghosh, Handbook of optical constants of solids (Academic press, 1998), Vol. 3.

31. I. Liberal and N. Engheta, “Near-zero refractive index photonics,” Nat. Photonics 11(3), 149–158 (2017). [CrossRef]  

32. S. Campione and F. Capolino, “Composite material made of plasmonic nanoshells with quantum dot cores: loss-compensation and epsilon-near-zero physical properties,” Nanotechnology 23(23), 235703 (2012). [CrossRef]  

33. S. Campione, M. Albani, and F. Capolino, “Complex modes and near-zero permittivity in 3D arrays of plasmonic nanoshells: loss compensation using gain [Invited],” Opt. Mater. Express 1(6), 1077–1089 (2011). [CrossRef]  

34. M. A. Vincenti, D. de Ceglia, and M. Scalora, “Nonlinear dynamics in low permittivity media: the impact of losses,” Opt. Express 21(24), 29949–29954 (2013). [CrossRef]  

35. M. Kamandi, C. Guclu, T. S. Luk, G. T. Wang, and F. Capolino, “Giant field enhancement in longitudinal epsilon-near-zero films,” Phys. Rev. B 95(16), 161105 (2017). [CrossRef]  

36. J. Luo, P. Xu, H. Chen, B. Hou, L. Gao, and Y. Lai, “Realizing almost perfect bending waveguides with anisotropic epsilon-near-zero metamaterials,” Appl. Phys. Lett. 100(22), 221903 (2012). [CrossRef]  

37. B. Edwards, A. Alù, M. G. Silveirinha, and N. Engheta, “Reflectionless sharp bends and corners in waveguides using epsilon-near-zero effects,” J. Appl. Phys. 105(4), 044905 (2009). [CrossRef]  

38. L. V. Alekseyev, E. E. Narimanov, T. Tumkur, H. Li, Y. A. Barnakov, and M. A. Noginov, “Uniaxial epsilon-near-zero metamaterial for angular filtering and polarization control,” Appl. Phys. Lett. 97(13), 131107 (2010). [CrossRef]  

39. S. Feng and K. Halterman, “Coherent perfect absorption in epsilon-near-zero metamaterials,” Phys. Rev. B 86(16), 165103 (2012). [CrossRef]  

40. K. Halterman, S. Feng, and V. C. Nguyen, “Controlled leaky wave radiation from anisotropic epsilon near zero metamaterials,” Phys. Rev. B 84(7), 075162 (2011). [CrossRef]  

41. M. A. Vincenti, M. Kamandi, D. de Ceglia, C. Guclu, M. Scalora, and F. Capolino, “Second-harmonic generation in longitudinal epsilon-near-zero materials,” Phys. Rev. B 96(4), 045438 (2017). [CrossRef]  

42. E. Sachet, C. T. Shelton, J. S. Harris, B. E. Gaddy, D. L. Irving, S. Curtarolo, B. F. Donovan, P. E. Hopkins, P. A. Sharma, A. L. Sharma, J. Ihlefeld, S. Franzen, and J.-P. Maria, “Dysprosium-doped cadmium oxide as a gateway material for mid-infrared plasmonics,” Nat. Mater. 14(4), 414–420 (2015). [CrossRef]  

43. M. Vincenti, D. de Ceglia, V. Roppo, and M. Scalora, “Harmonic generation in metallic, GaAs-filled nanocavities in the enhanced transmission regime at visible and UV wavelengths,” Opt. Express 19(3), 2064–2078 (2011). [CrossRef]  

44. L. Rodríguez-Suné, M. Scalora, A. S. Johnson, C. Cojocaru, N. Akozbek, Z. J. Coppens, D. Perez-Salinas, S. Wall, and J. Trull, “Study of second and third harmonic generation from an indium tin oxide nanolayer: Influence of nonlocal effects and hot electrons,” APL Photonics 5(1), 010801 (2020). [CrossRef]  

45. M. Scalora, J. Trull, D. de Ceglia, M. A. Vincenti, N. Akozbek, Z. Coppens, L. Rodríguez-Suné, and C. Cojocaru, “Electrodynamics of conductive oxides: Intensity-dependent anisotropy, reconstruction of the effective dielectric constant, and harmonic generation,” Phys. Rev. A 101(5), 053828 (2020). [CrossRef]  

46. M. Scalora, M. Vincenti, D. de Ceglia, N. Akozbek, V. Roppo, M. Bloemer, and J. Haus, “Dynamical model of harmonic generation in centrosymmetric semiconductors at visible and UV wavelengths,” Phys. Rev. A 85(5), 053809 (2012). [CrossRef]  

47. M. Scalora, M. A. Vincenti, D. de Ceglia, V. Roppo, M. Centini, N. Akozbek, and M. J. Bloemer, “Second- and third-harmonic generation in metal-based structures,” Phys. Rev. A 82(4), 043828 (2010). [CrossRef]  

48. J. E. Sipe, V. C. Y. So, M. Fukui, and G. I. Stegeman, “Analysis of second-harmonic generation at metal surfaces,” Phys. Rev. B 21(10), 4389–4402 (1980). [CrossRef]  

49. M. A. Vincenti, S. Campione, D. de Ceglia, F. Capolino, and M. Scalora, “Gain-assisted harmonic generation in near-zero permittivity metamaterials made of plasmonic nanoshells,” New J. Phys. 14(10), 103016 (2012). [CrossRef]  

50. D. de Ceglia, M. Scalora, M. A. Vincenti, S. Campione, K. Kelley, E. L. Runnerstrom, J.-P. Maria, G. A. Keeler, and T. S. Luk, “Viscoelastic optical nonlocality of low-loss epsilon-near-zero nanofilms,” Sci. Rep. 8(1), 9335 (2018). [CrossRef]  

51. J. F. García de Abajo, “Nonlocal Effects in the Plasmons of Strongly Interacting Nanoparticles, Dimers, and Waveguides,” J. Phys. Chem. C 112(46), 17983–17987 (2008). [CrossRef]  

52. M. Scalora, M. A. Vincenti, D. de Ceglia, N. Akozbek, M. J. Bloemer, C. De Angelis, J. W. Haus, R. Vilaseca, J. Trull, and C. Cojocaru, “Harmonic generation from metal-oxide and metal-metal boundaries,” Phys. Rev. A 98(2), 023837 (2018). [CrossRef]  

53. D. Maystre, M. Neviere, and R. Reinisch, “Nonlinear polarisation inside metals: A mathematical study of the free-electron model,” Appl. Phys. A: Mater. Sci. Process. 39(2), 115–121 (1986). [CrossRef]  

54. A. Benedetti, M. Centini, C. Sibilia, and M. Bertolotti, “Engineering the second harmonic generation pattern from coupled gold nanowires,” J. Opt. Soc. Am. B 27(3), 408–416 (2010). [CrossRef]