## Abstract

Ultrashort pulse propagation at the epsilon-near-zero spectral point is numerically investigated using the finite difference time-domain technique for the pump-probe experiment. Free carriers’ population dynamics in the conduction band for large intensities of the pump pulse and the transient response for rapidly varying pulses in two-level media are calculated. The auxiliary differential equation finite-difference time domain method was used to numerically investigate ultrashort probe pulse propagation in 300 nm of the AZO/ZnO metamaterial. Results show a dramatic change in shape for the probe pulse modulated using pump pulses of various duration (100-500 fs) and amplitude (10^{6}−10^{10} V/m).

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Epsilon-near-zero (ENZ) metamaterials, materials with permittivity |*ε*>| << 1, resulted in various intriguing phenomena [1]. Since ENZ condition requires the real permittivity to turn zero, a very large optical loss due to large imaginary permittivity can limit its effect. If the optical loss in the ENZ regime is minimized, achieving near-zero-index (NZI) materials with not only |*ε*| << 1 but also refractive index |*n*| << 1 is possible [2]. ENZ materials with a low optical loss have a very low group velocity of electromagnetic wave propagation [3]. This phenomenon of low group velocity bears some analogy to the manipulation of the speed of light using various optical systems [4,5]. The low group velocity leads to a modified light propagation effect which could enhance the frequency conversion efficiency in the ENZ material. The advantages of ENZ materials, such as the availability of the high local fields and the fact that frequency mixing could be done without phase-matching, allows for increased efficiency of second [6] and third harmonic generation [7].

There are many approaches to designing ENZ [1] and NZI [2] metamaterials. Transparent conducting oxides (TCO) are the most widely studied ENZ materials due to the broad tunability and low losses at the ENZ spectral point [8,9]. A wide range of established material fabrication methods for AZO/ZnO multilayered structures are available, including atomic layer deposition [10,11], pulsed laser deposition [12,13], and RF sputtering. In this manuscript, we consider an ENZ multilayered structure comprised of AZO/ZnO as it allows for precise engineering of the optical properties at the ENZ spectral point [10–13]. In our previous study, we found that varying fabrication parameters of multilayered AZO/ZnO structures (AZO period and thickness, deposition temperature, total sample thickness) allow for tuning the plasma frequency, *ω _{p}* and the damping frequency,

*γ*[13].

_{p}ENZ metamaterials based on TCO materials provide opportunities for enhancing nonlinear optical interactions due to low group velocity and optical field confinement. For example, the enhancement of the second [14] and third [15] harmonic generation in indium tin oxide (ITO) has been demonstrated. In addition, exceptionally large and ultrafast intensity- dependent index nonlinearities at ENZ spectral point has been demonstrated in various TCOs including ITO [16] and AZO [17]. Beam deflection measurements were used to investigate the nondegenerate nonlinear refraction of ITO at the ENZ wavelength for several different pump and probe polarizations [18].

Our numerical study of ultrashort pulse propagation in multilayered AZO/ZnO showed a “soliton-like” propagation of a 100 fs Gaussian pulse at the ENZ spectral point [19]. This unique propagation regime could be intuitively understood as the chromatic dispersion terms, β_{2} and β_{1}, acting together to prevent a temporal and spectral broadening. However, an ultrashort pulse propagating through distances larger than ∼300 nm experiences a temporal pulse distortion. Adaptive pre-shaping has been suggested [20] to compensate the effect of the enormous second-order dispersion, which could potentially be a bottleneck for applications in super resolution imaging. It is important to note that nonlinear effects become important for ultrashort pulse propagation in ENZ materials due to high optical confinement in ENZ photonic devices.

Several recent experimental and numerical studies have shown frequency shift and spectral distortions as the probe pulse is propagating through a medium with a large time-varying refractive index. A wavelength shift of 15 nm is obtained for near-infrared (1250 nm) probe pulses propagating in the Al-doped-ZnO layer driven by an intense ultrashort pulse at 785 nm [21]. The non-reciprocity of a modified form of Snell’s law was studied, and possible approaches for designing optical isolators using time-varying metasurfaces were proposed [22]. The investigation of the adiabatic frequency shift process in ENZ materials shows that a reduced group velocity enables the adiabatic frequency shift to be achieved in thin samples with relatively low pump energy [23]. A broadband frequency shift (up to 14.9THz) of a light beam using a time-varying subwavelength-thick indium tin oxide sample in its ENZ spectral range has been demonstrated [24].

