Light-driven nanoscale vectorial currents

Controlled charge flows are fundamental to many areas of science and technology, serving as carriers of energy and information, as probes of material properties and dynamics1 and as a means of revealing2,3 or even inducing4,5 broken symmetries. Emerging methods for light-based current control5–16 offer particularly promising routes beyond the speed and adaptability limitations of conventional voltage-driven systems. However, optical generation and manipulation of currents at nanometre spatial scales remains a basic challenge and a crucial step towards scalable optoelectronic systems for microelectronics and information science. Here we introduce vectorial optoelectronic metasurfaces in which ultrafast light pulses induce local directional charge flows around symmetry-broken plasmonic nanostructures, with tunable responses and arbitrary patterning down to subdiffractive nanometre scales. Local symmetries and vectorial currents are revealed by polarization-dependent and wavelength-sensitive electrical readout and terahertz (THz) emission, whereas spatially tailored global currents are demonstrated in the direct generation of elusive broadband THz vector beams17. We show that, in graphene, a detailed interplay between electrodynamic, thermodynamic and hydrodynamic degrees of freedom gives rise to rapidly evolving nanoscale driving forces and charge flows under the extremely spatially and temporally localized excitation. These results set the stage for versatile patterning and optical control over nanoscale currents in materials diagnostics, THz spectroscopies, nanomagnetism and ultrafast information processing.

such approaches have also remained limited to single laser spots or larger micrometer-scale structured light fields.
Plasmonic systems can be utilized to overcome these spatial limitations by concentrating light down into deeply sub-diffractive nanometer scales.Such systems are already known to provide exceptional control over energy flow between photonic, electronic, and thermal degrees of freedom on ultrafast timescales 20 .A better understanding of momentum flow in plasmonexcited hot carrier distributions is also emerging 21 , as determined by spatially-tailored plasmonic hot spots and directionality imposed by nanoscale geometry.Although such effects have primarily been explored via nonlinear photoemission into free space [22][23][24] , linear photocurrent responses under continuous-wave excitation have been observed recently in hybrid plasmonic systems serving as bias-free mid-infrared photodetectors 25,26 as well as spin-valley polarized valleytronic transistors 27 .
Here, we show that plasmonic metasurfaces offer far broader capabilities for harnessing charge flow at nanometer spatial scales and femtosecond time scales.Photocurrents are a universal manifestation of broken inversion symmetry 3 , and we find that asymmetric gold nanoantennas on graphene exhibit strong light-driven directional responses.Local directionality is determined by the orientation of individual nanoantennas, with general implications for global spatially-varying and optically-controlled photocurrents.As an immediate application, we demonstrate that these vectorial optoelectronic metasurfaces serve as efficient and versatile sources of ultrafast THz radiation, including broadband THz vector beams.Electrostatic gating and multiphysics modeling reveal a local photothermoelectric driving mechanism and elucidate previously unexplored dynamics occurring at the intersection of femtosecond excitation and nanoscale localization.

Generation of Ultrafast Directional Photocurrents
Our optoelectronic metasurface concept is illustrated in Fig. 1a, where inversion symmetrybroken gold nanoantennas act as lightning rods 28 with strong resonantly-enhanced light fields at their sharp tips (~15 nm radius; Fig. 1b,c insets).Hot carrier excitation at these nanoplasmonic hot spots drives local current density (, ) within the underlying graphene.The broken inversion symmetry leads to net current along the orientation of the nanoantennas, which decays on subpicosecond timescales due to momentum-relaxing scattering with optical phonons, impurities, and substrate phonons [29][30][31] .These transient photocurrents radiate free-space THz waves,  ∝ −/, which are utilized as a contact-free probe of the ultrafast current dynamics.Additionally, the time-averaged DC photocurrent is read out directly via electrical contacts for unambiguous verification of the overall charge flow behaviors.
We first investigate the photocurrent response in metasurfaces with uniform nanoantenna orientation and thus a globally preferred directionality.The plasmonic resonance is broadly tunable from visible to infrared wavelengths based on the nanoantenna structure (shape and size), material, and dielectric environment.Structural tuning is exploited here for resonances at 800 nm and 1550 nm, as clearly observed in the measured and simulated transmission spectra (Fig. 1b).
Corresponding enhancements of the ultrafast THz fields and DC photocurrents are shown in Fig. 1c, in good agreement with the simulated field intensity enhancements, which underscores the central role of the (asymmetric) plasmonic fields in the photocurrent generation process.
Single-cycle THz pulses are observed from both the 800 nm and 1550 nm (Fig. 1d) metasurfaces via free-space electro-optic sampling (Methods), indicating the photocurrent rise and decay on few-hundred femtosecond timescales.These pulses disappear when the exposed graphene is etched away without removing the nanoantennas (Supplementary Note 1), verifying that the THz radiation originates from graphene photocurrents rather than optical rectification directly at the nanoantennas 32 .Remarkably, with little optimization (Supplementary Note 2), the 30-nm-thick metasurfaces yield THz fields that are comparable to those from widely-utilized 1mm-thick 110 ZnTe nonlinear crystals under the same excitation conditions.This indicates the potential utility of these optoelectronic metasurfaces as new THz sources with versatile responses arising from sample orientation and patterning as well as incident polarization and wavelength, as explored further below.
A linear dependence on incident fluence (Fig. 1e) is observed in both the THz field amplitude and photocurrent readout for the 800 nm metasurface below ~0.8 µJ cm −2 , which can result from various possible mechanisms at metal-graphene junctions, such as photovoltaic 14,33 or photo-Dember 34 effects.However, photothermoelectric effects are generally found to prevail, including under continuous-wave excitation at the nanoscale 25,26 and ultrafast excitation at mesoscopic metal-graphene junctions 35,36 .We will later show that photothermoelectric effects remain dominant under simultaneous ultrafast, nano-localized excitation, but we first explore the general implications of artificially patterned linear photocurrent responses, which can be extended to other hybrid material systems and are independent of the specific mechanism.

