Extended high-harmonic spectra through cascade resonance in confined quantum systems

The study of high-harmonic generation in confined quantum systems is vital to establishing a complete physical picture of harmonic generation from atoms and molecules to bulk solids. Based on a multilevel approach, we demonstrate how intraband resonances significantly influence the harmonic spectra via charge pumping to the higher subbands and, thus, redefine the cutoff laws. As a proof of principle, we consider the interaction of graphene nanoribbons, with zigzag as well as armchair terminations, and resonant fields polarized along the cross-ribbon direction. Here, this effect is particularly prominent due to many nearly equi-separated energy levels. In such a scenario, a cascade resonance effect can take place in high-harmonic generation when the field strength is above a critical threshold, which is completely different from the harmonic generation mechanism of atoms, molecules and bulk solids. We further discuss the implications not only for other systems in a nanoribbon geometry, but also systems where only a few subbands (energy levels) meet this frequency-matching condition by considering a generalized multilevel Hamiltonian. Our study highlights that cascade resonance bears fundamentally distinct influence on the laws of harmonic generation, specifically the cutoff laws based on laser duration, field strength, and wavelength, thus unraveling new insights in solid-state high-harmonic generation.

The study of high-harmonic generation in confined quantum systems is vital to establishing a complete physical picture of harmonic generation from atoms and molecules to bulk solids. Based on a multilevel approach, we demonstrate how intraband resonances significantly influence the harmonic spectra via charge pumping to the higher subbands and, thus, redefine the cutoff laws. As a proof of principle, we consider the interaction of graphene nanoribbons, with zigzag as well as armchair terminations, and resonant fields polarized along the cross-ribbon direction. Here, this effect is particularly prominent due to many nearly equi-separated energy levels. In such a scenario, a cascade resonance effect can take place in high-harmonic generation when the field strength is above a critical threshold, which is completely different from the harmonic generation mechanism of atoms, molecules and bulk solids. We further discuss the implications not only for other systems in a nanoribbon geometry, but also systems where only a few subbands (energy levels) meet this frequency-matching condition by considering a generalized multilevel Hamiltonian. Our study highlights that cascade resonance bears fundamentally distinct influence on the laws of harmonic generation, specifically the cutoff laws based on laser duration, field strength, and wavelength, thus unraveling new insights in solid-state high-harmonic generation.
Depending on the intensity and frequency of the applied electric field, together with details of the band structure, even non-perturbative mechanisms can be responsible for the solid-state HHG [52]. In this case, solid-HHG is believed to have contributions from dynamical intraband and interband processes involving a k-space motion of Bloch electrons, typically described by a threestep model [12,13]. Nevertheless, a clear understanding For a confined quantum system, due to lack of translational symmetry, the spectrum consists of only a discrete set of energy levels instead of bands. It is known that the HHG spectra in quantum dots can be influenced by confinement conditions such as size and/or coupling parameters in coupled quantum dots [54,55]. One of the distinct advantages of confined quantum systems is the tunability of its energy spectra and wavefunctions by tailoring its size [56,57], or by external parameters, such as gate voltage or magnetic field [58][59][60].
A rather interesting situation is realized in partially confined systems, such as quasi-one-dimensional (quasi-1D) systems which are confined in one dimension but periodic in other resulting in a series of subbands in the band structure [61]. Therefore, such systems possess properties of both bulk and finite systems. The most notable example is perhaps the graphene nanoribbon (GNR) [62][63][64][65][66][67]. Depending on the polarization of the electric field, different aspects of HHG, corresponding to the bulk and/or the confined quantum systems can be explored. Arguably, the phenomenology of HHG in such a confined quantum system, therefore, forms a bridge between atoms, molecules and bulk solids. A detailed understanding of the underlying mechanism in such quantum systems may not only shed light on the HHG phenomenology in bulk solids but also unravel new arXiv:2112.08790v1 [cond-mat.mes-hall] 16 Dec 2021 phenomenology and HHG response in confined systems.
Understanding the HHG phenomenology in such systems necessitates going beyond the two-level (equivalent to the commonly employed two-band studies of the bulk solid-state HHG), thereby raising a number of conceptual questions: For a laser field polarized along the confining direction, what would happen if the laser frequency ω 0 approximately matches the subband gap, i.e., ω 0 ≈ ∆ sg ? Does it cause a resonance resembling Rabi flopping in a two-level system [68,69]? If many such (approximately) equi-separated subbands are present, does the resonance involves multiple subbands? If so, could this resonant excitation affect the harmonic process remarkably? In addition, Hansen et al. indicate that the observed cutoff law on HHG transits from atomic to solid-state type with increasing system size in a model 1D chain [70]. Does this phenomenon hold even in the presence of resonance?
Here, by considering a multilevel model and solving the time-dependent Schrödinger equations, we demonstrate that the subband resonance can lead to remarkable effects in the HHG spectra. We find that this effect is especially prominent in GNRs due to the presence of nearly equiseparated bands, To this end, as a proof of principle, we consider GNRs interacting with a resonant field polarized along the cross-ribbon (confined) direction where it leads to cascade resonance, a multilevel phenomena involving almost all the valence and conduction subbands, and an enhanced HHG spectra well beyond the current cutoff laws.
Specifically, the plateaus of the harmonic spectra are broadened significantly when the laser frequency matches the subband gap. The cascade resonance effect causes the excited electrons to gradually accumulate near the Dirac points also in higher-energy subbands (charge pumping), eventually leading to the extended HHG spectra. Subsequently, we establish the conditions for occurrence of cascade resonance in GNRs. Our analysis indicates that the occurrence of cascade resonance requires a threshold field strength in addition to satisfying the frequencymatching condition. In addition, we formulate the dependence of harmonic cutoff in GNRs on laser duration, field strength and wavelength, when the cascade resonance occurs. These dependencies are significantly different from those of HHG in atoms, molecules, and bulk solids currently studied. Finally, we discuss the possibilities for the cascade resonance in other materials by considering a general multilevel model and implications for HHG spectra.