Optical properties of AZO can be experimentally controlled on a sub-picosecond timescale by a combination of the interband and intraband nonlinear effects, which are driven by two different wavelength pump pulses (∼ 262 nm and ∼ 787 nm) [25]. This pump-probe experiment showed co-existence of both intraband and interband nonlinearities and dynamic control of the bandwidth of a probe optical pulse. The large modulation of the AZO refractive index with terahertz bandwidths is a new photonic capability that may lead to applications in optical metrology and ultrafast optics. Therefore, a numerical approach to investigating ultrashort probe pulse propagation in the presence of high intensity ultrashort pump pulse is needed for future efforts on probe pulse modulation at the ENZ spectral point.

When an electromagnetic wave propagates in a medium, the field induces a time varying dipole moment. The effects of the medium on ultrashort pulse propagation can be modeled by suitably defining the polarization vector of the medium and solving the equation for the macroscopic polarization along with Maxwell’s equations. Then the well-known finite-difference time-domain (FDTD) approach [26] can be used to model dispersive media [27]. This approach, known as the auxiliary differential equation finite-difference time domain (ADE-FDTD) method has also been used to incorporate nonlinear effects [28].

In general, simulating short pulse propagation in the presence of optical pump using the finite difference time-domain (FDTD) technique is a complicated task, as many variables need to be evaluated at the same time. A general multilevel carrier kinetics model is explored to simulate the saturation of absorption in the time-domain [29]. However, the challenge remains to describe ultrashort pulse propagation at the ENZ spectral point in the presence of ultrashort pump pulse.

Here, we propose to use the algorithm which is based on separating the calculations for rate equations and the auxiliary differential equation (ADE) method to investigate the ultrashort pulse propagation at the ENZ spectral point in the presence of time-varying permittivity due to strong UV pump. While the ADE-FDTD method has also been described previously in the context of dispersive media [28], the usage of the ADE-FDTD method for AZO/ZnO provides a unique insight into pump-probe ultrashort pulse dynamics at the ENZ spectral point. In particular, the proposed approach aims to investigate the pulse shape in time domain which is important for future applications in optical metrology and ultrafast optics. The resulting probe pulse shape in time domain can be characterized in the future by using either autocorrelation technique or frequency-resolved optical gating (FROG) method.

## 2. Methods

#### 2.1 Pump-probe approach

The pump–probe experiment is schematically shown in Fig. 1(a). In this case, the ultrashort pump pulse of various durations (100-500 fs) and amplitudes propagated in AZO/ZnO metamaterial. The ultrashort probe pulse (Fig. 1(b)) of various duration (80-100 fs) with ENZ central frequency was propagated with a delay *τ* (0-2 ps). The delay τ is the difference in time between pump and the probe pulses. We assumed that both pump and probe pulse initially have a zero phase in the time domain. The propagation length, *L*, was chosen to be 300 nm for all numerical simulations which is well within fabrication capabilities [11,13]. Normal incidence for both pump and probe beam is assumed in the numerical simulations. The central frequency of the pump pulse was in the UV/blue spectral region (1.14×10^{15} Hz).

The multilayered AZO/ZnO metamaterial in-plane and out-of-plane optical permittivity was modeled using effective medium approximation (EMA). The parameters for the EMA model (AZO/ZnO period thickness: 10 nm, AZO fill factor: 0.4, number of layers: 30) are based on our previous experimental study [13]. To model the optical permittivity, *ε*(*ω*), for AZO, we used the Drude model described below (Eq. (4)). The Adachi model was used as a model of the optical permittivity for ZnO [30].

One of the advantages of using multilayered metamaterial is the possibility of using both TE and TM polarization with a dramatically different response. In particular, the real part of the out-of-plane optical permittivity is *Re* (*ε*_{⊥}) = 0 and the in-plane permittivity is *Re* (*ε*_{∥}) = 2.87 at 1778 nm. Therefore, in the near infrared spectral region, the epsilon-near-zero dispersion is observed only in the case of TM-polarized (p-polarized) incident light because only in this polarization is the electric field component of the incident light able to probe both the in-plane and out-of-plane components of the dielectric function. In our numerical simulations normal incidence was assumed for both pump and probe beam. However, in the experiment, the angle of incidence could be varied to enhance the ENZ phenomena due to out-of-plane permittivity for AZO/ZnO multilayered metamaterial.

