Hot carrier spatio-temporal inhomogeneities in ultrafast nanophotonics

Light-induced hot carriers in nanostructures and their corresponding optical nonlinearity have been extensively examined during the last decades. However, nonlinear optical effects dictated by the spatio-temporal evolution of out-of-equilibrium electrons at the nanoscale represent a much more recent research focus. Here we theoretically discuss the role of spatial inhomogeneities that energetic electrons feature across individual nanoantennas in metasurface configuration upon illumination with femtosecond laser pulses. As exemplary cases, we consider two-dimensional geometries of gold meta-atoms having either a high aspect ratio or a tapered cross-section and model their ultrafast optical response. A comparison with numerical results obtained either neglecting or accounting for spatial effects indicates that deep sub-wavelength spatio-temporal transients of carriers may have a significant impact on the dynamics of the all-optically modulated signal, with major quantitative corrections up to predicted changes in sign. Our results present hot-electron local inhomogeneities as an emerging subject with potentially relevant applications in various ultrafast nanophotonic configurations.


Introduction
Recent advances in manufacturing nanostructured optical materials have inspired an extensive research aimed at engineering their properties and implementing functionalities down to the nanoscale, both with metals [1][2][3] and dielectrics [4,5]. In particular, to promote the interaction between matter and light, great efforts have been devoted to the study of nanopatterned periodic structures in quasi two-dimensional (2D) configuration, opening the way to the development of flat optics [6,7]. Indeed, these planar materials, known as metasurfaces, have been shown to serve as flexible platforms for the manipulation of light [8], as well as for nonlinear nanophotonics purposes, offering substantial advantages with respect to conventional set-ups [9][10][11] thanks to their compactness and efficiency. By exploiting strong field enhancements and localised resonances from the individual nanoantennas (meta-atoms), the system optical response can be significantly intensified and readily tailored, opening unprecedented opportunities to achieve advanced functionalities within thicknesses of a few hundreds of nanometers [12][13][14][15][16][17]. Interestingly, because of the peculiar and easy-to-control mechanisms of interaction between light and metasurfaces, the latter have been proven to be especially suitable in all-optical schemes, namely when light is used to both induce and experience changes in the system optical properties [18,19]. Moreover, when femtosecond light pulses are employed to excite such nanostructures, their optical response can be transiently modified on ultrashort time scales [20,21], offering thus novel routes to modulate light and reconfigure optical functionalities with unprecedented speed (see e.g. references [22][23][24][25][26][27][28][29]).
Although results reported so far have explored a variety of approaches, the ultrafast all-optical schemes employed to regulate the metasurface response rely on the onset of giant third-order nonlinear optical effects, whose efficiency is strongly enhanced compared to unstructured thin films thanks to the photonic resonances induced by nanoscale patterning [30][31][32].
Depending on the material under consideration, the third-order optical nonlinearity is governed by different physical processes. As an example, for semiconductor nanostructures, both instantaneous two-photon absorption and delayed nonlinearities due to photoinduced free carriers are involved [26][27][28][33][34][35][36]. For noble metal nanostructures, the χ (3) arises from both (i) free electrons excitation; and (ii) plasmon dephasing and resulting generation [37][38][39] of energetic 'hot' carriers, with the latter mechanism dominating (typically by two orders of magnitude) and exhibiting a delayed and spectrally dispersed behaviour with both real and imaginary contributions [40][41][42]. Details on the nonlinear character of the all-optical modulation of light from metallic structures and related χ (3) formulation are given in reference [41]. Very briefly, the dynamical properties of the permittivity modulation are described by a thermodynamic rate-equation model [43] detailing the temporal evolution of the out-of-equilibrium electrons in terms of their excess energy compared to the equilibrium state and of the main thermalization and scattering processes presiding over their relaxation [44]. Since an excited energy occupation probability of electrons affects the optical properties of the material (mainly its interband absorption), energetic out-of-equilibrium electrons are ultimately responsible for the ultrafast modulation of the Au permittivity triggered by photoexcitation.
Despite the large number of studies on the topic, the temporal evolution of the modulated optical signal driven by out-of-equilibrium carriers has been almost exclusively described upon the approximation that those are photogenerated uniformly on the spatial scale of the individual meta-atoms. This is especially true for plasmonic meta-atoms, in view of the high mobility offered by the metallic medium.
To the best of our knowledge, only a limited number of recent works have focused, either theoretically [45][46][47][48][49] or with experimental evidence as well [50][51][52][53][54][55], on the role of hot-carrier nanoscale spatial inhomogeneities in ultrafast nanophotonics. In particular, it has been shown that hot carrier diffusion, typically requiring hundreds of femtoseconds to lead to homogenization of the excitation, may play a major role in sub-ps light manipulation schemes. Overlooking thus the role of these inhomogeneities in the signal modulation, especially during the first picoseconds, could prevent from accurately predicting the effects on modulation.
In this work, we further develop this argument by showing how the spatio-temporal distribution of photogenerated hot carriers can strongly affect the temporal dynamics of the ultrafast optical response. Some exemplary Au-based structures are investigated upon illumination with ultrashort (100 fs) laser pulses by means of an inhomogeneous version of a wide-spread thermodynamic model for the energy relaxation processes in plasmonic nanoparticles (commonly referred to as the three-temperature model, 3TM [43]). A systematic comparison of the metasurface response predicted when spatial effects are either included or neglected in the model demonstrates that accounting for a non-uniform pattern of the excitation across single meta-atoms may lead to substantial differences in the temporal dynamics of the signal (both in transmission and reflection). Major quantitative discrepancies, up to a sign change, are predicted.
This study indicates that the effects of deep sub-wavelength spatial inhomogeneities featured by out-of-equilibrium carriers may lead to different temporal dynamics of the modulated signal even in rather common structures, with the most significant discrepancies for meta-atoms exhibiting strong field confinement effects. Our results suggest that the approach we pursued could provide the most suitable description of the optical behaviour of nanostructures featuring electromagnetic hot spots and localised field enhancements, especially on the early times (a few ps) following excitation, i.e. the time scale of most relevance for ultrafast nanophotonics.