II. HHG IN GRAPHENE NANORIBBONS
GNR is a quasi-1D material extending in two directions-ribbon (x) and cross-ribbon (y) directions [62][63][64][65][66][67]. GNR with armchair edges (AGNR) on both sides is classified by the number of dimer lines (N a ) across the ribbon width as shown in Fig. 1(a). The unit cell of AGNR consist of two chains P and Q. Likewise, GNR Band structure of (c) 10-AGNR and (d) 10-ZGNR. Red (blue) curves stand for the subbands belonging to the conduction (valence) band. ∆sg denotes the subband gap. Ka = 0 and Kz = 2π/3dz, 4π/3dz are Dirac points, respectively, for the armchair and zigzag GNRs when the boundary conditions of the transverse (y) direction are periodic. wa (wz) is the width of armchair (zigzag) nanoribbon. da (dz) is the distance of the armchair (zigzag) unit cell.
with zigzag shaped edges (ZGNR) on both sides is classified by the number of the zigzag chains (N z ) across the ribbon width [ Fig. 1(b)]. The chain in unit cell of ZGNR is labelled by C. We refer to GNR with N a armchair dimer lines as N a -AGNR and GNR with N z /2 zigzag chains as N z -ZGNR. Tight-binding model of GNR is presented in Appendix A.
Since the ribbon is macroscopically large along the x (longitudinal) direction, continuous band structure can be obtained across the Brillouin zone (BZ) shown in Figs. 1(c) and 1(d). By contrast, in the y (transverse) direction, transverse confinement gives rise to a discrete set of subbands which are one of the typical features for electronic structure of nanoribbons. Following previous convention, we adopt J(s) as index notation for the ribbon subbands, see Fig 1(c) and 1(d), where J = 1, 2, . . . , N is the band number and s is the band type with "c" and "v" representing conduction and valence bands, respectively. Different from subbands of ZGNR, the subband indices of the AGNR in Fig. 1(c) are classified into two groups labelled by {J} and {J } respectively, where subbands J(s) and J (s) merge at the boundaries of BZ. Atomic units are used throughout the paper unless otherwise indicated.
Recently, there has been a growing interest in study-ing HHG from GNRs or similar nanoribbons [71][72][73][74][75], owing to diverse electronic properties of GNRs, which arise particularly from different edge geometries, viz. zigzag and armchair edges. When the laser field is polarized along the ribbon direction, the bulk aspects of GNRs are reflected in the harmonic spectra. For example, the edge states of ZGNR enhance the emitting efficiency of low-order harmonics [72]; the onsite potential, breaking the mirror symmetry, causes the perpendicular harmonic emission [74,75]. It is easy to know that the cutoffs of longitudinal harmonic spectra scale linearly to the field strength and wavelength as shown in Appendix D, which is similar to the HHG in bulk graphene. However, the finite-size effects [70,76] in GNR along the cross-ribbon direction has not yet been fully examined. In GNRs, the HHG generated by transverse field can arguably be more interesting than generated by longitudinal field. Quantum confinement effects reflect on HHG when the applied laser field is polarized along crossribbon direction. Specifically, the nearly equal-energy spacing subbands play a role and induce a resonance excitation over the subbands. In general, however, optical transitions between two subbands are not always allowed. The optical selection rules for GNRs is a result of the wave function parity factor (−1) J , where J is the subband index which has been mentioned before. A detailed derivation and discussion of the optical selection rules is outlined in Appendix B.
A. Size dependence of the transverse HHG spectra The subbands near the Dirac points are approximately equidistant. We denote the gap between nearestneighbor subbands as ∆ sg , shown in Figs. 1(c) and 1(d). It is characterized by ∆ sg ∼ N −1 . We can see in Figs. 2(a) and 2(b) that both (∆ max sg ) −1 and (∆ ave sg ) −1 for AGNR and ZGNR are approximately linear with number of sites N , where ∆ max sg is the maximum subband gap, and ∆ ave sg is the average value of all subband gaps. Both ∆ max sg and ∆ ave sg are presented since the subband gaps at the Dirac points are not precisely equal.
Figures 2(c) and 2(d) display the size dependence of harmonic spectra for AGNR and ZGNR, respectively, while the driver field is polarized along the transverse direction. The laser frequency and intensity are fixed at ω 0 = 0.41 eV (0.0152 a.u.) and I 0 = 5.04 × 10 10 W/cm 2 , respectively. In the spectrum of AGNR-HHG as shown in Fig. 2(c), the harmonic cutoff is drastically extended for widths N a = 40, 128 and 220 (ribbon width w a ≈ 4.8 nm, 15.6 nm and 26.9 nm), forming three peaks p 1 , p 2 and p 3 respectively. These three peaks are depicted in Fig. 2(a) as red stars which abscissas correspond to ∆ sg = ω 0 , ω 0 /3 and ω 0 /5. For the first peak, at N a = 40, it is found in Fig. 2(a) that ∆ ave sg (40) ω 0 ∆ max sg (40), i.e., the laser frequency matches subbband gap [ω 0 ≈ ∆ sg (40)]. Therefore, the harmonic spectrum is affected significantly when the resonance condition (ω 0 = ∆ sg ) is met. Two addi- tional peaks at N a = 128 and 220 reveal that the resonant excitation occurs not only when ∆ sg (40) ≈ ω 0 , but also when ∆ sg (N a ) ≈ ω 0 /3 and ω 0 /5. This can be understood by the optical selection rule of AGNR in Appendix B when subbands are from the same group. The intraband (interband) transitions are permitted for such two subbands which band index difference ∆J = odd (∆J = even). Their energy gaps have such a relation: E Js − E J s ≈ (2j − 1)∆ sg , where j is an integer. Thus the resonance condition for AGNR can be written as ∆ sg (N a ) ≈ ω 0 /(2j − 1). Since intraband (interband) transitions are not permitted for ∆J = even (∆J = odd), we are finally not able to observe such peaks at N a = 84 and 172 satisfying ∆ sg (N a ) ≈ ω 0 /2j. Moreover, from p 1 to p 3 , the peaks gradually widen and their cutoff progressively shrinks, as shown in Fig. 2(c). The former can be attributed to the fact that the subband gap is not very sensitive to the change in system size when the size is large enough. The latter stems from the increase in detuning between the driving frequency ω 0 and the excitation gap (2j − 1)∆ sg , deviating from the resonance condition. Analogous to AGNR-HHG, the ZGNR-HHG spectrum shows similar features in Fig. 2(d) and can be explained in a same way.
To gain deeper insights into the mechanism of cutoff extension on the harmonic output, we will study the electron dynamics in the presence of the nearly resonant field. We will show how the electrons are excited to the highest conduction subbands emitting high-energy photons.

