Maximal terahertz emission in high harmonic generation from 3D Dirac semimetals

Unlike conventional semiconductor platforms, 3D Dirac semimetals (DSMs) require relatively low input laser intensities for efficient terahertz (THz) high harmonic generation (HHG), making them promising materials for developing compact THz light sources. Here, we show that 3D DSMs’ high nonlinearity opens up a regime of nonlinear optics where extreme subwavelength current density features develop within nanoscale propagation distances of the driving field. Our results reveal orders-of-magnitude enhancement in HHG intensity with thicker 3D DSM films, and show that these subwavelength features fundamentally limit HHG enhancement beyond an optimal film thickness. This decrease in HHG intensity beyond the optimal thickness constitutes an effective propagation-induced dephasing. Our findings highlight the importance of propagation dynamics in nanofilms of extreme optical nonlinearity. To increase the output intensity of conventional THz harmonic generation platforms, the medium is normally made continuously thicker. Here, 3D Dirac semimetal films are predicted to increase the THz output by orders of magnitude, but only up to an optimal thickness dictated by its nonlinear response.

T he terahertz (THz) regime has attracted much attention due to its broad range of potential applications, including electron acceleration 1-3 , imaging [4][5][6] , controlling ultrafast processes in materials [7][8][9] , and next-generation communications [10][11][12] . These emerging technologies have motivated the quest to realize compact solid-state THz high-harmonic generation (HHG) platforms. The efficiency of HHG, which involves light emission at integer multiples of the input laser frequency, favors materials with strong optical nonlinearity. Three-dimensional Dirac semimetals (3D DSMs) [13][14][15][16][17][18] , whose massless charge carriers result in extreme nonlinearity in the THz regime, is one such example. Recent experiments 19,20 and theory [21][22][23][24] have verified the promise of 3D DSMs in realizing compact, highly efficient THz HHG sources. In particular, even for moderate driving fields of ≤10 MV/m, highly efficient generation of the 3rd 19 and up to the 7th 20 harmonic have been demonstrated using the 3D DSM Cd 3 As 2 , with a theoretical study predicting efficient HHG up to the 31st order 23 . In contrast, conventional solid-state THz HHG platforms require driving fields ≥1 GV/m to generate comparable THz HHG intensities 25,26 .
Although 2D DSM graphene monolayers have also been shown in both theory [27][28][29][30][31] and experiments [32][33][34][35] to exhibit extremely high nonlinearities, its atomically thin nature results in a relatively small interaction volume, thereby limiting HHG conversion efficiencies. The finite film thickness of 3D DSMs, which results in a substantially enhanced interaction volume compared to graphene, presents itself as an additional degree-of-freedom through which the THz HHG output can be enhanced-a prospect that remains unexplored.
Here, we show that orders-of-magnitude enhancements of the output THz HHG intensity from 3D DSMs is attainable by increasing the propagation length. Specifically, our results reveal that for a 10 MV/m driving field, the output intensities of the 3rd and 31st harmonics can be enhanced by factors of 144 and 28, respectively in going from a film thickness of 50 nm to 1500 nm. We also show that highly nonlinear materials like 3D DSMs open up a regime of nonlinear optics in which an extremely subwavelength phase-flip in the induced current density appears over nanoscale propagation distances-an effect not seen in conventional nonlinear materials. This subwavelength phase-flip results in an optimal thickness for HHG in 3D DSM nanofilms, beyond which the output HHG intensity falls rapidly falls-an effect which can be understood as an effective propagation-induced dephasing mechanism. Additionally, for fixed field strengths, the optimal film thickness for all harmonics falls within a narrow range of values, indicating that many harmonics can be simultaneously optimized through a single choice of film thickness. Importantly, our result suggests that alternative strategies 36 are required to further enhance the HHG output beyond what is possible using 3D DSM films of the optimal thickness. Our findings highlight the importance of accounting for light propagation dynamics in highly nonlinear nanofilms (even beyond 3D DSMs) and pave the way to the development of efficient, solidstate THz light sources and optoelectronics based on 3D DSMs.

