Terahertz emission from laser-driven gas-plasmas: a plasmonic point of view

We disclose an unanticipated link between plasmonics and nonlinear frequency down-conversion in laser-induced gas-plasmas. For two-color femtosecond pump pulses, a plasmonic resonance is shown to broaden the terahertz emission spectra significantly. We identify the resonance as a leaky mode, which contributes to the emission spectra whenever electrons are excited along a direction where the plasma size is smaller than the plasma wavelength. As a direct consequence, such resonances can be controlled by changing the polarization properties of elliptically-shaped driving laser pulses. Both, experimental results and 3D Maxwell consistent simulations confirm that a significant terahertz pulse shortening and spectral broadening can be achieved by exploiting the transverse driving laser beam shape as an additional degree of freedom.


DEFINITIONS: FOURIER TRANSFORMS AND FAR-FIELD SPECTRA
In the following, we define the temporal Fourier transform f (r, ω) of a function f (r, t) bŷ Furthermore, we define the longitudinal spatial Fourier transformf (r ⊥ , k z , ω) of a functionf (r ⊥ , z, ω) by f (r ⊥ , z, ω) = f (r ⊥ , k z , ω)e ik z z dk z .
Note the difference in sign of the exponent for temporal and spatial transforms, which is common practice in the optical context. We introduce the spectral poynting fluxeŝ andS (r ⊥ , k z , ω) = 2/µ 0 {Ȇ(r ⊥ , k z , ω) ×B (r ⊥ , k z , ω)} , where µ 0 is the magnetic permeability and denotes the real part. The first expression is used to compute the angularly integrated far-field power spectrum in Maxwell-consistent 2D/3D simulations by integration over a closed surface around the plasma. To compute the angularly integrated far-field spectrum P x in Sec. 3 of the main article, we integrateS x along k z .

MODELING THE IONIZATION CURRENT MECHANISM FOR MICROPLASMAS
The THz generation by two-color laser pulses considered here is driven by the so-called ionization current (IC) mechanism [1]. A comprehensive model based on the fluid equations for electrons describing THz emission has been derived in [2]. In this framework, the IC mechanism naturally appears at the lowest order of a multiple-scale expansion. Besides the IC mechanism, this model is also able to treat THz generation driven by ponderomotive forces and others. They appear at higher orders of the multiple-scale expansion and are not considered here. In the following, we briefly summarize the equations describing the IC mechanism. The electromagnetic fields E and B are governed by Maxwell's equations in vacuum The plasma and electromagnetic fields are coupled via the conductive current density J governed by where electron-ion collisions lead to a damping of the current. The collision frequency is determined by [3] where λ ei is the Coulomb logarithm and E elec denotes the thermal and kinetic electron energy. The value λ ei = 3.5 turned out to match the results obtained by more sophisticated calculations with Particle-In-Cell (PIC) codes in [2]. The densities of Z times charged ions are determined by a set of rate equations for Z = 1, 2, 3, . . . , K, and the initial neutral density is n (0) ion (t = −∞) = n a . The tunnel ionization rate W (Z) in quasi-static approximation creating ions with charge Z is taken from [4,5]. Thus, W (Z) is a function of the modulus of the electric field E. The atoms can be at most K times ionized and thus W (K+1) = 0. The electron density is determined by the ion densities (S12) The electron energy density E = n e E elec is governed by The loss current accounting for ionization losses in the laser field reads where I Z p is the ionization potential for creation of an ion with charge Z. Even though ionization losses as well as higher order ionization (Z = 2, 3, . . .) are negligible in the framework of the present study, both are kept for completeness.
The model is implemented in the code ARCTIC, that solves Eqs. (S7)-(S8) by means of the Yee scheme [6]. We have previously benchmarked the code ARCTIC by the PIC code OCEAN [7] accounting for full kinetics of the plasma in [8].