Local Vectorial Photocurrents and Global Responses
For metasurfaces with uniform spatial patterning as shown in Fig. 2a,b, global current and THz emission measurements reveal the local light-matter interaction symmetries and photocurrent responses at the nanoscale unit cell level.The oriented metasurface (pm wallpaper group, Fig. 2a) exhibits a particularly simple cos () dependence on the incident linear polarization angle with respect to the nanoantenna principal axis of broken symmetry (Fig. 2c), which is characteristic of the linear dependence on the projected field intensity.The photocurrent directionality remains along this axis, regardless of the incident light polarization.Under circularly-polarized incident light, the overall nanoantenna-oriented photocurrent across the metasurface is revealed by the spatially-uniform linear polarization of the radiated THz beam, which is continuously rotated by sample rotation (Extended Data Fig. 1).
Additional symmetry can be introduced into the metasurface unit cell with multiple nanoantennas oriented in different directions.In general, a unit cell with Cn rotational symmetry will exhibit an isotropic linear response if n > 2, with noncentrosymmetric configurations (odd n) required for net directional photocurrent.We demonstrate this with a Kagome lattice metasurface (Fig. 2b), which has three-fold rotational symmetry and three planes of mirror symmetry (p3m1 wallpaper group).Global and corresponding local photocurrent directionality is continuously rotated by varying the incident linear polarization angle (Fig. 2d,e), with a linear combination of hot spot excitations (Fig. 2b) maintaining a constant photocurrent magnitude.This omnidirectional current control is described in a Cartesian basis by  =  ∑ cos ( −  ) cos( ) and  =  ∑ cos ( −  ) sin( ) , where  is the laser polarization angle and  = 0, 2/3, and 4/3 for the three nanoantenna orientations within the Kagome unit cell.
These prototypical examples already illustrate a broad design space for nanoscale vectorial photocurrents, which can be driven in arbitrary directions in space and/or time under different incident light fields.Localized vortical charge flows with no net global current can be realized under circularly-polarized excitation, for instance, by rotating the nanoantenna orientations within the Kagome metasurface unit cell (Fig. S11).The design space grows further with the introduction of different structures and resonant wavelengths.In such systems, the light-matter interaction symmetry may even be different than the underlying structural symmetry and tunable with wavelength, as read out by the polarization-and wavelength-dependent photocurrent response.
Such possibilities are demonstrated in Extended Data Fig. 2, where the Kagome metasurface unit cell contains three resonators of different dimensions.Although this perturbed Kagome lattice evidently exhibits the lowest possible degree of symmetry (p1 wallpaper group), an approximate mirror symmetry in the light-matter interaction appears either along a dominant resonance axis or between two equally excited resonances at intermediate wavelengths 24 .The general linear photocurrent response for any metasurface with uniform spatial patterning is given by  (, ) ∝ ∑  () cos ( −  ) cos( ) and  (, ) ∝ ∑  () cos ( −  ) sin( ) , where  () are the wavelength-dependent plasmonic fields at each tip hot spot.This shows a high degree of adaptability, with responses that are not rigidly fixed by metasurface patterning.These results also directly demonstrate important capabilities for bias-free polarization-resolved (Fig. 2d,e) and wavelength-resolved (Extended Data Fig. 2) photodetection for integrated ultrafast optoelectronic and information science applications 37 .

Spatially-Varying Vectorial Photocurrents
Nonuniform spatial patterning can be implemented for nearly arbitrary control over vectorial current distributions out to macroscopic scales, as demonstrated with radial (Fig. 3a) and azimuthal (Fig. 3e) metasurfaces.Corresponding radial and azimuthal transient charge flows generated upon uniform circularly-polarized excitation are expected to be evident in the far-field THz radiation patterns.We thus utilize a high-throughput THz emission spectroscopy system for hyperspectral imaging (Fig. 3d), scanning across the THz beam and reading out the full time, frequency, and polarization information as a function of position (Methods).These maps reveal clear radially-polarized (Fig. 3b,c) and azimuthally-polarized (Fig. 3f,g) THz fields, unambiguously confirming the vectorial photocurrents.For the radial metasurface device, the radial current and resonantly-enhanced excitation are also verified by direct photocurrent readout measurements (Extended Data Fig. 3).
While the mapped THz fields serve to demonstrate the spatially-varying vectorial currents, they also directly illustrate a new capability for direct generation of broadband vector beams in the THz frequency range (out to ~4 THz; Fig. 3d).These cylindrical THz vector beams hold significant promise for THz imaging, spectroscopy, and other applications 16 .Although various schemes have emerged recently for generating THz vector beams 16 , they often rely on bulky, narrowband, and/or multi-stage systems.By contrast, direct generation of arbitrary THz vector fields from spatially-patterned photocurrents offers an ultrathin, broadband, and single-stage source that can also be actively manipulated by incident optical fields or by electrostatic gating, as examined next.