A. Electron dynamics in GNRs under cascade resonance
Taking 36-ZGNR as an example, we study the HHG process combined with its electron dynamics under the resonance. To rule out the effect of pulse envelope on quantum paths, a 10-cycle trapezoid envelope with an one-cycle linear ramp is applied. Figure 3(a) displays the harmonic spectrum when laser frequency matches the subband gap. It can be seen that the cutoff can reach at the 27th order, whereas such a high-order harmonic cannot be observed if the laser is polarized along the ribbon direction. The time-frequency distribution of the corresponding time-dependent current is presented in Fig. 3(b), demonstrating the harmonic emission with subcycle temporal resolution. The maximum order of the emitted harmonics increases with the time evolution at the beginning 3.5 cycles (shadow area) and then keeps at 27th order for the following 6.5 cycles. We thus refer to the first 3.5 cycles as the rise stage and the next 6.5 cycles as the oscillation stage.
As shown in Figs. 3(c1)-(c2), after being excited to the conduction band, the electrons move further up to the higher subbands with increasing vector potential |A(t)|. When |A(t)| starts to decrease, however, the electrons near the Dirac points seem to do not deexcite to the lower subbands as those electrons governed by Bloch acceleration theorem do. Most of them tend to stay in place and wait for the next half cycle of the laser pulse. Once |A(t)| increases again, these (deposited) electrons will be excited to the higher subbands. As the electrons are excited and then deposited, accumulation zones are formed gradually near the K z and K z . The rise region of the harmonic emission in Fig. 3(b) is exactly ascribed to the cascade excitation process near the Dirac points. When all the subbands are involved (at about t = 3.5 T ), see Fig. 3(c3), the electrons are no longer excited to the higher subbands, but jump up and down in the accumulation zones, thus leading to the oscillation region in Fig. 3(b). We thus refer to it as cascade resonance. In AGNR, the resonant dynamics is similar. We provide the dynamic images of the conduction band population of 40-AGNR and 36-ZGNR in the supplementary material [77], so as to gain more insight into the electron resonant dynamics. Note that the electron dynamics in valence band is not shown since the distribution of the holes on the valence band is symmetric to the electron distribution on the conduction band.
With a long enough laser pulse, it is not difficult to foresee that the electrons can be driven to the highest subbands, generating high-order harmonics which energy equals the maximum band gap. This explains the significant cutoff extension on the size-dependent harmonic

B. Analysis of cascade resonance in simplified model
In order to analyze the cascade resonance near the Dirac points of GNRs and generalize the conclusions to other confined systems, we introduce a N -level model here which is widely used to simulate electron dynamics and HHG in quantum dots [54,78,79]. The Hamiltonian In the dipole approximation, the time-dependent Hamiltonian of the laser field with the N -level quantum dot is written as where electric field E(t) = E 0 cos(ωt)ŷ, and d ij is the dipole matrix element between states | i and | j . Rabi frequency is defined as Ω ij = E 0 · d y ij . In principle, the Hamiltonian Eq. (1) can be used to simulate electron dynamics and HHG in any confined system including GNRs along transverse direction, as long as the energy levels and dipole matrix elements are known.
In our simulation, energy levels E 1 , E 2 , . . . , E N/2 are involved as the initial state. Solving the time-dependent Schrödinger equation, the dipole moment can thus be

