Crystal-momentum-resolved contributions to multiple plateaus of high-order harmonic generation from bandgap materials

We study the crystal-momentum-resolved contributions to the high-order harmonic generation (HHG) in bandgap materials, and identify the relevant initial crystal momenta for the first and higher plateaus of the HHG spectra. We do so by using a time-dependent density-functional theory model of one-dimensional linear chains. We introduce a self-consistent periodic treatment for the infinitely extended limit of the linear chain model, which provides a convenient way to simulate and discuss the HHG from a perfect crystal beyond the single-active-electron approximation. The multi-plateau spectral feature is elucidated by a semiclassical k-space trajectory analysis with multiple conduction bands taken into account. In the considered laser-interaction regime, the multiple plateaus beyond the first cutoff are found to stem mainly from electrons with initial crystal momenta away from the Gamma point (k = 0), while electrons with initial crystal momenta located around the Gamma point are responsible for the harmonics in the first plateau. We also show that similar findings can be obtained from calculations using a sufficiently large finite model, which proves to mimic the corresponding infinite periodic limit in terms of the band structures and the HHG spectra.

The solid-state HHG manifests interesting features that are quite different from the atomic and molecular cases. For example, the high-energy harmonic cutoff in a bulk crystal was found to scale quasi-linearly with the driving field strength rather than quadratically as in gasphase atoms and molecules [1]. According to the present understanding of HHG from the class of solids, which can be adequately described by band theory (see, e.g., reviews [11,12]), two types of transition mechanisms are involved: (i) interband transitions, i.e., the electron transitions between different bands at specific crystal momentum k, and (ii) intraband transitions which involve kspace electron motions within specific bands. The latter type of transitions is unique to solids and determined by the shape of band structures. According to Refs. [21,22], the intraband and interband HHG exhibit a fundamentally different wavelength dependence. The interplay between intraband and interband mechanisms can possibly be used to control the HHG in solids [23]. Capturing the idea of the simple-man model for atomic HHG [24], a generalized three-step model [11,21] was proposed for the interband HHG in solids. It involves the following steps: (i) near the peak of the laser field, an electron tunnels from the valence band to the conduction band, producing a hole in the valence band; (ii) the electron and hole move in the corresponding bands, driven by the external field; and (iii) they recombine and emit HHG radiation with a frequency corresponding to the band energy difference at the crystal momentum at which the recombination occurs.
In Ref. [49], a finite-system 1D TDDFT model was in-arXiv:2007.08931v1 [physics.atom-ph] 17 Jul 2020 troduced to establish the essentials of HHG by linearly polarized laser pulses in solids. With this approach, one can go beyond the single-electron models and take the electron-electron interactions into account in the TDDFT sense. The importance of accounting for the contributions from all the orbitals in a fully-occupied band, i.e., from multiple k points, was pointed out in Ref. [49]. The response of a range of system sizes, from atomic to bulk solids, was studied with this model to reveal how the harmonic cutoff changes from atomic to solid-state HHG [50]. Due to its flexibility in constructing self-consistent model systems, this methodology was applied to investigate topological edge effects [51,52] and various types of imperfection effects such as doping [53], disorder [52,54] and vacancies [55]. We note that recently such a finitesystem model was also used for studying the carrierenvelope-phase effects [60], for examining the phase invariance of the SBEs [61] and for comparing semiclassical trajectory models [62]. For investigating HHG from an ideal periodic crystal lattice, however, it could be advantageous to extend the finite-system TDDFT model to the infinite periodic limit, as we will do in the present work.
Similar to previous studies [49,50,53,55], here we employ the 1D TDDFT model to study HHG from a linear chain of atoms exposed to linearly polarized laser pulses. In this work, we extend the previous finite-chain model to the corresponding infinitely-extended limit with a periodic treatment, and demonstrate that a sufficiently large finite model indeed mimics the infinite periodic limit in terms of band structures and HHG spectra. In terms of physical findings, we focus on the multi-plateau structure in the HHG spectra, and explain this feature with a semiclassical trajectory analysis accounting for the contributions of many independent electrons, i.e., from many different k points. The infinite periodic extension of the TDDFT model offers an efficient way to simulate HHG in a perfect crystal, and makes the crystal momentum k a good quantum number. We therefore investigate the contributions from electrons with different k values in the HHG processes, which further support the many-electron interpretation of the multiple plateaus and provide an opportunity to identify the relevant k regions for different parts of the HHG spectra. We find that electrons with initial crystal momenta located close to the Γ point are responsible for the first plateau, while the multiple plateaus beyond the first cutoff stem from electrons with initial crystal momenta away from the Γ point. In addition, we show that similar findings can also be obtained from the finite-system calculations.
This paper is organized as follows. In Sec. II, we describe the theoretical model and methods used in this work. In Sec. III, the results of our calculations are presented and discussed. Finally, we conclude with a brief summary in Sec. IV. The Appendix provides some details on the numerical approach used to simulate the periodic system. Atomic units (a.u.) are used throughout unless stated otherwise.

