Dynamic Nanophotonics in Epsilon‐Near‐Zero Conductive Oxide Films and Metasurfaces: A Quantitative, Nonlinear, Computational Model

The promise of dynamic nanophotonic technologies relies on the confinement and spatiotemporal control of light at the nanoscale. Confinement via plasmonics, dielectric resonators, and waveguides can be complemented with materials whose optical properties can be controlled using nonlinear effects. Transparent conducting oxides (TCOs) exhibit strong optical nonlinearities in their near‐zero permittivity spectral region, on the femtosecond timescale. Harnessing full spatiotemporal control over the nonlinear response requires a deeper understanding of the process. To achieve this, a self‐consistent multiphysics time‐domain model for the nonlinear optical response of TCOs is developed and implemented into a 3D finite‐difference time‐domain code. Simulations are compared and tuned against recently published experimental results for intense laser irradiation of thin indium tin oxide (ITO) films, achieving good quantitative agreement; the time‐domain dynamics of the nonlinear response and the phenomenon of time‐refraction are also explored. Finally, by simulating intense laser irradiation of a plasmonic particle on an ITO film, the applicability of the approach to time‐varying metasurfaces is demonstrated. As expected, significant enhancement of the nonlinear response of an ITO‐based metasurface over bare ITO thin films is found. This work thus enables quantitative nanophotonics design with conductive oxides in their epsilon‐near‐zero region.

The promise of dynamic nanophotonic technologies relies on the confinement and spatiotemporal control of light at the nanoscale. Confinement via plasmonics, dielectric resonators, and waveguides can be complemented with materials whose optical properties can be controlled using nonlinear effects. Transparent conducting oxides (TCOs) exhibit strong optical nonlinearities in their near-zero permittivity spectral region, on the femtosecond timescale. Harnessing full spatiotemporal control over the nonlinear response requires a deeper understanding of the process. To achieve this, a self-consistent multiphysics timedomain model for the nonlinear optical response of TCOs is developed and implemented into a 3D finite-difference time-domain code. Simulations are compared and tuned against recently published experimental results for intense laser irradiation of thin indium tin oxide (ITO) films, achieving good quantitative agreement; the time-domain dynamics of the nonlinear response and the phenomenon of time-refraction are also explored. Finally, by simulating intense laser irradiation of a plasmonic particle on an ITO film, the applicability of the approach to time-varying metasurfaces is demonstrated. As expected, significant enhancement of the nonlinear response of an ITO-based metasurface over bare ITO thin films is found. This work thus enables quantitative nanophotonics design with conductive oxides in their epsilon-near-zero region.
The nonlinear optical response of TCOs has been extensively studied in their ENZ spectral ranges, [20,23] where a large nonlinear response has been experimentally observed. This is attributed to intraband transitions occurring in the nonparabolic region of the conduction band of these materials. [32,33] For the case of ITO, different models accounting for the temperature dependence of the plasma frequency, [20,21,34,35] the effect of free-carrier scattering, [36] other carrier kinetics effects, [37] and nonlocality [38,39] have been proposed. However, the temporal nature of the nonlinearity is either not accounted for or, when it has been considered, the modeling of the transient response of electron dynamics is separate from the electrodynamics modeling of optical propagation. [20,21] Currently, no model can account for the temporal dynamics of the electron distribution and the optical properties simultaneously and self-consistently as the incident pulse propagates through and interacts with the material. This poses a limitation to understanding the nonlinear response of complex time-varying nanophotonic structures that include ENZ materials and limits their design.
In this work, a self-consistent multiphysics numerical model is presented that accounts for the temporal modification of the optical properties and electron distribution in TCO-based ENZ materials, as well as fully accounting for electrodynamic propagation. We use a two-temperature model (TTM) to describe the electron excitation and thermalization and account for its effect on the material's optical properties with a temperature-dependent plasma frequency. The implementation of this numerical model in a finite-difference time-domain (FDTD) electrodynamics solver provides access to local optical properties in the time domain. This enables the simulation of the material's nonlinear optical response in space and time and the design and optimization of complex nonlinear time-varying photonic structures incorporating TCO-based ENZ materials. We validate this model against several experimental measurements, make predictions, and gain insight into several phenomena including the temporal dynamics of the intensity-dependent refractive index, time refraction in transmission and reflection, and the nonlinear response of plasmonic metasurfaces interacting with TCO-based ENZ media.
The structure of the article is as follows. We introduce our spatiotemporal, multiphysics material model for the optical nonlinearity of TCOs in Section 2 and provide an intuitive explanation of the physical processes behind the time-varying optical properties. In Section 3, we refine and validate our model by reproducing experimental results from experiments published in ref. [20], exploring the effect of the pump pulse's central wavelength, incident angle, and intensity on the nonlinear optical response of ITO thin films. Comparisons between the simulation results and experimental data demonstrate quantitative agreement. In Section 4, we investigate the time-domain nature of the nonlinearity and compare our findings to pump-probe experiments published in ref. [20]. There, we also present movies of our simulations that show how the time-dependent, nonlinear response of the ENZ material modifies the propagation of the electric field amplitude. In Section 5, we use our computational approach to demonstrate time refraction, the change in the central frequency of a pulse while propagating through a medium with time-varying permittivity. [40] Our results predict a larger frequency shift for reflected light than for transmitted light; only the latter has been measured in experiment. [2] In Section 6, we simulate an ITO-based time-varying plasmonic metasurface, demonstrating the versatility and power of the FDTD implementation of our model. We find that the presence of the plasmonic particle enhances the effective refractive index change by a factor of 5 over a bare ITO thin film. Finally, in Section 7, we conclude our work.