Results
When an x-polarized THz pulse impinges upon a 3D DSM film at normal incidence (Fig. 1a), nonlinear currents are induced in the thin film, resulting in the emission of high harmonics. In momentum space (Fig. 1a, inset), the driving laser field induces carrier oscillations within and transitions between the conduction and valence bands, giving rise to intraband and interband currents, respectively, which emit light peaked at integer multiples of the input frequency.
We model the interaction between the externally incident THz pulse and the 3D DSM film by solving Maxwell's equations using a 1 + 1 D finite-difference time-domain (FDTD) method (see "Methods" section), and determine the induced nonlinear current density J x (z, t) via nonperturbative quantum theory (Supplementary Note 1). Our approach fully accounts for reflection and transmission at interfaces, pump depletion, and all other effects arising from nonlinear propagation of the combined pump and HHG fields. To confirm that the propagation-induced effects we observe are not a consequence of temperature or scattering, we consider the low-temperature limit (T → 0 K) and no carrier scattering (τ → ∞), although our model can fully account for these effects. Under these conditions, the current density is given by: 23 where we have considered only the intraband current, which dominates when ℏω <<2Ɛ F (ℏ is the reduced Planck's constant, Ɛ F is the Fermi energy, ω is the driving angular frequency), which is the case we study. Additionally, experiments 19,20 have shown that the intraband nonlinearities are dominant for high doping values in the THz regime. Specifically, we consider an incident pulse centered at ω 0 = 2π × (1 THz) and a Fermi energy Ɛ F = 60 meV, which satisfy ℏω ≪ 2Ɛ F . In Eq. (1), g = 4 is the product of the spin and valley degeneracies, e > 0 is the elementary charge, v i is the 3D DSM's Fermi velocity along the direction i∈{x,y,z}, and a x z; t ð Þ ¼ Àe Àt=τ R t À1 e t 0 =τ E x ðz; t 0 Þdt 0 is the modified vector potential 30 . When scattering is neglected (τ→∞), which is the case for our simulations, a x (z,t) reduces to the standard definition of the vector potential: A x z; t ð Þ ¼ À R t À1 E x ðz; t 0 Þdt 0 . We compute the HHG spectra as radiation from the current distribution J x (z,t). This produces the same results as computing the HHG spectra from the reflected and transmitted electric fields (Supplementary Note 2). We consider the experimentally measured Cd 3 As 2 Fermi velocities: ðv x ; v y ; v z Þ ¼ ð1:28; 1:30; 0:33Þ 10 6 ms À113 . Our choice of Ɛ F = 60 meV is attainable through chemical doping, and the film thicknesses considered have been realized using molecular beam epitaxy 19,20,[37][38][39][40][41] . Few-cycle THz pulses with peak field strengths similar to those we consider are readily accessible using compact sources [42][43][44][45][46] . Unless otherwise specified, we use these parameters throughout our work. We emphasize that our method fully incorporates the dispersive induced refractive index through the current density given by Eq. (1). Figure 1b shows our results for a Cd 3 As 2 film driven by an incident 2 ps-long (intensity full-width-at half-maximum) 1 THz pulse of peak field 10 MV/m (in free space). In going from a film thickness of 50 nm to 1500 nm, we find that the 3rd and 31st harmonics are enhanced by factors of 144 and 28, respectively. However, beyond the optimal thickness of ≈1.5 μm, the output harmonic intensities rapidly diminish with increasing film thickness.
To investigate why the output harmonic intensities decrease beyond a certain film thickness, we plot the 3rd harmonic of the current density, defined as ReJ Fig. 1c, d for film thicknesses of 1 and 5 μm respectively. Herẽ J x z; ω ð Þ is the Fourier transform of J x (z,t). For a film thickness of 1 μm-thinner than the optimal thickness-we see from Fig. 1c that the current profile undergoes negligible phase shift with propagation distance. Hence, at any given time, the radiation emitted by the current density at different propagation distances constructively interferes. However, for film thicknesses way beyond the optimal value (Fig. 1d), the current density undergoes a π-phase shift across the film. Consequently, the radiation emitted by currents at different z destructively interferes, resulting in reduced output intensity. It is noteworthy that this phase-flip takes place at an extremely subwavelength propagation distance of 1 μm -100 times smaller than the driving wavelength. Furthermore, this phase flip also occurs within the skin depth of the 3rd harmonic (about 1.4 μm for scattering times τ∼150 fs at T = 0 K), which emphasizes the importance of considering these extreme nonlinear propagation dynamics, even for films that are thin enough such that losses from attenuation do not dominate.
While we consider the τ→∞ and T→0 K limits here, we show in Supplementary Note 3 that our qualitative conclusions still hold when the effects of finite T and τ are included. In particular, we find that for faster τ down to tens of femtoseconds, the change in the phase of the 3rd harmonic current becomes more gradual as a function of z, and the difference in phase between current components on opposite ends of the film deviates further from π. While this leads to partial cancellation between radiation emitted by current components on opposite ends of the film, the reduction of output HHG intensity is still significant since the 3rd harmonic current still undergoes a net polarity reversal across the film. Additionally, we expect finite temperatures (up to T = 300 K) to have a less dominant effect on the nonlinear propagation dynamics relative to finite scattering times that fall within the range of experimentally measured values.
We find that this extreme subwavelength phase-flip only arises within a regime of nonlinear optics where highly nonlinear nanofilms, like 3D DSMs, are involved. We show this in Fig. 2, where the phase shift of the 3rd harmonic current density in a 2500 nm-thick, non-dispersive dielectric film driven by a peak field of 10 MV/m is plotted as a function of its linear and 3rdorder susceptibilities, χ (1) and χ (3) respectively. We see that a ≈ π-phase flip over subwavelength propagation distances only manifests for materials with χ ð3Þ ≥ 10 À12 m 2 =V 2 , such as 3D DSMs. In contrast, conventional materials, which possess relatively small 3rd-order susceptibility by comparison χ ð3Þ ≤ 10 À16 m 2 =V 2 À Á , see negligible phase shifts over subwavelength thicknesses. See "Supplementary Note 4" for details of the simulations used to plot Fig. 2. Figure 3 shows the output HHG intensity as a function of film thickness for various harmonics. While the optimal thickness generally differs across harmonics, they lie within a relatively narrow range of values, implying that the output intensity of all harmonics can be simultaneously optimized with a single choice of film thickness. Considering a peak field of 2 MV/m (Fig. 3a), the maximum output intensity for each harmonic is achieved with film thicknesses ranging from 150-300 nm. Considering a peak field of 10 MV/m (Fig. 3b), the optimal film thickness lies between 600 and 1600 nm. The blue shading in Fig. 3a, b indicate these ranges of optimal thicknesses.
We also observe that an increase in optimal film thickness generally accompanies stronger driving fields-a finding that seems unintuitive since larger fields are associated with stronger nonlinearity, which should in turn lead to a smaller optimal film thickness. However, it is important to note that stronger driving fields induce larger harmonic currents, which in turn emit more a An x-polarized laser pulse centered at angular frequency 1 THz impinges on a 3D DSM thin film at normal incidence, resulting in the radiation of higher harmonics. In momentum space p ¼ ðp ? ; p x Þ (inset), the driving field induces intraband carrier oscillations (solid purple arrows) within and interband carrier transitions (dotted purple arrow) between the valence and conduction bands of the Dirac cone, resulting in the emission of light peaked at odd-integer multiples of the input angular frequency ω 0 . b shows orders-of-magnitude enhancements in output intensity with increasing film thickness up to an optimal value of about 1500 nm, beyond which the harmonic output drastically falls. The existence of an optimal thickness arises due to the propagation-induced phase shift of the harmonic current (normalized 3rd harmonic current plotted in c, d) for 1 μm and 5 μm) across the film thickness. In thinner films c the phase shift in the current as a function of propagation distance z is insignificant. Thus, the emitted waves from different z are in-phase and add constructively. For films much thicker than the optimal thickness d a π-phaseflip occurs in the current density. The emitted radiation from opposite sides of this phase flip destructively interfere, resulting in the drastic decrease in HHG output with increasing film thickness seen in (b).
intense HHG. The more intense radiation emitted by the harmonic current components at each z, coupled with a phase flip in the harmonic current that occurs at larger z for stronger fields ( Supplementary Fig. S5 in Supplementary Note 5), implies that a larger propagation distance is required to significantly suppress output HHG intensity. This explains why a larger field strength results in a larger optimal film thickness.
Our investigations also reveal that the decrease in output HHG with increasing film thickness can be understood as a propagationinduced dephasing mechanism. We verify this by computing an effective scattering time τ eff in the absence of propagation. In Fig. 4, we show the excellent agreement between the exact spectrum obtained from Maxwell's equations and the spectrum obtained from Eq. (1) fitted to an effective scattering time τ = τ eff while neglecting propagation. Considering an incident pulse with 2 MV/m peak field (Fig. 4a-c), we find that the output spectra for film thicknesses of 800 nm, 1 μm, and 2.5 μm are effectively captured by τ eff = 16.5 fs, 12.8 fs, and 2.7 fs, respectively. As expected, thicker films (beyond the optimal thickness) correspond to faster scattering times-a trend that holds for a peak field of 10 MV/m (Fig. 4d-f): the spectra for film thicknesses 2 μm, 2.5 μm, and 5 μm are captured by τ eff = 7.5 fs, 5.1 fs, and 0.2 fs, respectively. Note that these values of τ eff are much shorter than the typical scattering times of τ ≈ 150 fs 19 in cleaner samples of Cd 3 As 2 , indicating that propagation-induced dephasing is the dominant dephasing mechanism for thicker films. Neglecting propagation effects-to show that the decrease in HHG output can be captured by τ eff -involves assuming a driving field throughout the film that has the same field profile as the irradiated surface. Our procedure finds τ eff such that the area under the spectrum from 2ω 0 to 32ω 0 is the same as the area under the exact propagated numerical spectrum over the same frequency range.

