Geometry Defines Ultrafast Hot‐Carrier Dynamics and Kerr Nonlinearity in Plasmonic Metamaterial Waveguides and Cavities

Hot carrier dynamics in plasmonic nanorod metamaterials and its influence on the metamaterial's optical Kerr nonlinearity is studied. The electron temperature distribution induced by an optical pump in the metallic component of the plasmonic metamaterial leads to geometry‐dependent variations of the optical response and its dynamics as observed in both the transmission and reflection properties of the metamaterial slab. Thus, the ultrafast dynamics of a metamaterial's optical response can be controlled via modal engineering. Both the transient response relaxation time and magnitude of the nonlinearity are shown to depend on the modal‐induced spatial profile of the electron temperature distribution and the hot‐electron diffusion in nanorods. The nonlocal effects, depending on the excitation‐induced losses in the metal, are shown to dictate the modal structure of the metamaterial slab and the associated dynamics of its nonlinear response. The opportunity of controlling the electron temperature profile induced in the plasmonic nanorods by changing the metamaterial's geometry and/or excitation conditions paves the way to achieve controllable dynamics of the nonlinear optical response for free‐space as well as integrated nanophotonic applications involving nonequilibrium electrons.


Introduction
The dynamics of excited, nonequilibrium carriers in plasmonic materials determines many interesting and important phenomena, such as hot-electron-induced photochemical plasmonic metamaterial can be dynamically modified by control light illumination, achieving modulation and switching of signal light. [4,5] Two-photon absorption was identified as the mechanism driving the nonlinear response of a gold film patterned with split-ring resonators. This metamaterial has been reported to achieve a modulation of up to 40% with a fluence of the optical control pulse of 270 µJ cm −2 when both the signal and control wavelengths are around the plasmonic resonance of the metamaterial. [6] Au thin films nanostructured with gratings can provide ultrafast modulation of surface plasmon polariton (SPP) excitation using interband electron transitions. [7] Nonlinear effects in ultrathin long-range plasmonic waveguides have been demonstrated in the spectral range of intraband absorption. [8] Due to its inherent nature related to the hotelectron distribution within a nanosctructured conductor in both geometrical space and within the conduction band, the dynamics of free-electron nonlinearity depends on the excitation pulse duration as well as on the geometry of the nanostructured system. In a simple smooth-film geometry, these dependences can be described via the nonlinear Schrodinger equation, [8] and thus inherently depend on the absorption of the system.
The optical properties of plasmonic metamaterials with hyperbolic dispersion are strongly sensitive to the changes of the refractive index of their constituents, thus providing an excellent platform for studies of Kerr-type nonlinear functionalities and their dynamics. [2,3c,9] Hyperbolic dispersion is the signature of anisotropic metamaterials in the spectral range where the components of the effective permittivity tensor have opposite signs. In this frequency range, the metamaterial supports bulk plasmon polaritons with large wavevectors. [10] If all the components of the permittivity tensor are positive, conventional elliptic isofrequency surfaces are observed when nonlocal spatial dispersion effects can be neglected. [11] The transition from elliptic to hyperbolic dispersion takes place through the so-called epsilon-near-zero (ENZ) regime when the real part of one of the permittivity tensor components vanishes. The ENZ regime occurs in the spectral range near the effective plasma frequency of the metamaterial and depends on both dimensional parameters and material composition of the metamaterial. [10,12] Strong and ultrafast nonlinear response of a gold nanorod-based hyperbolic metamaterial under interband femtosecond excitation has been achieved, recording a change in the metamaterial's transmission of up to 80% with a subpicosecond recovery time. [2] Such an enhanced nonlinearity was identified as originating from nonlocal (spatial dispersion) effects associated with the free-electron dynamics related to cylindrical surface plasmons. [13,14] In this paper, we numerically analyze the hot-electroninduced nonlinear optical dynamics of plasmonic nanorods forming a metamaterial by following its Kerr-type nonlinear ultrafast response under interband and intraband excitations. We demonstrate how the modal properties of these metamaterials, including waveguided and cavities modes, can influence the nonlinear behavior of the medium as well as its temporal response in different dispersion regimes, despite being governed by the same electron scattering processes in gold. In particular, we characterize the interplay between the spatial distribution of the mode excited by the control light field, the associated microscopic electron temperature profile in the nanorods, and the nonlinear optical dynamics of the metamaterial as probed in transmission and reflection measurements. The methodology and results of the present study are also relevant to the design of hybrid photonic devices, for example based on two-dimensional materials, such as, e.g., graphene and MoS 2 , combined with plasmonic films or nanostructures to increase light-matter interaction strength. [15] The spatially resolved, ultrafast response of the plasmonic metamaterial described here can be used to engineer ultrafast, strong optical modulation and optimize light harvesting when coupled with such monoatomic layers. The possibility of spatially shaping the electron temperature profile by controlling the metamaterial's geometry at the nanoscale opens up an additional degree of freedom for designing the strength and dynamics of the nonlinear response of integrated photonic components and nanodevices but also has implications for hot-electron extraction and applications in photochemistry and optoelectronics.