Model for TCO Time-Varying Nonlinear Optical Response
TCOs are wide-bandgap, degenerate semiconductors that possess a metallic optical response in the NIR due to an abundance of free electrons in the conduction band. Therefore, the frequency-domain permittivity of TCOs is well described in the NIR spectral range by the Drude model. [35] εðωÞ where ω p ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the plasma frequency, γ is the damping coefficient, ε ∞ is the high-frequency permittivity, ω is the angular frequency of the electromagnetic field, m* is the effective mass of the conduction electrons, e is the charge of an electron, N is the free electron density, and ε 0 is the permittivity of free space. According to Equation (1), the wavelength at which the real part of the permittivity equals zero (the zero-crossing wavelength) is given by λ ZC ¼ 2π ffiffiffiffiffi ffi ε ∞ p c=ω p , and this identifies the ENZ region.
The properties of the conduction electrons are dictated by the conduction-electron energy band. Typically, the parabolic band approximation can be applied, such that m* does not change upon interaction with light. This is generally not the case in TCOs exposed to intense, ultrafast laser pulses. Due to their high doping concentration (and thus high carrier density), a correction term is required in the dispersion relation [41] where E is the electron energy, ℏk is the magnitude of the crystal momentum, ℏ is the reduced Planck's constant, m Ã 0 is the effective mass at the conduction band minimum, and C is a nonparabolicity factor. This nonparabolic band structure results in a plasma frequency that depends on conduction band electron temperature, T e . An expression for ω p ðT e Þ is derived in Section S1, Supporting Information, where we use the semiclassical Boltzmann equation in the relaxation-time approximation, assuming an isotropic conduction band. The resulting formula is where f FD ¼ ½expððE À μÞ=k B T e Þ þ 1 À1 is the Fermi-Dirac distribution, μ is the chemical potential (discussed later), and k B is the Boltzmann constant. This temperature dependence of the plasma frequency is a consequence of the energy dependence of the electronic effective mass in the conduction band. [35] An alternate derivation of Equation (3) is given in ref. [33], resulting in the same equation.
In Figure 1, we plot both the plasma frequency (evaluated via numerical integration of Equation (3)) and the resulting real component of the Drude permittivity (Equation (1)) as a function of electron temperature (T e ) for several wavelengths, using parameters given in Section 3. The plasma frequency decreases monotonically with T e , leading to a dramatic increase in the relative permittivity. For wavelengths between 1200 and 1300 nm, room-temperature (300 K) ITO is an ENZ material. In fact, for λ > 1240 nm, room-temperature ITO behaves as a metal since Re(ε) < 0, but as the ITO heats, Re(ε) becomes positive, and it transitions to a dielectric. We discuss the dynamics of this metal-to-dielectric transition further in Section 4. Before introducing our model for the electron thermalization dynamics, we first qualitatively describe the effect on the conduction band electron distribution upon irradiation by an intense pulse with a center angular frequency of ω pump . As with other degenerate semiconductors, the Fermi energy E F is dependent on the doping concentration, though it typically lies in the conduction band. Defining the zero-energy level to be at the conduction band edge, the Fermi level of ITO is on the order of 1 eV inside the conduction band. As illustrated in Figure 2, when a pump pulse is incident on TCOs initially in equilibrium at an ambient temperature T e,0 , conduction band electrons with energies E between E F À ℏω pump and E F are excited above the Fermi level, leading to an athermal (nonequilibrium) distribution of electrons (Figure 2b) that is perturbed from the original Fermi-Dirac distribution (Figure 2a). [42] Subsequent electronelectron scattering results in the athermal energy being rapidly redistributed in the electron gas, resulting in a heated Fermi distribution with temperature T e,1 > T e,0 (Figure 2c). Now we turn to deriving our quantitative model. The thermalized electron and lattice subsystems are treated using an extended TTM [15,43] given by C e ðT e Þ ∂T e ðt, rÞ ∂t ¼ Àg ep ðT e , T l ÞðT e ðt, rÞ À T l ðt, rÞÞ þ Uðt, rÞ τ ee ðT e Þ (4b) where T l is the lattice temperature, C e is the electron heat capacity, C l is the lattice heat capacity, g ep is the electron-phonon coupling coefficient, τ ee is the electron-electron scattering time, τ ep is the electron-phonon scattering time, and P is the electromagnetic power absorbed by the conduction electrons. In Equation (4), the electron and lattice subsystems are treated separately, and energy is exchanged between the electron and lattice subsystems via electron-phonon coupling g ep . This version is referred to as the extended TTM since it accounts for the athermal electron energy density U(r, t) (J m À3 ), a scalar field which acts as a thermal energy source to the conduction electrons.
Due to the relatively low Fermi energy of ITO, the electronelectron scattering time τ ee is on the order of a few femtoseconds (becoming subfemtosecond for a hot electron gas) as we describe in Section S2, Supporting Information. This is advantageous in active nanophotonics as it means that the material response is ultrafast and limited mainly by the exciting pulse width. Because time scales of interest are typically on the order of hundreds of femtoseconds (much larger than τ ee ), the effect of including U in our modeling is negligible. For simplicity we thus use a simplified version of the TTM.
The parameters C e , C l , and g ep are all temperature dependent; expressions for calculating them are given in Section S3, Supporting Information, along with a method for calculating the temperature-dependent chemical potential μ.
where E is the electric field, H is the magnetic field, and J is the conduction electron current density that accounts for intraband transitions. For constant ω p , Equation (8) is equivalent to the standard Drude model (Equation (1)). As there are no interband transitions in the NIR for ITO, we use constant permittivity to describe the effects of the bound electrons that are absorbed into ε ∞ , though this can be easily generalized. The TTM and electrodynamics equations are coupled via two terms: 1) the absorbed power P ¼ J ⋅ E, given by Poynting's theorem, and used as the source in Equation (5a), and 2) the plasma frequency ω p ðT e Þ, calculated via Equation (3) and (5), and used as the source of the nonlinearity in the electrodynamics equations. Solved together, Equation (3), (5), (9)-(8) describe the spatiotemporal, thermally driven optical nonlinearity in TCOs.
We focus the implementation of our numerical model to account for the third-order nonlinearity responsible for the temperature-dependent permittivity. A more accurate treatment, especially for geometries with length scales on the order of 10 nm or less, would involve including hydrodynamic terms in Equation (8) and a nonlinear model for the bound electrons as in refs. [38,39]. These extra terms would broaden the functionality of our model by taking into account even-ordered nonlinearities such as second-harmonic generation. After some numerical experiments using our existing hydrodynamics solver, [44] we found that including hydrodynamic terms has very little effect on the permittivity change when solved in conjunction with the TTM. In fact the hydrodynamic terms are nonlocal, adding an unnecessary computational load, [45] and thus we chose to ignore these terms for the present work.
We use an in-house parallel 3D FDTD code [45,46] to solve Equation (9)-(8). In the same code, we also solve the TTM (Equation (5)) using the Runge-Kutta method. [47] At each step in the FDTD simulation, the electric field, magnetic field, current density, and the electron and lattice temperatures are updated self-consistently.