Discussion
Recent experiments 19,20 and theory 23 have shown the bulk nature of 3D DSMs enables the output harmonic intensity to exceed the Fig. 2 Extreme subwavelength phase-flip of the induced current in highly nonlinear materials. The colormap shows the phase shift in the 3rd harmonic current for a 2500 nm-thick, non-dispersive dielectric film as a function of its linear and third-order susceptibilities, χ (1) and χ (3) respectively, driven by an external field of 10 MV/m. We find that the extreme subwavelength phase-flip of the induced current occurs within a regime of nonlinear optics accessible only by highly nonlinear materials like 3D Dirac semimetals (DSMs); no such phase-flip occurs at extreme subwavelength propagation distances in conventional nonlinear materials, which have weaker nonlinearity (χ (1) and χ (3) values 48 of conventional materials plotted as white crosses). We use the same THz pulse parameters as Fig. 1.   Fig. 3 Optimizing film thickness for HHG in 3D Dirac semimetals (DSMs). a For a 3D DSM film driven by a pulse of peak amplitude 2 MV/m and central angular frequency ω 0 = 2π × 1 THz, the optimal thickness for all harmonics falls between 150 and 300 nm. For a larger driving field of 10 MV/m b the optimal thickness of all harmonics lies between 600 and 1600 nm. The optimal thicknesses for all harmonics from the 3rd to the 31st lie within the blue shaded areas. The optimal thickness generally shifts to larger values when driven by stronger fields. In both panels, the Nth curve from the top has been translated downwards by 10 4(N − 1) for visual clarity; the color of the data points represent their actual magnitudes in units of photons/1% bandwidth (BW). We consider the same parameters as Fig. 1.
THz HHG output 32 from monolayer graphene by over ten times. Importantly, we find that increasing the film thickness beyond previously considered values 19,20,23 is indeed a promising strategy to further enhance the HHG output intensity by orders of magnitude. However, our studies also reveal that enhancing HHG output intensity through increased film thickness is only possible up to an optimal thickness. We find that using nanofilms of highly nonlinear materials, of which 3D DSMs are only one example, places us within a regime of nonlinear optics in which the induced current undergoes as much as a π-phase-shift over extreme subwavelength propagation distances. This phase-flip results in the existence of an optimal thickness, beyond which the HHG output intensity falls. We show that this decrease in HHG intensity can be captured by an effective propagation-induced dephasing time τ eff . It is noteworthy that this propagationinduced dephasing occurs even over extremely subwavelength propagation lengths of ∼100 nm to ∼1 μm-about 100-1000 times shorter than the central driving wavelength. In contrast, a recent theoretical study 47 has shown that propagation-induced dephasing in conventional materials, which have weaker nonlinearity, only becomes significant when the film thickness far exceeds the central driving wavelength. We emphasize that their result is entirely consistent with our findings presented in this work: in conventional nonlinear materials, negligible phasechange of the induced current occurs over extreme subwavelength propagation distances, which in turn results in negligible propagation-induced dephasing on subwavelength scales. These observations are in line with conventional nonlinear optical theory, which tells us that nonlinear propagation effects become significant over a characteristic length scale that falls inversely with the nonlinear susceptibility: 48 L / 1=ðω 0 χ ð3Þ Þ. Additionally, as a result of a large χ (3) , frequency up-conversion is extremely efficient, leading to strong high-harmonic fields, which themselves accumulate significant nonlinear effects over even shorter propagation lengths (since L / 1=ω 0 ). This explains why propagation-induced dephasing in conventional nonlinear materials is only expected for film thicknesses far exceeding the central driving wavelength. Our study thus establishes the important role of light propagation effects on extreme subwavelength distances in this regime of nonlinear optics where highly nonlinear materials-such as emerging topological semimetals beyond 3D DSMs 49-58 , for instance-are considered.
Importantly, our results suggest that restrictions on the optimal film thickness may be circumvented by using various metamaterial configurations. For instance, a recent experiment 36 demonstrated orders-of-magnitude enhancement in the THz harmonic output from a graphene-photonic grating metamaterial relative to the output from a graphene monolayer-a result enabled by the field enhancement from the grating. Platforms that combine 3D DSM films with similar photonic structures could be a promising method allowing far greater enhancements of THz HHG output beyond what is fundamentally permitted by a single 3D DSM layer. An alternative metamaterial configuration would be a superlattice of 3D DSM nanofilms (e.g., 1D array of 3D DSM films) interleaved with other materials (e.g., metals or dielectrics). This could be useful since 3D DSMs nanostructured in this manner could also be used to enhance the output HHG intensity even beyond the optimal film thickness in a single-layer unstructured 3D DSM film. These 3D DSM superlattices would also introduce additional versatility to the output HHG spectrum, for instance, to selectively amplify or attenuate specific harmonics, enabling on-demand THz light-shaping on a chip-scale platform.

