Generation and control of localized terahertz fields in photoemitted electron plasmas

Abstract

Dense micron-sized electron plasmas, such as those generated upon irradiation of nanostructured metallic surfaces by intense femtosecond laser pulses, constitute a rich playground to study light–matter interactions, many-body phenomena, and out-of-equilibrium charge dynamics. Besides their fundamental interest, laser-induced plasmas hold great potential for the generation of localized terahertz radiation pulses. However, the underlying mechanisms ruling the formation and evolution of such plasmas are not yet well understood. Here, we develop a comprehensive microscopic theory to predictably describe the spatiotemporal dynamics of laser-pulse-induced plasmas. Through detailed analysis of electron emission, metal screening, and plasma cloud interactions, we investigate the spatial, temporal, and spectral characteristics of the so-generated terahertz fields, which can be extensively controlled through the metal morphology and the illumination conditions. We further describe the interaction with femtosecond electron beams to explain recent ultrafast electron microscopy experiments, whereby the position and temporal dependence of the observed electron acceleration permits assessing the associated terahertz field. Besides its potential application to the design of low-frequency light sources, our work contributes fundamental insight into the generation and dynamics of micron-scale electron plasmas and their interaction with ultrafast electron pulses.

---

1 Introduction

Terahertz (THz) radiation—a part of the electromagnetic spectrum sandwiched between microwaves and infrared light—has recently attracted significant attention because of its potential application in areas such as spectroscopy,^1–3 sensing,⁴ imaging,^5–7 and communication technologies.⁸ In this context, nanophotonics constitutes a suitable arena to test and capitalize on some unique properties of THz radiation such as the ability to penetrate through optically opaque materials⁹ and a high sensitivity to chemical composition,^5,10 both of which can be manipulated through material nanostructures. However, the efficient generation of THz light remains a challenge, and even more so when aiming for nanoscale sources.

The production of THz fields typically relies on nonlinear optical phenomena such as wave mixing,^11,12 optical rectification,^13–15 and frequency conversion.^16,17 These methods involve simple setups fed by high-frequency optical sources such as lasers, but they generally have low efficiencies and are limited by the availability of suitable nonlinear crystals.

Electron plasmas have emerged as an appealing alternative for the generation of intense THz fields.^18–20 Such plasmas can be obtained by extracting electrons from metal surfaces upon intense laser-pulse irradiation via multiphoton photoemission and thermionic emission.²¹ If the intensity of the ionizing laser is large enough, a short-lived electron-plasma plume of a few picoseconds in duration can be formed, characterized by a high density of emitted electrons that are eventually reabsorbed by the surface or escaping away from the metal. The associated charge motion gives rise to intense, transient localized THz fields, but the precise underlying mechanisms are not yet fully understood.²²

Dense electron-plasma plumes undergo a complex spatiotemporal dynamics ruled by the collective interaction among many electrons in the presence of screening by the metallic structure, thereby posing an important challenge for a comprehensive theoretical description. Nevertheless, besides their potential for application in THz technologies, the study of this phenomenon bears interest as a source of fundamental insight into the ultrafast dynamics of complex nanoscale systems, as revealed by recent experimental results obtained by employing ultrafast electron microscopes,^23,24 whereby electron beam (e-beam) pulses are made to interact with the plasma at controlled delay times relative to the laser pulses. In fact, high-energy electrons are ideal probes for ultrafast and localized phenomena such as charged plasmas^25–29 due to their ultraconfined nature, enabling a spatial/temporal resolution down to sub-nanometer/femtosecond scale, combined with a high sensitivity to electromagnetic interactions.³⁰

In this work, we study the formation and evolution of electron plasma produced upon irradiation of metal nanostructures by intense laser pulses, as depicted in Fig. 1(a and b), through a parameter-free theoretical formalism that incorporates a quantitative description of electron emission, metal screening, and cloud dynamics, including electron reabsorption and the generation of localized THz fields. The process also involves substantial heating of conduction electrons in the metal, triggering ultrafast thermal dynamics that needs to be accounted for to formulate accurate predictions on the behavior of the plasma. The present model has been successfully used to explain recent experimental results on the ultrafast nanoscale spatiotemporal dynamics of electron plasmas probed by femtosecond electron pulses.²³ We provide a comprehensive description of the formalism and apply it to study the so-generated transient THz fields, whose duration, spatial distribution, and spectral composition are strongly dependent on the metal morphology and illumination conditions. The latter provide suitable knobs to control the THz field characteristics for potential applications.

Fig. 1 Electron plasma formation by laser-pulse irradiation. (a and b) Schematic representation of (a) a metallic sample irradiated by a high-fluence laser pump pulse at time t = t₀ and (b) the ensuing formation of an electron plasma evolving at t > t₀. A synchronized e-beam probe pulse is used to study the spatiotemporal dynamics of the plasma by scanning the e-beam position and delay time relative to the laser pulse. (c) Scheme of the specific system under consideration, consisting of a metallic wedge that is translationally invariant along y and characterized by its length a, angle α, and tip radius r, under illumination by an external field polarized along y. The surface contour is parameterized by s, ranging counterclockwise from 0 to 1. (d) Dynamics of the electron temperature T_e at the surface of a copper wedge as a function of position s and time. The vertical dashed lines mark the positions of the blue dots in (c) along the surface of the wedge, with the one at s = 0.5 corresponding to the tip apex (point C). The wedge parameters in this specific calculation are a = 10 μm, r = 1 μm, and α = 20°, while the laser pulse has a fluence F_pump = 200 mJ cm⁻² and its maximum is arriving at time t₀ = 0.25 ps.

2 Theoretical description of laser-induced electron plasmas