To model ultrashort-pulse propagation in AZO/ZnO at the ENZ spectral point in the presence of the UV pump, the permittivity needs to be modeled to include the dynamic change of the carriers. One of the methods to model dynamic carriers’ change is the use of the rate equations which is described below. The carrier dynamics is determined by the rate of change of population inversion. We assume in our simulations that the propagation of the ultrashort pulse is a one-dimensional problem. Generalization to the two-dimensional problem could potentially provide insights into the spatial modulation effects for the probe pulse.

#### 2.2 Rate equations for two-level system

AZO possesses large ultrafast interband and intraband nonlinearities. Carrier dynamics in this material include a variety of the phenomena such as interband transitions, intraband relaxation, and non-radiative processes. AZO exhibits an excitonic peak in the ultraviolet region and a defect related peak in the visible region. The photoluminescence (PL) peak in the UV region is attributed to the radiative recombination of excitons. The wide visible band, centered at ∼500 nm, is due to its intrinsic defects. The un-doped ZnO film has an intense ultraviolet peak together with a broad defect-related PL peak. The defect-related broad peak of ZnO is suppressed in AZO.

In our previous work [31], we developed a complete model for electron dynamics based on the set of rate equations, which includes the formation and recombination of excitons, the creation and recombination of photons, and the electron–hole recombination. Our numerical results [31] for low carrier concentration (∼10^{17} cm^{−3}) show that at low pump intensity (0.01 mJ/cm^{2}) visible emission is dominant in the emission spectrum and, as the pump intensity increases (∼1 mJ/cm^{2}), the UV emission becomes dominant. Ultrafast dynamics study shows that using shorter pump pulses (<1ns) leads to on order of magnitude increase of UV intensity. In this manuscript, high pump intensities (>1-1000 mJ/cm^{2}) and shorter pump pulses were considered. The electrons in the valence band are pumped into the conduction band using an ultrashort UV pulse. Therefore, in the numerical simulations we use the simplified two-level system with only the conduction and valence band, as shown in Fig. 1(b)).

The rate equations describe the time evolution of the two-level system populations under the inﬂuence of applied pump pulse. The conduction and valence bands represent upper and lower energy levels with populations *N _{2}* and

*N*respectively. The electric polarization in two-level transitions can be described by the following equation [32]:

_{1,}**P**is the nonlinear macroscopic polarization,

**E**is the electric field,

*ΔN*

**=**

*N*is the difference in the number of electrons between the upper and lower levels,

_{2}-N_{1}*ω*is the transition frequency, and

_{a}*Δω*is the bandwidth of the transition,

_{a}*κ*= 6

*πε*

_{0}c^{3}/

*t*

_{0}ω^{3}, where

*t*is the decay time from the upper level to the lower level.

_{0}The relationship between population of the conduction and valence band are described by the rate equations. An ideal two-level system can be described as [33]:

*N*is upper level population,

_{2}*N*is lower level population,

_{1}*τ*is upper level lifetime. The effect of doping was included into simulations by limiting the minimum value of the carries in the conduction band at any given time

_{21}*N*to a realistic Al-doping level,

_{2}*N*=3×10

_{2}^{20}cm

^{−3}. Our previous research [11,13] showed that the number of the free carriers defines the position of the ENZ wavelength for investigated AZO/ZnO metamaterial. We found that the ENZ wavelength was

*λ*=1778 nm for this investigation. The transition frequency for AZO,

*ω*, is in the UV spectral region (wavelength 375 nm). The transition bandwidth,

_{a}*Δω*, was estimated to be 50 nm by measuring the width of the UV photoluminescent peak in [34]. The decay time or lifetime of the upper level,

_{a}*τ*, was estimated from our previous work to be 0.5 ns [31]. Nonlinear effects require using high pulse energies. The highest pump pulse electric field values (10

_{21}^{9}−10

^{10}V/m) correspond to the 1-100 TW/cm

^{2}pump pulse intensities. The estimate for the damage threshold of the bulk AZO is ∼19 J/cm

^{2}(corresponds to ∼190 TW/cm

^{2}for a pulse of 100 fs duration) [35]. Therefore, the highest pump pulse amplitude was chosen relatively close to the AZO damage threshold. Although, it is important to note that the optical surface damage threshold of any material in an experiment might depend on the number of parameters including pulse duration, wavelength, and laser repetition rate.

This approach has been used to investigate numerically wave propagation in gain media [32], saturation of absorption in the time-domain [29], and dynamic anisotropic gain media [36].

