All-dielectric nanostructured materials have been extensively applied to enhance and confine electric and magnetic fields at subwavelength scales for nonlinear optical applications [1–4]. Several strategies were proposed to increase the second- and third-order nonlinear conversion efficiency in semiconductors, based on Mie resonances in high-refractive-index nanostructures [5,6], as well as resonant metasurfaces [7–9]. Recently, a new paradigm of harmonic generation beyond the perturbative regime by strong field excitation has been demonstrated in semiconductors, opening new frontiers for previously non-accessible enhanced high-order harmonic spectra, consisting of odd and even harmonics and stretching far beyond third-order nonlinearities [10–15], and ultrafast all-optical modulation [16–19]. Comprehensive understanding of these coherent nonlinear excitation processes in semiconductors is still lacking, and self-consistent models are highly challenging due to the inclusion of the dynamical material response on a quantum level beyond the perturbative regime. Nonlinear conversion efficiency depends strongly on multiple involved phenomena such as intra-band and inter-band excitation dynamics [10,11,17], absorption induced by free carriers limiting the harmonics yield above the material bandgap [7] and shifting the nanostructure resonance positions [4,20], and nonlinear propagation effects [21,22].

The nature of odd and even harmonic generation in semiconductors has been established only recently, underlying the pronounced role of quantum interference and indirect transitions between valence and conduction bands [11,23,24]. The latter effects can be described by a quantum approach, as provided by semiconductor Bloch equations (SBEs) under inclusion of more than two mutually dipole coupled bands [25]. This momentum-resolved approach evaluates microscopic coherent carrier densities and polarizations within the non-equilibrium electronic system, considering the full Brillouin zone (BZ). Typically, a realistic model for the electronic band structure for a particular semiconductor is provided by density functional theory (DFT) simulations. In contrast to perturbative models, SBEs automatically include nonlinearities of any order as well as the full dynamic evolution of electronic material excitations. Furthermore, the non-perturbative nonlinearities were shown to play a significant role in off-resonant strong-field excitation and high-harmonic generation in semiconductors, attributed to dynamic intra-band oscillations that can be captured only by employing a dynamic microscopic model such as SBEs [26,27]. These non-perturbative scaling laws for high harmonics are influential in gallium arsenide (GaAs) for laser intensities of at least 1GW/cm2 [21,26] and can lead to an increase in the high-harmonic yield.

To simulate the nonlinear spectrum evolution upon ultrashort laser pulse propagation through bulk material, the SBEs were further coupled to unidirectional or 1D/2D propagation codes by macroscopic polarization currents [22,28,29]. To date, only a few works have focused on the self-consistent coupling between full-vector Maxwell equations and SBEs [30–32], albeit not for harmonic generation. A 3D modeling is essential here not only to reproduce the exact geometry and resonant nature of the nanostructures, but also to take account the inhomogeneous distribution of laser-induced carriers [17,18,33,34].

In the current work, we demonstrate the feasibility and robustness of a coupled Maxwell–Bloch approach for investigating the high-order harmonic spectra for finite slab and 3D sub-wavelength nanoparticle geometries in the non-perturbative regime of ultrashort laser pulse excitation. Driven by either a Fabry–Perot-like standing wave resonator or magnetic dipole resonances, inhomogeneously photo-induced carriers are localized in the semiconductor material, producing enhanced even and odd harmonics.

The approach consists of two coupled parts: Maxwell equations with nonlinear polarization P→ and current J→ sources responsible for inter-band and intra-band excitations, respectively, and dynamic SBEs, here utilizing a three-band structure model of a direct-gap non-centrosymmetric semiconductor using material parameters representative of GaAs. The SBEs are sourced by the electric field E→ providing the non-perturbative material responses in the form of macroscopic polarization and current for Maxwell input. The SBEs are solved in each spatial point of the nonlinear material.

