High harmonic generation in monolayer and bilayer of transition metal dichalcogenide

In transition metal dichalcogenides (TMDCs), charge carriers have spin, pseudospin, and valley degrees of freedom associated with magnetic moments. The monolayers and bilayers of the TMDCs, in particular, MoS$_2$, lead strong couplings between the spin and pseudospin effects. This feature have drawn attention to TMDCs for their potential use in advanced tech devices. Meanwhile, high-order harmonic generation (HHG) has recently been applied to the characterization of the electronic structure of solids, such as energy dispersion, Berry-curvature, and topological properties. Here, we show theoretical results obtained with the `philosophy' of using HHG to investigate the structural effects of the monolayer and bilayers of MoS$_2$ on nonlinear optical emission. We use a simple model for MoS$_2$ in the 3R AB staking. We find that the pseudospin and valley indexes (the Berry curvature and the dipole transition matrix element) in TMDC driven by circularly polarized laser (CPL) can encode in the high energy photon emissions. This theoretical investigation is expected to pave the way for the ultrafast manipulation of valleytronics and lead to new questions concerning the spin-obit-coupling (SOC) effects on TMDC materials, Weyl Semimetals, and topological phases and transitions in topological insulators.

The process of High-order Harmonic Generation (HHG) is a highly nonlinear phenomenon in the sense that an incoming middle infrared or infrared (MIR or IR)hω 0 ∼ 0.309 eV (with wavelength λ 0 = 4 µm) interacting with the solid produces a new spectrum, as a consequence of the optical responses of the lattice to the strong laser field (illustrated in Fig. 1) [6][7][8]. The outgoing emission contain high energy photons, nhω 0 , regarding the fundamental driving laser,hω 0 , where n is typically an integer [6,8,9] (for a mathematically formal description, See Appendixes A and B). This paper explores the high-order harmonics emission from monolayer and bilayer of a typical TMDC such as MoS 2 . The monolayer and bilayer crystalline structures of MoS 2 are illustrated in Figs. 1(a) and 1(b), respectively [10][11][12]. In the case of the bilayer, the stacking is featured by the AB type stacking which is so called 3R polytype [12]. This means that the 'S1' atom of Layer1 interacts with an 'empty space' (dashed black link in Fig. 1(b)) of Layer2; and, the 'Mo' atom of Layer1 interacts with an orbital of 'S2' atom in Layer2. This is the nearest neighbour interaction between inter-layers. The interactions between the 'Mo' and 'S' atoms in different layers are mediated by the hopping parameter t 11 [10,11] which is described in Section 2 and 3. The manipulation of monolayers, bilayer up to bulk TMDCs has been studied by several experimental and theoretical efforts using Angular Resolved Photoelectron-Energy Spectrum arXiv:2111.10986v2 [cond-mat.mtrl-sci] 25 Nov 2021 (ARPES) techniques [13,14]. However, its complete characterization remains yet a downright challenging research field. Returning to the key physical insight of HHG in solids, the simplified description of the HHG mechanism is depicted in Figs. 1(c) and 1(d) for both monolayer and bilayer MoS 2 [11]. The accelerated carriers in the conduction (valence) band have two intrinsic mechanisms [7]: intraband current and interband current. Note that our theory is subjected to the Keldysh approximation in which 0 ≥hω 0 , where 0 is the band gap of the crystal (See Appendix B). Thus, the interband mechanism of HHG can be interpreted in terms of three simple steps [15][16][17]: (i) at a time t when the driving laser-field, E(t), reaches its maximum, an electron-hole (e) and (h) pairs is created (see Figs. 1(c)-1(d)) by the action of the oscillating external electric field E(t), (ii) between t and t, the e and h are accelerated or propagated under the quasi-classical action S(k, t , t) Eq. (A2) in the bands of the MoS 2 (See Appendixes A and B.2), gaining a 'considerable' amount of energy from the laser and bands, and simultaneously encoding rich information of the lattice structure in the crystal (see Fig. 1(a) and Fig. 1(b)) [3,7,15,16,18]. 'Finally', (iii) after these e and h are accelerated on the Brillouin zone (BZ), the e and h can find a time t at which they annihilate each other. Through the conservation of energy [15], this physical picture converts the accumulated e-h energy into a high energy photon with respect to the fundamental-laser as depicted in Figs. 1(c) and 1(d). Surprisingly, this HHG process encodes rich information about the atomic structural, electron dynamics in the TMDCs. This opens the door to investigate fundamental questions and technological development such as symmetries of the Hamiltonian in the solids, the fermiondynamics, the effects of the Berry curvature, pseudospin, spin-valley, topological phases and transitions [16,[18][19][20][21][22][23]. Additionally, these features are relevant to the study of the fundamental protections of the symmetries, carrier measurements, Bloch oscillations at the natural timescale of the electrons, sub-femtoseconds (10 −18 s) [24]. In addition, it has been demonstrated that in TMDCs, the number of layers and the polytype geometrical features modify the electron structure of the TMDCs and their transport features [10,11]. This turns MoS 2 from a simple semiconductor to a material with a strong-spin orbit coupling (SOC) effect [11,25,26].
In this paper, we use HHG spectroscopy to explore the effects of interlayer coupling in bilayer MoS 2 and to investigate how the pseudospin and valley-index can differentiate the monolayer from the bilayer. Furthermore, we pursue to understand whether or not the Berry Curvature and Dipole Transition Matrix Elements (DTMEs) have an impact on the harmonic emissions [25]. Pioneering studies of HHG spectra from monolayer and bilayer graphene have shown qualitative differences [27,28] and, in addition, shown an attractive un-typical enhancement of the emissivity as a function of the ellipticity of the driven laser, particularly when the fundamental field is elliptically polarized. Similar effects have also been theoretically observed by Tancogne-Dejean in bulk MgO [20] by the fully time-dependent density functional theory [29]. Here, we calculate the HHG spectra, using the time-dependent density matrix and analyze them using a tight-binding model (TBM), described in Sections 2 and 3. The Berry curvature around the K/K under the k · p approach yields [10,11,25,30]: The parameters in Eq. (1) are described in Ref. [11] for the monolayer of TMDC. Note, M 0 is related to the onsite potential of the lattice ( Fig. 1(a)), t 1 is nearest neighbor (NN), and a 0 is lattice constant and τ the valley index for this simple monolayer model. this is the Berry curvature for the monolayer system MoS 2 of our TBM up to a linear expansion of the field-free HamiltonianĤ 0 (k) obtained from Refs. [18,19]. According to Kormányos et al. in Ref. [10] the bilayer Berry curvature about the K/K', for the valence band, is written as The parameters τ, δE ll , are defined according to the spin-valley index and the onsite energies Ref. [10]. λ 1 and λ 2 are related to the hoping strength between the layers of MoS 2 .
Interestingly the bilayer Berry curvatures shows dependence of the interlayer interaction strength. Suggesting the main hypothesis of this paper: "whether or not HHG can encode information of the monolayer peudospin and valley indexes once the system is subjected to strong middle infrared (MIR) lasers." We numerically examine the HHG spectra in terms of the Berry Curvature and the DTMEs for both monolayer and bilayer MoS 2 . The results show an interesting interconnection between the selection rules as well as a particular enhancement in the HHG spectrum for a few harmonic orders (HOs). We define HO as the ratio between the emitted photon frequency ω and the freq. of fundamental driver ω 0 , i.e., HO = ω/ω 0 . Further, The analysis of angular rotation, pattern of emission reveals a slight difference between mono and bilayers for the high-order harmonics. However, for the harmonic about the bandgap (Fig. A1), we find an interesting difference which might be attributable to the Berry Curvature and the DTMEs. These quantities contain pseudospin and valley indexes information of Mo 2 for monolayer or bilayer. This suggest that the angular rotation of the harmonics can be linked to DTMEs. Surprisingly, in our HHG produced by RCP light, we notice that the HOs in the plateau region exhibit a layer difference with respect to the Berry Curvature and DTMEs. We take advantage of these quantities to analyze our numerical results, the selection rules that these quantities impose and the accumulating phase of the electron and hole wavepackets at each bands (Appendixes B). In the next, we will review the total currents on the Bloch basis and split them into interband and intraband currents. Additionally, we will describe the time-dependent density matrix formalism to compute the density at the time and k-space.

Current and Time-dependent Density Matrix
The microscopic electron-charge current j(k, t) in a periodical crystalline structure subjected to an external oscillating laser E(t), and the macroscopic 'measurement' is calculated by integrating the k-elementary-microscopic currents in the BZ (atomic units are used throughout this paper unless otherwise indicated): Here we define the elementary-microscopic current j(k, t) as, This is the expectation value of the velocity operatorv = −i Ĥ (t),x whereĤ(t) is the full time-dependent Hamiltonian in the length-gaugeĤ(t) =Ĥ 0 +x · E(t). The current j(k, t) is defined in terms of the density matrixρ =ρ(k, t), and the momentum matrix element P m,n (k) = u m,k |∂ kĤ0 (k)|u n,k [18,[31][32][33]. Eq. (3) is conventionally split into two different contributions, the interband current, J er (t), (associated with the momentum or dipole matrix element transition between the valence and conduction bands) and the intraband current, J ra (t), which show how the electron-wave and hole-wave move in each band (See Fig. 1(c) and 1(d), respectively). Mathematically, those contributions are written as follows [18,19,31]: which obviously lead to the total current, To further evaluate the time-propagation density matrixρ(k, t), and the interband and intraband currents, we numerically solve, in the time-k space, the density matrix given by the Liouville-von Neumann equation, whereρ(k, t) is evaluated in the electro-magnetic length-gauge representation ofĤ(t). We will also consider the effects of the dissipation of scattering electrons in the lattice and potential thermal bath influences byD 2 . As a state-of-the-art method in high-order harmonics, the phenomenological dephasing time, T 2 , is included in the off-diagonal of the 'dissipation' D 2 . We adopt our numerical method, following Refs. [19,34,35] in the moving frame of BZ = K − A(t), where K and A(t) denote the crystal canonical momentum and the vector potential of the laser-field, respectively. The vector potential is computed according to E(t) = −∂ t A(t) [7,15]. Givenρ(k, t) matrix and the transition 'canonical' momentum matrix element, P nm (k), we compute the HHG spectra by Fourier transforming the total current (Eq. (7)) and taking its absolute square. This procedure allows us to compare the nonlinear optical response from monolayer with that from bilayer MoS 2 .

Hamiltonian model for monolayer and bilayer MoS 2
According to the tight-binding model (TBM), the laser-free field HamiltonianĤ 0 of the monolayer [25] and bilayer of TMDC, under the staking 3R-polytype symmetry for the bilayer [10,36,37] (Figs. 1(a) and 1(b)) is written as: where summatory indices n, i, j , and i, j are number of layers, the intra-layer nearestneighbour (NN) and the next nearest-neighbour (NNN) indices. n is the onside potential of each atom in the lattice.ĉ † jn Andĉ in are the creation operator at the j th lattice site and the annihilation operator at the i th atom of the n th layer (Figs. 1(a) and 1(b), respectively). t 1 , t 2 and t 11 are the hoping parameters of the NN, NNN intra-layers, as well as the inter-layer, respectively [10,25,36].

Harmonics from mono-layers and bi-layers
In the following sub-sections, we will demonstrate the high-order harmonics emitted from monolayer and bilayer MoS 2 subjected to a strong MIR laser-field. In our simulations, we use both a linearly polarized laser (LPL) and a circularly polarized laser-field (CPL). We also analyze the emitted HHG spectra as a function of the interlayer strength t 11 .

Response to linear polarized lasers
We begin presenting our numerical results by showing simulated HHG spectra from monolayer and bilayer MoS 2 driven by a linearly (or circularly) polarized laser-fields (LPL or CPL) with respect to the Γ-K direction. The results are shown in Figs. 2(a) and 2(b).
The HHG spectra from monolayer MoS 2 shows a few differences with respect to the normalized (normalization of the current with respect to the number of layers) HHG from bilayer MoS 2 . In particular, there is no substantial difference in the intensity of HHG emissions below the band gap (harmonic orders (HOs) ≤ 0 /ω 0 ). This observation is indicated in blue and red triangles. This result is surprising, since one can intuitively expect that the bilayer will enhance the emission due to more channels for tunneling excitation and interband transition. By contrast, only the 13 th order shows an enhancement of more than one order of magnitude for the MoS 2 bilayer. This can be interpreted to indicates a complex e-h recombination. In addition, Figure 2(a) shows the same cut-off regions for both monolayer and bilayer geometrical configurations. (a) HHG spectra from monolayer (blue) and 3R phase bilayer (red) MoS 2 driven by linearly polarized laser along Γ-K direction. The HHG intensity is normalized with respect to the number of layers. (b) HHG spectra with the right handed circularly polarized driving pulse for monolayer (blue) and bilayer (red). (c) and (d) Interband and intraband components of HHG spectrum for monolayer and bilayer MoS 2 , respectively, driven by linearly polarized pulse. The laser parameter used for these simulations are: λ 0 = 3.8 µm (hω 0 = 0.3263 eV), the carrier-envelop phase is fixed to zero and the number of optical cycles N cy = 10 at FWHM under a gaussian envelope. The intensity of the laser is I 0 = 0.315 × 10 12 W · cm −2 , The dephasing time is fixed to T 2 = 5.2 fs. In the spectra the harmonic order (HO) corresponding to the band gap are localized about HO6 th . According to our formalism, the e (h) can be driven in several region of the energy dispersion of the BZ for the monolayer and bilayer MoS 2 (See Figure A1).
The 'chaotic' intensity behavior observed in Fig. 2(a) is rationalized via the selection rules and accumulated phase of the radiation emission, and of course, via the breaking inversion symmetry of the Hamiltonian. These selection rules are imposed by the DTMEs, d mn (k) = i u m,k |∂ k |u n,k , and the quasi-classical action phase, S(k, t , t) (Appendix B), between conduction and valence bands for the monolayer system with respect to the bilayer system [18,19,31]. We think that the different coherences in the density matrixρ(k, t) lead to larger effects on the 13 th -order, from the interband and intraband contributions take place (Figs. 2(c) and 2(d)) [7,15]. While the electron (hole) wavepacket is driven in the conduction (valence) bands and its corresponding dipolar selection rules shown in Fig. A2 and Fig. A3 (Appendix), these "complex propagation" effects can easily lead to constructive or destructive interfering channels.
These results exhibit even and odd HOs and illustrate the selection rules because the interband currents are governed by the dipole matrix element [35] which breaks the inversion symmetry of the laser-field free Hamiltonian,Ĥ 0 (k). In case of the intraband current, the HHG spectra show a complete symmetry behaviour in the sense that only odd HOs show up in the spectra.

Response to circularly polarized lasers
The HHG produced by a contra-wise clock of the circularly polarized laser are shown in Fig.  2(b). For the different geometrical-phases, i.e., monolayer and bilayer, HHG spectra exhibit interesting differences.

1.
A surprising enhancement by almost one order of magnitude in the harmonics emitted around the plateau region is observed between monolayer and bilayer (See blue triangles), 2.
All the harmonics follows the 3-fold symmetry of the system with the co-rotating (3n + 1) and contra-rotating (3n + 2) harmonic orders in the plateau region, as has already been experimentally observed in Ref. [38] in solids with a 3-fold rotational symmetry. Our results fully capture these selection rules. As such, it is clear that circularly-polarized drivers have significant potential in bringing out the signatures between the monolayers and bilayers configurations, their differences and in particular, the valley index and pseudo-spin information as a function of the layers [10].
As the selection rules in the HHG spectra are well described for the LPL and CPL, inversion symmetry breaking is manifested in the HHG spectra. This implies that the Berry curvature and the transition dipoles have an interesting effects on harmonic emission [16,19,39] for both monolayer and bilayer MoS 2 . The results shown in Fig. 2(b) can be attributed to Berry curvature effects about the K'-point in the BZ and the 'rotation of DTMEs' coupled to the CPL and the action phase S(k, t , t) (Appendix B.1). Since the dipolar selection rules mn (k) also govern the photon emission via P (±) mn (k) (the momentum matrix element is proportional to the dipole transition matrix element d mn (k), obviously the indexes m = n [19,35]), the RCP or LCP will excite the e-h about the K' and K point respectively. Note, however, that RCP-light will only excite electrons around the K' points of the BZ as clearly indicated in Figs. A2-A5. The mathematical definition of CPL selection rule weighs the Berry curvature around the excitation e (h) wave packet and the recombination of the electron and hole. Thus, we find the potential physical picture which explains approximately the enhancement of the harmonics in the plateau region in monolayer MoS 2 compared to the bilayer.

Interlayer strength in high-order harmonics
Here, we discuss the role of the inter-layer parameter t 11 in the HHG spectra. Fig. 3 reports the HHG spectra as a function of the inter-layer interaction t 11 . We find small differences in the range of t 11 = [0, 0.55] eV. This range of t 11 is logically expected to be on the order of the magnitude of the NN parameter t 1 . Naturally, the coupling between the atoms 'Mo' of Layer1 and 'S' of Layer2 is not strong enough within our approximation to observe dramatic modifications on the HHG spectrum for the bilayer MoS 2 (See Figs. 3(a), 3(b) and 3(c), as well as their corresponding comparisons between the approximated interaction strength of t 11 ∼ 0.27, 0.46 and 0.05 eV, respectively). However, note that once the interaction strength t 11 increases up to 2-3 eV, dramatic differences are observed between monolayer and bilayer MoS 2 . We did not present or address this result, since it is unlikely that the interplay hopping interaction will reach such a large energy value of 2 or 3 eV.

Structural angular rotation of the high-order harmonics
Now, We explore the angular rotation of HHG emission for low HOs (those below the band-gap), the plateau and the cut-off region as a function of the laser wavelength. We also study the response of the HHG spectra as a function of the ellipticity of the laser for the monolayer and bilayers of MoS 2 . This interestingly links to the pseudo-spin and valley effects of our model into the HHG signal. The angular rotation of the harmonic orders is another interesting quantity which might contain information of the band energy structure, the dipole matrix element and the Berry curvature as well [2,16,24,40]. Therefore, we calculate the angular rotation of the harmonics using the following procedure: (1) the electric field, E(t), of the driving laser is linearly polarized at an angle of θ 0 with respect to the Γ-K crystal orientation (x−direction on the k x , k y crystal momentum space), (2) we project the resulting charge current components J x (t) and J y (t) with respect to the laser orientation. As a result, we focus our attention on the parallel J (ω) and perpendicular J ⊥ (ω) components with respect to the violet laser polarization [24].
We then compute the angular rotation of the harmonic orders for the parallel and perpendicular components as a function of the laser wavelengths for low HOs, bandgap and cutoff region of the HHG spectra. The results of the total HHG, I HHG (ω, θ 0 ) = |J | 2 (ω, θ 0 ) + |J ⊥ | 2 (ω, θ 0 ) are shown in Fig. 4. The 2 nd -order shown in Figs. 4(a,d,g) for different laser wavelengths can be understood in terms of DTMEs and its selection rules. The DTMEs show, in Figs. A2 and A3, an 'start'-like structure for mono and bilayer MoS 2 . The difference comes from the selection rules and the laser via the Bloch theorem in the BZ, i.e., k = k 0 + A(t) − A(t ), where, k 0 and t are the excitation crystal momentum and time, t the emission time. A(t) is the vector potential of the laser field E(t) = −∂ t A(t). Additionally, we can understand the difference easily by applying perturbative theory to our formalism. Note, however, that this is out of the scope of the current paper.
The most interesting aspect of this angular analysis, in between mono and bilayer MoS 2 , is: the angular rotation of the harmonic orders around the band gap depicted in Figs. 4(e) and 4(h), specially for driven wavelength λ 0 = 3.8 and 5.0 µm at HO6 th and HO7 th , respectively .

Ellipticity dependence of low-and high-order harmonics
In this section we present the result of the HHG spectra as a function of the ellipticity ε of the driving-MIR laser. This ε takes values between ε = [−1, 0, −1], meaning left-circularly polarized-laser (LCP), linear-polarized laser (LPL) and right-circularly polarized-laser (RCP).  Figure 5 show the calculated HHG spectra for co-rotating harmonic order (HO7 th ) and contra-rotating harmonic order (HO8 th ) between monolayer and bilayer MoS 2 . Two harmonic orders, HO7 th and HO8 th show interesting differences as a function of the ellipticity of the incoming laser for both the monolayer and bilayer MoS 2 . In Figs. 5(a,c) and 5(b,d), the ellipticity of the HO7 th with respect to the ellipticity of the laser is shown. We observe almost the same tendencies with the difference that the emission yield from the monolayer is larger with respect to the bilayer emission for the LPL. The harmonic yield of the LCP vs RCP are almost zero in both monolayer and bilayer cases (Figs. 5(a,c)).
The HO8 th , contra-rotating order with respect to the incident RCP (LCP), shows an impressive enhancement at ε = ±1 (Fig.5 (e)) for the emission from monolayers. In contrast, the yield of HO8 th for the bilayer MoS 2 shows a dramatic reduction of 67% in comparison to the monolayer case ( Fig.5 (g)).
Finally, the contra-rotating harmonic order exhibits totally different ellipticity from the co-rotating HO. This is expected, according to the selection rules (Figs. A4, A5, and A6). The quasi-classial action phase, S(k, t , t) (See Appendix A2), allows us to link this signature to pseudo-spin localized at K' or K valley. The valley index τ can be associated to the ellipticity analysis of the HOs of Fig. 5(f,h) with signs of the Berry curvatures for mono and bilayer MoS 2 at the K or K' points.

Conclusions and outlook
In summary, our theory shows that the high-harmonic generation (HHG) spectrum is capable to: • observe difference on the nonlinear optical emission about the band gap between monolayers and bilayers of TMDCs.
• describe a unique difference between the angular rotations and ellipticity dependence of the emitted harmonics as a function of the number of layers concerning the ellipticity of the laser, via the selection rule of the dipole matrix element and Berry Curvature to a linear and circularly polarized laser beam.
• is susceptible to breaking the inversion symmetries (IS, respectively) and thus sensitives to the Berry Curvature and its pseudospin character.
These features promise the applications of the TMDC into spin-valley manipulation of the information, thanks to its geometrical and pseudo metallic properties. The TMDC in 2D can especially exhibit topological insulating phases, making them attractive for quantum computing, storage of information, and other applications [25]. Finally, it is important to note that all of these ideas can be applied and tested in quantum simulators [41], in particular using ultracold atoms and lattice shaking [42], as well as Rydberg atoms, polaritons, or circuit QED [43][44][45][46][47][48]. The usage to study strong-field phenomena is an emerging application of the spintronic and valleytronic technological advances [49][50][51]. Since the pseudospin is related to magnetic properties and is naturally suited to the study of systems in ordered lattices, it should provide additional clarity to the role of material with pseudospin couplings and valley indexes in high-harmonic emissions.
the high-order harmonic spectrum are dominated by the inteband currents. [7,18]. We then can discover a unique physical feature by the interplay between the energy dispersion and the interband currents.

Appendix B. Keldysh approximation and quasi-classical analysis
Here we can gain more physical insight into the physics of the HHG spectra by applying the so-called Keldysh approximation, as discussed by Vampa et al. [7]. This approximation reads: ∑ m ρ m (k, t) ∼ 0, where m is an integer index that stands for the number of conduction bands. Thus, the Keldysh approximation is valid in the limit that a 0 E 0 ∼ 0 and 0 hω 0 . This essentially means that the population transferred to the conduction band is very small compared to that remaining in the valence band. This approximation is very similar to the one used in the Strong Field Approximation (SFA), which was originally developed for atoms and molecules [6,18,52,53] -we will thus hereafter term it SFA. We focus on the discussion of the inter-band current in this appendix.

Appendix B.1. Inter-band current
In this approximation, we will restrict ourselves to the harmonic radiation produced by the interband current according to Ref. [18], We obtain a closed form of the expression for the i th vectorial-component (i = x, y) within the Hamiltonian matter-gauge [34] for 2D materials, where S(K, t, t ) is the so-called quasi-classical action for the electron-hole, which is defined as Here, j = x, y indicates the component of the electric field and the transition-dipole product which depends on the polarization of the driving laser. Expressions (A1) and (A2) are direct analogues of the Landau-Dykhne formula for HHG in atoms, which was derived in Ref. [53], following the idea of the simple man's model [52] (See [6] for a recent review). Below we will analyze these expressions using the saddle point approximation over crystal momentum to derive the effects of the Berry curvature on the relevant trajectories.

Appendix B.2. Quasi-classical electron-hole pair quantum paths
Assuming that the exponentiated quasi-classical action e −iS = e −iS(K,t,t ) oscillates rapidly as a function of the canonical crystal momentum K, one can apply the saddle-point approximation to find the points K s at which the integrand's contributions to the inter-band current ( A1) concentrated. These are solutions to the saddle-point equation ∇ K S(K, t, t )| K s ≈ 0, which can be rephrased as Two different trajectories are identified from the last equation, the first one is related to the excited electron ∆x c (K s , t, t ) in the conduction band, whereas the second one involves the  trajectory ∆x v (K s , t, t ), followed by the hole in the valence band. We then obtain a general m th trajectory for the electron (m = c) and hole (n = v), which is expressed as  where (−1) m is the alternating sign (−1) c = +1 and (−1) v = −1, and the group velocity of the m th band is v gr,m = ∇ K ε m . Here we recognize the Berry curvature Ω m as well as the anomalous velocity v a,m , which are given by Ω m = ∇ K × ξ m and v a,m = E(t) × Ω m , respectively, for the electron-hole trajectories of Eq. (A4). The previous expression can be rewritten as These electron-hole pair trajectories, along with the saddle-point condition of Eq. (A3), should produce complex-valued solutions for K s , as is the case for HHG of gases. However, determining the solutions K s is not a trivial task, as it depends explicitly on the geometrical features, such as the Berry curvature and connection, and the phase of the dipole matrix elements. This gets further complicated as the eigenstates and eigenvalues that make up the energy bands are expected to exhibit branch cuts and branch points connecting the two bands [54,55]. Once the momentum is allowed to take on complex values (as it does in complex band structure theory [56]), this leads to a nontrivial geometrical problem with a high dimensionality whose analysis requires detailed attention. Nevertheless, these saddle points K s should have a component perpendicular to the driving laser field E(t) (in the case of linear drivers), which appears as a consequence of anomalous-velocity features and of the Berry curvature Ω m (k) in particular.

Appendix B.3. Berry Connection, Berry Curvature, and Chern number
To calculate the radiation-interaction features of these eigenstates, the DTME between each state of from conduction to valence bands or vice versa. This DMTE is unique, since it contains the full information of the selection rules for the emitted radiation. We define it in the basis of the Bloch function [57], Φ m,k |x|Φ n,k = −i∇ k δ mn δ(k − k ) where d mn (k) = i u m,k |∇ k |u n,k is a regular function which encodes the momentum gradient of the periodic part of the Bloch functions. For n = m, d mn (k) defines the DTME as also have been noticed in the text. The main diagonal of Eq. (A7) elements is the Berry connection of the m th -band, which is responsible for the parallel transport of wave function phase around the band. This parallel transport is measured by the Berry curvature, given by the gauge-invariant curl This parallel transport is observed in the quasi-classical action phase of the Eq. (A7).

Appendix B.4. Dipole moment and Berry curvature in the Haldane model (HM)
We now turn to the dipole transition matrix element, as initially defined in Eq. (A7), and its deep relationship with the Berry curvature which is given here as [18] On the one hand, we find that the cross product of the HM dipole Eq. (A7) yields Expanding the Berry curvature in the Haldane model, given by Eq. (A9), one obtains [18,30] v/c (k) = ∓ We therefore conclude that [18]: which demonstrates the close relationship between dipole matrix elements and the Berry curvature in the HM. This confirmation of the relationship between the dipole product and the Berry curvature is extremely important, as it leads to a direct connection of the inter-band transition current of Eq. (A1) and the topological invariant for this model.