#### 2.3 Number of free carriers and permittivity

When optical pump pulse energy is greater that the bandgap, the dominant contribution to nonlinearity results from interband transitions. As the pump propagates through material, electrons that reside in the valence band are moved into the conduction band. This results in the increase of the total free-carrier density. The excited carrier recombination at a later time can be described via rate equations which include excitons formation and recombination, photons creation and recombination, and electron-hole recombination [31]. The increase of the total free carriers will affect the plasma frequency, ω_{p}, according to:

*e*is the charge of an electron,

*ε*is the vacuum permittivity, and

_{0}*m**is the reduced mass of an electron for AZO (0.256m

_{e}according to [37]). The change of electrons due to optical pumping are incorporated using n, which is equal to the number of electrons in the conduction band (

*N*) at any given time.

_{2}Since UV pump pulse affects plasma frequency *ω _{p}*, the permittivity

*ε*(

*ω*) will be affected as well. This change to permittivity is described by the Drude model:

*ε*

_{∞}= 3.0 is the permittivity due to bound electrons,

*f*is the pole strength, and

_{p}*γ*= 10

_{p}^{11}Hz is the damping rate. Only pole strength

*f*=0.4 was used. The model parameters were chosen to match experimental data from our previous study [13].

_{0}To simulate ultrashort probe pulse propagation near ENZ spectral point in a Drude media, the auxiliary differential equation (ADE) method was used. This method for implementing FDTD models of dispersive materials uses time-domain auxiliary differential equations linking the polarization and the electric flux density. These equations are numerically solved in time domain in synchronism with Maxwell's curl equations [27]. ADE method uses a phasor polarization current that includes ω_{p} by:

The goal of the ADE technique is to develop a simple time-stepping scheme for phasor polarization current, **J**, which can be updated synchronously with Ampere's law in the time domain. This can be incorporated into the FDTD simulation for any number of Drude poles.

## 3. Results and discussion

Numerical simulations were used to investigate ultrashort probe pulse propagation after one pass through the material in the presence of time dependent permittivity. In our numerical simulations “zero time” is when we start the pump. Home-built FDTD algorithm was created using Python. The time step was defined *dt* = 5×10^{−18} s. The space step was set up *dx* = *dt×c* ≃1.5 nm, where *c* is a speed of light. Stability was estimated and confirmed by the Courant-Friedrichs-Lewy stability limit. First, FDTD simulations were used to investigate the time-dependent optical permittivity ε due to change in free carriers in conduction band for the two-level system at varying intensities of the pump pulse (Sec. 2.1.). Then, the auxiliary diﬀerential equation (ADE) method for ultrashort probe pulse propagation was used for ultrashort pulse modulation investigation for various pump pulse electric fields strengths. (Sec. 3.2). The home built FDTD algorithm was incorporated using Python both for the two-level system and the ADE method.

#### 3.1 Two-level system population and permittivity

Pump pulse propagation in a two-level system was incorporated into the analysis through the polarization response of the medium (Eq. (2)). The home built FDTD algorithm was tested and compared with the known results from [32]. The boundary conditions were specified in the following way. At the *x* = 0 boundary, the incident electric field begins to propagate through the media using the staggered Yee cell. The pulse starts to propagate inside the multilayer at *x*=0. At *x *= *L*, the field energy must dissipate into free space. Computationally, this is done by introducing a reﬂectionless absorbing material, commonly called perfectly matched layers, to keep the energy from reﬂecting back.

The main goal for using FDTD approach for rate equations was to model the transient response of the medium as well as large signal effects. The results of the time-dependent permittivity for various amplitudes of the pump pulse are shown in Fig. 2. We found that electric field amplitudes less than 2×10^{6} V/m do not change optical permittivity significantly (Fig. 2(a)). The results of the optical permittivity change due to pump pulses of the 10^{8} V/m amplitude are comparable within an order of magnitude with experimental results from [38]. As the amplitude increases, the carriers’ population in the conduction band increases and optical permittivity is modulated as a result (Fig. 2(b,c)). Optical permittivity is saturated and resemble a step function as the pump pulse amplitude reaches 2×10^{10} V/m.