High aspect ratio nanowire metasurface
To highlight how the spatially non-uniform photogeneration of hot carriers can affect the ultrafast dynamics of the optical response modulation in a metallic metasurface, we first designed a meta-atom with the following characteristics: (i) resonant plasmonic modes with locally well-structured spatial patterns, so to induce significant inhomogeneities of the near fields across the individual nanostructure, and (ii) static extinction enabling to achieve, in the metasurface configuration, non-negligible transmission and reflection. Among relatively simple plasmonic resonators meeting these requirements, a suitable structure is represented by metallic nanowires with high aspect ratio cross-section. Indeed, high aspect ratio nanostructures are known to support higher order plasmon resonances [56] featuring inhomogeneous distributions and high localisation of the induced fields which, as such, can change significantly with the wavelength [51,57]. We therefore investigated the array depicted in figure 1(a), consisting of Au nanowires (lying on a CaF 2 substrate in air) with 1D subwavelength periodicity p = 150 nm, cross-sectional height and width of H = 150 nm and W = 15 nm, respectively. Under p-polarised plane wave illumination (i.e. in-plane electric field, refer to the figure 1(a) inset) at 45 • incidence, the optical behaviour of the system is reported in figure 1(b), showing the unperturbed transmission T 0 (green), reflection R 0 (blue) and absorption A 0 (black) as a function of wavelength. The structure exhibits a well-defined plasmonic resonant feature in transmission, manifesting as a dip at a wavelength of ∼565 nm, which relates to a peak in absorption, blue-shifted by ∼20 nm where A 0 exceeds 60%. Static reflection also features a peak at ∼595 nm, reached after a gradual increase with the wavelength from few percentages to its maximum value of approximately 40%. In general, the static response of the metasurface suggests the spectral region comprised between 530 nm and 630 nm as a suitable working window for all-optical modulation, since both T 0 and R 0 exhibit appreciably high values.
Furthermore, since the absorption pattern of radiation driving the excitation is the element introducing a space dependence in the modulation process, inspecting its distribution as a function of the illumination wavelength is particularly meaningful as the ultimate connection between the electromagnetic problem and the optical perturbation. Starting thus from the local ohmic losses Q loss (r, λ) = 1 2 Re {E · j * D } for an electric field E inducing a displacement current density j D , we computed the local electromagnetic power dissipation density ρ(r, λ) as Q loss per unit cross-sectional surface (in 2D), normalised to its integral over the meta-atom and scaled by a factor equal to the static absorption A 0 (λ), and then varied the excitation wavelength. Figures 1(c) and (d) report ρ(r, λ) (normalised to its maximum) evaluated at a wavelength of λ = 560 nm and λ = 460 nm respectively, revealing indeed major differences.
The former configuration (figure 1(c), for λ = 560 nm sitting on the red edge of the plasmon) exhibits a strong localisation in the lower half of the nanowire as a consequence of the resonant near fields in the same region. This excitation condition has indeed been identified as an optimum balance between the need for high static absorption and significantly inhomogeneous pattern in space. When computing instead ρ(r, λ = 460 nm), shown in figure 1(d), the obtained pattern is substantially more uniform than the one at 560 nm and exhibits much less pronounced effects of localised enhancement across the meta-atom. As such, for a comparable value of absolute A 0 , the spatio-temporal inhomogeneities at 460 nm are expected to affect the response modulation to a lesser degree.
To determine the role of the spatial pattern of ρ(r, λ) on the all-optical modulation of the metasurface, we focussed then on its transient response. We thus implemented a wide-spread description of the ultrafast relaxation processes at the nanoscale following photoexcitation which relies on a well-established rate-equation model, the aforementioned three-temperature model [43] (see also references [58,59]). To elucidate on the role of the spatio-temporal evolution of the hot carriers population, we specifically compared theoretical predictions employing the 3TM either in its homogeneous [59] or inhomogeneous (I3TM, [48,53]) formulation, so to identify the fingerprint of inhomogeneities in the temporal evolution of the modulated signal. The model solves for three internal energetic degrees of freedom within the plasmonic nanostructure to account for the dynamics of hot carriers, which in turn drive ultrafast changes in the metal permittivity. The resulting modulation of the optical signal, governed by the third-order optical nonlinearities of Au (refer to section 5 for details), can then be retrieved in a pump-probe-like interaction: an intense ultrashort light (pump) pulse is shined on the structure to bring it out of equilibrium and, after the excitation, its response dynamical evolution is monitored by the transient absorption of a delayed broadband weaker (probe) pulse as a function of the pump-probe time delay. The three quantities solved for in the model to describe the out-of-equilibrium state of the metal illuminated by the pump pulse are N, the excess energy density stored in a non-thermal fraction of the electron population directly coupled to photoabsorption, Θ E , the temperature of thermalised electrons excited by non-thermal carriers via electron-electron scattering, and Θ L , the lattice temperature modified because of electron-phonon interactions. Its inhomogeneous formulation reads as follows: with a, b the coefficients governing non-thermalised electrons relaxation either via electron-electron or electron-phonon scattering respectively, C E , C L the thermalised electron and lattice heat capacities, κ E , κ L their thermal conductivities, G a thermal electron-phonon coupling term. For a detailed discussion on the coefficients in equations (1)-(3), derived from direct comparisons with more fundamental calculations [60][61][62], refer e.g. to references [53,59] and references therein. Note that here each degree of freedom is a function of position and time (dropped for sake of conciseness). The source term P abs (r, t) includes the electromagnetic power dissipation density previously defined evaluated at the pump wavelength and reading where F is the pump pulse fluence, p the periodicity of the 1D array analysed, S the cross-sectional area of the metasurface unit cell and g(t) the normalised Gaussian intensity profile of the exciting pulse, g(t) = 4 ln 2/πΔt 2 exp −4 ln 2(t − t 0 ) 2 /Δt 2 , centred at t 0 with time duration (full width at half maximum in intensity) Δt. In our study, fluence was set to F = 2 mJ cm −2 and Δt = 100 fs. To retrieve the homogeneous formulation of the 3TM [59], the diffusion terms for temperatures drop and, in the uniform drive term P abs (t), the quantity ρ(r, λ p ) reduces to A 0 (λ p ) after replacing Q loss by its integral over the meta-atom. Note that for the (I)3TM to provide a reliable description of the photoexcitation, the following assumptions are required to be fulfilled: (i) the pump pulse temporal duration is long enough to contain a large number of optical cycles (as it is for ultrashort pulses in the visible range of wavelengths); (ii) the pump fluence is low enough to disregard nonlinear intra-pulse self-actions and to model its absorption in the linear regime. By comparison with self-consistent (yet approximated) models, this assumption has been shown [63] to introduce quantitative discrepancies in the pump absorbed power, with no dramatic changes in the excitation dynamics, mostly because of the typically limited delay effects induced in plasmonic nanostructures.
The main results of our numerical simulations are reported in figure 2, comparing the temporal dynamics of the optical signal obtained with either the inhomogeneous (dark tones) or homogeneous (light tones) formulation of the model. Temporal traces of the differential transmission (green colour scale) , with T(λ, τ ) the probe transmission from the metasurface excited by the pump, T 0 (λ) its static unperturbed value, and reflection (blue colour scale, similarly defined), , are shown at selected wavelengths. The latter have been selected in correspondence with both high signal modulation and pronounced corrections introduced by the I3TM. In figures 2(a)-(f) a pump wavelength of λ p = 560 nm was considered, corresponding to the highly inhomogeneous absorption pattern shown in figure 1(c). The modulated signal is then observed up to a pump-probe time delay τ of 1.5 ps at the probe wavelengths λ of 620 nm (figure 2(a)), 630 nm (figure 2(b)) and 640 nm (figure 2(c)) in transmission, and of 530 nm (figure 2(d)), 540 nm (figure 2(e)) and 550 nm (figure 2(f)) in reflection. For all traces, non-negligible corrections are shown to be introduced by the I3TM (dark curves) with respect to the homogeneous modelling approach (light curves), mainly in terms of amplitude of the modulated signal. This is especially true in reflection (figures 2(d)-(f)), where major quantitative corrections are predicted, with values differing by more than ∼10% for the three wavelengths around the peak of ΔR(λ, τ )/R 0 (λ), up to ∼20% ( figure 2(d)).
To ascertain that these discrepancies can be ascribed to the degree of spatial non-uniformity in the hot carrier population, we performed the same numerical simulations by changing the pump wavelength. Indeed, impinging on the metasurface at λ p = 460 nm induces an electromagnetic absorption pattern driving the excitation (shown in figure 1(d)) which is more homogeneous across the meta-atom. This should have no effects on the temporal dynamics of the signal when the modulation is described with a homogeneous approach, while corrections introduced by a space-dependent description of the optical modulation promoted by hot carriers should be more moderate because of the increased degree of spatial non-uniformity. Results, shown for ΔR(λ, t)/R 0 (λ) in figures 2(g)-(i) for the same probe wavelengths as in figures 2(d)-(f), corroborate our intuition. The signal dynamics obtained either with the 3TM (light curves) or its inhomogeneous variant (dark) are comparatively more similar, and the quantitative discrepancies predicted for the same probe wavelengths when λ p = 560 nm (figures 2(d)-(f)) are significantly reduced. The obtained signal trend with λ p , combined with the computed absorption pattern of the drive term in the 3TM, conclusively confirms that the observed mismatch between model predictions is ultimately due to the spatio-temporal evolution of the hot carrier population. In fact, note that by varying the pump wavelength, slight differences in the differential reflection at fixed probed wavelengths are retrieved also when a homogeneous approach is pursued (refer to light curves in figures 2(d), (g), (e), (h), (f) and (i) respectively). This behaviour however, utterly unrelated to the local pattern of radiation dissipation, can be understood by considering that, in the most general case, different pump wavelengths correspond to: (i) a different amount of energy released by the light pulse and absorbed by the system, which is dictated by the value of A 0 (λ p ); (ii) a different modulation of the interband permittivity term arising from the non-thermal portion of out-of-equilibrium electrons, N, the spectral dispersion of which depends on the pump photon energy (i.e. on λ p ) [64].