To demonstrate the ability of electron plumes to generate and control localized THz fields, as well as to reveal fundamental insights into their origin, we introduce a theoretical framework that describes the generation and evolution of laser-pulse-generated plasmas. The theory here presented can be applied to any translationally invariant morphology, while a generalization to fully arbitrary shapes is straightforward. For concreteness, we provide numerical simulations for a system composed of a translationally invariant wedge (along y) with the cross-sectional geometry shown in Fig. 1(c), characterized by a length a along the x direction, an angle α, and a tip radius r. The wedge surface is described by a parameter s that varies counterclockwise from 0 to 1, starting from the surface point opposite to the tip apex, as indicated in Fig. 1(c). We study this structure under uniform illumination from the side by a high-fluence laser pulse with an external field polarized along y (i.e., the direction of translational symmetry), as indicated by the red arrow in Fig. 1(c). In the simulations presented below, we consider a laser central wavelength λ = 800 nm, a laser pulse duration Δt = 60 fs, and a fluence F_pump = 200 mJ cm⁻².

2.1 Sample temperature dynamics

The absorption of laser energy by the metal structure raises the temperature of its conduction electrons, whose thermal distribution depends on the local field enhancement, and this in turn on the surface geometry and the illumination conditions. We model the evolution of the electron temperature T_e(r, t) as a function of time t and position r on the metal surface using the two-temperature model described in Sec. 5.1 in Methods.

An example of the surface electron temperature dynamics is presented in Fig. 1(d) for a copper wedge with geometrical parameters a = 10 μm, r = 1 μm, and α = 20°. The temperature is first experiencing a fast increase within ∼100 fs, with a hotspot concentrated near the tip apex (s = 0.5), where it reaches T_e ∼ 8000 K. This is followed by a slow cooldown lasting for ∼5 ps, after which the conduction electrons return to ambient temperature. The temporal profile of the temperature along the A–E cuts in Fig. 1(d) are presented in Fig. S1 of the ESI,† together with analogous results for a similar wedge with smaller tip radius r = 0.1 μm.

2.2 Plasma emission

During the time over which the metallic sample is being illuminated, as well as the subsequent period in which conduction electrons near the sample surface remain hot, substantial electron emission can occur primarily due to thermionic emission and n-photon photoemission (with n = 3 in the present study, see Sec. 5.2 in Methods). The former occurs while the electron temperature remains elevated (roughly ≳1 ps), allowing some high-energy electrons to overcome the work function Φ of the metal. In contrast, 3-photon photoemission (3PPE) occurs exclusively during the pumping period (<100 fs) via the absorption of three photons by one electron on the metal surface. These mechanisms, which are discussed in Sec. 5.2 in Methods, yield a total emission rate given bywhere the electron emission rates corresponding to each mechanism, and , are given by eqn (8) and (11), respectively. Notably, these rates are strongly dependent on the local electron temperature, which is in turn a function of position and time.

Using eqn (1) together with the data in Fig. 1(d), we can retrieve the density of electrons emitted from a sample surface position r_s at time t, which is here normalized per unit length along the direction y of translational invariance in the sample. An example of total emission (i.e., integrated along the transverse surface profile) from each of the two mechanisms considered is shown in Fig. 2(a) for the copper wedge introduced in Fig. 1(d). We note that the 3PPE process is dominant during the pumping period (centered around t₀ = 0.25 ps). In contrast, thermionic emission becomes dominant at a later time [when the electron temperature reaches a maximum, see Fig. 1(d)], and remains dominant during a comparatively longer time.

Fig. 2 Plasma-induced THz electric field. (a) Electron emission rate due to thermionic emission (red), 3PPE (blue), and both channels combined (black) as a function of time. (b) Temporal evolution of the total number of emitted electrons (purple) and the number of electrons remaining in free space (i.e., those that are not yet reabsorbed; green) under the same conditions as in (a). (c–j) Electron plasma density ρ_e and image-charge surface density eσ_s (c, e, g and i), along with the corresponding electric field maps (d, f, h and j), at four different times t = t₀ + τ referred to the time t₀ of maximum incident pump laser intensity [indicated in panels (a and b) by an orange dashed line]. Color maps in (d, f, h and j) represent the electric-field amplitude as a function of position, supplemented by field lines (black) tangent to . All panels are calculated for a copper wedge with parameters a = 10 μm, r = 1 μm, and α = 20°, as well as a pump of fluence F_pump = 200 mJ cm⁻² arriving at t₀ = 0.25 ps.

2.3 Plasma dynamics

When electrons are ejected in large numbers from the surface, they form an electron plasma characterized by a volume density ρ_e(r, t) that evolves rapidly due to the electromagnetic interaction among the different electrons and the effect of screening by the sample. The latter is driven by the accumulation of positive image charges along the sample surface, with a hole surface density σ_s(r_s, t) depending on the distribution of electrons outside the metal. After an initial fast expansion of the plasma, the attractive interaction between the negatively charged electron cloud and the positively charged surface produces a deceleration in the outgoing motion of the ejected electrons and, eventually, partial reabsorption of plasma electrons. Incidentally, metal screening has a characteristic time³¹ ≲1 fs (over which the plasma changes negligibly), and the penetration of the surface screening charge is of the order of the Thomas–Fermi screening length ≲1 nm; consequently, we model screening in the perfect-metal approximation (see Sec. 5.4 in Methods).

To realistically describe the spatiotemporal plasma dynamics, we develop a theoretical approach capable of simulating the emission and subsequent evolution of the dense electron plasma cloud, as we describe in Sec. 5.3 in Methods. Using this procedure, the simulated spatiotemporal evolution of the plasma density ρ_e(r, t) and the associated surface hole density σ_s(r_s, t) are both represented in Fig. 2(c, e, g and i) as snapshots for selected time delays τ relative to the time t₀ = 0.25 ps at which the incident pump laser intensity is maximum.

Right after pumping (τ ≈ 0 ps), the emission is dominated by the 3PPE channel, and therefore, the plasma accumulates heavily near the sample surface. This occurs because electrons emitted via this process carry a comparatively small kinetic energy K ≈ 3ħω − Φ ∼ 10's meV (see Methods), and consequently, they remain relatively close to the surface. As a result, the vast majority of them are rapidly reabsorbed due to image attraction by the metal surface. This behavior is revealed by the sharp peak observed in the green curve of Fig. 2(b), representing the number of electrons that remain in free space (i.e., those that have not been reabsorbed) as a function of time, which we compare to the total number of emitted electrons (purple curve). After a few 100's fs, most of the initially emitted electrons have already been reabsorbed.