Equi-spaced N-level model
First of all, we consider an ideal situation: the N -level system with equal separations where i = 1, 2, . . . , N − 1, and ∆ sg is the gap between energy levels. The dipole matrix elements (in atomic unit) take the form of It is constructed based on the features of dipole matrix elements in GNRs, for example, Eq. (B2). For N = 2, the model is reduced to the classical Rabi model which has been investigated extensively before. In the resonant case, ω = ∆ sg , the parameter is commonly adopted to identify different interaction regimes, where Ω R = Ω 12 = E 0 · | d 12 |. In the weak coupling regime γ R 1, the Rabi flopping can be observed that the population inversion oscillates between −1 and 1 at frequency Ω R periodically [68], and thus the area theorem is valid [80]. While transforming into the strong coupling regime γ R 1, the area theorem breaks down [81][82][83]. A chaotic oscillation mode displaces the periodic mode since the contribution of counterrotating terms becomes prominent, which is known as the carrierwave Rabi flopping [84][85][86].
In a system with N > 2, we need to reanalyze its resonant dynamics, because the electron dynamics of the two-level system no longer applies. For comparison, the concepts of weak coupling, strong coupling, and coupling parameter are borrowed. In multilevel systems, we define the coupling parameter as where . . . denotes the average value, d ij is the dipole matrix element between the two energy levels which meet the resonance condition in frequency, and Rabi frequency Ω R = | d ij | · E 0 . We focus on the case ω = ∆ sg in the following because other resonance cases [ω = (2j−1)·∆ sg , j > 1] are similar. Taking N = 20 as an example, the resonant dynamics on conduction levels is shown in Figs. 4(a1)-(a2). We set the energy separation ∆ sg = 1 a.u., and the driver frequency ω = 1 a.u.. The valence population is not shown because it is symmetric to the conduction population and satisfies Fig. 4(a1), the excitation process can be divided into two stages: rise and oscillation, similar to what has occurred in the laser-driving GNR. In the rise region (shadow area), we can find that the electrons are excited to the conduction subbands in a cascaded way. Correspondingly, the total conduction population (P c = N i=N/2+1 | C i | 2 ) exhibits a monotonous increase shown by the black dashed line in Fig. 4(a1). Thus the rise time (T r ) is defined as the period from the beginning of the evolution to the first time that P c reaches the maximum. After entering in the oscillation region, the population of conduction levels starts to collectively oscillate up and down with a roughly uniform period T osc ≈ 35 T , where T = 2π/ω. Nevertheless, this periodicity is not strict and the population does not return to the initial value even for long enough pulses. Hence the area theorem fails for cascade resonance even in the weak coupling regime. For strong coupling γ R = 2 (E 0 = 2 a.u.), the population dynamics is unaltered compared to the case of γ R = 0.2 a.u., as shown in Figs. 4(a2). The only significant change is that the oscillation period T osc becomes 3.5 T . Therefore, the regime transformation in the equispaced system does not fundamentally shift the behavior of resonant excitation. This is completely different from the two-level resonance.
Based on the above analysis, it is found that the cascade resonance can take place in equi-spaced N -level systems as long as the frequency conditions are met, no matter how weak the electric field is. Due to all the energy levels are involved in the cascade resonance, the cutoffs of harmonic spectra can thus extend to maximum energy gap (E 20 − E 1 ). Therefore, we can see from where δ is the scale factor. This equation indicates a fact that the weaker the field strength and the more the energy levels, the longer the rise time. Reflected in the HHG process of the GNR, it implies that the larger the size and the weaker the field strength, the longer the rise region in harmonic emission.

ZGNR-like N-level model
Now we turn to a more realistic situation taking into account the detuning between the driver frequency and the energy separation. Simulating the subband structure of ZGNR at the Dirac points, energy levels have the form: where f and g are the parameters which determine the nearest-neighbor energy difference. In the following simulation, we set N = 20, ω = ∆ sg = 1 a.u., f = 1.15, and g = 1.1. The energy difference near the Fermi level is E 11 − E 10 = 1.15 ω, slightly deviating from the driver frequency. The dipole matrix elements we used the same as in Eq. (3).
In this ZGNR-like system, we focus on the population dynamics in different coupling regimes as before. In the weak coupling regime (γ R = 0.2), see Figs. 4(b1), the cascade resonance disappears but is displaced by a behavior close to the classical Rabi flopping. This is called near resonance. Only four of conduction levels near the Fermi level are involved, which implies that only a small proportion of the electrons can be driven to the highest energy level. Therefore, the corresponding HHG spectrum (red dashed curve) in Fig. 4(b3) terminates at 7th order which energy is much lower than the maximum gap [(E 20 − E 1 )/ω = 15.37]. However, in the strong coupling regime (γ R = 2), the collective resonance reappears in the ZGNR-like system, as shown in Figs. 4(b2). The rise and oscillation stages constitute the entire excitation process again. Correspondingly, the cutoff of harmonic spectrum can thus extend to 15th order shown as the blue curve in Fig. 4(b3). Therefore, the cascade resonance takes place merely for γ R 1 when the field detuning exists. This is different from the population dynamics in the equispaced model.

C. Resonance conditions for laser-driving GNRs
Combining the analyses of resonant dynamics in the 36-ZGNR (Sec. III A) and in the simplified models (Sec. III B) together, it is natural to infer that one of the condition for the occurrence of cascade resonance in the laser-driving GNRs is γ R 1. We verify the coupling parameter at p 1 in the ZGNR-HHG spectrum [ Fig. 2(d)], and find that γ R (N z = 36) ≈ 1.3 > 1 at the Dirac points.
Thereby, the resonance conditions for a laser-driving GNR can be summarized as where j is an integer, and ∆J is the index difference of resonant subbands. Equation (8a) is the frequency condition we have obtained in Sec. II A. It determines which driver frequency can lead to cascade resonance in a GNR with size of N . Then equation (8b) states that the occurrence of cascade resonance requires the laserdriving GNR to be in the strong coupling regime, and γ R = 1 is the critical point. Further, we define a concept -critical field strength: which is the minimum field strength resulting in the cascade resonance in a N -GNR. Substituting Eq. (8a) into Eq. (9), the critical field strength can be evaluated by where (α a , β a , γ a ) ≈ (0.025, 0.039, 0.014), and (α z , β z , γ z ) ≈ (0.021, 0.04, 0.014). It exhibits a dependence of E cri (N ) ∼ N −2 in the large N limit. It clearly indicates that the cascade resonance induced harmonic cutoff extension is more easily observed in large size GNRs.

