Efficient generation of extreme terahertz harmonics in 3D Dirac semimetals

Frequency multiplication of terahertz signals on a solid state platform is highly sought-after for the next generation of high-speed electronics and the creation of frequency combs. Solutions to efficiently generate extreme harmonics (up to the $31^{\rm{st}}$ harmonic and beyond) of a terahertz signal with modest input intensities, however, remain elusive. Using fully nonperturbative simulations and complementary analytical theory, we show that 3D Dirac semimetals (DSMs) have enormous potential as compact sources of extreme terahertz harmonics, achieving energy conversion efficiencies beyond $10^{-5}$ at the $31^{\rm{st}}$ harmonic with input intensities on the order of $10$ MW/cm$^2$, over $10^5$ times lower than in conventional THz high harmonic generation systems. Our theory also reveals a fundamental feature in the nonlinear optics of 3D DSMs: a distinctive regime where higher-order optical nonlinearity vanishes, arising as a direct result of the extra dimensionality in 3D DSMs compared to 2D DSMs. Our findings should pave the way to the development of efficient platforms for high-frequency terahertz light sources and optoelectronics based on 3D DSMs.

Frequency multiplication of terahertz signals on a solid state platform is highly sought-after for the next generation of high-speed electronics and the creation of frequency combs. Solutions to efficiently generate extreme harmonics (up to the 31 st harmonic and beyond) of a terahertz signal with modest input intensities, however, remain elusive. Using fully nonperturbative simulations and complementary analytical theory, we show that 3D Dirac semimetals (DSMs) have enormous potential as compact sources of extreme terahertz harmonics, achieving energy conversion efficiencies beyond 10 −5 at the 31 st harmonic with input intensities on the order of 10 MW/cm 2 , over 10 5 times lower than in conventional THz high harmonic generation systems. Our theory also reveals a fundamental feature in the nonlinear optics of 3D DSMs: a distinctive regime where higher-order optical nonlinearity vanishes, arising as a direct result of the extra dimensionality in 3D DSMs compared to 2D DSMs. Our findings should pave the way to the development of efficient platforms for high-frequency terahertz light sources and optoelectronics based on 3D DSMs.
High-harmonic generation (HHG) is a nonlinear process involving the emission of light at integer multiples of the driving laser frequency. HHG in gaseous media is a well-established method of generating high-frequency light and broadband frequency combs, which provide attosecond resolution in the study of atomic and materials phenomena [1][2][3][4]. More recently, HHG has been demonstrated in solids [5][6][7][8][9][10][11], revealing its potential for realizing chip-integrable optoelectronics and compact radiation sources. However, solid-state HHG typically requires strong driving fields exceeding 1 GV/m (i.e., on the order of TW/cm 2 peak power) [5][6][7][8][9][10][11][12]. With rapidly growing interest in the terahertz (THz) regime as the answer to next-generation computation, imaging and communications, the ability to perform efficient THz HHG using modest input intensities on a solid-state platform is a highly sought-after solution.
Here, we show that the recently discovered class of 3D Dirac semimetal (DSM) materials is a promising candidate for solid-state THz HHG, achieving high efficiencies with input field intensities on the order of 10 MW/cm 2 , more than 10 5 times weaker than required for conventional THz HHG. This applies not only to the generation of lower-order harmonics [13][14][15][16][17], but also to the generation of extreme harmonics, up to the 31 st order and beyond. In particular, by considering a Cd 3 As 2 thin film of 250 nm thickness (readily obtained using methods like molecular beam epitaxy [14,16,[18][19][20][21][22]), we predict that conversion efficiencies exceeding 10 −5 can be achieved as high as the 31 st harmonic. Notably, we find that the 31 st harmonic peak can be brought within 5 orders of magnitude of the driving harmonic using an input field of 10 MV/m and below, well within reach of existing table-top THz sources [23][24][25][26][27]. In contrast, conventional solid-state THz HHG requires driving fields approaching 10 GV/m -both in theory and in experiment -just to bring the 22 nd harmonic peak within 8 orders of magnitude of the driving harmonic [6].
Our studies also reveal a novel regime where the higher order nonlinearities of the intraband current in 3D DSMs completely vanish. As a result, the efficiency of HHG beyond the 3 rd harmonic is greatly suppressed in this regime. As we show, this is a fundamental feature of 3D DSMs directly arising from their extra dimensionality compared to 2D DSMs, which have no such regime of suppression. This breaks the common notion that 2D DSMs -well known for their nonlinear optical response [28][29][30][31][32][33][34] -share the same essential physics as 3D DSMs. We identify a physical quantity -termed the critical field strength -that divides the regime of nonlinear suppression from the regime of extreme nonlinearity (and hence highly efficient HHG) in 3D DSMs. We show that in the latter regime, 3D DSMs have orders-of-magnitude performance enhancements compared to existing solid-state HHG platforms.
RESULTS. A laser pulse impinging on a DSM material induces carrier oscillations within (intraband current) and carrier transitions between (interband current) the upper and lower bands of the Dirac cone band structure. These carrier oscillations and transitions in turn emit light peaked at multiples of the driving laser frequency (Fig. 1a).
The scenario we consider is illustrated in the inset of Fig. 1b: An external laser pulse is normally incident on a DSM material, which emits multiple harmonics of the incident pulse. Maxwell's equations are used to model the electrodynamics of the pulse and its interaction with the DSM material. The properties of the DSM material are in turned determined using nonpertubative time-domain quantum simulations. Full details of the treatment are FIG. 1. Highly efficient generation of extreme harmonics (up to 31 st harmonic) in 3D DSM Cd3As2 at modest driving field strengths. (a) HHG in a 3D DSM occurs when a driving laser pulse (red-arrowed waveform) induces carrier oscillations (current density depicted at one instant in time on the Dirac cone graphic) and transitions which lead to the emission of high harmonic light (multi-colored output arrows). Driving a Cd3As2 thin film with a linearly polarized pulse of central frequency 1 THz and peak field strength 10 MV/m ((b) inset) produces the emitted spectrum shown in (b), where we see that harmonics up to the 31 st order and beyond can be generated at energy conversion efficiencies well in excess of 10 −5 . (c)-(g) show the change in the output energy spectrum and conversion efficiency as a function of externally incident driving field strength. These results reveal that the generation of extreme harmonics continue to remain relatively efficient even at much lower field strengths. The markers in (c)-(g) are the result of our numerical simulations, whereas the connecting lines are visual guides. We considered a 250 nm thick Cd3As2 thin film of radius 1 mm and Fermi energy EF = 60 meV, uniformly illuminated by a 2 ps long pulse of 1 THz peak frequency.
given in Methods, with the key points as follows.
To model the dynamics in time t of an electron in a Dirac cone band structure, we start with the effective Hamiltonian where is the reduced Planck constant, v j are the Fermi velocities along Cartesian directions j, σ j are the Pauli matrices, and p j are the initial electron momenta, with j ∈ {x, y, z} and j ∈ {x, y} for 3D and 2D DSMs, respectively. We include the field via a modified minimal coupling substitution p p p → π π π(t) = p p p+ea a a(t), where e is the elementary charge, and a a a(t) is the modified vector potential inside the DSM, which is related to the electric field E E E(t) by the relation a a a(t) = −e −t/τ t −∞ E E E(t )e t /τ dt [35], where the inelastic intraband scattering time is τ . In the absence of intraband scattering, i.e., τ → ∞, a a a(t) is exactly the vector potential A A A(t). From Eq. (1), we obtain the total induced current J i (t), in terms of the popula-tion inversion N p p p (t) and the interband coherence Γ p p p (t), as where g = 4 is the combined valley and spin degeneracy and the integral extends over all momentum space. For 3D DSMs, we take n = 3 and i ∈ {x, y, z}; for 2D DSMs, we have n = 2 and i ∈ {x, y}. We also define the instantaneous energy E(t) = i v 2 i π 2 i (t), (Λ x , Λ y , Λ z ) = (cos θ cos φ, cos θ sin φ, − sin θ), and (∆ x , ∆ y , ∆ z ) = (− sin φ, cos φ, 0), where θ(t) = arccos[v z π z (t)/E(t)], and φ(t) = arctan[v y π y (t)/v x π x (t)]. For 2D DSMs, we set θ = π/2 and J z = v z = 0, obtaining an expression that reduces to the special case of isotropic 2D DSMs in [31,[35][36][37][38] when we further set v x = v y = v F . In Eq.
(2), the first term of the integrand represents the contri-bution of intraband current. The remaining second and third terms of the integrand represent the contribution of interband current. The effect of finite temperature and interband carrier scattering are taken in account in the solutions to N p p p (t) and Γ p p p (t).
We evaluate Eq. (2) to obtain a fully-closed-form, nonperturbative expression for the total current that is valid when the intraband current dominates the material response, as is the case in our regimes of interest. In 3D DSMs, the form of the total current depends on the value of the parameter , i ∈ {x, y, z}, and E F is the Fermi level. We refer to the field strength that exactly satisfies eΦ/E F = 1 as the critical field strength. When eΦ(t)/E F < 1, we obtain the x-component of the current as J 3D,sub where i ∈ {x, y, z}. When eΦ(t)/E F > 1, the expression becomes Equations (3) and (4) were obtained assuming that the driving frequencies ω satisfy the relation ω 2E F (intraband current dominantes), and that the initial Fermi distribution f D (E) is well-approximated by f D (E) = 1 when E < E F , and f D (E) = 0 when E > E F (a condition that holds exactly when T = 0 K). Details of the derivation are presented in Supplementary Information (SI) Section IA. Throughout this work, we will term eΦ max /E F < 1 the subcritical regime, and eΦ max /E F > 1 the supercritical regime, where Φ max is the maximum amplitude of Φ(t). It is noteworthy that that the first-and third-order conductivities of 3D DSMs extracted from Eq. (3) agree with known results [15,39]. At the same time, we note that the full, non-perturbative response of 3D DSMs has never been presented until now in equations (3) and (4).
The HHG energy spectrum presented in Fig. 1b indicates that extremely efficient intraband HHG up to the 31 st order and beyond can be generated with significant efficiencies at modest driving laser energies. We consider a 250 nm thick Cd 3 As 2 thin film (readily grown using molecular beam epitaxy [14,16,[18][19][20][21][22]) of Fermi level E F = 60 meV, driven by a normally incident, 2 ps-long pulse of peak frequency 1 THz and peak field amplitude E 0,x = 10 MV/m (Fig. 1b inset). We used experimentally verified values of the Fermi velocities for Cd 3 As 2 : (v x , v y , v z ) = (1.28, 1.3, 0.33) × 10 6 m/s [41]. We find that the 31 st harmonic of the output energy spectral density lies within 5 orders of magnitude of the fundamental harmonic. This compares favorably with the performance of conventional solid-state THz HHG, where the 22 nd harmonic lies within 8 orders of magnitude of the fundamental harmonic, when using input intensities 10 5 times stronger than what we consider [6]. As we operate in the supercritical regime eΦ(t) E F , Eq. (4) shows that the intraband current approaches J 3D,sup , which describes a square-wave temporal profile containing only odd-ordered frequency components. In Fig. 1, we consider T = 0 K and the absence of carrier scattering. As shown in SI Section II, the high output intensities at higher harmonics persist for higher temperatures and non-zero carrier scattering. Figures 1c-1g, which show the output intensities of various harmonic orders generated by the same Cd 3 As 2 thin film when the incident peak field strength is varied between 2 MV/m and 10 MV/m, reveal that efficient HHG performance at extreme harmonics can still be accessed even at more modest field strengths. We define energy conversion efficiency as U N /U in where U in is the incident laser energy (considering only the part of the pulse that interacts with the sample) and U N is the output energy of the N th harmonic, obtained by integrating the energy spectrum over a bandwidth of ω 0 about the harmonic peak at N ω 0 . For each harmonic peak, we generally observe a rise followed by a fall in the energy spectral density as driving electric field strengths increase, indicative of the saturation of lower harmonics at larger driving fields.
Our studies also reveal the existence of a novel regime in 3D DSMs where where higher order harmonic emission is instead suppressed. As can be seen directly from Eq. (4), the intraband current in the subcritical regime (eΦ max /E F < 1) is made up purely of the 1 st and 3 rd harmonics. As such, the emission due to the intraband current -which we term intraband emission -also consists only of the 1 st and 3 rd harmonics (Fig. 2a). This is a noteworthy feature of 3D DSMs, especially since 2D DSMs exhibit no such vanishing of higher order harmonics under any condition (Fig. 2b). To further illustrate this, the emission spectra of 3D and 2D DSMs in the subcritical regime is shown in Figs. 2c,e and 2d,f respectively. This difference in 3D and 2D DSM behavior can also be seen by directly comparing the intraband current for 2D DSMs, given by Supplementary Eq. (S22) -which contains every order of nonlinearity -in SI Section IB with the intraband current for 3D DSMs in Eq. (4).
The reason for the vanishing of higher-order nonlinearities in 3D DSMs lies in the extra dimensionality 3D DSMs possess compared to 2D DSMs. Due to the extra dimensionality, the integral in Eq. (2) takes place over 3 dimensions in momentum space for 3D DSMs -resulting in the spherical region of integration shown in Fig. 2g as opposed to just 2 dimensions in momentum space for 2D DSMs (Fig. 2h). It is also possible to mathemati-FIG. 2. Zero higher-order intraband emission from 3D DSMs at subcritical field strengths, leading to strongly suppressed nonlinearities (including HHG) in this regime. Near or below the critical field strength, 3D DSMs (a) cannot efficiently generate harmonics beyond the third order, unlike 2D DSMs (b). As noted in the text, however, 3D DSMs can achieve energy conversion efficiencies orders of magnitude beyond 2D DSMs and other conventional solid-state HHG platforms. The regime of nonlinear suppression is thus an important fundamental aspect of 3D DSMs, but does not prevent 3D DSMs from emitting high harmonics efficiently under the right conditions (as we see in Fig. 1). We show the numerically computed spectra in (c) and (d) for a 3D DSM and a 2D DSM respectively. As the driving field strength is less than the critical field strength, the intraband current in 3D DSMs contain no nonlinearities beyond the 3 rd harmonic, as we see from the output spectrum in (c). In contrast, the output from a 2D DSM (d) contains all harmonics. That this phenomenon holds for all field strengths below the critical field strength is illustrated in (e) and (f). The interband current in 3D DSMs does give rise to output beyond the 3 rd harmonic, but its contribution is very weak and has been verified to fall below the range of plotted intensities. In (e) and (f), we also see that the results of our fully closed-form, nonperturbative expressions (solid curves) are in excellent agreement with rigorous numerical simulations (solid circles). The reason for the vanishing of higher-order nonlinearities in 3D DSMsalthough no such phenomena occurs in 2D DSMs -lies in the extra dimensionality 3D DSMs possess compared to 2D DSMs. Although 3D and 2D DSMs share the same expression for current Eq. (2), the required integration over 3D momentum space (g) for 3D DSMs as opposed to integration over 2D momentum space (h) for 2D DSMs, leads to very different behavior in these materials. For all cases in this figure, we consider Fermi velocities vx = vy = vz = 10 6 m/s (same as graphene's [40]), Fermi energy EF = 250 meV at temperature T = 0 K, and no carrier scattering.
cally illustrate this vanishing of all higher orders in the 3D DSM subcritical regime using a Legendre polynomial expansion, which we do in SI Section III.
In Fig. 2, we consider temperature T = 0 K and the absence of carrier scattering. However, we find that even when finite temperatures and carrier scattering effects are considered, the suppression of higher-order light emission remains significant (SI Section IV). Furthermore, we note that this suppression of HHG in the subcritical regime of 3D DSMs is a phenomenon that is robust against variations in driving field polarization and phase, and Fermi velocity anisotropy (SI Section V).
In Fig. 3, we explore HHG in both subcritical and supercritical regimes for experimentally realized 3D DSMs Cd 3 As 2 [41] and Na 3 Bi [42]. We used experimentally obtained Fermi velocities for Na 3 Bi: (v x , v y , v z ) = (4.17, 3.63, 0.95) × 10 5 m/s [42], and, unless otherwise specified, the same parameters as in Fig. 2. Once again, we observe excellent agreement between our numerical spectra (markers) and closed-form expressions given by equations (3) and (4) (curves). When driven by fields larger than the critical field strength, we observe that 3D DSMs behave in a qualitatively similar manner as graphene [35,36,43]. Note, however, that a comparison between Cd 3 As 2 and (single-layer) graphene shows that Cd 3 As 2 can surpass the power output of graphene by over 2 orders of magnitude under the same input conditions, especially at higher harmonics. Using a hypothetical multilayer graphene-dielectric structure to achieve performance comparable with Cd 3 As 2 would require at least 9 layers of graphene-dielectric, a structure that has yet to be realized (SI Section VI). Our results show that DSMs with higher Fermi velocities along the direction of the driving field polarization are able to access the supercritical regime at lower field strengths, as expected from the fact that the critical field strength scales proportionally with Φ = i v 2 i A 2 i . As a result, we see in Figs. 3a-3c that higher harmonics rapidly emerge at a lower field strength for Cd 3 As 2 as compared to Na 3 Bi and graphene since Cd 3 As 2 has the highest Fermi velocity in the direction of the driving laser polarization.
DISCUSSION. Our work reveals two distinct operation regimes in 3D DSMs: The supercritical regime where optical nonlinearities are strong and HHG efficiently generates harmonics up to the 31 st order and beyond; and the subcritical regime where nonlinearities beyond the 3 rd order completely vanish from the intraband current, resulting in greatly diminished higher harmonic output. While much work has risen around 3D DSMs as bulk versions of 2D DSMs [14,[44][45][46], our findings break the common notion that these two systems share the same essential physics. As we show, the extra dimensionality in 3D DSMs could lead to much weaker nonlinear response from 3D DSMs in the subcritical regime compared to 2D DSMs, even though 3D DSMs seem to have a larger interaction volume. Choosing parameters that put us in the supercritical regime, however, we see that 3D DSMs can efficiently generate HHG up to the 31 st harmonic and more, well beyond recent experiments that demonstrated frequency upconversion with 3D DSMs up to the 5 th [14] and 7 th harmonics [16]. Furthermore, the extreme THz HHG we study is performed at modest driving intensities, over 10 5 times lower than in conventional THz high harmonic generation systems.
Our studies also reveal the important role of anisotropy in HHG from 3D DSMs. Because the critical field strength scales as i v 2 i A 2 i , and the amplitude of the intraband current is proportional to 1/v y v z for an xpolarized input field in the deep supercritical regime, we see that the value of the Fermi velocity v v v has significant influence over the nonlinear optical properties of 3D DSMs. This motivates the development of bandstructure engineering methods that can give us greater control over the Fermi velocity of a 3D Dirac cone, for instance, to enhance the HHG intensity through an appropriate choice of Fermi velocity (as we show through an example in SI Section VII).
Our findings are also relevant to the broader study of 3D DSMs, even beyond high harmonic generation. In particular, we note that the full, non-perturbative response of 3D DSMs has never been presented until now in Eqs. (3)-(4), where we present them in fully closedform, analytical expressions. These analytical expressions also allow very convenient implementation in numerical electrodynamics solvers to study the electromagnetic response of 3D DSMs, as we do using an FDTD algorithm. With growing interest in the optical properties of unconventional topological bandstructures [13,[47][48][49], the theory and numerical implementation we present here could also potentially prove useful in studying the nonperturbative light-driven dynamics of low-energy electrons in various systems.
In summary, we show the ability of 3D DSMs to efficiently generate extreme THz harmonics up to the 31 st harmonic and beyond, using modest laser intensities on the order of 10 MW/cm 2 , which are 10 5 times weaker than in conventional solid-state HHG [5][6][7][8][9][10]. Our studies reveal two distinct operation regimes for 3D DSMs: the supercritical regime, where the extreme nonlinearity of the 3D DSM response leads to HHG energy conversion efficiencies exceeding 10 −5 at the 31 st harmonic; and the subcritical regime, where intraband emission vanishes and HHG output is greatly diminished. We further show that the vanishing of intraband emission in the subcritical regime is linked to the extra dimensionality of 3D DSMs compared to 2D DSMs, breaking the common notion that 3D DSMs are simply bulk versions of 2D DSMs. Our studies at different temperatures, scattering times, incident field polarization and material anisotropy show that our conclusions are robust under a broad range of parameters. Our work fills a vital gap in the understanding of nonlinear physics in 3D DSMs, paving the way to the development of highly efficient, chip-integrable THz light sources and optoelectronics.

Nonperturbative massless electron dynamics.
Here, we derive the ordinary differential equations (ODE) governing the coupled behavior of the interband and intraband electron dynamics in a 3D DSM. Our derivation here follows from main text Eq. (1). Using the minimal coupling substitution, the low-energy Hamiltonian reads where j ∈ {x, y, z} runs over Cartesian directions, is the reduced Planck's constant, v j are the Fermi velocities associate with direction j, σ j are the Pauli matrices, π j (t) = p j + ea j (t) are the components of the modified minimal coupling momentum, p j is the unperturbed electron momentum, e is the elementary change, a j (t) is defined as [35] a where τ is the inelastic intraband scattering time, and E j (t) the time-varying, spatially-uniform electric field. In the absence of intraband scattering, i.e., τ → ∞, a j (t) is exactly equivalent to the unmodified vector po- The fields should include the self-consistent response of the material, but for simplicity we approximate it as the external field, neglecting induced-field effects. Following previous studies on graphene [35][36][37], we write the normalized, timedependent eigenstates in the adiabatic limit: The "c" and "v" subscripts denote the conduction and valence band states respectively. We define The most general wave function that satisfies Eq. (5) can be constructed by superposing the above eigenstates: The complex coefficients C p p p,c (t) and C p p p,v (t) describe how the probability of finding the electron in either state evolves. By using Eq. (8) as an ansatz to the Schödinger Eq., i.e., i ∂ t Ψ p p p (t) =Ĥ p p p Ψ p p p (t), and defining the population inversion (difference of electron population between the valence and conduction bands), N p p p (t) = |C p p p,c (t)| 2 − |C p p p,v (t)| 2 , and the interband coherence, Γ p p p (t) = [C p p p,v (t)] * C p p p,c (t)e −2iΩ(t) , where * denotes the complex conjugate, we can derive a set of ODEs, which describe how these quantities vary in time: Here, inelastic interband damping has been introduced phenomenologically through scattering time τ . While τ in Eqs. (9a) and (9b) can be different from τ in Eq. (6), we chose to use the same value for both scattering times. At some initial time, t 0 , long before the laser pulse interacts with the electrons, we have is the Fermi-Dirac distribution for a chemical potential µ(T ) at temperature T , and k B is the Boltzmann constant. Equations (9a) and (9b) describe how a massless electron moves within (intraband) and transitions between (interband) energy bands under the influence of driving fields, unlike the semiclassical Boltzmann equation which only describes intraband motion. In deriving these ODEs, we have made no approximations beyond a constant scattering time and the massless electron limit. The equivalent equations for the 2D case, which have been successful in describing infinitely-extended graphene interacting with a laser [35][36][37], are obtained by neglecting the zcomponents and setting θ = π/2, effectively limiting interactions to within the p x -p y plane. The induced current due to a single momentum value p p p is computed as j j j p p p (t) = −eΨ † p p p (t)∇ π π πĤp p p Ψ p p p (t), where ∇ π π πĤp p p = (v x σ x , v y σ y , v z σ z ) is the group velocity operator and the † supercript denotes Hermitian conjugate. The individual current components are j p p p,x = −ev x N p p p sin θ cos φ − 2 cos θ cos φRe Γ p p p − 2 sin φIm Γ p p p (11a) j p p p,y = −ev y N p p p sin θ sin φ − 2 cos θ sin φRe Γ p p p + 2 cos φIm Γ p p p (11b) j p p p,z = −ev y N p p p cos θ + 2 sin θRe Γ p p p .
The total induced current is obtained by integrating over all momentum space: where g is the electron state degeneracy (e.g., number of Dirac cones in the first Brillouin zone) and the integral extends over the entire momentum space. When performing the integral in Eq. (12), the N p p p terms in Eqs. (11a)-(11c) must replaced by N p p p + 1, as part of a regularization procedure to prevent divergences at large momenta due to the assumption of an infinitely-extended Dirac cone [35][36][37]50]. Eq. (12) is equivalent to Eq. (2) in the main text.
Nonperturbative time-domain quantum simulations. Equations (9a) and (9b) are discretized on a scaled momentum space grid (q i = v i p i , where i ∈ {x, y, z}). We then numerically integrate them over time using the Dormand-Prince adaptive step solver (Boost C++ library). The single-electron currents are integrated over all momentum space using the trapezoidal rule.
While simulations of graphene and other 2D DSMs are manageable on a discretized momentum grid, the simulations quickly become intractible due to the additional dimension in 3D Dirac semimetals. When the incident laser is linearly polarized, we can exploit cylindrical symmetry by aligning the polarization axis parallel to q z = v z p z , reducing the problem to an effectively 2D one. Equations (9a) to (9b) then becomė The time derivatives of the angular terms becomeφ = 0 andθ(t) = ev z E z (t)q ρ /E(t) 2 , where E z (t) is the electric field component along z and q ρ is the radial coordinate in the scaled momentum space. Note that the above equations and the initial conditions have no angular dependence, which is also the case for the single-electron current: Since both the single-electron current and Eqs. (13a) and (13b) are rotationally-invariant about the axis of polarization, we have effectively reduced the 3D problem to a 2D one. Thus, for linearly polarized illumination, we compute the macroscopic current along z, integrated over all momentum-resolved contributions, as Modelling ultrafast THz pulses. High intensity THz pulses from compact sources are frequently realized in the few-to-single cycle regime. To simulate a pulsed plane wave that does not contain spurious non-zero DC components, we chose a Poisson power spectrum, which accurately describes ultrafast pulses with durations down to single laser period [51]: Here, ω is angular frequency, φ 0 is a phase constant, s is a real, positive parameter that determines the pulse duration, ω 0 is the peak angular frequency of the pulse, Γ is the gamma function, and Θ is the Heaviside step function. In the narrow bandwidth limit (i.e., the manycycle limit), Eq. (16) approaches a Gaussian spectrum of central angular frequency ω 0 . The electric field whose spectrum is given by Eq. (16) is where the peak amplitude is E E E 0 = (E x,0 , E y,0 , E z,0 ) and t pk is the time where the pulse reaches its peak. The intensity full width at half maximum (FWHM) τ FWHM is related to parameter s via τ FWHM = 2s 2 log(2)/(2s − 1)/ω 0 . Throughout this work, all pulse durations refer to the intensity FWHM. For τ FWHM = 2 ps, the corresponding shape factor is s ≈ 56.4. The vector potential A A A(t) = (A x (t), A y (t), A z (t)) at time t is given by THz interaction with DSM thin film To model the interaction of a linearly polarized driving THz pulse with a finite 3D DSM thin film at normal incidence, we employ a 1 + 1D finite difference time domain (FDTD) routine which fully incorporates all intraband nonlinearities induced by the incident field. At the starting location z = z in , we have the input fields E E E in = (E in,x , 0, 0) and H H H in = (0, E in,x /η 0 , 0), where η 0 is the free space impedence and E in,x is given by The wavenumber is k 0 = ω 0 /c. We define z = z in − z pk , where z pk is the initial pulse peak location. Maxwell's equations are given by The total current is J x is the free current and J P x = ∂ t P x is the current arising from the dielectric polarization P x . To evaluate J x , we use our analytical solutions for the intraband current, given by Eqs. (3) and (4).
We discretize these equations on a Yee grid, which is staggered in both time and space for the evaluation of both the electric and magnetic fields. Defining normalized time and spatial steps ∆(ω 0 t) and ∆(k 0 z) respectively, we obtain We have assumed the material is linear and homogeneous in magnetic field response: B y = µ 0 µ r H y (µ r = 1 through this work). The upper index, n, denotes the time step index. The lower index, j, denotes the spatial grid index. We also implement Mur absorbing boundary conditions [52] at both ends of our FDTD grid. Note that J x | n+1 j makes this scheme implicit since the current depends on the field E x | n+1 j at the current time step. To obtain the correct J x | n+1 j and E x | n+1 j , we employ a fixed-point interation method where J x | n+1 j = 0 is used as an initial guess to compute E x | n+1 j . We use this first pass solution to obtain a refined approximation of J x | n+1 j . This procedure is iterated until the error between two consecutively refined values of E x | n+1 j are within a specified tolerance. We varied this tolerance, the values of ∆(k 0 z), the Courant number ∆(ω 0 t)/∆(k 0 z), z pk , and simulation box width until convergence was achieved for the overall simulation. We assume free space on either side of the 3D DSM film, although this algorithm can be readily adapted to account for the presence of adjoining complex dispersive media (by incorporating the implementation in [53]), or adjoining nonlinear media of other kinds whose behavior can be captured by current J x as a function of E x like above.
We note that the use of a plane wave input pulse and a 1+1D FDTD -as opposed to a focused laser pulse and a 3+1D FDTD -ignores effects arising from beam diffraction and wavefront curvature. However, it should be noted that for a weakly focused laser pulse (even one of high field strengths), a plane wave pulse is a reasonable approximation.
Details of the implementation of the 1+1D FDTD algorithm for 2D DSMs (such as graphene) are given in SI Section VIII.
Computing HHG spectra For an x-polarized plane wave pulse E x propagating in z and therefore impinging on the DSM thin film at normal incidence, we can express the induced current density at position z within the film as a linear combination of harmonic components: J x (z, t) = Re 1 2π J x (z, ω)e iωt dω .
The energy radiated per unit solid angle Ω per unit angular frequency in the far-field is given by (derivation in SI Section IX) where A is the area of the sample. We take A = πR 2 , where R = 1 mm throughout this paper.F (θ, ω) is given by the Fourier transform of F (θ, t) = 1 2π ike iωt−ikr J 1 (kR sin θ) kR sin θ × D 0J x (z , ω)e ikz cos θ dz dω.
The wavenumber corresponding to the angular frequency component ω is k = ω/c, the first-order Bessel function is J 1 , and the thickness of the DSM thin film is D (we choose D = 250 nm throughout our work). For Fig. 1, we numerically evaluate the integral over z in Eq. (24) based on our FDTD results, followed by numerically integrating the intensity spectrum, given by Eq. (23), over all solid angles in the forward emission direction to get the energy spectral density dU/dω. To obtain the spectrum in units of photons per 1% bandwidth (BW), we divide dU/dω by a factor 100 . For Fig. 2 and 3, we present our results using Eq. (23) evaluated at θ = 0 (forward emission), and assume that the current is spatially uniform throughout the entire sample volume.
Computing energy efficiency The energy conversion efficiency is defined as the ratio of the energy of the N th harmonic, U N , to the incident pulse energy U in . The incident pulse energy is computed as where E in,x (t) is the temporal profile of the incident driving electric field and A is the area of the sample as defined in Eq. (23). To get the energy conversion efficiency of the N th harmonic, we integrate the energy spectral density dU/dω over the frequency domain from (N − 1)ω 0 to (N + 1)ω 0 .