PLASMONIC RESONANCES IN THE PLASMA SLAB MODEL
To elaborate the origin of the laser-polarization-dependence of the THz-emission spectra from elliptically shaped 2C-laserinduced plasmas, we consider a simplified system as sketched in Fig. 5(a) of the main article: A plasma slab of thickness d in x direction with time-invariant electron density n 0 and collision frequency ν ei . Both quantities are translationally invariant in y and z. Above and below the slab we assume semi-infinite vacuum with a constant (relative) permittivity v = 1. By writing the total electric field as the sum of the laser field E L and the field due to laser-plasma interactionẼ, Eq. (2) in the main article reads as with the source term It is important to note that we assume that the source term ι depends on the product of the laser electric field and the timedependent electron density n e (t), as it is produced during the laser gas interaction. Only in the description of the response of the plasma slab we make the simplification of a time-invariant density n 0 .
Let us now use Eq. (S15) and Maxwell's equations to determine the response of the system. In frequency space (see Sec. 1 for definition), they can be rewritten for angular frequency ω = 0 as where for sake of brevity we introduced the source term The dielectric permittivity reads = p in the plasma slab (|x| ≤ d/2), and = v = 1 in the vacuum (|x| > d/2). The complex dielectric permittivity of the plasma is given by where the plasma frequency ω p = n 0 q 2 e m e 0 involves the time independent density n 0 of the slab. The time-independent collision frequency ν ei can be disregarded (ν ei = 0) and is kept here just for completeness.
Same as for , we consider an excitation that is translational invariant in y, that is, ∂ yι = 0 and ∂ yQ = 0. Therefore, we can set all the y-derivatives to zero and Eqs. (S17)-(S18) separate into two sets of equations. The translational invariance of the slab in z allows to write down these two sets of equations in the spatial Fourier domain with respect to z (∂ z → ik z ) giving where "˘" indicates the (ω, k z )-domain according to the definition in Sec. 1. The 1 st set of equations (S21)-(S23) is associated with the so-called S polarization. In the waveguide context, it is also often termed the transverse electric (TE) mode, because the only electric field component E y is polarized in the transverse translational invariant direction. Here, the only fields different from zero are (B x ,Ȇ y ,B z ). The 2 nd set of equations (S24)-(S26) is associated with the so-called P polarization and describes the evolution of (Ȇ x ,B y ,Ȇ z ). In the waveguide context it is also often termed the transverse magnetic (TM) mode.
In the following, we will consider two different configurations: Note that (i) corresponds to the THz generation by y-polarized "elliptical beams" and (ii) by x-polarized "elliptical beams" as investigated in the main article. Before computing the response of the plasma slab to an excitation by ι, it is interesting to investigate the resonances of the homogeneous system (ι = 0). To this end, it is convenient to solve the reflection transmission problem and to consider the reflection coefficient. Singularities of the reflection coefficient describe resonant modes, i.e., modes which provide non-zero fields in the slab without an incoming wave and without a source term. When then exciting the system by the source ι, these modes are most likely excited and thus of great interest.

A. Reflection coefficient of a slab
In the following, we consider the problem of computing the reflection coefficient for a plane wave arriving from one of the vacuum half-spaces and interacting with a plasma slab (see Fig. S1). To this end, we will operate in the (ω, k z )-space. We will detail the P polarization case, because it is of most interest. The S polarization case follows analogously with a small modification which will be highlighted in the coming derivation.
We denote byB i the amplitude of the magnetic field of the incident P polarized wave and byB r the amplitude of the reflected wave. Without loss of generality, the incident wave arrives from the positive half-space. Then, the magnetic field in the vacuum Fig. S1. Illustration of the plasma slab. For the magnetic field components, see text.
Denoting its x-derivative byB y = ∂ xBy , this implies the relation where we introduced the transformation matrix K.
For P polarization, only the transverse fieldB y is continuous at an interface while for the derivative we can only use that −1 ∂ xBy is continuous. The situation is easier for S polarization where bothȆ y andȆ y are continuous at an interface. Thus, to handle the P polarization and later on the S polarization case in the same manner, we introduce α P = 1/ p and α S = 1. Then, we relate the fields at the vacuum plasma interface at where we introduced the interface transition matrix T. The magnetic field in the plasma (|x| < d/2) can be decomposed into a forwardB + and a backwardB − running wave asB where Λ p = k 2 z − p ω 2 /c 2 . Using this equation, a short computation relates fields at the two interfaces inside the plasma In complete analogy to Eq. (S29), we obtain the fields in vaccum at the lower interface at x = −d/2 by We denote byB t the amplitude of the transmitted wave. Because there is no incoming wave from the negative half-space, the field in the vacuum for x ≤ d/2 reads In summary, we therefore obtain After performing all the matrix multiplications, we end up two equations relating the field amplitudesB i ,B r , andB t . We can thus eliminate the transmitted field amplitudeB t and get the complex valued reflection coefficient This derivation can be straightforwardly repeated for the S polarization case, where we governȆ y instead ofB y and just have to replace α P by α S . Then, we obtain the intensity reflectivity ρ S/P := |r S/P | 2 as This is exactly Eq. (2) of the main article.