D. Electron dynamics and HHG away from resonance condition in GNRs
A natural question at this stage is how sensitive is the cascade resonance to the frequency-matching and fieldstrength conditions.
In order to examine the frequency-matching condition, we choose the 80-ZGNR (∆ sg ≈ 0.2 eV) but use the laser field with wavelength of 3 µm (ω ≈ 0.41 eV) which causes cascade resonance in 36-ZGNR (N z = 36) to study the HHG process. In Figs. 5(a) and 5(b), the plateau of harmonic spectrum merely extends to 17th order, which is evidently shorter than that in 36-ZGNR (27th order). In Figs. 5(c1)-(c3), it is found that the electron dynamics is far from the cascade resonance. On the subbands with energies above 3.5 eV, the population is not as significant as in the case of a cascade resonance. This is because the relation between subband gap and laser frequency (∆ sg ≈ ω/2) does not agree with the frequency condition in Eq. (8a). Nevertheless, the subband gaps are not uniform. The gaps near the Fermi energy are larger than those higher, which can satisfy the frequency condition to some extent. Therefore, we observe that subbands close to the Fermi energy are involved in a near resonance. The cascade resonance for 80-ZGNR, however, appears if we employ a laser with wavelength λ = 6 µm (ω ≈ 0.21 eV) satisfying the frequency condition (∆ sg ≈ ω). The accumulation zones of electrons span from 0 to 6.5 eV and harmonic cutoff exceeds 11 eV (not shown).
When the laser field does not reach the threshold strength [Eq. (10)], the electron dynamics in GNRs is similar to what has been explored in the ZGNR-like model, see Sec. III B. The electrons are not able to be excited to high-energy subbands like what has been presented in Figs. 5(c1)-(c3). The intensity of harmonic spectrum is lower and the plateau is shorter than those in the strong laser field. In next section, the study of cutoff law on the field strength systematically shows the variation of harmonic spectrum from weak coupling regime to strong coupling regime.

IV. CUTOFF LAW OF HHG IN GRAPHENE NANORIBBONS
The cutoff law of HHG is a fundamental issue in strongfield physics. In GNRs, it has been shown that the energy cutoff increases linearly with increasing field strength and wavelength for fields along the ribbon direction, as also shown in the Appendix D. This is consistent with the generic observations in bulk solids. When the laser polarization turns to the cross-ribbon direction, however, the accelerated motion of electrons on the energy bands switches to the electron dynamics for finite systems like cascade resonance, and thus leading to alteration of the cutoff law as well. In the following, we explore the cutoff law of GNR-HHG on the laser duration, field strength and wavelength.
Duration. - Figure 6(a) shows the HHG spectra in the 80-AGNR subjected to resonant laser pulses with different durations. In the calculation, laser pulses in form of Eq. (A10) are employed, so the number of laser cycles (n cyc ) can be used to indicate pulse duration. We note that the plateau of harmonic spectrum is gradually stretched with increasing the laser duration. In Fig. 6(b), the harmonic cutoff clearly exhibits linear scaling with pulse duration for n cyc 5 and reaches the maximum attainable energy at n cyc = 11. Between n cyc = 6 and n cyc = 10, the cutoff varies slowly with increasing the laser duration. This cutoff law can be interpreted qualitatively by the two-stage excitation process in GNRs under the cascade resonance as elaborated in Secs. III A and III B. The linear increase and the slow variation of the cutoff corresponds to the rise and oscillation stages in the resonant excitation, respectively.
Field strength. -In Fig. 6(c), we study the field strength dependence of harmonic emissions in the 40-AGNR interacting with resonant pulses. The blue line represents the critical field strength for the 40-AGNR in which E cri = 8 × 10 −4 a.u. (i.e., I cri = 2.24 × 10 10 W/cm 2 ), and the red line is located at the field strength (E sec = 4×10 −4 a.u.) in which harmonic spectrum starts to possess marked second plateau. Given these two lines, the spectrum is intuitively divided into three regions: (i) γ R 1, (ii) γ R 1, and (iii) γ R > 1. (i) In the weak coupling regime (γ R 1), the energy cutoff presents a linear dependence on the field strength, see Figs. 6(c) and 6(d). (ii) When entering the transition regime (γ R 1), we can observe two-plateau structures on HHG spectra like the red dashed curve in Fig. 6(e), which result from the difference of population on lower and higher subbands. This difference can be understood by the similar electron dynamics in the ZGNR-like system, that is, the lower subbands participate in collective resonance, whereas the higher subbands are in near resonance or even in non-resonance. As Fig. 6(d) shows, the cutoff of the first plateau (E (1) cut ) preserves the linear dependence in region (i), but the second cutoff (E (2) cut ) approaches the attainable maximum tardily. (iii) As the field strength increases up to E cri , the blue curve in Fig. 6(e) shows that the first and the second plateaus merge thogether since all subbands are involved in the collective resonance. For this reason, the energy cutoff in the strong coupling regime (γ R > 1) tends to saturate and remains at the maximum band gap. We have thus revealed that the cutoff law of HHG in GNRs depends strongly on the interaction regime to which the system is subjected, the cutoff neither following the linear scaling in electric field of the bulk solids nor the quadratic dependence in electric field of gas. Wavelength.
- Figure 6(f) shows the HHG spectrum as a function of driver wavelength in a 40-AGNR (∆ sg ≈ 0.4 eV). It can be found that the cutoff frequency is extended abruptly around wavelength λ = 3 µm (ω ≈ 0.41 eV) where the frequency condition is satisfied exactly. It is because the cascade resonance appears as the wavelength approaches 3 µm. When the driver wavelength increases away from 3 µm, the laser-driving AGNR begins to deviate from the resonance point, and thus the plateau of the harmonic spectrum shrinks drastically.