The effect of the pump pulse duration on the modulated plasma frequency is shown in Fig. 3. For high pump amplitudes, the plasma frequency is increased over an order of magnitude as pump amplitude is doubled for both pump pulse durations. The time-dependent plasma frequency modulation occurs at around ∼850 fs with various rise time intervals for various pump pulse durations. Shorter pump pulse duration (Fig. 3(a)) results in much faster transition for the plasma frequency at any amplitude compared to longer pump pulse duration (Fig. 3(b)). This dependence could intuitively be understood using Eq. (3) which indicates that plasma frequency is directly related to the number of the electrons in the valence band.

It is important to note that for shorter pump pulse durations and high pump pulse amplitudes the optical permittivity changes over the time period which is comparable with probe pulse duration. This effect of the modulated plasma frequency is included via Eq. (4) and incorporated into FDTD simulations for ultrashort probe pulse propagation described in Sec. 3.2 below.

#### 3.2 Ultrashort pulse propagation through modulated resonance

The ADE-FDTD algorithm for modeling a Drude medium, three-step fully explicit procedure, was used next. The time-domain basis of ADE makes modeling of arbitrary nonlinear dispersive media particularly attractive. The results for the change in carrier concentration (Fig. 2) were incorporated into Drude model (Eq. (4)). The initial probe pulse spectrum is centered at the ENZ wavelength (1778 nm). The time delay *τ* between pump and probe pulse was varied. The time delay *τ* = 850 fs resulted in the strongest probe pulse modulation. As expected, longer and shorter delays did not have significant effect on the pulse since no time-dependent plasma frequency modulation has occurred. It is important to note that the average position in time of the time-dependent plasma frequency modulation will depend on the transition bandwidth *Δω _{a}* for AZO in the UV spectral region and the decay time of the upper level.

The results of the ADE-FDTD algorithm for ultrashort probe pulse propagation (probe pulse duration 100 fs) in 300 nm of the AZO/ZnO metamaterial with time delay between pump and probe pulse *τ* = 850 fs is presented in Fig. 4. The results show that increasing full-width-at-half-maximum (FWHM) pump pulse duration from 100 fs and 500 fs leads to a dramatic change in output pulse behavior. Longer pump pulse (500 fs) leads to much stronger absorption of the probe pulse compared to shorter pump pulse (100 fs) as the pump pulse amplitude is increased. For the highest investigated pump pulse amplitude (limited by damage threshold), the output probe pulse modulated by shorter pump pulse has strongly asymmetric time domain profile (Fig. 4(e)). This asymmetric shape change can intuitively be understood as follows: a high intensity pump pulse produces a time-dependent optical permittivity change (Fig. 2). For a shorter pump pulse, optical permittivity change happens on the time scale comparable to the probe pulse duration which results in a strong influence of the higher order chromatic dispersion terms on the probe pulse. Additional numerical study shows that higher damping frequency (*γ _{p}* = 5×10

^{12}Hz) for AZO/ZnO material does not change the pulse modulation.

The results of Fig. 4 show that a certain threshold of the amplitude of the pump pulse is required to observe modulation of the probe pulse due to time-changing optical permittivity. The threshold behavior is related to the excitation of sufficient number of free carriers which is necessary to observe a noticeable change to the optical permittivity (Fig. 2). We found that the incident UV pump pulse does not experience any significant distortions during propagation, only a slight reduction in amplitude, so it was not included in this discussion. This is expected because the pump pulse spectrum is centered at the UV spectral region which is far from ENZ spectral point (1778 nm), and optical permittivity at the bandgap is not sensitive to changes in the plasma frequency.

It is important to note that the effect of the probe pulse modulation manifests itself in various ways depending on the pump pulse duration. This could intuitively be understood as follows: the shorter pump (100 fs) leads to a high *dω _{p}*/

*dt*over a short period of time which results in a sharp optical permittivity change (Fig. 3(a)). As a result, for the highest investigated pump pulse amplitude (limited by damage threshold), the output probe pulse modulated by a shorter pump pulse has a strongly asymmetric time domain profile (Fig. 5(c)) and slight reduction of the peak amplitude. A longer pump pulse (500 fs) leads to a smaller

*dω*/

_{p}*dt*(Fig. 3(b)) compared to the shorter pump pulse which results in a strong absorption of the probe pulse with pulse shape remaining relatively un-changed in time domain (Fig. 5(d)). In particular, the peak of the electric field amplitude for the output probe pulse modulated by shorter (100 fs) pump pulse is ∼56 times higher compared to the peak of the amplitude for probe pulse modulated with longer pump pulse (500 fs) (Fig. 5(c, d)).