Metamaterial Geometry
We consider a plasmonic metamaterial geometry consisting of an assembly of gold nanorods embedded in a porous alumina (AAO) matrix and supported by a transparent substrate (Figure 1a). The nanorods, vertically aligned in the z-direction in the Cartesian frame of Figure 1a, form a uniaxial metamaterial that can be represented within the effective medium theory (EMT) by a diagonal permittivity tensor of the form ε = diag[ε x , ε y = ε x , ε z ], where the spectral dependence is implicit and the z-direction is along the optical axis of the metamaterial. Expressions for the components of the permittivity tensor can be obtained via the Maxwell-Garnet effective medium approximation as 2 is the nanorod concentration, with r and d being the nanorod radius and the average separation between nanorods respectively, ε Au and ε h are the permittivities of Au and the host medium (anodized aluminium oxide, AAO), respectively. [16] This model is valid, in general, away from the Brillouin zone edges of the nanorod array of period d, but fails to accurately reproduce the optical response of the metamaterial when losses are "small," in which case the nonlocal corrections to the effective medium theories need to be considered. [11,13] The spectral dependence of the effective permittivity components of the studied metamaterial is shown in Figure 1b,c for a metamaterial with a nanorod radius of r = 15 nm and an average interrod distance of d = 90 nm. Here, the permittivity of Au was modeled as gold intra i nter ε ε ε = + , where ε intra is a Drude contribution, describing the response from free-electrons in the 5sp conduction band, while ε inter accounts for electronic processes involving interband transitions from the 5d to the 5sp bands as a Lorentz-type contribution. [17] The intraband contribution to the permittivity was corrected to take into account the restricted mean free path R of electrons in the electrochemically grown nanorods with respect to the mean free path in bulk Au (R bulk = 37.5 nm). [13a] Depending on metamaterial postprocessing conditions, such as thermal annealing, the restricted mean free path can take typical values ranging from R = 3 nm to R = 30 nm.
The real part of the effective permittivity component Re(ε x,y ), for light polarized normal to the nanorod axis, is similar to that of a resonant dielectric while the permittivity tensor component Re(ε z ), for light polarized along the nanorod axis, varies monotonically from positive values at short wavelengths to negative values at longer wavelengths, crossing zero at around 830 nm, defining a transition in the metamaterial's dispersion between the regions of elliptic dispersion where Re(ε x,y,z ) > 0 and hyperbolic dispersion where Re(ε x,y ) > 0 and Re(ε z ) < 0. The cross-over frequency from the elliptic to the hyperbolic regime is defined by the metamaterial's effective plasma frequency z (Re( ) 0) p eff ω ω ε = = [10] around which the metamaterial behaves as the so-called ENZ medium. Nonlocal spatial dispersion effects become especially important in the description of the optical response of the metamaterial in both linear and nonlinear regimes as Im(ε z ) → 0. [2][3][4][5][6][7][8][9][10][11] In particular, when modal losses are low enough, the metamaterial supports the so-called additional (nonlocality-enabled) electromagnetic waves, with a significant contribution to the optical response of the system as is discussed below. [3b,11,13] These additional modes provide highwavevector waves when ε h + ε Au < 0, leading to strong local fields supported by the metamaterial. It should be noted that these spatial-dispersion effects are automatically taken into account in the numerical simulations of the composite nanorod metamaterial which do not rely on effective metamaterial parameters. Figure 1a illustrates a schematic of the pump-probe experiment reproduced in the simulations. The metamaterial is illuminated with an ultrashort control pulse and a white-light continuum probe, both transverse-magnetic (TM) polarized. The probe pulse is used to characterize the modifications of the optical properties of the metamaterial induced by the control pulse either through reflection or transmission spectra at increasing time delays following the excitation by the control pulse.