Conclusions
In summary, we have shown that the output intensity of THz HHG in 3D DSMs can be enhanced by orders of magnitude with increased propagation length. Additionally, we find that using nanofilms made of highly nonlinear materials like 3D DSMs opens up a regime of nonlinear optics where extremely subwavelength features develop in the induced current density, an effect that is negligible in conventional nonlinear materials. The phase-flip in current density results in the destructive interference of emitted radiation from different parts of the film. We verify that this phenomenon constitutes propagation-induced dephasing of the HHG process by computing an effective dephasing time τ eff for light generation in 3D DSMs. We also show that the optimal film thickness for all output harmonics lies within a relatively narrow range for a given driving field, and that this optimal film thickness generally increases for stronger fields. Our findings highlight the importance of including pulse propagation effects for highly nonlinear materials, including 3D DSMs and other bulk topological semimetals. Furthermore, our findings For all panels, we plot the exact output HHG spectrum computed via Maxwell's equations as solid red curves, and the output HHG spectrum computed via Eq. (1) fitted to an effective scattering time τ eff -while neglecting propagation effects-as blue dashed-dotted curves. a, b, c show the spectra when 3D DSM films of thicknesses 0.8 μm, 1 μm, and 2.5 μm, respectively are driven by a peak field of 2 MV/m. d, e, f show the spectra for film thicknesses of 2 μm, 2.5 μm, and 5 μm respectively, driven by a peak field of 10 MV/m. Our fitted values of τ eff are in good agreement with the exact spectra, implying that the reduced HHG can be captured by an effective dephasing time and is thus due to propagationinduced dephasing effects. As expected, larger film thicknesses correspond to shorter τ eff . We consider the same parameters as in Figs. 1 and 3. Here, we normalize all angular frequencies of the output ω against the fundamental driving frequency ω 0 .
suggest that appropriate nanostructuring could enable stillgreater THz HHG output efficiency to be achieved. Our work paves the way to efficient solid-state THz light sources and optoelectronics based on highly nonlinear material platforms like 3D DSMs.