V. DISCUSSION OF CASCADE RESONANCE IN GENERAL CONFINED SYSTEMS
We now turn our attention to the possibility of observing the cascade resonance induced HHG spectra in general confined systems. The key ingredients for cascade resonance, viz., the frequency-matching and field strength conditions, can be approximately met in a variety of systems in a nanoribbon geometry [61], for example in the Kagmoe lattice [87] and α − T 3 lattice [88,89] systems. Additionally, quantum dots offer a high degree of tunability in terms of the the energy spectra via the synthesis route, size dependence and/or the external parameters. In general, however, for most materials/systems not only the frequency-matching condition over a large number of subbands may not be possible, but also the fundamental gap (E g ) could be significantly different from the (average) subband gap (∆ sg ).
To explore such systems for possible cascade resonance induced HHG phenomena, we resort to the simplified models akin to the ones introduced in Sec. III B, but the dipole matrix element takes the form of We consider E g = ∆ sg and laser frequency matched to either the fundamental gap or the average sub(band) gap. Therefore, we need two parameters denoting the coupling strength for fundamental gap and sub(band) gap, respectively. We additionally consider cases where ∆ sg is non-uniform. To illustrate these cases, we also consider a 20-level model. Case I: E g < ∆ sg -For brevity, we consider the fundamental gap E g = 1 a.u. while all other energy levels are equally-spaced in energy with a larger gap: When the laser frequency matches with the fundamental gap E g , the system exhibits near resonance irrespective of the coupling strength. The electron dynamics in the strong coupling (E 0 = 2 a.u., γ E = 2) is shown in Fig. 7(a1). Significant electron population and oscillation merely occur at the energy levels close to the Fermi energy. On the other hand, when the laser frequency matches the sub(band) gap ∆ sg , the cascade resonance can be achieved, in principle, even for the weak coupling between the sub(band) levels, as shown in Fig. 7(a2) for γ ∆ = 0.125 (E 0 = 2 a.u.). The population on each energy level is comparable. The corresponding HHG spectra for different tuning of the laser frequency are shown in Fig. 7(a3) We clearly see that the HHG spectrum for the cascade resonance when ω = ∆ sg is much intense than for the near resonance when ω = E g . The frequency cutoff of resonant harmonic spectrum reaches 72 a.u. which is close to maximum band gap (73 a.u.).
Case II: E g > ∆ sg -We set the fundamental gap and sub(band) gap as otherwise. (16b) Figure 7(b1) shows the population dynamics for the laser frequency matching the fundamental gap E g . The field strength is set by E 0 = 1.0 a.u. (γ E = 0.625). The near resonance which has been observed in Fig. 7(a1) occurs in the energy levels near the Fermi level. Therefore, in the corresponding HHG spectrum, see Fig. 7(b3), the intensity is relatively lower and the plateau structure is not clear. On the other hand, when the laser frequency matches the sub(band) gap, a cascade resonance can occur only when γ ∆ 1. Figure 7(b2) shows the electron dynamics when E 0 = 2 a.u. (γ ∆ = 2). The cascade resonance takes place. The harmonic cutoff, blue curve in Fig. 7(b3), can thus reach maximum energy gap at 22 a.u. (E 20 − E 1 ). We further note that the cascade resonance depends crucially on the ratio of E g /∆ sg . For much larger fundamental gap values, the electron density in the higher subbands decreases as the excitation of electrons to the first conduction level is probabilistically low. The cascade resonance then cannot take place. In the specific context of GNRs, this can be achieved by adding staggered onsite potentials which open or enlarge the fundamental gap.
In Appendix E, we show the HHG spectrum and population dynamics in 36-ZGNR with a 1.0 eV gap. We can still observe the cascade resonance near the Dirac points but the lower conduction population in comparison to that of gapless 36-ZGNR. The intensity of corresponding HHG spectrum is lower than gapless ZGNR-HHG spectrum. However, the cascade resonance disappears when the fundamental gap is larger than 1.5 eV (not show).
A more realistic scenario in materials is when the sub(band) gaps are non-uniform. To account for these effects, we consider a model where ∆ sg is chosen randomly (uniform probability distribution) in a given energy range, while E g > ∆ sg , where . . . denotes the average value. The energy levels are: where Ran(n, m) means a random number between n and m. The laser frequency matches the average sub(band) gap (ω = ∆ = 1 a.u.). Figure 7(c1) shows the electron dynamics in the weak coupling regime, where E 0 = 0.2 a.u. (γ ∆ ≈ 0.2). Because of the large fundamental gap (E g = 4 a.u.), the electrons can hardly be excited to conduction levels. The population of the first conduction level is even less than 0.01. In comparison, in the strong coupling regime (E 0 = 2 a.u.), shown in Figure 7(c2), all the higher lying energy levels have sizable population. The details of the energy-level distribution seem immaterial to the cascade resonance phenomena. These features are also clearly reflected in the HHG spectra. In Fig. 7(c3), the intensity of harmonic spectrum in strong coupling regime is much higher than that in weak couling regime.
Consequently, the cascade resonance can be anticipated in real nanoribbon materials where a few subbands are nearly equi-spaced in energy while other bands are randomly distributed. When laser frequency matches subbband gap, the resonant excitation allows more electrons to occupy higher conduction subbands, which eventually enhances and broadens the harmonic spectra in experimental measurements.