Model Validation and Refinement in ITO Thin Films
To validate the model presented in the previous section, we compare it against experimental results obtained by femtosecond laser irradiation of ITO thin films. Our choice of ITO is due to the wide availability of experimental studies on its nonlinear optical response under various scenarios and the string of recent interest on this material as an ENZ medium. Indeed, our model can be applied to any highly doped semiconductor with an ultrafast thermalization time, as long as the model parameters are known.
Despite the many experimental studies in recent years, [2,20,21,32,35,37,38] due to the sensitivity of ITO's optical properties to manufacturing specifications, we focus first on reproducing the results of the article by Alam, De Leon, and Boyd, [20] in which the intensity-dependent complex refractive index of an ITO thin film in the ENZ spectral region is investigated using several tests. First, the authors measure the effect of the pump central wavelength and incident angle on the nonlinear refractive index of ITO. Second, they show how the reflection, transmission, and absorption shift as a function of the pump intensity. Third, they measure the time-domain dynamics of the nonlinearity via pump-probe measurements of the transmission. We use our model to reproduce the first two of these tests in this section and tackle the third (time-domain dynamics) in the next section. Through these comparisons, we are able to further refine our model to most accurately match this particular set of experimental data. We also perform an optimization procedure The pump pulse acts by exciting electrons below the Fermi level in the range from (E F À ℏω pump ) to E F (red shading) to states above the Fermi level, from E F to (E F þ ℏω pump ) (blue shading). This can be described as a perturbation added to the Fermi-Dirac distribution δf ðEÞ. c) Through electron-electron collisions, the distribution thermalizes to a hot Fermi-Dirac distribution with temperature T e;1 > T e;0 .
www.advancedsciencenews.com www.adpr-journal.com to determine optimal pulse parameters for maximizing the nonlinearity.
We consider the structure shown in Figure 3a. It consists of a 310 nm-thick ITO layer on a glass substrate, bounded by air on top, as in Ref. [20]. The structure is illuminated from the air side with p-polarized light incident at angle θ inc . To model the 1D thin-film geometry (along y), we apply Bloch periodic boundary conditions to the x and z directions of our 3D simulation domain enabling us to use oblique-incident pulses. [48,49] We use the room-temperature Drude model parameters for ITO reported in Ref. [20], given by ω p ¼ 2π Â 473 THz, γ ¼ 0.0468 ω p , and ε ∞ ¼ 3.8055, resulting in a zero-crossing wavelength of Ref. [20] quantifies the nonlinear response by assuming it is Kerr-like with an intensity-dependent refractive index n ¼ n 0 þ n 2 I. The real part of the effective nonlinear index, Re(n 2 ), and the effective nonlinear attenuation constant, β ¼ ð4π=λÞ Im(n 2 ), are extracted via z-scan measurements. It was observed that the maximum Re(n 2 ) and minimum β are found when the pump wavelength is centered at 1240 nm and the angle of incidence is 60°with respect to the ITO surface normal.
Treating the thermal nonlinearity described in Section 2 as Kerr-like is a rough approximation for two primary reasons: 1) there is no reason to assume that higher odd-order terms do not contribute to the nonlinearity and 2) the nonlinearity is not instantaneous. Although our numerical implementation allows us to calculate the nonlinear polarization (or current density) field using more formal definitions, to compare with the experimental results, we calculate an effective n 2 from the simulation data. Though there are several ways to go about this, we use the following computationally efficient method. We calculate the mean conduction electron temperature (T e ) within the ITO over space and time. This is used to obtain an effective plasma frequency and thus, by applying Equation (1), an effective, temperature-dependent, refractive index nðT e Þ ¼ Re( ffiffi ffi ε p ), which differs from the room-temperature refractive index n 0 by Δn ¼ nðT e Þ À n 0 . We then extract an effective n 2 using where I is the time-averaged pulse intensity. The true refractive index n r, t ð Þ will vary in space and time during pulse irradiation and therefore Δn r, t ð Þ ¼ n r, t ð ÞÀ n 0 also varies in space and time. While the time and space dependence of n 2 is not measured in Ref. [20], we do explicitly calculate Δn r, t ð Þ in Section 6 and Section S4, Supporting Information.
We carry out a 2D parameter sweep of pump incident angle (in 10°increments) and center wavelength (in 50 nm increments) and plot the effective n 2 and β values in Figure 3b,c, respectively. Our results indicate that n 2 is maximized at a wavelength of 1250 nm and an angle of 60°, in agreement with the experimental results. Further, our values of n 2 ¼ 0.16 cm 2 /GW and β ¼ À13 000 cm/GW are in reasonable quantitative agreement with the values obtained in Ref. [20], which are n 2 ¼ 0.11 cm 2 GW À1 and β ¼ À7500 cm GW À1 .  www.advancedsciencenews.com www.adpr-journal.com Next, we determine the effect of the incident pulse intensity on the reflectance, transmittance, and absorptance for a central pump wavelength of 1240 nm and an incident angle of 30°as in Ref. [20]. We show the results of our calculations in Figure 3d where we see an increase in the transmittance from T ¼ 0.13 to T ¼ 0.4 as the intensity increases from 0 to 250 GW cm À2 and a corresponding decrease in the reflectance from R ¼ 0.19 to R ¼ 0.02. Our results agree quantitatively with the experimental measurements, [20] which show an increase in the transmittance from T ¼ 0.12 to T ¼ 0.4 and a corresponding decrease in reflectance from R ¼ 0.2 to R ¼ 0.04 for the same intensity range.
A feature of note in Figure 3d is the saturation of the nonlinearity for high intensities. For simulations in which we assume that the Drude damping coefficient γ is temperature independent, this saturation occurs for intensities larger than 250 GW cm À2 . In order to produce Figure 3d, and reach agreement with experiments, we needed to consider a model wherein the damping coefficient increases with temperature. As a firstorder approximation, we use where γ 0 is the low-intensity Drude damping coefficient and T 0 ¼ 15 000 K is a fitting parameter with units of temperature. The increase in the Drude damping coefficient is reasonable because we would indeed expect a decrease in electron mobility due to an increase in electron-phonon, electron-impurity, and grain boundary collisions. This effect was also found experimentally in Ref. [35], where an almost-linear relationship was found between an effective γ and the pulse intensity (and thus, indirectly temperature), by fitting the measured spectra to the Drude model. The fact that γ depends on the temperature highlights the benefit of time-domain modeling as we can extract simple, yet effective models for optical parameters that otherwise would require semiclassical or quantum mechanical calculations. Finally, we optimize properties of the pump pulse (center wavelength, incident angle, and linear chirp for a given intensity) in order to maximize the change in the refractive index Δn(r, t). Using a Nelder Mead algorithm, we find that for lower intensities, Δn is largest when the center wavelength is tuned to the ENZ wavelength. Interestingly, for higher intensities, a larger center wavelength is required to maximize Δn. For all intensities, the optimal incident angle is θ inc ¼ 30°and the linear chirp has negligible effect. These results are presented and further explained in Section S4, Supporting Information.