B. S polarization with transverse excitation in y
In the following, we consider the full inhomogeneous set of Eqs. (S21)-(S26) including the excitation source term ι. While in the general case these equations have to be solved numerically, for the plasma slab they can be solved even analytically as it will be shown in the following. Let us start with case (i), that is,ι =ι y e y . Firstly, the general solutions inside the plasma slab and the neutral gas are computed, and secondly, the continuity of the transverse fields is used to determine the entire solution.
Equations (S21)-(S23) give the evolution equation for the transverse fieldȆ y inside the plasma and the neutral gas For sake of simplicity, we consider thatQ y is constant inside the plasma slab and zero outside of it. In particular, this choice makesQ y symmetric with respect to x, andȆ y has to obey the same symmetry. Then, the general symmetric solution in the plasma (|x| ≤ d/2) reads In the vacuum at x > d/2 the general solution reads For k 2 z ≥ ω 2 /c 2 , we use the "−" sign since for the "+" sign the field would grow exponentially when x → ∞. For k 2 z < ω 2 /c 2 , we use the "+" sign to obtain only outgoing propagating waves (along x). Note that solution in the vaccum at x < d/2 follows from symmetry considerations.
According to Maxwell's interface conditionsȆ y ,B z and thus ∂ xȆy have to be continuous at the plasma-vacuum interface. These two conditions determine the yet unknown amplitudes A p and A v to with common denominator Finally, Eqs. (S21)-(S22) determine the magnetic field components asB Equations (S40) and (S46) were used in Sec. 3 of the main article to compute the transverse Poynting fluxes via Eq. (S6).

C. P polarization case with transverse excitation in x
Next, case (ii) withι =ι x e x is considered. Equations. (S24)-(S26) give the evolution equation for the transverse fieldB y inside the plasma and neutral gas respectively, In analogy to the S polarization case in the previous section, we obtain in the plasmȃ and in upper vacuum The difference to the S polarization case appears when applying the interface conditions at the plasma-air interface:B y andȆ z are continuous but according to Eq. (S26) ∂ xBy is not, because changes at the interface. Applying these P polarization interface conditions determines A p and A v to with common denominator Finally, Eqs. (S25)-(S26) determine the electric field components asȆ

THZ EMISSION FROM ELLIPTICALLY SHAPED MICROPLASMAS
In the main article we investigate terahertz (THz) emission from elliptically shaped two-color (2C)-laser-induced gas-plasmas.
In such a configuration, the free electron density profile with a small transverse size along x and a large transverse size along y is created. For the considered microplasmas, the transverse plasma profile is strongly elliptical, that is, along x direction the plasma size is less than 1 µm, whereas along y direction the plasma is approximately 10 µm wide. Thus, by rotating the linear laser electric field polarization with respect to the transverse plasma profile, we can select the strength of the electron density gradients along the excitation direction.
As it can be seen in the 3D angularly integrated far-field spectra presented in Fig. S2, when exciting the plasma along strong electron density gradients (x-polarized laser associated with P polarization), the THz spectrum is broadened up to about 50 THz (dark red solid line), which corresponds to the maximum plasma frequency ν max p . In contrast, when exciting the plasma along the weak electron density gradients (y-polarized laser associated with S polarization) no such broadening is found (light gray solid line).
Results of corresponding 2D simulations (∂ y = 0) are shown as dashed lines in Fig. S2. Here, we find a similar behavior as in 3D: no broadening if the laser electric field is oriented in the now translationally invariant y direction (S), and broadening up to ν max p if the laser electric field points in the direction of the strong electron density gradient, that is, along the x direction (P). Treating the problem in 2D geometry, i.e., assuming translational invariance in one transverse direction (here y), the electromagnetic fields separate into two cases: S polarization case that governs the fields B x , E y , B z for a y-polarized driving laser pulse and the P polarization case that governs the fields E x , B y , E z for an x-polarized driving laser pulse. Any other polarization state in 2D can be written as the superposition of these two cases, if coupling through the laser generated plasma is ignored. For example, an incoming laser pulse that is linearly polarized under 45 • in the xy plane will give an electric field solution that can be written as E = E S / √ 2 + E P / √ 2, where E S and E P are the solutions for an S and a P polarized driving laser pulse, respectively. Therefore, if the 3D elliptical beam can be indeed approximated by the idealized 2D case, and no detrimental nonlinear coupling occurs, this property should hold. We checked this by comparing the angularly integrated THz-far-field power spectrum for a simulation with polarization at 45 • (black solid line) and the result for the superposed fields (black dashed line) in Fig. S2. Both overlap almost perfectly, which further justifies our analysis of the 2D configuration. Moreover, the possibility of superposing the THz fields which are produced by an xand y-polarized laser pulse could be important for applications, because it implies that the THz emission spectrum can be tuned by rotating the linear polarization of the incoming elliptical laser pulse.