VI. CONCLUSION & OUTLOOK
In conclusion, we demonstrated theoretically that the cascade resonance provides a systematic way to extend and enhance the harmonic spectrum of a confined quantum system. The cascade resonance can pump electrons to higher subbands, resulting in enhanced higher-order harmonic emissions, which is, therefore, fundamentally distinct from the high-harmonic generation mechanism of atoms, molecules and bulk.
Based on the study of size-dependent GNR-HHG and dynamic analysis in a multilevel model, the resonance conditions for laser parameters like frequency and field strength are established. The cascade resonance occurs when the laser frequency matches the subband gap and the field exceeds the threshold strength. While large deviations from the ideal frequency-matching condition leads to disappearance of the cascade resonance, for small deviations, strong field strengths can still induce cascade resonance although the strength of the HHG spectra becomes relatively weak.
These predictions are well-within the experimental reach, and relatively straight forward to verify in present experimental setups. These ideas can also be applicable to other confined materials/systems which, in general, may not meet the frequency matching condition perfectly. Perhaps, the most interesting situation is when only a part of the subbands meet the frequency matching condition, where the cascade resonance of few levels may take place by carefully tuning the laser frequency and intensity, for example in other two-dimensional materials in a nanoribbon geometry withstanding, or in quantum dots. Among these, the nanoribbons based on kagome lattice and α − T 3 lattice systems are particularly appealing due to the presence of flat band together with the linearly dispersing Dirac bands. A detailed study of HHG phenomena in such system is, however, beyond the scope of this work.
Looking forward, our theoretical framework for attosecond physics of confined systems establishes cascade resonance as a powerful tool for obtaining intense simple ultraviolet/extreme-ultraviolet X-ray sources. Attosecond technology has previously focused on the modulation of ultrafast processes by the intensity, wavelength, polarization, and time delay (two pulses) of the incident laser. The resonance mechanism will provide a more diverse manipulation approach because of the sensitivity to material size and drive duration, which may finally allow the establishment of highly tunable solid-state XUV sources. Furthermore, the relation between HHG and cascade resonance in the nanoribbon system will provide a new platform and idea for the study of carrier-wave Rabi flopping.

ACKNOWLEDGEMENTS
Xiao Zhang acknowledges fruitful discussions with Jinbin Li. The authors acknowledge supported from the National Natural Science Foundation of China (NSFC) (Grants No. 11834005, No. 11874030, No.11904146, No. 12047501).

Appendix A: Methods
The atomic structures of GNRs with armchair edges and zigzag edges are presented in Fig. 1(a) and 1(b), respectively. The distance of unit cell for AGNR and ZGNR are d a = √ 3a and d z = a, respectively, where the lattice constant a = 2.46Å.