The propagation model is based on full-vector Maxwell equations with macroscopic currents to include linear dispersion and inter-band and intra-band excitation upon ultrashort pulse excitation as follows: {∂E→∂t=∇×H→ϵ0ϵb−∑iJ→dispi+J→+∂P→∂tϵ0ϵb∂H→∂t=−∇×E→μ0,(1)where E→ and H→ are electric and magnetic fields, and ϵ0 and μ0 are the permittivity and permeability of free space, respectively; ϵb is the background dielectric constant. The dispersion relation for our GaAs model is applied from Ref. [35] through the current J→dispi in Eq. (1). J→ and ∂P→∂t correspond to inter- and intra-band currents, respectively.

Throughout the paper, a linearly polarized plane wave excitation source is used. The electric field source is chosen to be an x-polarized Gaussian pulse with a full width at half maximum of θ=100fs. Light propagates along the z direction. Maxwell equations are solved by a finite-difference time-domain (FDTD) approach with convolutional perfectly matched layers (CPMLs) at the boundaries [36]. For the 1D problem, a spatial step of 10 nm and temporal step of ≈2 as were chosen. For the 3D problem, the corresponding resolution was 5 nm spatially and ≈4 as temporally. Generally, temporal steps smaller than the ones defined by the Courant–Friedrichs–Lewy condition were required to resolve high-order harmonics.

Our GaAs band structure parameters are obtained via DFT assuming (111) crystal orientation. In our simulations, we include the mutually dipole coupled conduction e, split-off valence band h1, and light hole valence band h2. The dipole coupling strength dll′ for the corresponding electronic transitions ll′ and energy dispersion of each band εl, shown in Figs. 1(a) and 1(b), are derived as a fit to the results of our DFT calculation using the Vienna ab initio simulation package [37,38]. The accurate reproduction of the bandgap is realized by employing the Heyd–Scuseria–Ernzerhof hybrid functional [39] customized for GaAs. To keep our evaluations as simple as possible, we neglect the dipole phases and Berry phases, which is an acceptable approximation for the chosen excitation conditions [24]. We use a momentum grid size [−kmax;kmax] with 1200 points, where kmax=3π/a, and a=5.663Å is the lattice constant. Our consideration of a single path through the BZ along Γ-L direction is justified by the use of a linearly polarized source along this direction. To avoid nonphysical reflections at the boundaries of the BZ, a third-order polynomial fit is applied to the energy boundary points to conserve both the energy values and their derivatives. The resulting energy gradients ∇kϵl are shown in Fig. 1(c), smoothly decaying to zero at the edges of the grid.

The macroscopic intra-band current and inter-band polarization are calculated through complex microscopic polarizations pkll′ and carrier distributions fkl as follows: {Px=∑k,l,l′dkll′Re(pkll′),Jx=eℏ∑k,l(∇kεkl)fkl,(2)where l,l′ stand for the electron e or hole h1,h2 bands, and the summation is performed for all momentum points k. SBEs with sourced linearly x-polarized electric fields [25] are solved to define the microscopic polarizations and carrier densities as follows: {iℏ∂∂tpkhie=[εke+εkhi+ieEx(t)∇k−iℏ/T2]pkhie+[dkhjhipkhje−dkehjpkhihj−dkehi(1−fke−fkhi)]Ex(t),iℏ∂∂tpkh1h2=[εkh1−εkh2+ieEx(t)∇k−iℏ/T2]pkh1h2+[dkeh1pkeh2−dkh2epkh1e+dkh2h1(fkh1−fkh2)]Ex(t),ℏ∂∂tfke=−2Ex(t)∑i[Im(dkehipkehi)]+eEx(t)∇kfke,ℏ∂∂tfkhi=−2Ex(t)[Im(dkhjhipkhjhi)+Im(dkehipkehi)]+eEx(t)∇kfkhi,(3)where i=1,2, j≠i, T2=20fs is the dephasing time, and Re(⋅) and Im(⋅) stand for the real and the imaginary parts of the microscopic polarization, respectively. Since the Coulomb interaction was shown to play only a minor role for off-resonant strong field excitation, it is sufficient for the present purposes to include it only via the dephasing term [11]. The microscopic carriers in the conduction band are summed up and averaged through N=1200 momentum points as follows: Ne=∑kNfke/N.

To solve the SBEs, a fourth-order Runge–Kutta method is applied for the ensemble of microscopic polarizations pkll′ and carrier distributions fkl. Macroscopic quantities Px and Jx are then evaluated according to Eq. (2) after each four-step procedure. A fully symmetric fourth-order finite-difference approximation is implemented for microscopic gradient values ∇kpkll′ and ∇kfkl. Periodic boundary conditions are then applied for the microscopic polarizations and distributions at the borders of the BZ. For an uncoupled SBE model, the electric field is sourced directly to the SBE, and the Fourier transforms of the local emission spectra |iJx+Pxω|2 are plotted for the laser wavelengths of interest in Fig. 1(d). The results with N=1200 momentum points have been tested for convergence, giving identical results for the larger number of points N=2400. The local emission spectra indicate a complex behavior for high-order odd and even harmonics. The odd harmonics are mainly the consequence of direct transitions from the valence to conduction band, which can be illustrated by a perturbative expansion of the contributions to the material responses from SBEs. Equation (3) for microscopic polarizations gives rise to odd-order terms that involve only the coupling between two bands dkhie(1−fke−fkhi)Ex∝dhieE+dhie3E3+…, whereas the even order contributions appear in mixed terms of the expansion (dkhihjpkhje−dkhjepkhihj)Ex∝dhihjdhjeE2+(dhihjdhje)2E4+…, which show a coupling of three bands and can be attributed to indirect transition pathways. As a result, at least a three-band structure is required to fully describe both even and odd harmonics in our model.

In the framework of the SBE model, the major changes in both real and imaginary refractive indices are considered via the complex intra-band and inter-band currents. Band-filling dynamics is included via the Pauli blocking factor (1−fke−fkh), which reduces the oscillation strength as the plasma density increases. Complementary contributions might be expected by the bandgap shrinkage due to carrier–carrier scattering in Coulomb interactions, which are not included in the current model and could be further incorporated in the semiconductor Bloch formalism as proposed in Ref. [25]. It was shown in that reference that for strongly off-resonant excitation, Coulomb effects do not lead to additional features, but influence only spectral details and do so only in materials with exceptionally high Coulomb interaction such as transition metal dichalcogenide monolayers (TMDCs). For the GaAs system investigated here, Coulomb-induced renormalizations are almost two orders of magnitude smaller than in TMDCs and should be of minor importance. Including Coulomb effects would dramatically increase the numerical effort. Therefore, we include only the Coulomb scattering effects through the polarization dephasing term but neglect the bandgap and field renormalization.

Our numerical method for the coupled Maxwell Eq. (1) and SBE Eq. (3) evaluates additionally the nonlinear electric field component Ex inside a fourth-order Runge–Kutta procedure applied for the microscopic polarization and carrier densities. The macroscopic quantities at each spatial position affect the carrier dynamics and nonlinear propagation at each subsequent temporal step. Our computing algorithm uses spatially parallel implementation based on graphics processing units (GPUs).

The field is sourced in vacuum at a distance from the nonlinear material. The time history of the electric field evolution is recorded in space positions before and after hitting the semiconductor material to obtain separate reflected and transmitted field signals, which are then further analyzed by Fourier transforms. A harmonic and decaying window function is applied at the boundaries of the temporal domain [35] to reduce the noise floor due to the finite size of the used Gaussian pulse and the accumulation effects related to the numerical dispersion.

For our study, there are two specific regimes of interest for harmonic generation in semiconductor nanomaterials: a propagative regime in finite structures on the distances of order or larger than laser wavelength and a field confinement regime in subwavelength nanostructures, where the field can be significantly enhanced and the photo-generated carriers are localized inside the nanostructure.

First, we investigate the ultrashort laser pulse propagation through finite slabs of different thicknesses, effectively representing examples of a Fabry–Perot resonator due to the refractive index differences between vacuum and semiconductor material. Light waves entering the slab, being smaller than the optical pulse length d<cθ≈30μm, undergo multiple partial reflections, constructively interfering at the resonant frequency. Apart from creating multiple transmitted sub-pulses, nonlinear propagation is accompanied by the generation of electron densities and absorption losses upon propagation. For our study, an ultrashort pulse with a fundamental λ=10.6μm wavelength is chosen exhibiting a pronounced interplay of direct and indirect inter-band transitions contributing to high-order harmonic generation.

The spatiotemporal evolution of the electric fields for slab thicknesses of 12 and 24 μm is shown in Figs. 2(a) and 2(c) with the laser pulse propagating from down upwards. We see that for the thicker slab, the first and secondary reflection signals are separated in time, whereas for the thinner slab, they overlap. The maximum electron density distributions established at the end of the pulse are detailed in Figs. 2(b) and 2(d), representing periodic concentrated regions of increased charge carriers spaced by λ/n≈3.1μm, created by standing waves as a result of an interference between the transmitted and reflected waves from the slab boundaries (t>1ps). Strong spatial density gradients are formed on the back side of the slab, potentially contributing to an enhanced harmonic yield.

In Fig. 3(a), we show spectra of the transmitted fields as a function of the propagation distance in comparison to the local emission spectrum |iJx+Pxω|2 obtained by the uncoupled SBE model. Here, we refer to the uncoupled point SBE model by indicating that the electric field is sourced by a plane wave source inside the material without solving Maxwell equations, i.e., excluding propagation effects. We note that the harmonic yield tends to increase for a larger slab thickness, keeping approximately a similar shape to the spectrum of the point model. In addition, the appearance of the sub-harmonics is also apparent between sixth and seventh harmonics attributed to frequency mixing in nonlinear propagation, as well non-trivial behavior for higher-order harmonics with increasing propagation distance.

Reflection and transmission spectra are compared for 12 and 48 μm slabs in Figs. 3(b)–3(e) [Fourier transforms in Figs. 3(b) and 3(c), and corresponding temporal signals in Figs. 3(d) and 3(e)].

Nonlinear field effects can be significantly enhanced in subwavelength nanostructures of high-refractive-index materials supporting Mie resonances. We investigate here the ultrashort laser pulse excitation of a single spherical GaAs-like nanoparticle at λ=1.5μm near-infrared wavelength (refractive index n=3.3817, transparency range). In this excitation regime, the harmonics starting from the second lie above the material bandgap (Eg=1.42eV). Nevertheless, as the considered subwavelength radii in the range 150–240 nm are much smaller than the attenuation length of the photo-induced free carriers, the generated harmonics can be efficiently extracted, even if lying within the absorption region.

Extinction cross sections as a function of the nanoparticle radii at 1.5 μm wavelength are calculated by Mie theory [40] and plotted in Fig. 4(a). The sizes of choice are marked by colored points, corresponding to different positions respective to a magnetic dipole resonance (“MD” in the figure) and contrasting electric field distributions inside the nanoparticles. The intensities corresponding to an incident field E=108V/m are plotted in Figs. 4(c)–4(f) for the sizes of 150, 180, 210, and 240 nm, respectively. For the smallest nanoparticle within the Rayleigh range, the fields inside the material are weak and almost homogeneous. In contrast, larger sizes exhibit toroidal field distribution with the strongest field enhancement for radius 210 nm close to the magnetic dipole resonance. With increasing size, the intensity decreases following the extinction cross section. The resulting 1D spatial carrier distributions along the propagation direction (left to right) are compared for different radii in Fig. 4(b). The highest values and inhomogeneity of charge densities are obtained for the nanoparticle of 210 nm radius.

Figures 5(a) and 5(b) show the computed transmission and reflection spectra for nanoparticles of different radii, indicating the excitation of second-, third-, and fourth-order harmonics. The field spectra are collected at (0, 0, 1) and (0,0,−1)μm, respectively, at 1 μm distance from the nanoparticle center along propagation axis z. The harmonic yield for the nanoparticle R=210nm close to the magnetic dipole resonance exceeds by several orders of magnitude the results for off-resonant particles of smaller sizes R=150nm and R=180nm, as well as at least by two orders the nonlinear efficiency of third- and fourth-order harmonics for a larger nanoparticle of R=240nm, close to an electric dipole resonance.

In conclusion, we have introduced a microscopic Maxwell–SBE model that allows us to fully analyze high-harmonic generation in semiconductor nanostructures. This approach can be fully customized to accurately model a specific semiconductor material and explicitly consider different crystal orientations and direction-dependent effects. Time-dependent carrier dynamics and individual band contributions are fully accessible, allowing us to pinpoint the influences of band structure features. Additionally, phase-dependent effects could be incorporated, including but not limited to the influence of Berry phases for topologically non-trivial materials. The model is also flexible to include different features of the driving field in the analysis, e.g., carrier–envelope phase or even a combination of resonant and off-resonant pulses. This specificity allows for a highly improved quantitative comparison with experimental measurements implying different polarization states without relying on experimental or phenomenological input. The coupled SBEs and Maxwell equations approach can be applied to complex geometries and metasurfaces, benefiting from the periodic boundary conditions for electric and magnetic fields in the Maxwell solver.

We apply our model to the regimes of ultrashort pulse propagation in a finite slab, supporting Fabry–Perot resonances, and in a subwavelength nanoparticle, supporting magnetic dipole Mie resonance. Both regimes are characterized by inhomogeneous carrier distribution inside the structures and generation of even and odd harmonics of a higher yield for resonant conditions. Overall, implementation of our self-consistent approach to arbitrary geometries has a strong potential for optimal design of semiconductor optical metadevices, with emerging applications for high-resolution imaging and attosecond science.

Acknowledgment. The authors are grateful to Prof. Miroslav Kolesik for fruitful discussions.