Analysis of Time-Domain Dynamics
In the previous section, we used spectral results to validate our model. However, part of the interest in ITO's strong nonlinearity stems from its ultrafast response time, limited mainly by the pulse width (since the athermal relaxation time is small), allowing precise control of the instantaneous refractive index. Thus, we are particularly interested in modeling the time-varying optical response of ITO.
Degenerate pump-probe measurements in ITO thin films have shown that the rise time is less than 200 fs for a 150 fs laser pulse and the fall time was %360 fs. [20] Here we present corresponding pump-probe simulations using a two-step numerical process. A primary FDTD-TTM simulation computes the temperature fields as a function of space and time as the material is excited by a pump pulse. These temperature fields are stored and subsequently used in separate probe simulations, where it is assumed that the probe intensity is low enough that it has no effect on the time-and space-varying temperature field.
We plot the transmission at λ ¼ 1240 nm as a function of probe pulse delay in Figure 4, which shows a rise time of 170 fs and a fall time of 385 fs. As in Ref. [20] these are defined as the time it takes for the transmittance to rise (fall) from 10% (90%) to 90% (10%) of its maximum change. In order to obtain Figure 4, we needed to adapt a model for the electron-phonon coupling strength g ep that was originally derived for metals (Equation S34, Supporting Information) and contains a free parameter. Since the fall time of the optical response is controlled by g ep , we compared numerical calculations to the experimental   results to obtain the free parameter; this is presented in Section S3, Supporting Information.
In Figure 5a, we show the corresponding effective refractive index of the ITO film as a function of time, which is extracted using an effective medium theory where we consider the thin film to be an effective (uniform) medium. [50] We see that the real part of the refractive index increases by %0.5 and the imaginary part decreases by %0.2. Such a dramatic refractive index change suggests a transition from metallic (Re[ε] < 0) to dielectric (Re[ε] < 0) behavior in the epsilon-near-zero region, which also explains the large increase in transmittance in Figure 4.
We illustrate this further in Figure 5b,c, which shows snapshots of the electric field amplitude during the simulation of a 150 fs pump-pulse incident on an ITO film on glass. We see in Figure 5b that, early in the simulation, a large amount of the incident light is reflected, as evidenced by the interference patterns, indicating that ITO behaves like a metal. Figure 5c shows the electric field amplitude later in the same pulse where the ITO layer is hotter and there is much less reflection (no discernible interference pattern), indicating that ITO behaves as a dielectric.
Our results show that the coupled FDTD-TTM model for ITO agrees remarkably well with experiments, indicating the strength of the FDTD implementation. Because our modeling is in the time domain, we can track the nonlinear dynamics as a function of both space and time. In the Supporting Information, we present 10 sets of movies showing the evolution of the pump electric field amplitude for different incident angles, center wavelengths, and intensities. Also included are the corresponding videos of the electronic temperature field evolution. These are briefly discussed in Section S5, Supporting Information.
Here we highlight one example set of movies. Movie 1 shows the electric field amplitude evolution of a 250 GW cm À2 intensity pump, centered at 1240 nm, incident at 30°, and Movie 1-T shows the corresponding electronic temperature field evolution. In Movie 1, we see that at low temperatures (before the pump peak), the ITO film exhibits metallic properties as it reflects much of the incident pulse and absorbs much of the transmitted light, resulting in an evanescent wave in the ITO layer. Indeed, the angle of transmission in the ITO is small due to the low refractive index. As the pulse grows in intensity (and T e increases), we see the skin depth of the ITO increase (Im(Δn) < 0), the angle of transmission in the ITO increase (Re(Δn) > 0), and a corresponding increase in transmission into the glass substrate. At this point ITO behaves like a dielectric. Finally, near the end of the pulse, we see the interference pattern in vacuum diminish significantly due to decreased reflectance.