Electron Gas Dynamics Modeling
The control pulse has a Gaussian envelope with a timedependent power P(t) = P 0 exp(-(t -t 0 ) 2 /t 2 p ). The pulse is centered at t 0 = 100 fs and has half-width t p = 50 fs. Due to absorption in gold, some of the pulse energy is deposited in the electron gas in the nanorods. The excitation process is followed by fast electron gas thermalization, establishing a Fermi-Dirac energy distribution of the electrons, characterized by an elevated temperature T e . [2,18] This increased electron temperature is then mapped to a corresponding change of the Au permittivity through a random phase approximation (RPA) model, allowing modelling of transient optical properties of the metamaterial. [19] Figure 1b,c exemplifies the modifications of the effective permittivity tensor components for two different electron temperatures. In general, the component ε z experiences an increase in both its real and imaginary parts in most of the frequency range except in a narrow range near the interband transitions. The effective plasma frequency p eff ω of the metamaterial decreases as the condition Re (ε z ) = 0 red-shifts with the increasing electron temperature and, at the same time, the metamaterial experiences increased losses, which may contribute to changes in its nonlocal behavior.
The two-temperature model (TTM) is used to describe dynamics of the thermalized electron gas in gold, governed by electron-phonon and phonon-phonon scattering processes. Subsequently to the thermalization phase mainly achieved within the first few hundred femtoseconds through electron-electron scattering, electron-phonon scattering becomes the main relaxation pathway for times up to several picoseconds while phonon-phonon scattering dominates on typical timescales of tens to hundreds of picoseconds. Therefore, the simulation of nonlinear optical processes with the THz response rates, relevant to recent experimental observations and practical applications, can neglect the influence of both electron-electron and phonon-phonon scattering on the metamaterial dynamics to a good approximation. While this affects the overall simulated dynamics of the system, it does not impact the interrelation between field distributions and ultrafast dynamics at the nanoscale, which is of interest to us here, but rather greatly simplifies numerical simulations and interpretations. With these assumptions, the temporal evolution of the electron temperature can be described by a TTM, [20] which is based on the supposition that the excited hot electrons reach thermal equilibrium instantaneously, so that the electron gas can always be described by a Fermi-Dirac distribution with a well-defined equilibrium temperature T e . The time-space evolution of electrons temperature then follows the coupled differential heat equations [20] where C e (T e ) = 67.69 J m −3 k −2 × T e and C p = 3.5 J m −3 k −2 are the volumetric heat capacity of electrons and phonons, respectively, T p is the phonon temperature, k(T e ) = k 0 (T e /T p ) = 318(T e /T p ) W m −1 K −1 is the electron thermal conductivity, G = 2 × 10 16 W m −1 K −1 is the electron-phonon coupling constant, and s r t  ( , ) is the incident light power density absorbed by the metamaterial. During the first few picoseconds following the optical excitation, the electron temperature increase exceeds that of phonons by up to several orders of magnitude due to the much higher heat capacity of phonons, providing good ground to assume T p being constant during the time of interest to us. The second equation of the TTM is thus neglected in the simulations, keeping the phonon temperature constant at 300 K. Equation (1) is solved numerically for the electron equilibrium temperature in the nanorod geometry, so that a temperature-dependent permittivity distribution in the nanorods is monitored at various time delays following the excitation.
The time-dependent modifications of the effective permittivity tensor obtained using the TTM are plotted in Figure 2 as the relative changes of the tensor components, assuming uniform electron temperature in the nanorods. The strongest relative changes of the real parts of effective permittivities are observed at different wavelengths for different light polarizations, with Re(ε z ) varying strongest near the effective plasma frequency, while Re(ε x,y ) and the imaginary parts of the tensor components experience the biggest change in the range of interband transitions. Both real and imaginary parts of the permittivity components are modified, with the imaginary part being affected strongest in relative terms. While effective parameters provide interesting insight and guidance for intuitive nonlinearity understanding, how the modifications of the effective permittivity are reflected in the changes of the optical properties depends on the modal structure of the metamaterial slab, and, generally, considerations beyond the local effective medium theory are required to adequately describe the processes involving electromagnetic fields inside the metamaterial. [11] In the following, we will make use of the effective medium properties solely to understand the results of numerical simulations of the nonlinear behavior of the metamaterial, which is based on direct solutions of vectorial Maxwell equations for a nanorod array.