Tapered cross section nanowire metasurface
The analysis performed on high aspect ratio meta-atoms has illustrated how spatial inhomogeneities at the nanoscale have an impact on the temporal dynamics of the optical response of a metasurface based on nanostructures with rather common geometry, as the one presented in figure 1. Furthermore, results suggest that the effect benefits from local hot spots and strong field enhancements within the individual meta-atoms, supporting non-uniform patterns of electromagnetic power absorption and subsequent hot carrier generation. Thus, having in mind the relevance of confinement effects in single nanoresonators, we designed a suitable meta-atom to further highlight the role of excitation inhomogeneities at the nanoscale. When resorting to plasmonic systems supporting highly localised field enhancements, the choice may fall on tapered structures. Indeed, tapered plasmonic waveguides are well known to promote the generation of surface plasmon polaritons, which are slowed down when propagating towards the tip as a result of an adiabatic compression, causing energy concentration and a high local field at the tip. The pioneering work reporting such phenomenon [65] was followed by several investigations of the nanofocusing effect from conical structures [66][67][68], enabling energy delivery to a strongly localised region in space [69][70][71][72], of great interest for hot electron exploitation [73]. Moreover, recent studies have demonstrated that illumination of nanocones promotes the generation of temperature gradients at the nanoscale thanks to the unique inhomogeneities featured by the pattern of light absorption [74,75], disclosing the relation between optical excitation and spatial inhomogeneities in tapered structures.
In order to exploit the strong nanofocusing and high field localisation at the tip of cone-shaped nanostructures, the metasurface investigated is a subwavelength 1D array of Au nanowires with tapered cross-section. Height H and base width W of the wire cross-section are respectively 525 nm and 150 nm, while the array periodicity (centre-to-centre distance between neighbouring plasmonic wires) is set to 200 nm, so to avoid higher diffraction orders in the system response, in accordance with the metasurface paradigm. A tilted optical excitation (at 45 • ) with a linearly p-polarised plane wave (electric field in plane) is considered, since an angle of incidence different from zero (normal incidence) entails an electric field with odd symmetry with respect to the wire central axis, as required to excite the slow plasmonic waves capable of high field enhancement at the tip. Figure 3(a) reports the spectra of transmission (green), reflection (blue) and absorption (black) of the structure in unperturbed conditions over a broad spectral range. A typical plasmonic peak in absorption is observed at ∼580 nm, corresponding to reflection approaching zero and eventually vanishing, before increasing up to a broad peak centred at around ∼670 nm. The metasurface transmission, becoming negligible at wavelengths shorter than 550 nm, increases then with wavelength, up to ∼80% at 850 nm. Generally, the three quantities feature relatively high values over the broad spectral region from 610 nm to 850 nm, where effective all-optical modulation may be envisaged in a pump-probe spectroscopy scheme.
To better underlining the effect of spatial inhomogeneities on the metasurface response modulation, we aimed at identifying the most inhomogeneous pattern of absorption across the individual meta-atoms, to select accordingly the pump wavelength. In the spectral range under consideration, a highly non-uniform distribution is predicted precisely at the plasmonic resonance peak, λ = 580 nm, which moreover ensures the highest absorption of the impinging pulse (A 0 exceeds 80%, see figure 3(a)), a further reason to set the pump at 580 nm. The obtained spatial distribution of absorption at the selected wavelength is depicted (normalised) in figure 3(b), also reporting a sketch of the illumination conditions for both the pump (same specifics as in figure 2) and the probe pulses. The spatial pattern of ρ(r, λ) exhibits indeed the expected highly localised enhancement of absorption close to the tip of the wire, thanks to the field confinement in such a small spatial region within the meta-atom.
Following the numerical procedure previously discussed, we could then determine the evolution in space and time of the internal energetic degrees of freedom detailing the energy relaxation processes in the case of the tapered meta-atom. In particular, figures 3(c) and (d) show the temporal evolution of N (figure 3(c)) and ΔΘ E (figure 3(d)) evaluated at three exemplary points across the plasmonic nanostructure. As a result of the strongly inhomogeneous absorption of the pump, the electronic population is brought out of equilibrium with a non-uniform pattern, with the highest excitation levels close to the tip, as confirmed by the (red) curves for both non-thermal carriers (figure 3(c)) and thermalised electrons (figure 3(d)), substantially higher than the curves evaluated at the middle (brown) and bottom (yellow) points of the meta-atom. Moreover, the space-dependent variables computed with the I3TM can be directly compared to the solution of the homogeneous formulation of the nonlinear dynamical model, which provides instead a single-valued time-varying quantity for the whole structure, reported as black curve in figure 3(c) for N and figure 3(d) for ΔΘ E respectively. For the case of N, quantitative corrections are predicted when the I3TM is employed. A more complex scenario should instead be considered for the electronic temperature, non-uniformly excited by N and subsequently undergoing a diffusion process in space across the meta-atom. Indeed, by comparing the four curves in figure 3(d), note that traces do not simply differ by a scaling factor, but follow different dynamics, dictated by the local space-time evolution of carriers in the plasmonic nanostructure.
By then determining the spatio-temporal optical perturbation driven by out-of-equilibrium carriers and interrogating the metasurface with a probe pulse, p-polarised (in-plane electric field), impinging with an angle of 45 • at a time delay τ with respect to the pump, the corresponding ultrafast modulation of the system optical response is retrieved, either accounting for or neglecting spatial effects in the 3TM. To quantify the discrepancy between the two model formulations, we introduced the quantities δ T , δ R and δ A as a measure of the correction introduced by the I3TM with respect to the more standard homogeneous approach, namely δ T (λ, τ ) = T inhom. (λ, τ ) − T homog. (λ, τ ), with T inhom.(homog.) the dynamical transmittance ensuing from the inhomogeneous (homogeneous) simulation. Likewise, δ R (λ, τ ) and δ A (λ, τ ) are similarly defined for dynamical reflection and absorption.
To assess the impact of a non-uniform description of photoexcitation on the optical signal in terms of these quantities, figure 4(a) shows temporal sections of the three δ up to a pump-probe delay of 1.5 ps for some selected wavelengths within the probe spectrum, λ 1 = 610 nm, λ 2 = 630 nm and λ 3 = 730 nm.
The so-defined correction terms, as well as being wavelength-dependent, are different both in amplitude and dynamics (compare curves referring to the same wavelength in figure 4(a)), however they all reach values as high as some percent for all cases, thus affecting T, R and A to the same extent. Moreover, the fact that traces evolve over time and change in sign demonstrates that the I3TM does not simply predict a constant off-set in the system response compared to the homogeneous 3TM. Rather, by providing a non-uniform description of the optical modulation, the I3TM includes mechanisms that are dynamically involved in the relaxation process. This manifests more notably at early times following pump absorption, i.e. at the onset of the modulation (see curves in figure 4(a) during the first hundreds of femtoseconds time delay), suggesting a dominant role of non-thermal carriers to the predicted discrepancy.
The weight of the corrections shown in figure 4(a) can be better appreciated when normalized to the static optical response of the metasurface, as detailed in figure 4(b), showing the map of δ T (λ, τ )/T 0 (λ). Similarly, figure 4(c) reports the same quantity evaluated in reflection, i.e. δ R (λ, τ )/R 0 (λ). Inspecting the maps of the normalised quantities reveals indeed that inhomogeneous predictions of modulation differ from homogeneous results over a broad range of wavelength (spectra are shown from 610 nm to avoid spurious effects due to the drop of T 0 and R 0 below ∼605 nm, as shown in figure 3(a)). Also, note that, even for the moderate pump fluences under consideration, the discrepancy between inhomogeneous and homogeneous modelling, measured by the δ T /T 0 and δ R /R 0 figure of merits, reaches peak values as high as ∼10% (in modulus), which is comparable to the ΔT/T 0 and ΔR/R 0 transient optical signal.
To highlight this issue, figure 4(d) shows the temporal dynamics of the ΔR(λ, τ )/R 0 (λ), evaluated at three wavelengths (the same as in figure 4(a)) and obtained by implementing either the I3TM (dark shade) or its homogeneous formulation (light shade). Pump-probe traces unveil notable aspects in the reflection modulation. While the mismatch between models is almost uniquely quantitative for λ 3 (compare bottom curves in figure 4(d)), the ΔR/R 0 at λ 2 = 630 nm (middle curves in the panel) follows different dynamics according to the model employed, and traces are not a scaled version of each other. Besides the difference in amplitude, the rise time of the modulated signal in the case of an inhomogeneous approach is almost doubled (from ∼400 fs to ∼800 fs). Such delay effect in reaching the peak is indeed assigned to spatio-temporal inhomogeneities taking place at the nanoscale [55] and explained as the fingerprint of a homogenization process undergone by excitation across the individual meta-atoms, which modifies the dynamics of the optical excitation. However, the most striking difference between model predictions is observed at λ 1 = 610 nm (top curves in figure 4(d)). Consistently with the sign and values of δ R (refer to figures 4(a) and (b)), the homogeneous approach underestimates the I3TM results throughout the entire time scale analysed, and a major quantitative mismatch in the reflection modulation is obtained. Yet, in this case the ΔR/R 0 computed via the I3TM undergoes a change in sign, from positive (during the first hundreds of femtoseconds, up to ∼370 fs) to negative values (at longer time delays up to complete relaxation of the system), an utterly absent feature in the signal temporal trace retrieved by the homogeneous model. Such an unexpected behaviour can be ascribed to the contribution to the transient signal arising from non-thermal carriers, which is reported in figure 4(e). Note that at λ 1 this contribution is completely overlooked by the homogeneous model because of the dilution of the excitation over the entire domain of the meta-atom (light trace in the top panel of figure 4(e), corresponding to the 3TM simulation, is orders of magnitude lower than the dark trace, corresponding to the I3TM). However, spectral dispersion of the permittivity modulation also plays a key role in this picture, because when considering a different probe wavelength, e.g. λ 2 , the homogeneous 3TM turns out to overestimate the contribution from non-thermal carriers (light trace in the middle panel of figure 4(e), corresponding to the homogenous 3TM, dominates over the dark trace, corresponding to the I3TM calculation). At λ 3 the two models retrieve similar estimation of the ΔR/R 0 from non-thermal carriers (light and dark traces in the bottom panel of figure 4(e) are comparable to each other), and thus the corrections arising from spatial inhomogeneities are only quantitative, with no substantial impact on the temporal dynamics.