Various fabrication methods [11,13] offer ways to tune spectral position of the ENZ point for out-of-plane permittivity by adjusting the doping concentration and/or thickness of the AZO/ZnO metamaterial during fabrication. Our previous work [20] indicates that spectral tuning of the initial Gaussian pulse with FWHM duration of 100 fs propagating through AZO/ZnO sample will significantly affect output pulse temporal profile. Figure 5 illustrate the results of the spectral tuning of the initial probe pulse for various durations of the pump pulse. The initial probe pulse is tuned in such a way that the central frequency is shifted from below ENZ (positive optical permittivity (Fig. 5(a,b))) to beyond ENZ (negative optical permittivity) spectral point (Fig. 5(e-f)). The time-varying optical permittivity for Fig. 5(a, c, e) and Fig. 5(b, d, f) is defined by plasma frequency in Fig. 3(a) and Fig. 3(b) respectively. The results (Fig. 5) show that the output probe pulse experiences a slight reduction in amplitude for probe pulse with initial central frequency beyond ENZ.

It is important to note that the numerical calculations for the high pump pulse intensities could possibly require expanding the existing model to include the effect of the band filling. While ZnO has a large bandgap of ∼3.3 eV, the Burstein-Moss effect [39,40] results in an increase in bandgap. In our recent work [41], we investigated the interplay between band gap renormalization and band filling (Burstein-Moss effect). We found that in some cases the energy shift due to Burstein-Moss compensates the red-shift from the bandgap renormalization. Many experiments suggest that the Burstein-Moss effect and band gap renormalization for semiconductors with high doping will also strongly depend on the sample’s fabrication method [42,43]. Therefore, the modification of the band filling for the photoluminescence model can only be done in relationship to the particular fabrication method of the Al:ZnO/ZnO which will be a subject of our future research.

## 4. Conclusions

The results of the FDTD numerical study for pump-probe ultrashort pulse modulation in AZO/ZnO metamaterial at the epsilon near zero spectral point show that pump amplitude and FWHM duration are the main factors affecting ultrashort probe pulse propagation. The results of the ADE-FDTD numerical simulations show that a certain threshold of the amplitude of the pump pulse is required to observe modulation in the probe pulse. The effect of the probe pulse modulation varies depending on the pump pulse duration. The shorter pump (∼ 100 fs) pulse leads in a sharp optical permittivity change over a short period of time and results in the strongly asymmetric shape of the output probe pulse for the highest investigated pump pulse amplitude (limited by damage threshold). A longer pump pulse (500 fs) leads to the optical permittivity change over a much longer time interval which results in the much stronger absorption of the probe pulse. We found that the incident UV pump pulse does not experience any significant distortions during propagation. The presence of ENZ properties in near infrared for AZO/ZnO multilayered metamaterial, combined with enhanced interband nonlinearities, could lead to the development of optically controlled telecommunications devices.

## Funding

San Diego State University (San Diego State University UGP Grant (242518)); National Science Foundation (Graduate Research Fellowship Program 1321850).

## Acknowledgment

Priscilla Kelly gratefully acknowledges the financial support from National Science Foundation (NSF) (Graduate Research Fellowship Program 1321850).

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

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

**2. **N. Kinsey, C. DeVault, A. Boltasseva, and V. M. Shalaev, “Near-zero-index materials for photonics,” Nat. Rev. Mater. **4**(12), 742–760 (2019). [CrossRef]

**3. **M. Javani and M. Stockman, “Real and Imaginary Properties of Epsilon-Near-Zero Materials,” Phys. Rev. Lett. **117**(10), 107404 (2016). [CrossRef]

**4. **R. Boyd and D. Gauthier, “"Slow'‘ and “fast” light,” Prog. Opt. **43**, 497–530 (2002). [CrossRef]

**5. **J. Khurgin, “Slow light in various media: a tutorial,” Adv. Opt. Photonics **2**(3), 287–318 (2010). [CrossRef]

**6. **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 (2015). [CrossRef]

**7. **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]

**8. **M. A. Noginov, L. Gu, J. Livenere, G. Zhu, A. K. Pradhan, R. Mundle, M. Bahoura, Yu. A. Barnakov, and V. A. Podolskiy, “Transparent conductive oxides: plasmonic materials for telecom wavelengths,” Appl. Phys. Lett. **99**(2), 021101 (2011). [CrossRef]

**9. **G. V. Naik, J. Kim, and A. Boltasseva, “Oxides and nitrides as alternative plasmonic materials in the optical range,” Opt. Mater. Express **1**(6), 1090 (2011). [CrossRef]