Tight-binding model for graphene nanoribbon
The tight-binding Hamiltonian for AGNR takes the form wherep † (q † ) andp (q) are creation and annihilation operators corresponding to the chain P (Q), i (j) labels the site of an atom in the x (y) direction, γ = 3.03 eV is the hopping integral, and N x is the number of unit cells in the x direction. Using the similar method, the ZGNR Hamiltonian can be expressed bŷ Since we have periodic boundary conditions in the x direction, the Fourier transform can be made bŷ where k x is the quasi-momentum of the x direction, x i is the atomic position in the x direction, and kx∈BZ is the summation over the Brillouin zone (BZ where hopping integrals γ a . We obtain the AGNR Hamiltonian in the form ofĤ a = kx Φ a † kx H a kx Φ a kx using bases of Φ a kx = (p kx,1 ,q kx,1 ,p kx,2 ,q kx,2 , · · · ,p kx,Na ,q kx,Na ) T , where the matrix notation of AGNR Hamiltonian H a Likewise, we can write the ZGNR Hamiltonian in terms of Φ z kx = (ĉ kx,1 ,ĉ kx,2 , · · · ,ĉ kx,Nz ) T , in the from ofĤ z = kx Φ z † kx H z kx Φ z kx , where H z kx is written as The Hamiltonian H a kx is a 2N a × 2N a matrix because the AGNR contains 2N a sites per unit cell. For the ZGNRs, however, Hamiltonian H z kx only has N z orthogonal eigenstates, because each unit cell is composed of N z sites.
We solve time-independent Schrödinger equation obtaining the eigenstates S i (k x ) and eigenvalues E i (k x ) of the system. In Figs. 1(c) and 1(d), we show the energy bands of AGNR and ZGNR in the case of N a = 10 and N z = 10, which are calculated by plugging Eqs. (A6) and (A7) into Eq. (A8), respectively.

Coupling to an external laser field and numerical calculations
In order to simulate the interaction of laser field and GNR, we use the time-dependent Schrödinger equation where H(t) is the matrix form of time-dependent Hamiltonian, which has been coupled to the external field A(t).
Here, the vector potential A(t) = (0, A y ) T reads in the domain 0 t n cyc T , where E 0 is the amplitude of the electric field, ω is the fundamental frequency, T = 2π/ω is the time period, and n cyc is the total number of laser cycles. For a tight-binding GNR, the driving fields shift the quasi-momentum k x like because of the periodic boundary condition in the x direction [90]. The y component of the drivers will induce phase factors on the hopping terms, which makes the replacement of where γ ij represents the hopping integral between site i and site j [90]. Thereby, according to the above statement, the matrix elements of the time-dependent AGNR Hamiltonian H a (t) are expressed by where the phase factor γ a y (t) = e iaAy(t)/2 and H a ij (k x ) correspond to the elements in Eq. (A6). Similarly, the elements of ZGNR Hamiltonian H z (t) read where γ z y,1 (t) = e and H z ij (k x ) correspond to the elements in Eq. (A7). It is important to note that only the upper triangular elements of H(t) are given in Eqs. (A13) and (A14), and one can easily obtain the values of lower triangular elements by using H ji = (H ij ) * . In this simulation, the initial states Ψ nkx (0) are determined by the half lowest eigenstates of the Hamiltonian.
The wave functions are numerically propagated with the Crank-Nicolson method [91]: where I is the identity matrix, and ∆t is the discrete time step for temporal evolution. Then we can calculate the value of induced electric current as where n labels the state, N is the size of the Hamiltonian matrix, α = x, y, and Ψ † nkx = Ψ * nkx T . Ultimately, the harmonic spectrum along α direction can be evaluated from the current by Note that the derivative of current is multiplied by a Blackman window before the Fourier transform.

Appendix B: Optical selection rules in GNRs
Since the selection rules of GNR-subbands [92][93][94][95][96] are crucial for understanding the size-dependent HHG spectra in Sec. II A, it is necessary for us to clarify them in this section.
The dipole matrix element between subbands J 1 (s) and J 2 (s ) reads wherer α is position operator along α direction,p α is the momentum operator along α direction, s and s are the band type indices, J 1 and J 2 are the subband numbers, and E J is the band dipersion of subband J. These indices have been introduced in Sec. A 1. In what follows, only the y component of the dipole matrix element d y J1s,J2s is discussed, because we mainly focus on the optical transition while the GNR is driven by y-polarized light. In Eq. (B2), we give a example of dipole matrix (modulus) of 10-ZGNR at Dirac point K z ,  in the basis (|5v , |4v , . . . |1v , |1c , . . . |5c ) T . We can see some zero matrix elements which result in the interesting selection rule in ZGNRs.
Firstly, let us see the selection rule in AGNRs. According to the analytical derivations in Refs. [94][95][96] and our numerical calculations, the dipole matrix elements of intraband from the same group, i.e., s = s, and (J 1 , J 2 ) ∈ {J} or {J }, possess the following feature: (B4) That is to say the intraband (interband) transition between those two subbands is forbidden whenever the difference in corresponding indices (∆J) is an even (odd) number. Then let us move to another case, considering the transition between the subbands in different groups, that is J 1 ∈ {J} and J 2 ∈ {J }. The intraband (s = s ) dipole moments are shown as d y J1s,J 2 s = 0, ∆J = J 2 − J 1 ∈ even d y J1s,J 2 s = 0, ∆J = J 2 − J 1 ∈ odd. (B5) For the interband transition (s = s ), the situation is just opposite to Eq. (B5), that is, the optical transition is forbidden whenever the index difference ∆J is even, shown as d y J1s,J 2 s = 0, ∆J = J 2 − J 1 ∈ even d y J1s,J 2 s = 0, ∆J = J 2 − J 1 ∈ odd. (B6) Secondly, we consider the selection rule of ZGNR which is simpler than that of AGNR, because the ZGNR only has one group of subbands. Based on the similar calculation and analysis, we find that the features of dipole matrix elements for intraband and interband transitions are just same as which has been shown in Eqs. (B3) and (B4). The intraband transition is allowed when the index difference ∆J is odd, whereas the interband transition is allowed when the index difference ∆J is even.
Focusing on the intraband transitions near the Dirac points, we compute the dipole matrix elements d y J,J+1 for AGNR and ZGNR of nearest-neighbor subbands at Dirac points using Eq. (B1). It shows an interesting characteristic that the moduli of d y J,J+1 (K) are almost the same with different band numbers J, for example, d y ∆J=1 (K z ) in Eq. (B2). So we plot | d y ∆J=1 (K a ) | and | d y ∆J=1 (K z ) | , the average values of dipole matrix elements of nearest-neighbor subbands at Dirac points, in Figs. 8(a) and 8(b), respectively. For AGNR, the results of matrix elements from the second group {J } are not presented in Fig. 8(a). We can see clearly that the dipole matrix elements | d y ∆J=1 | , for both AGNR and ZGNR, vary linearly with the number of sites: This indicates that the coupling between nearestneighbor subbands becomes stronger while increasing the ribbon width.
Appendix C: Rise time depending on total number of energy levels and coupling strength in the equi-spaced model We know that in equi-spaced model the cascade resonance occurs when frequency condition is satisfied no matter how weak the electric field is. The electrons are excited to conduction levels in a cascade way from lower to higher. We can foresee that the weaker the field strength and the more the energy levels, the longer it takes for the electrons to be excited to the highest energy level. Therefore, the rise time should depend on the total number of energy levels and the coupling strength. From the perspective of HHG, this relation indicates that the manipulation of the harmonic cutoff using laser duration is easier to achieve in a large-size nanoribbon material.
Appendix D: HHG in GNRs by applying driving fields along the ribbon direction We study the harmonic generation in GNRs by applying the driver field to the ribbon direction. Figures 10(a) and 10(b) respectively show the harmonic spectra for 40-AGNR and 40-ZGNR as a function of the field strength, and the wavelength is 3 µm. It is found that the cutoffs of harmonic spectra for both AGNR and ZGNR scale linearly with the field strength, see the white lines. It is because the dispersion of subbands of GNRs near the Fermi level is almost linear with k x . In Figs. 10(c) and 10(d), we also see the linear dependence between the cutoff energy and the laser wavelength, which agrees well with the general behavior of HHG in bulk solids. The reason is also attributed to the linear band dispersion.

Appendix E: HHG in GNRs with onsite potential
In principle, we can open a fundamental gap in GNRs by adding staggered onsite potential on the two nearestneighbor sites respectively. Figure 11 shows the HHG spectrum and conduction dynamics in 36-ZGNR with a 1.0 eV gap. Since a large fundamental gap between valence and conduction bands, the excitation rate decreases, and fewer electrons are excited to the first conduction band in Figs. 11(c1)-(c3). However, one can still observe the cascade resonance, and the electrons still accumulate near the Dirac points although the total population becomes lower in comparison to that of gapless 36-ZGNR. Correspondingly, the intensity of the HHG spectrum in the Figs. 11(a) and 11(b) is lower than that in the Figs. 3(a) and 11(b).