A. Finite chain model
Let us first revisit the finite chain model and the (TD)DFT approach used in previous works [49,50,53,55]. A theoretical treatment of the corresponding infinitely extended system will be presented in Sec. II B. We consider a linear chain of N equally charged nuclei with a separation a and located at The ionic potential reads where Z is the charge of each ion and is a softening parameter which smoothens the Coulomb singularity in the 1D treatment. In this work we set Z = 4, a = 7 and = 2.25, as in Refs. [49,50,53,55]. For sufficiently large N , this model can qualitatively represent a bulk solid with two fully occupied valence bands. Note that the considered model system is charge and spin neutral. Thus the number of electrons (with spin index ↑, ↓) is N e↑ = N e↓ = 2N for the present case of Z = 4. For brevity, we use the spin-restricted scheme and drop the spin index in the following formulation. According to the Kohn-Sham (KS) theory [63], the field-free electronic state is described by a set of KS orbitals, determined by with the static KS potential With N occ denoting the number of occupied spatial orbitals, the total density is n(x) = 2 and the exchange-correlation potential is treated in a local density approximation (LDA) Following previous works [49,50,[53][54][55], here we use the LDA exchange for the three-dimensional (3D) electron gas as a simple choice to construct our model, which aims to capture the essentials of HHG from a 3D bulk system driven by linearly polarized laser fields rather than to solve an exact 1D system.
The interaction of our model system with laser fields is simulated using TDDFT in the velocity gauge and dipole approximation. For the driving laser pulse linearly polarized along the x axis, we consider the vector potential with ω 0 the angular frequency (photon energy) and n cyc the number of cycles. The laser-driven system is governed by the time-dependent KS equations is determined by the time-dependent density n(x, t) = 2 Nocc j=1 |ϕ j (x, t)| 2 . The numerical calculations are performed on an equidistant grid with spacing ∆x = 0.1, and the number of grid points is set to 28000 for the considered system size N = 200. We propagate the KS orbitals [Eq. (8)] using the Crank-Nicolson method with a predictor-corrector step for updating the KS potential [63,64]. A fixed time step size ∆t ≈ 0.011 (i.e., 1/25000 optical cycle for 2µm wavelength to be considered below) is used in this work. The convergence is further checked against time propagation using the Arnoldi-Lanczos method [65] with adaptive Krylov subspaces. The initial conditions for the TDDFT calculations, i.e., the field-free ground-state KS orbitals are found via imaginary time propagation with orthogonalization in each time step [64,66]. We note that as in previous studies [49,50,53,55], the considered laser interaction does not cause significant changes to the density; therefore we also perform time propagation with the KS potential frozen to its ground-state form, and find that the frozen KS approach captures basically the same HHG features as the dynamic KS approach. For simplicity, all the HHG spectra presented in this paper (including those for the infinite periodic system) refer to the frozen KS simulations, implying an independent-electron picture for which the contributions from different orbitals can be uniquely identified for further analysis (see below).
To calculate the HHG spectra of this model, we compute the total time-dependent current with J j (t) the contribution from the jth spatial orbital: where the factor of 2 accounts for the spin degeneracy. The HHG spectral intensity is then evaluated as the modulus square of the Fourier-transformed current, i.e., S(ω) ∝ dt W (t)J tot (t)e −iωt 2 , where W (t) is a window function introduced to improve the signal-to-noise ratio. A variety of window functions can be chosen for harmonic analysis [67]. For the HHG spectra shown in this paper, we use Blackman windows with the full temporal width set equal to the laser pulse duration.
B. Infinitely extended chain model The above-described finite model, using a selfconsistent many-electron approach, offers a straightforward way to simulate a bulk system by simply increasing the system size. In spite of its straightforwardness, the finite-size model simulations become computationally demanding for large N . If one aims to simulate an ideal crystal, the numerical difficulty of the finite-size model calculations can be avoided by considering an infinitely extended system with periodicity, which is commonly used to model a bulk solid. Moreover, a periodic treatment makes the crystal momentum a good quantum number [68], which could bring convenience for exploring the underlying physics. The infinite chain limit (N → ∞ with the nuclear charge Z and the internuclear separation a fixed) can be achieved by real-space calculations in a single unit cell with the periodic boundary condition imposed.
In the periodic treatment, we introduce the crystal momentum k and the band index m for the KS orbitals. According to Bloch's theorem, we express the field-free orbitals as with u m,k (x) the periodic part u m,k (x + a) = u m,k (x). The specific model considered here has two valence bands fully occupied (N b = 2), therefore the total density in this case is n(x) = (2π/a) −1 π/a −π/a dk 2 is the static periodic effective potential, and ε m,k gives the band structure of the mth band. Note that the ionic potential and the Hartree potential correspond to the infinitely extended limit of Eq. (2) and Eq. (5), which are unbounded with opposite signs as the system size increases, N → ∞. Therefore in practical calculations, we only evaluate their , which has finite values. In our implementation, the ions are simply placed at the center of each unit cell, e.g., a unit cell [−a/2, a/2] contains an ion at x = 0. Due to the periodicity and the charge neutrality, we can directly construct with a convergent infinite sum (for each pair of points x and x ), which can be numerically determined before time propagation of the KS equations.
When including laser interactions in the velocity gauge, the periodicity of the Hamiltonian in space is preserved and the crystal momentum k remains a good quantum number. With the time-dependent KS orbitals expressed as ϕ m,k (x, t) = e ikx u m,k (x, t), the time evolution of the periodic part u m,k (x, t) follows (13) and Eq. (16) can be solved in a single unit cell [−a/2, a/2] with the periodic boundary condition imposed; see Appendix for details. Compared with the finite-system approach described in Sec. II A, the periodic treatment for the infinitely-extended model allows us to restrict the numerical calculations to a single unit cell, which significantly reduces the computational cost. In the periodic case, we calculate the time-dependent current over the first Brillouin zone (BZ) as with J k (t) the crystal-momentum-resolved current again the factor of 2 accounts for the spin degeneracy. For the real-space discretization, we use the same grid spacing as in the finite-size model calculations, i.e., a single unit cell for the considered model (a = 7) is discretized with 70 grid points. Note that the Hamiltonians in Eq. (13) and Eq. (16) explicitly depend on the crystal momentum k. Thus we use the Lanczos time propagator (with adaptive Krylov subspaces) for solving Eq. (13) and Eq. (16), which guarantees that all the KS orbitals are propagated with the same level of high accuracy. This is an advantage over some other commonly-used propagators such as Crank-Nicolson and split-operator schemes, for which the numerical error would be k dependent. We use 400 equally spaced k points to sample the BZ [−π/a, π/a], and check the convergence against calculations using 800 BZ sampling points. All the results discussed below are fully converged.

A. Finite vs Periodic: Band structures and HHG
The finite chain model considered here has been used in previous works for discussing many-electron effects [49], finite-size effects [50] and influence of lattice defects [53,55]. It has been demonstrated to behave like a bulk solid when the system size is sufficiently large. Since it has not been explicitly shown to what extent the finitesize system captures the features of the corresponding infinitely-extended system, here we begin our discussion with a comparison between the finite-system model and its infinite periodic limit.
First of all, we present the band structures of the model in Fig. 1, including two valence bands (VB1 and VB2) and four conduction bands (from CB1 to CB4). Conceptually, the band structures are well defined for a periodic system. For a sufficiently large finite system, the k-space distribution of the KS orbitals can be seen as an approximate representation of the band structures. As shown in Fig. 1, the band structures for the N = 200 system and the infinite periodic limit are basically identical, except that some finite-system features, i.e., the free-space dispersion k 2 /2 and several edge states with "out-of-band" energies, are absent for the periodic system. For typical laser parameters considered so far, it has been demonstrated that the orbitals in the highestoccupied band VB2 dominantly contribute to the total HHG spectra [49]. For the considered finite system with N = 200, the orbitals with indexes from 203 to 400 belong to VB2, and the orbitals with indexes 201 and 202 are actually a pair of edge states with energy slightly below VB2. Similarly, the orbitals with indexes 401 and 402 are also a pair of edge states with energy slightly below CB1. Thus in the periodic treatment, the bandgap between VB2 and CB1 is found to be 0.239, which is slightly larger than the previously-reported value (0.235) [49,53,55] obtained from the energy difference between the lowest-unoccupied and highest-occupied orbitals for the finite system. Nevertheless, as we will see below, the above-mentioned finite-system features play negligible roles in the HHG processes for the considered model.
Having shown that the finite system indeed mimics the infinite periodic limit in terms of the band structures, we now compare the HHG spectra obtained from the finite system and the infinitely-extended periodic system. Band structures for the considered finite-system model with N = 200, obtained from k-space distribution of the KS orbitals. The markers are band structures for the corresponding infinite periodic system in the first Brillouin zone. The plot includes two valence bands (VB1 and VB2, red circles) and four conduction bands (from CB1 to CB4, green squares). While the free-space (FS) dispersion k 2 /2 naturally appears for the finite system, this finite-box feature is absent for the infinite periodic system. Since the response per atom is considered in the periodic treatment, for the sake of comparison, we evaluate the laser-driven current per atom for the finite system, and correspondingly scale the total HHG spectrum of the Natom finite system by N −2 . As one can see from Fig. 2, the HHG spectra obtained from the finite system and the infinite periodic system are in very good agreement, which justifies both approaches. For the considered laser frequency ω 0 = 0.0228 (cor-responding to a wavelength of ∼2µm), a transition from VB2 to CB1 requires at least 11 photons. This regime is similar to the one considered in experiments, where the minimum number of photons required for an excitation from the valence to the conduction band is typically between 10 and 20: A minimum number of 9 photons was required for ZnO [1], 13 for solid Kr [7], 15 for solid Ar [7], and 14 to 19 for bulk GaSe [4,6]. For our model, the harmonics of orders ≥11 are therefore in the abovebandgap regime, where the interband mechanism plays an important role. As shown in Fig. 2, four plateaus can be clearly identified in the spectra (with their cutoffs marked by vertical arrows), corresponding to the recombination transitions from the four conduction bands (CBi, i = 1, 2, 3, 4) to VB2. A detailed analysis of the underlying dynamics will be presented below.
B. Multiple plateaus from a many-electron reciprocal-space perspective The multi-plateau spectral structure, often observed in various numerical simulations, is a signature of the presence of multiple conduction bands in a bulk solid. So far the harmonics beyond the first cutoff have been experimentally measured in solid Ar and Kr [7]. There have been several theoretical explanations for the multiple plateaus and it seems that a consensus has not yet been reached. For example, the experimentally observed multiple plateaus were interpreted from a multilevel perspective by considering only the Γ-point (k = 0) contribution [27] in a model with a direct bandgap at the Γ point. Also restricting to the dynamics of a single Bloch electron at the Γ point, a quasiclassical model was proposed, assuming that the population of higher conduction bands is pumped from the lower band step by step [29]. Accounting for many independent electrons, a k-space semiclassical model was proposed, in which the band-climbing processes are also treated in a step-by-step manner [69]. It should be noted that the Γ-point-only single-electron approximation was demonstrated to disagree with many-electron calculations [49,57]. This point was further emphasized in a recent paper [45], where a preacceleration step was introduced for the electrons with initial crystal momenta away from the Γ point. For the multi-plateau structure observed in Fig. 2, we find that the many-electron k-space semiclassical model introduced in Ref. [69] can indeed provide useful insights. For completeness, we first recall the basic assumptions in this model.
(a) The laser-driven intraband motion of a Bloch electron is described by the semiclassical acceleration theorem, i.e., the instantaneous crystal momentum of an electron is k 0 +A(t), with k 0 the initial crystal momentum and A(t) the vector potential related to the electric field F (t) according to A(t) = − t dt F (t ).
(b) Tunneling excitation to an upper band is considered to occur at a bandgap, e.g., at k = 0 for the VB2-CB1 transition and at k = ±π/a for the CB1-CB2 transition. This is a simplifying approximation relying on the fact that the highest tunneling probability is typically obtained when the energy gap is smallest [70]. In fact transitions in the vicinity of the bandgap are also possible, which we will come back to later.
(c) The interband harmonics might be emitted at any time t as long as there are electrons already excited into conduction bands, and the photon energy is given by the electron-hole energy difference at the instantaneous crystal momentum k(t). This is in contrast to the atomic case where recombination happens when the electron returns to the parent ion in real space. The former assumption about emission times to some extent reflects the spatiallydelocalized nature of electrons in solids. In the semiclassical picture, the electron and hole wavepackets are localized in k space and spread in position over many lattice sites [68], their real-space overlaps will therefore be imperfect more often than not. Note that this description does not completely conflict with the generalized three-step model for HHG in solids [11,21], where electron-hole recollision in real space is usually considered as a condition for the emission of harmonics. In fact the real-space electron-hole recollision assumption will simply select the dominant contributions in some cases; more generally, the electron-hole recombination could happen even in case of an imperfect overlap of the electron and hole wavepackets, which was quantified very recently [46]. Thus the description of harmonic emission in this k-space semiclassical model, albeit simple, can well capture the interband harmonics, as will be shown below.
According to this k-space semiclassical model, electrons can climb up to higher bands in a step-by-step manner, as illustrated in Fig. 3. For example, an electron with k 0 may undergo a VB2-CB1 transition when k 0 + A(t 1 ) = 0, and at a later time it continues to climb onto CB2 when k 0 + A(t 2 ) = π/a. With additional optical cycles, climbing onto higher bands is possible, as long as k 0 + A(t) = 0 or ±π/a can be fulfilled. Note that the vector-potential amplitude A 0 considered in this work fulfills 0.5π/a<A 0 <π/a, for which a part of electrons will reach higher conduction bands via the band-climbing processes and the HHG spectra will then show clear multiplateau structures. By scanning the tunneling time (denoted by t s ) for the VB2-CB1 transition and accounting for all the possible paths in k space, we obtain a full set of trajectories. With the assumption of bandgap transitions, a specific electron with initial k 0 = −A(t s ) is automatically selected when scanning the VB2-CB1 tunneling time t s , which implies a many-independent-electron picture. We note in passing that in this semiclassical model, the laser-induced energy shifts of bands are neglected for simplicity, i.e., the considered laser-interaction regime leads to only a small portion of electrons excited into the conduction bands, causing no significant changes to the band structures during the laser pulse.
To examine the correctness of this k-space trajectory analysis, we compare the temporally resolved harmonic The band-climbing processes illustrated by considering the central optical cycle (approximately treated as a sinusoidal form for simplicity). In the presence of the laser field, the instantaneous crystal momentum of an electron is k0+A(t), with k0 the initial crystal momentum and A(t) the vector potential (see text). The considered vector-potential amplitude A0 fulfills 0.5π/a<A0<π/a, which applies to all the results presented in this paper. Tunneling excitation to an upper band is considered to occur at a bandgap, e.g., at k(t) = 0 (±π/a) for the VB2-CB1 (CB1-CB2) transition. Due to the form of k(t), there are typically two moments for a bandgap transition within one optical cycle. Similar processes can happen among higher conduction bands with additional optical cycles, e.g., after being excited into CB2, the electron may undergo a CB2-CB3 transition when k(t) reaches 0 again. Inset: the vector potential in one optical cycle, with the circles (squares) indicating the moments of the VB2-CB1 (CB1-CB2) transitions for the depicted paths. emission of the semiclassical trajectories with the timefrequency profile of the HHG spectrum. The latter, for the finite system and the infinite periodic system, is extracted by performing a Gabor transform of the timedependent current Eq. (10) and Eq. (17), respectively: where the width of the time window τ is chosen to be 4π (a.u.). Figure 4 shows the Gabor-transformed spectrum |G cell (ω, t)| 2 for the periodic-system current. The timefrequency profile of the finite-system HHG is very similar (not shown). In the semiclassical trajectory analysis, the band structures of VB2 and four conduction bands (CBi, i = 1, 2, 3, 4) are taken into account, which are sufficient to describe the four plateaus observed in the spectra. The step-by-step band-climbing processes are indeed reflected in the time-frequency profile of the HHG, i.e., the higher plateaus emerge later than the lower one.
The temporally resolved harmonic emission (i.e., the harmonic energy as a function of time) described by each k-space trajectory follows the instantaneous energy difference between the conduction and valence band structures, E cv [k(t)]. Note that in the semiclassical trajectory analysis with multiple conduction bands taken into account, the conduction band involved in E cv [k(t)] is the one on which the electron moves at time t. As shown in Fig. 4, the k-space trajectories, without any constraint on "real-space recollision" for the electron-hole pairs, are found to agree well with the time-frequency profile. The electron-hole "real-space recollision" time (denoted by t r ) is typically determined by [11,21,37] Although this equation is sometimes interpreted as "realspace recollision" of the electron-hole pairs, it actually stems from the saddle-point method on the basis of spatially-delocalized Bloch waves [11,21,37]. It therefore seems appropriate to see this as a constraint for selecting some significant contributions, rather than a conventional real-space picture as in the atomic and molecular cases. Applying the additional constraint of Eq. (20) results in a subset of trajectories covering much fewer harmonic emission events, as also shown in Fig. 4. Therefore Eq. (20) seems to be a too restrictive condition, as also pointed out in Ref. [71]. In the semiclassical description, as noted above, the spatially-delocalized nature of the electron and hole wavepackets will inevitably cause electron-hole recombination due to imperfect wavepacket overlaps, which are not captured by Eq. (20). Since the k-space trajectory analysis (without the additional "realspace recollision" constraint) works well for illustrating the multi-plateau HHG, here we do not go into details of the trajectories constrained by Eq. (20). A detailed analysis of the "real-space recollision" trajectories was presented in a recent paper [62], with three bands (VB2, CB1 and CB2) considered. Nevertheless, it should be emphasized that the so-called "real-space recollision" model is simply a subset of the k-space trajectories. Note that a statement in Ref. [62] that the k-space trajectories disagree with the time-frequency profile turns out to be invalid: the k-space trajectories shown in Fig. 3(c) of Ref. [62] did not even include the "real-space recollision" trajectories in Fig. 3(b) of Ref. [62] as they should, implying that the k-space trajectory analysis presented in Ref. [62] seems to have missed some essential ingredients.
Using the same TDDFT model, we show in Fig. 4 that a properly-implemented k-space trajectory analysis captures the time-frequency profile of the multiple plateaus.
We also see from Fig. 4 that the k-space trajectories do not cover some relatively weaker harmonic emission before the center of the pulse, e.g., do not cover the harmonics beyond the first cutoff emitted during the 7th optical cycle. This disagreement is due to the aforementioned assumption of the bandgap transitions: according to this simple approximation, the initial crystal momen- tum of an electron tunneling into CB1 during the early part of the pulse lies in a limited range around k 0 = 0, then the corresponding instantaneous crystal momentum k(t) cannot reach the BZ edge (±π/a) to tunnel into CB2 for the considered A(t). Accounting for tunneling transitions in the vicinity of the bandgap in the semiclassical model could therefore lead to a better description. Also, we can infer that the electrons with different crystal momenta behave differently in the band-climbing processes and may contribute to different parts in the HHG spectra. This will be further discussed in the following subsection.

C. Contribution of different orbitals to HHG
With the insights gained from the k-space semiclassical model, we study how the electrons with different initial crystal momenta contribute to the HHG spectra. Recently crystal-momentum-resolved contributions were investigated for the first-plateau regime of a Kronig-Penney model in Ref. [72]. In the present work, we perform this type of analysis for the multi-plateau spectra of our model. Correspondingly the laser vector-potential amplitude considered here fulfills 0.5π/a<A 0 <π/a, for which it is possible to find a range of crystal momenta responsible for the high-order harmonics beyond the first cutoff.
The considered 15-cycle pulse is relatively long such that the carrier-envelope-phase effects play a very minor role. The k-space semiclassical model implies that the electrons tunneling into CB1 should have initial crystal momenta k 0 in an interval [−A 0 , A 0 ]. In pursuit of a better description, here we consider an extended region [−A 0 −0.01π, A 0 +0.01π], with ±0.01π a rough estimate of the small k range for tunneling transitions near the bandgap. Similarly, for tunneling into higher conduction bands, we have an additional condition |k 0 |≥(π/a−A 0 −0.01π). Having the step-by-step band-climbing picture in mind, we expect the region R II = {k 0 | (π/a−A 0 −0.01π)≤|k 0 |≤(A 0 +0.01π)} to capture the second and higher plateaus, and the region R I = {k 0 | |k 0 |≤(π/a−A 0 −0.01π)} to dominantly contribute to the first plateau. As shown in Fig. 5, such an analysis is supported by the crystal-momentum-resolved profile of the periodic-system HHG spectrum, which is obtained from the Fourier transform of Eq. (18). Note that the notation k used in Sec. II B stands for the initial crystal momentum (k 0 ) and Eq. (18) naturally defines the contribution from a specific initial crystal momentum. The total spectrum is a coherent sum of all the contributions, however, the spectral intensity of individual contributions [ Fig. 5(a)] clearly shows that the second and higher plateaus mainly originate from the region R II . This is further confirmed in Figs. 5(b) and (c) by comparing the spectra obtained from the current calculated within the above-defined k 0 regions [i.e., replacing the integral in Eq. (17) π/a −π/a dk → R I(II) dk]. Indeed, R II captures all the harmonics beyond the first cutoff and R I roughly captures the first plateau. We note that while the region R II nicely reproduces the spectrum beyond the first cutoff, it also contributes to a small part of the first plateau (see, e.g., the "tail" near the first cutoff at order ∼25). This is because in the step-by-step band-climbing picture, before an electron reaches higher conduction bands, it must undergo some motions within the first conduction band, partially contributing to the first plateau by CB1-VB2 recombination.
As another demonstration of the k 0 -resolved contributions, we present in Fig. 6 the crystal-momentumresolved profile for a different set of laser parameters A 0 = 0.3 and ω 0 = 0.01824, with the peak intensity unchanged and the wavelength increased to ∼2.5µm. The definition of k 0 regions R I and R II also applies when varying the laser parameters. A 0 turns out to be an important parameter which governs the many-electron dynamics: increasing A 0 in the range 0.5π/a<A 0 <π/a makes the region R I (R II ) smaller (larger). It could be interesting to note that A 0 is also directly linked with the Keldysh parameter [70], which typically classifies the laser-matter interaction regimes. When increasing the wavelength at a fixed field strength, the Keldysh parameter decreases and the interband excitation becomes more deeply into the tunneling regime. We mention in passing that increasing the wavelength while keeping A 0 = 0.24 fixed leads to very similar observation as Fig. 5 except for lower HHG signals due to the weaker laser field (not shown).
In addition, we mention that as an analogue to the crystal-momentum-resolved contributions for the infinite periodic model, the so-called orbital-index profile of the finite-system HHG introduced in Ref. [55] provides similar information when using the finite-system model. Recall that the band structures can be approximately represented by the k-space distribution of the KS orbitals for a sufficiently large finite system (Fig. 1), therefore the above-defined k 0 regions can be mapped, according to the valence band structure, to the orbital-index ranges for the finite-system model. The orbital-index profile of the finite-system HHG spectrum, obtained from the Fourier transform of Eq. (11), is presented in Fig. 7. Here we only focus on the orbital contributions in VB2 (i.e., orbital indexes from 203 to 400), since the HHG contributions from other orbitals are almost negligible. Similar to the k 0 regions defined for the infinite periodic system, the orbital-index ranges R I and R II for the finite system are found to be responsible for harmonics in the first plateau and the higher plateaus, respectively.
Identifying relevant orbital contributions for the HHG spectral structures is helpful for understanding the underlying processes from a many-electron perspective. It could also lead to a reduction of computational cost, since one may skip the time propagation of the less important orbitals in some cases, which would be rather advanta- geous for 2D or 3D models. Regarding the latter point, however, we note that while a set of significant orbitals can be defined for reproducing the interband harmonics as shown in Figs. 5-7, a correct description of the intraband harmonics in the below-bandgap regime still requires the contributions from all the occupied orbitals to be coherently taken into account, as also found in Ref. [43].

IV. CONCLUSION
Using a self-consistent TDDFT approach, we numerically studied the HHG from a model solid with a bandgap irradiated by linearly polarized laser pulses, in particular the multi-plateau spectral structures. We supplemented the previous model used in Refs. [49,50,53,55] with an infinite periodic treatment, and demonstrated that a sufficiently large finite model indeed mimics the infinite periodic limit in terms of the band structures and the HHG spectra. On the basis of a k-space semiclassical trajectory analysis, the multi-plateau spectral feature was well explained from an independent-electron perspective, where electrons with different crystal momenta contribute differently to the total HHG spectra. It may be worthwhile to point out that the concept of semiclassical electron-hole recollisions for HHG in solids should only be considered with the spatially-delocalized nature of the Bloch waves kept in mind: an exact coincidence of the electron and hole wavepacket centers in real space could be a too restrictive condition, which cannot capture the harmonic emission in the scenario of imperfect electron-  Fig. 5 to the VB2 orbital energy ranges R I and R II , according to the band structure. (b) Orbital-index profile, corresponding to the finite-system HHG spectrum shown in Fig. 2 (without the scaling factor N −2 ). The dashed and solid lines indicate the VB2 orbital energies corresponding to the k0 values described in Fig. 5. For this finite system with N = 200, the regions R I and R II are defined by indexes from 323 to 400 and indexes from 281 to 322, respectively. (c) and (d) HHG spectra obtained from the current calculated within the orbital-index ranges R I and R II , respectively. The total spectrum is also shown for comparison, with the vertical arrows indicating the cutoffs of the four plateaus as in Fig. 2. hole wavepacket overlaps.
With the infinite periodic extension of the TDDFT model, we investigated the crystal-momentum-resolved contributions in the HHG processes, and identified the relevant initial-crystal-momentum (k 0 ) regions for different parts of the HHG spectra. Specifically, for the considered laser-interaction regime, the multiple plateaus beyond the first cutoff mainly originate from electrons with k 0 away from the Γ point, and the first-plateau harmonics are predominantly contributed by electrons with k 0 around the Γ point. Such particular k 0 ranges can be controlled by the vector potential. Similar findings were also obtained from the finite model calculations, keeping in mind that the band structures can be ap-proximately represented by the k-space distribution of orbitals for a sufficiently large finite system. This work highlights the importance of many-electron contributions in the solid-state HHG processes. Moreover, it shows the possibility of isolating the significant contributions to some particular spectral features based on a semiclassical independent-electron picture, which lays the foundation for possible control of the HHG spectra by laser parameters to be explored in future studies. In practical calculations, the wavefunctions are discretized as finite-size vectors and all the operators are expressed in terms of matrices. In general, a detailed implementation of boundary conditions depends on the specific discretization scheme used for numerical calculations. For completeness, here we briefly present how the 1D periodic boundary condition is implemented in this work.
The periodic treatment means that we only deal with periodic functions in practice, e.g., the periodic part of a Bloch wave which fulfills u(x) = u(x + a). To account for such a periodicity, we use a uniform n-point grid of size a with spacing ∆x = a/n. Due to the periodicity, the grid-point positions can be chosen with an arbitrary shift; here we simply put them symmetrically in a unit cell centered at the origin [−a/2, a/2]: x i = [i − (n + 1)/2]∆x, (i = 1, · · · , n). Imagine that we apply this procedure to all the unit cells, the entire x axis would then be discretized by an infinitely-extended uniform grid. As a straightforward approach, the wavefunctions and potentials can be represented by their values on grid points; the remaining task is therefore to approximate the first-and second-order derivatives in discretized form.
Our implementation uses the matrix Numerov method [73], which is a finite-difference scheme suitable for Schrödinger equations. To introduce the details, let us begin with the simplest central-finite-difference method, in which the first-and second-order derivatives at the ith point u i and u i are expressed as When restricting to the grid in a single unit cell, one should pay attention to the boundary grid points x 1 and x n . We denote the grid point to the left of x 1 by x 0 and the one to the right of x n by x n+1 , respectively. Clearly, the periodic ansatz implies u 0 = u n and u n+1 = u 1 . Thus with the periodic boundary condition, one can express the first-and second-order differential operators as n×n matrices: The finite-difference method described above is for illustrative purpose only, which shows how the periodic boundary condition is imposed by setting the corner elements of the matrices. In our practical calculations, we achieve higher accuracy [O(∆x 4 )] by using an implicit three-point stencil in the spirit of the Numerov method [64]. This elegant idea has been underlying the basics of some TDSE solvers [74][75][76]. Here we apply it to the periodic case; correspondingly, the matrix-form operators are modified as