Conclusions
In summary, this study has illustrated the effects that the spatio-temporal evolution of light-driven hot carriers may have on the dynamics of the ultrafast all-optical modulation in plasmonic metasurfaces. Indeed, although great interest has been devoted in the last years to explain how the optical signal varies over time following excited electron relaxation, the role of spatial inhomogeneities on the hot-carrier driven optical nonlinearity represents a much more recent research topic. Our analysis focussed on two exemplary Au-based metasurfaces featuring either a high aspect ratio or a tapered cross-section of the meta-atoms, and accounted for the transient photoexcitation via the well-established three-temperature model by proposing an inhomogeneous formulation (the I3TM) to include space-related effects in the processes of hot carrier excitation and subsequent permittivity modulation within the individual meta-atoms, i.e. on a deep sub-wavelength scale. We theoretically predicted major quantitative corrections, being of the same order of magnitude as the transient optical signal, and even changes in sign, due to hot electrons spatial inhomogeneities, when compared to a homogeneous modelling approach, disclosing the impact of the spatio-temporal evolution of the hot carrier population on the dynamics of pump-probe traces. Moreover, our comparative study on the two selected metasurfaces seems to disclose the crucial role of the absorption pattern experienced by the pump pulse, imprinting a characteristic spatial distribution of the hot carrier population which then evolves in time following diffusion processes across the individual meta-atom. More pronounced discrepancies between predictions from homogeneous versus inhomogeneous model are indeed retrieved when the metasurface supports local hot spots and strong field confinement in the metallic domain. This condition, as shown by the structures investigated, can be achieved by either properly tuning the excitation wavelength or carefully designing the meta-atom. Our results indicate that local inhomogeneities of hot electrons in nanostructures can indeed introduce substantial corrections in extended systems featuring localised field enhancement and non-uniform modal patterns under illumination. This is particularly true at early times following photoexcitation, representing the most relevant time window for applications in ultrafast nanophotonics, for which an accurate description of the evolution of the hot-carrier driven optical perturbation is crucial.