Conversely, after this short period, thermionic emission takes over as the dominant emission mechanism, and on average, the so-emitted electrons have a much higher kinetic energy K ∼ 1 eV, producing a noticeable expansion of the plasma cloud up to a few microns from the surface, as revealed by the density maps in Fig. 2(c, e, g and i). This evolution is accompanied by fast electron scattering along the wedge side [cf.Fig. 2(e) and (g)] due to the transverse asymmetry of the charge distribution, which is strongly concentrated near the tip apex. For this reason, the newly emitted electrons feel a weakened charge barrier when compared to smoother geometries,²³ thus resulting in a faster expansion of the plume. Nevertheless, the intense attractive force of the metal image charges, represented by green lines in Fig. 2(c, e, g and i), progressively produces a deceleration and subsequent reversion of the cloud expansion, so that most electrons are eventually reabsorbed. After 10 ps, only ∼10% of the emitted electrons remain in free space, most of which can escape. This amounts to a relatively high portion of escaping electrons, a fact that we explain by the aforementioned small charge barrier effect of this particular geometry.

In the calculations presented in this work, we assume the metal structure to be electrically isolated, such that the system maintains charge neutrality (i.e., the number of electrons in the plasma is fully compensated by the number of holes distributed on the metal surface). For grounded samples, additional charges should partially refill the holes, therefore reducing electron reabsorption and affecting the plasma dynamics. We expect this effect to be small for smooth surfaces characterized by large curvature radii, which should produce a screening that is locally approaching the limit of a planar metal surface (i.e., simultaneously meeting the conditions of charge neutrality and a vanishing surface potential).

3 Results and discussion

The charged electron cloud generates an intense electric field in its interior and vicinity. We show in Fig. 2(d, f, h and j) the electric field generated by the electron and surface-charge distributions plotted in Fig. 2(c, e, g and i), respectively. These maps display a strong concentration of the field amplitude in a region extending up to ∼5–10 μm around the tip surface. The field reaches a maximum value right after pumping (i.e., when a large electron pileup is found close to the wedge surface) and slowly dies out as the electrons scatter mainly due to interactions with other electrons in the plasma. Similar dynamics can be observed in a wedge with smaller tip radius (see Fig. S2 in the ESI†), in which the electrons are emitted from a smaller surface area, thus producing a more concentrated spatial field distribution.

3.1 Frequency decomposition of the generated fields

We perform a spectral analysis of the field produced by the laser-induced plasma by Fourier transforming the time-dependent field [shown in Fig. 2(d, f, h and j) for selected instants]. Representative examples of the resulting frequency-domain field are shown in Fig. 3(a) and (b) for wedges with a tip radius of 1 and 0.1 μm, respectively, and for the specific frequency f = ω/2π = 7 THz. In Fig. 3(c) and (d), we represent the full spectral decomposition of the field at selected positions lying at a distance of 0.5 μm from the metal surface, as indicated by color-coordinated dots in Fig. 3(a) and (b), respectively. Upon comparison of the field maps in Fig. 3(a) and (b), we conclude that a smaller tip radius produces a stronger spatial concentration of the near field in the vicinity of the sample, which we explain as the result of electron emission arising from a smaller surface area. Nevertheless, in both cases the field intensity decreases rapidly with distance to the tip, as illustrated in Fig. 3(e) and (f) by plotting the spectral decomposition of the THz field at several positions separated by a distance of 10 μm from the sample surface.

Fig. 3 Plasma THz electric field in the frequency domain. (a and b) Spatial map of the f = 7.0 THz frequency component of the Fourier-transformed electric field near (a) the same copper wedge as in Fig. 2, with a tip radius r = 1 μm, and (b) a similar wedge but with r = 100 nm and all other parameters unchanged. The color map represents and we superimpose field lines (black) tangent to the vector . (c and d) Spectral decomposition of the electric field around the wedges in (a and b), respectively, at the specific spatial positions marked by color-coordinated dots in the insets [and also in (a and b)], situated at a distance of 500 nm from the metal surface. The vertical dashed lines mark the f = 7.0 THz component. (e and f) Same as (c and d), but for positions marked by the square dots in the insets, placed at a distance of 10 μm from the metal surface. We set the pump fluence to F_pump = 200 mJ cm⁻² in all panels. The arbitrary units used for the electric field amplitude are common to all plots.

In both wedges, the electric field is dominated by low THz frequency components with a mean frequency 〈f〉 ∼ 7 THz. The THz nature of these fields is commensurate with the time scale over which the electron cloud evolves: at a fixed spatial position, there is an initial fast variation in charge density (over the first ∼2 ps), responsible for the high-frequency components, followed by a longer period (∼10 ps) in which most electrons have been reabsorbed and the density displays a slow evolution [see Fig. 2(b)], giving rise to the sharp increase in low-frequency contributions observed in Fig. 3(c–f). Remarkably, fast electrons that are able to escape the sample leave a distinctive signature along their trajectory, as observed in Fig. 3(b) [see also Fig. S2(b, d, f and h)†]: due to the rapid variation of the time-dependent field at those positions, we observe high-frequency field components in the purple curve in Fig. 3(d).

The similarity between the average field frequency observed for the two wedges can be attributed to the comparable plasma charge density in the two scenarios (cf.Fig. 2(c, e, g, i) and S2(a, c, f, h) in the ESI†), which crucially influences the motion of plasma constituents and, therefore, the frequency composition of the generated field. This is, in turn, prompted by the tradeoff arising from a larger temperature enhancement over a smaller surface area for the smaller wedge tip. Nevertheless, the detailed spectral profile for each wedge depends on both position and surface geometry, suggesting that the latter, together with the illumination conditions, offers potentially useful means to control the generated THz field.

3.2 Electron probing of THz fields