**10. **C. T. Riley, T. A. Kieu, J. S. T. Smalley, S. H. A. Pan, S. J. Kim, K. W. Post, A. Kargar, D. N. Basov, X. Pan, Y. Fainman, D. Wang, and D. J. Sirbuly, “Plasmonic tuning of aluminum doped zinc oxide nanostructures by atomic layer deposition,” Phys. Status Solidi RRL **8**(11), 948–952 (2014). [CrossRef]

**11. **P. Kelly, M. Liu, and L. Kuznetsova, “Designing optical metamaterial with hyperbolic dispersion based on an Al:ZnO/ZnO nano-layered structure using the atomic layer deposition technique,” Appl. Opt. **55**(11), 2993–2997 (2016). [CrossRef]

**12. **G. V. Naik, J. Liu, A. V. Kildishev, V. M. Shalaev, and A. Boltasseva, “Demonstration of Al: ZnO as a plasmonic component for near-infrared metamaterials,” Proc. Natl. Acad. Sci. **109**(23), 8834–8838 (2012). [CrossRef]

**13. **P. Kelly, W. Zhang, M. Liu, and L. Kuznetsova, “Engineering the structural, plasmonic, and optical properties of multilayered aluminum-doped zinc oxide metamaterial grown by pulsed laser deposition,” Appl. Opt. **58**(21), 5681–5686 (2019). [CrossRef]

**14. **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]

**15. **T. S. Luk, D. de Ceglia, S. Liu, G. A. Keeler, R. O. 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]

**16. **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]

**17. **L. Caspani, R. P. M. Kaipurath, M. Clerici, M. Ferrera, T. Roger, J. Kim, N. Kinsey, M. Pietrzyk, A. Di Falco, V. M. Shalaev, A. Boltasseva, and D. Faccio, “Enhanced nonlinear refractive index in epsilon-near-zero materials,” Phys. Rev. Lett. **116**(23), 233901 (2016). [CrossRef]

**18. **D. J. Hagan, S. A. Benis, and E. W. Van Stryland, “Large, ultrafast induced index changes in ITO,” Proc. SPIE **11080**, 110800E (2019). [CrossRef]

**19. **P. Kelly and L. Kuznetsova, “Pulse shaping in the presence of enormous second-order dispersion in Al: ZnO/ZnO epsilon-near-zero metamaterial,” Appl. Phys. B **124**(4), 60 (2018). [CrossRef]

**20. **P. Kelly and L. Kuznetsova, “Adaptive pre-shaping for ultrashort pulse control during propagation in AZO/ZnO multilayered metamaterial at the epsilon-near-zero spectral point,” OSA Continuum **3**(2), 143 (2020). [CrossRef]

**21. **A. M. Shaltout, M. Clerici, N. Kinsey, R. Kaipurath, J. Kim, E. G. Carnemolla, D. Faccio, A. Boltasseva, V. M. Shalaev, and M. Ferrera, “Doppler-Shift Emulation Using Highly Time-Refracting TCO Layer,” in Conference on Lasers and Electro-Optics, OSA Technical Digest (online) (Optical Society of America, 2016), paper FF2D.6.

**22. **A. Shaltout, A. Kildishev, and V. Shalaev, “Time-varying metasurfaces and Lorentz non-reciprocity,” Opt. Mater. Express **5**(11), 2459–2467 (2015). [CrossRef]

**23. **J. B. Khurgin, M. Clerici, V. Bruno, L. Caspani, C. DeVault, J. Kim, A. Shaltout, A. Boltasseva, V. M. Shalaev, M. Ferrera, D. Faccio, and N. Kinsey, “Adiabatic frequency shifting in epsilon-near-zero materials: the role of group velocity,” Optica **7**(3), 226 (2020). [CrossRef]

**24. **Y. Zhou, M. Z. Alam, M. Karimi, J. Upham, O. Reshef, C. Liu, A. E. Willner, and R. W. Boyd, “Broadband frequency translation through time refraction in an epsilon-near-zero material,” Nat. Commun. **11**(1), 2180 (2020). [CrossRef]

**25. **M. Clerici, N. Kinsey, C. DeVault, J. Kim, E. G. Carnemolla, L. Caspani, A. Shaltout, D. Faccio, V. Shalaev, A. Boltasseva, and M. Ferrera, “Controlling hybrid nonlinearities in transparent conducting oxides via two-colour excitation,” Nat. Commun. **8**(1), 15829 (2017). [CrossRef]