Photothermoelectric Dynamics
We now investigate the rich photothermoelectric dynamics that occur within the goldgraphene metasurfaces under femtosecond optical excitation at the nanometric tip hot spots.The photocurrent driving mechanism is revealed by electrostatic gating of large-area metasurface devices (Methods), with a back-gate voltage Vg applied to tune chemical potential µ and corresponding electrical conductivity  in the bare graphene regions (Fig. 4a).By contrast, the chemical potential of the graphene beneath the nanoantennas remains pinned to that of gold 38 , leading to a spatially-varying Seebeck coefficient (Fig. 4c) described by the Mott relation at room temperature  ,  = − , where kB is the Boltzmann constant, e is the elementary charge, and || ∝  −  / with  ≈ 9 V at the charge neutrality point (Supplementary Note 1).
The photocurrent closely follows the nonmonotonic dependence of the difference Δ between the bare and gold-doped graphene regions (Fig. 4b), which is a clear signature of a photothermoelectric effect 39,40 and is inconsistent with photovoltaic (monotonic Vg dependence) or other linear responses that can occur at graphene junctions 14,33,34 .
We therefore combine electromagnetic, thermodynamic, and hydrodynamic modeling to understand the evolution of the far-from-equilibrium electronic temperature  (, ), the effective photothermoelectric driving field ∝ −(, )∇ , and the resulting nanoscale charge flow (, ).
The localized optical power absorption generates a peak Te approaching 3000 K at the nanoantenna tip (Fig. 4d) for an incident pulse fluence of 0.5 µJ cm −2 in the linear response regime.In-plane electronic thermal diffusion and scattering with optical phonons 29 along with out-of-plane scattering with substrate phonons 30 then cool the superheated carrier distribution within a few hundred femtoseconds (Supplementary Note 3).Before cooling, a net acceleration proportional to Δ∇ acts on the electronic system along the nanoantenna axis (Fig. 4e).We note that while ΔS exhibits a nonlinear temperature dependence at high Te (Supplementary Note 3), this is not expected to alter the gating behavior shown for Δ (Fig. 4b) beyond a scaling factor.The resulting charge flow is modeled via the time-dependent incompressible Navier-Stokes equation (Fig. 4f; Methods).Such a hydrodynamic framework [43][44][45] is based on the condition that  ≪  , where  is the electron-electron scattering time and  is momentum relaxation time associated with electron-impurity/phonon scattering.This treatment can be justified even for relatively lowmobility graphene ( ≈ 45 fs) by the electronic superheating, which dramatically increases the phase space for electron-electron scattering.Upon photoexcitation,  ∝  is driven from a few hundred femtoseconds at room temperature into the few-femtosecond range at several thousand Kelvin (Supplementary Note 3), supporting a hydrodynamic treatment and the possibility of lightinduced viscous electron flow 46 .
Indeed, our multiphysics modeling reveals the rapid evolution of Te,  , and corresponding local physical quantities (e.g., electron viscosity, thermopower, and heat capacity) that occurs over hundreds-of-femtosecond timescales and tens-of-nanometer spatial scales (Extended Data Fig. 4).
The calculations also capture the essential macroscopic behaviors measured for the various metasurfaces, including polarization, frequency, and intensity dependencies, with time-averaged photocurrents that agree to within a factor of 2 of experimental measurements (Fig. S9).
These robust model results elucidate the ultrafast nanoscale charge flows that produce net local and global vectorial photocurrents measured via THz emission and electrical readout.Beyond gold and graphene, the choice of materials for vectorial optoelectronic metasurfaces is expansive, with photocurrent responses expected due to photovoltaic, photoelectric, and photothermoelectric effects in semiconductor, semimetal, topological insulator, ferromagnetic and other material systems.Symmetry-broken nanoscale structuring may therefore also be imposed on new materials to probe unknown local responses and physical properties 3 .Furthermore, the introduction of planar chirality and thus net vorticity can lead to dynamic nano-magnetic moments (Supplementary Note 4).More generally, spatially-varying photocurrents can be utilized to generate complex magnetic fields 47 in free space or nearby materials.While we have directly demonstrated polarization/wavelength-resolved photodetection and ultrafast THz generation in symmetrybroken optoelectronic metasurfaces towards information science, microelectronics, and THz technologies, we anticipate that many more opportunities will emerge in other applications (e.g., nano-magnetism, materials diagnostics, analog computing, and energy harvesting) using the basic concepts introduced in this work.

Methods Metasurface Fabrication
Metasurface fabrication begins with large-area monolayer graphene grown via chemical vapor deposition (CVD) and transferred onto fused quartz substrates (Supplementary Note 1).
Substrates are spin-coated with bilayer PMMA (PMMA 495 bottom layer and PMMA 950 top layer) at 3000 rpm for 30 s, followed by 1 min prebake at 180 °C.A conductive polymer layer (DisCharge) is then applied to ensure adequate charge dissipation during electron beam exposure.
Arrays are then patterned via electron beam lithography using a 100 kV JEOL JBX-6300FS system.For < 30 nm features, interparticle proximity effect correction is utilized with an optimized total dosage around 1650 µC cm −2 .Under these conditions, writing of a 1 mm 2 array takes approximately 1 hr.Following a 30 s rinse in isopropanol (IPA) for DisCharge removal, samples are then developed in a 3:1 IPA:MIBK (methyl isobutyl ketone) solution for 1 min.A 30 nm gold layer is then deposited via electron-beam or sputter deposition without any adhesion layer.Liftoff is performed after overnight acetone soaking, using gentle acetone and IPA wash bottle rinsing to remove the residual PMMA/gold.For most samples, no particle delamination is observed.
Devices for simultaneous ultrafast THz emission and average photocurrent electrical readout studies are prepared using maskless photolithography (Heidelberg MLA 150).We utilize AZ 5214E photoresist in image reversal mode, precoating the samples with a hexamethyldisilazane (HMDS) layer to promote photoresist adhesion.Excess graphene is etched to obtain desired geometries via 2 min exposure to O2 plasma (100 W, 10 sccm; Anatech RIE).Electrodes consist of 50 nm Au with a 3 nm Ti adhesion layer.
For back-gated devices, the back-gate electrode with minimal overlap with the top electrodes (precluding any shorting through defects in the dielectric spacer layer) is prepared via photolithographic patterning, deposition of 3 nm/30 nm Ti/Pt, and a liftoff process.Various spacer layer materials (SiO2, Al2O3, and HfO2), deposition methods (atomic layer deposition and physical vapor deposition), and thicknesses (5-100 nm) are tested for optimal electrostatic gating of largearea devices.We find that 30 nm SiO2 offers a good trade-off between high gating capacitance and relatively modest nanoantenna resonance red-shifting due to the image charge oscillation within the Pt and dielectric environment of the spacer layer.Subsequent device preparation follows as described above.