Nonlinear Response Modeling
Transient spectra recorded in a pump-probe experiment can be reproduced in numerical simulations applying two independent steps. In a first step, the pump-induced changes in the Au nanorod permittivity are simulated at different times following optical excitation. This first step provides a set of time-dependent Au permittivities, which are introduced in the linear optical simulations in a second step to compute the transient transmission or reflection of the metamaterial as independent steady-state effective permittivities seen by the probe light. This approximation is valid as long as nonlinear interactions between pump and probe pulses can be neglected. In our study, we considered the AAO matrix to be nondispersive and with linear optical response in the frequency range and pump intensity range of interest. The effect of modal structure on the ultrafast dynamics of the metamaterial was further assessed by comparing it to the dynamics of the same system when the spatial profile of the electron temperature is neglected. In this geometry, the dynamics was simulated setting the electronic temperature to be homogeneous across the nanorods with a value set by the TTM without the diffusion term.

Linear Optical Properties of the Metamaterial Slab
In order to understand the modal structure of the metamaterial, both its linear reflection and extinction spectra were calculated with local and nonlocal EMTs (Figure 3). [11,13] Here, we restrict the study to the optical properties of the metamaterial under TM-polarized plane-wave excitation with which both elliptic and hyperbolic dispersion regimes can be assessed. The local EMT assumes that the metamaterial supports propagation of a single TM-polarized mode with the dispersion given by the EMT parameters in Figure 1b,c. The nonlocal EMT shows the existence of two TM-polarized modes simultaneously propagating in the metamaterial. Both models adequately describe the optical response of the nanorod metamaterial for relatively small angles of incidence (AOI). However, as the AOI increases, the accuracy of the local EMT diminishes. At the same time, the nonlocal EMT adequately describes the optical properties of the nanorod metamaterial for the whole range of frequencies and incident angles and is in excellent agreement with the full vectorial numerical simulations of the nanorod array composite (Figure 4a,c). The extinction dispersion, plotting angular and spectral dependence of OD = −log 10 (T), where T represents the transmission of the metamaterial, reveals two distinct resonances at around 525 nm (2.36 eV) and a double-peak resonance at 820 nm (1.51 eV) and 730 nm (1.7 eV) ( Figure 3a). The shortwavelength peak is essentially dispersion-less and corresponds to the resonance in the spectrum of ε x seen in Figure 1b,c which can be excited with an electric field component polarized normal to the nanorod axes. For light polarized along the nanorod axes, the extinction has strong angular dependence. In particular, a mode splitting is observed at around 30° near the effective plasma frequency, which corresponds to the existence of an additional TM wave in the metamaterial as compared to the local response. This is a signature of nonlocality originating from the cylindrical surface plasmons nature of the modes supported by the nanorods. [3b,13] This behavior is determined by the internal structure of the metamaterial, which cannot be recovered in a conventional, local effective medium theory.
The minima in the reflection dispersion correspond to the incident light coupling to cavity modes of the metamaterials slab (Figure 3b). Above the light line in the superstrate (air), these modes are unbound Fabry-Perot-type resonances of the planar metametarial slab accessible from both the superstrate and the substrate. Between the light line in air and the light line in the substrate, the modes are evanescent at the metamaterial-air interface and leaky in the substrate, therefore, only accessible to plane waves from the substrate side at the angles of incidence above the critical total internal reflection angle. [10] These modes have the unusual property of low negative group velocity, appearing relatively flat in the dispersion plots. [10,21] They can be identified from their characteristic field distributions with several lower order modes (q = 1, 2, 3, 4) shown in Figure 3b. The TM-guided modes above the effective plasma frequency originate from the additional wave in the nonlocal material and exist in the range of frequencies where ε h + ε Au < 0. The predictions of the nonlocal EMT agree with the full-wave solutions of Maxwell equations that take into account the composite structure of the material.
To underline the importance of nonlocality in the description of the observed guided modes, we have repeated the simulations with a conventional (local) EMT (Figure 3 c,d). Both extinction and reflection spectra fundamentally change as the result of this approximation. The double resonance observed in the nonlocal extinction disappears in the local regime, leaving only a single resonance at a wavelength of around 820 nm. In addition, the negative group velocity modes are only present below the effective plasma frequency where the local effective medium theory provides adequate description of the optical response of the metamaterial. Only TM1 and TM2 modes are predicted by the local EMT model. Higher order modes are eneabled by nonlocality and cannot be predicted by the local EMT. Interestingly, the increased losses in Au primarily affect the high-index nonlocal waves, thus driving the metamaterial response toward the predictions of the local EMT.
We will now study the potential of various modes of the metamaterial slab for optimization of the magnitude and dynamics of the nonlinear response of the metamaterial.