The above analysis suggests that fast electron pulses traversing the plasma with a controlled delay time relative to the laser pump can serve as excellent probes of the temporal, spectral, and spatial characteristics of the generated THz field. To illustrate this idea, we extend our theoretical formalism to incorporate the interaction with a fast probing electron, producing excellent results in comparison with experiments, as shown in recent publications for a different geometry.^23,24 For the wedge structure investigated in this work, we consider an electron passing at a distance b from the surface, with a velocity vector v_e making an angle θ relative to the positive x direction, as indicated by the green arrows in Fig. 4(a) for different values of θ in the α/2 ≤ θ ≤ π − α/2 range, such that only aloof electron trajectories are considered. Although the angle θ spans the range [0, 2π), we note that the spatiotemporal dynamics of the plasma is symmetric with respect to the y = 0 plane [see Fig. 1(c)], and hence, the trajectories described by θ and 2π − θ produce the same outcome. We assume the electron wavepacket to be well focused in the transverse e-beam direction and spanning a full-with-half-maximum (FWHM) temporal duration Δt_e. Furthermore, we neglect any changes produced by the interaction on the electron velocity v_e (nonrecoil approximation, see Sec. 5.6).

Fig. 4 Probing THz fields with e-beam pulses. (a) Scheme showing different electron trajectories (green arrows) passing at a distance b from the surface of the copper wedge considered in Fig. 2 and forming an angle θ with the x axis in the α/2 ≤ θ ≤ π − α/2 range, where α = 20°. We introduce a delay τ between the swift electron and laser pulse. (b) Frequency decomposition of the electric field along the electron trajectory as a function of frequency f and angle θ for three different values of τ (see labels). Green curves represent the peak frequency as a function of θ. (c) Net post-interaction electron energy change ΔE as a function of the trajectory angle θ. The vertical dashed lines indicate the θ limits α/2 and π − α/2. We take a pump fluence F_pump = 200 mJ cm⁻², an impact parameter b = 1 μm, an electron velocity v_e ≈ 0.7c, and an electron pulse duration Δt_e = 600 fs.

In Fig. 4(b), we show the frequency-domain electric field amplitude as a function of frequency f and electron trajectory angle θ for three selected values of the electron delay τ relative to the pump laser pulse. We set the electron velocity to v_e ≈ 0.7 c and the impact parameter to b = 1 μm. Analogous results for a wedge with smaller tip radius are presented in Fig. S3 of the ESI.† The maps in both figures confirm that the spectral landscape of the field experienced by the electron depends strongly on both the trajectory angle θ and the delay τ. In particular, the dominant frequency (i.e., that for which the field intensity is maximum for a given angle θ, represented by the green curves in each plot) can be varied within the ∼3–10 THz range by adjusting these parameters. Likewise, the electron velocity v_e and the impact parameter b are additional trajectory parameters that can be varied to explore the frequency landscape, which is also strongly dependent on sample geometry and illumination conditions. As previously reported^23,24 and discussed in Sec. 5.6 in Methods, these electric fields produce a net energy variation ΔE of the probe electron, which we plot in Fig. 4(c) as a function of the trajectory angle θ for the same values of the delay τ as in panel (b) (see also Fig. S3 in the ESI†). The dependence of ΔE on these parameters is complex because it is mediated by the plasma dynamics during the interaction time. For example, a net gain or loss is observed depending on the electron trajectory. State-of-the-art electron spectrometers can currently resolve energy changes down to <10 meV,³² rendering this approach highly sensitive to minuscule details in the plasma dynamics along the probe electron trajectory.

4 Conclusions

The electron clouds that arise upon irradiation of metallic surfaces with intense laser pulses act as sources of intense THz electromagnetic fields localized over micrometer-sized regions. We have shown that both the spatial extension and the spectral composition of these fields are extremely sensitive to the surface geometry and the characteristics of the pumping. The surface morphology and the illumination conditions are thus elements that can be engineered to control the generated THz fields. These sources could find application in sensing molecular vibrations at similarly low frequencies, with a spatial resolution well below the field wavelength gained by appropriately shaping the metal surface. The localized nature of the generated field is appealing to eliminate spurious signals coming from far regions away from the sample of interest. In addition, the produced THz radiation could be detected without any contamination from the incident laser field, which lies within a completely different spectral region.

The studied process involves the presence of metal holes due to electron ejection. Such holes are redistributed along the surface, acting as screening charges that strongly affect the plasma dynamics. In this work, we have assumed electrically isolated structures in which charge neutrality leads to a number of holes exactly compensating for the number of plasma electrons. An interesting scenario could be encountered when considering subsequent laser pulses (i.e., impinging on a previously ionized sample), for which we would expect different plasma dynamics under the influence of the net electrostatic potential landscape produced by previous pulses, and eventually, a stationary regime should be reached in which no electrons escape from the structure. A different behavior is also anticipated for grounded samples, in which additional electrons can refill the holes as electrons are ejected away from the surface region. Deviations in the performances of grounded and isolated structures are expected to be stronger in the presence of sharp surface profiles like those considered in this work. We envision the use of an externally controlled degree of insulation (e.g., through a variable resistor) to switch between these two scenarios, thus providing additional means of active control over the generated THz radiation. Nevertheless, when large structures are considered, regions away from the laser-irradiated area act as an effective ground, and therefore, this discussion should be more relevant in smaller samples.