Time Refraction in ITO Thin Films
In the previous sections, we refined and validated our model using experimental measurements published in Ref. [20]. We demonstrated how intense optical pulses can actively change the refractive index of ITO. However, the reverse is also true. A recent publication cleverly investigated the elusive timedomain response of ITO by harnessing a little-known optical effect known as adiabatic frequency translation or more colloquially, time refraction. [2] Adiabatic frequency translation occurs as the result of a wave traveling through a medium with time-varying permittivity. [3,40] The center frequency of an optical pulse leaving the time-varying medium is shifted from that of the incident pulse. Its colloquial title is inspired by space refraction where the wavevector changes due to spatial variations of the permittivity. Here, the frequency (i.e., photon energy) changes due to temporal variations in the permittivity.
It was recently demonstrated experimentally that ITO produces measurable frequency translation when pumped with a high-intensity pulse. [2] Degenerate pump-probe experiments, www.advancedsciencenews.com www.adpr-journal.com with a pulse width of 120 fs, were used to demonstrate time refraction in a 620 nm film of ITO bounded on either side by glass. By changing the pump-probe delay, they demonstrated a time-dependent shift in the probe's center frequency. Using our self-consistent modeling approach, we observe the experimental adiabatic frequency shifting in transmission of Ref. [2]. Moreover, we predict a stronger effect in reflection, where we find that, for certain pump-probe delays, two distinct reflected pulses are created. We show in Figure 6b,c the probe pulse transmission and reflection spectra, respectively, as a function of pump-probe delay for three different pump intensities, where we find that the reflection spectra show a much more pronounced frequency translation. For all pump simulations, we used the same excitation conditions reported in Ref. [2], that is, the pump field is p-polarized with an incidence angle of 15°. Normal incidence was used for the broadband, probe simulations.
We found that the probe pulse must be down chirped to show pronounced frequency translation in our simulations. An unchirped 120 fs pulse was found to contain too few frequency components to show appreciable translation. Requiring chirp is not explicitly mentioned in the main article of Ref. [2]. However, in the accompanying Supplementary Material, the frequencyresolved optical gating (FROG) measurement of their pulse is plotted and is shown to be greatly down chirped. In our numerical experiments, we find that a chirp parameter (as defined in Ref. [51]) of a ¼ À1.5 (where the negative indicates a down chirp) best replicates the experimental data. We believe this insight presents the chirp as a new degree of freedom to explore in future works. (c) Figure 6. Frequency translation results obtained through pump-probe simulations as a function of frequency shift (x-axis) and pump-probe delay (y-axis) for three pump intensities as indicated. a) The material structure consists of a 620 nm film of ITO bounded on both sides by glass. The pump field is p-polarized with an incidence angle of θ inc ¼ 15°. b) Transmission and c) reflection spectra for varying pump-probe delay.
www.advancedsciencenews.com www.adpr-journal.com As the reflected light was not measured in the experiments, our calculation that shows enhanced frequency translation in reflection over transmitted light represents a novel prediction. We see that the reflection spectrum is composed of two distinct peaks for some values of pump-probe delay. To explore this further, we plot in Figure 7 the time domain reflected (blue) and transmitted (red) electric fields for a pump-probe delay of 25 fs and an incident pump intensity of 268 GW cm À2 (right panel of Figure 6c). Such pump intensity has been used in experiments without material breakdown. [20,23] We observe that the reflected field is composed of two distinct envelopes in time, unlike the transmitted field.
Taking a windowed Fourier transform of the reflected timedomain field centered at different times, we obtain an effective spectrum as a function of time, and we plot its evolution as a contour plot in Figure 7b. We see that the spectral content of each distinct time-domain lobe is different. In Figure 7c, we plot the spectra at the times corresponding to the peaks of the two lobes (indicated by the vertical dotted lines in Figure 7b); we see that the two spectra are indeed shifted.
The shift itself is not surprising, since we used a chirped pulse in the simulation as described earlier. Indeed, when we repeated this simulation for the low-intensity limit where a constant temperature and thus constant plasma frequency is maintained (not shown), we also see a shift but in a single lobe. As the peak instantaneous wavelength shifts from 1300 to 1200 nm, the real part of the relative permittivity at the peak instantaneous wavelength increases from %À0.2 to þ0.4 as it crosses through the ENZ region; this is seen in Figure 1 for the case of low temperature.
The two-lobe reflected pulse only occurs for high intensity, as a direct consequence of the nonlinearity. As the temperature increases, the plasma frequency decreases, causing a corresponding increase of the relative permittivity up to a value of 2 or more. During the intense pulse irradiation, the effective index of ITO briefly comes to match that of the surrounding medium (SiO 2 ), leading to a null in reflection that does not occur in transmission.
In sum, the ability of this model to quantitatively replicate experimental results in the time and frequency domain gives us confidence in its accuracy and the novel phenomena it predicts. This allows us to explore even more complex scenarios, to which we now turn.