Transient Nonlinear Dynamics for Local and Average Electron Temperatures
The dynamics of the optical response of the metamaterial can be where O symbolizes the monitored optical properties, such as, e.g., reflection or transmission determined by the modal structure of the metamaterial, ε is the effective permittivity of the metamaterial, and T e is the electron temperature (Equation (1)). The first term on the righthand side depends on the modal properties of the metamaterial slab and incorporates the spectral dependence of the resonant response of the supported modes. The second term describes how the effective permittivity depends on an electron temperature, which is modelled via the RPA. The third term bears the intrinsic transient dynamics of the electron gas, which can be modeled via the TTM. Both the electron temperature diffusion term and the source term s r t  ( , ) in Equation (1) are of importance when considering nanostructured materials because they incorporate a spatial dependence of the electron temperature dynamics, being driven by the modal properties of the metamaterial. Indeed, the energy deposited by the excitation pulse in the metamaterial, which is expressed as α ω being the polarizability of mode i. The choice of the illumination conditions will, therefore, determine the energy transferred from the excitation pulse to the metamaterial and, through Equation (1), the subsequent dynamics of hot electrons and the associated nonlinear response. Figure 4 shows the correspondence between the selected metamaterial modes, which can be excited by plane waves, and the induced electron temperature. The volume-averaged temperature increase is of the order of 10 3 K for a 0.1 mJ cm −2 fluence of the excitation light at a 465 nm wavelength and is weakly angular dependent, with an absorption mainly determined by a small spectral overlap with the maxima of the imaginary parts of the effective permittivity (Figure 4a). In contrast, detuning the excitation conditions from both material and geometrical resonances, e.g., for a wavelength of 820 nm at normal incidence, leads to a relatively negligible temperature increases on the order of 10 K (Figure 4a). For this wavelength, a more significant ≈10 3 K increase of electron temperature is only observed at oblique incidence for which the excitation is resonant with a mode supported by the metamaterial slab (Figure 4b). In addition to the magnitude of the absorption, tuning the illumination conditions to excite specific metamaterial modes enables the control of the electron temperature spatial distribution (Figure 4). This is particularly remarkable when the excitation of the metamaterial in reflection (Figure 4c) clearly recovers the spatial field distribution of modes q = 2 through q = 5 in the electron temperature maps. Moreover, while the third term in Equation (1) allows to map the incident power to the position-dependent absorbed power, the heat diffusion term in Equation (1) establishes a direct link between the modal field distribution and the local electron relaxation dynamics. Assuming a temperature-independent electronic temperature diffusion rate k T k ( ) e ≈ , this contribution to the dynamics is governed by the curvature of the temperature profile ∇ 2 T e which together with the electron-phonon coupling term establishes an interrelationship between modal properties Figure 5. The dynamics of a) volume-averaged electron temperature in the nanorods for various excitation conditions. b) Position-dependent electron temperature at different points inside the nanorods for a TMpolarized excitation pulse (820 nm wavelength, 20° AOI, 0.1 mJ cm −2 fluence) at t = 200 fs plotted in the plane containing the symmetry axis of the nanorod. c) Relative instantaneous electron temperature dynamics  showing non-uniform, positiondependent electron cooling in the nanorod. and ultrafast time response. This enables an opportunity to use nanostructured systems and their intricate and tailorable resonant field distributions to design the ultrafast response at the nanoscale by geometrical means.
Let us now compare the dynamics associated with the local electron temperature effects to the volume-averaged dynamics of the electron gas. In the latter case, we assume that a uniform electron temperature distribution along the nanorods is established instantaneously. In this case, the time constant of the recovery of the electron temperature, as expressed by the TTM, is determined by the electron-phonon scattering rate for a given absorbed energy, since the difusion term is absent (Equation (1)). Figure 5a shows that the general trend is for this time constant to increase with the amount of energy deposited in the electron gas by the excitation pulse, as expected from the temperature dependence of the elecron scattering rates. However, Equation (1) tells us that the electron temperature relaxation dynamics is potentially more complex when considering resonant structures in the presence of the heat diffusion term.
To show the influence of the diffusion term in the electron cooling dynamics, Figure 5c (Figure 5c). In particular, while at some locations within the nanorods the electron relaxation is accelerated by the inhomogeneous electron temperature distribution by more than 40%, other locations show a slowed down relaxation process by more than 120%, with an average dynamics described by the TTM in the absence of the diffusion term.
The transient spectral maps ΔOD(λ, t) for the metamaterial with the linear optical properties as in Figure 3a are shown in Figure 6a for the excitation conditions as in Figure 5c,d when probed in transmission at an AOI of 20°. As mentioned above, the dynamics of the optical properties of the metamaterial at selected wavelengths is not necessarily that of the electronic comes into play when characterizing the dynamics of optical properties. In particular, the resonant response of the metamaterial adds to the complexity of the dynamics through this dependence, which is determined by both the absorptive and dispersive behavior of the resonances present in the transient extinction or reflection spectra. This will usually result in a non-monoexonential relaxation behavior of the nonlinear response at any given frequency of the resonant nanostructured material. Selected spectral cross-sections of the nonlinear dynamics (Figure 6a,b) show the absorptive behavior of the high-frequency resonance, which is governed by a transient broadening, and the more complex spectral behavior in the hyperbolic spectral range. The impact of the position-dependent absorption (and, thus, the position-dependent electron temperature) on the transient optical properties of the metamaterial is shown in Figure 6c, which plots the transient dynamics of Figure 6a relative to the dynamics of the metamaterial assuming a uniform electron temperature across the nanorods. The strongest difference in dynamics is observed at very short time scales (within 100 fs of the excitation) in the ENZ regime (Figure 6d). In fact, the relative transient transmission change for the nonuniform electron temperature exceeds by almost one order of magnitude that of the uniform temperature assumption at the wavelength of 820 nm (cf. Figure 7). This difference results from the strong spatial overlap between the excitation-induced electron temperture distribution and the field distribution of the probe as clearly revealed in the transient transmission for both the nonuniform and uniform electron temperature distributions (Figure 7a). It shows that for both metamaterial resonances accessible in transmission, at around 525 nm and 820 nm, an enhanced nonlinear response is recorded in the wavelength range corresponding to the maximum spatial overlap between the probe mode and the spatial distribution of the excitation-induced permittivity changes. Here, this maximum overlap is achieved when probing in the ENZ spectral range where the excitation energy is initially deposited. In addition to affecting the magnitude of the transient extinction in the early times following optical pumping, the dynamic response in the ENZ regime experiences a deacreased recovery time. In fact, the dynamics for the nonuniform temperature distribution shows two time constants of 0.1 ps and 1.4 ps, compared to a single time constant of 1.5 ps obtained for the uniform temperature (Figure 7b). The short sub-ps time constant is related to the modal mapping of the electron temperature distribution, governing the dynamics of the transmission at the time scales of heat diffusion within the nanorod. Both uniform and nonuniform temperture dynamics converge within a time frame of ≈400 fs, corresponding to the time in which the electron temperature has mostly homogenized within the nanorod, minimizing the contribution of the heat diffusion term in Equation (1) with respect to the electron-phonon scattering term (Figure 7c).