Basic Characterization
The Metasurface optical properties are verified via white light transmission spectroscopy using a tungsten-halogen white light source passed through a broadband polarizer for linear polarization control.Two 20× microscope objectives are utilized to focus the light through the sample and collect the transmitted light into both visible (Acton SP2300) and near-infrared (Ocean Optics NIRQuest) spectrometers.The transmittance is given by the ratio of signal (on the metasurface) to reference (off the metasurface but on the same substrate), with the ambient background collected with the white light source turned off and subtracted from each.

Terahertz Emission Spectroscopy
THz emission experiments are performed across several systems with different functionality, with overlapping datasets verifying reproducibility.All measurements except wavelength-dependent THz emission (see below) are performed in the low-fluence regime (< 1 µJ cm −2 ), utilizing a Ti:sapphire oscillator (Chameleon Vision-S) providing 800 nm, ~100 fs pulses at the sample location with 80 MHz repetition rate.As seen in Fig. 1e, these low fluences ensure linear responses and also preclude thermally-induced nanoantenna deformation or photochemical degradative effects 48 that can occur in the intense hot spot regions.The gradual onset of a sublinear regime beyond ~0.8 µJ cm −2 in Fig. 1e is attributed to competing temperature-dependent thermodynamic contributions (Supplementary Notes 3 and 4).For metasurfaces operating at 800 nm, The THz radiation is measured via electro-optic sampling 49 using a 1 mm 110 ZnTe crystal.
Three wire-grid polarizers are utilized to measure the x and y field components of the THz radiation, with the first polarizer at 0 or 90, second polarizer at 45, and third polarizer at 0.The THz beam path is contained within a dry-air purged environment (< 2% relative humidity).
Pump-wavelength-dependent THz emission experiments are performed using a noncollinear optical parametric amplifier to generate pump pulses across 725-875 nm (pulse width ≲ 50 fs) and an optical parametric amplifier to generate pump pulses across 1450-1650 nm (pulse width ≲ 200 fs).Both amplifiers are seeded by a 1 MHz fiber amplified solid state laser producing ~270 fs pulses at 1032 nm.In all of these measurements, we tune the pump wavelength while the incident fluence is maintained at a constant value using a neutral density attenuator.The direct output of the seed laser is used to gate the THz emission using electro-optic sampling with a 0.5 mm thick 110 GaP crystal.This enables wavelength-dependent pumping without affecting the detection.This system is also utilized to obtain the THz emission for the 1550 nm metasurface presented in Fig. 1d, right panel.
For spatial mapping of THz vector beams, we utilize a custom-developed Menlo Systems apparatus, consisting of a 100 MHz mode-locked 1560 nm erbium-doped all-fiber laser oscillator that seeds two separate amplifiers.The first line uses a high-power erbium fiber amplifier driven in the linear pulse propagation regime followed by free-space second harmonic generation, outputting 1 W of 780 nm, 125 fs pulses and serving as the pump beam for generation of THz radiation from metasurfaces.The second amplifier line is used to gate a broadband, fiber-coupled PCA detector 50 (TERA 15-RX-PC), which is mounted on a two-dimensional stage for automated xy spatial scanning.Along with the high-repetition integrated delay stage (> 45 Hz for a 20 ps temporal scanning range), this system enables rapid hyperspectral imaging of THz vector beams collimated and focused by two TPX lenses (50 mm focal length).The PCA is situated before the focal point of the second lens, with resulting deviations from a planar phase-front calibrated and corrected for in the results presented in Fig. 3.

Device Gating and Photocurrent Readout
Time-averaged photocurrents are electrically read out via a lock-in amplifier (LIA; for the highest signal-to-noise ratio) and picoammeter (to directly measure polarity).Gating studies are performed using a computer-controlled dual-channel sourcemeter (Keithley 2614B) for simultaneous application of a gate voltage and readout of either the device resistance (via sourcemeter) or photocurrent (via picoammeter or LIA).While nonlocal photocurrents on electrodes are often evaluated in terms of a Shockley-Ramo type response 1,51 , the square metasurface and electrode geometries here lead to a simplified, direct readout of the overall x and y photocurrent components.Wavelength-dependent photocurrent responses are measured by directly tuning the output wavelength of the 80 MHz Ti:sapphire laser while maintaining a constant power.While the electrode gold-graphene junctions can also contribute to the photothermoelectric current, the current induced by the nanoantennas is isolated via (i) the spatial position of the focused beam, (ii) the much stronger nanoantenna polarization dependence, and (iii) the lack of a similar response in nanoantenna-free devices (see Extended Data Fig. 5 for further details).