Numerical modelling
The numerical simulations of the processes involved in the photoexcitation of the plasmonic metasurfaces here considered have been conducted via a finite element method-based commercial software (COMSOL Multiphysics 5.6), employed both to describe the electromagnetic behaviour of the (either unperturbed or excited) structure and to determine the dynamical evolution of the degrees of freedom of the I3TM (the same numerical methods apply to the 3TM, after having adapted the expression of the drive term as outlined in the main text). To first compute the optical behaviour of the metasurface, full-wave simulations solve for the scattering problem in the frequency domain. A parametrically swept monochromatic plane-wave illumination is considered, and periodic boundary conditions in the port formalism are set in order to compute the near fields across the unit cell of a perfectly periodic array. The substrate (CaF 2 ) refractive index is set to 1.43, while the permittivity of gold in unperturbed conditions ε 0 Au (λ) is retrieved from the analytical model proposed in reference [76] fitted on experimental data [77]. To determine the transient response of the system, a segregated approach has been pursued [53], starting from the calculation of ρ(r, λ p ) in static conditions, to be used in the drive term of the I3TM, subsequently solved in the time domain. Once the dynamics of the three internal degrees of freedom have been computed, frequency-domain (i.e. treating time as a parameter) simulations are performed to evaluate the interaction of the perturbed metasurface with a probe pulse at a fixed time delay τ from the pump pulse. Formally, the numerical model is the same as the one implemented in the first step (static, unperturbed, conditions), for a modified permittivity in the metallic meta-atom ε 0 Au (r, τ , λ) = ε 0 Au (λ) + Δ ε(r, τ , λ), where the space-time modulation of permittivity Δε(r, τ , λ) is computed starting from the solution of the I3TM via a semi-classical model of the thermo-modulational nonlinearities in Au [78][79][80]. In brief, both N and ΔΘ E entail a modulation of the interband transitions in the metal, because of the modifications of the electronic distribution f(E) around the Fermi level. For the case on non-thermal electrons, the change in f(E) is a double-step-like function with width in energy set by the incoming photon energy, while an increased electronic temperature causes a smearing of the Fermi-Dirac distribution. Due to those changes in the energy distribution of out-of-equilibrium electrons, the probe beam impinging at a fixed delay time experiences a modified absorption profile, accounted for a modified imaginary part of the permittivity. The corresponding real part of the variation is retrieved by Kramers-Kronig relations. In addition, the intraband term of Au permittivity is also modified throughout photoexcitation because of changes in the Drude damping factor and plasma frequency, both mostly dictated by an increase in the Au lattice temperature, Θ L . A combination of those effects defines thus the total transient complex-valued permittivity modulation Δε(r, τ , λ). Further details on the processes involved can be found in [48,53,59,81].