Nonlinear Optical Response of the Metamaterial Slab
We now consider the full dispersion dynamics of both the extinction and reflection of the metamaterial after the excitation by a control pulse at a wavelength of 820 nm and an angle of incidence of 20°. This excitation configuration has been chosen as it induces the electron temperature changes which overlap the spatial distribution of various metamaterial modes, thus providing a pathway to assess their nonlinear optical properties. The restriction path was kept at R = 30 nm allowing the investigation of the role of a rich modal structure of the metamaterial slab and the role of nonlocal behavior. The nonlinear response of the metamaterial slab in the local EMT regime has been studied elsewhere. [9] The dispersions of linear extinction ΟD = −log 10 (T 0 ) and reflection R 0 are shown in Figure 8a,c together with the associated disperions of the transient extinction ΔΟD = −log 10 (T/T 0 ) and ΔR = (R − R 0 ) shown in Figure 8b,d at the time 200 fs after the excitation (here, T and R are the time-dependent transmission and reflection of the metamaterial, respectively). In the wavelength range of the elliptic dispersion, the strongest changes in the extinction are observed near the resonance of ε x , essentially determined by an increased damping of the resonance (Figure 8e) as seen in Figure 6a,b. This nonlinear response is not sensitive to the angle of incidence in this spectral range. The induced extinction changes in the hyperbolic dispersion regime are strongly angular dependent but also follow the extinction spectrum features. Similarly, the nonlinear response in reflection are also governed by the modal structure of the metamaterial slab (Figure 8d,f) with the induced changes being much stronger for the leaky modes below the light line than the Fabry-Perot-type modes above the light line.
Varying the excitation wavelength throughout the elliptic and hyperbolic dispersion range and different angles of incidence of the excitation pulse mostly influences the modulation depth and time response, while the shape of the spectral response remains mostly unchanged, completely determined by the linear modal structure of the metamaterial slab. Since the modulation depth is determined by the electron temperature increase induced by the illumination, it is stronger for the wavelength/angle of incidence combinations for which strong absorption is reached. The absolute modulation in reflection reaches, in absolute values, more than 10% for cavity modes and 30% for guided modes for all pump-probe wavelength and angles of incidence configurations for which these resonances are probed. Interestingly, due to the dominating nonlocal mode structure in the elliptic regime, the nonlinear changes are very strong on both sides of the effective plasma frequency of the metamaterial. When losses in Au are increased, e.g., in the electrodeposited, unannealed Au nanorods, the situation changes and nonlinearity in the elliptic regime is strongly suppressed.