Dynamical Multiphysics Model
Coupled electromagnetic, thermodynamic, and hydrodynamic equations are solved using the finite element method (COMSOL Multiphysics 6.1).Nanoscale optical fields are first simulated The full thermal evolution of the 2D graphene electron and lattice (specifically, optical phonon 29 ) systems is calculated through the coupled heat equations, where ce is the volumetric electronic specific heat described in Supplementary Note 3,  =   is the in-plane electronic thermal conductivity approximating Wiedemann-Franz law behavior (approximately preserved in the Fermi liquid regime of graphene 45 ),  = 3 MW m −2 K −1 is utilized as the out-of-plane electronic thermal conductance due to coupling with the SiO2 substrate phonons 30 (corresponding to a cooling length  / ≈ 150 nm at 3000 K), Top is the optical phonon temperature, cop is the optical phonon specific heat based on time-resolved Raman studies of graphite 29,55 , and  ≈ 1.7 × 10 W K −1 is the lattice thermal conductance based on thermal transport measurements of supported graphene 56 (1:1 device aspect ratio, here).The     Extended Data Fig. 3 | Radial metasurface device and resonant excitation.a, Optical micrograph of a radial metasurface, prepared as a device for photocurrent readout.Inset: Radially-outward-oriented nanoantennas with resonance at 800 nm and nanoantenna spacing ~500 nm.b, Measured DC photocurrent (dashed line) and simulated plasmonic hot spot field intensity (solid fill; same as in Fig. 1c), verifying the resonantly-coupled/enhanced radial photocurrent.
Extended Data Fig. 4 | Ultrafast time evolution of local thermodynamic and hydrodynamic quantities.Calculated electron temperature, photothermoelectric acceleration field, hydrodynamic flow profile, and relative scattering rates within the square metasurface unit cell, as a function of time with respect to the 100 fs laser pulse.All plots are normalized to the peak values.
Extended Data Fig. 5 | Electrostatically gated metasurface.a, Metasurface device, with false coloring indicating the 3 nm/50 nm Ti/Au (light yellow), nanoantenna array (darker yellow), graphene (dark gray), 30 nm Pt back-gate electrode with 30 nm silica spacer layer (light gray), and substrate (light blue).b, Measured (solid line) and simulated (solid fill) reflectance spectrum of device, showing resonant absorption around 875 nm.c, Gate-dependent photocurrent as a function of position, with locations indicated in (a).Beam spot ~25 µm.To isolate the electrode contributions, the left and right electrode measurements are performed on an array-free device, with little laser polarization dependence observed.A polarity reversal is observed between the left and right electrodes.By contrast, strong polarization dependence is observed on the array at the center of the sample.Both metasurface and array-free devices have 250 µm square graphene between the electrodes (metasurface device is same as in Fig. 4).Photocurrents are normalized to the peak positive or negative values.
1.2 Nanoantenna structure and effect of exposed graphene removal Gold nanostructures are fabricated with 15 nm radius of curvature tips, with good sticking between the gold and graphene obviating the need for a (dissipative) metal adhesion layer.To test the role of the graphene currents in THz generation, we remove the exposed graphene in a metasurface device via a minimally destructive 45 s plasma etch (100 W, 10 sccm O2; Fig. S2a,b).This causes a small blueshift attributed to slight tip deformation during the plasma treatment (Fig. S2c,d) but maintains the nanostructure integrity and strong plasmonic resonance (similar quality factor).Graphene can thus serve as an atomically-thin, minimally-damping adhesion layer between the gold and the underlying fused quartz substrate (or, in general, a variety of materials upon which graphene can be readily transferred).Most importantly for the present studies, the conductivity drops to zero after etching away the bare graphene and no photocurrents are observed.The disappearance of the THz radiation (Fig. S2e,f) under otherwise identical testing conditions illustrates the central role of the graphene photocurrents, as opposed to any local currents are optical rectification within, at the surface of, or beneath the gold nanostructures.

Supplementary Note 2: Electromagnetic response 2.1 Preliminary investigation of lattice density and coupling effects
The resonance of plasmonic nanoantennas arranged in a regular array-with spacings comparable to the incident wavelength-is influenced by intermediate-/far-field coupling between resonators.This is known as the surface lattice resonance 4 , which can be tuned to influence the resonance wavelength and Q factor compared with isolated nanostructures.Such effects are accounted for here in the electromagnetic simulations via periodic boundary conditions.Rather than iteratively optimize the resonator design, array geometry, and spacing for an 800 nm (or other desired) resonance, we simply optimize the resonator design in a 500 nm square lattice then test the effect of varying density (Fig. S3a).A shift in resonance position and quality factor is observed for the varying lattice constants, in good agreement with simulations (Fig. S3b).The 500 nm pitch metasurface exhibits the highest tip intensity enhancement (Fig. S3c) and THz emission (Fig. S3d) and is thus utilized for the present studies of the uniformly-oriented arrays.However, we note that significant further optimization is likely to be possible with different resonator designs, lattices, and nanoantenna spacings.
Because the same 800 nm nanoantenna design is utilized for the Kagome lattice, the surface lattice resonance is not optimized in these systems, leading to the lower field enhancements observed in Fig. 2 of the main text.Nevertheless, the linear response and properties of interest here are not affected.It is important to note that such lattice resonance effects will modify the plasmonic field enhancements but have little effect on the photocurrent distributions beyond the magnitude.In particular, unlike strong near-field coupling effects that can complicate the picture of linear superpositions of individual resonators, this intermediate/far-field coupling between elements of the sub-lattices do not appreciably modify the metasurface linear responses.