Time-Varying Plasmonic Metasurface on ITO
In the previous sections, we refined and validated our model using experimental measurements published in Refs. [2,20] using both spectral and time-domain results. It was found that the pump's wavelength and angle of incidence have a strong effect on the magnitude of the refractive index change. Intuitively, one might suppose that the strength of the nonlinear optical response can be further increased by enhancing the near fields within ITO films by including nanostructures. Plasmonic nanoparticles have been used to enhance optical processes of all sorts, including, nonlinear optical effects. [52] As our model is built into FDTD, we are not restricted to simulations of the simple planar geometries presented above, but we can simulate the refractive index shift in ITO in the presence of nanostructuring, including nearby plasmonic particles, where the 3D spatial structure and dispersion of the plasmonic material are taken into account self consistently.
As an illustrative example, we consider the plasmonic metasurface depicted in Figure 8a. It consists of a squared array of gold rectangular nanoparticles on a 310 nm-thick layer of ITO on glass. The array is characterized by a 600 nm Â 600 nm unit cell. The rectangular nanoparticle has a length of 500 nm, width of 300 nm, and thickness of 40 nm. This geometry was chosen because it has a very large absorption peak of A ¼ 0.9 within the www.advancedsciencenews.com www.adpr-journal.com wavelength range of interest; this peak occurs at λ ¼ 1350 nm. We use Bloch boundary conditions in our FDTD simulations, where a plane wave that is linearly polarized along the length of the particle is incident on an infinitely periodic metasurface at an angle of 30°. As in Section S4, Supporting Information, we use a 2 GW cm À2 pump pulse for ease of comparison to our results there. We save 3D temperature field snapshots within the ITO at various times during the pump pulse simulation-the one at the pump pulse peak is shown in the inset of Figure 8a. For each temperature field snapshot, we perform separate probe simulations. From the resulting probe reflection and transmission coefficient spectra, we use effective medium theory [50] to extract the effective refractive index of the metasurface as a function of wavelength and pump pulse simulation time. We plot in Figure 8b the effective index change of the metasurface as a function of pump simulation time at λ ¼ 1350 nm. We find that the maximum value of Δn is five times larger than for the bare ITO thin film (Section S4, Figure S4b, Supporting Information) irradiated at the same incident pump intensity. In Figure 8c we plot the corresponding change in reflectance and transmittance, ΔR and ΔT, respectively. We find that the transmittance decreases during the pulse irradiation (ΔT < 0), in contrast to the increase that was demonstrated for bare ITO thin films (ΔT > 0) previously (Figure 4 and S4c, Supporting Information). In addition to this change in sign, the magnitude of ΔT is also two times larger than in Figure S4c, Supporting Information. The large enhancement in ΔT here is due to the plasmonic particle that creates electromagnetic hotspots when excited near the localized surface plasmon resonance. As their names imply, these hotspots are regions of enhanced Joule heating, leading to a stronger nonlinearity being induced in the ITO film (see Figure 8a). By comparing Figure 8c and S4c, Supporting Information, we also observe a longer transient nonlinear response when the plasmonic particle is present. This occurs because the decay time of the plasmon resonance is longer than the decay time of the ITO nonlinear response.
Effective medium theory does not only output an effective refractive index n ¼ ffiffiffiffiffi εμ p but also an effective admittance Typically, for sufficiently small feature sizes, we expect these values to be equal as long as the material is nonmagnetic (μ ¼ 1). We have verified that this indeed remains true for bare ITO films considered in the previous sections. However, this is not the case for the metasurface investigated in this section, where applying effective medium theory suggests an effective magnetic response of the metasurface. In fact, we can see this effect in Figure 8 where we see an increased imaginary component of the effective index (Figure 8b) coupled with a decrease in the absorptance (Figure 8c); this can only occur if μ 6 ¼ 1. This effective magnetic behavior is due to complex response of the plasmonic particle, and will fundamentally change the ENZ behavior, because we must now take μ into account. This will be investigated in future work.
Previous experimental work has investigated ultrafast irradiation of ITO-based time-varying plasmonic metasurfaces. [21] There, it was also found that the transmittance decreases during  Figure 3a. The structure is illuminated with a p-polarized optical field incident at an angle of θ inc ¼ 30°. The inset shows a 3D time-domain snapshot of the temperature field in the ITO film at the peak of the laser pulse, within a unit cell of the plasmonic particle array. b,c) The change in refractive index and the change in reflectance and transmittance, obtained using the FDTD-TMM model, respectively.
www.advancedsciencenews.com www.adpr-journal.com pulse irradiation for the metasurface but increases (with a smaller line-width) for a bare ITO film. This further validates the results of this (and the previous) section. The magnitudes of ΔT and Δn here are not directly comparable to those of Ref. [21], as their ITO sample was different (with different film thicknesses and optical properties), and they applied a different pump intensity. However, the enhancement of Δn when a plasmonic particle is present versus a bare ITO film is comparable. We also present movies showing the time evolution of the electric field amplitude and temperature field of the plasmonic metasurface on ITO while being irradiated by an oblique incidence pump pulse. In Movie 11, we see the incident pulse exciting the plasmonic particle and creating electromagnetic hotspots at the particle edges. These result in electronic temperature hotspots, as shown in Movie 11-T and in the inset of Figure 8a. As the pump pulse peaks, we start to see strong interference between the incident and reflected pulse in the vacuum region due to the increased reflectance (as also indicated by the ΔR plot in Figure 8c). There is very little decrease in the transmittance to balance the reflectance (ΔT plot in Figure 8c). That, in Movie 11, is evident by the increased penetration of the field into the ITO film, suggesting an increased skin depth during intense field irradiation.

Conclusion
We introduced a self-consistent, multiphysics, spatiotemporal numerical model for the intensity-dependent, time-varying refractive index of TCO-based ENZ materials. The nonlinear optical response of TCOs is described by Drude-like permittivity where the plasma frequency is dependent upon the conduction band electron temperature. A TTM describing the electron temperature distribution in time and space as induced by irradiation with an ultrafast intense laser pulse was implemented directly into a 3D FDTD electrodynamics solver. Our method captures the full nonlinear spatiotemporal optical response of TCO-based ENZ materials, including those within complex geometries containing, for example, plasmonic nanostructures in close proximity. We validated our implementation via comparisons with experimental results for a bare ITO film, where we obtained quantitative agreement, as well as an ITO film decorated with plasmonic nanoantennas. We also used our model to demonstrate time refraction wherein the time-varying permittivity of the ITO film results in a shifting of the central frequency of an optical pulse. We show that the reflected light is predicted to have a larger and more noticeable frequency change than the transmitted light. The ability to model and understand the complex nonlinear response of TCOs will be instrumental for designing sophisticated nonlinear nanophotonic devices with femtosecond-scale response times, such as ultrafast reconfigurable metasurfaces for active switching, modulation or beam manipulation.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.