Methods
Numerical simulation of 3D DSM thin-film interaction with a THz pulse. To model the interaction of a driving THz pulse linearly polarized along x with a finite 3D DSM thin film at normal incidence, we employ a (1 + 1)D finite-difference time-domain (FDTD) routine that fully incorporates all nonlinear propagation effects of the combined pump and higher-harmonic field, as well as reflection and transmission at interfaces. At the simulation box edge z = z in , we have the input fields E in ¼ ðE in;x ; 0; 0Þ and H in ¼ ð0; E in;x =η 0 ; 0Þ, where η 0 is the free space impedance and E in,x is given by a plane wave with a Poisson power spectrum: 59 where k 0 ¼ ω 0 =c is the free-space light wavenumber, c is the speed of light in free space, ϕ 0 is a phase constant, and s > 0 is a real parameter that determines the pulse duration (we use s = 56.4 throughout our work, which corresponds to an intensity full-width at half-maximum (FWHM) pulse duration of ≈2 ps). We define z 0 ¼ z in À z pk , where z pk is the initial pulse peak location. We note that a plane wave input pulse and a (1 + 1)D FDTD neglects effects like beam diffraction and wavefront curvature (as opposed to a focused laser pulse and a (3 + 1)D FDTD). However, for a weakly focused laser pulse (even one of the high field strengths), a plane wave pulse is a reasonable approximation. We proceed to write Maxwell's equations in the form The total current is the sum J x is the free current and J P x ¼ ∂ t P x is the current arising from the dielectric polarization P x . Discretizing these equations on a Yee grid and defining normalized time and spatial steps Δ(ω 0 t) and Δ(k 0 z) respectively, we obtain We assume that the material is linear and homogeneous in magnetic field response: B y ¼ μ 0 μ r H y (μ r = 1 throughout our work). The upper index n denotes the time step. The lower index j refers to the spatial grid position. Note that J x nþ1 j makes this scheme implicit because the current depends on the field E x nþ1 j at the current time step. To obtain the self-consistent J x nþ1 j and E x nþ1 j , we adopt a fixedpoint interation method in which J x nþ1 j ¼ J x n j (i.e., the previous time step value) 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 is within a specified tolerance. We assume free space on either side of the 3D DSM film and implement Mur absorbing boundary conditions 60 at both ends of our FDTD grid. In the limit _ω 0 ( 2μ (where μ is the chemical potential) and finite temperatures, J x ðz; tÞ is computed using Eq. (S14) in Supplementary Note 1. When T→0 K, we compute J x z; t ð Þ using main text Eq. (1). The HHG spectrum in terms of the energy spectral density can then be computed from the numerical solution of J x z; t ð Þ using: where d 2 U/dωdΩ, which denotes the energy U radiated into a unit solid angle Ω and unit angular frequency ω in the far-field, takes the following form: 23 d 2 U dΩdω ¼ A 2 8π 3 ϵ 0 c cos 2 ϕ cos 2 θ þ sin 2 ϕ À Á ike Àikr J 1 ðkR sin θÞ kR sin θ Z D

0J
x z 0 ; ω ð Þe ikz 0 cos θ dz 0 Here, A = πR 2 is the area of the circular radiating region of radius R, θ is the polar angle (with respect to forward propagation directions), ϕ is the azimuthal angle, k = ω/c is the wavenumber, J 1 is the first-order Bessel function, r is the radial position to the observer, D is the film thickness, andJ x z 0 ; ω ð Þis the Fourier transform of J x (z′, t). The integral is taken with respect to the source coordinate z′. Note that the current J x is assumed to be uniform in the transverse (x,y) directions. We consider R = 1 mm for all cases presented in our work.

Data availability
The simulation results presented in this paper are available from the corresponding author on reasonable request.

Code availability
The custom computer code used to generate the results is the 1+1 D finite-difference time-domain scheme described in the "Methods" section.