Graphene layer dependence
The effect of graphene layer number is tested by comparing the THz emission from metasurfaces fabricated on monolayer, bilayer, 3-5 layer, and 6-8 layer graphene (ACS Material), with SEM micrographs shown in Fig. S4a.An overall decreasing trend of THz emission with layer number is observed (Fig. S4b,c).This can be at least partially explained by increased damping with increase layer number, which leads to less absorbed power per layer (Fig. S4d) and thus less electronic heating, although the total absorbed power across all layers remains similar (Fig. S4e).Furthermore, the effect of metal doping will decrease away from the interfacial layer, leading to less acceleration of the hot carriers excited within deeper layers.Thus, the overall decrease in the current and THz emission signal is expected (simulation result in (Fig. S4c).Nevertheless, further work is needed to understand the various relevant contributions to the layer dependence, as well as ruling out systematics or sample-to-sample variations in the fabrication process for different numbers of layers.The different dispersion around the Dirac point for monolayer versus multilayer graphene are expected to have a small effect here, given the highly environmentally p-doped graphene (chemical potential far from the Dirac point, well into the Fermi liquid regime) and high photon energies (1.55 eV here).For the present studies, therefore, monolayer graphene appears to be the optimal choice.Supplementary Note 3: Thermodynamics

Two-temperature modeling
The full spatiotemporal evolution of electron and lattice temperatures in the metasurface unit cell requires a numerical treatment, as described in the main text (see Fig. 4, Methods, and Extended Data Fig. 5).However, the temperature evolution within bare graphene and individual nanostructures can be treated via simplified two-temperature modeling, neglecting the spatial degrees of freedom, so long as the absorption cross-sections are known.This offers additional insight on the role of the nanolocalized hot spot heating compared with uniform heating due to the 2.3% absorption in bare monolayer graphene.It also offers insight into the internal metal nanostructure heating, and whether these energy sources ought to also be accounted for in the hybrid system.Thus, here we implement the following simplified two-temperature kinetic model for the electronic ( ) and phonon ( ) temperatures, evaluated separately for the graphene and gold systems: Coupling with the substrate is neglected here, and it is approximated that the absorbed power from the optical pulse ( ()) is instantaneously transferred to the thermalized electron gas, bypassing the tens-of-femtosecond electron-electron thermalization times in both systems.
For graphene, we focus only on the energy transfer to the strongly-coupled optical phonons 5 , with the electronic specific heat  ( ) described below (note 3.3), absorbed power  () =  () (Methods), energy relaxation coupling constant  = 2 × 10 WK -1 m -2 (see note 3.2), and optical phonon specific heat  taken from previous work 5 .The results are shown in Fig. S5a for 0.5 µJcm -2 incident fluence and 100 fs pulse duration.Upon comparing the peak  ≈ 360 K with the few-thousand Kelvin peak  increase in the hybrid metasurface system (Fig. S8) under the same incident fluence, the role of the concentrated and enhanced power within the plasmonic hot spot in driving the photothermoelectric dynamics is further underscored.Note that at very short times before the excited carrier distribution has thermalized,  is not strictly defined, but we maintain the instantaneous carrier-carrier thermalization approximation here.Meanwhile, in both the bare graphene and metasurface systems, the lattice (here meaning strongly-coupled nanoantenna volume, and  = 4 × 10 nm 2 the simulated nanoantenna absorption cross section for resonant 800 nm excitation),  = 2.5 × 10 W K -1 m -3 is the electron-phonon coupling constant 6 , and  = 2.4 × 10 6 J K m for gold.The results are shown in Fig. S5b for 0.5 µJ cm - 2 incident fluence.For this system, the electronic temperature increase is an even more modest ~30 K and the lattice temperature increase is negligible.Thus, the gold may serve as a thermal sink (unaccounted for in the present work) for the local hot carrier dynamics within the graphene around the tip hot spot but would not be sufficiently hot to contribute as a source in any of the dynamics.

Ultrafast dynamics measured via transient reflectivity
We perform transient reflectivity measurements to characterize the coupling constant ( ) between the electronic and lattice systems in our graphene samples.For these measurements we utilize collinear, cross-polarized 800 nm + 800 nm pump-probe microscopy to measure transient (sub-picosecond) changes in the graphene reflectivity.Pump and probe pulses from a Coherent Vitesse oscillator (80 MHz repetition rate) are focused through a Mitutoyo 20× objective onto the sample under high vacuum (~10 -6 Torr; Janis ST500 optical cryostat).The reflected probe beam is collected on a Si photodiode, with pump filtered out via a linear polarizer.A lock-in amplifier (SR830) is utilized to read out the modulated probe signal induced by the 2.5 kHz-chopped pump beam.The transient reflectivity traces (Fig. S6a) reveal two well-known general behaviors of graphene 7,8 : a fast (~400 fs) decay due to electron-optical phonon scattering and a slower (~1.1 ps) decay due to thermalization with acoustic phonons (acoustic phonon bottleneck).The slower thermalization with the acoustic phonon system can proceed via optical-acoustic phonon-phonon coupling or directly via disorder-assisted electron-acoustic phonon scattering 9 .
For suitably high signal-to-noise measurements on the monolayer graphene (with peak ~ 0.5 × 10 ), we work at high fluences up to 100 µJ cm -2 .These values also approach the locally plasmon-enhanced fluences in the metasurface system, but with stronger overall heating without the diffusive cooling pathway.The coupling parameter,  , is tuned for best agreement between the measured fluence-dependent electron-optical phonon thermalization times and twotemperature models (Fig. S6b).We estimate  ≈ 2 × 10 WK -1 m -2 based on these model assumptions and the high-temperature behavior of graphene discussed in the next section.

High-temperature properties of graphene
Under ultrafast exposure, the electronic system is transiently far from equilibrium with the lattice, reaching temperatures of several thousand Kelvin while the lattice remains near room temperature (thus avoiding issues such as electrode melting that would occur at high temperatures in studies around equilibrium).Expressions commonly utilized at low temperatures or up around room temperature will no longer be valid as  →  ( = 3600 K for  = −300 meV).We thus include a discussion on the behavior of several important quantities in the high-temperature regime.
The  -dependent electronic specific heat is critical for determining the heating of the electronic system, and can be written generically as The result is plotted in Fig. S7b, revealing  ( ) ~  dependence.This behavior is responsible for locally driving the electronic system around the nanoantennas into an apparently hydrodynamic regime, where  ≪  , influencing the spatially-varying viscosity,  ~  .The room temperature value of  was determined to be ~200 fs in previous work 14 .
The Seebeck coefficient at high temperatures can be determined from the generalized form of the Mott relation 15,16 , The measured electrical conductivity is fit to an idealized functional form 17,18 ,  () =  1 + where  is the minimum conductivity and  = 100 meV is the width of the charge neutrality region (Fig. S7c).This removes extrinsic contributions that do not influence the local conductivity at the nanoantenna tips.We approximate  (which also appears in the Wiedemann-Franz expression for  ) to be independent of  , as Umklapp scattering is suppressed in graphene.At the high temperatures relevant to the present studies, Eq.S9 yields a sublinear and even nonmonotonic dependence on  , plotted in Fig. S7d for both the bare and gold-doped graphene regions.For the bare graphene, Eq. 6 (Fig. S7a) is utilized for  ( ), while we approximate  ≈ −50 meV as a constant pinned by the gold chemical potential.The resulting Δ( ) is then responsible for the photothermoelectric acceleration of charge.

Thermodynamic energy flow
The energy cascade in graphene following ultrafast optical excitation generally involves as (i) rapid carrier-carrier thermalization within ~100 fs, (ii) thermalization with optical phonons within a few hundred femtoseconds, and (iii) full thermalization with acoustic phonons on > 1 ps timescales 7,8,10,11,19,20 .As mentioned above, thermalization with the acoustic phonon system can be accelerated via a disorder-assisted (so-called supercollision) cooling channel, directly between the electrons and acoustic phonons 9,10 .If photon energy > 2 such that interband transitions are allowed (as is the case here), separated Fermi-Dirac distributions for the electron and hole gases can be established on tens-of-femtosecond timescales before a single electronic Fermi-Dirac Fig. S8 | Electronic and optical phonon temperature evolution in graphene.Spatial temperature distributions at various times relative to the incident laser pulse, assuming instantaneous transfer of energy into a locally thermalized hot carrier distribution. .distribution is established on < 150 fs timescales 11,19 .In the modeling performed here, we approximate an instantaneous transfer of energy between the optical pulse and the fully thermalized electronic distribution.While the convolution of these ~100 fs dynamics with the ~100 fs optical pulse duration can be expected to skew the short-timescale behaviors, such effects will have little bearing on the overall thermodynamic energy and hydrodynamic momentum flows examined here.
In some previous works 5,19 , the strong coupling between the electronic and optical phonon systems has led to rapid optical phonon heating up to > 1000 K.These works, however, employ 30-300-fold higher fluences than typically utilized in the present studies (< 1 µJ cm -2 ) on bare graphene.The total absorbed power here is modest, with electronic superheating only occurring due to the highly concentrated plasmonic field distribution.As a result of this relatively low overall heat load, fast electronic thermal diffusion, and the delayed energy transfer to the optical phonon system, the peak optical phonon temperature within the unit cell is calculated to be merely 390 K (Fig. S8).Coupling to the acoustic phonon system on longer timescales is neglected here.The strongly heated electronic system (with dramatically faster  ) and near-room-temperature phonon system (with little change to  ) thus supports the view of a transient light-induced hydrodynamic phase.The resulting charge flows are discussed in the main text and examined further in the next section.
~10 nm via Gaussian convolution over the two domains.(iii) The effects of nonuniform charge density.(iv) High-temperature behaviors in the thermodynamic and hydrodynamic parameters.Despite these sources of uncertainty, the simulated flow behaviors remain robust to changes in the exact functional forms and values of the heat capacities, thermal conductivities, thermopower, etc.Therefore, the good agreement between experiment and modeling-based on the best approximations described above and using no ad hoc parameter variation/optimization-suggests that such simulations are a reasonable starting point for understanding and predicting the nanoscale flow behaviors occurring in our system.
Calculated polarization dependence is also in good agreement with experiments (and the analysis based on a simple linear response), as shown in Fig. S10.Compared with the simple linear responses and superposition of currents from the different nanoantennas described in the main text, the hydrodynamic simulations also account for coupling that occurs due to the incompressibility constraint.Even so, the symmetry of the light-matter interaction evidently prevails.

Chirality and magnetization
Most of the metasurfaces studied here have been achiral, with the exception of the perturbed Kagome lattice (Extended Data Fig. 2) with local planar chirality, and the azimuthal metasurface (Fig. 3) with global planar chirality.Transient spatially-varying magnetic fields will accompany all of the current flow profiles, though the introduction of structural chirality (i.e., with no mirror symmetry planes perpendicular to the metasurface) will introduce chiral orbital nanocurrents and net orbital magnetization.A simple example of this is shown in Fig. S11.Such systems may prove useful for controlling magnetic interactions within various symmetry-broken metasurfaces or in nearby materials, although a simple estimate limits the transient magnetic fields to the µT or low-mT ranges in graphene.This limit is imposed by the current density,  =  , where  <  (and ≪  under most experimentally accessible conditions) and  ≲ 10 cm -2 by electrostatic doping.197-200 (2005).
quality of metasurface fabrication is characterized via optical microscopy, scanning electron microscopy (SEM), and white light transmission spectroscopy.Optical bright and darkfield imaging under 50× magnification is sufficient for resolving individual nanostructures to verify successful liftoff and large-scale sample uniformity, while also resolving the monolayer graphene edge in etched devices.A more detailed view of the graphene and nanostructure morphologies (Figs.S1 and S2) is provided via SEM micrographs collected using FEI Magellan, Nova NanoSEM 450, and Nova NanoLab 600 systems.The fabricated nanoantennas closely reproduce the design profile down to the tens-of-nanometer scale, with extra-sharp ~15 nm radius of curvature tips.
using a 3D rectangular metasurface unit cell domain consisting of the air superstrate, gold nanoantennas (dielectric function given by Johnson and Christy 52 ), conductive graphene boundary layer, and fused quartz substrate 53 .Periodic (Floquet) boundary conditions are applied in the transverse x and y directions, with input/output ports along the z direction backed by perfectly matched layers.The graphene optical conductivity is approximately constant in the near-infrared frequency range 54 , with real component  = 6.1 × 10  −1 and imaginary component  = −2.1 × 10  −1 .Resonant transmission spectra are obtained from the simulated S-parameters under plane wave excitation at normal incidence, consistent with the experiments.The same nanoantenna geometry file is used for both lithography and simulations.

Fig. 4 ,
Fig.4, with further discussion of the hydrodynamic modeling provided in Supplementary Note 4.

Fig. 2 |
Fig. 2 | Polarization-dependent local responses and omni-directional control.a,b, SEM images of uniformlyoriented (a) and Kagome (b) metasurfaces with 800 nm resonances.Scale bars, 500 nm.Insets: Simulated resonant plasmonic field enhancements for different incident linear polarization angles (black double arrows), with the calculated net current direction indicated (blue single arrows).c,d, Measured x (red) and y (blue) components of the radiated THz field (solid lines with data markers) and photocurrent (dashed lines) for the uniformly-oriented (c) and Kagome (d) metasurfaces with respect to the incident linear polarization angle.Calculated linear responses (solid fills) are shown for comparison, with  signs indicating lobe polarity.For clarity, a small residual y component is not shown in (c).e, The Kagome metasurface exhibits nearly constant photocurrent magnitude (purple) and continuously-rotatable direction (gray), consistent with analytic predictions.

Fig. 3 |
Fig. 3 | Global vectorial currents and THz vector beams.a, SEM image of the central region of a radial vector metasurface, with white arrows illustrating photocurrent in the radial direction upon circularly-polarized excitation.Scale bar, 1 µm.b, Spatial mapping of radial THz vector field.c, Measured (top) versus ideal (bottom) Hermite-Gaussian modes for the x and y field components of the radial THz vector beam.d, Lineout from the radial Ex image, showing the transient THz field as a function of x position, along with example THz time-domain waveforms illustrating the polarity reversal across the beam.Inset: The corresponding THz field amplitude spectrum in the frequency-domain on a logarithmic scale.e-g, same as (a-c) but for the azimuthal vector metasurface and THz beam.

Fig. 4 |
Fig. 4 | Local photothermoelectric driving force and resulting nanoscale charge flow.a, Measured resistivity of a metasurface device as a function of gate voltage.Insets: Chemical potentials of the gold-pinned and bare graphene regions for three gate voltages as indicated in (b).b, Measured gate-dependent current (dashed line) compared with calculated S0 (solid fill).Inset: Illustration of the gated device configuration, with tip-localized heating and corresponding net carrier motion.c, Spatial distribution of S0. d, Snapshot of the calculated Te distribution at the time of the incident pulse peak (0 fs). e, Photothermoelectric x acceleration field (along the principal nanoantenna axis) at  = 0 V and 0 fs.f, Hydrodynamic velocity field for the hole fluid at  = 0 V and 0 fs.

Fig. S2 |
Fig. S2 | Effect of bare graphene removal on metasurface resonance and THz emission.a, Scanning electron micrograph before and, b, after removal of exposed graphene between nanoantennas.Insets: Individual nanoantennas with 15 nm radius circle illustrating tip radius.c, White light transmission spectrum before and, d, after graphene removal.e, THz time trace before and, f, after graphene removal.

Fig. S3 |
Fig. S3 | Effect of metasurface density.a, Electron micrographs of metasurfaces with 250 nm (bottom) to 750 nm (top) square lattice constant.b, Experimental (top) and simulated (bottom) relative transmittance spectra, correcting using the bare graphene/substrate as a reference.c, Simulated tip hot spot field intensity.d, Experimental (solid-dot line) THz field amplitude as a function of lattice constant, measured at 800 nm incident wavelength, compared with the simulated value from the density-scaled tip intensity (solid fill).

Fig. S4 |
Fig. S4 | Effect of graphene thickness.a, Electron micrographs of metasurfaces fabricated on monolayer, bilayer, 3-5 (~4) layer, and 6-8 (~7) layer graphene.b, THz time domain traces for different layer numbers.c, Summary of the measured THz field amplitude (dash-dot line) compared with the peak simulated tip field intensity (solid fill) from panel d.d, Simulated field intensities at the nanoantenna tip, a few nanometers above the graphene interface.With graphene defined as a boundary layer, the layer number () is accounted for via modification of  →  .e, Same as panel d, but scaled by .
optical phonon) temperature rise is minimal at only a few tens of Kelvin.See note 3.4 for calculations of the time-evolving  within the metasurface unit cell.Now considering only the internal gold nanostructure thermal evolution, the relevant parameters are  =  2  B 2   2 F   (from Sommerfeld theory;  = 5.53 eV and  = 5.9 × 10 cm -3 ),  () =  () (with  = 377 Ω the impedance of free space,  = 3 × 10 nm 3 the

Fig. S5 |
Fig. S5 | Two-temperature modeling for bare graphene and gold nanoantennas.a, Electron and optical phonon temperature evolution for bare graphene (no metasurface) under 100 fs pulse excitation centered around 0 delay time, neglecting acoustic phonon coupling.b, Electron and lattice (optical and acoustic phonon) temperature evolution for gold nanoantenna. .

Fig. S7 |
Fig. S7 | High-  calculations for graphene.a, Chemical potential (top) and specific heat (bottom).b, Scattering time relative to  value compared with standard  Fermi liquid theory behavior.c, Measured conductivity for graphene device and idealized fit excluding extrinsic effects, utilized for calculations of the Seebeck coefficient.d, Seebeck coefficient calculated for the bare and metal-doped graphene regions based on Eq.S9, compared with the linear- behavior that prevails at low temperatures. .

Fig. S10 |
Fig. S10 | Comparison between analytic and hydrodynamic polarization dependence.a, Flow profile for the unit cell f the square lattice.b, Polarization dependence of the photocurrent (x component; y component is negligibly small) calculated via hydrodynamic model.c, Polarization dependence of the photocurrent calculated analytically for a linear response (x component; y component is identically zero by symmetry).d-f, Same as panels a-c but for the Kagome lattice unit cell (now including y photocurrent component). .

Fig. S11 |
Fig. S11 | Local vortical flow and net magnetization in a chiral unit cell.An example of a chiral nanostructure arrangement with corresponding counter-clockwise hydrodynamic flow, serving as a local magnetic dipole. .