Conclusion
We have studied the hot-electron dynamics and diffusion in plasmonic nanorods forming a metamaterial composite. A microscopic model allows us to take into account the nonuniform profile of the electron temperature distribution within the nanorods, which depends on the modal distribution at the excitation wavelength and determines the magnitude and the dynamics of the transient nonlinear response. This opens up the possibility to study nonequlibrium electron dynamics in plasmonic nanostructures via the observation of their nonlinear optical response and, in turn, to engineer the magnitude and dynamics of nonlinear transients in plasmonic metamaterials via reshaping the excitation absorption profile, and, thus, the electron temperature profile, by engineering of the metamaterial's modal structure. For the metamaterial under consideration, the superposition of the control light-induced absorption profile with that of the probe beam in the ENZ spectral region allows to obtain the fastest transient behavior with a time constant of around 100 fs governed only by electronphonon scattering and electron diffusion. The rich spectrum of the modes of the metamaterial slab provides opportunities to engineer the sign and the temporal response of the nonlinearity in transmission and reflection or for waveguided modes by selectively coupling signal and control light to the modes of the metamaterial slab.
The relaxation time of the optical response monitored through the observation of the transient reflection and extinction, while linked to the relaxation time of the electron temperature distribution within the nanorods, may differ from the latter. The spatial distribution of the electron temperature within the nanorod has an important influence on the metamaterial dynamics at relatively short times following the excitation, when the diffusion term in the TTM is not zero, either accelerating or slowing the relaxation process at different locations within the metamaterial, compared to the uniform electron temperature approximation. The results show the opportunities to control the electron dynamics and, thus, the time response of nanostructured plasmonic systems by tailoring the modal properties of the metamaterial. The engineering of the relaxation dynamics via nanostructuring of plasmonic systems may have many implications in nonlinear optics as well as for hotelectron effects in photophysics and optoelectronics.