**26. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**(3), 302–307 (1966). [CrossRef]

**27. **A. Taflove and S. C. Hagness, * Computational Electrodynamics: The Finite-Difference Time-Domain Method*. (Artech House, Norwood, MA, 2005), 3rd ed.

**28. **P. M. Goorjian and A. Taflove, “Direct time integration of Maxwell’s equations in nonlinear dispersive media for propagation and scattering of femtosecond electromagnetic solitons,” Opt. Lett. **17**(3), 180–182 (1992). [CrossRef]

**29. **S. I. Azzam and A. V. Kildishev, “Time-domain dynamic of saturation of absorption using multilevel atomic systems,” Opt. Mater. Express **8**(12), 3829 (2018). [CrossRef]

**30. **Y. Yoshikawa and S. Adachi, “Optical constants of ZnO,” Jpn. J. Appl. Phys. **36**(Part 1, No. 10), 6237–6243 (1997). [CrossRef]

**31. **B. Campbell, E. Zarate, P. Kelly, and L. Kuznetsova, “Three-level system for numerical modeling of ultraviolet and visible photoluminescence of Aluminum-doped Zinc Oxide,” J. Opt. Soc. Am. B **36**(4), 1017 (2019). [CrossRef]

**32. **A. S. Nagra and R. A. York, “FDTD Analysis of Wave Propagation in Nonlinear Absorbing and Gain Media,” IEEE Trans. Antennas Propag. **46**(3), 334–340 (1998). [CrossRef]

**33. **A. E. Siegman, * Lasers*. (Univ. Sci. Books, 1986).

**34. **Y. Liu and J. Lian, “Optical and electrical properties of aluminum-doped ZnO thin films grown by pulsed laser deposition,” Appl. Surf. Sci. **253**(7), 3727–3730 (2007). [CrossRef]

**35. **Y. Xu, Y. Lu, Y. Zuo, F. Xu, and D. Zuo, “Z-scan measurements of nonlinear refraction and absorption for aluminum-doped zinc oxide thin film,” Appl. Opt. **58**(22), 6112 (2019). [CrossRef]

**36. **A. A. Al-Jabr, D. P. San Roman Alerigi, B. S. Ooi, and M. A. Alsunaidi, “An FDTD algorithm for simulating light propagation in anisotropic dynamic gain media,” Proc. SPIE **9134**, 91341U (2014). [CrossRef]

**37. **M. A. Bodea, G. Sbarcea, G. V. Naik, A. Boltesseva, T. A. Klar, and J. D. Pedarnig, “Negative permittivity of ZnO thin films prepared from aluminum and gallium doped ceramics via pulsed-laser-deposition,” Appl. Phys. A **110**(4), 929–934 (2013). [CrossRef]

**38. **N. Kinsey, C. DeVault, J. Kim, M. Ferrera, V. Shalaev, and A. Boltesseva, “Epsilon-near-zero Al-doped ZnO for ultrafast switching at telecom wavelengths,” Optica **2**(7), 616 (2015). [CrossRef]

**39. **E. Burstein, “Anomalous optical absorption limit in InSb,” Phys. Rev. **93**(3), 632–633 (1954). [CrossRef]

**40. **T. S. Moss, “The interpretation of properties of indium antimonide,” Proc. Phys. Soc., London, Sect. B **67**(10), 775–782 (1954). [CrossRef]

**41. **B. Campbell, P. Kelly, and L. Kuznetsova, “The interplay between band gap renormalization and Burstein-Moss effect for multilayered Al:ZnO/ZnO metamaterial,” presented at SPIE Optics and Photonics (August, 2020).

**42. **C. S. Granerød, S. R. Bilden, T. Aarholt, Y.-F. Yao, C. C. Yang, D. C. Look, L. Vines, K. Magnus Johansen, and Ø. Prytz, “Direct observation of conduction band plasmons and the related Burstein-Moss shift in highly doped semiconductors: a STEM-EELS study of Ga-doped ZnO,” Phys. Rev. B **98**(11), 115301 (2018). [CrossRef]

**43. **K. G. Saw, N. M. Aznan, F. K. Yam, S. S. Ng, and S. Y. Pung, “New insights on the Burstein-Moss shift and band gap narrowing in indium-doped zinc oxide thin films,” PLoS One **10**(10), e0141180 (2015). [CrossRef]