As a way to characterize plasma dynamics in this context, we have shown that a passing electron beam pulse with a controlled trajectory can selectively probe specific frequency components, thus offering a unique way to map the spatiotemporal evolution of laser-pulse-induced microplasmas. Probing the ultrafast out-of-equilibrium dynamics of charged-carrier clouds is a challenging problem, whose solution bears interest from both fundamental and applied perspectives. By using the present theory, we have explained recent experiments of spatiotemporal plasma mapping in the context of ultrafast electron microscopy,²³ which have served as a testbed to elucidate the ingredients that play a relevant role in such a complex process, involving different scales of time (from sub-femtosecond metal screening to picosecond plasma evolution), length (from a few nanometers in electron emission and surface charge dynamics to microns in plasma plume dynamics), and energy (from a few electronvolts needed to eject electrons from the metal surface to 100's keV probe electron energies). The effects produced on the probing electron suggest the possibility of designing a disruptive type of micron-sized electron optics component, whereby the wave function associated with free electrons is manipulated by subjecting them to a sizeable and widely controllable interaction with plasma plumes.

5 Methods

5.1 Two-temperature model

We describe the temperature dynamics in a metallic sample irradiated by ultrafast laser pulses through the two-temperature model (TTM), in which the electron and lattice temperatures within the material (T_e and , respectively) are taken as independent variables. For the pulse fluences here considered, the variation in lattice temperature can be neglected (, where T₀ is the ambient temperature) and the electron temperature thus obeys the differential equation³³where c_e is the electron heat capacity, κ_e is the electron thermal conductivity, p^abs is the power density absorbed from the laser, and G describes electron-phonon coupling.

We calculate the electronic heat capacity of the metal c_e = ∂Q_e(T_e)/∂T_e from the derivative of the temperature-dependent electronic heat density,

where ρ(E) is the density of states (DOS), Θ is the step function, f_μ,Te(E) = {exp[(E − μ)/(k_BT_e)] + 1}⁻¹ is the Fermi–Dirac distribution, and μ is the chemical potential. The latter depends on temperature as determined by the conditionexpressing the conservation of the number of electrons in the system.

In eqn (2), the absorbed power density at position r and time t is given by

where ω is the pump frequency, ε(ω) is the metal permittivity, and is the optical electric field (including scattering by the metal structure). We write the latter as , where is the pump field amplitude expressed in terms of the fluence F_pump, t₀ marks the time of maximum pulse intensity, we define with Δt standing for the FWHM of the intensity, and is the local field enhancement that we calculate using the boundary-element method³⁴ (BEM). The latter is presented in Fig. S1(a) of the ESI† for the wedge geometries considered throughout this work.

By assuming that the material surface has a smooth profile characterized by a local curvature radius that is large compared with the light wavelength λ = 2πc/ω, we solve the TTM locally as a 1D model in which any lateral heat diffusion (i.e., along directions parallel to the surface) is neglected and only diffusion along the local direction perpendicular to the surface is considered. This approximation largely simplifies the problem, so that the evolution of the temperature T_e(ζ, t) (with ζ standing for the distance from the metal surface towards its interior) can be readily determined from eqn (2) using a standard partial differential equation solver.

In this work, we apply this procedure to copper structures, using tabulated data for the DOS of this material³⁵ and a Fermi energy E_F ≈ 9.5 eV³⁵ corresponding to the chemical potential at T_e = 0. In addition, we adopt experimental values for the thermal conductivity³⁶κ_e and the electron-phonon coupling coefficient³⁷G. In the present simulations, we set the light wavelength to λ ≈ 800 nm, for which the copper permittivity is ε ≈ − 25.07 + 2.54i, and consider pulses with a duration Δt = 60 fs. Using these parameters to feed the BEM and the TTM, we find the temperature dynamics illustrated in Fig. 1(d) and S1 (see ESI).†

5.2 Photothermal electron emission mechanisms

Under the illumination conditions considered in this work, we assume that electron emission from the metal surface is dominated by two different mechanisms: (1) direct thermionic emission and (2) n-photon photoemission. The first of these mechanisms takes place while the metal surface remains hot (for ≳1 ps), such that the elevated electron temperature promotes electrons from lower-to higher-energy states according to the Fermi–Dirac distribution, thus dramatically increasing the electron population for energies above the potential barrier and resulting in electron escape, as depicted in Fig. 5(a). In contrast, n-photon photoemission occurs only during the pumping time (<100 fs) and is driven by the absorption of n photons by one electron, providing it with enough energy to overcome the potential barrier, as depicted in Fig. 5(b). We describe each of these two mechanisms below, as well as alternative emission processes.

Fig. 5 Dominant electron emission mechanisms. Schematic representation of (a) thermionic emission and (b) n-photon photoemission. In both panels, the dashed orange curve represents the density of states, the filled solid orange curve shows the density of occupied states ρ_occ for electrons at room temperature T_e ≈ T₀, and z stands from the vacuum distance away from the metal surface. The barrier height V₀, the Fermi level E_F, and the work function Φ are the same in (a and b). The energy origin E = 0 is chosen at the barrier top. Electron emission (red arrows) is assisted by an elevated electron temperature T_e in (a), where the filled purple curve shows ρ_occ for T_e ≫ T₀, and by absorption of n photons of energy ħω in (b) (n = 3 in the scheme).

5.2.1 Thermionic emission. We study this process for an infinite planar surface normal to the z direction, under the approximation of a smooth surface profile (see above). Assuming that conduction electrons inside the metal are confined to a potential well of width L along z, the probability per unit area of ejecting an electron across the barrier can be written aswhere A is the surface area, the factor of 2 accounts for spin degeneracy, is the electron wave vector, E = ħ²k²/2m* − V₀ is the electron energy, E_z = E cos²θ is the electron energy along z, θ is the emission angle with respect to the z direction,is the transmittance across the surface energy barrier of height V₀, τ = 2L/v is the average time interval separating two consecutive electron collisions against the potential barrier, is z component of the electron velocity, and m* is the effective electron mass. Transforming the sums in eqn (6) into integrals through the substitutions and , we obtainwhererepresents the probability of thermionic emission of an electron of energy E along an angle θ when the electron surface temperature is T_e. From here, it follows that electrons are primarily emitted around the surface normal (average emission angle θ_av = 0) with average energycorresponding to an emission velocity .

5.2.2 n-Photon photoemission. Under the illumination conditions considered in this work (photon energy ħω ≈ 1.55 eV, copper work function³⁸Φ ≈ 4.65 eV), we have 4 > Φ/ħω ≳ 3, so that this emission channel is dominated by n = 3 processes. The corresponding photoemission rate is calculated using the well-known Fowler–Dubridge model, according to which the emission probability is given by³⁹where is the Richardson constant, I_abs is the absorbed power density, Φ_Te = Φ + E_F − μ is the temperature-corrected work function, is the Fowler function,³⁹ and a₃ ∼ 5 × 10⁻³⁶ cm⁶ A⁻³ (ref. 23) represents the likelihood of the emission. Finally, the average energy distribution of the photoemitted electrons can be calculated aswhile the average emission angle is again θ_av = 0 due to symmetry.

5.2.3 Alternative emission mechanisms. Under strong field illumination, conduction electrons could escape from the metal via tunneling into the vacuum due to the periodic lowering of the potential barrier by the incident laser electric field. According to the Keldysh criterion, this mechanism is negligible compared to n-photon photoemission if γ ≪ 1, whereis the dimensionless Keldysh parameter,^21,40 defined in terms of the metal work function Φ and the ponderomotive energy , which is in turn expressed in terms of the light frequency ω and the electric field amplitude . For the copper samples considered in this work, under illumination by a 800 nm laser with a peak of intensity of ∼300–1500 GW cm⁻², we have γ ∼ 10–20 ≫ 1, and consequently, we neglect tunneling emission. Another possible emission mechanism is the escape of nonthermal high-energy electrons during a short period right after pumping when the system is strongly out of equilibrium. However, we expect this contribution to only amount to a small correction in the total emission, and thus, we neglect it as well.

5.3 Plasma dynamics

We now describe our numerical implementation to simulate the emission and spatial evolution of the electron plasma, starting with a discretization of the surface through a set of positions s_j, and also the time intervals t_i at which electrons have been emitted, where i and j are discretization indices. Electrons within each set of (i, j) indices are evolved independently, taking into consideration the interaction with both surface charges and other sets. We are interested in the evolution of the corresponding densities of emitted electrons ρ_ij(t), and further represent the dynamics of each (i, j) set with a time-dependent average velocity v_ij(t). The initial population of every (i, j) set is determined by the emission probabilities in eqn (1), which are described in Sec. 5.2. To alleviate the computational burden, we consider all electrons to be ejected normally to the local surface with a velocity determined by the corresponding average energy E_av (see above), which is, in turn, dependent on (i, j) through the local field amplitude and electron temperature.

To compute the dynamical evolution of each (i, j) set, we need to calculate the force acting on the electrons at any given time t ≥ t_i:

where refers to the electron–electron (ee) interaction with the rest of the previously emitted (i′, j′) sets, while is the contribution of surface charges [electron–hole (eh) interaction; see Sec. 5.4] summed over surface positions . Neglecting magnetic interactions due to the small drift velocity of the emitted electrons, the ee force component is given bywhere r_ij = (x_ij, 0, z_ij) and we introduce the retarded time of interaction t_r = t − |r_ij − r_i′j′|/c. We find that this retardation correction affects the results because of the large extension of the plume, which is not negligible compared with the wavelength associated with the generated THz field. It should be noted that the y and y′ integrals are needed because ρ_e (and also the surface charge density, see below) is defined per unit length along that direction and we assume translational invariance in both the geometry and the pump. To connect with experiments, in which the pump laser beam has a finite lateral extension, we have introduced a parameter D_y accounting for an effective length along y, within which we approximate the density of plasma electrons (and also surface charges, see below) to be constant. We set D_y = 25 μm in the present calculations. Analogously, the eh contribution readswhere is the density of holes per unit length along y within the surface interval represented by the point of coordinates , and .

The corresponding acceleration that this force exerts on electrons in the (i, j) set is given by a_ij(t) = F_ij(t)/M_ij(t), where we assimilate M_ij(t) = m_eD_yρ_ij(t) to the total mass of a uniform charged line with extension D_y along y and mass density m_eρ_ij(t). The velocity and position are then updated at each time step δt according to Newton's equations as v_ij(t + δt) = v_ij(t) + a_ij(t)δt and r_ij(t + δt) = r_ij(t) + v_ij(t)δt, respectively. Simultaneously evolving all (i, j) sets, we construct the time-dependent electron density , from which the surface hole distribution is also updated at each time t using the method described in Sec. 5.4.

When electrons move back to the surface, such that r_ij(t + δt) is located inside the metal at time t + δt (but outside at time t), we introduce the effect of electron reabsorption and partial reflection by considering that a fraction of the arriving electrons is specularly reflected. This is done by inverting the sign of the surface-normal component of the velocity and making , where is the normal electron energy and is a transmittance coefficient given by eqn (7). Since , this procedure produces a depletion in the number of plasma electrons (i.e., reabsorption).

5.4 Surface screening charge

We approximate the metal as a perfect conductor because screening has a characteristic time and length of³¹ ≲1 fs and ≲1 nm, much smaller than the spatiotemporal scales involved in the formation and evolution of the plasma. The screening charge is then obtained by an adaptation of the boundary-element method for perfect conductors.⁴¹ Taking a structure that is translationally invariant along y, we consider the transverse profile r_s = (x_s, 0, z_s), parametrized by a variable s that evolves linearly from 0 to 1 as we go around the perimeter length P. We now consider a line of charge aligned along y, placed at a transverse position r_c = (x_c, 0, z_c), and having a charge density q per unit length. The presence of the external charge places the metal surface at a potential V, which is uniform in the limit of a perfect conductor. Then, the distribution of induced surface charges eσ_s (charge per unit of surface area) is determined by the condition
which needs to be satisfied at all positions s. The y integral of each of the fractions in this expression produces a logarithmic divergence at large distances. However, the overall divergence cancels due to the neutrality of the total external plus induced charges, and thus, upon integration, we find
Incidentally, we normalized the arguments of the logarithms by dividing by P² to obtain dimensionless quantities, but any normalization factor in these functions cancels because of charge neutrality. We solve this equation by discretizing the transverse surface profile through a set of N equally spaced points corresponding to , with , leading to the linear equation M·σ = b, where we define an N × N matrix of components , as well as the N-vectors and .

For a biased structure, the potential V is taken as a parameter (e.g., V = 0 for grounded samples) and the screening charge in the presence of such potential is just obtained upon inversion of the aforementioned N × N linear system of equations. Conversely, for electrically isolated metal structures, charge neutrality must be imposed through the equation . Then, the potential V is no longer a parameter, but rather a variable determined by the new equation, and the above matrix and vectors need to be supplemented with additional components , , M_N,N = 0, b_N = −q, and σ_N = V, thus defining an enlarged (N + 1) × (N + 1) system. In both scenarios (grounded and isolated structure), we need to invert the corresponding system of linear equations to find the surface hole distribution σ_s for a line charge q placed at r_c. At each time t along the evolution of the plasma, we then calculate the total surface charge as the superposition of those generated by all plasma elements q = −eρ_ij(t) at positions r_c = r_ij(t).

We note that the difference between the grounded and isolated scenarios is significant in structures of characteristic size comparable with the average plasma-sample separation distance. In contrast, large grounded structures in which electron emission occurs mainly from a comparatively small surface area (as is the case in this work) can be approximately taken as isolated because of the very large number of free electrons available in distant parts of the sample, which act as an effective ground reservoir. In the present work, we adopt this approximation because it allows us to consider a relatively small tip length (L = 10 μm) –thus alleviating the computational burden of our simulations– without introducing artificial effects arising from the distribution of the potential around the reduced-size sample.

5.5 Electric field calculation

The electric field generated by the evolving plasma at position r and time t is given by , where ϕ = ϕ_e + ϕ_h is the potential generated by the plasma electrons (ϕ_e) and the induced surface charges (ϕ_h). We ignore the effect of the vector potential due to the low drift velocity of the emitted electrons. Correspondingly, the electric field can be separated into the contributions arising from the emitted plasma electrons,and the associated induced surface charges,where r = (x, 0, z), r′ = (x′, 0, z′), r_s = (x_s, 0, z_s) runs over surface positions parameterized by s, P is the perimeter of the metal cross section, ρ_e and σ_s are the densities of emitted electrons and surface holes, respectively [see Fig. 2(c, e, g and i) in the main text], and t_r is the retarded time defined as t − |r − r′|/c in eqn (16a) and t − |r − r_s|/c in eqn (16b). The frequency-domain electric field can be obtained from the time-domain counterpart by means of a Fourier transform.

In the numerical implementation of the calculation of the electric field, the electron plasma distribution at each time t is discretized through a uniform grid of element size δx × δz, constructed such that all electrons placed inside each grid element are assimilated to a single effective charge, with a linear (along y) density given by the sum of those associated with the enclosed electrons. Analogously, the surface is also discretized with elements of equal length δs along the transverse surface profile, each of them containing an effective positive image charge. This procedure runs smoothly when evaluating the field (and the induced forces) at large distances by simply placing the effective charges at the center of the grid or surface elements. However, extra care needs to be taken at short distances, for which we consider each grid element to be uniformly charged. Eqn (16) are then corrected by performing the transformations

andwhere is a surface-tangent vector at the s-dependent position r_s. To evaluate these integrals, we can safely neglect the effect of the grid size on the terms and in eqn (16a) and (b), respectively, as we have and D_y ≫ Pδs; therefore, we only need to compute the integralsandwhere we define x_± = x ± δx/2, z_± = z ± δz/2, ,

g(x,z) = x arctan(z/x) + (z/2)log(x² + z²), (19a)

andfinally, this correction is applied by using eqn (18a) and (b) to replace the terms (r − r′)/|r − r′|² and (r − r_s)/|r − r_s|² in eqn (16a) and (16b), respectively.

We additionally note that, due to the finite time step δt adopted in the dynamical evolution of our model, the frequency-domain electric field can be resolved up to a maximum frequency ω_c = π/δt.⁴² To avoid significant aliasing effects, δt must be chosen small enough such that the field spectral components are concentrated in the ω ≪ ω_c range.

5.6 Energy variation of the probe electron

When an electron with energy E₀ and velocity v_e passes by the vicinity of the metal structure, it interacts with the plasma and the sample, thus undergoing a variation in energy by an amountwhere is the electric field on the trajectory of the electron r_e(t), calculated using eqn (16), we define Γ(t) as the electron energy variation rate, and we adopt the nonrecoil approximation (i.e., v_e is constant and |ΔE| ≪ E₀). The classical energy change represented by eqn (20) is a good approximation even when considering electrons as quantum wavepackets, as shown in ref. 23 and 24. The e-beam energy variation rates can equally be separated into the corresponding contributions as . To account for the finite FWHM of the electron wavepacket Δt_e, we correct eqn (16a) and (b) by performing a Gaussian convolution with the same FWHM, such that Γ(t) is replaced by with .

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been supported in part by the European Research Council (Advanced Grant 789104-eNANO and Starting Grant 851780-NanoEP), the European Commission (Horizon 2020 Grant 964591-SMART-electron), the Spanish MICINN (PID2020-112625GB-I00 and Severo Ochoa CEX2019-000910-S), Google Inc., the Catalan CERCA Program, and Fundaciós Cellex and Mir-Puig.

References

1. P. U. Jepsen, D. G. Cooke and M. Koch, Laser Photon. Rev., 2011, 5, 124–166.
2. T. Kampfrath, K. Tanaka and K. A. Nelson, Nat. Photonics, 2013, 7, 680–690.
3. R. Ulbricht, E. Hendry, J. Shan, T. F. Heinz and M. Bonn, Rev. Mod. Phys., 2011, 83, 543.
4. M. Nagel, M. Först and H. Kurz, J. Phys.: Condens. Matter, 2006, 18, S601.
5. Y. Watanabe, K. Kawase, T. Ikari, H. Ito, Y. Ishikawa and H. Minamide, Appl. Phys. Lett., 2003, 83, 800–802.
6. A. Dobroiu, C. Otani and K. Kawase, Meas. Sci. Technol., 2006, 17, R161.
7. S. Nakajima, H. Hoshina, M. Yamashita, C. Otani and N. Miyoshi, Appl. Phys. Lett., 2007, 90, 041102.
8. T. Nagatsuma, G. Ducournau and C. C. Renaud, Nat. Photonics, 2016, 10, 371–379.
9. V. P. Wallace, E. MacPherson, J. A. Zeitler and C. Reid, J. Opt. Soc. Am. A, 2008, 25, 3120–3133.
10. C. Yan, B. Yang and Z. Yu, Anal. Lett., 2013, 46, 946–958.
11. P. Zhao, S. Ragam, Y. J. Ding and I. B. Zotova, Appl. Phys. Lett., 2011, 98, 131106.
12. Y. Jiang, K. Vijayraghavan, S. Jung, F. Demmerle, G. Boehm, M. C. Amann and M. A. Belkin, J. Opt., 2014, 16, 094002.
13. X.-C. Zhang, X. F. Ma, Y. Jin, T.-M. Lu, E. P. Boden, P. D. Phelps, K. R. Stewart and C. P. Yakymyshyn, Appl. Phys. Lett., 1992, 61, 3080–3082.
14. A. Rice, Y. Jin, X. F. Ma, X.-C. Zhang, D. Bliss, J. Larkin and M. Alexander, Appl. Phys. Lett., 1994, 64, 1324–1326.
15. J. A. Fülöp, L. Pálfalvi, G. Almási and J. Hebling, Opt. Express, 2010, 18, 12311–12327.
16. I. G. Savenko, I. A. Shelykh and M. A. Kaliteevski, Phys. Rev. Lett., 2011, 107, 027401.
17. Z. Fang, H. Wang, X. Wu, S. Shan, C. Wang, H. Zhao, C. Xia, T. Nie, J. Miao, C. Zhang, W. Zhao and L. Wang, Appl. Phys. Lett., 2019, 115, 191102.
18. H. Hamster, A. Sullivan, S. Gordon, W. White and R. W. Falcone, Phys. Rev. Lett., 1993, 71, 2725.
19. W. P. Leemans, C. G. R. Geddes, J. Faure, C. Tóth, J. Van Tilborg, C. B. Schroeder, E. Esarey, G. Fubiani, D. Auerbach, B. Marcelis and M. A. Carnahan, Phys. Rev. Lett., 2003, 91, 074802.
20. L. Zhang, A. Tcypkin, S. Kozlov, C. Zhang and X.-C. Zhang, Ultrafast Sci., 2021, 2021, 9892763.
21. P. Dombi, Z. Pápa, J. Vogelsang, S. V. Yalunin, M. Sivis, G. Herink, S. Schäfer, P. Groß, C. Ropers and C. Lienau, Rev. Mod. Phys., 2020, 92, 025003.
22. G. Liao, Y. Li, H. Liu, G. G. Scott, D. Neely, Y. Zhang, B. Zhu, Z. Zhang, C. Armstrong, E. Zemaityte, P. Bradford, P. G. Huggard, D. R. Rusby, P. McKenna, C. M. Brenner, N. C. Woolsey and W. Wang, Proc. Natl. Acad. Sci., 2019, 116, 3994–3999.
23. I. Madan, E. J. C. Dias, S. Gargiulo, F. Barantani, M. Yannai, G. Berruto, T. LaGrange, L. Piazza, T. T. A. Lummen, R. Dahan, I. Kaminer, G. M. Vanacore, F. J. García de Abajo and F. Carbone, ACS Nano, 2023, 17, 3657–3665.
24. M. Yannai, R. Dahan, A. Gorlach, Y. Adiv, K. Wang, I. Madan, S. Gargiulo, F. Barantani, E. J. C. Dias, G. M. Vanacore, N. Rivera, F. Carbone, F. J. García de Abajo and I. Kaminer, ACS Nano, 2023, 17, 3645–3656.
25. J. Vogelsang, G. Hergert, D. Wang, P. Groß and C. Lienau, Light: Sci. Appl., 2018, 7, 55.
26. G. Hergert, A. Woste, J. Vogelsang, T. Quenzel, D. Wang, P. Gross and C. Lienau, ACS Photonics, 2021, 8, 2573–2580.
27. A. Ryabov and P. Baum, Science, 2016, 353, 374–377.
28. S. Sun, X. Sun, D. Bartles, E. Wozniak, J. Williams, P. Zhang and C.-Y. Ruan, Struct. Dyn., 2020, 7, 064301.
29. M. Centurion, P. Reckenthaeler, S. A. Trushin, F. Krausz and E. E. Fill, Nat. Photonics, 2008, 2, 315–318.
30. F. J. García de Abajo and V. Di Giulio, ACS Photonics, 2021, 8, 945–974.
31. F. J. García de Abajo and P. M. Echenique, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 13399–13407.
32. O. L. Krivanek, N. Dellby, J. A. Hachtel, J.-C. Idrobo, M. T. Hotz, B. Plotkin-Swing, N. J. Bacon, A. L. Bleloch, G. J. Corbin, M. V. Hoffman, C. E. Meyer and T. C. Lovejoy, Ultramicroscopy, 2019, 203, 60–67.
33. R. Yu, Q. Guo, F. Xia and F. J. García de Abajo, Phys. Rev. Lett., 2018, 121, 057404.
34. F. J. García de Abajo and A. Howie, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 115418.
35. The Electronic Band Structure of Copper, https://lampx.tugraz.at/hadley/ss2/bands/dft/calculations/Cu.php, Accessed: 2022-06-06.
36. Thermal Conductivity of Copper, https://www.efunda.com/materials/elements/TC_Table.cfm?Element_ID=Cu, Accessed: 2022-03-09.
37. Z. Lin, L. V. Zhigilei and V. Celli, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 075133.
38. S. O. Kasap, Principles of Electronic Materials and Devices, McGraw-Hill, New York, 2006, vol. 2.
39. G. Ferrini, F. Banfi, C. Giannetti and F. Parmigiani, Nucl. Instrum. Methods Phys. Res., 2009, 601, 123–131.
40. L. V. Keldysh, Sov. Phys., 1965, 20, 1307–1314.
41. S. Thongrattanasiri, I. Silveiro and F. J. García de Abajo, Appl. Phys. Lett., 2012, 100, 201105.
42. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes, Cambridge University Press, New York, 1992.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3na00168g

---